info  MACRO
) File containing macros for least squares nonlinear fitting
) including fitting ARIMA models to time series by unconditional
) least squares.
) Version of 030813
)
) Copyright (C) 1999, 2000, 2001, 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                                     *
) *********************************************************************
)
) This file serves as its own help file.
)
) With the standard configuration you can get help and usage on a topic,
) say hannriss, by arimahelp(hannriss) and arimahelp(hannriss,usage:T)
)
) You can get a list of all topics by arimahelp(index:T) or simply
) by arimahelp()
)
) Otherwise you can use help(file:"Arima.mac",arima),
) usage(file:"Arima.mac",hannriss) and help(file:"Arima.mac",arima_index)
)
)                 SIGN CONVENTIONS FOR COEFFICIENTS
) macros arima, hannriss, innovest, acfarma, and all recognize keywords
) keywords 'arsign' and 'masign' with values +1 or -1
)   arsign has default -1 or variable ARSIGN if it exists
)   masign has default -1 or variable MASIGN if it exists
) These affect the interpretation of autoregressive coefficients phi
) and moving average coefficients theta
) The AR operator is
)     Y[t] + Arsign*(phi[1]*Y[t-1] + phi[2]*Y[t-2] + ...)
) The MA operator is
)     Z[t] + Masign*(theta[1]*Z[t-1] + theta[2]*Z[t-2] + ...)
)
) Macros included are
) arima*        Unconditional least squares and MLE estimation of
)               ARIMA model or linear regression with ARIMA errors
) hannriss*     ARIMA fitting using Hannan-Rissannen algorithm
) innovest      ARIMA fitting using innovations algorithm
) neg2logLarma* Computes -2*log(L) and other stuff for ARIMA model
)               with given coefficients;
) arimares      Computes residuals from ARIMA model;
) innovations   Computes "innovations" algorithm to compute the
)               coefficients for one step prediction in terms of previous
)               one-step prediction errors.  Used by innovest
) moveoutroots  Fixes up coefficients for a MA or AR operator so that
)               all the roots are outside the unit circle in the
)               comples plane
) acfarma*      Compute autocovariance function of ARMA model
) specarma*     Compute spectrum of ARMA model
) rhatvar       Bartlett's formula for variances of estimated
)               autocorrelations
) rhatcovar     Bartlett's formula for variances or covariances of
)               estimated autocorrelations or entire variance matrix
) _polish       Macro used by hannriss and innovest to gave a final
)               adjustment to estimates
)
) * = affected by MASIGN and ARSIGN
)
)) Version of 990112
)) Version of 990127 arima and arimares allow MLE estimation
)) Version of 990129 Added hannriss, innovest, innovations and _polish
))   and made bug fixes.
)) Version of 990131 Regularize checking for existence of macros
)) Version of 990202 Added moveoutroots, modified innovest to use it
))   added some checks for stationarity/invertibility, and changed
))   the way the determinant of the covariance matrix is computed in
))   arimares.
)) Version of 990206 Minor changes to hannriss, innovest and _polish
)) Version of 990207 Modified innovations to use cholesky().  This
))   resulted in a huge speed up
)) Version of 990211 Pure AR and MA models forced to be stationary/
))   invertible before polishing in hannriss and innovest
)) Version of 990304 Fixed bugs in arimares; removed " in $S" from
))   calls to error()
)) Version of 000909 Stripped $$, updated use of keyvalue() and
))   enabled alternative sign conventions of coefficients using
))   keywords 'arsign' and 'masign'
)) Version of 000914, modified use of arsign and masign so that
))   defaults are -1 and -1
)) Version of 001003, added copy of arimahelp() from MacAnova.mac
)) Version of 001010, fixed bug in arima
)) Version of 001016, modified innovest help
)) Version of 001019, fixed bug in arimares
)) Version of 001117, added neg2logLarma and fixed another arimares bug
)) Version of 001130, modified MLE iteration in arima and changed
)     computation of logdet and fixed bugs in arimares
)) Version of 001207 fixed bug in arima
)) Version of 001213 added subtopics to help
)) Version of 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
)) Version of 010505 corrected bug in arima
)) Version of 010724 corrected bug in computing HII in arima with mle:T
))   also XTXINV and hessian do not include the first row of residuals
))   computed by levmar()
))   The hessian returned by levmar() is now computed from the jacobian
))   that is returned
)) Version of 030320 removed levmar(), _cgrad() and _lmout() (moved to
))   math.mac; added more keys to help entries
))   Also removed nlreg(), linear(), asymptot(), drapSmthFunc(), testfun(),
))   SandCT19.8.1, DandS_T10.2 (moved to regress.mac)
)) Version of 030428 Minor change to innovest()
)) 030814 added subtopic titles to help
%info%

===> arimahelp <===
arimahelp     MACRO
) Macro to obtain help on the macros in file arima.mac
) Usage:
)    arimahelp(topic1 [, topic2 ...] [,scrollback:T])
)      prints help on topics in file arima.mac
)    arimahelp(topic1 [, topic2 ...], usage:T [,scrollback:T])
)      prints usage information on topics in file arima.mac
)    arimahelp(index:T [,scrollback:T])
)      prints index of topics in the file arima.mac
)) Version 990929
# usage $S(topic1 [, topic2 ...] [help keywords])
#       $S(topic1 [, topic2 ...], usage:T) with no other keywords
#       $S(index:T)
if(!ismacro(_gethelp)){
	getmacros(_gethelp,quiet:T,printname:F)
}
___HELPFILE_ <- "arima.mac"
___INDEXTOP_ <- "arima_index"
___MACRO_ <- "$S"
_gethelp($0)
%arimahelp%

===> arima <===
arima         MACRO DOLLARS
) Least squares and MLE ARIMA model fitting using nonlinear least squares
) arima(y [,pdq:pdq] [,PDQ:PDQ,seasonal:period, [,x:x,fitmean:T,\
)         start:b0,active:active,cast:n, cycles:m,mle:T,mlecycles:m1,\
)         masign:Masign, arsign:Arsign,\
)         maxit:itmax, minit:itmin,crit:crvec,print:T,keep:T,quiet:T])
) y       REAL vector containing the time series
) pdq     vector(p,d,q), p, d, q nonnegative integers
) PDQ     vector(P,D,Q), P, D, Q nonnegative integers
) period  integer > 1; ignored when PDQ:PDQ is not an argument
) x       Optional REAL vector or matrix whose columns are linear
)         predictors; nrows(x) = nrows(y)
) fitmean When T, mean will be fitted.  Default is T when d = 0,
)         F, otherwise
) b0      Optional REAL vector of starting values
) active  LOGICAL vector the same length as b0
) cast    Value is how far back casting will be done in computing
)         residuals
) cycles  Non-negative integer specifying number of forecasting/
)         backcasting cycles (default 0)
) mle     When T, use approximate Maximum likelihood estimation
) mlecycles Number of updates of sigmahatsq allowed in final
)         stage of MLE (only with mle:T).  Default is 0, meaning
)         sigmahatsq as estimated using unconditional LS estimates
)         is used while including log determinant in objective
)         function
) Masign  -1 (default) or +1, affects definition of MA parameters.  For
)         MA(2), for example, the model is
)          Y(t) = Z(t) + Masign*theta[1]*Z(t-1) + Masign*theta[2]*Z(t-2)
) Arsign  -1 (default) or 1, affects definition of AR parameters.  For
)         AR(2), for example, the model is
)          Y(t) = Z(t) - Arsign*phi[1]*y(t-1) - Arsign*phi[2]*y(t-2)
) itmax   Maximum number of iterations allowed (default 30, 0 is ok)
)         With mle:T, when itmax < 0, no iterations of least squares
)         fitting is done and up to abs(itmax) iterations of
)         MLE using sigmahat based on coefficients in b0.  This
)         is helpful in re-running arima with mle:T when convergence
)         was not achieved.
) itmin   Minimum number of iterations carried out (default 0)
) crvec   vector(numsig, nsigsq, delta), convergence criteria
)               Converge when                             Default
)         Relative change in all coefs < 10^-numsig     numsig =  5
)         Relative change in RSS < 10^-nsigsq           nsigsq = 10
)         ||gradient|| < delta (ignored when delta = 0)  delta =  0
) print   When T, partial results printed at each iteration
) keep    When T, arima returns the structure returned by levmar()
)         with components, coefs, hessian, jacobian, gradient,
)         residuals, nobs, npar, pdq, PDQ, seasonal, active, iter, iconv,
)         rss, neg2logL, aicc
) quiet   When F (default, unless keep:T), no summary results are
)         printed, although side effect variables are computed.
)
) If variable ARSIGN and/or MASIGN exist with value +1 or -1 they are used as
) defaults for 'arsign' and 'masign' instead of -1 and -1.
)
) arima creates the following side effect variables
)   COEF = vector(phihat,thetahat)
)   XTXINV = analogue of solve(X' %*% X) matrix in regression, with rows
)    and columns corresponding to inactive parameters set to 0
)   ALLRESIDUALS = residuals from fitted model including back cast
)    residuals
)   RESIDUALS = ALLRESIDUALS without back cast and incomplete
)    residuals
)   HII, REAL vector of leverages
)   RSS = sum(ALLRESIDUALS^2)
)   NEG2LOGL = -2*log(likelihood)
)   NPAR = number of active parameters + 1 for the mean if estimated
)     as sample mean as it will be when active[1] is F and a mean is
)     estimated.
)   JACOBIAN = REAL matrix of derivatives of ALLRESIDUALS with respect to
)    active parameters
)   GRADIENT = REAL gradient vector with an element for each active
)    parameter
)   AICC = value of AICC
)
)) Other macros used (should be loaded automatically if available)
))  arimares()
))    may use acfarma()
))  levmar()
))    uses _cgrad(), may use _lmout()
) Written by C. Bingham, January 1999
)
)) 000908 Stripped $$, updated keyword handling
)) 001010 corrected bug in using @signs inactive parameters
)) 001130 Modified final stage when log(det) is included in
))        objective function
)) 001207 fixed bug in computing log(det) with mle:F and either
))        Masign = 1 or Arsign = 1
)) 010103 minimum default for cast increased from 15 to 50
)) 010510 corrected bug by adding @x[1,] to call to arimares
)) 010724 bug fix; now trim of row 1 of jacobian when mle:T
))        also remove contribution of this row from hessian and
))        hence from xtxinv.  New side effect variable AICC
))        The change in xtxinv affects the estimated standard errors
))        and t-statistics.
# $S(y [,pdq:pdq] [,PDQ:PDQ,seasonal:period], [,x:x] [,fitmean:T,\
#         start:b0,active:active,cast:n, cycles:m,mlecycles:m1,\
#         masign:Masign, arsign:Arsign,\
#         maxit:itmax, minit:itmin,crit:crvec,print:T,keep:T,quiet:T])
# y       REAL vector containing the time series
# pdq     vector(p,d,q), p, d, q nonnegative integers
# PDQ     vector(P,D,Q), p, d, q nonnegative integers
# period  integer > 1; ignored when PDQ:PDQ is not an argument
# x       Optional REAL vector or matrix whose columns are linear
#         predictors; nrows(x) = nrows(y)
# fitmean When T, mean will be fitted.  Default is T when d = 0,
#         F, otherwise
# b0      Optional REAL vector of starting values
# active  LOGICAL vector the same length as b0
# cast    How far back casting will be done
# cycles  number of forecasting/backcasting cycles (default 0)
# mle     When T, use approximate Maximum likelihood estimation
# mlecycles Number of updates of sigmahatsq allowed in final
#         stage of MLE (only with mle:T).  Default is 0, meaning
#         sigmahatsq estimated using unconditional LS estimates
#         is used while including log determinant in objective
#         function
# Masign  -1 (default) or +1, affects definition of MA parameters.  For
#         MA(2), for example, the model is
#          Y(t) = Z(t) + Masign*theta[1]*Z(t-1) + Masign*theta[2]*Z(t-2)
# Arsign  -1 (default) or 1, affects definition of AR parameters.  For
#         AR(2), for example, the model is
#          Y(t) = Z(t) - Arsign*phi[1]*y(t-1) - Arsign*phi[2]*y(t-2)
# itmax, itmin Limits on number of iterations (defaults 30 and 0)
# crvec   vector(numsig, nsigsq, delta), 3 criteria for convergence
# print   When T, partial results printed at each iteration
# keep    When T, $S returns the structure returned by levmar() with the
#         the addition of component 'edf' = error df
# quiet   When F (default, unless keep:T), no summary results printed
if ($v != 1){
	error("$S must have exactly 1 non-keyword argument",macroname:F)
}

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

@keys <- if ($k > 0) {
	structure($K)
} else {
	structure(notakey:NULL)
}

@pdq <- keyvalue(@keys,"pdq","nonneg integer vector", default:rep(0,3))
if (length(@pdq) > 3){
	error("length(pdq) > 3")
}
@pdq <- padto(@pdq, 3)
@p <- @pdq[1]
@d <- @pdq[2]
@q <- @pdq[3]

@PDQ <- keyvalue(@keys,"PDQ","nonneg integer vector", default:rep(0,3))
if (length(@PDQ) > 3){
	error("length(PDQ) > 3")
}
@PDQ <- padto(@PDQ, 3)
@P <- @PDQ[1]
@D <- @PDQ[2]
@Q <- @PDQ[3]

@fitmean <- keyvalue(@keys,"fitmean","TF",default:@d + @D == 0)

@seasonal <- keyvalue(@keys,"season*","positive count")
if (alltrue(!isnull(@seasonal), @seasonal == 1)){
	error("value for 'seasonal' must be > 1")
}
if (isnull(@seasonal)) {
	@seasonal <- 0
}

if (@P + @D + @Q > 0){
	if (@seasonal == 0){
		error("'seasonal:period' required with non-zero 'PDQ'")
	}
}elseif (@seasonal > 0){
	print("WARNING: value of 'seasonal' ignored with zero 'PDQ'")
	@seasonal <- 0
}

@mle <- keyvalue(@keys,"mle","TF", default:F)
@mlecycles <- keyvalue(@keys,"mlecycle*","count")
if (!@mle) {
	if (!isnull(@mlecycles)){
		print("WARNING: keyword 'mlecycles' ignored without 'mle:T'",\
			macroname:T)
	}
	delete(@mlecycles)
} else {
	@mlecycles <- if (isnull(@mlecycles)) {
		1
	} else {
		@mlecycles + 1
	}
}

@param <-\
	structure(pdq:@pdq,PDQ:@PDQ,seasonal:@seasonal,fitmean:@fitmean)

@cast <- keyvalue(@keys,"cast","count", default:\
	if (@p + @P == 0){
		@q + @seasonal*@Q
	}else{
		max(50, @q + @seasonal*@Q + 4*(@p + @seasonal*@P))
	})

@param <- strconcat(@param, cast:@cast,\
				  cycles:keyvalue(@keys,"cycle*","count", default:0))

@maxit <- keyvalue(@keys,"maxit", "integer scalar", default:30)

@minit <- keyvalue(@keys,"minit","count", default:0)

 #  Default         Converge when
 #  crvec[1] =  5   Relative change in all coefs < 1e-5 = 10^-crvec[1]
 #  crvec[2] = 10   Relative change in RSS < 1e-10 = 10^-crvec[2]
 #  crvec[3] =  0   ||gradient|| < crvec[3] (ignored when crvec[3] = 0)
@crit <- keyvalue(@keys,"crit*","nonmissing real vector",\
	default:vector(5, 10, 0))

if (length(@crit) > 3){
	error("value of 'crit' not vector with length <= 3")
}
@crit <- padto(@crit,3)

@b0 <- keyvalue(@keys,"start*","nonmissing real vector")

@active <- keyvalue(@keys,"active","nonmissing logical vector")

@print <- keyvalue(@keys,"print","TF",default:F)

@keep <- keyvalue(@keys,"keep","TF",default:F)

@quiet <- keyvalue(@keys,"quiet","TF",default:@keep)

@sign <- if (alltrue(isscalar(MASIGN, real:T),!anymissing(MASIGN))) {
	MASIGN
} else {
	-1
}
@masign <- keyvalue(@keys,"masign","number", default:@sign)
if (abs(@masign) != 1) {
	error("MASIGN or value for 'masign' != +- 1")
}
@sign <- if (alltrue(isscalar(ARSIGN, real:T),!anymissing(ARSIGN))) {
	ARSIGN
} else {
	-1
}
@arsign <- keyvalue(@keys,"arsign","number", default:@sign)
delete(@sign)
if (abs(@arsign) != 1) {
	error("ARSIGN or value for 'arsign' != +- 1")
}
@x <- keyvalue(@keys,"x","real nonmissing matrix")

delete(@keys)

@signs <- if (@fitmean) {
	1
} else {
	NULL
}
if (!isnull(@x)) {
	if (nrows(@x) != nrows(@y)){
		error("nrows(x) != nrows(y)")
	}
	@signs <- vector(@signs,rep(1,ncols(@x)))
} else {
	@x <- 0
}

if (@p > 0) {
	@signs <- vector(@signs, rep(-@arsign,@p))
}
if (@q > 0) {
	@signs <- vector(@signs, rep(-@masign,@q))
}
if (@P > 0) {
	@signs <- vector(@signs, rep(-@arsign,@P))
}
if (@Q > 0) {
	@signs <- vector(@signs, rep(-@masign,@Q))
}

delete(@arsign, @masign)
if (isnull(@signs)) {
	error("no coefficients to fit")
}

@npar <- length(@signs)

if (isnull(@b0)){
	@b0 <- rep(0, @npar) #starting values
}elseif (length(@b0) < @npar){
	@b0 <- padto(@b0, @npar)
}elseif (length(@b0) > @npar){
	error("more starting values than there are coefficients")
}

@b0 <-* @signs
if (@fitmean && @b0[1] == 0){
	@b0[1] <- ?
}

if (isnull(@active)){
	@active <- rep(T,@npar)
}elseif(length(@active) != @npar){
	error("length of value of active != number of coefficients")
}
@signs1 <- @signs[@active]
@nactive <- sum(@active)
if (@nactive == 0){
	error("all elements of 'active' are False")
}

@ybarcomp <- ismissing(@b0[1]) || @fitmean && @npar == 1

@N <- nrows(@y)
@nobs <- @N - @d - @D*@seasonal # nobs after differencing

@minNobs <- @fitmean + @p + @q + (@P + @Q)*(@seasonal + 1)
if (@nobs <= @minNobs){
	error(paste("you need at least",@minNobs,"observations for this model"))
}else{
	delete(@minNobs)
}

