Help file for macro file Tser.mac (C) 1999, 2000, 2001, 2002, 2003 by
Gary W.  Oehlert and Christopher Bingham
Created 950330
Last update 030814 CB

???? Starting marker for list of up to 32 comma/newline separated keys
ARIMA models
Autocorrelation
Autocovariance
Complex numbers
Fourier transforms
Frequency domain
General
Plotting
Spectrum analysis
Time domain
????
Ending marker for keys

Note: There are no topics with names longer than 15 characters so as to
allow 5 column output from 'help()' with width >= 80.

Each topic starts with '====topicName' in column 1, followed by an
unlimited number of lines of information.  It is helpful, but not
required, that topics be in alphabetical order.  Topic names should, if
possible, be no longer than 12 characters, since otherwise they must be
quoted in a help() command.  The convention is used that when referring
to functions, but not macros, by name, '()' is appended.

Each topic name that refers to a function or a pre-defined macro has
'()' appended to it.  This does not affect output but is intended to
useful in semi-automatic generation of the reference manual from the
help file.

Names of topics that are not to be included in the Reference Manual are
followed by '*'.  This does not affect help() output.

Names of topics whose usages are not to be included in the Reference
Manual are followed by '%'.  This does not affect help() output.

Keys separated by commas can follow the topic name (and trailing '()',
'*' or '%') if separated from it by '#'.  These should be selected from
up to 32 keys listed between lines starting with '????' somewhere before
the first help topic.  Case is ignored.

The file should be terminated by a line starting '_E_O_F_'

See file MacAnova.hlp for a more complete description of the format of
a help file.

See topic 'tser_index' for a list of topics.

====arspectrum()#spectrum analysis,frequency domain,autoregression,arima models
%%%%
arspectrum(Y,P [,nfreq:nfreq] [,nospec:T]), Y a REAL vector, integer P >
  0, integer nfreq > 0.
%%%%
@@@@usage#Usage
arspectrum(Y,P) estimates the spectrum of Y considered as a discrete
parameter AR(P) (order P autoregressive) time series.  Y must be a REAL
vector with no MISSING elements and P > 0 an integer.

The value returned is a structure(phi:Phi,var:Var, spectrum:Sy)
   Phi    REAL vector of length P containing estimated AR coefficients
   V      REAL scalar containing the estimated variance of the residuals
          (innovations)
   Sy     REAL vector of length Nfreq (see below) containing the
          estimated spectrum computed at frequencies 0, 1/Nfreq,
          2/Nfreq, ..., (Nfreq-1)/Nfreq cycles per Delta_t, the
          interval between observations.

The default value for Nfreq is determined as follows:
  1.  If variable S is defined and is an integer > 2, Nfreq = S.  It is
  an error if S has a prime factor > 29.

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

arspectrum(Y,P,nfreq:Nfreq) or arspectrum(Y,P,Nfreq), where Nfreq > 0 is
an integer, does the same, computing the spectrum at Nfreq frequencies.
It is an error if Nfreq has a prime factor > 29.

@@@@nospec_keyword#Keyword nospec
arspectrum(Y,P,nospec:T) computes phi and var but not spectrum,
returning structure(phi,var).

@@@@method#Method
Estimated AR coefficients Phi are computed from the first P sample
autocorrelations by solving the Yule-Walker equations.  See
yulewalker().  The variance V = gammahat(0)*prod(1 - pacf^2),
where gammahat(0) = sum((y - ybar)^2)/n and pacf is the vector of
P partial autocorrelations associated with Phi.

@@@@see_also#Cross reference
See also burg(), yulewalker().
@@@@______

====autocor()#time domain,autocorrelation
%%%%
autocor(y [, nlags [, nfreq]] [,full:T] [center:T] [,degree:d]), y a REAL
  vector or matrix, nlags > 0 and nfreq > 0 integers, d an integer
  scalar or vector
%%%%
@@@@usage#Usage
r_y <- autocor(y, nlags) computes sample autocorrelations acf(h,y) (y a
REAL vector) or acf(h,y[,j]) (y a REAL matrix) for lags h = 1, 2, ...,
nlags.

The columns of y, which may have no MISSING elements, are interpreted as
time series defined at equally spaced time points.

For a vector y of length N, the sample autocorrelations are defined as
   acf(h,y) = acvf(h,y)/acvf(0,y),
where
   acvf(h,y) = sum((y[i] - ybar)*(y[i+h] - ybar),i=1,...,N-h)/N
is the sample autocovariance.  ybar = sum(y)/N is the sample mean of y.

When ny = ncols(y) = 1, r_y is a vector of length L = nlags with r_y[k]
= acf(k,y).  When ny > 1, r_y is a L by ny matrix with r_y[k,j] =
acf(k,y[,j]).

autocor(y) is the same as autocor(y, nrows(y) - 1), computing all
non-degenerate sample autocorrelations.

@@@@degree_keyword#Keyword degree
You can specify that columns of y be detrended by subtracting a
polynomial in time using keyword phrase 'degree:d', where d is an
integer scalar or vector of length ny = ncols(y).  A scalar d is
interpreted as rep(d,ny).

r_y <- autocor(y [,nlags] degree:d) computes autocorrelations of the
residuals from polynomial trends fit by least squares to each column of
y.  The values of the fitted polynomial replace ybar in the definition
of acvf().

d[j] is the degree of the polynomial fit to column j.  When d[j] = 0,
only the sample mean is subtracted (default).  When d[j] < 0, nothing is
subtracted.

@@@@keywords#Keywords
You can modify the behavior of autocor() using keyword phrases 'full:T'
(compute acf for negative lags, with acf(-h,y) = acf(h,y) in row
2*nlags+2-h), and 'center:T' (same, but with acf(0,y) in row nlags+1).
Note that the output includes acf(0,y) which is not included in the
default output.  See crosscor() for details.

@@@@computation#Computation
Autocovariances acvf(h,y) are computed using crosscov() which makes use
of discrete Fourier transforms (DFTs) of length goodfactors(N+nlags),
the smallest integer >= N + nlags which has no prime factors > 29.  See
topic 'fourier'.

autocor(y, nlags, M) where M is an integer, does the same, using DFTs of
length M.  It is an error if M < N + nlags or if M has a prime factor >
29.  See topic 'fourier'.

@@@@plotting_autocorrelations#Plotting autocorrelations
You can use macro tsplot() to plot the computed autocorrelations.  Any
of the following might be appropriate to plot the first 50
autocorrelations of the columns of y.

  Cmd> tsplot(autocor(y, 50), 1) # line plot
  Cmd> tsplot(autocor(y, 50), 1, impulse:T) # impulse plot
  Cmd> tsplot(autocor(y, 50), 1, impulse:T, lines:T) # line and impulse
  Cmd> tsplot(autocor(y, 50), 1, char:charVec)#plot lines & symbols
  Cmd> tsplot(autocor(y,50,center:T),-50, impulse:T)

On each of these it would usually be appropriate to include keyword
phrases ymin:-1 and ymax:1 since autocorrelations are between -1 and 1.

@@@@see_also#Cross reference
See also autocov(), crosscov(), crosscor(), tsplot().
@@@@______

====autocov()#time domain,autocovariance
%%%%
autocov(y [, nlags [, nfreq]] [,full:T] [center:T] [,degree:d]), y a REAL
  vector or matrix, nlags > 0 and nfreq > 0 integers, d an integer
  scalar or vector.
%%%%
@@@@usage#Usage
c_y <- autocov(y, nlags) computes sample autocovariances acvf(h,y) (y a
REAL vector) or acvf(h,y[,j]) (y a REAL matrix) for lags h = 0, 1, ...,
nlags.

The columns of y, which may have no MISSING elements, are interpreted as
time series defined at equally spaced time points.

For a vector y of length N, the sample autocovariances are defined as
    acvf(h, y) = sum((y[i] - ybar)*(y[i+h] - ybar),i=1,...,N-h)/N
where ybar = sum(y)/N is the sample mean of y.

When ny = ncols(y) = 1, c_y is a vector of length L = nlags + 1 with
c_y[k] = acvf(k-1,y).  When ny > 1, c_y is a L by ny matrix with
c_y[k,j] = acvf(k-1,y[,j]).

autocov(y) is the same as autocov(y, nrows(y) - 1), computing all
non-degenerate sample autocovariances.

@@@@degree_keyword#Keyword degree
You can specify that columns of y be detrended by subtracting a
polynomial in time using keyword phrase 'degree:d', where d is an
integer scalar or vector of length ny = ncols(y).  A scalar d is
interpreted as rep(d,ny).

c_y <- autocov(y [,nlags] degree:d) computes autocovariances of the
residuals from polynomial trends fit by least squares to each column of
y.  The values of the fitted polynomial replace ybar in the definition
of acvf().

d[j] is the degree of the polynomial fit to column j.  When d[j] = 0,
only the sample mean is subtracted (default).  When d[j] < 0, nothing is
subtracted.

@@@@keywords#Keywords
You can modify the behavior of autocov() using keyword phrases 'full:T'
(compute acvf for negative lags, with acvf(-h,y) in row 2*nlags+2-h), and
'center:T' (same, but with acvf(0,y) in row nlags+1).  See crosscov() for
details.

@@@@computation#Computation
Autocovariances are computed using crosscov() which makes use of
discrete Fourier transforms (DFTs) of length goodfactors(N+nlags), the
smallest integer >= N + nlags which has no prime factors > 29.  See
topic 'fourier'.

