Next: detrend() Up: Time Series Macros Help Previous: crosscov()   Contents


crsspectrum(y, len [,reps]), y a REAL matrix with ncols(y) >= 2, len > 0
  and reps > 0 integers.

Keywords: frequency domain, spectrum analysis
crsspectrum(y, len) computes smoothed periodogram and cross-periodogram
estimates of the spectrum and cross-spectrum of the multivariate time
series in the columns of REAL matrix y.

Smoothing is accomplished by 4 successive convolutions of a "boxcar"
smoother of length len.  That is, the smoothing weights are
rep(1/len,len).  The total number of frequencies involved in each
spectrum estimate is 4*len - 3.  When len = 1 no smoothing is done.

crsspectrum(y, len, reps) does the same, except that reps convolutions
with a boxcar are used, where reps > 0 is an integer scalar.  When reps
is odd, then len must also be odd.

See topic 'bandwidth' for information on the approximate bandwidth and
equivalent degrees of freedom associated with these estimates.

Suppose p = ncols(y) and N = nrows(y).  Then, after Sy <-
crsspectrum(y, len [,reps]), Sy is a N by Q matrix, where Q = p +
p*(p-1)/2 = p*(p+1)/2.

Sy[,j], j = 1,...,p are the estimated spectra of y[,j], in Real form.

Sy[,p+1], Sy[,p+2],...,Sy[,Q] contain the estimated cross spectra of
y[,i] and y[,j], i < j, in Hermitian form.  The order is (i,j) = (1,2),
(1,3), ..., (1,p), (2,3),..., (2,p), ... (p-1,p)

Specifically, the estimated cross spectrum of y[,i] and y[,j] is in
column i*(p - (i + 1)/2) + j.

For example, when x and y are vectors, crsspectrum(hconcat(x,y), len)
returns hconcat(Sxx, Syy, Sxy), where Sxx and Syy are estimated spectra
and Sxy is the estimated cross spectrum.  Sxx and Syy are in Real form
and Sxy is in Hermitian form.

The column means are subtracted before computing the periodograms and
cross-periodograms, but no tapering or other detrending is done.

The spectra and cross spectra are computed at the Nfreq Fourier
frequencies 0, 1/Nfreq, 2/Nfreq, ..., (Nfreq-1)/Nfreq cycles per unit
time, were Nfreq is determined as follows:

  When variable S is defined and is a positive integer, Nfreq = S.  It
  is an error if S has a prime factor > 29.

  Otherwise, Nfreq = smallest integer >= 2*nrows(y) which has no prime
  factors > 29, that is Nfreq = goodfactors(2*nrows(y)).

Most of the time, a better choice for estimating cross spectra is macro
compfa() which includes tapering and detrending and for which the amount
of smoothing is specified as a bandwidth or an equivalent degrees of
freedom.  See compfa() for details.

Gary Oehlert 2003-01-15