\section*{Examples of ARIMA fitting}
The purpose of these notes is to take a set of data (river data from
the trinity river) and attempt a first pass at ARIMA modelling.  The
tools we will use are arima.mle (to do the mle estimation), arima.diag
(to do the diagnostic checking), acf (to get ACF and PACF plots), diff
(to do differences of the series), and tsplot(to plot things).
In order to help a little, I'll include the help pages for three of
these commands here
{\small
  \begin{verbatim}

Estimate    Autocovariance,    Autocorrelation     or     Partial
Autocorrelation


DESCRIPTION:
       Estimates and displays autocovariance, autocorrelation  or
       partial autocorrelation functions.

USAGE:
       acf(x, lag.max=<<see below>>, type="correlation", plot=T)

REQUIRED ARGUMENTS:
x:     a univariate or multivariate time series or  a  vector  or
       matrix.  Missing  values are allowed only at the beginning
       or end of  series.  If  a  matrix,  rows  are  treated  as
       timepoints and columns as univariate series.

OPTIONAL ARGUMENTS:
lag.max:    the maximum number of lags at which to  estimate  the
       autocovariance.   If  this is not supplied, it is a number
       proportional to the logarithm of the length of the series.
type:      a  character  string:  "covariance"  to  estimate  the
       autocovariance    function,    "correlation"    for    the
       autocorrelation function, or  "partial",  if  the  partial
       autocorrelation  function is desired.  The start of one of
       the strings will suffice.
plot:     If TRUE, the autocovariance or autocorrelation function
       between  pairs  of univariate series will be plotted in an
       array of at most 5 by  5  plots  per  page.   If  type  is
       "correlation"  or  "partial"  approximate  95%  confidence
       limits are drawn on the plots.

VALUE:
       a list with the following components:
acf:   a three-dimensional array containing the autocovariance or
       autocorrelation  function  estimates.   acf[i,j,k]  is the
       covariance (or correlation) between  the  j-th  series  at
       time t and the k-th series at time t+1-i.
lag:   an array the same shape as acf  containing  the  lags  (as
       fractions   of  the  sampling  period)  at  which  acf  is
       calculated.  If j > k  and  i  >  1,  then  lag[i,j,k]  is
       negative.
n.used:   the number of observations in which no  missing  values
       occur.
type:     a character string indicating  the  type  of  function,
       "covariance", "correlation" or "partial".
series:   the name of x, including transformations.

SIDE EFFECTS:
       if plot is TRUE, a high density  plot  (type="h")  of  the
       resulting  function  will  be  produced  including the 95%
       confidence limits about zero.

DETAILS:
       The autocovariance function is estimated  by  summing  the
       lagged  products and dividing by the length of the series.
       For the  autocorrelation  function,  all  covariances  are
       further divided by the geometric mean of the corresponding
       variances.  For the partial autocorrelation function,  the
       Levinson-Durbin recursion is used to fit AR(p) models to x
       successively for p = 1, ..., lag.max,  and  from  the  AR-
       coefficients  the  partial autocorrelation function values
       are derived.

BACKGROUND:
       Both  the  autocorrelation  function  (acf)  and   partial
       autocorrelation  function (pacf) are useful tools in ARIMA
       model   identification.    The   autocovariance   function
       describes  the  serial  dependence  structure  of  a  time
       series.  The  time  series  is  generally  assumed  to  be
       stationary before these tools are used for Box and Jenkins
       modeling purposes.

REFERENCES:
       The Time Series chapter of the S-PLUS User's Manual.

SEE ALSO:
       acf.plot, ar, arima.mle.

EXAMPLES:
       acf.lynx <- acf(lynx, 36, "correlation")
       pacf.lynx <- acf(lynx, 36, "partial")

  \end{verbatim}
}

