CSSWD

   more...
Estimates the nonnormalized cross‑spectral density of two stationary time series using a spectral window given the time series data.
Required Arguments
X — Vector of length NOBS containing the first stationary time series. (Input)
Y — Vector of length NOBS containing the second stationary time series. (Input)
F — Vector of length NF containing the frequencies at which to evaluate the cross‑spectral density estimate. (Input)
The units of F correspond to the scale specified by IFSCAL. The elements of F must be in the range (π/TINT, π /TINT), inclusive, for IFSCAL = 0 and
(1/(2 * TINT), 1/(2 * TINT)), inclusive, for IFSCAL = 1.
ISWVER — Option for version of the spectral window. (Input)
ISWVER
Action
1
Modified Bartlett
2
Daniell
3
Tukey‑Hamming
4
Tukey‑Hanning
5
Parzen
6
Bartlett‑Priestley
Refer to the “Algorithm” section for further details.
M — Vector of length NM containing the values of the spectral window parameter M. (Input)
For the Parzen spectral window (ISWVER = 5), all values of the spectral window parameters M must be even.
CPM — (N/2 + 1) by 10 matrix containing a summarization 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).
Note N = NOBS + NPAD.
CSMNF by (NM * 7 + 2) matrix containing a summarization of the cross‑spectral analysis. (Output)
The k‑th element of the j‑th column of CSM is defined as
Col.
Description
1
Frequency, F(k).
2
Period, pk where pk = 2π/F(k) for IFSCAL = 0 and pk = 1/F(k) for IFSCAL = 1. If F(k) = 0, pk is set to missing.
3
spectral density estimate at F(k) using the spectral window parameter M(1).
4
Y spectral density estimate at F(k) using the spectral window parameter M(1).
5
Cospectrum estimate at F(k) using the spectral window parameter M(1).
6
Quadrature spectrum estimate at F(k) using the spectral window parameter M(1).
7
Cross‑amplitude spectrum estimate at F(k) using the spectral window parameter M(1).
8
Phase spectrum estimate at F(k) using the spectral window parameter M(1).
9
Coherence estimate at F(k) using the spectral window parameter M(1).
 
