====info             MACRO
) This file contains macros useful in time series analysis
) Copyright (C) 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003
) by Christopher Bingham
) *********************************************************************
) * Some of these macros use features introduced with Release 3       *
) * of MacAnova 4.11 and will not run with earlier versions.          *
) * You should not install them unless you have a version of MacAnova *
) * dated after December 10, 2000                                     *
) *********************************************************************
) The macros included are
) gettsmacros()   Eases retrieval of macros from tser.mac
) ffplot()        plot a frequency function against frequency
) tsplot()        plot time series against time
) detrend()       remove polynomial trend
) autocor()       compute autocorrelation function
) autocov()       compute autocovariance function
) crosscor()      compute auto- and crosscorrelation functions
) crosscov()      compute auto- and crosscovariance functions
) testnfreq()     test for absence of prime factor > 29
) spectrum()      compute a smoothed periodogram, detrending, no tapering
) crsspectrum()   compute smoothed cross periodogram, detrending, no taper
) costaper()      compute a cosine taper, percentage specified
) compza()        compute a cosine tapered Fourier transform, detrending
) compfa()        compute smoothed modified periodogram, cos taper, detrend
) arspectrum()    compute autoregressive spectrum solving Yule-Walker equations
) burg()          compute autoregressive spectrum using Burg's algorithm
) dpss()          compute discrete prolate spheroidal sequences
) multitaper()    compute multitaper spectrum estimate using dpss()
) evalpoly()      evaluate polynomial with real coefficients
) Version of 020701
))        Changed cat() to vector(), makestr() to structure()
))        Added many uses of delete(var,return:T) to save memory
))        Used keyvalue() to decode keywords
))        Used real:T and logic:T on isxxxx() to check arguments
))        Converted macros to new style format
))        Unavailable macros now all retrieved by getmacros()
)) 990129 added keyword phrase nospec:T to burg()
)) 990315 Added macro autocor()
)) 990324 Some macros use argvalue()
))        Removed "ERROR: " from error messages.
))        Added macro testnfreq()
)) 990331 Updated tsplot() to allow character plotting
)) 001214 Updated tsplot() to better handle impulse plotting; with
))        impulse:T, default for 'lines' is F and vice versa.
))        Also uses variable START as a default starting time
)) 001226 tsplot() and ffplot() recognize timeunit:Unit, with default
))        taken from variable TIMEUNIT, used in labelling
)) 001230 autocov() and autocor() now ignore S, automatically pick
))        number of frequencies
)) 010104 modified handling of Fourier transform lengths
))        It is an error if argument or keyword value Nfreq has a prime
))        factor > 29.  When S is used as default for Nfreq, it is an
))        error if it has a prime factor > 29
)) 020701 Added new macros crosscor() and crosscov() and rewrote
))        autocor(), autocov() and detrend() with enhancement.
))        detrend() now uses orthogonal polynomials and can now use a
))        different degree polynomial for each column
)) 030318 Modified tsplot() to allow usage tsplot(y,time [,...]), where time is
))        a vector with nrows(time) = nrows(y)
%info%

====> gettsmacros <====
gettsmacros      MACRO   DOLLARS
) Macro to retreive time series macros from file
) The file name is the value of TSMACROS if it exists and is a CHARACTER
) scalar or "tser.mac"
) Since it searches only one file it should be faster than getmacros
) usage: gettsmacros(macro1,macro2,... [,quiet:T]).
) Version of 000910; stripped $$
#gettsmacros(macro1,macro2,... [,quiet:T]) retrieves macros from TSMACROS
if($v==0){
	error("must specify at least one macro name")
}
@tsmacros <- if(isscalar(TSMACROS,char:T)){
	TSMACROS
}else{
	"tser.mac"
}
@args <- $A
for(@i,run($v)){
	@macname <- @args[@i]
	<<@macname>> <- if($k==0){
		macroread(@tsmacros,@macname)
	}else{
		macroread(@tsmacros,@macname,$K)
	}
};;
%gettsmacros%

====> autocov <====
autocov         MACRO   DOLLARS
) Macro to compute autocovariance function
) Usage:
)   acvf <- autocov(y [, nlags] [,full:T] [,center:T] [,degree:d])
)   y          REAL vector or matrix of data with no MISSING values
)              considered as ny = ncols(y) time series of length
)              N = nrows(y)
)   nlags      integer > 0, number of lags wanted; when omitted, the default
)              is N - 1. 
)   d          integer scalar or integer vector of length ny;  a scalar d
)              is equivalent to rep(d,ny).  The default is d = rep(0,ny).
)              With d[j] = 0, only the mean of y[,j] will be
)              subtracted before computing sums of lagged products.
)              Otherwise, the least squares degree d[j] polynomial trend
)              will be subtracted from y[,j] before computing sums of 
)              lagged products.  When d[j] < 0, nothing will be subtracted.
)   acvf       Without full:T or center:T, a REAL nlags+1 by ny matrix of
)              auto- and cross-covariances with acvf[k+1,j] = cyy(k,j,j),
)              k=0,...,nlags, where
)                  cyy(k,j,j) = sum(y1[l+k,j]*y1[l,j],l=1,...,N-k)/N,
)              with y1[,j] = detrend(y[,j],d[j]).
)
)              With full:T, a REAL 2*nlags+1 by ny matrix with rows 1 to
)              nlags+1 as before and rows 2*nlags+1, 2*nlags, ..., nlags+2
)              containing cyy(-1,j,j) = cyy(1,j,j), cyy(-2,j,j) = cyy(2,j,j),
)              ..., cyy(-nlags,j,j) = cyy(nlags,j,j).
)
)              With center:T, a REAL 2*nlags+1 by ny matrix with the auto-
)              covariance functions centered with lag 0 in row nlags+1.
)              
) Computation uses Fourier transforms of length >= N + nlags
)
) acvf <- autocov(y, nlags, nfreq ... ) does the same, except the number of
) frequencies used is nfreq, an integer.
) It is an error when nfreq < N + nlags or when nfreq has a prime
) factor > 29.
)
)) Version 020531  Completely new version using macro crosscov()
# usage: acf <- autocov(y [, nlags])
if ($v == 0 || $v > 3) {
	error("usage: acf <- $S(y [,nlags])")
}
if (!ismacro(crosscov)) {
	getmacros(crosscov, silent:T)
}
@center <- keyvalue($K,"center","TF",default:F)
@full <- keyvalue($K,"full","TF",default:@center)
@d <- keyvalue($K,"degree","integer vector",default:0)
@C <- if ($v == 1) {
	crosscov($1, full:@full, center:@center, auto:T, degree:@d)
} elseif ($v == 2) {
	crosscov($1, $2, full:@full, center:@center, auto:T, degree:@d)
} else { #nfreq an argument
	crosscov($1, $2, full:@full, center:@center, auto:T, degree:@d,\
		nfreq:argvalue($3,"nfreq","positive count"))
}
delete(@d)
delete(@C,return:T)
%autocov%

====> autocor <====
autocor         MACRO   DOLLARS
) Macro to compute autocorrelation function
) Usage:
)   acf <- autocor(y [, nlags] [,full:T] [,center:T] [,degree:d])
)   y          REAL vector or matrix of data with no MISSING values
)              considered as ncols(y) time series of length N = nrows(y)
)   nlags      integer > 0, number of lags wanted; when omitted, the default
)              is N - 1. 
)   d          integer scalar or integer vector of length ny;  a scalar d
)              is equivalent to rep(d,ny).  The default is d = rep(0,ny).
)              With d[j] = 0, only the mean of y[,j] will be
)              subtracted before computing sums of lagged products.
)              Otherwise, the least squares degree d[j] polynomial trend
)              will be subtracted from y[,j] before computing sums of 
)              lagged products.  When d[j] < 0, nothing will be subtracted.
)   acf        Without full:T or center:T, a REAL nlags by ny matrix of
)              autocorrelations with acf[k,j] = ryy(k,j,j), k=1,...,nlags,
)              where ryy(k,j,j) = cyy(k,j,j)/cyy(k,0,0) with
)              cyy(k,j,j) = sum(y1[l+k,j]*y1[l,j],l=1,...,N-k)/N,
)              with y1[,j] = detrend(y[,j],d[j]).
)
)              With full:T, a REAL 2*nlags+1 by ny matrix with rows 1 to
)              nlags+1 as before and rows 2*nlags+1, 2*nlags, ..., nlags+2
)              containing cyy(-1,j,j) = cyy(1,j,j), cyy(-2,j,j) = cyy(2,j,j),
)              ..., cyy(-nlags,j,j) = cyy(nlags,j,j).
)
)              With center:T, a REAL 2*nlags+1 by ny matrix with the auto-
)              covariance functions centered with lag 0 in row nlags+1.
)              
) Computation uses Fourier transforms of length at least N + nlags.
)
) acf <- autocor(y, nlags, nfreq) does the same, except the number of
) frequencies used is nfreq, an integer >= N + nlags.  It is an error
) when nfreq < N + nlags or when nfreq has a prime factor > 29.
)
) Macro crosscov() is used to compute the autocorrelation function
)) Version 020531  Completely new version using macro crosscov()
)) Maintained only for backward compatibility
# usage: acf <- autocor(y [, nlags])
if ($v == 0 || $v > 3) {
	error("usage: acf <- $S(y [,nlags])")
}
if (!ismacro(crosscov)) {
	getmacros(crosscov, silent:T)
}
		  
