Estimates the nonnormalized cross-spectral density of two stationary time series using a weighted cross periodogram given the spectral densities and cross periodogram.
N Number of
observations in each of the appropriately centered and padded
time series X and Y. (Input)
N must be
greater than or equal to two.
SX Vector of length NF containing the estimate of the spectral density of the first time series X. (Input)
SY Vector of length NF containing the estimate of the spectral density of the second time series Y. (Input)
CPREAL Vector of length
⌊N/2⌋ +1 containing
the real part of the cross periodogram between X and Y.
(Input)
The real part of the cross periodogram evaluated at (angular)
frequency
ωk = 2πk/N is given by CPREAL(k + 1),
k = 0, 1,
, ⌊N/2⌋.
CPIMAG Vector
of length ⌊N/2⌋. + 1 containing
the imaginary part of the cross periodogram between X and Y.
(Input)
The imaginary part of the cross periodogram evaluated at (angular)
frequency
ωk = 2πk/N is given by CPIMAG(k + 1),
k = 0, 1,
, ⌊N/2⌋.
F Vector of length NF containing the (angular) frequencies at which the spectral density is estimated. (Input)
WT Vector of
length NWT
containing the weights used to smooth the periodogram. (Input)
The actual weights are the values in WT normalized to sum
to 1 with the current periodogram ordinate taking the middle weight for NWT odd or the weight
to the right of the middle for NWT even.
COSPEC Vector of length NF containing the estimate of the cospectrum. (Output)
QUADRA Vector of length NF containing the estimate of the quadrature spectrum. (Output)
CRAMPL Vector of length NF containing the estimate of the cross-amplitude spectrum. (Output)
PHASE Vector of length NF containing the estimate of the phase spectrum. (Output)
COHERE Vector of length NF containing the estimate of the coherence. (Output)
NF Number of
(angular) frequencies. (Input)
F must be greater than
or equal to one.
Default: NF = size (F,1).
NWT Number of
weights. (Input)
NWT must be greater
than or equal to one.
Default: NWT = size (WT,1).
Generic: CALL CSWEP (N, SX, SY, CPREAL, CPIMAG, F, WT, COSPEC, QUADRA, CRAMPL, PHASE, COHERE [, ])
Specific: The specific interface names are S_CSWEP and D_CSWEP.
Single: CALL CSWEP (N, SX, SY, CPREAL, CPIMAG, NF, F, NWT, WT, COSPEC, QUADRA, CRAMPL, PHASE, COHERE)
Double: The double precision name is DCSWEP.
Routine CSWEP estimates the nonnormalized cross-spectral density function of two jointly stationary time series using a fixed sequence of weights given the modified cross-periodogram and spectral densities of the appropriately centered and padded data
and
for t = 1, , N. The routine CPFFT may be used to compute the modified periodograms
and cross-periodogram
over the discrete set of nonnegative frequencies
(Here, ⌊a⌋ means the greatest integer less than or equal to a.) Either routine SSWP or routine SWEP may be applied to the periodograms to obtain nonnormalized spectral density estimates
over the set of frequencies
ω = i, i = 1, , n
where nf = NF. These frequencies are in the scale of radians per unit time. The time sampling interval Δt is assumed to be equal to one. Note that the spectral window or weight sequence used to compute
may differ from that used to compute
The cross-spectral density function is complex-valued in general and may be written as
The cospectrum is estimated by
and the quadrature spectrum is estimated by
where
and k(ω) is the integer such that ωk,0 is closest to ω. The
sequence of m = NWT
weights {wj} for
j =
−⌊m/2⌋,
, (m
− ⌊m/2⌋ − 1) satisfies
Σjwj = 1. These weights are
fixed in the sense that they do not depend on the frequency ω at which to
estimate hXY(ω). Usually, m
is odd with the weights symmetric about the middle weight w0. If m is even,
the weight to the right of the middle is considered w0. The argument WT
may contain relative weights since they are normalized to sum to one in the
actual computations. The nonnormalized cross-spectral density estimate is
computed over the same set of frequencies as the nonnormalized spectral density
estimates. However, the particular weight sequence used to compute
need not correspond to either the weight sequence or spectral window used to compute either
An equivalent representation of hXY(ω) is the polar form defined by
The cross-amplitude spectrum is estimated by
and the phase spectrum is estimated by
Finally, the coherency spectrum is estimated by
The coherence or squared coherency is output.
1. The periodograms of X and Y and cross periodogram between X and Y may be computed via the routine CPFFT. The spectral densities of X and Y may then be estimated using any of the routines SSWD, SWED, SSWP, or SWEP. Thus, different window types and/or weight sequences may be used to estimate the spectral and cross-spectral densities given either the series or their periodograms. Note that use of the modified periodograms and modified cross periodogram ensures that the scales of the spectral and cross-spectral densities and their estimates are equivalent.
2 The time sampling interval, TINT, is assumed to be equal to one. This assumption is appropriate for discrete parameter processes. The adjustment for continuous parameter processes (TINT > 0.0) involves multipication of the frequency vector F by 1/TINT and multiplication of the spectral and cross-spectral density estimates by TINT.
3. To convert the frequency scale from radians per unit time to cycles per unit time, multiply F by 1/(2π).
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 years 1770 through 1869. Application of routine CSWEP to these data produces the following results.
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER LDCPM, LDCSM, LDRDAT, N, NDRDAT, NF, NOBS, &
NPAD, NWT
PARAMETER (LDRDAT=100, NDRDAT=4, NF=10, NOBS=100, &
NWT=7, LDCSM=NF, NPAD=NOBS-1, N=NOBS+NPAD, &
LDCPM=N/2+1)
!
INTEGER I, IPVER, NROW, NVAR
REAL COHERE(NF), COSPEC(NF), CPIMAG(LDCPM), &
CPM(LDCPM,10), CPREAL(LDCPM), CRAMPL(NF), &
CSM(LDCSM,9), F(NF), FLOAT, PHASE(NF), PI, PX(LDCPM), &
PY(LDCPM), QUADRA(NF), RDATA(LDRDAT,NDRDAT), &
SX(NF), SY(NF), WT(NWT), X(NOBS), Y(NOBS)
CHARACTER CLABEL1(5)*24, CLABEL2(6)*16, FMT*8, RLABEL(1)*6, &
TITLE1*32, TITLE2*40
INTRINSIC FLOAT
!
EQUIVALENCE (X(1), RDATA(1,2)), (Y(1), RDATA(1,3))
EQUIVALENCE (PX(1), CPM(1,3)), (PY(1), CPM(1,6))
EQUIVALENCE (CPREAL(1), CPM(1,9)), (CPIMAG(1), CPM(1,10))
!
DATA WT/1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0/
DATA FMT/'(F12.4)'/
DATA CLABEL1/'%/%/ k', '%/%/Frequency', '%/%/Period', &
'Spectral%/Estimate%/of X', 'Spectral%/Estimate%/of Y'/
DATA CLABEL2/'%/ k', '%/Cospectrum', '%/Quadrature', &
'Cross%/Amplitude', '%/Phase', '%/Coherence'/
DATA RLABEL/'NUMBER'/
DATA TITLE1/'Results of the Spectral Analyses'/
DATA TITLE2/'%/Results of the Cross-Spectral Analysis'/
! Initialization
PI = 2.0*ASIN(1.0)
DO 10 I=1, NF
F(I) = PI*FLOAT(I)/FLOAT(NF)
CALL SCOPY (NF, F, 1, CSM(1:,1), 1)
CSM(I,2) = 2.0*FLOAT(NF)/FLOAT(I)
10 CONTINUE
! Robinson data
CALL GDATA (8, RDATA, NROW, NVAR)
! 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)
! Estimate the spectral densities
CALL SWEP (N, PX, F, WT, SX)
CALL SWEP (N, PY, F, WT, SY)
! Estimate the cross-spectral density
CALL CSWEP (N, SX, SY, CPREAL, CPIMAG, F, WT, COSPEC, &
QUADRA, CRAMPL, PHASE, COHERE)
! Print results
!
! Copy results to output matrices
CALL SCOPY (NF, SX, 1, CSM(1:,3), 1)
CALL SCOPY (NF, SY, 1, CSM(1:,4), 1)
CALL SCOPY (NF, COSPEC, 1, CSM(1:,5), 1)
CALL SCOPY (NF, QUADRA, 1, CSM(1:,6), 1)
CALL SCOPY (NF, CRAMPL, 1, CSM(1:,7), 1)
CALL SCOPY (NF, PHASE, 1, CSM(1:,8), 1)
CALL SCOPY (NF, COHERE, 1, CSM(1:,9), 1)
! Call printing routines
CALL WRRRL (TITLE1, CSM, RLABEL, CLABEL1, NF, 4, FMT=FMT)
CALL WRRRL (TITLE2, CSM(1:,5), RLABEL, CLABEL2, NF, 5, FMT=FMT)
!
END
Results of the Spectral
Analyses
Spectral
Spectral
Estimate
Estimate
k
Frequency
Period of
X of
Y
1
0.3142 20.0000
116.9550
1315.8370
2
0.6283 10.0000
1206.6086
1005.1219
3
0.9425
6.6667 84.8369
317.2589
4
1.2566
5.0000 55.2120
270.2111
5
1.5708
4.0000 46.5748
115.6768
6
1.8850 3.3333
12.4050
250.0125
7
2.1991
2.8571
7.0934
82.6773
8
2.5133
2.5000
3.4091
62.3267
9
2.8274
2.2222
5.6828
12.8970
10
3.1416
2.0000
4.0346
17.6441
Results of the Cross-Spectral
Analysis
Cross
k Cospectrum
Quadrature
Amplitude
Phase
Coherence
1
94.0531 -254.0125
270.8659
1.2162
0.4767
2
702.5118
21.9823
702.8557
-0.0313
0.4073
3
70.2379
-31.4431
76.9547
0.4209
0.2200
4
-1.8715
-36.1639
36.2123
1.6225
0.0879
5
36.6366
-18.5925
41.0843
0.4696
0.3133
6
32.7071
-6.6569
33.3776
0.2008
0.3592
7
9.4887
-9.1692
13.1950
0.7683
0.2969
8
9.2534
-0.3000
9.2583
0.0324 0.4034
9
-0.5568
2.9455
2.9977
-1.7576
0.1226
10
1.7640
-1.8321
2.5433
0.8043 0.0909
PHONE: 713.784.3131 FAX:713.781.9260 |