NM * 7 + 2
Coherence estimate at F(k) using the spectral window parameter M(NM).
where k = 1, , NF.
Optional Arguments
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 cross periodogram and cross‑spectral density estimate based on a specified version of a spectral window for a given set of spectral window parameters.
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.
NF — Number of frequencies at which to evaluate the cross‑spectral density estimate. (Input)
Default: NF = size (F,1).
TINT — Time interval at which the series are sampled. (Input)
For a discrete parameter process, usually TINT = 1. For a continuous parameter process, TINT > 0. TINT is used to adjust the cross‑spectral density estimate.
Default: TINT = 1.0.
NM — Number of spectral window parameters M used to compute the cross‑spectral density estimate for a given spectral window version. (Input)
NM must be greater than or equal to one.
Default: NM = size (M,1).
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).
LDCSM — Leading dimension of CSM exactly as specified in the dimension statement of the calling program. (Input)
LDCSM must be greater than or equal to NF.
Default: LDCSM = size (CSM,1).
FORTRAN 90 Interface
Generic: CALL CSSWD (X, Y, F, ISWVER, M, CPM, CSM[])
Specific: The specific interface names are S_CSSWD and D_CSSWD.
FORTRAN 77 Interface
Single: CALL CSSWD (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, NF, F, TINT, ISWVER, NM, M, CPM, LDCPM, CSM, LDCSM)
Double: The double precision name is DCSSWD.
Description
Routine CSSWD estimates the nonnormalized cross‑spectral density function of two jointly stationary time series using a spectral window 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 modified periodogram of
is estimated by
where
and
represent the
cosine and sine transforms, respectively, and K is the scale factor equal to 1/(2πn). The modified periodogram of
is estimated 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
The routine PFFT is used to compute the modified periodograms of both
The computational formula for the cross periodogram is given by
where
and
The routine CPFFT is used to compute the modified cross periodogram between
The nonnormalized spectral density of Xt is estimated by
and the nonnormalized spectral density of Yt is estimated by
where the spectral window Wn(θ) is specified by argument ISWVER. The following spectral windows Wn(θ) are available.
Modified Bartlett
where FM(θ) corresponds to the Fejér kernel of order M.
Daniell
Tukey
where DM(θ) represents the Dirichlet kernel. The Tukey‑Hamming window is obtained when
a = 0.23 and the Tukey‑Hanning window is obtained when a = 0.25.
Parzen
where M is even. If M is odd, then M + 1 is used instead of M in the above formula.
Bartlett-Priestley
The argument NM specifies the number of window parameters M and, hence, corresponds to the number of spectral density estimates to be computed for a given spectral window. Note that the same spectral window Wn(θ) and set of parameters M are used to obtain both
The above spectral density formulas assume the data {Xt} and {Yt} correspond to a realization of a bivariate discrete‑parameter stationary process observed consecutively in time. In this case, the observations are equally spaced in time with interval Δt = TINT equal to one. However, if the data correspond to a realization of a bivariate continuous‑parameter stationary process recorded at equal time intervals, then the spectral density estimates must be adjusted for the effect of aliasing. In general, the estimate of hX(ω) is given by
and the estimate of hY(ω) is given by
The nonnormalized spectral density is estimated over the set of frequencies
ω = fi,      i = 1, …, nƒ
where nƒ = NF. These frequencies are in the scale specified by the argument IFSCAL but are transformed to the scale of radians per unit time for computational purposes. The frequency ω of the desired spectral estimate is assumed to be input in a form already adjusted for the time interval Δt.
The cross‑spectral density function is complex‑valued in general and may be written in the following form:
The cospectrum is estimated by
and the quadrature spectrum is estimated by
Note that the same spectral window Wn(θ) and window parameter M used to derive
are also used to compute
The nonnormalized cross‑spectral density estimate is computed over the same set of frequencies as the nonnormalized spectral density estimates with a similar adjustment for Δt.
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.
Comments
1. Workspace may be explicitly provided, if desired, by use of C2SWD/DC2SWD. The reference is:
CALL C2SWD (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, NF, F, TINT, ISWVER, NM, M, CPM, LDCPM, CSM, LDCSM, CX, COEF, WFFTC, CPY)
The additional arguments are as follows:
CX — Complex work vector of length N. (Output)
COEF — Complex work vector of length N. (Output)
WFFTC — Vector of length 4N + 15.
CPY — Vector of length 2N.
2. The centered and padded time series are defined by
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 normalized cross‑spectral density estimate is obtained by dividing the nonnormalized cross‑spectral density estimate in matrix CSM by the product of the estimated standard deviation of X and the estimated standard deviation of Y.
Example
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 CSSWD to these data produces the following results:
 
USE UMACH_INT
USE GDATA_INT
USE CSSWD_INT
USE WRRRL_INT
 