autocov(y, nlags, M) where M is an integer, does the same, using DFTs of
length M.  It is an error if M < N + nlags or if M has a prime factor >
29.  See topic 'fourier'.

@@@@plotting_autocovariances#Plotting autocovariances
You can use macro tsplot() to plot the computed autocovariances.  Any of
the following might be appropriate to plot the first 50 autocorrelations
of the columns of y.

  Cmd> tsplot(autocov(y, 50), 0) # line plot
  Cmd> tsplot(autocov(y, 50), 0, impulse:T) # impulse plot
  Cmd> tsplot(autocov(y, 50), 0, impulse:T, lines:T) # impulse and lines
  Cmd> tsplot(autocov(y, 50), 0, char:charVec) # plot lines & symbols

@@@@see_also#Cross reference
See also autocor(), crosscov(), crosscor(), tsplot().
@@@@______

====bandwidth%#spectrum analysis,frequency domain
%%%%
Type help(bandwidth) for information about the approximate bandwidth and
equivlent degrees of freedom of a smoothed periodogram.
%%%%
@@@@introduction#Introduction
Suppose a spectrum estimate is computed by circular convolution (see
MacAnova topic convolve()) of weights w[1], w[2], ..., w[nfreq] with a
periodogram computed from an untapered time series at nfreq equally
spaced frequences 0,1/nfreq, 2/nfreq, ..., (nfreq-1)/nfreq cycles per
delta_t where delta_t is the time interval between observations.

An important characteristic of the estimate is its bandwidth, the
effective frequency range over which appreciable smoothing occurs.  The
greater the bandwidth, the more stable (smaller variance) is the
resulting estimate, but the greater the potential for bias because of
the smoothing.  Conversely, the smaller the bandwidth, the smaller is
the bias, at the cost of increased variance of the estimate.

A common definition of the bandwidth associated with a periodogram
computed at nfreq frequencies and smoothed by w[i] is
          B = {sum(w[i])^2/sum(w[i]^2)}/(nfreq*delta_t).

This is in units of cycles per unit time.

When {w[i]} is a "box car" of length m, that is w[i] = 0 except for m
consecutive values of i for which w[i] is constant, then B =
m/(nfreq*delta_t).  In the common situation when nfreq ~= 2*N, where N
= nrows(y), B = .5*m/(N*delta_t).

@@@@convolution_weights#Convolution weights
A common type of smoothing weights are obtained by convolving
(non-circularly) a box car with itself P times.  This is sometimes
called the P-th convolution power of the box car.

For P = 2, the a graph of the weights is a triangle; for higher P the
graph becomes bell shaped.  For large P the graph approximates a normal
curve.  MacAnova macros spectrum(), crsspectrum() and compfa() all use
the 4-th convolution power of a box car as the smoothing weights.

When the weights are computed as the P-th convolution power of a box
car of length m, sum(w^2)/sum(w)^2 and sum(w)^2/sum(w^2) are as follows.
    P             sum(w^2)/sum(w)^2                  sum(w)^2/sum(w^2)
    1     1/m                                               m
    2     (2*m^3+m)/(3*m^4)                       1.5*m-.75/m+O(m^-3)
    3     (11*m^5+5*m^3+4*m)/(20*m^6)            1.818*m-.826/m+O(m^-3)
    4     (151*m^7+70*m^5+49*m^3+45*m)/(315*m^8) 2.086*m-.967/m+O(m^-3)

Thus, when nfreq ~= 2*N, the default smoothing using the 4th convolution
power of a boxcar smoother of length m, has bandwidth approximately
     B = (2.09*m - .97/m)/(nfreq*delta_t) ~= 1.043*m/(N*delta_t)

@@@@equivalent_degrees_of_freedom
                  Equivalent degrees of freedom (EDF)
One common way to summarize the stability of a spectrum estimator is by
its equivalent degrees of freedom or EDF.  Large EDF means smaller
relative standard deviation and thus greater stability.

For a positive random variable W, EDF[W] = 2*E[W]^2/Var[W] = 2/CV[W]^2,
where CV[W] = SD[W]/E[W] is the coefficient of variation.

When W = k*X^2, where X^2 is distributed as chisq(f) (chi-squared on f
degrees of freedom), EDF[W] = f.

Even when W is not a multiple of chi-squared, its distribution may
often be well approximated by that of (E[W]/EDF[W])*chisq(EDF).

When the spectrum is sufficiently smooth, the approximate EDF of a
smoothed periodogram with bandwidth B at a frequency f between B/2 and
.5/deltat_t - B/2 is approximately 2*B*N*delta_t.  At 0 and .5/deltat_t
the EDF is B*N*deltat_t.  For the default smoother, EDF ~=
(4.17*m-1.93/m)*N.  Expressing m in terms of EDF you have the
approximate relationship m ~= EDF/4.17 + 1.93/EDF

@@@@using_a_taper_or_data_window
                      Using a Taper or Data Window
In most situations it is desirable to "taper" the time series before
computing the periodogram.  If yr[i] = y[i] - muhat[i], where muhat[i]
is an estimate of E[y[i]], the tapered time series is h[i]*yr[i],
where (usually) h[i] tapers smoothly to 0 near i = 1 and i = N.  The
sequence h[1],...,h[N] is called a taper or a data window.

Tapering reduces "leakage" -- bias at a given frequency arising from
variation a distant frequencies.  However, it also reduces the EDF of a
smoothed periodogram by a factor of approximately
             R = sum(h^2)^2/{N*sum(h^4)} <= 1
thus increasing the variance.  The greater the amount of tapering, the
smaller is R.

@@@@edf_for_cosine_taper#EDF for cosine taper
Macro compfa(), which computes estimated spectra and, optionally, cross
spectra, uses a "cosine" taper which tapers approximately A*N
observations on each end, where 0 <= A <= .5.  For this taper, the
factor by which the EDF is reduced is

 R = (1 - 5*A/8)^2/(1 - 93*A/128) = 1 - 0.5234*A + 0.0103*A^2 + O(A^3).

Therefore, for the default smoothing, with 100*A percent tapering on
each end and nfreq ~= 2*N, EDF ~= (1 - 0.5234*A)*1.043*m*N.

@@@@see_also#Cross reference
See costaper() for details on the form of a cosine taper.
@@@@______

====burg()#spectrum analysis,frequency domain,time domain,autoregression
%%%%
burg(Y, P [,degree:d, nfreq:Nfreq]), Y a REAL vector, P an integer > 0,
  Nfreq an integer > 0 (length of DFT used).
%%%%
@@@@usage#Usage
burg(Y,P,degree:D,nfreq:Nfreq) estimates the spectrum of Y considered
as a discrete parameter AR(P) (order P autoregression) time series.

The value returned is a structure with the following components:
   phi       REAL vector of length P containing estimated AR
             coefficients
   var       REAL scalar containing the estimated variance of the res-
             iduals (innovations)
   spectrum  REAL vector of length Nfreq (see below) containing the
             estimated spectrum computed at frequencies 0, 1/Nfreq,
             2/Nfreq, ..., (Nfreq-1)/Nfreq cycles per Delta-t, the
             interval between observations.

Y must be a REAL vector and P an integer > 0.  The keyword phrases are
optional.  See below for the default values of D and Nfreq.

Before estimating the spectrum, Y is detrended by subtracting a degree D
polynomial in time fit by least squares.  D = 0 corresponds to
subtracting the sample mean and D < 0 directs that no detrending is
to be done, not even subtracting a mean.

Nfreq must be an integer >= nrows(Y) + P and must have no prime factors
> 29.

If degree:D is omitted, the default value for D is 0 (subtract the
mean).

@@@@default_number_of_frequencies#Default number of frequencies
If nfreq:Nfreq is omitted, the default value for Nfreq = S if S is a
positive integer variable; it is an error if S has a prime factor >
29.

When such an S does not exist, Nfreq = smallest integer >= N + P that
has no prime factors > 29 otherwise, that is Nfreq = goodfactors(N+P),
where N = nrows(Y).

@@@@algorithm#Algorithm
The estimated AR coefficients are computed using an algorithm due to
Burg which does not involve computing the sample autocorrelations.  The
method is sometimes called the maximum entropy method, although that can
equally well describe any method of estimating a spectrum by fitting an
autoregressive model.

@@@@see_also#Cross reference
See also arspectrum(), detrend(), getmacros().
@@@@______

====compfa()#frequency domain,spectrum analysis,cross spectrum
%%%%
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
%%%%
@@@@usage#Usage
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.

@@@@cross_spectrum_usage#Cross spectrum usage
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.

@@@@amount_of_smoothing#Amount of smoothing
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.

@@@@default_keyword_values#Default keyword values
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)).

@@@@detrending#Detrending
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.

@@@@amount of tapering#Amount of tapering
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.

@@@@comparison_with_spectrum_and_crsspectrum#Comparison with spectrum() and crsspectrum()
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.

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

@@@@see_also#Cross reference
See also 'complex_data'.
@@@@______

====complex_data%#complex numbers,frequency domain
%%%%
Type help(complex_data) for information on how MacAnova represents
  complex data in REAL vectors and matrices
%%%%
@@@@representation_of_complex_data
                MacAnova representation of complex data
An unrestricted complex series (vector) of length N is
represented as a REAL N by 2 matrix with the real and imaginary parts in
columns 1 and 2.