@center <- keyvalue($K,"center","TF",default:F)
@full <- keyvalue($K,"full","TF",default:@center)
@d <- keyvalue($K,"degree","integer vector",default:0)
@R <- if ($v == 1) {
	crosscov($1, full:@full, center:@center, auto:T, cor:T, degree:@d)
} elseif ($v == 2) {
	crosscov($1, $2, full:@full, center:@center, auto:T, cor:T, degree:@d)
} else { #nfreq an argument
	crosscov($1, $2, full:@full, center:@center, auto:T, cor:T, degree:@d,\
		nfreq:argvalue($3,"nfreq","positive count"))
}
delete(@d)
if (!delete(@full,return:T) && !delete(@center,return:T)){ # remove lag 0
	@R <- @R[-1,]
}
delete(@R,return:T)	
%autocor%

====> crosscov <====
crosscov        MACRO   DOLLARS
) Macro to compute auto and cross covariance or correlation functions
) Usage:
)   c_yy <- crosscov(y [,nlags] [,degree:d] [,cor:T])
)   c_yy <- crosscov(y [,nlags] [,degree:d], auto:T [,cor:T])
)   c_yy <- crosscov(y [,nlags] [,degree:d], full:T [,auto:T] [,cor:T])
)   c_yy <- crosscov(y [,nlags] [,degree:d], center:T [,auto:T] [,cor:T])
)   c_yy <- crosscov(y, i, j [,nlags] [,degree:d] [,center:T] [,cor:T])
)   y          REAL vector or matrix of data with no MISSING values con-
)              sidered as ny = ncols(y) time series of length N = nrows(y)
)   i, j       positive integers <= ny
)   nlags      optional positive integer < N specifying the number
)              of lags for which covariances or correlations are to be
)              computed; default is nlags = N - 1
)   d          integer scalar or integer vector of length ny;  a scalar d
)              is equivalent to rep(d,ny).  The default is d = rep(0,ny).
)              With d[j] = 0, only the mean of y[,j] will be
)              subtracted before computing sums of lagged products.
)              Otherwise, the least squares degree d[j] polynomial trend
)              will be subtracted from y[,j] before computing sums of 
)              lagged products.  When d[j] < 0, nothing will be subtracted.
)
)   Without full:T, c_yy will be a nlags+1 by ny by ny array with
)   c_yy[k+1,i,j] = ccvf(k,i,j), k=0,...,nlags, where
)          ccvf(k,i,j) = sum(y1[l+k,i]*y1[l,j],l=1,...,N-k)/N,
)   with y1[,j] = detrend(y[,j],d[j]).  Note that y[,j] lags k behind y[,i].
)   Covariances where y[,j] leads k ahead of y[,i] are in c_yy[k+1,j,i]
)
)   Note the divisor is N, not N-k or N-k-1, in the definition of ccvf(k,i,j)
)
)   With full:T, covariances will be computed for both positive and negative
)   lags with the negative lags (leads) running backward from the last row
)   c_yy will be a 2*nlags+1 by ny by ny array with
)   c_yy[k+1,i,j] = ccvf(k,i,j), k=0,...,nlags, and
)   c_yy[2*nlags+1-k,i,j] = ccvf(-k,i,j) = ccvf(k,j,i), k = 1,...,nlags
)
)   When i and j are arguments, c_yy will be vector(ccvf(0,i,j),ccvf(1,i,j),
)   ..., ccvf(nlags,i,j),ccvf(-nlags,i,j),...,ccvf(-2,i,j), ccvf(-1,i,j)),
)   with length 2*nlags + 1.  With full:F, ccvf(k,i,j) is computed only
)   for k >= 0.
)
)   With center:T, the rows of the result will be "rotated" so that
)   in c_yy, lags run from -nlags in row 1 to + nlags in row 2*nlags+1
)   with lag 0 in row nlags+1.  center:T implies full:T and is illegal
)   with full:F
)
)   With auto:T (illegal on crosscov(y, i, j ...) with i != j), only
)   autocovariances are computed.  The result is a nlags+1 by ny
)   matrix (without full:T or center:T) or a 2*nlags+1 by ny
)   matrix (with full:T or center:T).  crosscov(y,i,i,auto:T) returns
)   the same as crosscov(y,i,i).
)
)   With cor:T, auto- and cross-correlations will be computed instead of
)   auto- and cross-covariances.
)
)   Computation is by means of discrete Fourier transforms (DFTs) of
)   length >= N + nlags
)
)   c_yy <- crosscov(y [,i, j] [,nlags] [keywords], nfreq:nf) does
)   the same, except the length of the DFTs used is nfreq, an
)   integer >= N + nlags with no prime factors > 29.
)
))020531 New macro.  Written by C. Bingham, kb@stat.umn.edu
# usage: ccvf <- $S(y [,i,j] [,nlags] [,full:T or F] [,center:T][, cor:T])
@X <- argvalue($1, "argument 1", "real matrix nonmissing")
@N <- nrows(@X)
@nx <- ncols(@X)

@nlags <- if ($v == 2) {
	argvalue($2, "nlags","positive count")
} elseif ($v == 4) {
	argvalue($4, "nlags","positive count")
} else {
	@N - 1
}
@auto <- keyvalue($K,"auto*","TF",default:F)
@cor <- keyvalue($K,"cor*","TF",default:F)
if ($v > 2) {
	@i <- argvalue($2,"i","positive count")
	@j <- argvalue($3,"j","positive count")
	if (max(@i,@j) > @nx) {
		error("both i and j must be <= ncols(x)")
	}
	@X <- if (@i == @j) {
		@nx <- 1
		vector(@X[,@i])
	} else {
		if (@auto) {
			error("$S(x, i, j [,nlags], auto:T ...) illegal when i != j")
		}
		@nx <- 2
		@X[,vector(@i,@j)]
	}
}
@degree <- keyvalue($K,"degree*","integer vector",default:0)
if(!isscalar(@degree)) {
	if (length(@degree) != @nx) {
		error("value of 'degree' not scalar but of wrong length")
	}
	if (sum(@degree != @degree[1]) == 0) { # all elements the same
		@degree <- @degree[1]
	}
}
@center <- keyvalue($K,"center*","TF",default:F)
@full <- keyvalue($K,"full","TF",default:$v > 2  || @center)
if (@center && !@full) {
	error("'full:F' with 'center:T' illegal")
}
if ($v > 4) {
	error("usage: $S(x [,i,j] [,nlags] [,full:T or F] [,center:T][, cor:T])")
}
@good <- goodfactors(@N+@nlags)
@nfreq <- keyvalue($K, "nfreq*", "positive count",default:@good )
if (@nfreq < @N + @nlags){
	error(paste("Number of frequencies < n + nlags =",@N + @nlags))
}
if (@nfreq < @good){
	error(paste("Number of frequencies",@nfreq,\
			"has prime factor > 29; try", @good))
}
if (!ismacro(detrend)) {
	getmacros(detrend,silent:T)
}

if (isscalar(@degree)) {
	if (@degree >= 0) {
		@X <- detrend(@X,@degree)
	}
} else {
	for(@i, 1, @nx) {
		if (@degree[@i] >= 0) {
			@X[,@i] <- detrend(@X[,@i], @degree[@i])
		}
	}
}
if (@cor) {
	@X <-/ sqrt(sum(@X^2)/@N)
}
@X <- rft(padto(@X,@nfreq))
@J <- if (!delete(@full,return:T)) {
	run(@nlags+1)
} else {
	vector(run(@nlags+1),@nfreq - run(@nlags-1,0))
}
if (@center) {
	@J <- rotate(@J, @nlags)
}
if (@auto && @nx > 1) {
	@ccvf <- hft(hconj(hprdhj(@X)),divbyt:T)[@J,]/@N
} elseif ($v <= 2 && @nx > 1) {
	@ccvf <- array(rep(0,length(@J)*@nx^2),length(@J),@nx,@nx)
	for (@i, 1, @nx) {
		@ccvf[,@i,] <- hft(hconj(hprdhj(@X[,@i],@X)),divbyt:T)[@J,]/@N
	}
} elseif (@nx == 1) {
	@ccvf <- hft(hconj(hprdhj(@X)),divbyt:T)[@J]/@N
} else {
	@ccvf <- hft(hconj(hprdhj(@X[,1],@X[,2])),divbyt:T)[@J]/@N
}
if (@cor && ($v < 3 || @nx == 1)) {
	# make sure lag 0 autocorrelations are exactly 1
	@i <- @center*@nlags + 1
	if (!@auto) {
		for(@j,1,@nx) {
			@ccvf[@i,@j,@j] <- 1
		}
	} else {
		@ccvf[@i,] <- 1
	}
}
delete(@i,@j,silent:T)	
delete(@nfreq,@nlags,@N,@nx,@J, @X, @good,@center)
delete(@ccvf,return:T)
%crosscov%