{\small
  \begin{verbatim}

ARIMA Modeling via Gaussian Maximum Likelihood


DESCRIPTION:
       Returns a  list  representing  a  univariate  ARIMA  model
       estimated by Gaussian maximum likelihood.

USAGE:
       arima.mle(x, model, n.cond=<<see below>>, xreg=NULL,  ...)

REQUIRED ARGUMENTS:
x:     a univariate time series  or  a  vector.   Missing  values
       (NAs) are allowed.
model:    a list specifying an ARIMA model.  If  a  simple  ARIMA
       model  is specified (i.e., only one component), then model
       should be a list with component names from: order, period,
       ar,  ma,  ndiff,  ar.opt,  ma.trans,  and ma.opt.  See the
       DETAILS section for further details.

OPTIONAL ARGUMENTS:
n.cond:   the number of observations on which  to  condition  the
       likelihood.   n.cond  must  be at least P+D (the default),
       where  P  and  D  are   the   orders   of   the   expanded
       autoregressive and differencing polynomials.
xreg:     a univariate or multivariate time series or a vector or
       matrix with univariate time series per column.  These will
       be used as additive regression variables.
...:   additional arguments can be passed to nlmin.

VALUE:
       a list representing the result of the fitted model:
model:     the  same  as  the  input  model  with  the  estimated
       coefficients substituted in the components ar and ma.
var.coef:   the variance-covariance matrix for the autoregressive
       and moving average coefficients.
loglik:   -2 times the log likelihood  of  the  model  (up  to  a
       constant factor).
aic:   Akaike's information criteria,  which  is  loglik  plus  2
       times the number of parameters fit.
sigma2:   the estimated innovations variance.
n.used:     the  number  of  observations  used  to  compute  the
       likelihood.
n.cond:   the number of observations on which the  likelihood  is
       conditioned.
reg.coef:   the estimated regression  coefficients  (returned  if
       xreg is provided).
converged:      if  the  optimizer  has  apparently  successfully
       converged  to a minimum, then converged is TRUE; otherwise
       converged is FALSE.
conv.type:       a  character  string  describing  the  type   of
       convergence.
series:   a character string giving  the  name  of  x,  including
       transformations.
reg.series:    a  character  string  giving  the  name  of  xreg,
       including   transformations  (only  returned  if  xreg  is
       provided).

DETAILS:
       SPECIFYING A MODEL.  The components of the model  argument
       should  be  some or all of:  order, period, ar, ma, ndiff,
       ar.opt, ma.trans, and ma.opt.  order is a vector of length
       3 specifying the order of an ARIMA(p,d,q); the elements of
       the vector are (p, d, q), where p and q are the  order  of
       the  autoregressive  and moving average operators and d is
       the number of differences.  period specifies the period of
       the ARIMA operators; e.g., a seasonal operator for monthly
       data would have a period of 12 (the default is 1).  ar and
       ma  are  vectors  of initial values for autoregressive and
       moving average coefficients respectively (if not provided,
       the  initial values to the optimizer are set to 0).  ndiff
       is  an  alternative  way  to   specify   the   number   of
       differences.   A model can be specified using either order
       or ar, ma, and ndiff, or both (but they must agree if both
       are  specified).  The model doesn't explicitly allow for a
       non-zero mean of the  series,  though  the  estimation  is
       invariant to the mean of the data for many of the possible
       models; if you estimate an AR model, you must subtract the
       mean  yourself  or  use  the xreg argument in order to get
       sensible estimates.

       If ma.trans is TRUE  (the  default),  the  moving  average
       coefficients  will  be  transformed before passing them to
       the  optimizer,  ensuring  invertibility  of  the   model.
       ar.opt  and  ma.opt  are logical vectors of length p and q
       respectively.  If  the  ith  element  is  TRUE,  then  the
       optimizer  will  optimize  over  the ith autoregressive or
       moving average coefficient.  By default, the  vectors  are
       TRUE.   This  option is useful to fit models in which some
       low order coefficients are set to 0.

       If a multiplicative ARIMA model is specified  (i.e.,  more
       than  one component), then model is a list of lists.  Each
       list represents a component of  the  multiplicative  ARIMA
       model, and is specified in the same manner as above (i.e.,
       the components of each list  are  ar,  ma,  ndiff,  order,
       period, ma.trans, ar.opt, ma.opt).

       THE ALGORITHM.  The likelihood is  conditioned  on  n.cond
       observations.   If  no  missing  values  are  present, the
       likelihood   is   computed   by   forming   the   Choleski
       decomposition of the covariance matrix of the process (see
       Ansley, 1979).  If missing values are  present,  then  the
       likelihood  is computed using the Kalman filter applied to
       a state space representation of the likelihood.  The state
       space  representation is the same as that used by Kohn and
       Ansley (1986).  The  method  of  initializing  the  Kalman
       filter  recursions  is  that  given  by  Bell  and Hillmer
       (1987).  The regression parameters are concentrated out of
       the  likelihood,  as  in  Kohn  and  Ansley  (1985).   The
       autoregressive   and   moving   average   parameters   are
       transformed to ensure stationarity and invertibility as in
       Jones (1980).

WARNING:
       Unlike the ar function, the mean of the series will not be
       estimated  by  arima.mle unless you use the xreg argument.
       arima.mle assumes a zero mean series.

REFERENCES:
       Ansley,  C.  F.  (1979).   An  algorithm  for  the   exact
       likelihood   of   a  mixed  autoregressive-moving  average
       process.  Biometrika 66, 59-65.

       Bell, W. and Hillmer, S. (1987).  Initializing the  Kalman
       filter   in   the  nonstationary  case.   Research  Report
       CENSUS/SRC/RR-87/33, Statistical Research Division, Bureau
       of the Census, Washington, DC, 20233.

       Jones, R. H. (1980).  Maximum likelihood fitting  of  ARMA
       models   to   time   series   with  missing  observations.
       Technometrics 22, 389-395.

       Kohn, R. and Ansley, C. F. (1985).   Efficient  estimation
       and   prediction   in   time   series  regression  models.
       Biometrika 72, 694-697.

       Kohn,  R.  and  Ansley,   C.   F.   (1986).    Estimation,
       prediction,   and  interpolation  for  ARIMA  models  with
       missing  data.   Journal  of  the   American   Statistical
       Association 81, 751-761.

       The Time Series chapter of the S-PLUS User's Manual.

SEE ALSO:
       ar, arima.diag, arima.filt, arima.forecast, nlmin, sabl.

EXAMPLES:
       # Simulate an MA(2) series and fit the series using Gaussian
       # Maximum Likelihood
       ma <- arima.sim(model=list(ma=c(-.5, -.25)))
       ma.fit <- arima.mle(ma, model=list(ma=c(-.5, -.25)))

       # Fit a Box-Jenkins (0,1,1)x(0,1,1)12 Airline model to the ship data
       # Use zeros as the starting values for the optimizer
       al.mod <- list(list(order=c(0,1,1)), list(order=c(0,1,1), period=12))
       fit <- arima.mle(ship, model=al.mod)

  \end{verbatim}
}
{\small
  \begin{verbatim}

Compute Diagnostics for ARIMA Model


DESCRIPTION:
       Computes diagnostics for an ARIMA model.  The  diagnostics
       include the autocorrelation function of the residuals, the
       standardized residuals, and the  portmanteau  goodness  of
       fit test statistic.

USAGE:
       arima.diag(z, acf.resid=T, gof.lag=10, lag.max=<<see below>>, resid=F,
                  std.resid=T, plot=T, ...)

REQUIRED ARGUMENTS:
z:     a list like the output from arima.mle.

OPTIONAL ARGUMENTS:
acf.resid:     logical flag: if TRUE, the autocorrelation of  the
       residuals will be returned.
lag.max:    the maximum number of lags at which to  estimate  the
       autocovariance.   If  this  is  not  supplied,  it  is the
       maximum  between  gof.lag  plus  the   number   of   model
       parameters  and  a number proportional to the logarithm of
       the length of the series.
gof.lag:    if gof.lag>0, then gof.lag plus the number  of  model
       parameters  is the number of lags to use for computing the
       Portmanteau goodness of fit statistic.  If gof.lag=0, then
       the statistic will not be computed
resid:    logical flag: if  TRUE,  then  the  residuals  will  be
       returned.
std.resid:     logical  flag:  if  TRUE,  then  the  standardized
       residuals will be returned.
plot:     logical flag: if TRUE, the diagnostics will be  plotted
       using the function arima.diag.plot
...:    additional  arguments  may  be  passed  to  the  function
       arima.diag.plot

VALUE:
       a list (which is returned invisibly when  plot=TRUE)  with
       the following elements:
acf.list:   a list representing the autocorrelation  function  of
       the residuals.  See acf for details.
gof:   a  list  representing  the  Portmanteau  goodness  of  fit
       statistics  computed  for  a  range of lags.  The list has
       four elements: lag, statistic,  df,  p.value.   lag  is  a
       vector   of  the  number  of  lags  used  to  compute  the
       statistics.   statistic  is  the  vector   of   statistics
       corresponding  to  each  lag  used.   df  is the number of
       degrees of freedom the test statistics have under the null
       hypothesis that the model is correct.  p.value is a vector
       of the p-values for the  statistics  using  a  Chi-Squared
       distribution with the appropriate degrees of freedom.
resid:    the residuals or innovations for the process.
std.resid:     the standardized  residuals.   The  residuals  are
       standardized  to  have  unit variance under the assumption
       that the model is correct and the process is Gaussian.
series:   the name of x, including transformations.

SIDE EFFECTS:
       if plot is TRUE, the diagnostics will be plotted using the
       function arima.diag.plot

DETAILS:
       The residuals (both standardized  and  raw)  are  computed
       using   the   function  arima.filt.   The  autocorrelation
       function of the residuals is computed using  the  function
       acf.   The  portmanteau test statistic is derived from the
       autocorrelation function of the residuals  (see  the  Time
       Series chapter in the S-PLUS User's Manual for details).

REFERENCES:
       Box, G. E. P. and Jenkins,  G.  M.  (1976).   Time  Series
       Analysis:  Forecasting  and Control.  Holden-Day, Oakland,
       Calif.  Chapter 8.

       The Time Series chapter of the S-PLUS User's Manual.

SEE ALSO:
       arima.diag.plot arima.filt, arima.mle.

EXAMPLES:
       # compute and plot diagnostics for simulated AR(1) series with mean 5
       x <- arima.sim(model=list(ar=.9)) + 5
       xreg <- rep(1,100)
       fit <- arima.mle(x,model=list(ar=.9), xreg=xreg)
       diag <- arima.diag(fit)

       lynx.arma11 <- arima.mle(lynx, model=list(ar=0, ma=0))
       arima.diag(lynx.arma11)

  \end{verbatim}
}

The data we with to look at have the following first few lines
\begin{verbatim}
TRINITY RIVER AT LEWISTON, CALIFORNIA, OCTOBER 1912 - SEPTEMBER 1960
83.42 100.71 86.5 106.15 107.77 78.79 
85.67 110.47 107.49 76.07 90.48 105.36 
74.58 99.29 89.93 113.55 101.06 59.08 
97.34 105.89 108.67 83.38 82.43 104.72 
\end{verbatim}

We read in the data as (with skip to skip the first line of rubbish).
\begin{verbatim}
trinity <- as.ts( scan("trinity.dat", skip=1))
\end{verbatim}

The first thing we want to do is a number of plots of the data, and
acf, pacf plots.
Notice that the seasonality is easier to see in the partial data
(second plot) and that the 12 month effects show up strongly in the
ACF and PACF
\begin{verbatim}
> tsplot(trinity)       
> title("tsplot(trinity)")
> tsplot(trinity[1:96])   
> title("tsplot(trinity[1:96])")
> acf(trinity)                  
> acf(trinity, type="p")   
\end{verbatim}
\centerline{\psfig{figure=trin.eps,angle=-90,width=6in}}

We have a lot of reason to think about 12th order difference of the
data, so lets try that.
\begin{verbatim}
> delme <- diff(trinity, lag=12)
> tsplot(delme)
> title("tsplot(delme)")
> tsplot(delme[1:96])   
> title("tsplot(delme[1:96])")
> acf(delme)
> acf(delme, type="p")
\end{verbatim}
\centerline{\psfig{figure=trin2.eps,angle=-90,width=6in}}

Those plots now look quite good.  One might temporarily guess that a
MA(2) model might do the trick, or maybe an ARMA(1,1).   Lets try
those and see where we get.
First fit the model that has 12th order differences, and is otherwise
an MA(2) model.
\begin{verbatim}
>  mymod <- arima.mle(trinity, model=list( list(ndiff=1, period=12), 
+ list(order=c(0,0,2))))                                            
> mymod$loglik
[1] 4271.951

\end{verbatim}
Also try the same except an ARIMA(1,0,1) model.
\begin{verbatim}
>  mymod2 <- arima.mle(trinity, model=list( list(ndiff=1, period=12),
+ list(order=c(1,0,1))))  
 > mymod2$loglik                                                      
[1] 4256.779
\end{verbatim}
Just for kicks we'll try the raw AR(1) and raw MA(1) and see how they
do.
\begin{verbatim}
> arima.mle(trinity, model=list( list(ndiff=1, period=12),           
+ list(order=c(1,0,0))))$loglik                           
[1] 4307.937
> arima.mle(trinity, model=list( list(ndiff=1, period=12),
+ list(order=c(0,0,1))))$loglik                           
[1] 4305.263
\end{verbatim}
Either of these is demonstrably works (likelihood ratio test) than the
combined models.
Lets quickly look at the diagnostics for the two candidate models.
\begin{verbatim}
> arima.diag(mymod)
> arima.diag(mymod2)
\end{verbatim}
I could not not figure out how to automatically get the pacfs, so here is
the manual approach to that
\begin{verbatim}
> resid2 <- arima.diag(mymod2,std=T, plot=F)$std.resid
> resid <- arima.diag(mymod,std=T, plot=F)$std.resid  
> acf(resid,type="p")
> acf(resid2,type="p")
\end{verbatim}
The plots are very similar, but I (again) prefer the second one.
The summary of our fit is
\begin{verbatim}
$model:
$model[[1]]:
$model[[1]]$ndiff:
[1] 1

$model[[1]]$period:
[1] 12

$model[[1]]$order:
[1] 0 1 0


$model[[2]]:
$model[[2]]$order:
[1] 1 0 1

$model[[2]]$ar:
[1] 0.5485699

$model[[2]]$ndiff:
[1] 0

$model[[2]]$ma:
[1] 0.8451971

\end{verbatim}
\centerline{\psfig{figure=trin3.eps,angle=-90,width=5in}}
\centerline{\psfig{figure=trin4.eps,angle=-90,width=5in}}
\centerline{\psfig{figure=trin5.eps,angle=-90,width=5in}}


% Local Variables: 
% mode: latex
% eval: (LaTeX-math-mode t)
% eval: (my-latex-mode)
% eval: (read-abbrev-file "~/.abbrev_defs")
% abbrev-mode: t
% TeX-master: "class"
% End: 