if (@ybarcomp){
  # if you estimate mu by ybar, count it as an estimated parameter
	if (!@active[1]){
		@nactive <-+ 1
	}
	if (@d + @D == 0){
		@b0[1] <- sum(@y)/@nobs
		if (@npar == 1){
			@residuals <- @y - @b0[1]
			@rss <- sum(@residuals^2)
		}
	}else{
		@dif <- @y
		#compute ordinary differences
		if (@d > 0){
			for(@j,1,@d){
				@dif <- movavg(1,@dif)[-1]
			}
		}
		#compute seasonal differences
		if (@D > 0){
			@S <- reverse(padto(1,@seasonal))
			for(@j,1,@D){
				@dif <- movavg(@S,@dif)[-run(@seasonal)]
			}#compute seasonal differences
			delete(@S)
		}
		delete(@j)
		@b0[1] <- sum(@dif)/@nobs
		if (@npar == 1){
			@residuals <- @dif - @b0[1]
			@rss <- sum(@residuals^2)
		}
		delete(@dif)
	}
}# if (@ybarcomp)

if (@npar > 1 || !@fitmean){
	if (!ismacro(levmar)){
		getmacros(levmar,quiet:T,printname:F)
	}
	if (!ismacro(arimares)){
		getmacros(arimares,quiet:T,printname:F)
	}

	@maxit1 <- if(@maxit < 0){
		0
	}else{
	  	@maxit
	}
	@minit1 <- if(@minit < 0){
		0
	}else{
	  	@minit
	}

	@result <- levmar(@b0,@x,@y,@param,resid:arimares,\
					  crit:@crit, active:@active,\
					  maxit:@maxit1, minit:@minit1,print:@print)

	@iter <- @result$iter

	if (@mle && (@result$iconv > 0 || @maxit <= 0)){
		@maxit1 <- if (@maxit > 0){
			max(1,@maxit - @iter)
		}else{
			-@maxit
		}
		@maxit <- abs(@maxit)
		# On each cycle @sigmahat is kept fixed and levmar() is allowed
		# two iterations except on the last one
		# The number of cycles is controled by keyword 'mlecycles'

		@sigmahat <- sqrt(sum(@result$residuals^2)/@nobs)

		@param <- strconcat(@param,sigmahat:0)
		for (@i, 1, @mlecycles) {
			@param[ncomps(@param)] <- @sigmahat
			@maxit2 <- if (@i < @mlecycles){
				2
			} else {
				@maxit1
			}
			@result <- levmar(@result$coefs,@x,@y,@param,resid:arimares,\
							  crit:@crit, active:@active,\
							  maxit:@maxit2,minit:@minit,print:@print)
			@iter <-+ @result$iter
			@sigmahat <- sqrt(sum(@result$residuals[-1]^2)/@nobs)
		}
		delete(@mlecycles, @maxit2, @i)
		@maxit <- -@maxit # signal we reached MLE part of iteration
	}# if (@mle && @result$iconv > 0)
	delete(@maxit1, @minit1)
}else{ # if (@npar > 1 || !@fitmean)
	@iter <- 0
	@result <- structure(coefs:@b0,\
		hessian:@nobs,\
		jacobian:rep(-1, @nobs),\
		gradient:if(@ybarcomp){0}else{?},\
		rss:@rss,\
		residuals:delete(@residuals,return:T),\
		nobs:@N,\
		iter:if(@ybarcomp){1}else{0},\
		iconv:1)
}
delete(@b0, @y)

ALLRESIDUALS <- @result$residuals

if (@mle && isdefined(@sigmahat)){
	@logdet <- (ALLRESIDUALS[1]/@sigmahat)^2
	ALLRESIDUALS <- ALLRESIDUALS[-1]
	@jacindex <- match("jacobian",compnames(@result),0)
	@hesindex <- match("hessian",compnames(@result),0)
	@result[@jacindex] <- @result[@jacindex][-1,]
	@result[@hesindex] <- @result$jacobian %c% @result$jacobian
	delete(@jacindex,@hesindex)
}else{
	@logdet <- arimares(@result$coefs,@x[1,],0,\
		structure(pdq:@pdq,PDQ:@PDQ,seasonal:@seasonal,cast:@cast,\
		fitmean:@fitmean,cycles:0,sigmahat:-@N))^2
	# sigmahat:-N means just compute sqrt(logdet)
}

COEF <- @result$coefs * @signs

@xtxinv <- solve(@result$hessian)
XTXINV <- dmat(@npar, 0)
XTXINV[@active,@active] <- (@signs1 * @xtxinv) * @signs1'

@edf <- @nobs - @nactive

@sigmahat <- sqrt(sum(ALLRESIDUALS^2)/@nobs)

@twoLogL <- -2*@nobs*log(@sigmahat) - @logdet - @nobs -\
		@nobs*log(2*PI)

AICC <- -@twoLogL + 2*@nactive * @nobs/(@nobs - 1 - @nactive)
NEG2LOGL <- -delete(@twoLogL,return:T)
NPAR <- delete(@nactive,return:T)
RSS <- sum(ALLRESIDUALS^2)

@nx <- if(!isscalar(@x)){
	ncols(@x)
}else{
	0
}

JACOBIAN <- @result$jacobian * @signs1'
GRADIENT <- @result$gradient * @signs1

@extra <- nrows(ALLRESIDUALS) - @nobs

if (@extra > 0){
	RESIDUALS <- ALLRESIDUALS[-run(@extra)]
	@jacobian <- JACOBIAN[-run(@extra),]'
}else{
	RESIDUALS <- ALLRESIDUALS
	@jacobian <- JACOBIAN'
}
HII <- vector(sum(@jacobian* (@xtxinv %*% @jacobian)))
delete(@extra,@jacobian)


@iconv <- @result$iconv
@hessian <- @signs1 * (@result$hessian * @signs1')
@residuals <- @result$residuals
delete(@result)

if (!@quiet){
	@mse <- sum(RESIDUALS^2)/@edf

	@method <- if(@mle){
		@objfun <- "-2logL"
		"Maximum likelihood"
	}else{
		@objfun <- "RSS"
		"Unconditional least squares"
	}

	@what <- if (@d + @D + @P + @Q == 0){
		paste("ARMA(",@p,",",@q,")",sep:"")
	}elseif (@P + @Q + @D == 0){
		paste("ARIMA(",@p,",",@d,",",@q,")",sep:"")
	}elseif (@p + @q + @d == 0){
		paste("Seasonal ARIMA(",@P,",",@D,",",@Q,\
			")<",@seasonal,">",sep:"")
	}else{
		paste("Seasonal ARIMA(",@p,",",@d,",",@q, ")x(",\
			@P,",",@D,",",@Q,")<",@seasonal,">",sep:"")
	}

	print(paste(@method, @what))
	delete(@what,@method)

	@fmt <- getoptions(format:T)
	@xtrwid <- 2*(@seasonal > 1)
	@width <- floor(vecread(string:@fmt,silent:T))
	print(paste(charwidth:5+@xtrwid," ",charwidth:@width,\
		"Coef","StdErr","t","P Value",justify:"r"))

	@se <- sqrt(@mse*diag(@xtxinv))
	@tstat <- COEF[@active]/@se

	@ix <- @fitmean
	@ip <- @ix + @nx
	@iq <- @ip + @p
	@iP <- @iq + @q
	@iQ <- @iP + @P

	@j <- @k <- 0

	for(@i,1,@iQ + @Q){
		@j <-+ 1
		@k <-+ @active[@j]
		if(@active[@j] || COEF[@j] != 0){
			@lab <- if(@j > @iQ){
				@l <- @iQ
				paste("MA",@seasonal,intwidth:2,sep:"")
			}elseif(@j > @iP){
				@l <- @iP
				paste("AR",@seasonal,intwidth:2,sep:"")
			}elseif(@j > @iq){
				@l <- @iq
				"MA"
			}elseif(@j > @ip){
				@l <- @ip
				"AR"
			}elseif(@j > @ix){
				@l <- @ix
				"B"
			}else{
				@l <- 0
				"Mu"
			}
			@tmp <- if (@j == @fitmean){
				paste(charwidth:5+@xtrwid,@lab,format:@fmt,COEF[@j])
			}else{
				paste(charwidth:2+@xtrwid,@lab,format:"2.0f",\
					@i-@l,format:@fmt,COEF[@j])
			}
			if (@active[@j] && @edf > 0){
				@tmp <- paste(@tmp,format:@fmt,\
					@se[@k], @tstat[@k],\
					2*(1-cumstu(abs(@tstat[@k]),@edf)))
			}
			print(@tmp)
		}#if(@active[@j] || COEF[@j] != 0)
	}# for(@i,1,@iq + @q)
	delete(@i, @j, @k, @l)

	print(paste(rep("-",9+@xtrwid+4*@width),sep:""))

	@tmp <- paste("MSE: ", @mse, ", DF: ", @edf)
	@tmp <- paste(@tmp, ", -2*log(L): ", NEG2LOGL,", AICC: ",AICC,\
			sep:"")
	print(@tmp)

	@tmp <- paste("Complete RSS:",RSS)
	print(@tmp)
	@tmp <- paste("N: ",@N, ", Backcasting horizon: ",@cast, sep:"")

	if (@d + @D > 0){
		@tmp <- paste(@tmp, ", N after differencing: ",	@nobs, sep:"")
	}
	print(delete(@tmp, return:T))

	if(@iconv <= 0){
		if (@maxit != 0){
			@msg <- paste("Did not converge in",@iter,"iterations")
			if (@mle){
				@part <- vector("LS","ML")[(@maxit < 0) + 1]
				@msg <- paste(@msg,"during",delete(@part,return:T),\
					"part of iteration")
			}
		}else{
			@msg <- "No iteration done"
		}
		print(delete(@msg,return:T))
	}elseif (@iconv <= 3){
		@msg <- if (@iconv == 1){
			paste("relative change in all coefs <", 10^-@crit[1])
		}elseif (@iconv == 2){
			paste("relative change in",@objfun, "<", 10^-@crit[2])
		}else{
			paste("norm of gradient =",sqrt(sum(GRADIENT^2)),"<",\
				@crit[3])
		}
		print(paste("Converged with",@msg,"in", @iter,"iterations"))
		delete(@msg)
	}else{
		print(paste("Halving step did not reduce",@objfun, "on",\
			@iter,"iteration"))
	}
	if (@p > 0 && @cast > 0){
		if (sum(@residuals[@mle + run(5)]^2)/\
			sum(@residuals^2) > 1e-3*5/length(@residuals)){
			print("WARNING: backward residuals not dying out rapidly enough")
		}
	}
	if (@p > 0){
		if (max(creal(cpolar(polyroot(COEF[@ip + run(@p)])))) > 1){
			print("WARNING: AR part of model is not stationary")
		}
	}
	if (@q > 0){
		if (max(creal(cpolar(polyroot(COEF[@iq+run(@q)])))) > 1){
			print("WARNING: MA part of model is not invertible")
		}
	}
	if (@P > 0){
		if(max(creal(cpolar(polyroot(COEF[@iP+run(@P)])))) > 1){
			print("WARNING: seasonal AR part of model is not stationary")
		}
	}
	if (@Q > 0){
		if (max(creal(cpolar(polyroot(COEF[@iQ+run(@Q)])))) > 1){
			print("WARNING: seasonal MA part of model is not invertible")
		}
	}
	delete(@ix, @ip, @iq, @iP, @iQ)
	delete(@fmt,@se,@tstat,@param,@npar, @width,@fitmean)
}# if (!@quiet)

delete(@p,@d,@q,@cast,@maxit,@mle,@objfun,@residuals,@signs,@signs1)

if (delete(@keep,return:T)){
	@result <- structure(coefs:COEF,\
		hessian:delete(@hessian,return:T),\
		jacobian:JACOBIAN,\
		gradient:GRADIENT,\
		residuals:ALLRESIDUALS,\
		nobs:delete(@nobs,return:T),\
		npar:NPAR,\
        pdq:delete(@pdq,return:T),\
        PDQ:delete(@PDQ,return:T),\
        seasonal:delete(@seasonal,return:T),\
		active:delete(@active,return:T),\
		iter:delete(@iter,return:T),\
		iconv:delete(@iconv,return:T),\
		rss:RSS,\
		neg2logL:NEG2LOGL,\
		aicc:AICC)
	delete(@result, return:T)
}else{
	delete(@edf, @nactive, @active, @pdq, @PDQ, @seasonal, silent:T)
}
%arima%

===> hannriss <===
hannriss      MACRO DOLLARS
) Experimental version of macro to compute the Hannan-Rissanen
) initial estimates of the parameters of an ARIMA series.
) Usage:
)   hannriss(x, pdq:vector(p,d,q) [,degree:degree] [,maxlag:M] \
)            [,arsign:1] [,masign:1] [,useburg:T] [,polish:T, cycles:nc])
)    p >= 0, d >= 0, q >= 0, M >= p + q, nc >= 0,degree are integers
)    masign:1 (default -1) alters definition of MA parameters.  For
)      MA(2), e.g., y(t) = Z(t) + theta[1]*Z(t-1) + theta[2]*Z(t-2)
)      instead of y(t) = Z(t) - theta[1]*Z(t-1) - theta[2]*Z(t-2)
)    arsign:1 (default -1) alters definition of AR parameters.  For
)      AR(2), e.g., y(t) = Z(t) - phi[1]*y(t-1) - phi[2]*y(t-2)
)      instead of y(t) = Z(t) + phi[1]*y(t-1) + phi[2]*y(t-2)
)
) If variable ARSIGN and/or MASIGN exist with value +1 or -1 they are used as
) defaults for 'arsign' and 'masign' instead of -1 and -1.
)
) The first step computes yulewalker estimates based on M autocorrelations
) unless q = 0 when only p autocorrelations are used.  If useburg:T,
) is an argument the Burg estimates are computed instead.
) The default for M is 20 + p + q
)
) The possibly differenced data is detrended with a polynomial of
) order degree.  degree < 0 means nothing is subtracted.
) In the absense of degree:degree, the default value for degree is 0
) when d = 0 and -1 when d > 0
)
) With polish:T, it includes the extra step described on page 155
) of Brockwell & Davis.  With cycles:nc the extra step is iterated nc times.
)
) Returns structure(phi:phihat,theta:thetahat,xtxinv,nobs:n,npar:np,
)   rss:residSS,neg2logL:-2log(likelihood),aicc:AICC),
) where np = p + q + degree1 + 1 with degree1 = degree if degree > 0
) and = -1 otherwise, residSS is sum of all squared residuals,
) including backcast residuals, likelihood is the normal likelihood of
) the fitted model
)
) theta is defined so Theta(B) is 1 - theta[1]*B - theta[2]*B^2 ...
)
) When p = q = 0, phi, theta, xtxinv are NULL
)
) hannriss creates the following side effect variables
)   COEF = vector(phihat,thetahat) (NULL when p = q = 0)
)   XTXINV = analogue of solve(X' %c% X) matrix in regression
)          = NULL when p = q = 0
)   ALLRESIDUALS = residuals from fitted model including back cast
)    residuals
)   RSS = sum(ALLRESIDUALS^2)
)   NEG2LOGL = -2*log(likelihood)
)   NPAR = p + q + degree1 + 1 = number of coefficients estimated
)
)) Other macros used (should be loaded automatically if available)
))  detrend
))  autocov (without burg:T)
))  burg    (with burg:T)
))  _polish
))  arimares (copied to macro _resid)
))    may use acfarma
) There is currently no provision for seasonal ARIMAs
) Version 000908
# $S(x,p ,q [,maxlag:m] [,useburg:T] [,polish:T])
if ($v != 1){
	error("usage: $S(x, pdq:vector(p,d,q) [,degree:d], [,maxlag:m] [,polish:T,ncycles:nc)",\
		macroname:F)
}
@x <- argvalue($1,"$1","nonmissing real vector")

@keys <- if ($k > 0) {
	structure($K)
} else {
	structure(notakey:NULL)
}

@pdq <- keyvalue(@keys,"pdq","nonnegative integer vector", default:rep(0,3))
if (length(@pdq) > 3){
	error("length of value of 'pdq' > 3")
}
@pdq <- padto(@pdq,3)
@p <- @pdq[1]
@d <- @pdq[2]
@q <- @pdq[3]
delete(@pdq)

@npar <- @p + @q

@m <- keyvalue(@keys,"maxlag","positive count", default:\
	if(@q == 0){
		max(@p,1)
	}else{
		max(20, @npar)
	})

@polish <- keyvalue(@keys,"polish","TF", default:F)

@cycles <- keyvalue(@keys,"cycle*","positive count")
if (isnull(@cycles)){
	@cycles <- 1
}elseif(!@polish){
	print("WARNING: value of 'cycles' ignored without 'polish:T'",macroname:T)
}

@degree <- keyvalue(@keys,"degree","integer scalar", default:-@d)

@useburg <- keyvalue(@keys,"useburg","TF", default:F)

if (@d > 0){
	@diffop <- padto(1, @d+1)
	for(@i,1,@d){
		@diffop <- movavg(1,@diffop)
	}
	@diffop <- -@diffop[-1]
	@x <- movavg(@diffop,@x)[-run(@d)]
}
@nobs <- length(@x)

if (@degree >= 0){
	#either degree specified or @d > 1
	if(!ismacro(detrend)){
		getmacros(detrend,quiet:T,printname:F)
	}
	@x <- detrend(@x, @degree) # remove mean or trend
}

@sign <- if (alltrue(isscalar(MASIGN, real:T),!anymissing(MASIGN))) {
	MASIGN
} else {
	-1
}
@masign <- keyvalue(@keys,"masign","number", default:@sign)
if (abs(@masign) != 1) {
	error("MASIGN or value for 'masign' != +- 1")
}
@sign <- if (alltrue(isscalar(ARSIGN, real:T),!anymissing(ARSIGN))) {
	ARSIGN
} else {
	-1
}
@arsign <- keyvalue(@keys,"arsign","number", default:@sign)
delete(@sign)
if (abs(@arsign) != 1) {
	error("ARSIGN or value for 'arsign' != +- 1")
}

@maxpq <- max(@p,@q)
@offset <- keyvalue(@keys,"offset","count", default:@maxpq)
delete(@keys)

@n1 <- @nobs - @m - @offset