====> crosscor <====
crosscor    macro  dollars
) Macro to compute auto and cross correlation functions
) Usage:
)   ccf <- crosscor(x [, nlags] [,degree:d])
)   ccf <- crosscor(x [, nlags] [,degree:d], auto:T)
)   ccf <- crosscor(x [, nlags] [,degree:d], full:T [, auto:T])
)   ccf <- crosscor(x [, nlags] [,degree:d], center:T [, auto:T])
)   ccf <- crosscor(x, i, j [, nlags] [,center:T])
)   x          REAL vector or matrix of data with no MISSING values con-
)              sidered as nx = ncols(x) time series of length N = nrows(x)
)   i, j       positive integers <= nx
)   nlags      optional positive integer < N specifying the number
)              of lags for which correlations are to be computed; default
)              is nlags = N - 1
)   d          integer scalar or integer vector of length nx;  a scalar d
)              is equivalent to rep(d,nx).  The default is d = rep(0,nx).
)              With d[j] = 0, only the mean of x[,j] will be
)              subtracted before computing sums of lagged products.
)              Otherwise, the least squares degree d[j] polynomial trend
)              will be subtracted from x[,j] before computing sums of 
)              lagged products.  When d[j] < 0, nothing will be subtracted.
)
)   Without full:T ccvf will be a nlags+1 by nx by nx array with
)   ccf[k+1,i,j] = corxx(k,i,j) = cxx(k,i,j)/sqrt(cxx(k,i,i)*cxx(k,j,j)),
)   k=0,...,nlags, where 
)          cxx(k,i,j) = sum(x1[l+k,i]*x1[l,j],l=1,...,N-k)/N,
)   with x1[,l] = detrend(x[,l],d[l]).  Note that x[,j] lags k behind x[,i].
)   Correlations where x[,j] leads k ahead of x[,i] are in ccf[,j,i]
)
)   With full:T, correlations will be computed for both positive and negative
)   lags with the negative lags (leads) running backward from the last row
)   ccf will be a 2*nlags+1 by nx by nx array with
)   ccf[k+1,i,j] = corrxx(k,i,j), k=0,...,nlags, and
)   ccf[2*nlags+1-k,i,j] = corrxx(-k,i,j) = corrxx(k,j,i), k = 1,...,nlags
)
)   With i and j as arguments, ccf will the 2*nlags+1 length vector
)   vector(corxx(0,i,j),corxx(1,i,j),...,corxx(nlags,i,j),
)   corxx(-nlags,i,j),...,  corxx(-2,i,j), corxx(-1,i,j)).  With full:F,
)   cxx(k,i,j) is computed only for k >= 0.
)
)   With center:T, the rows of ccf will be "rotated" so that
)   lags run from -nlags in row 1 to +nlags in row 2*nlags+1 with
)   lag 0 in row 0.  center:T implies full:T and is illegal with full:F
)
)   With auto:T (illegal on crosscor(x, i, j ...)), only the auto-
)   correlations are computed.  ccf will be a nlags+1 by nx
)   matrix (without full:T or center:T) or a 2*nlags+1 by nx
)   matrix (with full:T or center:T).
)
) The computation uses Fourier transforms of length >= N + nlags
)
)   ccf <- crosscor(x [,i, j] [,nlags] [keywords], nfreq:nf) does
)   the same, except the length of the Fourier transforms used is nfreq,
)   an integer, which must be >= N + nlags and have no prime factors > 29.
)
) crosscor() is equivalent to crosscov() with 'cor:T'
)
))020531 New macro.  Written by C. Bingham, kb@stat.umn.edu
# usage: ccf <- crosscor(x [, nlags] [,degree:d] [,full:T or center:T] [, nfreq:nf])
if (!ismacro(crosscov)) {
	getmacros(crosscov, silent:T)
}
crosscov($0, cor:T)
%crosscor%

====> detrend <====
detrend         MACRO   DOLLARS
) Macro to remove a polynomial trend from the columns of a matrix
) Usage:
)  y <- detrend(x [,degree] [,time:t])
)   x       REAL matrix with no MISSING elements considered as
)           nx = ncols(x) time series of length N = nrows(x)
)   t       REAL vector of length N with no MISSING elements; default
)           is run(N)
)   degree  integer scalar or vector specifying degree(s) of
)           polynomial trend(s) removed. Default is degree = 1.
)           When degree is not a scalar it must be a vector of length
)           nx.  A scalar is equivalent to rep(degree,nx)
)           When degree[j] = 0, the mean of x[,j] will be subtracted
)           from x[,j]
)           When degree[j] = 1 (default) the least squares linear
)           trend in t is subtracted from x[,j]
)           When degree[j] > 1 the least squares degree[j] polynomial
)           in t trend will be subtracted from x[,j]
)           When degree[j] < 0, nothing will be subtracted from x[,j]
)
)   y       REAL matrix the same size and shape as x.
)
)) 020601 New version written by C. Bingham, kb@stat.umn.edu.
))        Argument degree can now be a vector
))        orthopoly() used instead of swp()
))        Values for time epochs can be specified
# usage: $S(y [,degree] [,time:tt]), real matrix y, integer scalar or
#        vector degree; REAL vector tt
if ($k > 1 || $v == 0 || $v > 2){
	error("usage: $S(x [,degree] [, time:t])", macroname:F)
}
@degree <- if($v == 2){
	argvalue($02,"detrending degree","integer vector")
}else{
	1
}

@Y <- argvalue($1,"argument 1","real matrix nonmissing")
@N <- nrows(@Y)
@ny <- ncols(@Y)
if ($k > 0) {
	@t <- keyvalue($K,"time","nonmissing vector")
	if (isnull(@t)){
		error(paste("'",nameof($K),"' is not a legal keyword",sep:""))
	}
} else {
	@t <- run(@N)
}
if (!isscalar(@degree) && length(@degree) != @ny) {
	error("nonscalar degree and length(degree) != ncols(y)")
}
if (length(unique(@degree)) == 1) {
	@degree <- @degree[1]
}

if (nrows(@t) != @N){
	error("number of time values must be nrows(y)")
}
@dims <- dim(@Y)
@maxd <- max(@degree)

if (@maxd >= 0) { # some detrending to do
	@Y <- matrix(@Y)
	@J <- if (isscalar(@degree)){
		run(@ny)
	} else {
		@degree >= 0
	}
	@Y[,@J] <- @Y[,@J] - sum(@Y[,@J])/@N
	delete(@J)

	if (@maxd > 0){
		if (!ismacro(orthopoly)){
			getmacros(orthopoly,silent:T)
		}
		@X <- orthopoly(@t, @maxd, d)[,-1]/sqrt(@N)
		if (isscalar(@degree)) {
			@Y <-- @X %*% (@X %c% @Y)
		} else {
			for(@j,1,@ny) {
				if (@degree[@j] > 0) {
					@XJ <- @X[,run(@degree[@j])]
					@Yj <- @Y[,@j]
					@Y[,@j] <- @Yj - @XJ %*% (@XJ %c% @Yj)
				}
			}
			delete(@j,@XJ,@Yj)
		}
		delete(@X)
	}
}
delete(@degree,@maxd,@N,@ny)
array(delete(@Y,return:T),delete(@dims,return:T))
%detrend%

====> testnfreq <====
testnfreq       MACRO   DOLLARS
) Macro to test for the presence of a factor > 29
) Usage:
)  ok <- testnfreq(nfreq)
)   nfreq         positive integer scalar or vector
)  ok is a LOGICAL vector the same length as nfreq with
)  ok[i] = T if and only if nfreq[i] has no prime factors > 29
)
) Typical usage:
)  if (!testnfreq(nfreq)){error("bad nfreq")}
)) 990324 Created
)) 001214 rewritten using goodfactors()
@nf <- argvalue($1,"nfreq","positive integer vector")
@result <- goodfactors(@nf) == @nf
delete(@nf)
delete(@result,return:T)
%testnfreq%

====> spectrum <====
spectrum        MACRO   DOLLARS
) Macro to do spectrum analysis by smoothing the periodogram of a time
) series
) Usage:
)  Sy <- spectrum(y, len [,reps])
)  y      REAL vector or matrix with no MISSING values.  y is considered
)         to contain ncols(y) time series of length nrows(y).  The result
)         is a vector or matrix with nfreq rows, where nfreq is the number
)         of frequencies at which the spectrum is estimated
)  len    positive integer
)  reps   positive integer, default 4.  reps must be even when len is even.
)
)  Value for Nfreq = number of frequencies used
)   When S is defined and is a positive integer, Nfreq = S.  It is an
)    error when S has a prime factor > 29
)   Otherwise, Nfreq = goodfactor(2*nrows(y)).  No warning is printed
)    when nrows(y) has a prime factor > 29.
)   For both choices, if necessary Nfreq is increased to the next
)    integer with no prime factors > 29.  
)
) The smoothing weights are computed as reps (default 4) convolutions of
) "boxcar" weights of length len.  When len = 1, no smoothing is done,
) so spectrum(y,1) computes the periodogram.
)
) Usually macro compfa is preferable, since it can do tapering,
) detrending, and the amount of smoothing can be specified in terms
) of effective degrees of freedom or bandwidth.
)) 000910 stripped $$
))  uses argvalue(), tests for legal Nfreq.
)) 001230 now always increases Nfreq when necessary, printing warning message.
# usage: spectrum(y,len[,reps])  at S or 2*nrows(y) frequencies if S not defined
if($v < 2 || $v > 3){
	error("usage is $S(y,len[, reps])", macroname:F)
}
@len <- argvalue($2, "argument 2", "positive count")
@reps <- if ($v <= 2){
	4
}else{
	argvalue($3, "argument 3", "positive count")
}
@y <- matrix(argvalue($1, "time series", "nonmissing real matrix"))
@shift <- @reps*(@len-1)/2
if(@shift != floor(@shift)){
	error("Arg 3 to $S must be even when Arg 2 is", macroname:F)
}
@N <- nrows(@y)

