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.

CSM — NF 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