An unrestricted N by p complex matrix (p unrestricted complex series)
is represented by a REAL N by 2*p matrix, with the real and imaginary
parts of the j-th series in columns 2*j-1 and 2*j.  A N by 2*p-1 matrix
can be considered to represent p-1 complex series and 1 real series in
the last column.

When the series represents a frequency function evaluated at Nfreq
frequencies, row k usually contains the values for frequency (k-1)/Nfreq
cycles.

@@@@hermitian_form#Hermitian form
A complex vector with Hermitian symmetry (a Hermitian series) {x(j)} of
length N is represented as a REAL vector, say, hx, of length N with
  hx[1] = x(0)
  hx[j+1] = Re(x(j)), 1 <= j < N/2
  hx[N+1-j] = Im(x(j), 1 <= j < N/2.
  hx[N/2+1] = x(N/2) (only when N is even).
p Hermitian series are represented by the columns of a N by p matrix.

When the series represents a frequency function evaluated at Nfreq
frequencies, row 1 contains the value for frequency 0 and rows j and
Nfreq+1-j contain the real and imaginary parts of the function at
frequency (j-1)/Nfreq cycles.

@@@@see_also#Cross reference
See also topic 'hermitian'.
@@@@______

====complex_fun%#complex numbers,frequency domain,fourier transforms
%%%%
Type help(complex_fun) for information on representing complex data in
  MacAnova and on functions and macros for working with complex data.
%%%%
Here is a brief summary of MacAnova functions and macros for working
with complex data.

@@@@notational_conventions#Notational conventions
  rx is a REAL vector or matrix whose columns represent real series
  hx is a REAL vector or matrix whose columns represent complex series
    with Hermitian symmetry (see topics 'hermitian' and 'complex_data')
  cx is a REAL matrix with successive pairs of columns representing
    unrestricted complex series.  If ncols(cx) is odd, the last column
    represents a real series.

@@@@fourier_transforms#Fourier transforms
rft(rx) returns in Hermitian form the DFT of real series in the columns
of rx.

hft(hx) returns as real series the DFT of Hermitian complex series in
hx.

cft(cx) returns as unrestricted complex series the DFT of unrestricted
complex series in cx.

@@@@complex_conjugation#Complex conjugation
hconj(hx) and cconj(cx) compute the complex conjugate of Hermitian
series or unrestricted complex series in hx or cx.

@@@@real_and_imaginary_parts#Real and imaginary parts
hreal(hx) and creal(cx) return the complete real parts of Hermitian
series or unrestricted complex series in hx or cx.

himag(hx) and cimag(cx) return the complete imaginary parts of Hermitian
series or unrestricted complex series in hx or cx.

The results of hreal(), creal(), himag() and cimag() are all
unretricted Real series.

cmplx(rx1,rx2) returns the unrestricted complex form of the series
rx1+i*rx2, where i = sqrt(-1).

@@@@conversion_to_polar_form#Conversion to polar form
hpolar(hx) and cpolar(cx) return the polar forms of Hermitian series or
unrestricted complex series in hx or cx.  The amplitudes (moduli) are
returned in the real parts and the phases (arguments) are returned in
the imaginary parts.

The output of hpolar() is in Hermitian form.  The phases are in radians,
cycles, or degrees depending on the value of option 'angles'.  See
regular topic setoptions() for details.  By default, the phases are
"unwound" to minimize discontinuities if they wrap around the circle.

@@@@conversion_from_polar_form#Conversion from polar form
hrect(hx) and crect(cx) return the usual representation in terms of real
and imaginary parts of the polar forms of Hermitian series or
unrestricted complex series in hx or cx.  The phases of the polar form
are assumed to be in radians, cycles or degrees depending on the value
of option 'angles'.

@@@@converting_to_and_from_hermitian_form#Converting to and from hermitian form
htoc(hx) returns the unrestricted complex form of the Hermitian series
in hx.

ctoh(cx) returns the Hermitian symmetrization of the unrestricted
complex series in cx.  If cx has Hermitian symmetry, then htoc(ctoh(cx))
returns cx.

@@@@complex_multiplication#Complex multiplication
hprdh(hx1, hx2) returns the Hermitian form of the elementwise complex
product of the Hermitian series in hx1 and hx2.  hprdh(hx) returns
hprdh(hx,hx).

hprdhj(hx1,hx2) and hprdhj(hx) are equivalent to hprdh(hx1,hconj(hx2))
and hprdh(hx,hconj(hx)), respectively.

cprdc(cx1, cx2) returns the unrestricted complex form of the elementwise
complex product of the complex series in cx1 and cx2.  cprdc(cx) returns
cprdc(cx,cx).

cprdcj(cx1,cx2) and cprdcj(cx) are equivalent to cprdc(cx1,cconj(cx2))
and cprdc(cx,cconj(cx)), respectively.

@@@@complex_division#Complex division
hdivh(hx1, hx2) returns the Hermitian form of the elementwise complex
ratio of the Hermitian series in hx1 and hx2.  hdivh(hx) returns
hdivh(hx,hx).  When hx is a vector and no represented complex elements
are 0, cdivh(hx,hx) is the Hermitian form of rep(1,nrows(hx)).

hdivhj(hx1,hx2) and hdivhj(hx) are equivalent to hdivh(hx1,hconj(hx2))
and hdivh(hx,hconj(hx)), respectively.

cdivc(cx1, cx2) returns the unrestricted complex form of the elementwise
complex product of the complex series in cx1 and cx2.  cdivc(cx) returns
cdivc(cx,cx).  When cx has 1 or 2 columns, with no zero rows,
cdivc(cx,cx) is the fully complex form of rep(1,nrows(cx)).

cdivcj(cx1,cx2) and cdivcj(cx) are equivalent to cdivc(cx1,cconj(cx2))
and cdivc(cx,cconj(cx)), respectively.

@@@@macros_for_complex_matrices#Macros for complex matrices
In the following a and b are REAL vectors or matrices representing
complex matrices A and B in unrestricted complex form.  All the macros
return their result in unrestricted complex form.

cmatmultc(a,b [,op:"%*%"), cmatmultc(a,b,op:"%c%") and cmatmultc(a,b,
op:"%C") return matrix products of A and B.

ctrace(a) returns the trace of A.

cdiag(a) returns the complex diagonal of square A.

ctranspose(a) returns the transpose of A.

cjtranspose(a) is equivalent to ctranspose(cconj(a)).

csubscr(A), csubscr(A,i), csubscr(A,i,j) simulate A[i], A{i,], and
A[i,j].  i and or j can be empty.

csolve(a) returns the inverse of non-singular square A.

ceigen(a) returns a structure containing the real eigenvalues and complex
eigenvectors of a Hermitian matrix A (ctranspose(a) = cconj(a)).

@@@@______

          Additional functions useful in time series analysis

@@@@autoreg#autoreg()
autoreg(phi,x) applies an autoregressive operator specified by REAL
vector phi of length p to the columns of REAL matrix x.  Keywords
'reverse', 'limits' and 'start' allow for reversing the direction,
operating on a subset of rows and providing starting values.

Uses for autoreg() include generation of autoregressive (AR(p)) time
series (x random), computing residuals from a moving average (MA(p)
time series with coefficients phi, computing cumulative sums (phi = 1)
and solving difference equations.

The sign convention used corresponds to ARSIGN = -1.  See arimahelp
topic 'MASIGN'.

@@@@convolve#convolve()
convolve(wts,x) returns the circular convolutions of REAL vector wts and
the columns of REAL matrix x.

convolve(wts,x,reverse:T) returns the sums of circularly lagged products
of wts with the columns of x.  With keyword phrase 'decimate:n' as an
argument it returns only rows 1, n+1, 2*n+1, ... of the convolutions.

@@@@movavg#movavg()
movavg(theta,x) applies an moving average operator specified by REAL
vector theta of length q to the columns of REAL matrix x.  Keywords
'reverse', 'limits' and 'start' allow for reversing the direction,
operating on a subset of rows and providing starting values.

Uses for movavg() include generation of moving average (MA(q)) time
series (x random), computing residuals from a autoregressive (AR(q))
time series with coefficients theta, computing first differences (theta
= 1), and computing non-circular convolution.

The sign convention used corresponds to MASIGN = -1.  See arimahelp
topic 'MASIGN'.

@@@@padto#padto()
padto(x,n) returns matrix x padded with n additional rows of zeros.  If
n < nrows(x), the result consists of the first n rows of x.

padto() is useful in computing Fourier transforms at Nfreq frequencies
where Nfreq > nrows(x), for example, rft(padto(x,2*nrows)) computes the
Fourier transforms of the columns of x at Nfreq = 2*nrows equally spaced
frequencies.

@@@@partacf#partacf()
phi_kk <- partacf(rho) computes the partial autocorrelations
corresponding to autocorrelations in the column or columns of REAL
vector or matrix rho.

rho <- partacf(phi_kk,reverse:T) computes the autocorrelations
corresponding to partial autocorrelations in the column or columns of
REAL or matrix vector phi_kk.

@@@@polyroot#polyroot()
polyroot(coefs) returns a REAL matrix with 2*ncols(coefs) columns
containing the complex zeros of real polynomials whose coefficients are
specified by the columns of REAL vector or matrix coefs.

polyroot() is useful in determining whether a MA operator is invertible
or a AR operator is stationary.

REAL vector theta defines an invertible MA operator if and only if
   max(creal(cpolar(polyroot(theta)))) < 1.

REAL vector phi defines a stationary AR operator if and only if
   max(creal(cpolar(polyroot(phi)))) < 1.

The sign convention used corrsponds to MASIGN = -1 and ARSIGN = -1 and
is not affected by options 'masign' and 'arsign.  See help topic
'MASIGN'.

@@@@reverse#reverse()
reverse(x) returns a vector or matrix whose rows are the same as those
of REAL or LOGICAL vector or matrix x but in reverse order.

@@@@unwind#unwind()
unwind(theta) returns a vector or matrix derived from the columns of
REAL vector or matrix theta interpreted as angles in units of radians,
cycles or degrees depending on option 'angles'.  Multiples of 2*PI
radians, 1 cycle or 360 degrees are added or subtracted to remove large
row to row jumps.  You can specify the size of jumps using keyword
'crit'.

@@@@yulewalker#yulewalker()
phi <- yulewalker(rho) computes a vector or matrix whose columns are the
autoregression coefficients that solve the Yule-Walker equations based
on autocorrelations in the corresponding column of REAL vector or matrix
rho.

rho <- yulewalker(phi,reverse:T) computes a vector or matrix
autocorrelations that satisfies phi = yulewalker(rho).
@@@@______

                           Informational Topics
@@@@see_also#Cross reference
See topic 'hermitian' for the definition of Hermitian symmetry.

See topic 'fourier' for information on Fourier transforms.

See topic 'complex_data; for information on representing complex data.

See regular MacAnova help topics for full information on cft(), hft(),
rft(), creal(), hreal(), cimag(), himag(), cmplx(), cpolar(), hpolar(),
crect(), hrect(), ctoh(), htoc(), cprdc(), cprdcj(), cdivc(), cdivcj(),
hprdh(), hprdhj(), hdivh(), hdivhj(), autoreg(), movavg(), convolve(),
padto(), partacf(), polyroot(), reverse(), unwind() and yulewalker().
@@@@______

====compza()#frequency domain,spectrum analysis
%%%%
compza(y [,degree:D, alpha:A, nfreq:Nfreq]), y a REAL vector or matrix,
  D and Nfreq != 0 integers, 0 <= A <= .5
%%%%
@@@@usage#Usage
compza(y,nfreq:Nfreq,degree:D,alpha:A) computes scaled Fourier
transforms of the columns of REAL matrix y after detrending and
tapering them.  They are scaled by dividing by sqrt(Ka) where Ka =
sum(taper^2).

The number of frequencies at which the Fourier transform is computed is
abs(Nfreq).  It is an error if Nfreq > 0 and has a prime factor > 29.

Integer D = degree of polynomial used to detrend the columns of y.  D <
0 means no detrending, not even subtracting a mean.

Scalar A >= 0 is the proportion of tapering on each end of the time
series.

See below for defaults for Nfreq, A and Nfreq.

The value is structure(za:Za, n:nrows(y), ka:Ka, alpha:A, degree:D)
  Za        The Nfreq by ncols(y) matrix of Fourier transforms of the
            detrended and then tapered columns of y, in Hermitian form
  Ka        sum(W^2), where vector W are the tapering factors

@@@@default_keyword_values#Default keyword values
The values of the optional keyword phrase arguments 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       Integer != 0 with abs(Nfreq)
                                       not having prime factors > 29

When nfreq:Nfreq is not an argument and positive integer variable S
does not exist, Nfreq is the smallest integer >= 2*nrows(y) with no
prime factor > 29, that is Nfreq = goodfactors(2*nrows(y)).

When S does exist and is a positive integer, nFreq = S.  It is an error
if S has a prime factor > 29.

@@@@detrending#Detrending
When D > 0, the columns of y are detrended with a polynomial in time 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.

@@@@tapering#Tapering
When A = 0, no tapering is done.

When A > 0, after detrending, each column of y will be multiplied by a
cosine taper which modifies approximately A*nrows(y) values on each end
of the series, approximately 2*A*nrows(y) values in all.  Thus alpha:.5,
say, means that the all the observations, except the middle one, when N
is odd, are tapered.  Similarly, alpha:.1 directs that 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.

@@@@method#Method
After detrending followed by tapering, enough zero rows are added to the
end to bring the total number of rows to Nfreq and then rft() is used to
compute the Hermitian form of the Fourier transforms of each column.
These are sometimes known as "modified Fourier transforms".  Finally
they are divided by sqrt(Ka), where Ka = sum(taper^2); see above.

@@@@other_macros_used#Other macros used
compza() uses macros detrend() and costaper().  If they have not
previously been loaded, they are read from file tser.mac.

@@@@see_also#Cross reference
See also detrend(), costaper(), 'hermitian', 'complex_fun',
'complex_data', rft(), 'complex'.

@@@@______

====costaper()#tapering,frequency domain,spectrum analysis,time domain
%%%%
costaper(N, alpha), N > 0 an integer, 0 <= alpha <= .5 a REAL scalar.
%%%%
@@@@usage#Usage
costaper(N, alpha) returns a vector h containing a 100*alpha percent
cosine taper (data window) of length N, tapering approximately N*alpha
elements on each end.  N must be a positive integer and alpha a REAL
scalar with 0 <= alpha <= .5.

When alpha = 0, h = rep(1,N), a "taper" that does no tapering.

For 0 < alpha <= .5 and L = ceiling(alpha*N),
   h[j] = .5*(1 - cos(PI*(j-.5)/L) = sin(.5*PI*(j-.5)/L)^2 , 1 <= j <= L
   h[j] = 1, L + 1 <= j <= N - L
   h[j] = h[N-j+1], N - L + 1 <= j <= N

A taper of length N is used to multiply a (usually) detrended time
series of length N before computing its discrete Fourier transform.

@@@@interpretation_of_tapering_proportion#Interpretation of tapering proportion
NOTE: The interpretation of tapering proportion alpha differs from that
used by some practitioners, for whom the tapering proportion is the
proportion of the entire series modified that is tapered.  To compute a
cosine taper to modify 100*P percent of a series of length N, use
costaper(N, P/2).
@@@@______

====crosscor()#autocorrelation,crosscorrelation,autocovariance,crosscovariance,time domain
%%%%
ccf <- crosscor(x [,nlags] [,degree:d] [,auto:T] [,full:T or center:T]\
  [,nfreq:M]), REAL matrix x, integer nlags > 0, integer scalar or vector
  d, integer M > 0
ccf <- crosscor(x, i, j [,nlags] [,degree:d] [,full:T  or  center:T]\
  [,nfreq:M]), integers i > 0, j > 0
%%%%
@@@@usage#Usage
r_yy <- crosscov(y, nlags) computes sample auto- and crosscorrelation
functions ccf(h,y[,i],y[,j]) for lags h = 0, 1, ..., nlags.

The columns of y, a N by ny REAL matrix with no MISSING elements, are
interpreted as time series defined at equally spaced time points.

r_yy will be a nlags+1 by ny by ny REAL array, with r_yy[h+1,i,j] =
ccf(h,y[,i],y[,j]), h = 0, ..., nlags, where
  ccf(h,y[,i],y[,j]) =
    ccvf(h,y[,i],y[,j])/sqrt(ccvf(h,y[,i],y[,i])*ccvf(h,y[,j],y[,j])),
with
   ccvf(h,y[,i],y[,j]) =
       sum((y[l+h,i]-ybar[i])*(y[l,j]-ybar[j]),l=1,...,N-h)/N.
Note the divisor is N, not N-h or N-h-1.  ybar[j] = sum(y[,j])/N is the
sample mean of y[,j].

Note that y[,j] lags h behind y[,i].  Correlations where y[,j] leads h
ahead of y[,i] are in r_yy[h+1,j,i].

r_yy <- crosscor(y) does the same, except nlags = N-1.

In contrast to the output of autocor(), r_yy contains lag 0 correlations.

r_yy <- crosscor(y, i, j [,nlags]) computes ccf(h,i,j) for h =
-nlags,...,0,...,nlags as a vector of length 2*nlags+1.  Without
center:T (see below), results for lags 0, ..., nlags are in rows 1
through nlags+1 and those for lags -1, ..., -nlags are in rows
2*nlags+1, 2*nlags, ..., nlags+2.  If you want ccf(h,i,j) only for h >=
0, use full:F as an argument.

@@@@keywords#Keywords
You can also use crosscov() keyword phrases 'degree:d' (to control
detrending), 'auto:T' (to restrict computations to autocorrelations),
'full:T' (to compute results for negative lags), 'center:T' (to center
lag 0 between negative and positive lags) and 'nfreq:M' (to specify the
length of Fourier transforms used).  See crosscov() for details.

@@@@computation#Computation
crosscor() is implemented as a macro which uses crosscov() with keyword
phrase 'cor:T' to compute auto- and cross-correlations using discrete
Fourier transforms.  See subtopic 'crosscov:"computation" for details
and for description of the use of keyword 'nfreq'.

@@@@see_also#Cross reference
See also crosscov(), autocor(), autocov().
@@@@______

====crosscov()#crosscovariance,crosscorrelation,autocovariance,autocorrelation,time domain
%%%%
ccvf <- crosscov(y [,nlags] [,degree:d] [,auto:T] [,cor:T]\
  [,full:T or center:T] [,nfreq:M]), REAL matrix y, integer nlags > 0,
  integer scalar or vector d, integer M > 0
ccvf <- crosscov(y, i, j [,nlags] [,degree:d] [,cor:T] \
  [,full:T  or  center:T] [,nfreq:M]), integers i > 0, j > 0
%%%%
@@@@usage#Usage
c_yy <- crosscov(y, nlags) computes sample auto- and crosscovariance
functions ccvf(h,y[,i],y[,j]) for lags h = 0, 1, ..., nlags.

The columns of y, a N by ny REAL matrix with no MISSING elements, are
interpreted as time series defined at equally spaced time points.

c_yy will be a nlags+1 by ny by ny REAL array, with c_yy[h+1,i,j] =
ccvf(h,y[,i],y[,j]), h = 0, ..., nlags, where
   ccvf(h,y[,i],y[,j]) =
       sum((y[l+h,i]-ybar[i])*(y[l,j]-ybar[j]),l=1,...,N-h)/N.
Note the divisor is N, not N-h or N-h-1.  ybar[j] = sum(y[,j])/N is the
sample mean of y[,j].

Note that y[,j] lags h behind y[,i].  Covariances where y[,j] leads h
ahead of y[,i] are in c_yy[h+1,j,i].

crosscov(y) is the same as crosscov(y,nrows(y)-1), computing all
non-degenerate sample autocovariances.

c_yy <- crosscov(y, i, j [,nlags]) computes ccvf(k,i,j) for k =
-nlags,...,0,...,nlags as a vector of length 2*nlags+1.  Without
center:T (see below), results for lags 0, ..., nlags are in rows 1
through nlags+1 and those for lags -1, ..., -nlags are in rows
2*nlags+1, 2*nlags, ..., nlags+2.  If you want ccvf(k,i,j) only for k >=
0, use full:F as an argument.

@@@@auto_keyword#Keyword auto
c_yy <- crosscov(y [,nlags], auto:T) computes only autocovariances.  The
result is a matrix with ny columns with c_yy[h+1,j] = ccvf(h,
y[,j],y[,j]), h = 0,...,nlags, j = 1,...,ny.

crosscov(y, i, j [,nlags], auto:T) is legal only when i = j.

@@@@full_and_center_keywords#Keywords full and center
c_yy <- crosscov(y [,nlags] full:T) does the same, except auto- and
cross-covariances are also computed for negative lags.

c_yy will be a 2*nlags+1 by ny by ny array, with results for lags 0
through nlags in c_yy[1,,], c_yy[2,,],..., c_yy[nlags+1,,], and results
for lags -1 through -nlags in c_yy[2*nlags+1,,], c_yy[2*nlags,,],...,
c_yy[nlags+2,,].

c_yy <- crosscov(y [,nlags] center:T) does the same as with full:T,
except the auto- and cross-covariance functions are centered at row
nlags+1.  Specifically, results for lags -nlags, -nlags+1, ..., -1, 0,
1, ..., nlags are in c_yy[1,,], c_yy[2,,], ..., c_yy[2*nlags+1,,]

With auto:T, with either full:T or center:T, c_yy will be a 2*nlags+1 by
ny matrix.

@@@@degree_keyword#Keyword degree
You can specify that columns of y be detrended by subtracting a
polynomial in time using keyword phrase 'degree:d', where d is an
integer scalar or vector of length ny = ncols(y).  A scalar d is
interpreted as rep(d,ny).

c_yy <- crosscov(y [,nlags], degree:d) computes auto- and cross-
covariances of the residuals from polynomial trends fit by least squares
to each column of y.  The values of the fitted polynomial replace ybar
in the definition of ccvf.

d[j] is the degree of the polynomial fit to column j.  When d[j] = 0,
only the sample mean is subtracted (default).  When d[j] < 0, nothing is
subtracted.

You can use keyword phrases 'full:T', 'center:T' and 'auto:T' with
'degree:d'

@@@@correlations#Correlations
r_yy <- crosscov(y [,i,j] [,nlags], cor:T ...) computes auto- and cross-
correlations instead of auto- and cross-covariances.  You can use any of
the other keywords with 'cor:T'.

Note that, unlike autocor(), the result includes lag 0 correlations.

@@@@computational_method#Computational method
The necessary sums of lagged products are computed using discrete
Fourier transforms (DFTs) of length >= N + nlags, the actual length
chosen being M = goodfactors(N + nlags).  See goodfactors().

When working with long time series, it may be possible to find M1 > M
with more small factors resulting in faster DFT computation.  For
example, when N + nlags = 24388, M = goodfactors(24388) = 24389 = 29^3,
but a DFT of length M1 = 24576 = 3*2^13 takes only about 40% as long to
compute as a DFT of length M.

c_yy <- crosscov(y[,i,j] [,nlags], nfreq:M1 ...) uses DFTs of length M1.

@@@@see_also#Cross reference

See also crosscor(), autocov(), autocor().
@@@@______

====crsspectrum()#frequency domain,spectrum analysis,cross spectrum
%%%%
crsspectrum(y, len [,reps]), y a REAL matrix with ncols(y) >= 2, len > 0
  and reps > 0 integers.
%%%%
@@@@usage#Usage
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.

@@@@form_of_result#Form of result
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.

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

@@@@length_of_fourier_transform#Length of fourier transform
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)).

@@@@comparison_with_compfa#Comparison with compfa()
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.
@@@@______

====detrend()#time domain
%%%%
detrend(x [,Degree]), REAL vector or matrix x, integer scalar or vector
  Degree
%%%%
@@@@usage#Usage
y1 <- detrend(y, k, time:tt) detrends each column of REAL vector or
matrix by subtracting a degree k polynomial P_k(tt).  Each column of y,
which can have no MISSING values, is considered as a time series
observed at time points tt[1], tt[2], ..., tt[N], where N = nrows(y).
tt is a REAL vector of length N and with no MISSING values.  y1 has the
same size and shape as y.

k is usually an integer scalar but can be a length ny = ncols(y) vector
of integers, in which case column j is detrended by a degree k[j]
polynomial.

When k = 1, a linear trend is removed.  When k = 0, only the sample mean
is subtracted.  When k < 0, nothing is subtracted and y1 = y.

Each polynomial is fit by least squares.

You can find the trend values themselves by y - detrend(y,i,time:tt).

detrend(y, time:tt) is equivalent to detrend(y, 1, time:tt), subtracting
least squares regression lines a + b*tt from columns of y.

detrend(y [,k]) is equivalent to detrend(y [,k], time:run(nrows(y))).
This is what you should use with discrete parameter time series observed
at equally spaced time points.

Caution: Some early versions of detrend() used k = 0 as the default,
so that detrend(y) removed only the column means.  It was changed
because it seemed inconsistent with the usual meaning of "detrend".

@@@@examples#Examples
Compute residuals from linear trend:

  Cmd> detrended <- detrend(y) # or detrend(y,1,time:run(nrows(y)))

Plot a discrete parameter monthly time series starting January 1995
together with cubic trend:

  Cmd> tsplot(hconcat(y,y-detrend(y,3)),1995,1/12,title:"Time series")

Detrend y[,1] linearly, y[,2] using a cubic polynomial and y[,3] not at
all:

  Cmd> y1 <- detrend(y,vector(1,3,-1))

@@@@see_also#Cross reference
See also regress(), orthopoly(), tsplot().

@@@@______

====dpss()#tapering,spectrum analysis,frequency domain,tapering
%%%%
dpss(N, W, K [,First]), N >= 1, 1 <= K <= N, 1 <= First <= N - K + 1
  integers, 0 < W < .5
%%%%
@@@@usage#Usage
dpss(N, W, K) computes Discrete Prolate Spheroidal Sequences (DPSS) with
half bandwidth W for orders 0, 1, ..., K-1.  The order L DPSS is the
L+1-th orthonormal eigenvector of a certain tridiagonal symmetric matrix
of order N (see below).  The value returned is N by K REAL matrix, with
the first K eigenvectors in columns 1, 2, ... K in order of decreasing
eigenvalue.

DPSS sequences are used as tapers in multi-taper spectrum estimation.
See below.

dpss(N, W, K, J) does the same, except the DPSS have orders J-1, J, ...,
J+K-2, that is they are eigenvectors J, J+1, ..., J+K-1 of the
tridiagonal matrix.  It returns a N by K REAL matrix, with eigenvectors
J, J+1,..., J+K-1 in columns 1, 2, ... K.

N, K and J must be integers satisfying N >= 1, 1 <= K <= N, 1 <= J <= N
- K + 1.

W must be a REAL scalar, 0 < W < 0.5.

@@@@use_as_tapers#Use as tapers
Discrete Prolate Spheroidal Sequences (DPSS) are used as tapers (data
windows) in multi-taper spectrum estamation.  Their continuous Fourier
transforms are very highly concentrated in low frequencies with a very
sharp cutoff near frequencies W and -W cycles.  Because they are
eigenvectors of a symmetric matrix, they are orthogonal.

The diagonal and subdiagonal of the tridiagonal matrix are
    d = cos(2*PI*W)*(.5*run(-N+1,N-1,2))^2
and
    e = (run(N-1)*run(N-1,1))/2

The DPSS are computed by trideigen(d,e,J,J+K-1,values:F), followed by
certain sign changes.  See regular topic trideigen().
@@@@______

====evalpoly()#
%%%%
evalpoly(coef,z), coef and z REAL matrices with ncols(z)=2*ncols(coef)
evalpoly(coef,z,T), coef and z REAL matrices with ncols(z)=ncols(coef)
%%%%
@@@@usage#Usage
evalpoly(coef,z) evaluates polynomials with REAL coefficients in coef
for complex arguments specified by z.

Specifically, when coef is a n by p REAL matrix and z is a N by 2*p
REAL matrix considered as representing a N by p complex matrix Z,
evalpoly(coef,z) evaluates
  Z^n - coefs[1,]*Z^(n-1) - ... - coef[n-1,]*Z - coef[n,],

The result is a N by 2*p REAL matrix representing a N by p complex
matrix.

By definition, evalpoly(coef,polyroot(coef)) should be zero within
rounding error.

It is an error when ncols(z) != 2*ncols(coef).

evalpoly(coef,x,T) does the same except x is considered to be a real
rather than complex matrix with ncols(x) = ncols(coef).  The result
is a matrix with the same dimensions as x and is the same as
creal(evalpoly(coef,cmplx(x)))

@@@@see_also#Cross reference
See also topics 'complex_data', cmplx(), creal() and polyroot().
@@@@______

====ffplot()#frequency domain,plotting,complex numbers
%%%%
ffplot(G [, Range[, delta_t]] [,timeunit:Unit] [,plotting keywords]),
  G a REAL matrix, Range a real scalar or vector of length 2, REAL
  scalar delta_t > 0, CHARACTER scalar Unit
%%%%
@@@@usage#Usage
ffplot(G, Range, delta_t) plots the columns of REAL matrix G considered
as functions of frequency f for values of f specified by Range.

If N = nrows(G), the i-th row of G is associated with frequency (i-1)/N
cycles per delta_t time units, or (i-1)/(N*delta_t) cycles per unit time
so the full range of frequencies in G is assumed between 0 and
((N-1)/N)/delta_t.

When Range is vector(f1, f2), G is plotted for all Fourier frequencies
between f1 and f2, inclusive.

Range = f, where f is a non-zero scalar, is equivalent to Range =
vector(0,f).

Range = 0 is equivalent to Range = vector(0,.5/delta_t).

You can provide a default value for delta_t by setting variable DELTAT
appropriately.  You can also provide a time unit to be used in
constucting the default x-axis label by setting variable TIMEUNIT.  See
below.

You can use the usual graphic keywords, including 'title', 'xlab',
'ylab', 'xmin', 'xmax', 'ymin', 'ymax', and 'linetype'.

@@@@timeunit_keyword#Keyword timeunit
ffplot(G,Range,delta_t,timeunit:Unit), where Unit is a CHARACTER scalar
such as "year", does the same, the default x-axis label will be, say,
"Frequency (cycles/year)".

Unit should be specified consistantly with delta_t.  For example, with
monthly data and delta_t = 1/12, Unit should be "year", while with
delta_t = 1, Unit should be "month".

@@@@defaults_for_arguments#Defaults for arguments
You can omit delta_t or both Range and delta_t, and provide a default
time unit.

When you omit argument delta_t (ffplot(G, Range)), the default for
delta_t is variable DELTAT if it is a positive scalar and is 1
otherwise.

When you omit both arguments Range and delta_t (ffplot(G), the defaults
for Range and delta_t are vector(0,.5/DELTAT) and DELTAT, when DELTAT
is a positive scalar, and vector(0, .5) and 1 otherwise

Without keyword 'timeunit', the default x-axis label will be
constructed from the value of CHARACTER scalar TIMEUNIT, if it exists
and differs from "".

Macro tsplot() also uses variables DELTAT and TIMEUNIT, as well as
variable START, to construct a default title and x-axis label.  You
can set them once and forget about them.  For example
  Cmd> START <- 1991; DELTAT <- 1/12; TIMEUNIT <- "year"

to ensure that the x-axis label for both frequency and time domain
plots will be informative.

@@@@plotting_outside_normal_range#Plotting outside normal range
If some or all of Range is outside of the interval (0,.5/delta_t), the
values plotted are the periodic extension (with period N) of each column
of G.  Thus ffplot(G,vector(-.5,.5), 1) is legal and plots the a full
cycle of each column from frequency -.5 to .5, with the values with
frequencies < 0 coming from rows i with i > N/2.

@@@@hermitian_argument#Hermitian argument
If the columns of G are complex in Hermitian form, ffplot(G, Range) will
produce the same plot as ffplot(hreal(G), Range) as long as the range
specified is contained in the interval (0, .5/DELTAT).  See
'complex_data'.

@@@@see_also#Cross reference
See also tsplot().
@@@@______

====fourier%#frequency domain,fourier transforms,complex numbers
%%%%
Type help("fourier") for information on discrete and continuous Fourier
  transforms
%%%%
@@@@discrete_fourier_transform
                    Discrete Fourier Transform (DFT)
Let {x(t)} = {x(0), x(1), ..., x(N-1)} = {x(t),0<=t<=N-1}} be a finite
complex or real series, that is, a series of N complex or real numbers.
Then the DFT (discrete Fourier transform) of {x(t)} is the finite
complex series of length N {X<N>(k),0<=k<=N-1}, where

   X<N>(k) = sum(x(t)*exp(-i*2*PI*t*k/N),0<=t<=N-1), k = 0,1,...,N-1

In the argument to exp(), i = sqrt(-1) is a pure imaginary number.

For any integer M > 0 we define
     X<M>(k) = sum(x(t)*exp(-i*2*PI*t*k/M),0<=t<=N-1), k = 0,1,...,M-1
When M > N, this can be viewed as the DFT of the series of length M
obtained by "padding" {x(t),0<=t<=N-1} with M - N 0's.

There are two aspects of this notation:
 The use of an upper case X for the Fourier transform of lower case x.
 The use of suffix <M> to specify that f = k/M is used in the complex
 exponential exp(-i*2*PI*t*f).

The alternative notation DFT<M>{x}(k) = X<M>(k) is sometimes useful.

The definition of x<M>(k) works for any integer k < 0 and k >= M.  With
this extension, x<M>(k) is periodic with period M, that is for any k
   x<M>(k+M) = x<M>(k-M) = x<M>(k).

Note that in the definition of the DFT, indices start with 0 rather than
1 as is assumed in MacAnova.  In working with Fourier transforms in
MacAnova, this correspondence must be kept in mind.  The first element
in a series is always considered to be x(0) or X(0), but to access it
you will need to use x[1].

@@@@inversion_formula
                           Inversion formula
If {X<N>(k)},k = 0, 1, ..., N-1 is the DFT of {x(t)},t = 0, 1, ..., N-1,
then the following inversion formula holds:

  x(t) = (1/N)*sum(X<N>(t)*exp(+i*2*PI*t*k/N),0<=k<=N-1)
       = (1/N)*conj(DFT<N>{conj(DFT<N>{x})}(t)), t  = 0, 1, ..., N-1,

where conj(x) is the complex conjugate of the complex number z.

More generally, when M > N

  x(t) = (1/M)*sum(X<M>(t)*exp(+i*2*PI*t*k/M),0<=k<=M-1)
       = (1/M)*conj(DFT<M>{conj(DFT<M>{x})}(t)), t = 0, 1, ..., N-1,

For t = N,...,M-1, the sum is 0.

@@@@periodic_extensions
                          Periodic Extensions
When dealing with the DFT it is sometimes convenient to consider a
finite complex or real series as being a segment of length N from a
periodic infinite complex or real series {x<N>(t),-oo<t<oo} with period
N:
  x<N>(t) = x(t),      t  = 0, 1, ..., N-1
  x<N>(t) = x<N>(t-N), t = N, N+1, ...
  x<N>(t) = x<N>(t+N), t = -1, -2, ...

More generally, we can define x<M>(t) as the periodic extension with
period M of x(t) padded with M-N 0's:
  x<M>(t) = x(t),      t  = 0, 1, ..., N-1
  x<M>(t) = 0,         t = N,...,M-1
  x<M>(t) = x<M>(t-M), t = M, M+1, ...
  x<M>(t) = x<M>(t+M), t = -1, -2, ...

With this extended definition, for any integer t0

  X<M>(k) = sum(x<M>(t)*exp(-i*2*PI*t*k/M),t0 <= t <= t0+M-1)

that is, as a summation over an arbitrary complete period of x<M>(t).

@@@@continuous_fourier_transform
                   Continuous Fourier Transform (CFT)
The continuous Fourier transform of a finite real or complex series
{x(t),0<=t<=N-1} is the continuous function of the real variable f

        X^(f) = CFT{X}(f) = sum(x(t)*exp(-i*2*PI*t*f),0<=t<=N-1)

The argument f is the frequency at which X^(f) is evaluated and is in
units of cycles.

X^(f) is a periodic function of f with period 1, that is
  X^(f) = X^(f+k) = X^(f-k) for any integer k.

You can consider the DFT to be a "sampling" of the CFT, in the sense
that X<M>(k) = X^(k/M).

Because X<M>(k), k = 0, ..., M-1 is the DFT of x(t) "padded" with M - N
0's, by padding {x(t)} with enough additional zeros, you can use the DFT
to compute the CFT at an arbitrarily dense set of frequencies. For
purposes of computing spectra of a series of length N it is usually
desirable to compute Fourier transforms at approximately M = 2*N
frequencies.  Function padto() is useful for adding zero rows to a
vector or matrix.

You can define the continuous Fourier transform of an infinite real or
complex series {x(t), -oo < t < oo} in a similar manner as
  X^(f) = sum(x(t)*exp(-i*2*PI*t*f), -oo < t < oo)
when the infinite sum converges as it always will if only x(t) != 0 for
only a finite set of t's.

In another terminology, {x(t), -oo < t < oo} are the Fourier coeficients
of the periodic function X^(f).
@@@@______

====gettsmacros()#general
%%%%
gettsmacros(name1 [,name2 ... ]), name1, name2 ... unquoted of macro
  names on file TSMACROS or "tser.mac" if CHARACTER scalar TSMACROS does
  not exist
%%%%
@@@@usage#Usage
gettsmacros(Macro1,Macro2,...) retrieves macros Macro1, Macro2, ... from
a file.  The file name is the value of CHARACTER scalar TSMACROS if it
exists or "tser.mac" if it does not.  The macro names must not be
enclosed in quotes or be CHARACTER variables.

gettsmacros(Macro1,Macro2,...,quiet:T) retrieves the macros but
suppresses printing the descriptive comments associated with them.

If there is more than one copy of any of the named macros in the file,
getttmacros retrieves the first one found.

TSMACROS has no predefined value.

Example:  If TSMACROS does not exist or has value "tser.hlp",
  Cmd> gettsmacros(detrend, spectrum, costaper)
retrieves macros detrend(), spectrum() and costaper() from file
"tser.mac".

gettsmacros() should be faster than pre-defined macro getmacros()
because it searches only one file.  See MacAnova help topic getmacros().
@@@@______

====hermitian%#frequency domain,fourier transforms,complex numbers
%%%%
Type help(hermitian) for information on complex series with Hermitian
symmetry
%%%%
@@@@definition#Definition
A complex series {y(j),0<=j<=N-1} of length N that satisfies

         y(0) real and y(n-j) = conj(y(j)), j = 1, ..., N-1

is said to have Hermitian symmetry, or simply to be a Hermitian series.
Here conj(z) denotes the complex conjugate of complex number z.

This should not be confused with the Hermitian symmetry of a square
complex matrix A for which A' = complex conjugate of A.

When N is even, y(N/2) is real for Hermitan y.

If you define the periodic extension as the infinite complex series
  y<N>(j) = y(j),      j = 0, ..., N-1
  y<N>(j) = y<N>(j-N), j = N, N+1, ...
  y<N>(j) = y<N>(j+N), j = -1, -2, ...

then Hermitian symmetry is equivalent to

  y<N>(-j) = conj(y<N>(j)), j = 0, +-1, +-2, ...

If {x(t),0<=t<=N-1} is a real series, then its DFT, {x<N>(t)} is a
Hermitian series.  Conversely if {x(t),0<=t<=N-1} is Hermitian, then
{x<N>(t)} is real.  See topic 'fourier' for information on the DFT.

If {x(t),0<=t<=N-1} is an unrestricted complex series, its Hermitian
symmetrized form is the Hermitian series {y(t),0<=t<=N-1} where

    y(0) = Re(x(0)), y(t) =(1/2)(x(t)+conj(x(N-t))), t = 1, ..., N-1

Using the periodic extensions {x<N>(t)} and {y<N>(t)}, this is
equivalent to

          y<N>(t) = (x<N>(t)+x<N>(-t))/2, j = 0, +-1, +-2, ...

@@@@see_also#Cross reference
See topic 'complex_data' for information on how complex series are
represented in MacAnova.
@@@@______

====multitaper()#spectrum analysis,tapering,frequency domain
%%%%
multitaper(y, W, K [,degree:D,nfreq:Nfreq],deltat:delta_t,wts:wts]),
  REAL vector y, scalar W > 0, integer K >= 1, integer D, integer Nfreq
  >= length(y), REAL scalar delta_t > 0, REAL vector wts with wts[i] > 0
%%%%
@@@@usage#Usage
multitaper(y, W, K, degree:D, deltat:delta_t, nfreq:Nfreq, wts:Wts)
computes multitaper spectrum estimates using K discrete prolate
spheroidal sequences (DPSS) as tapers, with total bandwidth
approximately 2*W cycles per unit time.

The keyword phrases are optional.  See below for defaults.

y is a REAL matrix with no MISSING values whose columns are discrete
parameter time series observed at equally spaced times delta_t units
apart.

delta_t > 0 is a scalar; see below for default.

REAL scalar W specifies the half bandwidth W; 0 < W <.5/delta_t is
required.

Integer K > 0 is the number of DPSS tapers used.  Sensible values for K
should satisfy K < 2*nrows(y)*W*delta_t, but this is not checked.  At a
frequency where the spectrum is smooth, the approximate equivalent
degrees of freedom EDF (see topic 'bandwidth') of an estimate is
2*sum(Wts)^2/sum(Wts^2).  For equal weights, EDF = 2*K.

Wts is a REAL vector of length K with Wts[i] > 0.  Default is Wts =
rep(1/K, K), that is, equal weights.

The estimated spectrum is the weighted average, using weights in Wts,
of K periodograms of detrended columns of y, tapered using DPSS tapers.

@@@@reference#Reference
See Spectral Analysis for Physical Applications by D. B. Percival and
A. T. Walden (Cambridge Univ. Press, 1993) for details on multitaper
spectrum estimation.

@@@@deltat_keyword#Keyword deltat
Keyword phrase deltat:delta_t specifies the assumed interval between
observations.  It affects only the interpretation of W which is
interpreted as cycles per unit time.  Because of this, W must satisfy 0
< W < .5/delta_t.

For example, with hourly data and W in units of cycles per day, you
would use deltat:1/24 and W must satisfy 0 < W < 12.  See below for the
default value of delta_t.

@@@@degree_keyword#Keyword degree
Keyword phrase degree:D specifies that the columns of y are detrended
using a polynomial of degree D fit by least squares.  When D < 0 no
detrending is done.  When D = 0 column means are subtracted.

@@@@number_of_frequencies#Number of frequencies
When nfreq:Nfreq is an argument, estimated spectra will be computed at
Nfreq frequencies.  It is an error if Nfreq < nrows(y) or if Nfreq has
a prime factor > 29.

When nfreq:Nfreq is not an argument and S is a positive integer
variable, Nfreq = S.  It is an error if S has a prime factor > 29.
When such an S does not exist, Nfreq = smallest integer >= 2*nrows(y)
with no prime factor > 29, that is Nfreq = goodfactors(@*nrows(y))

@@@@defaults_for_keywords#Defaults for keywords
The default values of omitted keywords are as follows:
  Keyword    Default value
  degree     D = 0 (subtract only the mean).
  deltat     delta_t = variable DELTAT if it exists; otherwise 1
  nfreq      Nfreq = variable S if it exists or the smallest integer >=
             2*nrows(y) otherwise; it is an error if S has a prime
             factor > 29.
  wts        Wts = rep(1/K,K)

@@@@other_macros_used#Other macros used
multitaper() uses macros dpss() and detrend().  If they are not defined,
an attempt is made to read them from tser.mac using getmacros.
@@@@______

====news*#general
@@@@march_2003#March 2003
030227 tsplot(y,time [keyword phrases]) is a new usage for tsplot(),
where time is a REAL vector with length(time) = nrows(y).

@@@@july_2002#July 2002
020701 There are new macros crosscov() and crosscor() to compute auto-
and cross- covariances and correlations from the columns of a matrix,
possibly after polynomial detrending.

autocor() and autocov() have been rewritten to use crosscov() and have
been enhanced.  You can detrend data and compute autocovariances and
autocorrelations for negative lags.

@@@@january_2001#January 2001
010103 Many small changes in macros, chiefly having to do with the
number of frequencies used in Fourier transforms.

  Keyword 'S' is deprecated; use 'nfreq' instead.  'nfreq:Nfreq' is
  always an error when Nfreq has a prime factor > 29

  When variable S is a positive integer, it is used by several macros as
  the number of frequencies at which a Fourier transform is computed.
  In such a case, it is now always an error if S has a prime factor >
  29.  The error message suggests the next higher integer with no prime
  factor > 29.

010102 detrend() now aborts on tolerance check failure in least squares
fitting of a polynomial.

dpss() does more thorough checking of arguments.

@@@@december_2000#December 2000
001230 autocov() and autocor() now ignore variable S, automatically
picking the number of frequencies unless 'nfreq' is an argument.  In the
latter case, it is an error if nfreq < N + nlags or has a prime factor >
29.

001226 tsplot() and ffplot() recognize keyword phrase 'timeunit:Unit',
where Unit is a CHARACTER scalar to be used in the default xaxis label
and, for tsplot(), in the default title.  The default for Unit is
variable TIMEUNIT if it is defined and is a CHARACTER scalar, and ""
otherwise.

001214 tsplot() updated to handle impulse plotting better.  With
'impulse:T', 'lines:F' is assumed unless 'lines:T' is an argument.
tsplot() also uses variable START as a default for argument start.
Previously the default was always 0.

====spectrum()#spectrum analysis,frequency domain
%%%%
spectrum(y, len [,reps]), y a REAL matrix, and len > 0 and reps > 0
  integers
%%%%
@@@@usage#Usage
Sy <- spectrum(y, len) computes a smoothed periodogram estimate of the
spectrum of each column of REAL vector or 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.

spectrum(y, len, reps) does the same, except that reps convolutions
with a boxcar are used.  If reps is odd, then len must also be odd.

spectrum(y, 1) computes the unsmoothed periodogram.

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

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

@@@@number_of_frequencies#Number of frequencies
The spectrum is 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)).