if (alltrue(isscalar(S,real:T), S > 0 , S == floor(S))){
	if (S != goodfactors(S)){
		error(paste("S =",@S,"has prime factor > 29;",@nfreq,"try S =",\
			goodfactors(S)))
	}
	@nfreq <- S
}else{
	@nfreq <- goodfactors(2*@N)
}

@boxcar <- rep(1/@len,@len)
@w <- padto(@boxcar,2*@shift+1)
if(@reps > 1){
	for(@i,run(2,@reps)){
		@w <- convolve(@boxcar,@w)
	}
}
@y <- hreal(hprdhj(rft(padto(@y-sum(@y)/@N,@nfreq))))/@N
delete(@N, @nfreq, @boxcar, @reps)
rotate(convolve(delete(@w,return:T),delete(@y,return:T)),\
	-delete(@shift,return:T))
%spectrum%

====> crsspectrum <====
crsspectrum     MACRO   DOLLARS
) Macro to compute spectra and cross spectra of multivariate time series
) by smoothing the periodograms and cross-periodograms.
) Usage:
)  Sy <- crsspectrum(y, len [, rep])
)  y      REAL matrix with no missing values, considered as a multivariate
)         time series with ncols(y) components of length nrows(y)
)  len    positive integer
)  reps   positive integer, default 4.  reps must be even when len is even
)
) The smoothing weights are computed as reps (default 4) repeated
) convolutions of "boxcar" weights of length len.  When len = 1, no
) smoothing is done, so spectrum(y,1) computes the periodograms and
) cross periodograms.
)
) The estimated frequency functions are computed at the nfreq frequencies
) 0, 1/nfreq, 2/nfreq, ..., (nfreq-1)/nfreq, where nfreq is as follows:
)
) 1. When S is not defined, or S is not a positive integer, nfreq is
)    the smallest integer >= 2*nrows(y) that has no prime factors > 29
) 2. When S is a positive integer, nfreq = S;  it is an error if S
)    has a prime factor > 29.
)
) Result:
)  Sy is a nfreq by Q matrix with Q = p + p*(p-1)/2 = p*(p+1)/2 columns,
)  where p = ncols(y).
)  Cols. 1, 2, ..., p are the estimated spectra in real form.
)  Cols. p+1, p+2, ..., p + p*(p-1)/2 are the estimated cross spectra
)  of y[,i] and y[,j] in Hermitian form.  The order is (i,j) = (1,2), (1,3),
)  ..., (1,p), (2,3),..., (2,p), ... (p-1,p).  That is, 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, reps)
)  is hconcat(Sxx, Syy, Sxy), where Sxx and Syy are estimated spectra and
)  Sxy is the estimated cross spectrum.
)) 000910 stripped $$
)) 001230 modified choice of nfreq so that it has no prime factors > 29
))        if S is defined and has prime factors > 29, goodfactors(S)
)         is used, with a warning message
# usage: crsspectrum(y,len [,reps]), at S or about 2*nrows(y) frequencies
if($v < 2 || $v > 3){
	error("usage is $S(y,len[,reps])", macroname:F)
}
@len <- argvalue($2,"argument 2","positive count")
@reps <- if($v <= 2){
	4
}else{
	argvalue($3,"argument 3","positive count")
}
@shift <- @reps*(@len-1)/2
if(@shift != floor(@shift)){
	error("Arg 3 must be even if Arg 2 is")
}
@y <- matrix(argvalue($1, "multivariate time series", "nonmissing real matrix"))
@N <- nrows(@y)
@p <- ncols(@y)

if (alltrue(isscalar(S,real:T), S > 0 , S == floor(S))) {
	if (S != goodfactors(S)){
		error(paste("S =",@S,"has prime factor > 29;",@nfreq,"try S =",\
			goodfactors(S)))
	}
	@nfreq <- S
}else{
	@nfreq <- goodfactors(2*@N)
}


@boxcar <- rep(1/@len,@len)
@w <- padto(@boxcar,2*@shift+1)
if(@reps > 1){
	for(@i,run(2,@reps)){
		@w <- convolve(@boxcar,@w)
	}
}
delete(@reps, @len)
@y <- rft(padto(@y-sum(@y)/@N,@nfreq))/sqrt(@N)
@SP <- padto(rep(0,@p*(@p+1)/2)',@nfreq)
@SP[,run(@p)] <- rotate(convolve(@w,hreal(hprdhj(@y))),-@shift)
if(@p > 1){
	@k <- @p+1
	for(@i,run(@p-1)){
		for(@j,run(@i+1,@p)){
			@SP[,@k] <- ctoh(rotate(convolve(@w,\
			   htoc(hprdhj(@y[,@i],@y[,@j]))),-@shift))
			@k <-+ 1
		}
	}
}
delete(@p,@y,@N,@nfreq,@shift,@i,@boxcar,@w,@k,silent:T)
delete(@SP,return:T)
%crsspectrum%

====> tsplot <====
tsplot        MACRO   DOLLARS
) Macro to make graph of time series against time
) Usage:
)  tsplot(y, time, [,timeunit:Unit] [,symbols:symb] \
)   [,lines:F or T] [,impulse:T] [,labeling keywords])
)  tsplot(y [,start [,deltat]]]  [,timeunit:Unit] [,symbols:symb] \
)   [,lines:F or T] [,impulse:T] [,labeling keywords])
)
)  y         REAL vector or matrix whose columns are time series to
)            be graphed
)  time      REAL vector, nrows(time) = nrows(y)
)  start     nonmissing REAL scalar, default START if defined or 0
)  deltat    positive REAL scalar, default DELTAT if defined or 1
)  Unit      CHARACTER scalar, default TIMEUNIT if defined or "", used
)            in default title and x-axis label
)  symb      REAL or CHARACTER scalar, vector or matrix as expected
)            by chplot()
) In the first usage y[i] or y[i,j] is plotted against time[i]
) i = 1, 2, ..., N = nrows(y)
) 
) In the 2nd usage time[i] is calculated as start + (i-1)*deltat
)
) A line plot will be drawn unless suppressed by lines:F or impulse:T
)
) With impulse:T, an impulse plot will be drawn;  it will have lines
)   only with lines:T.
) 
) To get a pure character plot, use lines:F without impulse:T, with
)   or without symbols:symb.
)
) With symbols:symb, where symb is a REAL or CHARACTER vector or matrix,
)   the graph will be drawn by chplot() with symbols:symb.
)   Among the conventions of chplot() is that when symb is a vector
)   with length(symb) = ncols(y), each symbol is used for a column
)
) symbols:? is a special case.
)   When y is a vector, symbols 1, 2, ..., are used for cases 1, 2, ...
)   When ncols(y) > 1, symbols 1, 2, 3, ... are used for columns
)   1, 2, 3, ... .
)
)) Version 000910 stripped $$
)) 001214 impulse:T without lines:T suppresses lines and vice versa
)) 001226 timeunit:Unit (default CHARACTER scalar TIMEUNIT) used to
))        contruct default xlab and title
)) 030227 New usage tsplot(y, time ...)
# usage: $S(y, time [,symbols:symb] [,lines:F] [,impulse:T] \
#            [,labeling keywords]) 
# or:    $S(y [,start [,deltat]] [,symbols:symb] \
#            [,lines:F] [,impulse:T] [,labeling keywords])
@y <- matrix(argvalue($1,"argument 1","real matrix"))

if($v >= 2){
	@t0 <- argvalue($02, "$2", "real vector")
	if (!isscalar(@t0)) {
		if (length(@t0) != nrows(@y)) {
			error("length of time vector != number of cases")
		}
		if ($v > 2) {
			error("> 2 non-keyword arguments with time vector")
		}
		@J <- !ismissing(@t0)
		if (alltrue(sum(@J) > 1,min(movavg(1,@t0[@J])[-1]) <= 0)) {
			error("Time vector not increasing")
		}
		delete(@J)
	}
} else {
	@t0 <- if (alltrue(isscalar(START,real:T),!ismissing(START))) {
		START
	}else{
		0
	}
}

if (isscalar(@t0)) {
	@dt <- if($v >= 3){
		argvalue($03, "delta_t", "positive real scalar")
	}elseif(isscalar(DELTAT,real:T,positive:T)) {
		DELTAT
	}else{
		1
	}
	@t0 <- vector(@t0,delete(@dt,return:T)) # (start, deltat)
}

@timeunit <- keyvalue($K,"timeunit*","string")
if (isnull(@timeunit)){
	@timeunit <- if (isscalar(TIMEUNIT,char:T)){
		TIMEUNIT
	} else {
		""
	}
}
if (@timeunit != ""){
	@xlab <- paste("Time (",@timeunit,")",sep:"")
	@title <- paste("Plot of $1 vs",@timeunit)
} else {
	@xlab <- "Time"
	@title <- "Plot of $1 vs time"
}

@impulse <- keyvalue($K,"impulse*","TF")
@lines <- keyvalue($K,"lines*","TF")
if (isnull(@impulse) && isnull(@lines)) {
	@impulse <- F
	@lines <- T
} elseif (isnull(@impulse)){
	@impulse <- !@lines && isnull(@symbols)
} elseif (isnull(@lines)){
	@lines <- !@impulse
}

@symbols <- keyvalue($K,"symb*","matrix")
		
