Estimates the nonnormalized cross‑spectral density of two stationary time series using a weighted cross periodogram given the time series data.

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

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

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.

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.