IMPLICIT NONE
INTEGER LDCPM, LDCSM, LDRDAT, N, NDRDAT, NF, NM, &
NOBS, NPAD
PARAMETER (LDRDAT=100, NDRDAT=4, NF=10, NM=2, &
NOBS=100, LDCSM=NF, NPAD=NOBS-1, N=NOBS+NPAD,&
LDCPM=N/2+1)
!
INTEGER I, ISWVER, J, JPT, M(NM), NOUT, NRCOL, NRROW
REAL ASIN, CPM(LDCPM,10), CSM(LDCSM,NM*7+2), F(NF), FLOAT, &
PI, RDATA(LDRDAT,NDRDAT), TINT, X(NOBS), Y(NOBS)
CHARACTER CLABEL1(3)*9, CLABEL2(6)*16, FMT*7, RLABEL(1)*6, &
TITLE*80
INTRINSIC ASIN, FLOAT
!
EQUIVALENCE (X(1), RDATA(1,2)), (Y(1), RDATA(1,3))
!
DATA FMT/'(F10.4)'/
DATA CLABEL1/' k', 'Frequency', 'Period'/
DATA CLABEL2/'%/ k', '%/Cospectrum', '%/Quadrature', &
'Cross%/Amplitude', '%/Phase', '%/Coherence'/
DATA RLABEL/'NUMBER'/
! Initialization
CALL UMACH (2, NOUT)
PI = 2.0*ASIN(1.0)
DO 10 I=1, NF
F(I) = PI*FLOAT(I)/FLOAT(NF)
10 CONTINUE
! Robinson Data
CALL GDATA (8, RDATA, NRROW, NRCOL)
! Center on arithmetic means
! Frequency in radians per unit time
! Spectral window parameters
M(1) = 10
M(2) = 30
! Time interval for discrete data
! Compute cross-spectral density
! using the Parzen window
ISWVER = 5
CALL CSSWD (X, Y, F, ISWVER, M, CPM, CSM)
! Print results
TITLE = 'Cross-Spectral Analysis Using Parzen Window'
CALL WRRRL (TITLE, CSM, RLABEL, CLABEL1, NF, 2, FMT=FMT)
DO 20 J=1, NM
JPT = 7*(J-1) + 5
TITLE = '%/Results of the Cross-Spectral Analysis With '// &
'Spectral Window Parameter M = '
WRITE (TITLE(77:78),'(I2)') M(J)
CALL WRRRL (TITLE, CSM(1:,JPT:), RLABEL, CLABEL2, NF, 5, FMT=FMT)
20 CONTINUE
!
END
Output
 
Cross-Spectral Analysis Using Parzen Window
k Frequency Period
1 0.3142 20.0000
2 0.6283 10.0000
3 0.9425 6.6667
4 1.2566 5.0000
5 1.5708 4.0000
6 1.8850 3.3333
7 2.1991 2.8571
8 2.5133 2.5000
9 2.8274 2.2222
10 3.1416 2.0000
 
Results of the Cross-Spectral Analysis With Spectral Window Parameter M = 10
Cross
k Cospectrum Quadrature Amplitude Phase Coherence
1 463.5888 -65.9763 468.2600 0.1414 0.2570
2 286.5450 -75.0209 296.2029 0.2561 0.1710
3 150.1073 -57.8263 160.8604 0.3677 0.1438
4 52.9840 -32.3642 62.0866 0.5483 0.0998
5 21.5435 -15.0888 26.3020 0.6110 0.0794
6 21.4228 -9.8188 23.5658 0.4298 0.1716
7 15.7005 -5.3704 16.5936 0.3296 0.2112
8 8.0118 -1.8887 8.2314 0.2315 0.1272
9 2.7682 0.2007 2.7754 -0.0724 0.0446
10 0.5777 0.1008 0.5864 -0.1727 0.0091
 
Results of the Cross-Spectral Analysis With Spectral Window Parameter M = 30
Cross
k Cospectrum Quadrature Amplitude Phase Coherence
1 169.7542 -193.4384 257.3615 0.8505 0.1620
2 452.6187 32.3813 453.7755 -0.0714 0.2213
3 94.5221 -90.8159 131.0800 0.7654 0.2629
4 -0.2096 -6.1127 6.1163 1.6051 0.0019
5 27.4711 -22.1946 35.3166 0.6796 0.2492
6 29.1329 -4.0128 29.4080 0.1369 0.3170
7 11.2058 -9.3403 14.5881 0.6948 0.2594
8 8.0017 0.8813 8.0501 -0.1097 0.1928
9 -0.4199 2.2893 2.3275 -1.7522 0.0468
10 0.5570 -1.0767 1.2123 1.0934 0.0678