if (isnull(@symbols)){
	# when !@impulse and !@lines, lineplot() is like plot()
	lineplot(delete(@t0,return:T),delete(@y,return:T),\
		$K,xlab:@xlab, ylab:"Time series",title:@title,\
		lines:@lines, impulse:@impulse)
}else{
	if (isreal(@symbols)){
		if (isscalar(@symbols) && ismissing(@symbols[1])){
			@symbols <-\
				if(!isvector(@y)){
					run(ncols(@y))
				}else{
					run(nrows(@y))
				}
		}elseif(anymissing(@symbols)){
			error("MISSING values in value for 'symbols'")
		}elseif(min(vector(@symbols)) < 0 ||\
			max(vector(@symbols)) > 999||\
			sum(vector(@symbols!=round(@symbols))) != 0){
			error("REAL value for 'symbols' not integers between 0 and 999")
		}
	}elseif(!ischar(@symbols)){
		error("value for 'symbols' not integers or CHARACTER")
	}
	chplot(delete(@t0,return:T),delete(@y,return:T),\
		symbols:@symbols,$K,xlab:@xlab, ylab:"Time series",	title:@title,\
		lines:@lines, impulse:@impulse)
}
delete(@symbols,@xlab,@title, @lines, @impulse)
%tsplot%

====> ffplot <====
ffplot          MACRO   DOLLARS
) Macro for plotting periodic frequency functions vs frequency in cycles
) per delta-t
) Usage:
)  ffplot(fx [,range [,deltat] ][,graphics keyword phrases][,timeunit:Unit])
)
)  fx       REAL vector or matrix whose columns are interpreted as frequency
)           functions sampled at frequencies j/(N*deltat),j = 0, 1, ..., N-1
)           where N = nrows(fx)
)  deltat   positive scalar, default DELTAT if defined or 1
)  range    REAL scalar or REAL vector of length 2 specifying range of
)           frequencies plotted; default vector(0,.5/deltat)
)  Unit     CHARACTER scalar specifying time unit used in default frequency
)           axis label; default is TIMEUNIT if defined or "" (no time
)           unit in label)
)
)  When range is a vector it is vector(lowFreq,hiFreq), where lowFreq <
)    hiFreq are in units of cycles per deltat.
)  A scalar range = Freq is equivalent to vector(0,Freq) when Freq > 0
)    and to vector(Freq,0) when Freq < 0, Freq in cycles per deltat
)  range = 0 is equivalent to vector(0,.5/deltat), the default range
)) 990324 uses argvalue()
)) 000910 strip $$
)) 001226 timeunit:Unit (default CHARACTER scalar TIMEUNIT) used to
))        contruct default xlab
# usage: ffplot(fx[,range[,deltat]][,labeling keyword phrases])
@f <- matrix(argvalue($1,"argument 1", "real matrix"))
@nfreq <- nrows(@f) # number of frequencies
@f <- matrix(@f,@nfreq)

@dt <- if($v >= 3){
	argvalue($03, "delta_t", "positive number")
}elseif(isscalar(DELTAT,real:T,positive:T)){
	DELTAT
}else{
	1
}

@timeunit <- keyvalue($K,"timeunit*","string")
if (isnull(@timeunit)){
	@timeunit <- if (isscalar(TIMEUNIT,char:T)){
		TIMEUNIT
	} else {
		""
	}
}
if (@timeunit != ""){
	@xlab <- paste("Frequency (Cycles/",@timeunit,")",sep:"")
} else {
	@xlab <- "Frequency (Cycles)"
}

@rng <- if ($v >= 2){
	argvalue($02, "range","nonmissing real vector")
}else{
	0
}
if(length(@rng) > 2){
	error("length(range) not 1 or 2")
}
if(length(@rng) == 1 && @rng[1] == 0){
	@rng <- .5/@dt
}
if(length(@rng) == 1){
	@rng <- vector(0,@rng)
}
@J <- run(vector(floor(@nfreq*@dt*vector(min(@rng),max(@rng))),1))
delete(@rng)
lineplot(@J/(@nfreq*@dt),\
	delete(@f,return:T)[1+delete(@J,return:T) %% @nfreq,],$K,\
	xlab:@xlab, ylab:"Frequency Function")
delete(@nfreq, @dt, @xlab)
%ffplot%

====> arspectrum <====
arspectrum        MACRO   DOLLARS
) Macro to compute autoregresssive spectrum estimate using Yule-Walker
) estimates of autoregression coefficients.
)
) Usage:
)  arspectrum(y,p [,nospec:T])
)  arspectrum(y,p,nfreq:nfreq [,nospec:T])
)  arspectrum(y,p,nfreq [,nospec:T])
)   y           REAL vector with no MISSING elements
)   p > 0       Integer order of AR model fit (number of lags in
)               autoregression)
)   nfreq > 0   Integer number of frequencies at which to estimate
)               spectrum; nfreq must have no prime factors > 29
)
) Default for nfreq
)  When S is defined and is a positive integer, nfreq = S; it is an
)  error if S has a prime factor > 29
)  Otherwise, nfreq = goodfactors(2*nrows(y))
)
) Result:
)  structure(phi:ARCoefs,var:V,spectrum:sy)
)   ARCoefs    REAL vector of fitted coefficients
)   V          REAL scalar estimating the innovation variance
)   sy         REAL vector of length nfreq containing estimated spectrum
) With nospec:T, nfreq is ignored and component spectrum is omitted.
)
)) Version of 990324 added keyword phrase nospec:T; use argvalue()
)) Version 000910; stripped $$
)) Version 020701; added phrase nfreq:L as alternative to non-keywor
))    argument 3
#usage: $S(y,nlags [,nfreq] [,nospec:T]) has components phi, var, spectrum
if($v < 2 || $v > 3){
	error("usage: $S(y,p [,nospec:T]) or  $S(y,p,nfreq [,nospec:T])",\
		macroname:F)
}
@y <- argvalue($1,"argument 1", "nonmissin real vector")
@p <- argvalue($2,"argument 2","positive count")

@nospec <- keyvalue($K,"nospec*", "TF", default:F)
@nfreq <- keyvalue($K,"nfreq*","positive count")
if (@nospec){
	if($v > 2){
		print("WARNING: argument 3 to $S ignored with 'nospec:T'",macroname:T)
	}
}else{
	if ($v > 2) {
		if (!isnull(@nfreq)){
			error("you must have only 2 non-keyword arguments with 'nfreq:L'")
		}
		@nfreq <- argvalue($03,"argument 3","positive count")
		if (@nfreq != goodfactors(@nfreq)){
			error(paste("number of frequencies",@nfreq,\
				"has a prime factor > 29; try", goodfactors(@nfreq)))
		}
	} elseif (!isnull(@nfreq)) {
		if (@nfreq != goodfactors(@nfreq)){
			error(paste("nfreq = ",@nfreq," has prime factor > 29; ",\
				"try nfreq:", goodfactors(@nfreq),sep:""))
		}
	} elseif (alltrue(isscalar(S,real:T), S > 2, S == floor(S))) {
		if (S != goodfactors(S)){
			error(paste("S =", S,"has prime factor > 29; try S =",\
				goodfactors(S)))
		}
		@nfreq <- S
	}else{
		@nfreq <- goodfactors(2*nrows(@y))
	}
}

if(!ismacro(autocov)){
	getmacros(autocov,quiet:T,printname:F)
}
@y <- autocov(@y, @p)

@phi <- @y[-1]/@y[1]
@sigma <- @y[1]*prod(1 - partacf(@phi)^2)
delete(@y,@p)

@phi <- yulewalker(@phi)

@result <- if (@nospec) {
	structure(phi:@phi,var:@sigma)
} else {
	structure(phi:@phi,var:@sigma,\
		spectrum:@sigma*hreal(hprdhj(rft(autoreg(@phi,padto(1,@nfreq))))))
}
delete(@phi,@sigma,@nospec)
delete(@result,return:T)
%arspectrum%

====> evalpoly <====
evalpoly        MACRO   DOLLARS
) Macro to evaluate polynomials with given coefficients
) Usage:
)  v <- evalpoly(c,z)   or   v <- evalpoly(c,x,T)
)   coef     REAL n by p matrix with no MISSING elements, with c[,j]
)            defining real polynomial Pj(w) = w^n - c[1,j]*w^(n-1) -
)            . . . - c[n-1,j]*w - c[n,j]
)   z        REAL m by 2*p matrix with no MISSING elements,
)            considered as an m by p complex matrix whos
)   x        REAL m by p matrix with no MISSING elements considered
)            as a real matrix
)
) With two arguments, v will be a m by 2*p matrix considered as m by p
) complex matrix with v[,J] = Pj(z[,J]), J = vector(2*j-1,2*j)
) With three arguments, v is m by p with v[,j] = Pj(x[,j]) 
)
) evalpoly(coef,polyroot(coef)) should be zero within rounding error
)
)) 990324; uses argvalue()
)) 000910: stripped $$
#usage: $S(coef,z) or evalpoly(coef,x,T)
if ($v < 2 || $v > 3){
	error("usage: $S(coef,z) or $S(coef,x,T)", macroname:F)
}
@c <- argvalue($1,"argument 1", "nonmissing real matrix")
@x <- argvalue($2,"argument 2", "nonmissing real matrix")

@cmplx <- if ($v == 3){
	!argvalue($03,"argument 3", "TF")
}else{
	T
}