@@@@comparison_with_compfa#Comparison with compfa()
Most of the time, macro compfa() is a better choice than spectrum() for
estimating spectra.  compfa() tapers the series and allows polynomial
detrending.  Moreover, you specify the amount of smoothing to use by
either a bandwidth or an EDF (equivalent degrees of freedom).  See
compfa() for details.

Other, quite different, choices for estimating spectra are macros
arspectrum() and burg().

@@@@see_also#Cross reference
See also arspectrum(), burg(), crsspectrum(), compfa(), compza().
@@@@______

====testnfreq()#fourier transforms
%%%%
testnfreq(nfreq), nfreq a vector of positive integers
%%%%
@@@@usage#Usage
testnfreq(nfreq) returns a LOGICAL vector the same length as nfreq, a
vector of positive integers.  Element j of the result is True if and
only if nfreq[j] has no prime factors > 29.

testnfreq() is useful in macros which use one of the fast Fourier
transform functions, rft(), dft() and hft(), since these operate only on
vectors or matrices x for which nrows(x) has no prime factors > 29.  It
allows you to test whether this condition is true.

The use of testnfreq() is deprecated since goodfactors() can be used in
an equivalent test.  When N is a scalar, goodfactors(N) == N is True if
and only if testnfreq(N) is true.  Moreover, the value of goodfactors(N)
is the next integer >= N with no prime factors > 29.

