|
Computes the cross periodogram of two stationary time series using a fast Fourier transform.
X Vector of length NOBS containing the first stationary time series. (Input)
Y Vector of length NOBS containing the second stationary time series. (Input)
CPM (⌊N/2⌋ + 1) by 10
matrix containing a summarization of the results of the cross periodogram
analysis. (Output)
For k = 0, 1,
, ⌊N/2⌋, the (k
+ 1)-st element of the j-th column of CPM is defined as
Col. |
Description |
1 |
Frequency, ωk where ωk = 2πk/N for IFSCAL = 0 or ωk = k/N for IFSCAL = 1. |
2 |
Period, pk where pk = 2π/ ωk for IFSCAL = 0 and pk = 1/ ωk for IFSCAL = 1. If ωk = 0, pk is set to missing. |
3 |
X periodogram ordinate, IX(ωk) |
4 |
X cosine transformation coefficient, AX(ωk) |
5 |
X sine transformation coefficient, BX(ωk) |
6 |
Y periodogram ordinate, IY(ωk) |
7 |
Y cosine transformation coefficient, AY(ωk) |
8 |
Y sine transformation coefficient, BY(ωk) |
9 |
Real part of the XY cross periodogram ordinate IXY(ωk). |
10 |
Imaginary part of the XY cross periodogram ordinate IXY(ωk). |
NOBS Number of
observations in each stationary time series X and Y. (Input)
NOBS
must be greater than or equal to two.
Default: NOBS = size (X,1).
IPRINT Printing
option. (Input)
Default: IPRINT = 0.
IPRINT |
Action |
0 |
No printing is performed. |
1 |
Prints the periodogram, cosine and sine series, and the real and imaginary components of the cross periodogram. |
XCNTR Constant
used to center the time series X.
(Input)
Default: XCNTR = the arithmetic
mean.
YCNTR Constant
used to center the time series Y.
(Input)
Default: YCNTR = the arithmetic
mean.
NPAD Number of
zeroes used to pad each centered time series. (Input)
NPAD must be greater
than or equal to zero. The length of each centered and padded time series is
N = NOBS
+ NPAD.
Default:
NPAD = NOBS 1.
IFSCAL Option
for frequency scale. (Input)
Default: IFSCAL= 0.
IFSCAL |
Action |
0 |
Frequency in radians per unit time |
1 |
Frequency in cycles per unit time |
IPVER Option
for version of the periodogram. (Input)
Default: IPVER = 0.
IPVER |
Action |
0 |
Compute usual periodogram. |
1 |
Compute modified periodogram. |
Refer to the Description section for further details.
LDCPM Leading
dimension of CPM
exactly as specified in the dimension statement of the calling
program. (Input)
LDCPM must be greater
than or equal to ⌊N/2⌋ +
1.
Default: LDCPM = size (CPM,1).
Generic: CALL CPFFT (X, Y, CPM [, ])
Specific: The specific interface names are S_CPFFT and D_CPFFT.
Single: CALL CPFFT (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, IPVER, CPM, LDCPM)
Double: The double precision name is DCPFFT.
Routine CPFFT computes the cross periodogram of two jointly stationary time series given a sample of n = NOBS observations {Xt} and {Yt} for t = 1, 2, , n.
Let
represent the centered and padded data where N = NOBS + NPAD,
and
is determined by
Similarly, let
represent the centered and padded data where
and
is determined by
The periodogram of the sample sequence {Xt}, t = 1, , n computed with the padded sequence
is defined by
where
and
represent the
cosine and sine transforms, respectively, and K is the scale factor
The periodogram of the sample sequence {Yt}, t = 1, , n computed with the padded sequence
is defined by
where
and
represent the
cosine and sine transforms, respectively. Since the periodogram is an even function of the frequency, it is sufficient to estimate the periodogram at the discrete set of nonnegative frequencies
(Here, ⌊a⌋ means the greatest integer less than or equal to a). The routine PFFT is used to compute the periodograms of both
according to the version specified by the argument IPVER. The computational formula for the cross periodogram is given by
where
and
The real part of the (modified) cross periodogram represents the raw sample cospectrum and the negative of the imaginary part of the (modified) cross periodogram represents the raw sample quadrature spectrum (Priestley 1981, page 695). The relationship between the cross periodogram and its complex conjugate is given by
and may be used to recover the cross periodogram at negative frequencies.
1. Workspace may be explicitly provided, if desired, by use of C2FFT/DC2FFT. The reference is:
CALL C2FFT (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, IPVER, CPM, LDCPM, CX, COEF, WFFTC, CPY)
The additional arguments are as follows:
CX Complex work vector of length N.
COEF Complex work vector of length N.
WFFTC Work vector of length 4N + 15.
CPY Work vector of length 2N.
2. The centered and padded time series are defined by
I |
AOV(I) |
CX(j) = X(j) − XCNTR |
for j = 1, , NOBS |
CX(j) = 0 |
for j = NOBS + 1, , N |
and |
|
CY(j) = Y(j) − YCNTR |
for j = 1, , NOBS |
CY (j) = 0 |
for j = NOBS + 1, , N where N = NOBS + NPAD. |
3.
The cross periodogram IXY(ω) is complex valued in general. The
relation
IXY(−ω) =
conj(IXY(ω)) for w > 0.0 recovers the
cross periodogram for negative frequencies since real(IXY(−ω)) =
real(IXY(ω)) and
imag(IXY(−ω)) = −imag(IXY(ω)). The
periodogram I(ω) is an even function of the frequency ω. The relation
I(−ω) = I(ω)
for ω > 0.0 recovers the periodogram for negative frequencies.
4.
Since cos(ω) is an even function of ω and sin(ω) is an odd function of ω, the
cosine and sine transformations, respectively, satisfy
A(−ω) = A(ω)
and B(−ω) =
−B(ω) for
ω > 0.0. Similarly, the complex Fourier coefficients, stored in COEF, satisfy
COEF(−ω) = conj(COEF(ω)).
5. Computation of the 2 * NOBS − 1 cross-covariances of X and Y using the inverse Fourier transform of the cross periodogram requires NPAD = NOBS − 1.
Consider the Robinson Multichannel Time Series Data (Robinson 1967, page 204) where X is the Wφlfer sunspot number and Y is the northern light activity for the time period from 1770 through 1869. Application of routine CPFFT to these data produces the following results. Note that CPFFT sets CPM (1, 2) to the missing value code via routine AMACH.p<.STREF.DOC!MACH;;. The printing of CPM (1, 2) depends on the computer.
USE GDATA_INT
USE CPFFT_INT
USE WRRRL_INT
IMPLICIT NONE
INTEGER LDCPM, LDRDAT, NDRDAT, NOBS, NPAD
PARAMETER (LDRDAT=100, NDRDAT=4, NOBS=100, NPAD=NOBS-1, &
LDCPM=(NOBS+NPAD)/2+1)
!
INTEGER IPVER, NRCOL, NRROW
REAL CPM(LDCPM,10), FLOAT, RDATA(LDRDAT,NDRDAT), &
X(NOBS), Y(NOBS)
CHARACTER CLABEL1(6)*9, CLABEL2(6)*9, FMT*7, RLABEL(1)*6, &
TITLE*41
INTRINSIC FLOAT
!
EQUIVALENCE (X(1), RDATA(1,2)), (Y(1), RDATA(1,3))
!
DATA TITLE/'Results of the Cross Periodogram Analysis'/
DATA FMT/'(F10.3)'/
DATA CLABEL1/'k+1', 'w(k)', 'p(k)', 'IX(w(k))', 'AX(w(k))', &
'BX(w(k))'/
DATA CLABEL2/'k+1', 'IY(w(k))', 'AY(w(k))', 'BY(w(k))', &
'Real IXY', 'Imag. IXY'/
DATA RLABEL/'NUMBER'/
!
! Robinson Data
CALL GDATA (8, RDATA, NRROW, NRCOL)
! Center on arithmetic means
! Frequency in radians per unit time
! Modified periodogram version
IPVER = 1
! Compute the cross periodogram
CALL CPFFT (X, Y, CPM, IPVER=IPVER)
!
! Print results (First 10 rows)
CALL WRRRL (TITLE, CPM, RLABEL, CLABEL1, 10, 5, FMT=FMT)
CALL WRRRL ('%/', CPM(1:,6), RLABEL, CLABEL2, 10, 5, FMT=FMT)
!
END
Results
of the Cross Periodogram
Analysis
k+1
w(k) p(k)
IX(w(k)) AX(w(k)) BX(w(k))
1
0.000
NaN
0.000
0.000 0.000
2 0.032
199.000 184.159
3.742 -13.044
3 0.063
99.500 1364.408
35.457 -10.354
4 0.095
66.333 2433.933
29.411 39.610
5 0.126
49.750 1351.002
-21.749 29.631
6 0.158
39.800 140.421
-11.716 -1.773
7 0.189
33.167 44.117
-4.671 4.722
8 0.221
28.429 121.186
-11.003 -0.343
9 0.253
24.875 176.275
-4.782
-12.386
10
0.284 22.111
144.867 10.038
-6.642
k+1 IY(w(k))
AY(w(k)) BY(w(k)) Real IXY Imag.
IXY
1
0.000
0.000
0.000
0.000 0.000
2
1689.212 -37.480
-16.866 79.776
-552.014
3 4113.003
41.232 -49.122 1970.577
-1314.779
4 3255.785
44.214 36.068
2729.031 -690.474
5
1757.663 -8.162
41.122 1396.006 -652.513
6 1002.050
-30.107 9.778
335.410 -167.954
7
62.360 -6.825
3.972 50.636
13.678
8 1481.396
-38.096 5.487
417.288 -73.451
9
1274.161 -17.176
-31.291 469.704
-63.095
10 488.479
-12.442 -18.267
-3.570 -265.992
PHONE: 713.784.3131 FAX:713.781.9260 |