if(ncols(@x) != {
	if(@cmplx){
		2
	}else{
		1
	}
}*ncols(@c)){
	error("wrong number of columns in 2nd argument to $S")
}
@n <- nrows(@c)
if(@cmplx){
	@zero <- 0*@c[1,]
}
@s <- @x
for(@i,run(@n-1)){
	@s <- if(@cmplx){
		cprdc(@s-cmplx(@c[@i,],@zero),@x)
	} else {
		(@s - @c[@i,])*@x
	}
}
@s <-- if(@cmplx){
	cmplx(@c[@n,],@zero)
}else{
	@c[@n,]
}
delete(@c,@x,@n,@cmplx)
delete(@s,return:T)
%evalpoly%

====> compfa <====
compfa          MACRO   DOLLARS
) Compute smoothed periodogram and optionally smoothed cross
) periodogram using cosine data taper and polynomial detrending.
) Usage:
)   Sy <- compfa(y, edf [,nfreq:Nfreq] [,degree:D] [,alpha:A] [,cross:T])
) with
)   y         REAL vector of length N or REAL N by p matrix with no
)             MISSING values
)   edf >= 0  REAL scalar <= .5 or >= 1 specifying equivalent
)             degrees of freedom or bandwidtg
)   Nfreq > 0 integer with no divisors > 29, the number of frequencies at
)             which spectrum is to be estimated
)   D         integer scalar, specifying degree for polynomial
)             detrending; D < 0 means no detrending
)   A >= 0    REAL scalar between 0 and .5, specifying amount of
)             cosine tapering.  Approximately A*N values are tapered
)             at each end
)
) When edf >= 1, it specifies the approximate EDF (equivalent degrees
) of freedom) of the estimates,
) When 0 <= edf <= .5, edf specifies the approximate bandwidth in units
) of cycles/delta_t; the EDF is computed as EDF = 2*edf*N
)
) Estimates are computed at the Nfreq frequencies 0, 1/Nfreq, 2/Nfreq,
) ..., (Nfreq - 1)/Nfreq cycles per delta_t.  See below for the default
) value of Nfreq
)
) When y is a vector, Sy is a vector of length Nfreq containing the
) smoothed periodogram of y.  Otherwise, the value of Sy depends on
) whether or not cross:T is an argument.
)
) cross:T is an argument
)  Sy is a nfreq by Q matrix with Q = p + p*(p-1)/2 = p*(p+1)/2 columns,
)  where p = ncols(y).
)   Cols. 1, 2, ..., p are the estimated spectra in Real form.
)   Cols. p+1, p+2, ..., p + p*(p-1)/2 are the estimated cross spectra
)   of y[,i] and y[,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, Sy <- compfa(hconcat(x,y), edf)
)  computes Sy = hconcat(Sxx, Syy, Sxy), where Sxx and Syy are estimated
)  spectra and Sxy is the estimated cross spectrum.
)
) cross:T not an argument:
)  Sy is Nfreq by p matrix with Sy[,j] containing the smoothed
)  periodogram of y[,j].  y can also be a N by p1 by p2 by ... by pk
)  array, in which case Sy is a Nfreq by p1 by p2 by ... by pk array,
)  and Sy[,j1,j2,...,jk] contains the smoothed periodogram of
)  y[,j1,j2,...,jk]
)
) When nfreq:Nfreq is not an argument
)   Nfreq = S if S is defined and is a positive integer.  It is an error
)     if S has a prime factor > 29
)   Nfreq = goodfactors(2*dim(y)[1]) otherwise, that is the smallest
)     integer >= 2*dim(y)[1] with no prime factors > 29.
)
) Uses macro getmacros to retrieve compza if not present.
)) Version of 990204 fixed potential bug
)) 990324 converted to use argvalue()
))        explicit check for factors in nfreq
)) 000910 stripped $$
)) 001231 Added check on value of S
)) 030408 Added keyword 'silent' as a synonym for 'quiet'
# usage: $S(y, edf [, degree:D, alpha:A, S:nfreq [,cross:T]]) , y a
# REAL vector or matrix, nfreq > 0, D integers, edf >=0, 0 <= A <= .5
if($v != 2 || $k > 4){
	error("usage is $S(y, edf [,S:n] [,degree:D] [,alpha:a][,cross:T])",\
		macroname:F)
}
@edf <- argvalue($02,"edf", "nonnegative number")
if (@edf == 0){
	@edf <- 2
}
@bwidth <- (@edf <= .5)
if (!@bwidth && @edf < 2){
	error("edf not a REAL scalar, 0 <= edf <= .5 or edf >= 2")
}

@quiet <- keyvalue($K,"quiet","TF")
@silent <- keyvalue($K,"silent","TF")
if (alltrue(!isnull(@quiet), !isnull(@silent), @quiet != @silent)) {
	error("Inconsistent usage of 'quiet' and 'silent'")
}
if (!isnull(@silent)) {
	@quiet <- @silent
}
delete(@silent)
@cross <- keyvalue($K,"cross","TF", default:F)
@nfreq <- keyvalue($K, "nfreq*","integer scalar")
@nfreq1 <- keyvalue($K, "S","integer scalar")

@y <- argvalue($1,"argument 1",vector("real","matrix","nonmissing"))
if (@cross && isvector(@y)){
	error("with cross:T, argument 1 for $S must not be vector")
}

@n <- nrows(@y)
if (!isnull(@nfreq) && !isnull(@nfreq1)) {
	if (@nfreq != @nfreq1) {
		error("values of keywords 'S' and 'nfreq' differ; don't know which to use")
	}
}
if (isnull(@nfreq)){
	@nfreq <- @nfreq1
}
if (!isnull(@nfreq)){
	if (@nfreq < 0) {
		@nfreq <- abs(@nfreq) # don't check value
	} elseif (@nfreq == 0) {
		error("number of frequencies is 0")
	} elseif (@nfreq != goodfactors(@nfreq)) {
		error(paste("number of frequencies",@nfreq,\
					"has a prime factor > 29; try", goodfactors(@nfreq)))
	}
} elseif (alltrue(isscalar(S,real:T),S > 0 && S == floor(S))) {
	if (S != goodfactors(S)){
		error(paste("S =",@S,"has prime factor > 29;",@nfreq,"try S =",\
			goodfactors(S)))
	}
	@nfreq <- S
} else {
	@nfreq <- goodfactors(2*@n)
}

if(!ismacro(compza)){
	getmacros(compza,quiet:T,printname:F)
}

  # compute Fourier transform of tapered detrended data
@za <- if($k > 0){
	compza(@y,nfreq:-@nfreq,$K)
}else{
	compza(@y,nfreq:-@nfreq)
}
@n <- @za$n
@alpha <- @za$alpha
if (@bwidth){
	@edf <- .5*@edf*@n
}
@nfreq <- dim(@za$za)[1]
if (!@cross || isvector(@za$za)){
	# compute modified periodograms
	@ia <- hreal(hprdhj(delete(@za,return:T)$za))
}else{
	@p <- ncols(@za$za)
	@ia <- hconcat(delete(@za, return:T)$za,\
		padto(rep(0,@p*(@p-1)/2)',@nfreq))
    @k <- @p
	# compute modified cross periodograms
	for (@i,run(@p-1)){
		@l <- run(1,@p - @i)
		@ia[,	@k + @l] <- hprdhj(@ia[,@i],@ia[,@i + @l])
		@k <-+ length(@l)
	}
	# compute modified periodograms
	@ia[,run(@p)] <- hreal(hprdhj(@ia[,run(@p)]))
}
if (@edf != 2){ #do some smoothing
	@ra <- (1 - 1.25*@alpha)^2/(1 - 1.453125*@alpha)  # always >= 1
	@factor <- @edf/@ra # always <= @edf
	@len <- ceiling((.5*@nfreq/@n)*(@factor/2.086 + 0.967/@factor))
	if (@len > 1){
		@shift <- 2*(@len-1)
		@wts <- rep(1/@len,@len)
		@wts <- convolve(@wts,vector(@wts,rep(0,@shift/2)))
		@wts <- convolve(@wts,vector(@wts,rep(0,@shift)))
		@edf <- (2*@n/@nfreq)*@ra/sum(@wts^2)
		if (!@cross){
			@fa <- rotate(convolve(delete(@wts,return:T),\
				delete(@ia,return:T)),-@shift)
		}else{
			@ia[,run(@p)] <- rotate(convolve(@wts,@ia[,run(@p)]),\
				-@shift)
			@ia[,-run(@p)] <-\
			  ctoh(rotate(convolve(delete(@wts,return:T),\
				htoc(@ia[,-run(@p)])), -@shift))
			rename (@ia, @fa)
		}
		delete(@shift)
	}else{
		@edf <- 2
	}
	delete(@ra, @factor,silent:T)
}
if (@edf == 2){
	@len <- 1
	rename(@ia,@fa)
}
if(isnull(@quiet)){
	@quiet <- @edf == 2
}
if(!@quiet){
	print(paste("rep(1/",@len, ",", @len, ")^*4 smoother with ",\
		round(@edf,1), " edf, S = ",@nfreq,sep:""))
}
delete(@len,@quiet,@edf,@bwidth,@nfreq,silent:T)
delete(@fa,return:T)
%compfa%

====> compza <====
compza          MACRO   DOLLARS
) Macro to compute the Fourier transform of a cosing tapered time series
) Usage:
)   compza(y [,nfreq:Nfreq] [,degree:D] [,alpha:A])
)   y         REAL vector or matrix with no MISSING elements and N rows
)   Nfreq     integer != 0 such that abs(Nfreq) has no prime factors > 29,
)             the number of frequencies at which the Fourier transform is
)             to be computed
)   D         integer degree of polynomial used to detrend the columns
)             of y; D < 0 means nothing is subtracted. Default is D = 0
)   A         REAL scalar, 0 <= A <= .5 specifying the amount of cosine
)             tapering; approximately A*N observations are tapered at
)             both ends of y.  Default is A = 0
)
) With nfreq:Nfreq
)  When Nfreq < 0, no checking for factors is done and the result
)  is computed at abs(Nfreq) frequencies.
) Without nfreq:Nfreq
)  If S exists and is a positive integer, Nfreq = S; it is an error if S
)  as a prime factor > 29.  Otherwise Nfreq = goodfactors(2*N).
)
) Value
)  structure(za:Za, n:nrows(y), ka:Ka, alpha:A, degree:D)
)   Za is a Nfreq by ncols(y) REAL vector or matrix whose columns are
)   the Fourier transforms of the corresponding columns of y after
)   detrending and tapering.  They are scaled by dividing by sqrt(Ka)
)   Ka = sum(taper^2), where taper is the vector of tapering weights
)
)) Uses following macros
))  detrend, costaper, testnfreq
)) Version of 980430
))  Tests for factors > 29 in nfreq
))  Uses argvalue()
))  Permits degree < 0 meaning do no detrending at all
)) Version 000910: stripped $$
)) 001215 modified default number of frequencies
)) 010101 keyword 'nfreq' now synomymous with 'S' which is deprecated;
))        is an error if value has prime factors > 29
# usage: $S(y [,nfreq:Nfreq] [,degree:D] [,alpha:A])
if($v != 1){
	error("usage is $S(y [,nfreq:Nfreq] [,degree:D] [,alpha:A])")
}