@@@@example#Example
  if (!testnfreq(nrows(y))){
    error("nrows(y) has prime factor > 29")
  }else{
    yft <- rft(y)
  }

@@@@see_also#Cross reference
See also topics goodfactors(), primefactors(), rft(), hft(), cft(),
'fourier'

@@@@______

====tser_index*#
%%%%
Type help(tser_index) for a list of topics in file "tser.hlp".
%%%%
File tser.hlp contains help on the following topics:

arspectrum()  Macro to estimate spectrum of autoregression by solving
              the Yule-Walker equations
autocor()     Macro to compute sample autocorrelation function
autocov()     Macro to compute sample autocovariance function
bandwidth     Comments about the bandwidth and EDF of spectrum estimates
burg()        Macro to estimate autoregression coefficients using Burg's
              algorithm and optionally to compute the spectrum of the
              fitted moddel
compfa()      Macro to compute smoothed modified periodograms and,
              optionally, cross periodograms, using cosine tapering,
              with optional detrending
complex_data  Information on representing complex data and series in
              MacAnova
complex_fun   Summary of MacAnova functions for working with complex
              series
compza()      Macro to compute the Fourier transform of a cosine tapered
              series, optionally detrending
costaper()    Macro to compute a cosine taper with a specified amount of
              tapering