if (@npar > 0){
	if (@n1 <= @npar){
		error("sample size is too small for model specified")
	}
	if (@q > 0){
		if (!@useburg){
			if(!ismacro(autocov)){
				getmacros(autocov,quiet:T,printname:F)
			}
			@gamma <- autocov(@x,@m)
			@phi <- yulewalker(@gamma[-1]/@gamma[1])
			delete(@gamma)
		}else{
			if (!ismacro(burg)){
				getmacros(burg,quiet:T,printname:F)
			}
			@phi <- burg(@x,@m,nospec:T)$phi
		}
		@zhat <- movavg(@phi,@x)#estimated residuals
	}else{
		@zhat <- NULL
	}

	@J <- run(@n1) + @nobs - @n1

	 #build data matrix
	@X <- matrix(rep(0,@n1*(@npar+1)),@n1)
	@k <- 0
	if (@p > 0){
		# add columns of lagged x
		for(@i,1,@p){
			@k <-+ 1
			@X[,@k] <- @x[@J-@i]
		}
	}
	if (@q > 0){
		# add columns of lagged zhat
		for(@i,1,@q){
			@k <-+ 1
			@X[,@k] <- @zhat[@J-@i]
		}
	}

	@k <-+ 1
	@X[,@k] <- @x[@J] # add non-lagged x

	@XX <- swp(@X %c% @X, run(@npar))
	delete(@X)

	@phi <- if (@p > 0){
		vector(@XX[@k,run(@p)])
	}else{
		NULL
	}
	@theta <- if (@q > 0){
		-vector(@XX[@k,@p+run(@q)])
	}else{
		NULL
	}

	if (@p == 0){ # pure MA; make sure invertible
		if (!ismacro(moveoutroots)){
			getmacros(moveoutroots,quiet:T,printname:F)
		}
		@theta <- moveoutroots(@theta)
	}
	if (@q == 0){ # pure AR; make sure stationary
		if (!ismacro(moveoutroots)){
			getmacros(moveoutroots,quiet:T,printname:F)
		}
		@phi <- moveoutroots(@phi)
	}
	if (@polish){
		if (!ismacro(_polish)){
			getmacros(_polish,quiet:T,printname:F)
		}
		@polished <- _polish(@x,structure(phi:@phi,theta:@theta),\
			@p,@q,@cycles)
		@phi <- @polished$phi
		@theta <- @polished$theta
		@XX <- @polished$xtxswept
		delete(@polished)
	}# if (@polish)
}else{
	@phi <- @theta <- NULL
	@xtxinv <- NULL
}

 # compute residual SS by backcasting and also log(det)
if(!ismacro(arimares)){
	getmacros(arimares,quiet:T,printname:F)
}
@cast <- if (@p == 0){
	@q
}else{
	50
}
ALLRESIDUALS <- arimares(vector(@phi,@theta),0,@x,\
	structure(pdq:vector(@p,0,@q),seasonal:0,cast:@cast,fitmean:F,\
		cycles:0, sigmahat:1))

@logdet <- ALLRESIDUALS[1]^2
ALLRESIDUALS <- ALLRESIDUALS[-1]
@rss <- sum(ALLRESIDUALS^2)

@sigmahat <- sqrt(@rss/@nobs)
if (@degree >= 0){
	@npar <-+ @degree + 1
}
@neg2logL <- 2*@nobs*log(@sigmahat) + @logdet + @nobs +\
	@nobs*log(2*PI)
@aicc <- @neg2logL + 2*@npar * @nobs/(@nobs - 1 - @npar)
delete(@sigmahat, @logdet)

if (@p + @q > 0){
	@signs <- rep(vector(-@arsign, -@masign),vector(@p,@q))
	@J <- run(@p+@q)
	if (@p > 0) {
		@phi <-* -@arsign
	}
	if (@q > 0) {
		@theta <-* -@masign
	}
	XTXINV <- @xtxinv <- @signs * (@XX[@J,@J] * @signs')
	COEF <- vector(@phi, @theta)
	delete(@signs)
}else{
	XTXINV <- COEF <- NULL
}
delete(@arsign,@masign)
RSS <- @rss
NEG2LOGL <- @neg2logL
NPAR <- @npar

delete(@x,@zhat,@XX,@p,@d,@q,@J,@k,@i,@degree,\
	@maxpq,@cycles,silent:T)

structure(phi:delete(@phi,return:T),\
	theta:delete(@theta,return:T),\
	xtxinv:delete(@xtxinv,return:T),\
	nobs:delete(@nobs, return:T),\
	npar:delete(@npar, return:T),\
	rss:delete(@rss,return:T),\
	neg2logL:delete(@neg2logL,return:T),\
	aicc:delete(@aicc,return:T))
%hannriss%

===> innovest <===
innovest      MACRO DOLLARS
) Macro to compute the innovations preliminary estimate of
) coefficients for an ARMA(p,q) time series as described on pp 151-153
) of Brockwell and Davis
) Usage:
) innovest(x,pdq:vector(p,d,q) [,maxlag:M] [,degree:degree]\
)          [,arsign:1] [,masign:1] [,polish:T, cycles:nc]
)          [,checkroots:F] [,silent:T])
)  p >= 0, d >= q >= 0, M >= p+q, degree, c integers
) The default value for M is max(2*(p+q),17)
)  masign:1 (default -1) alters definition of MA parameters.  For
)    MA(2), e.g., y(t) = Z(t) + theta[1]*Z(t-1) + theta[2]*Z(t-2)
)    instead of y(t) = Z(t) - theta[1]*Z(t-1) - theta[2]*Z(t-2)
)  arsign:1 (default -1) alters definition of AR parameters.  For
)    AR(2), e.g., y(t) = Z(t) - phi[1]*y(t-1) - phi[2]*y(t-2)
)    instead of y(t) = Z(t) + phi[1]*y(t-1) + phi[2]*y(t-2)
)
) With degree:degree, a polynomial trend of order degree is removed
) possibly after differencing
) When degree < 0, nothing is subtracted not even a mean.
) The default value for degree is -d (mean removed when d = 0,
) nothing removed otherwise)
)
) With polish:T, does nc (default 1) repetitions of an approximate
) least squares iteration
)
) Unless checkroots:F is an argument, the coefficients phihat and
) thetahat are checked to see they define stationary/invertible lag
) operators.  If they do not, a warning message is printed, no
) 'polishing' step and rss, neg2logL, and aicc are set to MISSING.
) With 'silent:T', warning messages are suppressed.
)
) returns structure(phi:phihat, theta:thetahat, nobs:n, xtxinv:xtxinv,
)   npar:p+q+degree1+1,rss:residSS,neg2logL:-2*log(likelihood),
)   aicc:AICC),
)
) where degree1 = max(degree, -1) residSS = sum of squares of all
) residuals, including those back cast, likelihood is the likelihood
) of the fitted model
)
) Component xtxinv is NULL unless polish:T is an argument when it
) can be used as an approximate solve(X' %*% X) matrix in computing
) standard errors.
)
) theta is defined so Theta(B) is 1 - theta[1]*B - theta[2]*B^2 ...
) unless masign:1 is an argument or variable MASIGN exists and has value 1
)
) When p = q = 0, phi, theta, xtxinv are NULL
)
) innovest creates the following side effect variables
)   COEF = vector(phihat,thetahat)
)   ALLRESIDUALS = residuals from fitted model including back cast
)     residuals
)   RSS = sum(ALLRESIDUALS^2)
)   NEG2LOGL = -2*log(likelihood)
)   NPAR = p + q + degree + 1 = number of coefficients estimated
)
) There is currently no provision for seasonal ARIMAs
)
)) Other macros used (should be loaded automatically if available)
))  autocov
))  detrend
))  innovations
))  _polish
))  arimares
))    may use acfarma
) 0001210
)) p = q = 0 now permitted
# $S(x, pdq:vector(p,d,q) [,maxlag:m] [,degree:degree])
if ($v > 1){
	error("usage: $S(x, pdq:vector(p,d,q) [,maxlag:m] [,degree:degree])",\
		macroname:F)
}

@x <- argvalue($1,"$1",vector("real","vector","nonmissing"))

@keys <- if ($k > 0) {
	structure($K)
} else {
	structure(notakey:NULL)
}

@pdq <- keyvalue(@keys,"pdq","nonnegative integer vector")
if (isnull(@pdq)){
	error("pdq:vector(p,d,q) is a required argument to $S", macroname:F)
}
if (length(@pdq) > 3){
	error("length(pdq) > 3")
}
@pdq <- padto(@pdq,3)

@p <- @pdq[1]
@d <- @pdq[2]
@q <- @pdq[3]
delete(@pdq)

@npar <- @p + @q

@m <- keyvalue(@keys,"maxlag","positive count", \
		default: max(2*@npar,15 + @p + @q))
if (@m < @npar){
	error("value for 'maxlag' < p + q")
}

@polish <- keyvalue(@keys,"polish","TF", default:F)

@checkit <- keyvalue(@keys,"check*","TF",default:T)

@silent <- keyvalue(@keys,"silent","TF",default:F)

@cycles <- keyvalue(@keys,"cycle*","positive count")

if (isnull(@cycles)) {
	@cycles <- 1
} elseif(!@polish && !@silent){
	print("WARNING: value of 'cycles' ignored without 'polish:T'", macroname:T)
}

@sign <- if (alltrue(isscalar(MASIGN, real:T),!anymissing(MASIGN))) {
	MASIGN
} else {
	-1
}
@masign <- keyvalue(@keys,"masign","number", default:@sign)
if (abs(@masign) != 1) {
	error("MASIGN or value for 'masign' != +- 1")
}
@sign <- if (alltrue(isscalar(ARSIGN, real:T),!anymissing(ARSIGN))) {
	ARSIGN
} else {
	-1
}
@arsign <- keyvalue(@keys,"arsign","number", default:@sign)
delete(@sign)
if (abs(@arsign) != 1) {
	error("ARSIGN or value for 'arsign' != +- 1")
}
@signs <- rep(vector(-@arsign,-@masign),vector(@p,@q))

@degree <- keyvalue(@keys,"degree","integer scalar",default:-@d)

if (@d > 0){
	@diffop <- padto(1,@d+1)
	for(@i,1,@d){
		@diffop <- movavg(1,@diffop)
	}
	@diffop <- -@diffop[-1]
	@x <- movavg(@diffop,@x)[-run(@d)]
}

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

@nobs <- nrows(@x)

if (@npar > 0){
	if (!ismacro(autocov)){
		getmacros(autocov,quiet:T,printname:F)
	}
	@gamma <- autocov(@x, @m)

	if (!ismacro(innovations)){
		getmacros(innovations,quiet:T,printname:F)
	}
	@results <- innovations(delete(@gamma,return:T),lag:@m,final:T)
	@psihat <- vector(1, reverse(@results$theta)[run(@npar)])
	delete(@results)

	if (@p > 0){
		@origin <- max(@p - @q - 1, 0)
		if (@origin > 0){
			@psihat <- vector(rep(0,@origin), @psihat)
		}
		@J <- run(@p) + @q + @origin + 1

		@Psi <- matrix(rep(0,@p^2),@p)
		for(@i,1,@p){
			@Psi[,@i] <- @psihat[@J-@i]
		}

		@phi <- vector(solve(@Psi ,@psihat[@J]))
		@theta <- if (@q > 0){
			-movavg(@phi,@psihat[@origin+run(@q+1)])[-1]
		}else{
			NULL
		}
		@cast <- 30
		delete(@origin)
	}else{
		@phi <- NULL
		@theta <- -@psihat[1 + run(@q)]
		@cast <- @q
	}
	delete(@psihat)

	if (@p == 0){ # pure MA; make sure invertible
		if (!ismacro(moveoutroots)){
			getmacros(moveoutroots,quiet:T,printname:F)
		}
		@theta <- moveoutroots(@theta)
	}
	if (@q == 0){ # pure AR; make sure stationary
		if (!ismacro(moveoutroots)){
			getmacros(moveoutroots,quiet:T,printname:F)
		}
		@phi <- moveoutroots(@phi)
	}

	if (@polish){
		if (!ismacro(_polish)){
			getmacros(_polish,quiet:T,printname:F)
		}
		@polished <- _polish(@x,structure(phi:@phi,theta:@theta),\
			@p,@q,@cycles)
		@phi <- @polished$phi
		@theta <- @polished$theta
		@XX <- @polished$xtxswept
		delete(@polished)
	}# if (@polish)
	COEF  <- @signs * vector(@phi,@theta)
}else{
	COEF <- @phi <- @theta <- NULL
	@xtxinv <- NULL
}
NPAR <- @npar

if (@npar > 0){
	@ok <- T
	if (@checkit && (@polish || @p > 0 && @q > 0)){
		if(@p > 0){
			@roots <- polyroot(@phi)
			@ok <- max(hypot(creal(@roots),cimag(@roots))) < 1
			if (!@ok && !@silent){
				print("WARNING: Estimated AR operator is non-stationary",\
					macroname:T)
			}
		}
		if(@q > 0){
			@roots <- polyroot(@theta)
			if (max(hypot(creal(@roots),cimag(@roots))) > 1){
				@ok <- F
				if (!@silent){
					print("WARNING: Estimated MA operator is non-invertible",\
						macroname:T)
				}
			}
		}
	}# if (@checkit)
}else{# if (@npar > 0)
	@ok <- T
}
delete(XTXINV, silent:T)
if (@ok){
	 # compute residual SS by backcasting and also log(det)
	if (@p + @q > 0){
		@J <- run(@npar)
		if (@polish){
			@xtxinv <- @signs * (@XX[@J,@J] * @signs')
		}else{
			@xtxinv <- NULL
		}
	}else{
		@xtxinv <- COEF <- NULL
	}
	@cast <- if (@p > 0){
		30
	}else{
		@q
	}
	if(!ismacro(arimares)){
		getmacros(arimares,quiet:T,printname:F)
	}
	ALLRESIDUALS <- arimares(vector(@phi,@theta),0,@x,\
		structure(pdq:vector(@p,0,@q),seasonal:0,cast:@cast,fitmean:F,\
			cycles:0, sigmahat:1))

	@logdet <- ALLRESIDUALS[1]^2
	ALLRESIDUALS <- ALLRESIDUALS[-1]
	@rss <- sum(ALLRESIDUALS^2)

	@sigmahat <- sqrt(@rss/@nobs)
	if (@degree >= 0){
		@npar <-+ @degree + 1
	}
	@neg2logL <- 2*@nobs*log(@sigmahat) + @logdet + @nobs +\
		@nobs*log(2*PI)
	@aicc <- @neg2logL + 2*@npar * @nobs/(@nobs - 1 - @npar)
	delete(@sigmahat, @logdet, @degree)

	RSS <- @rss
	NEG2LOGL <- @neg2logL
}else{
	if (@polish){
		if (!@silent){
			print("WARNING: no 'polishing' step done")
		}
		@xtxinv <- matrix(rep(?,@npar^2), @npar)
	}else{
		@xtxinv <- NULL
	}
	RSS <- NEG2LOGL <- @rss <- @neg2logL <- @aicc <- ?
}

if (@polish){
	XTXINV <- @xtxinv
}

if (@p > 0) {
	@phi <-* -@arsign
}
if (@q > 0) {
	@theta <-* -@masign
}
delete(@x,@p,@q,@m,@arsign,@masign,@signs)

structure(phi:delete(@phi,return:T),\
	theta:delete(@theta,return:T),\
	xtxinv:delete(@xtxinv,return:T),\
	nobs:delete(@nobs, return:T),\
	npar:delete(@npar, return:T),\
	rss:delete(@rss,return:T),\
	neg2logL:delete(@neg2logL,return:T),\
	aicc:delete(@aicc,return:T))
%innovest%

===> neg2logLarma <===
neg2logLarma MACRO DOLLARS
) Macro to compute -2*log(L) and other quantities for an ARIMA model
) with specified parameters.  The calling sequence is very similar to
) that of arima, except you must provide values of the coefficients
) and there are keywords which specify what is returned.
)
) Usage:
)  neg2logLarma(y,coefficients, [,x:x], [pdq:vector(p,d,q)] \
)   [,PDQ:vector(P,D,Q),seasonal:s] [,fitmean:T or F] \
)   [,cast:m] [cycles:ncyc] [, sigmasq:sigmasq] [,neg2logL:F],\
)   [,residuals:T] [,logdet:T] [,sigmahatsq:T] [,all:T],\
)   [,masign:1 or -1] [,arsign:1 or -1])
) where
) y is a real nonmissing vector of length n (the response)
) x is a real nonmissing matrix with n rows (optional linear predictors)
) p, d, q are non-negative integers (specify non-seasonal ARIMA form)
) P, D, Q are non-negative integers (specify season ARIMA)
) s > 1 an integer; required if max(P,D,Q) > 0
) fitmean:T => model includes mean (to differences if d > 0 or D > 0)
) fitmean:F => model does not include mean
)   default for fitmean is T is d = D = 0; F otherwise
) m >= 0 an integer, specifying how far to backcast or forcast in
)   computing residuals
)   Without cast:m, an appropriate default is chosen.
) ncyc >= 0, the number of forecasting/backcasting cycles;
)   default is 0 (no forcasting, just backcasting)
) sigmasq real scalar  > 0 or MISSING; value for innovation variance
)   MISSING (default) means estimate from data and compute
)   concentrated likelihood
) What is returned is controlled by
)   'neg2logL', 'residuals', 'logdet', 'sigmahatsq'
) neg2logL      -2*log(L).  When sigmasq is MISSING, L is the
)               concentrated likelihood
) residuals     vector of residuals, including backcast ones
) logdet        log(det(Gamma)) where Gamma is covariance matrix
)               when innovation variance is 1
) sigmahatsq    The estimated innovation variance or sigmasq, when
)               sigmasq is not MISSING
) If more than one of these is T, the output is a structure with these
)   as component names; if only 1 is T, the output is a scalar or vector
) all:T implies neg2logL:T, residuals:T, logdet:T, sigmasq:T, but
)   a component can be suppressed by, say, residuals:F
) masign must be +1 or -1, default is MASIGN if it exists or -1
) arsign must be +1 or -1, default is ARSIGN if it exists or -1
)) 001117 modified from simpler earlier unreleased version
)) Version 001117 written by C. Bingham (kb@stat.umn.edu)
)) 010103 minimum default for cast increased from 15 to 50 as in arima
# neg2logLarma(y,coefficients, [,x:x] [,pdq:vector(p,d,q)] \
#   [,PDQ:vector(P,D,Q),seasonal:s] [,fitmean:T or F] \
#   [,cast:m] [,cycles:ncyc] [,sigmasq:sigmasq] [,neg2logL:F],\
#   [,residuals:T] [,logdet:T] [,sigmahatsq:T] [,all:T])
@y <- argvalue($1,"time series")
@n <- length(@y)
@b <- argvalue($2,"parameters") #(mu,beta,phi,theta,phiS,thetaS)
if (!isnull(@b)){
	@b <- argvalue(@b,"parameters","nonmissing real vector")
}
@keys <- if ($k > 0) {
	structure($K)
} else {
	structure(notakey:NULL)
}

