Next: complex_fun Up: Time Series Macros Help Previous: burg()   Contents

compfa()

Usage:
compfa(y, edf [, degree:D, alpha:A, nfreq:Nfreq, cross:T]), y a REAL
  vector or matrix, nfreq > 0, D integers, edf >=0 and 0 <= A <= .5



Keywords: frequency domain, spectrum analysis
Sf <- compfa(y, edf, degree:D, alpha:A, nfreq:Nfreq) computes spectrum
estimates which are smoothed periodograms of the detrended and tapered
data in the columns of REAL matrix y.

Sf is a Nfreq by ncols(y) matrix whose columns are the estimated
spectra of the detrended and tapered columns of y in Real form.

Nfreq must be a positive integer with no prime factors > 29.

More generally, y can be an array with dimensions n1, n2, n3, ... in
which case the result is also an array with dimensions Nfreq, n2, n3,
... with result[,i2,i3,...] containing the estimated spectrum of
y[,i2,i3,...].

See below for defaults for D, A and Nfreq.

Sf <- compfa(y, edf, cross:T, degree:D, alpha:A, nfreq:Nfreq), with y a
matrix with p > 1 columns and N rows, returns a matrix containing
estimated spectra (smoothed periodograms) and cross spectra (smoothed
cross periodograms) for the tapered detrended columns of y.

It is an error when ncols(y) == 1 (no cross spectrum is defined).

Sy is a Nfreq 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) so that
Sy[,i*(p-(i+1)/2)+j] is the estimated cross spectrum of y[,i] and y[,j].

For example, when x and y are vectors, compfa(hconcat(x,y), edf)
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.

When edf > 2 the amount of smoothing will be chosen so as to yield
spectrum estimates with approximatly edf equivalent degrees of freedom.
See 'bandwidth'.

When 0 < edf <= .5, edf is interpreted as a bandwidth B in cycles per
delta_t and is translated to a working edf = 2*B*N, where N = dim(y)[1].
delta_t is the time interval between successive rows of y.  See topic
bandwidth.

When edf = 0 or edf = 2, no smoothing of periodograms will be done.

The values of the optional keywords other than 'cross', are REAL
scalars with the following restrictions and default values.
 Keyword  Value       Default Value        Restrictions
 degree     D               0              Integer
 alpha      A               0              0 <= A <= .5.
 nfreq    Nfreq        See below           Positive integer with no
                                           prime factors > 29

Default for Nfreq when nfreq:Nfreq is not an argument:
  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 is the smallest integer >= 2*nrows(y) that has no
  prime factor > 29, that is Nfreq = goodfactors(2*nrows(y)).

When D > 0, the columns of y are detrended before any tapering with a
polynomial of degree D fit by least squares, that is the residuals from
the polynomial are tapered and analyzed.  When D = 0, the sample mean
is subtracted.  When D < 0, no detrending is done, not even subtraction
of the mean.

When A = 0, no tapering is done.

When A > 0, after detrending, each column of y will be multiplied by a
cosine taper.  The taper is chosen so as to modify approximately
A*nrows(y) values on each end of the series, approximately 2*A*nrows(y)
values in all.

For example alpha:.5, say, directs that all the observations, except
the middle one, when N is odd, are tapered.  Similarly, alpha:.1
directs that approximately 10% of the observations at the start and 10%
at the end, or 20% in all, will be tapered.

NOTE: The interpretation of tapering proportion A differs from that
used by some practitioners, for whom the tapering proportion is the
proportion of the entire series modified that is tapered.  To modify
100*P percent of the entire series, use alpha:P/2.

See topic costaper() for an exact definition of the taper used.

Generally compfa() is to be preferred to macros spectrum() and
crsspectrum() because it allows polynomial detrending and cosine
tapering which spectrum() and crsspectrum() lack.  In addition, you
specify the amount of smoothing in terms of bandwidth or EDF instead of
the length of smoothing weights.

compfa() uses macros compza(), costaper(), detrend(), read from file
tser.mac if they haven't yet been defined.

See also 'complex_data'.


Gary Oehlert 2003-01-15