crosscor()    Macro to  auto- and crosscorrelation functions
crosscov()    Macro to  auto- and crosscovariance functions
crsspectrum() Macro to compute smoothed periodgrams and cross
              periodogram with no tapering or detrending (deprecated in
              favor of compfa())
detrend()     Macro to remove a polynomial trend in equally spaced time
dpss()        Macro to compute discrete prolate spheroidal sequences
evalpoly()    Macro to evaluate a real polynomial of a complex variable
ffplot()      Macro to plot a frequency function against frequency
fourier       Information concerning Fourier transforms
gettsmacros() Macro to retrieve macros from tser.mac
hermitian     Information on complex series with Hermitian symmetry
multitaper()  Macro to compute multitaper spectrum estimates
news          Items summarizing changes to macros in tser.mac
spectrum()    Macro to compute smoothed periodograms with no tapering or
              detrending (deprecated in favor of compfa())
testnfreq()   Macro to test whether its argument has prime factors > 29
tser_index    This index
tsplot()      Macro to plot time series against time

====tsplot()#plotting,time domain
%%%%
tsplot(y, times [,symbols:Symb] [,lines:F] [,impulse:T]
  [,timeunit:Unit] [,graphics keywords]), y a REAL matrix, times a REAL
  vector with length(time) = nrows(y), Symb a REAL or CHARACTER scalar,
  vector or matrix, Unit a CHARACTER scalar