@pdq <- keyvalue(@keys,"pdq","nonneg integer vector", default:rep(0,3))
if (length(@pdq) > 3){
	error("length(pdq) > 3")
}
@pdq <- padto(@pdq, 3)
@p <- @pdq[1]
@d <- @pdq[2]
@q <- @pdq[3]

@PDQ <- keyvalue(@keys,"PDQ","nonneg integer vector", default:rep(0,3))
if (length(@PDQ) > 3){
	error("length(PDQ) > 3")
}
@PDQ <- padto(@PDQ, 3)
@P <- @PDQ[1]
@D <- @PDQ[2]
@Q <- @PDQ[3]

@fitmean <- keyvalue(@keys,"fitmean","TF",default:@d + @D == 0)

@seasonal <- keyvalue(@keys,"season*","positive count")
if (alltrue(!isnull(@seasonal), @seasonal == 1)){
	error("value for 'seasonal' must be > 1")
}
if (isnull(@seasonal)) {
	@seasonal <- 0
}

if (sum(@PDQ) > 0){
	if (@seasonal == 0){
		error("'seasonal:period' required with non-zero 'PDQ'")
	}
}elseif (@seasonal > 0){
	print("WARNING: value of 'seasonal' ignored with zero 'PDQ'")
	@seasonal <- 0
}

 # these control calculation
@x <- keyvalue(@keys,"x","real nonmissing matrix")
@cast <- keyvalue(@keys,"cast","count", default:\
	if (@p + @P == 0){
		@q + @seasonal*@Q
	}else{
		max(50,@q + @seasonal*@Q + 4*(@p + @seasonal*@P))
	})

@cycles <- keyvalue(@keys,"cycle*","count",default:0)

@sigmahat <- keyvalue(@keys,"sigmasq*","real scalar",default:?)
@concent <- ismissing(@sigmahat)
if (alltrue(!@concent, @sigmahat <= 0)){
	error("Value for 'var' not ? or positive")
}
if (@concent){
	@sigmahat <- 1
} else {
	@sigmahat <- sqrt(@sigmahat)
}
 # these control output
@all <- keyvalue(@keys,"all","TF",default:F)
@neg2logL1 <- keyvalue(@keys,"neg2log*","TF",default:T)
@resid1 <- keyvalue(@keys,"resid*","TF",default:@all)
@logdet1 <- keyvalue(@keys,"logdet*","TF",default:@all)
@sigmasq1 <- keyvalue(@keys,"sigmahat*","TF",default:@all)
delete(@all)

@sign <- if (alltrue(isscalar(MASIGN, real:T),!anymissing(MASIGN))) {
	MASIGN
} else {
	-1
}
@masign <- keyvalue(@keys,"masign","number", default:@sign)
if (abs(@masign) != 1) {
	error("MASIGN or value for 'masign' != +- 1")
}
@sign <- if (alltrue(isscalar(ARSIGN, real:T),!anymissing(ARSIGN))) {
	ARSIGN
} else {
	-1
}

@arsign <- keyvalue(@keys,"arsign","number", default:@sign)
delete(@sign)
if (abs(@arsign) != 1) {
	error("ARSIGN or value for 'arsign' != +- 1")
}
delete(@keys)

@signs <- if (@fitmean) {
	1
} else {
	NULL
}

if (!isnull(@x)) {
	if (nrows(@x) != nrows(@y)){
		error("nrows(x) != nrows(y)")
	}
	@signs <- vector(@signs,rep(1,ncols(@x)))
} else {
	@x <- 0
}

if (@p > 0) {
	@signs <- vector(@signs, rep(-@arsign,@p))
}
if (@q > 0) {
	@signs <- vector(@signs, rep(-@masign,@q))
}
if (@P > 0) {
	@signs <- vector(@signs, rep(-@arsign,@P))
}
if (@Q > 0) {
	@signs <- vector(@signs, rep(-@masign,@Q))
}
delete(@arsign, @masign, @p, @q, @P, @Q)

@npar <- length(@signs)

if (length(@b) != @npar) {
	error("number of coefficients doesn't match model")
}
if (@npar != 0) {
	@b <-* @signs  # adjust signs of parameters
}
if (!ismacro(arimares)){
	getmacros(arimares,quiet:T,printname:F)
}
if (@npar > 0) {
	@params <- structure(pdq:@pdq,PDQ:@PDQ,cast:@cast,\
		fitmean:@fitmean, seasonal:@seasonal,cycles:@cycles,\
		sigmahat:@sigmahat)

	@y <- arimares(@b,@x,@y,@params)
	@logdet <- @y[1]^2/@sigmahat^2
	@y <- @y[-1]
	delete(@params)
} else { # completely null model
	@logdet <- 0
}
@rss <- sum(@y^2)
@n <-- @d + @D*(@seasonal)
delete(@d,@D)

if (@concent) {
	@sigmasq <- @rss/@n
	@rss <- 1
} else {
	@sigmasq <- @sigmahat^2
	@rss <- (@rss/@n)/@sigmasq
}
if (@neg2logL1) {
	@neg2logL <- @n*(log(2*PI) + log(@sigmasq) + @rss) + @logdet
}
@ncomps <- 0
@result <- NULL
if (@neg2logL1) {
	@result <- strconcat(@result,neg2logL:@neg2logL)
}
if (@logdet1) {
	@result <- strconcat(@result,logdet:@logdet)
}
if (@sigmasq1) {
	@result <- strconcat(@result,sigmasq:@sigmasq)
}
if (@resid1) {
	@result <- strconcat(@result,residuals:@y)
}
delete(@b,@y,@pdq,@PDQ,@cast,@n,@rss,@concent,@sigmahat,@sigmasq,@logdet,\
	@logdet1,@resid1,@sigmasq1)
@result <- @result[-1]
delete(@result,return:T)
%neg2logLarma%

===> detarma <===
detarma   MACRO DOLLARS
) Macro to compute determinant of Toeplitz covariance associated with
) ARMA process.  At present it has no provision for seasonal processes.
) Usage:
)  detarma(phi, theta, n [,masign:Masign] [,arsign:Arsign] [,nfreq:Nfreq])
)   phi      non-MISSING REAL vector defining a causal AR operator
)   theta    non-MISSING REAL vector defining a MA operator
)   n        vector of positive integers
)   Masign and/or Arsign must be +1 or -1; they determine sign conventions
)            for interpreting phi and theta.  Type arimahelp(MASIGN)
)            for information.
)   Nfreq    Number of frequencies used in DFT to compute ACFV.  Default
)            is the smallest integer >= 2*max(max(n), 200) that has no
)            prime factors > 29
)
) Result is a REAL vector d, the same length as n with
)   d[i] = det(Sigma(n[i])), Sigma(n[i]) = Cov[(X(1),...,X(n[i]))]
)   when X(t) is ARMA with AR coefficients phi and MA coefficients theta
)) Written by C. Bingham (kb@stat.umn.edu)
)) Version 001129
)) 010103 S is no longer used as default for Nfreq
#$S(phi,theta, n [,masign:Masign] [,arsign:Arsign] [,nfreq:Nfreq])
@phi <- argvalue($1,"phi","nonmissing real vector")
@theta <- argvalue($2,"theta","nonmissing real vector")
@n <- argvalue($3,"n","positive integer vector")
@sign <- if(isscalar(MASIGN,real:T)) {
	MASIGN
} else {
	-1
}
@masign <- keyvalue($K,"masign","integer scalar",default:@sign)
if (abs(@masign) != 1) {
	error("value for 'masign' not +1 or -1")
}
@sign <- if(isscalar(ARSIGN,real:T)) {
	ARSIGN
} else {
	-1
}
@arsign <- keyvalue($K,"arsign","integer scalar",default:@sign)
if (abs(@masign) != 1) {
	error("value for 'arsign' not +1 or -1")
}
delete(@sign)
@nfreq <- keyvalue($K,"nfreq*","positive count")
if (isnull(@nfreq)) {
	# no longer uses S
	@nfreq <- goodfactors(2*max(@n, 200))
} elseif (goodfactors(@nfreq) > @nfreq) {
	error("

if (!ismacro(acfarma)) {
	getmacros(acfarma,quiet:T,printname:F)
}
if (isscalar(@theta) && @theta[1] == 0){
	@nlags <- length(@phi)
	@n <- @nlags + (@n - @nlags)*(@n <= @nlags)
} else {
	@nlags <- max(@n) - 1
}
@gamma <- acfarma(@phi,@theta,lag:@nlags,\
	arsign:@arsign,masign:@masign,nfreq:@nfreq)

@tmp <- vector(log(@gamma[1]),log(1 - partacf(@gamma[-1]/@gamma[1])^2))
delete(@phi,@theta,@gamma,@masign,@arsign,@nfreq,@nlags)

@det <- @n
for(@i,1,length(@n)) {
	@J <- run(@n[@i])
	@det[@i] <- exp(sum(reverse(@J)*@tmp[@J]))
}
delete(@i,@n,@J,@tmp)
delete(@det,return:T)
%detarma%

===> arimares <===
arimares      MACRO DOLLARS
) Compute residuals from ARIMA model with given coefficients
) Usage:
)   residuals <- arimares(b,x,y,params)
) b        REAL vector of coefficients
) x        vector or matrix of linear predictors; if a scalar x is
)          ignored
) y        the discrete time series being analyzed
) params   structure(pdq:vector(p,d,q),PDQ:vector(P,D,Q),seasonal:n1,
)          cast:n2,cycles:m3,fitmean:T or F [,sigmahat:s])
)  p,d,q   integers >= 0, AR order, difference, MA order
)  P,D,Q   integers >= 0, seasonal AR order, seasonal difference,
)          seasonal MA order
)  n1      integer > 1, length of seasonal
)  n2      integer >= - specifies how far to backcast
)  n3      integer >= 0 number of cycles of fore/backcasting
)  fitmean:T => mean will be fit
)  s       REAL scalar, estimate of residual stddev.  Its presence
)          signals that log(det(covariance matrix)) is to be
)          computed. If s < 0, this is all that arimares
)          does, returning the scalar sqrt(log(det)) for M x M
)          covariance matrix, where M = N - d - n1*D, N = -s.
) length(b) should be p + q + (P + Q)*seasonal + 1*fitmean +
)          ncols(x)*(!isscalar(x))
)
) NOTE: arimares assumes the parametrization implied by the functions
) autoreg(phi,Z) and movavg(theta,Z) which produce as output
) y[1] = Z[1], y[2] = Z[2] + phi[1]*y[1], y[3] = Z[3] + phi[1]*y[2] +
)   phi[2]*y[1] + ... (Arsign = -1)
) y[1] = Z[1], y[2] = Z[2] - theta[1]*Z[1], y[3] = Z[3] - theta[1]*Z[2] -
)   theta[2]*Z[1] + ... (Masign = -1)
) If a different convention is used, it should be handled in the controlling
) macro
) The value is
)   residual vector                             sigmahat is not provided
)   vector(sigmahat*sqrt(log(det)),residuals)   sigmahat > 0
)   sqrt(log(det))                              sigmahat < 0
)
)) Other macros used (should be loaded automatically if available)
))  acfarma
))
)) 001019 added arsign:-1 and masign:-1 in arg list of acfarma
)) 001116 Fixed bug
)) 001130 Fixed some bugs and corrected computation of determinant
) Version 001130
# usage: residuals <- $S(b,x,y,structure(pdq,PDQ,seasonal,cast,\
#        cycles,fitmean [,sigmahat]))
# If x is not scalar, it contains independent variables
@b <- $1 # starting values
@x <- $2 # matrix of exogenous variables if not a scalar
@y <- $3 # time series

@params <- $4 #structure(pdq,PDQ,seasonal,fitmean,cast,cycles[,sigmahat])

if (!isnull(@params$pdq)){
	@p <- @params$pdq[1]
	@d <- @params$pdq[2]
	@q <- @params$pdq[3]
}else{
	@p <- @d <- @q <- 0
}

@seasonal <- @params$seasonal
if (@seasonal > 0){
	@P <- @params$PDQ[1]
	@D <- @params$PDQ[2]
	@Q <- @params$PDQ[3]
}else{
	@P <- @D <- @Q <- 0
}

@fitmean <- @params$fitmean

@cast <- @params$cast # back and forecasting range

@cycles <- @params$cycles

@sigmahat <- keyvalue(@params, "sigmahat", "number")
@usedet <- !isnull(@sigmahat)
if (!@usedet){
	@sigmahat <- 1
}
@N <- if (@sigmahat < 0) {
	-@sigmahat
} else {
	length(@y)
}
delete(@params)

@nx <- if(isscalar(@x)){
	0
}else{
	ncols(@x)
}
@nobs <- @N - @d - @seasonal*@D
@n <- @nobs + @cast

@mu <- if (@fitmean){
	@b[1]
}else{
	0
}
@np <- @fitmean # keeps track of number of parameters

if(@nx > 0){
	if (@sigmahat > 0){
		@y <-- @x %*% @b[@np + run(@nx)]
	}
	@np <-+ @nx
}

@phi <- if (@p > 0){
	@b[@np + run(@p)]
}else{
	0
}
@np <-+ @p

@theta <- if (@q > 0){
	@b[@np + run(@q)]
}else{
	0
}
@np <-+ @q

if (@P > 0){
	@phiS <- rep(0,@seasonal*@P)
	@phiS[@seasonal*run(@P)] <- @b[@np + run(@P)]
	@phi <- if (@p > 0){
		@p <-+ @P*@seasonal
		-movavg(@phi,movavg(@phiS,padto(1,@p+1)))[-1]
	}else{
		@p <- @P*@seasonal
		@phiS
	}
	@np <-+ @P
	delete(@phiS)
}

if (@Q > 0){
	@thetaS <- rep(0,@seasonal*@Q)
	@thetaS[@seasonal*run(@Q)] <- @b[@np + run(@Q)]
	@theta <- if (@q > 0){
		@q <-+ @Q*@seasonal
		-movavg(@theta,movavg(@thetaS,padto(1,@q+1)))[-1]
	}else{
		@q <- @Q*@seasonal
		@thetaS
	}
	@np <-+ @Q
	delete(@thetaS)
}
delete(@P, @Q)

 #@p and @q are now orders of actual AR and MA operators in phi & theta
if (@sigmahat > 0){
	if (@d > 0){#compute differences
		for(@j,1,@d){
			@y <- movavg(1,@y)[-1]
		}
	}

	if (@D > 0){#seasonal differences
		for(@j,1,@D){
			@y <- movavg(1,@y,seasonal:@seasonal)[-run(@seasonal)]
 #			@y <-	movavg(@seasdif,@y)[-run(@seasonal)]
		}

	}
	@y <- vector(rep(0,@cast),@y - @mu)
} #if (@sigmahat > 0){}

if (@sigmahat < 0){# just computing sqrt(log(det))
	@sigmahat <- 1
	@res <- NULL
}elseif (@p + @q == 0){ # null model
	@res <- @y
}elseif (@cast > 0){ # do backcasting
	@pastlim <- vector(1,@cast)
	@past <- run(1,@cast)
	@lim2B <- vector(@cast+1,@n - @p)

	# compute backward residuals, then backcast
	@res <- if (@p > 0){
		movavg(@phi,@y,reverse:T,\
				limits:@lim2B,start:0*@y)#backres
	}else{
		@y
	}

	if (@q > 0){
		@lim2B[2] <- @lim2B[2] - @q
		@res <- autoreg(@theta,@res,reverse:T,\
			limits:@lim2B,start:@res)#backres
	}
	if (@q > 0){
		@res <- movavg(@theta,@res,reverse:T,\
			limits:@pastlim,start:@res)#backcast
	}
	if (@p > 0){
		@y <- autoreg(@phi,@res,reverse:T,\
			limits:@pastlim,start:@y)#backcast
	}else{
		@y[@past] <- @res[@past]
	}

	if (@p > 0){#forward resids
		@res <- movavg(@phi,@y)
	}else{
		@res <- @y
	}

	if(@q > 0){
		@res <- autoreg(@theta,@res)
	}#forward resids

	if(@cycles > 0){ # do alternating back and forecasting loop
		@res <- padto(@res,@n + @cast)
		@y <- padto(@y,@n + @cast)

		@lim2F <- vector(1, @n)
		@lim2B <- @lim2F + @cast

		@futlim <- @pastlim + @n #limits for forcast
		@future <- @past + @n

		for(@k,1,@cycles){#forecast/backcast loop
			# start forecast using forward residuals

			if (@q > 0){
				@res <- movavg(@theta,@res,\
					limits:@futlim,start:@res)#part of forecast of y
			}

			if (@p > 0){
				@y <- autoreg(@phi,@res,\
					limits:@futlim,start:@y)#forecast
			}else{
				@y[@future] <- @res[@future]
			}

			# Compute back residuals[-@pastlim], then back cast
			@y[@past] <- 0
			@res <- if (@p > 0){
				movavg(@phi,@y,reverse:T,limits:@lim2B,start:@y)#backres
			}else{
				@y
			}

			if (@q > 0){
				@res <- autoreg(@theta,@res,reverse:T,\
					limits:@lim2B,start:@res)#backres
			}

			# Now backcast @y[@pastlim]
			if (@q > 0){
				@res <- movavg(@theta,@res,reverse:T,\
							limits:@pastlim,start:@y)#backcast
			}

			if (@p > 0){
				@y <- autoreg(@phi,@res,reverse:T,\
					limits:@pastlim,start:@y)#backcast
				@y[@future] <- 0
			} else {
				@y <- @res
			}

			# compute forward residuals
			if (@p > 0){
				@res <- movavg(@phi,@y,limits:@lim2F,start:@y)
			}else{
				@res <- @y
			}

			if (@q > 0){
				@res <- autoreg(@theta,@res,limits:@lim2F,start:@res)
			}
			@res[@future] <- 0
		}# for(@k,run(@cycles))
		delete(@k, @futlim, @lim2F,@future)
		@res <- padto(@res, @n) # trim off excess
	}# if(@cycles > 0)

	delete(@pastlim,@lim2B,@past)
}else{
	@res <- movavg(@phi,autoreg(@theta, @y))[-run(@p + @q)]
}