@alpha <- keyvalue($K, "alpha", "nonnegative number", default:0)
if(@alpha > .5){
	error("value alpha must be between 0 and .5")
}
@degree <- keyvalue($K, "degree", "integer scalar", default:0)
@nfreq <- keyvalue($K, "nfreq*","integer scalar")
@nfreq1 <- keyvalue($K, "S","integer scalar")
if (!isnull(@nfreq) && !isnull(@nfreq1)){
	if (abs(@nfreq) != abs(@nfreq1)) {
		error("values of 'S' and 'nfreq' differ; don't know which to use")
	}
}

@y <- argvalue($1,"y", "nonmissing real matrix")

@n <- nrows(@y)
if (isnull(@nfreq)) {
	@nfreq <- @nfreq1
}
delete(@nfreq1)
if (!isnull(@nfreq)){
	if (@nfreq < 0) {
		@nfreq <- abs(@nfreq) # don't check value
	} elseif (@nfreq != goodfactors(@nfreq)) {
		error(paste("number of frequencies",@nfreq,\
					"has a prime factor > 29; try", goodfactors(@nfreq)))
	}
} elseif (alltrue(isscalar(S,real:T),S > 0 && S == floor(S))) {
	if (S != goodfactors(S)){
		error(paste("S =",@S,"has prime factor > 29;",@nfreq,"try S =",\
			goodfactors(S)))
	}
	@nfreq <- S
} else {
	@nfreq <- goodfactors(2*@n)
}

if (@degree >= 0){
	if(!ismacro(detrend)){
		getmacros(detrend,quiet:T,printname:F)
	}
	@y <- detrend(@y, @degree)
}
if(@alpha > 0){
	if(!ismacro(costaper)){
		getmacros(costaper,quiet:T,printname:F)
	}
	@taper <- costaper(@n,@alpha)
	@y <-* @taper
	@ka <- sum(@taper^2)
	delete(@taper)
}else{
	@ka <- @n
}
@y <-/ sqrt(@ka)

structure(za:rft(padto(delete(@y,return:T),delete(@nfreq,return:T))),\
	n:delete(@n,return:T),ka:delete(@ka,return:T),\
	alpha:delete(@alpha,return:T),degree:delete(@degree,return:T))
%compza%

====> costaper <====
costaper        MACRO   DOLLARS
) Compute cosine taper (data window)
) Usage:
)  taper <- costaper(n, alpha)
)   n > 0       integer
)   alpha >= 0  scalar <= .5
)
) The result is a vector of length n with approximately alpha*n elements
) on each end tapered and with the approximately (1 - 2*alpha)*n middle
) elements = 1.  alpha = 0.5 means the entire result is tapered.
)
) NOTE: This differs from another convention in which alpha is
)       the total proportion tapered, counting both ends.
)) 990324 uses argvalue()
)) 000910 stripped $$
# usage: $S(n, alpha), integer n > 0, scalar alpha between 0 and .5
if($v != 2){
	error("usage is $S(n, alpha)", macroname:F)
}
@n <- argvalue($01,"n","positive count")
@alpha <- argvalue($02,"alpha","nonnegative number")
if (@alpha > 0.5){
	error("alpha not a REAL scalar between 0 and .5")
}
if(@alpha > 0){
	@len <- ceiling(@alpha*@n - 1e-8)
	@taper <- vector(.5*(1-cos(.5*(run(.5,@len))/@len, cycles:T)),\
		rep(1,@n-@len))
	delete(@len)
	@taper <-* reverse(@taper)
}else{
	@taper <- rep(1,@n)
}
delete(@alpha, @n)
delete(@taper,return:T)
%costaper%

====> dpss <====
dpss           MACRO   DOLLARS
) macro to compute Discrete Prolate Spheroidal Sequences (dpss)
) usage:
)  dpss(N, W, K)   or   dpss(N, W, K, First)
)    N >= 1        integer length of tapers
)    1 <= K <= N   integer number of tapers computed
)    0 < W < .5    scalar half bandwidth
)    1 <= First <= N - K + 1,
)                  integer starting index for tapers
) This computes K dpss sequences of length N, starting with the
) First-th
)
)) Version of 990324
))  uses argvalue()
)) Version 000910: stripped $$
)) 010102 checks that K <= N, plus some cosmetic changes 
# usage: $S(N, W, K [,FirstVec]), 0 < W < .5, N,K,FirstVec integers > 0
if ($k > 0 || $v > 4 || $v < 3){
	error("usage: $S(N, W, K [,FirstVec])", macroname:F)
}
@N <- argvalue($01,"argument 1 (N)", "positive count")

@W <- argvalue($02,"argument 2 (W)", "positive number")
if (@W >= .5){
	error("argument 2 (W) not scalar between 0 and .5")
}
@K <- argvalue($03,"argument 3 (K)", "positive count")
if (@K > @N) {
	error("K = argument 3 > N = argument 1")
}

if ($v > 3){
	@first <- argvalue($4,"argument 4 (First)","positive count")
	if(@first > @N - @K + 1){
		error("argument 4 (First) not positive integer <= N-K+1")
	}
} else {
	@first <- 1
}

@d <- cos(@W, cycles:T)*(.5*run(-@N+1,@N-1,2))^2

@e <- .5*run(0, @N-1)*run(@N, 1)
@e <- trideigen(@d, @e, @first, @first+@K-1, values:F)

for(@i,run(@K)){
	@k <- @i + @first - 2
	@d <- @e[,@i]
	if ((@k %& 1) == 0){
		if (sum(@d) < 0){
			@d <- @e[,@i] <- -@d
		}
	}elseif (sum(run(@N-1,-@N+1,-2)*@d) < 0){
		@d <- @e[,@i] <- -@d
	}
}
delete (@d, @N, @W, @first, @K, @i, @k)
delete(@e,return:T)
%dpss%

====> burg <====
burg     MACRO   DOLLARS
) Macro to carry out "maximum entropy" spectrum estimation using the
) Burg algorithm
) Usage:
)  burg(y, p [,degree:D] [,nospec:T] [,nfreq:Nfreq])
)   y          REAL vector or matrix with no MISSING elements
)   p > 0      Integer specifying the order of the AR model to be fit
)   D          Integer specifying degree of polynomial detrending
)              D < 0 means nothing is subracted, not even a mean
)   Nfreq > 0  Integer with no prime factors > 29 specifying the
)              number of frequencies the spectrum will be estimated at
)
) Optional keyword phrase nospec:T suppresses the computation of
) a spectrum and keyword phrase nfreq:Nfreq is ignored if present.
)
) The result is structure(phi:Phi, var:Var, spectrum:Sy) or, with
) nospect:T, structure(phi:Phi, var:Var), where Phi is the vector of
) estimated AR coefficients, Var is an estimate of the innovation
) variance and Sy, each column of which is the estimated spectrum for
) the corresponding column of y.
)
) With nfreq:Nfreq
)  The spectrum will be estimated at Nfreq frequencies. It is an error
)  if Nfreq has a prime factor > 29.
)
) Without optional keyword phrase nfreq:Nfreq:
)  If S exists and is a positive integer, Nfreq = S; it is an error if S
)  as a prime factor > 29.  Otherwise Nfreq = goodfactors(2*N).
)
) The value for var is prod(1-phikk^2)*sum(detrend(y,degree:D)^2)/N
) where phikk is vector of partial autocorrelations and N = nrows(y)
) This differs from some other implementations of the Burg algorithm
)) 990115
))   new keywork nospec:T, uses argvalue()
)) 990324
))   Explicit test for absence of prime factor > 29 in nfreq
))   Slight change in definition of variance; divide by N, not N-1
)) 000910 stripped $$
)) 010101 added keyword nfreq to replace S and modified determination
))        of number of frequencies
# $S(y,nlags [,degree:m][,nospec:T][, nfreq:nf])
#   returns structure(phi, var [, spectrum])
if ($v != 2){
	error("usage is $S(y,p [,degree:m] [,S:nfreq])", macroname:F)
}