tsplot(y [,start [,deltat]] [,symbols:Symb] [,lines:F] [,impulse:T]
  [,timeunit:Unit] [,graphics keywords]), Start and delta_t > 0 REAL
  scalars
%%%%
@@@@usage#Usage
tsplot(Y, Times) does a line plot of the columns of REAL matrix y vs
Times, a REAL vector whose non MISSING values must be in increasing
order.

tsplot(Y, Start, Delta_t) does the same using Times = Start +
Delta_t*run(0,nrows(Y)-1).  That is, the plotting times are equally
spaced starting with Start and incrementing by Delta_t.

tsplot(Y, Start) does the same using a default value for Delta_t.

tsplot(Y) does the same using default values for Start and Delta_t.

Defaults for Start and Delta_t are taken from variables START and
DELTAT, if they exist and are non-missing REAL scalars with DELTAT > 0.
When they do not exist the defaults for Start and Delta_t are 0 and 1,
respectively.

You can provide a time unit to be used in constucting a graph title and
x-axis label by setting variable TIMEUNIT or by using keyword
'timeunit'.  See below.

You can use the usual graphic keywords, including 'title', 'xlab',
'ylab', 'xmin', 'xmax', 'ymin', 'ymax', 'linetype' and 'impulse'.

With 'impulse:T', no connecting lines are drawn unless 'lines:T' is also
an argument.