if (@usedet){
	if (@p + @q > 0){
		if (!ismacro(acfarma)){
			getmacros(acfarma,quiet:T,printname:F)
		}
		@nlags <- if (@q == 0){
			@p
		}else{
			@nobs-1
		}
		@gamma <- acfarma(@phi,@theta,lag:@nlags,masign:-1,arsign:-1)
		@logdet <- (@nlags+1)*log(@gamma[1]) +\
			sum(run(@nlags,1)*log(1 - partacf(@gamma[-1]/@gamma[1])^2))

		@logdet <- max(@logdet,0)

		@res <- vector(@sigmahat*sqrt(@logdet), @res)
		delete(@logdet,@gamma,@nlags)
	}else{
		@res <- vector(0, @res)
	}
}
delete(@p,@d,@q,@fitmean,@mu,@phi,@theta,@b,@n,@cast,@y,\
	@sigmahat,@usedet)
delete(@res,return:T)
%arimares%

===> innovations <===
innovations   MACRO DOLLARS
) innovations computes the "innovation" algorithm given on p. 71 of
) Brockwell, and Davis, computing a M by M matrix containing
) coefficients and prediction variances.  The result is not affected by
) variable MASIGN if it exists
)
) Usage:
)  innovations(gamma [,lag:m] [,final:T])
)   gamma          a covariance matrix or an autocovariance
)                  function (ACVF)
)   m > 0          integer with, default m = nrows(gamma) - 1.
)
) When gamma is an ACFV, innovations(gamma [,lag:m]) is the same as
) innovations(toeplitz(gamma) [,lag:m])
)
) The macro actually uses Cholesky decomposition rather than the
) algorithm in Brockwell and Davis
)
) Value:
)    structure(theta:Theta, v:V)                  without 'final:T'
)    structure(theta:vector(Theta[m,]),v:V[m+1])  with 'final:T'
)  where Theta and V are as follows:
)
) Suppose {X[1], X[2], ..., X[m+1]} is a sequence of random variables
) with covariance matrix gamma (toeplitz(gamma) in ACVF case).
)
)   Theta is an m by m matrix with the property that
)     Theta[m,m]*Z[m] + Theta[m,m-1]*Z[m-1] + ... + Theta[m,1]*Z[1]
)   is the best one step predictor of X[m+1] based on prediction
)   errors Z[1] = X[1], Z[2] = X[2] - Theta[1,1]*Z[1], ...,
)    Z[j+1] = X[j+1] - (Theta[j,j]*Z[j]+Theta[j,j-1]*Z[j-1]+...
)        +theta[j,1]*Z[1]),..., Z[m]
)
)   V is a length m+1 vector with V[1] = gamma[1] = Var[X[1]] and
)   V[j+1] = prediction error variance when X[j+1] is predicted using
)   Z[1], ..., Z[j].
)
)) Other macros used
))   None
) Version 000908
# $S(gamma [,lag:m] [,final:T])
if ($v != 1){
	error("usage: $S(gamma [,lag:m] [,final:T])", macroname:F)
}
@kappa <- argvalue($1,"argument 1 to $S","real matrix nonmissing")
@n <- nrows(@kappa)
if (!isvector(@kappa) && ncols(@kappa) != @n){
	error("argument 1 to $S not square matrix", macroname:F)
}
if (@n == 1){
	error("nrows(arg 1) < 2")
}
@keys <- if ($k > 0){
	structure($K)
} else {
	structure(notakey:NULL)
}
@lag <- keyvalue(@keys,"lag*","positive count",default:@n-1)
if (@lag >= @n){
	error("value of 'lag' >= nrows(arg 1)")
}

@final <- keyvalue(@keys,"final","TF",default:F)

delete(@keys)

if (@lag != @n - 1){
	@n <- @lag + 1
	@J <- run(@n)
	if (isvector(@kappa)){
		@kappa <- @kappa[@J]
	}else{
		@kappa <- @kappa[@J,@J]
	}
	delete(@J)
}
if (isvector(@kappa)){
	@kappa <- toeplitz(@kappa)
}

@r <- cholesky(delete(@kappa,return:T))
@v <- diag(@r)

@theta <- ((@r/@v)' - dmat(@n, 1))[-1,-@n]
@v <- @v^2

if (@final){
	@v <- @v[@n-1]
	@theta <- vector(@theta[@n-1,])
}
delete(@n, @lag, @final, @r)
structure(theta:delete(@theta,return:T),v:delete(@v,return:T))
%innovations%

===> moveoutroots <====
moveoutroots    MACRO DOLLARS
) moveoutroots(theta) where theta is a REAL vector defining the polynomial
) P(z) = 1 - sum(theta*z^run(length(theta))), returns a new vector
) theta1 defining a new polynomial P1(z) which has the same zeros as
) P(z) except that any zero z0 with |z0| < 0 is replace by 1/z0.
)
) For example, movoutroots(2) returns 1/2,
) moveoutroots(vector(2,-1.25)) returns vector(1.6, -.8)
) because 1 - 2*z + 1.25*z^2 has zeros 8 +- .4*i
) and 1 - 1.6*z + .8*z^2) has roots (1 +- .5*i) = 1/(.8 -+ .4*i)
) Thus moveoutroots(theta) are the coefficients of an invertible MA
) operator and, when no zero has modulus 1, the coefficients of
) a stationary AR operator
))
)) Note that polyroot computes the roots of the polynomial
)) z^n - coefs[1]*z^(n-1) - ... - coefs[n-1]*z - coefs[1] whose
)) zeros are the reciprocals of the roots of P(z)
))
)) 000908 Version by C. Bingham
#$S(coefficients)
@coefs <- argvalue($1,"argument 1 to $S",vector("real","vector","nonmissing"))

@degree <- length(@coefs)
if (@degree == 1){
	if (abs(@coefs) > 1){
		@coefs <- 1/@coefs
	}
}else{
	@roots <- polyroot(@coefs)
	if (max(hypot(creal(@roots),cimag(@roots))) > 1){
		@coefs <- padto(1,@degree+1)
		@k <- 0

		while(@k < @degree){
			@k <-+ 1
			@real <- creal(@roots[@k,])
			@imag <- cimag(@roots[@k,])
			if (@imag == 0){#real root
				if (@real > 1){
					@real <- 1/@real
				}
			}else{#pair of complex roots
				@k <-+ 1
				@modsq <- hypot(@real,@imag)^2
				if (@modsq > 1){
					@modsq <- 1/@modsq
					@real <-* @modsq
				}
				@real <- vector(2*@real,-@modsq)
			}
			@coefs <- movavg(@real,@coefs)
		}# while(@k < @degree)
		delete(@real,@imag,@modsq,@k,silent:T)
		@coefs <- -@coefs[-1]
	}# if (max(creal(cpolar(@roots))) > 1)
	delete(@roots)
}
delete(@degree)
delete(@coefs,return:T)
%moveoutroots%

===> acfarma <===
acfarma       MACRO DOLLARS
) Macro to compute the spectrum of a stationary ARMA time series
) Usage:
)  acfarma(phi,theta [, L] [,nfreq:Nfreq]\
)            [,arsign:Arsign] [,masign:Masign]),
)   phi       non MISSING REAL vector of AR coefficients
)   theta     non MISSING REAL vector of MA coefficients
)   L         positive integer, the number of lags
)   Nfreq     integer != 0 with no prime factors > 29
)   Arsign    +1 or -1 specifies sign convention for phi
)   Masign    +1 or -1 specifies sign convention for theta
) To omit part of model, use phi = 0 or theta = 0
)
) The default for Arsign is ARSIGN if it exists or -1 otherwise.
) The default for Masign is MASIGN if it exists or -1 otherwise.
)
) The ARMA model is defined as
)   (1 + sum(Arsign*phi*B^run(p)))X[t] =
)                 (1 + sum(Masign*theta*B^run(q)))Z[t].
) where B is the backshift operator and {Z[t]} is 0 mean white noise with
) variance = 1.
)
) The ACFV is computed using a discrete Fourier transform of length
) abs(Nfreq).  When Nfreq < 0, it is not checked for large prime factors.
) The default for Nfreq is goodfactors(max(500,L+1)), ensuring that it
) has no prime factors > 29
)
) For backward compatibility, you can use 'S' in place of keyword
) 'nfreq'
)
)) Other macros used
))   None
))
)) 001227 modified so can be used as acfarma(phi,theta,L ...) which
))        is equivalent to acfarma(phi,theta,lag:L ...)
)) 001231 fixed bug; nfreq now at least L+1
)) 010103 it is an error if Nfreq has a prime factor > 29
))        S is no longer used as a default for Nfreq
))
#$S(phi,theta [,lag:l] [,nfreq:nfreq] [,arsign:1] [,masign:1]),
#  phi, theta REAL vectors, integer l > 0
if ($v < 2 || $v > 3){
	error("usage: $S(phi,theta [,lag:nlags] [,nfreq:nfreq] [,arsign:1] [,masign:1])",\
		macroname:F)
}

@phi <- argvalue($1,"phi","nonmissing real vector")
if(alltrue(isscalar(@phi),@phi[1] == 0)){
	@phi <- NULL
} elseif (max(creal(cpolar(polyroot(@phi)))) >= 1) {
	error("AR coefficients do not define stationary operator")
}

@theta <- argvalue($2,"theta","nonmissing real vector")
if(alltrue(isscalar(@theta),@theta[1] == 0)){
	@theta <- NULL
}

@keys <- if ($k > 0){
	structure($K)
} else {
	structure(notakey:NULL)
}
@lag1 <- if ($v == 3){
	argvalue($3,"lag","positive count")
} else {
	NULL
}
@lag <- keyvalue($K,"lag*","positive count")

if (isnull(@lag) && isnull(@lag1)){
	@lag <- 100
}
if (alltrue(!isnull(@lag),!isnull(@lag1),@lag1 != @lag)){
	error("conflicting specification of number of lags")
}
if (isnull(@lag)){
	@lag <- @lag1
}

@S <- keyvalue($K, "S", "positive count")
@nfreq <- keyvalue($K,"nfreq*","integer scalar")
if (!isnull(@S) && !isnull(@nfreq)){
	if (@S != @nfreq) {
		error("'nfreq' and 'S' both used but with different values")
	}
} elseif(!isnull(@S)) {
	@nfreq <- @S
}
delete(@S)
if (!isnull(@nfreq)){
	if (@nfreq < 0) { # no checking done
		@nfreq <- abs(@nfreq)
	} elseif (goodfactors(@nfreq) != @nfreq){
		error(paste("No. of frequencies",@nfreq,\
			"has prime factor > 29; try",goodfactors(@nfreq)))
	}
} else {
	# no longer using S for default
	@nfreq <- goodfactors(max(500, @lag + 1))
}

@sign <- if (alltrue(isscalar(MASIGN, real:T),!anymissing(MASIGN))) {
	MASIGN
} else {
	-1
}
@masign <- keyvalue(@keys,"masign","number", default:@sign)
if (abs(@masign) != 1) {
	error("MASIGN or value for 'masign' != +- 1")
}
@sign <- if (alltrue(isscalar(ARSIGN, real:T),!anymissing(ARSIGN))) {
	ARSIGN
} else {
	-1
}
@arsign <- keyvalue(@keys,"arsign","number", default:@sign)
if (abs(@arsign) != 1) {
	error("ARSIGN or value for 'arsign' != +- 1")
}
delete(@sign,@keys)

@top <- if(isnull(@theta)){
	1
}else{
	@theta <-* -@masign
	hreal(hprdhj(rft(movavg(@theta,padto(1,@nfreq)))))
}

@bottom <- if(isnull(@phi)){
	1
}else{
	@phi <-* -@arsign
	hreal(hprdhj(rft(autoreg(@phi,padto(1,@nfreq)))))
}

@result <- @top*@bottom

@result <- if (isscalar(@result)){
	padto(@result,@lag+1)
}else{
	rft(@result,divbyt:T)[run(@lag+1)]
}

delete(@phi,@theta,@nfreq,@lag,@top,@bottom,@masign,@arsign)
delete(@result,return:T)
%acfarma%

===> specarma <===
specarma      MACRO DOLLARS
) Macro to compute the spectrum (in Real form) of an ARMA time series
) Usage:
)   specarma(phi ,theta [, nfreq:Nfreq] [,arsign:Arsign] [,masign:Masign])
)    phi       non MISSING REAL vector of AR coefficients
)    theta     non MISSING REAL vector of MA coefficients
)    Nfreq     positive integer scalar, number of equally spaced
)              frequencies, having no prime factors > 29
)    Arsign    +1 or -1 specifies sign convention for phi
)    Masign    +1 or -1 specifies sign convention for theta
)    To omit part of model, use phi = 0 or theta = 0
)
) The ARMA model is defined as
)   (1 + sum(Arsign*phi*B^run(p)))X[t] =
)                 (1 + sum(Masign*theta*B^run(q)))Z[t].
) where B is the backshift operator and {Z[t]} is 0 mean white noise with
) variance = 1.
)
) The result is a vector of length Nfreq containing the spectrum at
) frequencies 0, 1/Nfreq, 2/Nfreq, ..., (Nfreq-1)/Nfreq cycles per unit
) time.
)
) It is an error if nfreq:Nfreq is an argument and Nfreq has a prime
) factor > 29.
) The default value for Nfreq is S when S exists and is a
) positive integer, or 400 if not.  It is an error if S has a prime
) factor > 29.
)
) The default for Arsign is ARSIGN if it exists or -1 otherwise.
) The default for Masign is MASIGN if it exists or -1 otherwise.
)   Arsign = -1 and Masign = -1 correspond to the convention used by
)   Box and Jenkins.
)   Arsign = +1 and Masign = +1 correspond to the convention used by
)   Brockwell and Davis
)
)) Other macros used
))   None
)
) Version 010102 Nfreq must have no prime factors > 29
))
#$S(phi ,theta [, nfreq:nf] [,arsign:Arsign] [,masign:Masign])

if ($v != 2){
	error("usage: $S(phi ,theta [,nfreq:nf])", macroname:F)
}

@phi <- argvalue($1,"phi","nonmissing real vector")
if(alltrue(isscalar(@phi),@phi[1] == 0)){
	@phi <- NULL
} elseif (max(creal(cpolar(polyroot(@phi)))) >= 1) {
	print("WARNING: AR coefficients do not define stationary model",\
		macroname:T)
}

@theta <- argvalue($2,"phi","nonmissing real vector")
if(alltrue(isscalar(@theta),@theta[1] == 0)){
	@theta <- NULL
}

@keys <- if ($k > 0) {
	structure($K)
} else {
	structure(notakey:NULL)
}
@nfreq <- keyvalue(@keys,"nfreq","positive count")
@S <- keyvalue(@keys,"S","positive count")
if (!isnull(@S) && !isnull(@nfreq)){
	if (@nfreq != @S){
		error("values of 'S' and 'nfreq' differ; don't know which to use")
	}
} elseif (isnull(@nfreq)){
	@nfreq <- @S
}
delete(@S)
if(isnull(@nfreq)){
	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 <- 400
	}
} elseif (@nfreq != goodfactors(@nfreq)) {
	error(paste("No. of frequencies",@nfreq,\
		"has a prime factor > 29; try",goodfactors(@nfreq)))
}

@sign <- if (alltrue(isscalar(MASIGN, real:T),!anymissing(MASIGN))) {
	MASIGN
} else {
	-1
}
@masign <- keyvalue(@keys,"masign","number", default:@sign)
if (abs(@masign) != 1) {
	error("MASIGN or value for 'masign' != +- 1")
}
@sign <- if (alltrue(isscalar(ARSIGN, real:T),!anymissing(ARSIGN))) {
	ARSIGN
} else {
	-1
}
@arsign <- keyvalue(@keys,"arsign","number", default:@sign)
if (abs(@arsign) != 1) {
	error("ARSIGN or value for 'arsign' != +- 1")
}
delete(@sign,@keys)

@top <- if(isnull(@theta)){
	1
}else{
	@theta <-* -@masign
	hreal(hprdhj(rft(movavg(@theta,padto(1,@nfreq)))))
}

@bottom <- if(isnull(@phi)){
	1
}else{
	@phi <-* -@arsign
	hreal(hprdhj(rft(movavg(@phi,padto(1,@nfreq)))))
}
delete(@arsign,@masign)
@result <- if(isnull(@phi) && isnull(@theta)){
	rep(1,@nfreq)
}else{
	@top/@bottom
}
delete(@phi,@theta,@nfreq,@top,@bottom)
delete(@result,return:T)
%specarma%

===> rhatvar <===
rhatvar       MACRO DOLLARS
) Macro to compute variance of sqrt(n)*rhohat(k), using Bartlett's
) formulat
) Usage:
)  v <- rhatvar(rho [,lag:L])
)   rho     REAL vector of autocorrelations, starting with either lag 0
)           (rho[1] == 1) or lag 1 (rho[1] != 1)
)   L       positive integer
)
) The value returned is v = n*vector(var(rhohat(1)),var(rhohat(2)),...,
)   var(rhohat(L)))
) Default value of L is (maximum lag for rho)/2
) Version 001231
))
# usage: $S(rho [,lag:L]), rho REAL vector, L positive integer

@rho <- argvalue($1,"rho","vector real nonmissing")
if(@rho[1] != 1){
	@rho <- vector(1,@rho)
}

if(max(abs(@rho[-1])) >= 1){
	error("at least one element of rho is not between -1 and +1")
}

@m <- length(@rho) - 1

@lag <- keyvalue($K,"lag*","positive count",default:ceiling(@m/2))

@rho1 <- vector(rep(0,2*@m),reverse(@rho[-1]),@rho,rep(0,2*@m))

@var <- rep(0,@lag)
@J <- 3*@m + 1 + run(@m)
for(@i, 1, @lag){
	@var[@i] <-\
		sum((@rho1[@J + @i] + @rho1[@J - @i] - \
			2*@rho[@i+1]*@rho1[@J])^2)
}
delete(@rho,@rho1,@i,@m,@J,@lag)
delete(@var,return:T)
%rhatvar%

