CSWED
Estimates the nonnormalized cross‑spectral density of two stationary time series using a weighted cross periodogram 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.
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.
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). |
CSM — NF by 9 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 |
X spectral density estimate at F(k) using the specified relative weights contained in WT. |
|
4 |
Y spectral density estimate at F(k) using the specified relative weights contained in WT. |
|
5 |
Co‑spectrum estimate at F(k) using the specified relative weights contained in WT. |
|
6 |
Quadrature spectrum estimate at F(k) using the specified relative weights contained in WT. |
|
7 |
Cross‑amplitude spectrum estimate at F(k). |
|
8 |
Phase spectrum estimate at F(k). |
|
9 |
Coherence estimate at F(k). |
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 periodogram, cosine and sine transformations of each centered and padded time series, the real and imaginary components of the cross periodogram, and the cross‑spectral density estimate based on a specified weight sequence. |
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.0. For a continuous parameter process, TINT > 0.0. TINT is used to adjust the cross‑spectral density estimate.
Default: TINT = 1.0.
NWT — Number of weights. (Input)
NWT must be greater than or equal to one.
Default: NWT= size (WT,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 CSWED (X, Y, F, WT, CPM, CSM [, …])
Specific: The specific interface names are S_CSWED and D_CSWED.
FORTRAN 77 Interface
Single: CALL CSWED (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, NF, F, TINT, NWT, WT, CPM, LDCPM, CSM, LDCSM)
Double: The double precision name is DCSWED.
Description
Routine CSWED estimates the nonnormalized cross‑spectral density function of two jointly stationary time series using a fixed sequence of weights given a sample of n = NOBS observations {Xt} and {Yt} for t = 1, 2, …, n. Let
for t = 1, …, N represent the centered and padded data where N = NOBS + NPAD,
and
is determined by
Similarly, let
for t = 1, …, N represent the centered and padded data where
and
is determined by
The modified periodogram of
for t = 1, …, N 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 {Yt} for t = 1, …, N 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
(Here, ⌊a⌋ means the greatest integer less than or equal to a). 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
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 the spectral density. 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 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 equivalent 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
ω = ƒi, 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 as
hXY(ω) = cXY (ω) − iqXY(ω)
The cospectrum is estimated by
and the quadrature spectrum is estimated by
Note that the same sequence of weights {wj} used to estimate
is used to estimate
The nonnormalized cross‑spectral density estimate is computed over the same set of frequencies as the nonnormalized spectral density estimates discussed above 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 C2WED/DC2WED. The reference is:
CALL C2WED (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, NF, F, TINT, NWT, WT, CPM, LDCPM, CSM, LDCSM, CWK, COEFWK, WFFTC, CPY)
The additional arguments are as follows:
CWK — Complex work vector of length N. (Output)
COEFWK — Complex work vector of length N. (Output)
WFFTC — Vector of length 4N + 15.
CPY — Vector of length 2N.
2. 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 years 1770 through 1869. Application of routine CSWED 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=NPAD+NOBS, &
LDCPM=N/2+1)
!
INTEGER I, NRCOL, NRROW
REAL CPM(LDCPM,10), CSM(LDCSM,9), F(NF), FLOAT, PI, &
RDATA(LDRDAT,NDRDAT), WT(NWT), X(NOBS), Y(NOBS)
CHARACTER CLABEL1(5)*24, CLABEL2(6)*16, FMT*7, RLABEL(1)*6, &
TITLE1*32, TITLE2*40
INTRINSIC FLOAT
!
EQUIVALENCE (X(1), RDATA(1,2)), (Y(1), RDATA(1,3))
!
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)
10 CONTINUE
! Robinson data
CALL GDATA (8, RDATA, NRROW, NRCOL)
! Center on arithmetic means
! Frequency in radians per unit time
! Time interval for discrete data
! Compute the cross periodogram
CALL CSWED (X, Y, F, WT, CPM, CSM)
! Print results
CALL WRRRL (TITLE1, CSM, RLABEL, CLABEL1, NF, 4, FMT=FMT)
CALL WRRRL (TITLE2, CSM(1:,5), RLABEL, CLABEL2, NF, 5, FMT=FMT)
!
END
Output
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