@y <- argvalue($1,"$1","nonmissing real matrix")

@nlags <- argvalue($2,"$2", "positive count")

@m <- @n <- nrows(@y)
@ncols <- ncols(@y)

@degree <- keyvalue($K,"degree","integer scalar", default:0)

@nospec <- keyvalue($K,"nospec*", "TF", default:F)

@nfreq <- keyvalue($K,"S","positive count")
@nfreq1 <- keyvalue($K,"nfreq*","positive count")
if (@nospec) {
	@what <- if (!isnull(@nfreq)) {
		"S"
	} elseif (!isnull(@nfreq1)) {
		"nfreq"
	} else {
		NULL
	}
	if (!isnull(@what)) {
		print(paste("WARNING: value of '",@what,"' ignored with nospec:T",\
				sep:""))
	}
	delete(@what,@nfreq1)
} else {
	if (!isnull(@nfreq) && !isnull(@nfreq1)){
		if (@nfreq != @nfreq1) {
			error("values of 'S' and 'nfreq' differ; don't know which to use")
		}
	}
	if (isnull(@nfreq)){
		@nfreq <- @nfreq1
	}
	delete(@nfreq1)
	if (!isnull(@nfreq)){
		if (@nfreq != goodfactors(@nfreq)) {
			error(paste("number of frequencies",@nfreq,\
					"has a prime factor > 29; try", goodfactors(@nfreq)))
		}
	} elseif(alltrue(isscalar(S,real:T),S > 0 && S == floor(S))) {
		if (S != goodfactors(S)){
			error(paste("S =",@S,"has prime factor > 29;",@nfreq,"try S =",\
				goodfactors(S)))
		}
		@nfreq <- S
	} else {
		@nfreq <- goodfactors(2*@n)
	}
}

if (@degree == 0){
	@y <-- sum(@y)/@n
} elseif (@degree > 0) {
	if(!ismacro(detrend)){
		getmacros(detrend,quiet:T,printname:F)
	}
	@y <- detrend(@y, @degree)
}

@var <- sum(@y^2)/@n # 990324: changed (@n-1) to @n

@e <- @y[-1,]
@f <- @y[-@m,]
for(@k, run(@nlags)){
	@m <-- 1
	@phikk <- 2*sum(@e*@f)/(sum(@e^2)+sum(@f^2))
	@var <-* 1 - @phikk*@phikk
	@phi <- if (@k == 1){
		@phikk
	}else{
		vconcat(@phi - @phikk*reverse(@phi), @phikk)
	}
	if (@k < @nlags){
		@tmp <- (@e - @phikk * @f)[-1,]
		@f <- (@f - @phikk * @e)[-@m,]
		@e <- @tmp
	}
}
if (@nlags > 1){
	delete(@tmp)
}

if (!@nospec){
	@spectrum <- matrix(rep(0,@nfreq*@ncols),@nfreq)
	for (@k, run(@ncols)){
		@spectrum[,@k] <- vector(@var[@k]*hreal(hprdhj(rft(autoreg(\
			@phi[,@k],padto(1,@nfreq))))))
	}
	if (@ncols == 1){
		@spectrum <- vector(@spectrum)
	}
}

if (@ncols == 1){
	@phi <- vector(@phi)
}

delete(@k,@e,@f,@phikk,@y,@m,@n,@degree,@nlags)


@result <- if (@nospec) {
	structure(var:vector(delete(@var,return:T)), phi:delete(@phi,return:T))
} else {
	structure(var:vector(delete(@var,return:T)),phi:delete(@phi,return:T),\
		spectrum:delete(@spectrum,return:T))
}
delete(@nospec)
delete(@result,return:T)
%burg%

====> multitaper <====
multitaper      MACRO   DOLLARS
) macro to compute multitaper spectrum estimates using discrete
) prolate spheroidal sequences (dpss).  It requires function
) trideigen(), present in MacAnova versions released after 5/20/95
) Usage:
)  multitaper(y, W, K [,degree:D] [,nfreq:Nfreq] [,deltat:dt] [,wts:wts)
)   y         REAL vector or matrix with no MISSING elements, whose
)             whose columns are discrete parameter time series equally
)             spaced in time
)   W > 0     REAL scalar, half bandwidth in cycles per unit time
)   K > 0     integer scalar, the number of tapers to use
)   D         integer scalar, the degree of a polynomial to be fit for
)             detrending; D < 0 means nothing subtracted; D = 0 is default
)   Nfreq > 0 integer scalar with no prime factors > 29, the number of
)             frequencies at which the spectrum will be estimated; see
)             below for default
)   dt > 0    REAL scalar, the interval between observation times;
)             default is 1 or DELTAT
)   wts       REAL vector of length K with w[i] > 0; default is rep(1/K,K)
)
) The result is a nfreq by N = ncols(y) vector or matrix.
)
) Without deltat:dt, dt is taken to be variable DELTAT or 1 if DELTAT
) does not exist
)
) Half Bandwidth W is in units of cycles per unit time, not cycles
) per dt time units.  It must be between 0 and .5/dt.
)
) Without optional keyword phrase nfreq:Nfreq:
)  If S exists and is a positive integer, Nfreq = S; it is an error if S
)  as a prime factor > 29.  Otherwise Nfreq = goodfactors(2*N).
)
)) This macro uses macros dpss and detrend.  If they are not defined, an
)) attempt is made to read them from tser.mac using macro getmacros
)) 990324
))  use argvalue()
))  use testnfreq to check for prime factors > 29
)) 000910 stripped $$
)) 010101 modified choice of nfreq
# usage: multitaper(y, W, K [,degree:degree] [,S:nFreq])
if ($v != 3){
	error("usage is $S(y, W, K [,degree:degree] [S:nFreq])", macroname:F)
}
@y <- argvalue($01,"argument 1 (y)", "nonmissing real matrix")

@posint <- vector("positive","integer","scalar")
@W <- argvalue($02,"argument 2 (W)", "positive number")

@K <- argvalue($03, "argument 3 (K)", "positive count")

if(!ismacro(dpss)){
	getmacros(dpss,quiet:T,printname:F)
}

@degree <- keyvalue($K,"degree", "count", default:0)

@wts <- keyvalue($K,"wts","positive vector", default:rep(1,@K))
if (length(@wts) != @K){
	error(paste("value for wts not positive vector of length", @K))
}
@wts <-/ sum(@wts)

@n <- nrows(@y)
@p <- ncols(@y)

@nfreq <- keyvalue($K,"S","positive count")
@nfreq1 <- keyvalue($K,"nfreq*","positive count")
if (!isnull(@nfreq) && !isnull(@nfreq1)){
	if (@nfreq != @nfreq1) {
		error("values of 'S' and 'nfreq' differ; don't know which to use")
	}
}
if (isnull(@nfreq)){
	@nfreq <- @nfreq1
}
delete(@nfreq1)
if (!isnull(@nfreq)){
	if (@nfreq != goodfactors(@nfreq)) {
		error(paste("number of frequencies",@nfreq,\
				"has a prime factor > 29; try", goodfactors(@nfreq)))
	}
	if (@nfreq < @n) {
		error(paste("number of frequencies =",@nfreq,"< N =",@n))
	}
} elseif(alltrue(isscalar(S,real:T),S > 0, S == floor(S), S >= @n)) {
	if (S != goodfactors(S)){
		error(paste("S =",@S,"has prime factor > 29;",@nfreq,"try S =",\
			goodfactors(S)))
	}
	@nfreq <- S
} else {
	@nfreq <- goodfactors(2*@n)
}

@deltat <- keyvalue($K,"deltat", "positive number")
if (isnull(@deltat)){
	if (isscalar(DELTAT, real:T)){
		@deltat <- if (DELTAT > 0){
			DELTAT
		}else{
			1
		}
	}else{
		@deltat <- 1
	}
}

if(@W >= .5/@deltat){
	error("width $2 >= .5/deltat")
}

if (@degree >= 0){
	if(!ismacro(detrend)){
		getmacros(detrend,quiet:T,printname:F)
	}
	@y <- detrend(@y,@degree)
}

@h <- dpss(@n,@deltat*@W,@K)
delete(@n,@W,@degree,@deltat)

if (@p == 1){
	@result <- hreal(hprdhj(rft(padto(delete(@h,return:T)*\
					delete(@y,return:T),@nfreq)))) %*% @wts
}else{
	@result <- @wts[1]*hreal(hprdhj(rft(padto(@h[,1]*@y,@nfreq))))
	if (@K > 1){
		for(@j,2,@K){
			@result <-+ @wts[@j]*\
				hreal(hprdhj(rft(padto(@h[,@j]*@y,@nfreq))))
		}
	}
	delete(@y, @h, @j,silent:T)
}
delete(@p,@K,@nfreq,@wts)
delete(@result,return:T)
%multitaper%