===> rhatcovar <===
rhatcovar     MACRO DOLLARS
) Macro using Bartlett's formulat to compute n*covar{rhohat(i),rhohat(j)}
) or the covariance matrix of sqrt(n)*vector(rhohat(1),...,rhohat(L))
) Usage:
)  cov <- rhatcovar(rho,i, j)
)  V <- rhatcovar(rho, lag:L)
)   rho     REAL vector of autocorrelations starting with
)           either lag 0 (rho[1] == 1) or lag 1 (rho[1] != 1)
)   i, j    Positive integers specifying lags
)   L       Positive integer specifying maximum lag
)
) When lag:L is not an argument, cov = n*covar{rhohat(i),rhohat(j)}
)
) When lag:L is an argument, V is L by L and
)   V = Cov[sqrt(n)*rhat(1),..., sqrt(n)*rhat(L)]
)
) Version 000909
))
# usage: $S(rho,i, j) or $S(rho [,lag:L])

@usage <- "usage: $S(rho,i,j) or $S(rho [,lag:L])"
if ($v != 1 && $v != 3 || $k > 1){
	error(@usage)
}

@rho <- argvalue($01,"rho",vector("vector","real","nonmissing"))
if(@rho[1] != 1){
	@rho <- vector(1,@rho)
}
if(max(abs(@rho[-1])) >= 1){
	error("at least one element of rho is not between -1 and +1")
}

@single <- $v == 3

if (@single){
	@i <- argvalue($02,"argument 2","positive count")
	@j <- argvalue($03,"argument 3","positive count")
}
@m <- length(@rho) - 1

@lag <- keyvalue($K,"lag*","positive count")

if (@single){
	if (!isnull(@lag)){
		error(@usage)
	}
}else{
	if (isnull(@lag)){
		@lag <- ceiling(@m/2)
	}
}

if (!@single){
	@ilim <- if (@lag > 1) {
		run(@lag)
	}else{
		rep(1,2)
	}
	@j1 <- 1
}else{
	@ilim <- rep(max(@i,@j),2)
	@j1 <- min(@i,@j)
}

@result <- if(!@single){
	matrix(rep(0,@lag^2), @lag)
}else{
	0
}

@zeros <- rep(0,2*@m)
@rho1 <- vector(@zeros,reverse(@rho[-1]),@rho,@zeros)

delete(@lag,@zeros)

@J <- 3*@m+1+run(@m)
for(@i,@ilim){
	@left <- @rho1[@J + @i] + @rho1[@J - @i] -\
			2*@rho[@i+1]*@rho1[@J]
	for(@j,run(@j1,@i)){
		@right <- if (@i == @j){
			@left
		}else{
			@rho1[@J + @j] + @rho1[@J - @j] -\
			2*@rho[@j+1]*@rho1[@J]
		}
		@tmp <- sum(@right*@left)
		if (!@single){
			@result[@i,@j] <- @result[@j,@i] <- @tmp
		}else{
			@result<- @tmp
			break 2
		}
	}
}
delete(@rho,@rho1,@i,@j,@m,@left,@right,@j1,\
	@single, @ilim,@tmp,@J,silent:T)
delete(@result,return:T)
%rhatcovar%

===> _polish <===
_polish     MACRO DOLLARS
) Macro for carrying one or more "polishing" cycles as described
) on p. 155-156 of Brockwell & Davis.  This essentially computes
) approximated derivatives of the conditional sum of squares
) as described in Sec. 7.2.4 of Box and Jenkins.
) Usage:
)   _polish(x,structure(phi:phi,theta:theta), p, q, cycles)
)   where x is assumed to be differenced and detrended
)   cycles > 0 is the number of updating cycles
)  Return value: structure(phi,theta)
)) Other macros used
))   None
) Version 990206
# $S(x,structure(phi:phi,theta:theta),p,q,cycles)
@x <- $1
@b <- $2
@p <- $3
@q <- $4
@cycles <- $5

@phi <- @b$phi
@theta <- @b$theta
delete(@b)

@nobs <- nrows(@x)
@maxpq <- max(@p,@q)
@npar <- @p + @q
@limits <- vector(@maxpq + 1, @nobs)

@n1 <- @nobs - @maxpq
@J <- run(@n1) + @maxpq
@X <- matrix(rep(0,@n1*(@npar+1)),@n1)

for(@l,1,@cycles){
	# compute residuals
	if (@p > 0){
		@z <- movavg(@phi,@x)
	}else{
		@z <- @x
	}

	@z[run(@maxpq)] <- 0
	if (@q > 0){
		@z <- autoreg(@theta,@z,limits:@limits, start:@z)
	}

	# compute approximate derivatives
	if (@p > 0){
		@v <- autoreg(@phi,@z,limits:@limits,start:@z)
	}
	if (@q > 0){
		@w <- autoreg(@theta,@z,limits:@limits,start:@z)
	}

	@k <- 0
	if (@p > 0){
		for(@i,1,@p){
			@k <-+ 1
			@X[,@k] <- -@v[@J-@i]
		}
		delete(@v)
	}
	if (@q > 0){
		for(@i,1,@q){
			@k <-+ 1
			@X[,@k] <- @w[@J-@i]
		}
		delete(@w);
	}

	@k <-+ 1
	@X[,@k] <- @z[@J]

	@XX <- swp(@X %c% @X, run(@npar))
	@dphi <- if (@p > 0){
		vector(@XX[@k,run(@p)])
	}else{
		?
	}
	@dtheta <- if (@q > 0){
		vector(@XX[@k,@p + run(@q)])
	}else{
		?
	}

	if (@q > 0){
		@theta <-- @dtheta
	}
	if (@p > 0){
		@phi <-- @dphi
	}
}# for(@l,1,@cycles)

delete(@X,@maxpq,@npar,@l,@x,@p,@q,@k,@dphi, @dtheta)
structure(phi:delete(@phi,return:T),\
	theta:delete(@theta,return:T),\
	xtxswept:delete(@XX,return:T))
%_polish%

===> arimatest <===
arimatest MACRO
) This macro is in file Arima.mac
) It is here to make it easier to test which macro file is
) being read.
print("$S")
%arimatest%

_E_N_D_O_F_M_A_C_R_O_S_ Marker for end of macros and data
Help for macros in arima.mac
(C) 1999, 2000, 2001, 2002, 2003 by Gary W. Oehlert and Christopher
Bingham
Updated 030814 CB

!!!! Starting marker for message of the day
!!!! Ending marker for message of the day

???? Starting marker for list of up to 32 comma/newline separated keys
ARIMA models
Autocovariance
Complex numbers
Frequency domain
General
Nonlinear fitting
Preliminary estimation
Spectrum analysis
Time domain
???? Ending marker for keys

====acfarma()#arima models,time domain,autocovariance
%%%%
acfarma(phi,theta [,lag:L][,nfreq:Nfreq][,arsign:Arsign,\
  masign:Masign]), REAL vectors phi, theta, positive integer L, integer
  Nfreq != 0, Arsign and Masign either +1 or -1
%%%%
@@@@usage#Usage
Macro acfarma() computes the theoretical autocovariance function (ACVF)
of an ARMA time series with given coefficients and innovation variance
1.

gamma <- acfarma(phi, theta, lag:L), where phi and theta are REAL
vectors and L > 0 is an integer, returns the ACVF from 0 to L lags of
an ARMA time series with REAL coefficient vectors phi and theta.
gamma[1] is the variance and gamma[h+1] is the lag h autocovariance, h
= 1, ..., L.  The innovation variance is assumed to be 1.

To omit part of model, use phi = 0 or theta = 0.

'lag:L' can be omitted, in which case the default value is L = 100.

You can make an "impulse" plot of the first 50 lags of the ACFV by
  Cmd> tsplot(acfarma(phi, theta, lag:50),0,impulse:T,lines:F)

If you omit 'lines:F', lines are drawn between successive values.

@@@@sign_conventions#Sign conventions
The model assumed for series X is
  (1 + Arsign*sum(phi*B^run(p)))X[t] =
                    (1 + Masign*sum(theta*B^run(q)))Z[t].

where Arsign and Masign are either +1 or -1 and {Z[t]} is zero mean
white noise with standard deviation 1.

The default value for Arsign is variable ARSIGN if it exists and -1
otherwise.  The default value for Masign is variable MASIGN if it exists
and -1 otherwise.  Or you can include either or both keyword phrases
'arsign:Arsign' and 'masign:Masign' as arguments.  See topic 'MASIGN'
for an explanation.  If you want the convention used by Brockwell and
Davis, you should do the following before you start your analysis

  Cmd> MASIGN <- 1; ARSIGN <- -1

It is an error for ARSIGN, MASIGN or supplied values of Arsign and
Masign to be other than +1 or -1.

@@@@nfreq_keyword#Keyword nfreq
acfarma() uses the fast discrete Fourier transform to compute the ACVF.
The number of frequencies used is the smallest integer Nfreq >=
max(500,L+1) with no prime factors > 29.  Or you can specify Nfreq
by including nfreq:Nfreq as an argument, where Nfreq is non-zero
integer.  When Nfreq < 0, abs(Nfreq) is used without further checking;
when Nfreq > 0 it is an error when Nfreq has a prime factor > 29.

@@@@see_also#Cross references
See also specarma(), tsplot().
@@@@______

====arima()#arima models,time domain,nonlinear fitting
%%%%
arima(y [,pdq:vector(p,d,q)] [,PDQ:vector(P,D,Q),seasonal:period]\
  [,x:x,fitmean:T, start:b0,active:active,cast:n, cycles:m,mle:T,\
  mlecycles:m1, masign:Masign, arsign:Arsign, maxit:itmax, minit:itmin,\
  crit:vector(numsig, nsigsq, delta), print:T,keep:T,quiet:T]),
  REAL vector y, REAL variable x, nonnegative integers p, d, q, P, D, Q,
  period > 1, cast, n, m, m1, REAL vector b0, LOGICAL vector active with
  length(active) = length(b0), Masign and Arsign +1 or -1, integers
  itmax, itmin, numsig, nsiqsq, REAL scalar delta >= 0
%%%%
@@@@usage#Usage
arima(y [,pdq:vector(p,d,q)] [,PDQ:vector(P,D,Q),seasonal:Period]),
where y is a REAL vector with no MISSING values and p, d, q, P, D, Q
are non-negative integers and Period > 1 is an integer, uses
unconditional least squares to estimate the parameters of an
ARIMA(p,d,q)x(P,D,Q) time series model, with season length Period.

When d = D = 0, a non-zero mean is fit; otherwise no mean is fit.
Keyword 'seasonal' is required when keyword 'PDQ' is used.

arima(y, x:x [,pdq:vector(p,d,q)] [,PDQ:vector(P,D,Q),seasonal:Period])
does the same, except a linear regression with ARIMA errors of y on the
columns of REAL matrix x is carried out.  x must have no MISSING
elements and nrows(x) = nrows(y).  It should not include a constant
column.

arima(y [, x:x], pdq:pdq, PDQ:PDQ, mle:T) does the same, except maximum
likelihood estimation is used instead of unconditional least squares.

You can use keywords 'minit', 'maxit', 'cast', 'cycles', and 'mlecycles'
to control certain details of the estimation algorithm.  See below for
more information.

You can use keywords 'masign' and 'arsign' or define variables MASIGN
and ARSIGN to control the sign conventions assumed for moving average
and autoregressive coefficients.  See below and topic 'MASIGN' for
details.

@@@@model_coefficients#Model coefficients
The model coefficients estimated include some or all of the following:
   mu     mean (of differenced data when d > 0 or D > 0)
   beta   vector of coefficients of columns of x, if any
   phi    vector of AR coefficients
   theta  vector of MA coefficients
   phiS   vector of seasonal AR coefficients
   thetaS vector of seasonal MA coefficients.

They are grouped in a coefficient vector b = vector(mu, beta, phi,
theta, phiS, thetaS), omitting any coefficients not in the model.

You can provide starting values using keyword 'start'; see below.

To force inclusion of mu when d > 0 or D > 0, use 'fitmean:T' as an
argument.  To force exclusion of mu when d = D = 0, use 'fitmean:T'.
See below.

You can omit certain coefficients from being included in the estimation
algorithm using keyword 'active'; see below.  This is useful when
fitting models using a subset of lags.  Any starting values you provide
for them (non-zero starting value for mu) using 'start' are not changed
but will be used in computing residuals.

When mu is in the model, but is not an active parameter, it is estimated
by the sample mean of the possibly differenced data unless a non-zero
starting value is given.

@@@@output#Output
arima() prints a table of the estimated coefficients, their approximate
standard errors, t = coef/StdErr, and a nominal P-value based on the t
distribution.  It also prints the mean square error (MSE), its degrees
of freedom (DF), -2*log(L), where L = likelihood, a modification of
Akaike's information critirion (AICC), and the complete sum of squared
residuals, including backcast values, if any (RSS).  Printing can be
suppressed by either quiet:T or keep:T; see below.  Side effect
variables are always produced; see below.

DF = n - d - D*Period - (number of parameters being estimated).  When mu
is estimated by a sample mean, it counts as a parameter, even if it is
not active in the optimization (active[1] is false).

