Computes the cross periodogram of two stationary time series using a fast Fourier transform.

Routine CPFFT computes the cross periodogram of two jointly stationary time series 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 periodogram of the sample sequence {Xt}, t = 1, …, n computed with the padded sequence

is defined by

where

and

represent the

cosine and sine transforms, respectively, and K is the scale factor

The periodogram of the sample sequence {Yt}, t = 1, …, n computed with the padded sequence

is defined 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 periodograms of both

according to the version specified by the argument IPVER. The computational formula for the cross periodogram is given by

where

and

The real part of the (modified) cross periodogram represents the ’raw’ sample cospectrum and the negative of the imaginary part of the (modified) cross periodogram represents the ‘raw’ sample quadrature spectrum (Priestley 1981, page 695). The relationship between the cross periodogram and its complex conjugate is given by

and may be used to recover the cross periodogram at negative frequencies.

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 CPFFT to these data produces the following results. Note that CPFFT sets CPM (1, 2) to the missing value code via routine AMACH.The printing of CPM (1, 2) depends on the computer.