Without 'symbols:Symb', no symbols are plotted unless 'lines:F' is an
argument and 'impulse:T' is not an argument.  See below.

@@@@timeunit_keyword#Keyword timeunit
tsplot(Y, Times, timeunit:Unit) and tsplot(Y, Start, Delta_t,
timeunit:Unit)  do the same except that Unit is used in constructing
the default title and x-axis label.  Unit must be a CHARACTER scalar
such as "year" that is different from "".

Unit should be specified consistantly with Delta_t.  For example, with
monthly data and Delta_t = 1/12, Unit should be "year", while with
delta_t = 1, Unit should be "month".

Without keyword 'timeunit', the default title and x-axis label will be
constructed from the value of CHARACTER scalar TIMEUNIT, if it exists
and differs from "".

@@@@symbols_keyword#Keyword symbols
tsplot(Y, Times, symbols:Symb) and tsplot(Y [,Start ,[Delta_t]],
symbols:Symb) does the same, except that plotting symbols are taken from
REAL or CHARACTER scalar, vector or matrix Symb.  When Symb is a scalar,
it will be used for every point.  When Symb is a vector of length
ncols(Y), Symb[j] will be the plotting symbol for column j.  See
chplot() for further details.

'symbols:?'  is a special case.  It specifies that plotting symbols will
be 1, 2, ..., nrows(y) when Y is a vector and 1, 2, ..., ncols(y) for
each column of Y when ncols(Y) > 1.

@@@@impulse_keyword#Keyword impulse
tsplot(Y, Times, impulse:T, ...) and tsplot(Y [,Start ,[Delta_t]],
impulse:T, ...) do the same, except an "impulse" plot is drawn rather
than a line plot.  If you want both, also include 'lines:T' as an
argument.

@@@@lines_keyword#Keyword lines
tsplot(Y, Times, lines:F, ...) and tsplot(Y [,Start ,[Delta_t]],
lines:F, ...) do the same, except that no lines or impulses are drawn.
If symbols are not supplied, the symbols are the default symbols drawn
by plot().

@@@@graphics_keyword#Keyword graphics
All the usual graphic keywords can be used, including 'impulse',
'lines', 'title', 'xlab', 'ylab', 'xaxis', 'yaxis', 'xmin', 'xmax',
'ymin', 'add', 'linetype', etc.

In particular, 'impulse:T' draws an impulse plot without lines unless
'lines:T'is also an argument.  See regular help topic 'graph_keys' for
details.

@@@@examples#Examples
Suppose the columns of x contain 10 years of monthly data starting in
January, 1948.  Then
  Cmd> tsplot(x, 1948, 1/12, symbols:"\1",xlab:"Year")

will make a plot of the columns of x against time in years, using a
small triangle as plotting symbol.  If DELTAT has value 1/12, argument 3
can be omitted.  If also variable START is 1948, argument 2 can be
omitted.

 Cmd> tsplot(x,1948+run(0,119)/12, symbols:"\1",xlab:"Year")

makes the same plot, ignoring DELTAT and START.

Suppose rhohat contains autocorrelation functions for the columns of x,
starting with lag 1 month, perhaps computed as rhohat <- autocor(x,60).
Then
  Cmd> tsplot(rhohat,1/12,1/12,impulse:T,ymin:-1,ymax:1,\
       xlab:"Lag (Years)", ylab:"Autocorrelation",\
       title:"Autocorrelation functions for x")

makes an impulse plot of the autocorrelation functions.

If you wanted the lags in months, use tsplot(rhohat,1,1,...,xlab:"Lag
(months)",...)  After DELTAT <- 1/12, simply tsplot(rhohat,DELTAT,
impulse:T, ...) would produce the same plot.

@@@@see_also#Cross reference
See also ffplot(), autocor().
@@@@______

_E_O_F_#This should be the last line, an internal End Of File marker