@@@@side_effect_variables#Side effect variables
arima() creates the following side effect variables
  COEF = bhat = vector(muhat,betahat,phihat,thetahat,phiShat,thetaShat),
    omitting any coefficients not in the model.  When coefficient b[i]
    is inactive (seek keyword 'active' below), COEF[i] = the starting
    value, if any, or 0.
  ALLRESIDUALS = residuals from fitted model including incomplete and
    backcast residuals
  RESIDUALS = ALLRESIDUALS, omitting backcast and incomplete residuals
  RSS = sum(ALLRESIDUALS^2)
  NEG2LOGL = -2*log(likelihood) assuming a Gaussian series.  The likeli-
    hood includes a factor of (2*PI)^((n-d-D*Period)/2)
  NPAR = number of active parameters plus 1 for the mean if estimated
    as sample mean.  The active parameters correspond to rows k of
    XTXINV with XTXINV[k,k] != 0
  JACOBIAN = REAL matrix of derivatives of ALLRESIDUALS with respect to
   active parameters.
  XTXINV = analogue of solve(X' %*% X) matrix in regression computed
   from JACOBIAN, with rows and columns corresponding to inactive
   parameters set to 0
  HII, REAL vector of leverages of length n - d - D*Period, computed
   from JACOBIAN
  GRADIENT = REAL gradient vector (approximation to to partial
   derivatives sum of squares or log likelihood with respect to the
   coefficients) with an element for each active parameter

The innovation variances is estimated as MSE = RSS/(n-d-D*Period-NPAR).
MSE is used to compute estimated standard errors as
sqrt(MSE*diag(XTXINV))

@@@@choosing_starting_values#Choosing starting values
When there is no seasonal part to the model and no omitted lags, you may
be able to speed convergence by finding starting values using macros
hannriss() or innovest() and providing them to arima() using keyword
'start' (see below).

@@@@difference_from_other_programs#Difference from other programs
Other ARIMA estimation programs may differ in what is computed as the
estimated innovation variance.  Two alternatives to MSE as computed by
arima() are sum^2(ALLRESIDUALS)/DF and
sum(ALLRESIDUALS^2)/(n-d-D*Period), the maximum likelihood estimate.

Programs also may differ in sign conventions used for the autoregressive
and moving average coefficients.  In arima(), the signs are determined
either by variables MASIGN and ARSIGN (defaults are -1 and -1) or the
value of keywords 'arsign' and 'masign'.  See topic 'MASIGN'.

@@@@sign_conventions#Sign conventions
arima() allows you to choose among several conventions that have been
used for the signs of AR and MA coefficients.  These are specified by
numbers Arsign and Masign with values +1 or -1.  Either or both can be
set as values of keywords 'arsign' and 'masign'; see below.  If not set
by 'arsign', the default value for Arsign is the value of variable
ARSIGN if it exists, or -1 otherwise.  Similarly the default value for
Masign is MASIGN or -1.  For details, see topic 'MASIGN'.

To use the Brockwell and Davis convention, before you start your
analysis you should do the following:
  Cmd> MASIGN <- 1; ARSIGN <- -1

See topic 'MASIGN'.
@@@@______

                  Other Optional Keyword Phrases
  Keyword phrase               Explanation
  ---------------------------------------------------------------------
@@@@fitmean_keyword
  fitmean:T or F Include mu in model when T, omit mu otherwise. Default
                 is T with d = D = 0 and F otherwise.  When d > 0 or D
                 > 0, mu is the mean of the differenced series
@@@@active_keyword
  active:active  LOGICAL vector with length(active) = npars = length(b)
                 where b = vector(mu, beta, phi, theta, phiS, thetaS).
                 When active[i] is F, b[i] is "inactive" and does not
                 change from its starting value.  If you want to
                 estimate mu by the sample mean rather than least
                 squares or MLE, use active[1] = F and b0[1] = 0
@@@@start_keyword
  start:b0       REAL vector b0 = vector(mu0, beta0, phi0, theta0,
                 phiS0, thetaS0) of starting values values for the
                 iteration.  Default is rep(0,npar).  b0 should include
                 values for "inactive" parameters.  When the mean is in
                 the model, b0[1] = 0 is equivalent to b0[1] = sample
                 mean.
@@@@cast_and_cycles_keywords
  cast:nback     Residuals will be computed by "backcasting" for nback
                 >= 0 time points before the start of the series
  cycles:m       m (non-negative integer) cycles of forecasting/
                 backcasting cycles used in computing residuals; default
                 m = 1.
@@@@masign_and_arsign_keywords
  masign:Masign  Alters definition of MA parameters.  Masign must be +1
                 or -1.  As an example, a MA(2) model is
                  Y(t) = Z(t) + Masign*(theta[1]*Z(t-1) +
                       theta[2]*Z(t-2))
                 instead of
                  Y(t) = Z(t) - theta[1]*Z(t-1) - theta[2]*Z(t-2)
  arsign:Arsign  Alters definition of AR parameters.  Arsign must be +1
                 or -1.  As an example, an AR(2) model is
                  Y(t) = Z(t) - Arsign*phi[1]*Y(t-1) -
                       Arsign*phi[2]*Y(t-2)
                 instead of
                  Y(t) = Z(t) + phi[1]*Y(t-1) + phi[2]*Y(t-2)
@@@@maxit_and_minit_keywords
  maxit:itmax    Maximum number of iterations allowed (default 30, 0 is
                 ok).  With mle:T, when itmax < 0, no iterations of
                 least squares fitting is done and up to abs(itmax)
                 iterations of MLE using sigmahat based on coefficients
                 in b0.  When itmax = 0, no iteration is done at all,
                 but calculations are done using starting values.
  minit:itmin    Minimum number of iterations carried out (default = 0)
@@@@crit_keyword
  crit:vector(numsig, nsigsq, delta)
                 Integer numsig >= 1 is the number of accurate digits
                 wanted in coefficients.
                 Integer nsiqsq >= 1 is the number of accurate digits
                 wanted in the residual sum of squares.
                 When delta > 0, iterations stop when the norm of the
                 gradient <= delta
                 Iteration terminates the first time any of these
                 criteria is met.  When any of numsig, nsigsq or delta
                 are 0, that criterion is ignored.
                 Example: crit:vector(6,0,0) specifies iteration
                 continues until the relative change in all
                 coefficients is less than 1e-6.
@@@@print_keep_and_quiet_keywords
  print:T        Partial results are printed at each iteration
  keep:T         arima() returns a structure as value (see below).
  quiet:F        No summary results are printed, although side effect
                 variables are created.  Default is F except with keep:T

@@@@value_returned#Value returned
Without keyword phrase 'keep:T', arima() returns NULL as value.

With keyword phrase 'keep:T', arima() returns a structure with the
following components:
  Name           Contents
  coefs          vector(muhat,betahat,phihat,thetahat,phiShat,thetaShat)
                 with coefficients not in the model omitted.
  hessian        NPAR by NPAR matrix JACOBIAN %c% JACOBIAN
  gradient       Same as GRADIANT
  residuals      Same as ALLRESIDUALS
  nobs           length(y) - d - D*Period
  npar           number of active parameters
  pdq            zero padded value of keyword 'pdq' or rep(0,3)
  PDQ            zero padded value of keyword 'PDQ' or rep(0,3)
  seasonal       value of keyword 'seasonal' or 0
  active         logical vector the same length as coefs; active[i] is
                 True if and only if coefs[i] is an active parameter
  iter           Number of iterations to convergence or termination
  iconv          0: not converged by any criterion
                 1: converged by relative coefficient change criterion
                 2: converged by relative RSS or log L change criterion
                 3: changed by norm of gradient criterion
                 4: could not further reduce RSS or -2logL
  rss            Same as RSS
  neg2logL       Same as NEG2LOGL
  aicc           Same as AICC

@@@@see_also#Cross references
See also hannriss(), innovest().
@@@@______

====arimahelp()#general
%%%%
arimahelp(topic1 [, topic2 ...] [,usage:T] [,scrollback:T])
arimahelp(topic, subtopic:Subtopics), CHARACTER scalar or vector
  Subtopics
arimahelp(topic1:Subtopics1 [,topic2:Subtopics2 ...])
arimahelp(key:Key), CHARACTER scalar Key
arimahelp(index:T [,scrollback:T])
%%%%
@@@@usage#Usage
arimahelp(Topic1 [, Topic2, ...]) prints help on topics Topic1, Topic2,
... related to macros in file arima.mac.  The help is taken from file
arima.mac.

arimahelp(Topic1 [, Topic2, ...] , usage:T) prints usage information
related to these macros.

arimahelp(index:T) or simply arimahelp() prints an index of the topics
available using arimahelp().

arimahelp(Topic, subtopic:Subtopic), where Subtopic is a CHARACTER
scalar or vector, prints subtopics of topic Topic.  With subtopic:"?", a
list of subtopics is printed.

arimahelp(Topic1:Subtopics1 [,Topic2:Subtopics2], ...), where Suptopics1
and Subtopics2 are CHARACTER scalars or vectors, prints the specified
subtopics.  You can't use any other keywords with this usage.

In all the first 4 of these usages, you can also include help() keyword
phrase 'scrollback:T' as an argument to arimahelp().  In windowed
versions, this directs the output/command window will be automatically
scrolled back to the start of the help output.

arimahelp(key:key) where key is a quoted string or CHARACTER scalar
lists all topics cross referenced under Key.  arimahelp(key:"?") prints
a list of available cross reference keys for topics in the file.
@@@@______

arimahelp() is implemented as a predefined macro.

@@@@see_also#Cross reference
See help() for information on direct use of help() to retrieve
information from arima.mac.
@@@@______

====arimares()#arima models,time domain,nonlinear fitting
%%%%
residuals <- arimares(b,x,y,params), REAL vectors b, y, REAL
  vector or matrix x, params = structure(pdq:vector(p,d,q),
  PDQ:vector(P,D,Q),seasonal:n1,cast:n2,cycles:n3,fitmean:T or F
  [,sigmahat:s]) integers >= 0 p, d, q, P, D, Q, n1, n2, n3,
  integer seasonal > 0, REAL scalar sigmahat
%%%%
@@@@usage#Usage
residuals <- arimares(b,x,y,params) returns a REAL vector of residual
from a (p,d,q)x(P,D,Q) ARIMA model with specified coefficients.  The
model may optionally have one or more linear predictors.  Depending on
arguments it can do one or more forecast/backcast cycles to estimate
past residuals.  It is intended for use by other macros.

b = vector(mu,beta,phi,theta,phiS,thetaS) is a REAL vector of
coefficients.  mu is the mean, beta is a vector of coefficients of
linear predictors in the columns of x, phi, theta, phiS and thetaS are
vectors of AR, MA, seasonal AR and seasonal MA coefficients,
respectively.  Any of mu, beta, ..., thetaS are NULL if they are not in
the model.

x is either 0 (no linear predictors) or a nrows(y) by k REAL matrix of
linear predictors.

REAL vector y contains the time series being analyzed.

@@@@params_argument#Argument params
params = structure(pdq:vector(p,d,q),PDQ:vector(P,D,Q),seasonal:n1,
cast:n2,cycles:m3,fitmean:T or F [,sigmahat:s]) defines the form of the
model and details about the algorithm used.

@@@@value_returned#Value returned
The value is as follows
  residual vector                             sigmahat is not provided
  vector(sigmahat*sqrt(log(det)),residuals)   sigmahat > 0
  sqrt(log(det))                              sigmahat < 0

@@@@arima_model#ARIMA model
The ARIMA model is Phi(B)*PhiS(B)*((1 - B)^d*(1 - B)^D*Y[t] - mu) =
mu + Theta(B)*ThetaS(B)*Z(t), where B is the backshift
operator.  Here
  Phi(z) = 1 - phi[1]*z - phi[2]*z^2 - ... - phi[p]*z^p.
  Theta(z) = 1 - theta[1]*z - theta[2]*z^2 - ... - theta[q]*z^q
  PhiS(z) = 1 - phiS[1]*z^n1 - phiS[2]*z^(2*n1) - ... - phi[P]*z^(P*n1)
  ThetaS(z) = 1 - thetaS[1]*z^n1 - thetaS[2]*z^(2*n1) - ... -
     theta[Q]*z^(Q*n1)

When there are k linear predictors, there is also an term X %*% beta on
the right, where X is an nrows(y) by k matrix matrix of linear
predictors with coefficient vector beta

When another convention on the signs of the MA and/or AR coefficients is
wanted, the calling macro must make the adjustment.  That is, in the
notation used in describing arima, arimares() assumes Arsign = -1,
Masign = -1.

@@@@params_details#Details about argument params
The components of params are as follows:
 p,d,q   integers >= 0, AR order, difference, MA order
 P,D,Q   integers >= 0, seasonal AR order, seasonal difference,
         seasonal MA order
 n1      integer > 1, length of seasonal
 n2      integer >= - specifies how far to backcast
 n3      integer >= 0 number of cycles of fore/backcasting
 fitmean:T => mean will be fit
 s       REAL scalar, estimate of residual stddev.  Its presence
         signals that log(det(covariance matrix)) is to be
         computed. If sigmahat < 0, this is all that arimares()
         does, returning the scalar sqrt(log(det)).
length(b) should be p + q + (P + Q)*seasonal + 1*fitmean +
         ncols(x)*isscalar(x)
@@@@______

====arima_index*#
%%%%
Type help(arima_index) for a list of topics in this file
%%%%
This file contains help on the following topics.  At present most topics
are not in final form.
acfarma()*      Macro to compute autocovariance function of ARMA model
arima()*        Macro to do unconditional least squares and MLE estimation
                of ARIMA model or linear regression with ARIMA errors
arimahelp()     Macro to provide help on macros in this file
arimares()      Macro to compute residuals from ARIMA model; used by
                macros, arima(), hannriss(), innovest().
detarma()*      Macro to compute determinant of Toeplitz covariance matrix
hannriss()*     Macro to do ARIMA fitting using Hannan-Rissannen algorithm
innovations()   Macro which uses cholesky() to compute the coefficients
                for one step prediction in terms of previous
                one-step prediction errors.  Used by innovest().
innovest()*     Macro to do ARIMA fitting using innovations algorithm
MASIGN          Explanation about the choice of signs for coefficients
                in ARIMA models and use of keywords 'masign' and
                'arsign'
moveoutroots()  Macro to fix up coefficients for a MA or AR operator so
                that all the zeros are outside the unit circle in the
                complex plane
neg2logLarma()* Macro to compute -2*log(L) and other quantities for
                ARIMA model, with possible seasonal structure and linear
                predictors
news            Items summarizing changes to macros in arima.mac
rhatcovar()     Macro to compute variances and or covariances of sample
                autocorrelations or the entire variance matrix of the
                sample autocorrelation function using Bartlett's formula
rhatvar()       Macro to compute factor for variances of sample cor-
                relations using Bartlett's formula
specarma()*     Macro to compute spectrum of ARMA model
_polish()       Macro used by hannriss() and innovest() to adjust their
                estimates nearer unconditional least squares estimates.
* = affected by MASIGN and ARSIGN and recognize keywords 'arsign' and
'masign'

====ARSIGN%#arima models
%%%%
Type arimahelp(MASIGN) for information on modifying sign conventions of
coefficients in ARIMA models.
%%%%
@@@@description#Description
Various macros, including arima(), hannriss(), innovest(),
neg2logLarma(), acfarma() and specarma(), allow for different sign
conventions in the definition of autoregression and moving average
coefficients.  This is controlled by keywords 'arsign' and 'masign'
and/or scalar variables ARSIGN and MASIGN, all with values -1 or +1.

@@@@see_also#Cross reference
See topic 'MASIGN' for details.
@@@@______

====detarma()#arima models
%%%%
detarma(phi, theta, n [,masign:Masign] [,arsign:Arsign] [,nfreq:Nfreq]),
  non-MISSING REAL vectors phi and theta, positive integer vector n,
  Masign and Arsign +1 or -1, positive integer scalar Nfreq
%%%%
@@@@introduction#Introduction
Help is not complete but is taken from the comments associated with
detarma().

detarma() computes the determinant of the Toeplitz covariance associated
with ARMA process.  At present it has no provision for seasonal
processes.

@@@@usage#Usage
 detarma(phi, theta, n [,masign:Masign] [,arsign:Arsign] [,nfreq:Nfreq])
  phi      non-MISSING REAL vector defining a causal AR operator
  theta    non-MISSING REAL vector defining a MA operator
  n        vector of positive integers
  Masign and/or Arsign must be +1 or -1; they determine sign conventions
           for interpreting phi and theta.  Type arimahelp(MASIGN)
           for information.
  Nfreq    Number of frequencies used in DFT to compute ACFV.  Default
           is the smallest integer >= 2*max(max(n), 200) that has no
           prime factors > 29

@@@@value_returned#Value returned
The result is a REAL vector d, the same length as n, with d[i] =
det(Sigma(n[i])), where Sigma(n[i]) = Cov[(X(1),...,X(n[i]))] when X(t)
is ARMA with AR coefficients phi and MA coefficients theta

detarma() computes the autocovariance function as the inverse Fourier
transform of the spectrum computed at nfreq frequencies.  From this it
computes that partial autocorrelations which are used, together with the
variance to compute the determinant.

@@@@______

====hannriss()#arima models,time domain,preliminary estimation
%%%%
hannriss(x, pdq:vector(p,d,q) [,degree:d1] [,maxlag:m] \
   [,polish:T, cycles:nc] [,arsign:Arsign] [,masign:Masign]), integers
   p >= 0, d >= 0 q >= 0, m >= p + q, nc >= 0, d1, Arsign and Masign +1
   or -1
%%%%
@@@@usage#Usage
hannriss(x,pdq:vector(p,d,q)), where p >= 0, d >= 0 and q >= 0 are
integers computes preliminary estimates of ARIMA coefficients using the
Hannan-Rissanen method.

The first step computes yulewalker estimates based on m autocorrelations
unless q = 0 when only p autocorrelations are used.  The default for m
is 20 + p + q.

Its value is structure(phi:phihat,theta:thetahat,xtxinv,rss:rss,nobs:n)

@@@@maxlag_and_degree_keywords#Keywords maxlag and degree
hannriss(x,pdq:vector(p,d,q),maxlag:m), integer m >= p + q, does the
same using m autocorrelations.

hannriss(x,pdq:vector(p,d,q),degree:d1 [,maxlag:m]) does the same,
except the possibly differenced data is detrended with a polynomial of
order d1.  d1 < 0 means nothing is subtracted.  The default for d1 is 0
when d = 0 and -1 when d > 0.

@@@@side_effect_variables#Side effect variables
In addition to returning a value, hannriss() creates the following side
effect variables
  COEF = vector(phihat,thetahat)
  XTXINV = analogue of solve(X' %c% X) matrix in regression
  ALLRESIDUALS = residuals from fitted model including backcast
   residuals
  RSS = sum(ALLRESIDUALS^2)
  NEG2LOGL = -2*log(likelihood)
  NPAR = p + q + degree1 + 1 = number of coefficients estimated
    where degree1 = max(d1,-1).

@@@@sign_conventions#Sign conventions
phi is defined so the autoregressive operator is Phi(B) = 1 +
Arsign*phi[1]*B + Arsign*phi[2]*B^2 + ... + Arsign*phi[p]*B^p.

theta is defined so that the moving average operator is Theta(B) = 1 +
Masign*theta[1]*B + Masign*theta[2]*B^2 + ... + Masign*theta[q]*B^q,
where B is the backshift operator.

The default values for Arsign and Masign are -1 and -1, but you may
change them by keyword phrases 'arsign:Arsign' and 'masign:Masign' or by
creating variables ARSIGN and/or MASIGN with values +1 or -1.  See topic
'MASIGNS' for details.

@@@@polish_cycles_keywords#Keywords polish and cycles
You can also include the following keyword phrases as arguments
  polish:T      carry out one or more extra "polishing" steps that
                should move the estimates closer to the unconditional
                least squares estimates.
  cycles:nc   nc > 0 polishing cycles will be carried out; default is 1

@@@@limitation#Limitation
There is currently no provision for seasonal ARIMAs
@@@@______

====innovations()#arima models,time domain,preliminary estimation
%%%%
innovations(gamma [,lag:m] [,final:T]), REAL vector gamma, integer m >
  0
%%%%
@@@@introduction#Introduction
innovations() computes the "innovation" algorithm given on p. 71 of
Brockwell, and Davis, computing a M by M matrix containing
coefficients and prediction variances.  It actually uses Cholesky
decomposition rather than the algorithm as given in Brockwell and Davis.

@@@@usage#Usage
innovations(gamma [,lag:M] [,final:T]), returns structure(theta:Theta,
v:V) where Theta is a M by M REAL matrix and V is a vector of length
M+1.  gamma is either a REAL n by n covariance matrix or a REAL vector
of length n containing an autocovariance function (ACVF).  M < n is a
positive integer with default value n - 1.

When gamma is an ACFV, innovations(gamma [,lag:M]) is the same as
innovations(toeplitz(gamma) [,lag:M])

The result is not affected by variables MASIGN if it exists.

@@@@definition#Definition
Suppose {X[1], X[2], ..., X[M+1]} is a sequence of random variables
with covariance matrix gamma (toeplitz(gamma) in ACVF case).

Then Theta is the M by M matrix with the property that

    Theta[j,j]*Z[j] + Theta[j,j-1]*Z[j-1] + ... + Theta[j,1]*Z[1]

is the best one step predictor of X[j+1] based on prediction
errors Z[1] = X[1], Z[2] = X[2] - Theta[1,1]*X[1], ..., Z[j] =
Theta[j-1,j-1]*Z[j-1]+Theta[j-1,j-2]*Z[j-2]+...+Theta[j-1,1]*Z[1].

Note that Z[j] could be expressed as a linear combination of X[j],
X[j-1], ..., X[1] so that the prediction is also the best prediction of
X[j+1] as a linear combination of X[j], X[j-1], ..., X[1].

V is a length M+1 vector with V[1] = gamma[1,1] = Var[X[1]] and V[j+1]
= prediction error variance when X[j+1] is predicted using X[1], X[2],
..., X[j] (or Z[1], ..., Z[j]).

@@@@final_keyword#Keyword final
When final:T is an argument, the result is
   structure(theta:vector(Theta[m,]), v:V[m+1])

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

====innovest()#arima models,time domain,preliminary estimation
%%%%
innovest(x,pdq:vector(p,d,q) [,maxlag:m] [,degree:degree]\
  [,polish:T,cycles:nc] [,checkroots:F] [,silent:T] [,arsign:Arsign]\
  [,masign:Masign]), integers p>=0, d >= 0, q >= 0, maxlag >= p+q,
  degree, nc >= 1, Arsign and Masign = +1 or -1.
Value is structure(phi:phihat, theta:thetahat, xtxinv, nobs:n, npar,
  rss:residSS, neg2logL:-2*log(likelihood), aicc:AICC)
%%%%
@@@@introduction#Introduction
innovest() is a macro to compute the innovations preliminary estimate of
coefficients for an ARIMA(p,d,q) time series as described on pp 151-153
of Brockwell and Davis.  Keywords 'arsign' and 'masign' allow you to
specify the sign convention to be used in defining parameters.  See
below.

@@@@usage#Usage
innovest(x, pdq:vector(p,d,q) [,maxlag:M]), where x is a REAL vector
and p, d and q are nonnegative integers return a structure summarizing
the results of the preliminary innovations estimates for an
ARIMA(p,d,q) model fit to x.  M >= p + q is an integer with default
value p + q + max(15, p + q).  See below for the form of the results.

There is currently no provision for estimating seasonal ARIMA models

@@@@side_effect_variables#Side effect variables
innovest() creates the following side effect variables
  COEF = vector(phihat,thetahat), the estimated AR and MA coefficients
  ALLRESIDUALS = residuals from fitted model including backcast
    residuals
  RSS = sum(ALLRESIDUALS^2)
  NEG2LOGL = -2*log(likelihood) using the estimates found
  NPAR = p + q + degree + 1 = number of coefficients estimated

@@@@keywords#Optional keyword phrase arguments
There are several optional keyword phrases which affect what innovest()
does:
  degree:d1      A polynomial trend of order d1 is removed (after
                 differencing when d > 0).  Nothing is removed, not
                  even a mean, when d1 < 0.  The default for d1 is
                  -d (mean removed when d = 0, nothing otherwise).
  polish:T       The estimates are adjusted by one or more cylcels
                 of an approximate iteration in the direction of the
                  least squares estimates.
  cycles:nc      nc > 0 an integer, default 1 is the number of
                 "polishing" cycles.
  checkroots:F   suppresses checking for stationarity and invertability
  silent:T       suppress warning messages
  arsign:Arsign  +1 or -1; alters definition of AR paramaters; see below
  masign:Masign  +1 or -1; alters definition of MA paramaters; see below

@@@@value_returned#Value returned
innovest() returns as value
  structure(phi:phihat, theta:thetahat, nobs:n, xtxinv:xtxinv,
    npar:p+q+d1+1, rss:residSS, neg2logL:-2*log(likelhihood),aicc:aicc)

phihat and thetahat are the estimated AR and MA coefficients (NULL when
p or q is 0).

xtxinv is NULL without polish:T or when p = q = 0, Otherwise xtxinv is
an analogue of the regression solve(X' %*% X) derived from the final
polishing step.  Its diagonal elements can be used to compute
approximate standard errors of the coefficients.

residSS = sum of squares of all residuals, including those backcast.

The likelihood is the normal likelihood computed using backcasting and
includes a factor of (2*PI)^(-n/2) and aicc is a modification of
Akaike's information criteria (AIC).

Note that innovest() does not return a mean or other estimates of
detrending parameters.

@@@@check_on_operators#Check on operators
By default, estimated AR and MA coefficients are checked to see if they
define stationary (causal) and invertable operators, respectively.  If
they do not, a warning message printed and any "polishing" cycles (see
below) are suppressed and components rss, neg2logL and aicc of the
result are set to MISSING.

@@@@sign_conventions#Sign conventions
See topic 'MASIGN' for information on how Arsign (default -1 or the
value of ARSIGN if it exists) and Masign (default -1 or the value of
MASIGN if it exists) modify the definitions of the AR and MA parameters.

To use the convention used by Brockwell and Davis before you start the
analysis, you should use
  Cmd> MASIGN <- 1; ARSIGN <- -1

The convention used by Box and Jenkins (the default) corresponds to
Arsign = -1, Masign = -1.

@@@@see_also#Cross references
See also arima(), hannriss(), innovations().
@@@@______

====MASIGN%#arima models
%%%%
Type arimahelp(MASIGN) for information on modifying sign conventions of
coefficients in ARIMA models.
%%%%
@@@@sign_conventions#Sign conventions
Various macros, including arima(), hannriss(), innovest(),
neg2logLarma(), acfarma() and specarma(), allow for different sign
conventions in the definition of autoregression and moving average
coefficients.  The convention to be used may be determined by keyword
phrases 'arsign:Arsign' and 'masign:Masign' or the values of variables
ARSIGN or MASIGN, when they exist.

Some functions and macros, including movavg(), autoreg(), polyroot(),
and moveoutroots() *always* use a sign convention equivalent to Arsign
= Masign = -1 and are *not* affected by keywords 'arsign' and 'masign'
or the values of ARSIGN and MASIGN.

For the the affected macros, an ARMA model is assumed to have the form
 X[t] + Arsign*(phi[1]*X[t-1]+phi[2]*X[t-2]+...+phi[p]*X[t]) =
     Z[t] + Masign*(theta[1]*Z[t-1]+theta[2]*Z[t-2]+...+theta[q]*Z[t-q])
where {Z[t]} is white noise.

Arsign is -1 or +1.  Without 'arsign:Arsign', the default value of
Arsign is the value of variable ARSIGN when it exists or -1 when it
does not.

Masign is -1 or +1.  Without 'masign:Masign', the default value of
Masign is the value of variable MASIGN when it exists or -1 when it
does not.

@@@@examples_of_sign_conventions#Examples of sign conventions
Examples for an ARMA(2,2) model with zero mean and AR coefficients phi
and MA coefficients theta.

  For the defaults arsign:-1 and masign:-1, the model is
    Y[t] = phi[1]*Y[t-1] + phi[2]*Y[t-1] +
                Z[t] - theta[1]*Z[t-1] - theta[2]*Z[t-2]
  With arsign:-1 and masign:1, the model is
    Y[t] = phi[1]*Y[t-1] + phi[2]*Y[t-1] +
                Z[t] + theta[1]*Z[t-1] + theta[2]*Z[t-2]
  With arsign:1 and masign:-1, the model is
    Y[t] = -phi[1]*Y[t-1] - phi[2]*Y[t-1] +
                Z[t] - theta[1]*Z[t-1] - theta[2]*Z[t-2]
  With arsign:1 and masign:1, the model is
    Y[t] = -phi[1]*Y[t-1] - phi[2]*Y[t-1] +
                Z[t] + theta[1]*Z[t-1] + theta[2]*Z[t-2]

@@@@defaults#Defaults
The convention used by Box and Jenkins and many others corresponds to
Arsign = -1 and Masign = -1.  This is the convention that is implicit
in the definition of MacAnova functions autoreg(), movavg() and
polyroot().  These functions do not recognize 'arsign' and 'masign' and
are not affected by the values of ARSIGN and MASIGN.

The convention used by Brockwell and Davis corresponds to Arsign = -1
and Masign = +1.  If you want results consistent with this convention
you should do the following before you start your work:
  Cmd> MASIGN <- 1; ARSIGN <- -1

This is easier than using 'masign:1' whenever you are using one of the
affected macros.

@@@@unaffected_macros#Macros not affected by masign and arsign
Note that macros innovations() and arimares() are not affected by this
notation, but macros that use these compensate.

@@@@see_also#Cross references
See also autoreg(), movavg(), polyroot(), acfarma(), arima(),
arimares(), hannriss(), innovest(), neg2logLarma(), specarma().
@@@@______

====moveoutroots()#time domain,arima models
%%%%
moveoutroots(theta), REAL vector theta
%%%%
@@@@usage#Usage
moveoutroots(theta) where theta is a REAL vector defining the
polynomial P(z) = 1 - sum(theta*z^run(n)), n = length(theta), returns a
new vector theta1 defining a new polynomial P1(z) of the same form which
has the same zeros as P(z) except that any zero z0 with |z0| < 1 is
replaced by 1/z0.

Thus all the zeros of P1(z) are outside the unit circle in the complex
plane and moveoutroots(theta) are the coefficients of an invertible MA
operator and, when no zero has modulus 1, the coefficients of a
stationary (causal) AR operator.

Equivalently, if Q(z) = z^n - sum(theta*z^run(n-1,0)) and Q1(z) =
z^n - sum(theta1*z^run(n-1,0), any zero z0 of Q(z) with |z0| > 1 is
replaced by 1/z0.

@@@@sign_convention#Sign convention
moveoutroots assumes the sign convention used by Box and Jenkins and
implicit in movavg(), autoreg() and polyroot() (Arsign = -1, Masign =
-1).  When theta are MA coefficients using the sign convention of
Brockwell and Davis (Masign = +1), you should use moveoutroots(-theta).
See topic 'MASIGN'.

@@@@examples#Examples
Examples:
  Cmd> moveoutroots(2)
  (1)         0.5

This is correct because 1 - 2*z has zero .5 < 1 and 1 - .5*z has zero 2
> 1.

  Cmd> moveoutroots(vector(2,-1.25))
  (1)         1.6        -0.8

This is correct because 1 - 2*z + 1.25*z^2 has zeros z = .8 +- .4*i with
|z| = sqrt(.8) = 0.89443 < 1 and 1 - 1.6*z + .8*z^2) has roots z =
(1+-.5*i) = 1/(.8 -+ .4*i) with |z| = sqrt(1.25) = 1.118 > 1.

@@@@see_also#Cross references
See polyroot(), movavg(), autoreg().
@@@@______

====neg2logLarma()#arima models,time domain
%%%%
neg2logLarma(y,coefficients, [,x:x] [,pdq:vector(p,d,q)] \
  [,PDQ:vector(P,D,Q),seasonal:s] [,fitmean:T or F] \
  [,cast:m] [,cycles:ncyc] [,sigmasq:sigmasq] [,neg2logL:F],
  [,residuals:T] [,logdet:T] [,sigmahatsq:T] [,all:T]
  [,masign:1 or -1] [,arsign:1 or -1])
  nonMISSING REAL vector y, optional nonMISSING REAL matrix x,
  integers p, d, q, P, D, Q, m, ncyc, >= 0, integer s > 0,
  sigmasq REAL scalar > 0 or MISSING
%%%%
@@@@introduction#Introduction
Macro to compute -2*log(L) and other quantities for an ARIMA model with
specified parameters.  L is the likelihood or concentrated likelihood.

The calling sequence is very similar to that of arima(), except you must
provide values of the coefficients and there are keywords which specify
what is returned.

@@@@usage_arguments#Usage and arguments
 neg2logLarma(y,b, [,x:x], [pdq:vector(p,d,q)] \
  [,PDQ:vector(P,D,Q),seasonal:s] [,fitmean:T or F] \
  [, sigmasq:sigmasq] [,cast:m] [cycles:ncyc] [,neg2logL:F],
  [,masign:Masign] [,arsign:Arsign]\
  [,residuals:T] [,logdet:T] [,sigmahatsq:T] [,all:T])

y is a real nonmissing vector of length n (the response)
b = vector(mu, beta, phi, theta, phiS, thetaS) of coefficients (mu =
   mean of mean difference, beta = slopes of linear predictors, phi =
   non-seasonal AR coefficients, theta = non-seasonal MA coefficients,
   phiS = seasonal AR coefficients, thetaS = non-seasonal MA
   coefficients).  Any parts that are not needed are NULL or omitted.
x = optional non-MISSING REAL matrix of linear predictors, with
   nrows(x) = nrows(y)
p, d, q >= 0 are integers specifying the non-seasonal ARIMA component
P, D, Q >= 0 are integers specifying the seasonal ARIMA component
s > 1 an integer specifying the season length; required with any of P,
  D or Q non-zero
fitmean:T means mu is in the model; fitmean:F means mu is not in
  model.  The defualt is T when d = D = 0 and F otherwise
sigmasq a REAL scalar, either MISSING or positive.  If MISSING (the
  default) the innovation variance is estimated and the "concentrated"
  log likelihood is computed.  Otherwise, the log(L) is computed
  assuming innovation variance sigmasq
m >= 0 and ncyc >= 0 are integers controlling backcasting in computing
  residuals.  See arima() and arimares() for more information.
Arsign and Masign must be +1 or -1 and control what sign convention is
  used for coefficients.  See topics 'MASIGN' and 'ARSIGN' for details.

@@@@value_returned#Value returned
What is returned is controlled by keywords 'neg2logL', 'residuals',
'logdet', 'sigmahatsq' and 'all'.

  Keyword                Value returned
 neg2logL:T    -2*log(L).  When sigmasq is MISSING, L is the
                           concentrated likelihood
 residuals:T   vector of residuals, including backcast ones
 logdet:T      log(det(Gamma)) where sigmasq*Gamma is covariance matrix
 sigmahatsq:T  The estimated innovation variance when sigmasq is MISSING
               or sigmasq, otherwise
 all:T         Same neg2logL:T,residuals:T,logdet:T,sigmahatsq:T

If more than items is to be returned, the output is a structure with
the keywords as component names; if only one is T, the output is a
scalar or vector.  With all:T you can suppress a result by, say,
'residuals:F',

@@@@examples#Examples
To get the concentrated likelihood for an ordinary ARMA(p,q) model, use
  Cmd> result <- neg2logLarma(y,vector(mu,phi,theta), pdq:vector(p,0,q))

For an ARIMA model with d > 0 and zero mean differences, use
  Cmd> result <- neg2logLarma(y,vector(phi,theta), pdq:vector(p,d,q))

For a seasonal (p,q)x(P,Q) seasonal ARMA with no differences and, say
quarterly data, use
  Cmd> result <- neg2logLarma(y,vector(mu,phi,theta),\
    pdq:vector(p,0,q),PDQ:vector(P,0,Q),season:4)

For a seasonal (p,d,q)x(P,D,Q) seasonal ARIMA with no differences and,
say quarterly data, use
  Cmd> result <- neg2logLarma(y,vector(phi,theta),\
    pdq:vector(p,d,q),PDQ:vector(P,D,Q),season:4)

For an ARMA(p,q) model with mean that is cubic in time, try
  Cmd> predictors <- (run(n) - (n+1)/2)^run(0,3)'

  Cmd> result <- neg2logLarma(y,vector(beta,phi,theta),x:predictors,\
      pdq:vector(p,0,q), fit meanF) # beta a vector of length 4

@@@@see_also#Cross reference
Use arimahelp() to also see topic arima().
@@@@______

====news*#
@@@@july_2001#July 2001
010724 levmar() now returns a hessian computed from the jacobian after
the last iteration.  Previously it was computed from the jacobian at the
start of the last iteration.

arima() creates an additional side effect variable AICC.

Fixed arima() bugs affecting mle:T case.
  Now length(HII) = length(RESIDUALS) and nrows(JACOBIAN) =
  length(ALLRESIDUALS) (previously HII and JACOBIAN had one extra row).
  XTXINV is computed from the corrected JACOBIAN and standard errors of
  coefficients are computed from the corrected XTXINV.

@@@@january_2001#January 2001
010103 Many small changes in macros, especially in determining 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

  Variable S is no longer used as a default number of frequencies in
  macros such as neg2logLarma() that do not return a frequency function.

  When variable S is to be used as a default number of frequencies, it
  is an error if S has a prime factor > 29.

  The minimum default value for 'cast' in arima() and neg2logLarma() has
  been increased to 50 for models with an AR component.

@@@@december_2000#December 2000
001227 modified acfarma() so can be used as acfarma(phi,theta,L ...)
which is equivalent to acfarma(phi,theta,lag:L ...)

001207 Fixed arima() bug in computing log(det) with mle:F and either
Masign = 1 or Arsign = 1

@@@@november_2000#November 2000
001130 MLE iteration in arima() modified

Computation of logdet changed and fixed bugs in arimares().

001117  Added neg2logLarma() and fixed arimares() bug.

@@@@______

====rhatcovar()#time domain,autocorrelation
%%%%
rhatcovar(rho,i, j) or rhatcovar(rho,lag:L), REAL vector rho, integers
  i > 0, j > 0, L > 0
%%%%
@@@@usage#Usage
rhatcovar(rho [,lag:L]), where rho is a REAL vector of auto-
correlations, computes n*Var[rhohat], where Var[rhohat] is the large
sample covariance matrix of the sample autocorrelation function
rhohat, computed using Bartletts' formula.

The result is a REAL L by L matrix.  The default value of L is
(maximum lag for rho)/2.

When rho[1] == 1, rho[k] is assumed to contain the lag k-1 auto-
correlation; otherwise, rho[k] is assumed to contain the lag k
autocorrelation.  Thus if gamma is a REAL vector containing the
autocovariance function with gamma[1] = Var[X(t)], both
rhatvar(gamma/gamma[1]) and rhatvar(gamma[-1]/gamma[1]) return the same
result.

rhatcovar(rho,i,j) returns n*covar{rhohat(i),rhohat(j)}, that is,
rhatcovar(rho)[i,j].

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

====rhatvar()#time domain,autocorrelation
%%%%
rhatvar(rho [,lag:L]), REAL vector of correlations rho, integer L > 0
%%%%
@@@@usage#Usage
rhatvar(rho [,lag:L]), where rho is a REAL vector of autocorrelations
computes the vector of n*(var(rhohat(1)),var(rhohat(2)),...,
var(rhohat(lag)) using Bartlett's formula, where rhohat are sample
autocorrelations from a series with ACF rho.

When rho[1] == 1, rho[k] is assumed to contain the lag k-1 auto-
correlation; otherwise, rho[k] is assumed to contain the lag k
autocorrelation.  Thus if gamma is a REAL vector containing the
autocovariance function with gamma[1] = Var[X(t)], both
rhatvar(gamma/gamma[1]) and rhatvar(gamma[-1]/gamma[1]) return the same
result.

The default value of L is (maximum lag for rho)/2

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

====specarma()#frequency domain,time series,spectrum analysis
%%%%
specarma(phi ,theta [, nfreq:nf]), REAL vectors phi and theta, integer
  nf > 0
%%%%
@@@@usage#Usage
specarma(phi, theta [, nfreq:Nfreq]) computes the spectrum of an ARMA
time series with AR coefficients in phi and MA coefficients in theta and
innovation variance 1.

phi and theta are REAL vectors with no MISSING values. To omit part of
the model, use phi = 0 or theta = 0

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

The result is a vector of length nf containing the spectrum at
frequencies 0, 1/Nfreq, 2/Nfreq, ..., (Nfreq-1)/Nfreq cycles per unit
time.

Without nfreq:Nfreq, when no positive integer scalar S exists, the
default is Nfreq = 400.  When S does exist and is a positive integer
scalar, the default for Nfreq is S.  It is an error if S has a prime
factor > 29.

@@@@sign_conventions#Sign conventions
You can use keyword phrases arsign:Arsign and masign:Masign, where
Arsign and Masign are +1 or -1, to modify the interpretation of the
coefficients in the ARMA model.  The default for Arsign is variable
ARSIGN if it exists or -1 if not.  The default for Masign is variable
MASIGN if it exists or -1 if not.

The ARMA model assumed is defined as
  (1+Arsign*sum(phi*B^run(p)))X[t] = (1+Masign*sum(theta*B^run(q)))Z[t].
where B is the backshift operator and {Z[t]} is white noise with mean 0
and variance 1.

Arsign = -1 and Masign = -1 correspond to the convention used by Box
and Jenkins, and Arsign = -1 and Masign = +1 correspond to the
convention used by Brockwell and Davis

@@@@example#Example
You can plot the spectrum for the ARMA model fit by arima() by
  Cmd> result <- arima(y,pdq:vector(p,0,q), keep:T)

  Cmd> ffplot(specarma(result$phi, result$theta))

@@@@see_also#Cross references
See also arima(), acfarma(), ffplot().
@@@@______

====_polish()*#
%%%%
_polish(x,structure(phi:phi,theta:theta),p,q,cycles)
%%%%
@@@@usage#Usage
_polish(x,structure(phi:phi,theta:theta),p,q,cycles) carries out cycles
repetions of an approximate least squares algorithm to adjust the
values of phi and theta.  It is used by hannriss() and innovest() to
improve the first values computed.
@@@@______

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

