Partial Autocorrelation Function

Recall, in the dim past, that we saw the autocorrelation function (acf) for an AR(1) model looked just like an exponential decline (or an oscillatory exponential decline). What does the ACF of an AR(2) or AR(3) process look like? The following graph shows plots of two AR processes,

y_{t} = .7y_{t−1} + .2y_{t−2} + z_{t}

and

y_{t} = .7y_{t−1} + z_{t}

and their autocorrelation functions. Notice that the ACF’s have long lags and it is not easy to diagnose anything.

> tsplot(dd <- arima.sim(200, model=list(ar=c(.7, .2))))
> acf(dd)
> tsplot(dd <- arima.sim(200, model=list(ar=c(.7))))
> acf(dd)
>

What about a similar set of plots for a MA(1) and MA(2) process?

> tsplot(dd <- arima.sim(200, model=list(ma=c(.6))))
> acf(dd)
> tsplot(dd <- arima.sim(200, model=list(ma=c(.6, .4))))
> acf(dd)
>

The two models that we fit are

y_{t} = .6z_{t−1} + z_{t}

and

y_{t} = .6z_{t−1} + .4z_{t−2} + z_{t}

Notice that the ACF tails off after lag 1 and 2. That is the ACF provides a diagnostic of the order of the MA process. A simple (too simple!) conclusion is that if the ACF is zero beyond a certain lag, then that is evidence for an MA process of the appropriate order.

That turns out to be a useful observation, but it doesn’t help us with AR processes. We need to search for a diagnostic that will help us determine the order of an AR process—but we don’t have to search all that far. What we need is the partial autocorrelation function (PACF). The conceptual idea (not computational, but that doesn’t matter) is that we take a time series and fit AR models of orders 1, 2, 3, and so on. For the PACF at lag 1 we use the coefficient α_{1} of the fitted AR(1) model. At lag 2 we use the coefficient α_{2} of the AR(2) model, and so on. Clearly (I hope it is clear) if we really have an AR(k) model then all of the coefficients at lags greater than k should be zero—or at least they should be as close to zero as sampling error allows. It works really well. To actually calculate the PACF (and plot it in S) requires the use of the acf function with an extra argument.

acf(dd,type="p")

I have produced four plots, two of AR(2) processes and two of AR(1) processes. Notice that these plots look like the ACF plots of MA processes.

> acf(arima.sim(200, model=list(ar=c(.7, .2))), type="p")
> acf(arima.sim(200, model=list(ar=c(.7, .2))), type="p")
> acf(arima.sim(200, model=list(ar=c(.7))), type="p")
> acf(arima.sim(200, model=list(ar=c(.7))), type="p")

Maybe the PACF of a MA process will look like the ACF of an AR process? This sort of duality would be quite useful. Lets look and see. The simulations were generated by

> acf( arima.sim(200, model=list(ma=c(-.9))), type="p")
> acf( arima.sim(200, model=list(ma=c(.9))), type="p")
> acf( arima.sim(200, model=list(ma=c(.9, .6))), type="p")
> acf( arima.sim(200, model=list(ma=c(.9, .6))), type="p")

The following plots really do look like the ACF of AR’s. There is a mathematical duality between ACF’s and PACF’s, but the important thing to note is that each of them is quite good at helping to diagnose the order of a pure MA or a pure AR process.

241

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

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")

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)

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)

The data we with to look at have the following first few lines

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

We read in the data as (with skip to skip the first line of rubbish).

trinity <- as.ts( scan("trinity.dat", skip=1))

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

> tsplot(trinity)

> title("tsplot(trinity)")
> tsplot(trinity[1:96])
> title("tsplot(trinity[1:96])")
> acf(trinity)
> acf(trinity, type="p")

We have a lot of reason to think about 12th order difference of the data, so lets try that.

> 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")

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.

>  mymod <- arima.mle(trinity, model=list( list(ndiff=1, period=12),
+ list(order=c(0,0,2))))
> mymod$loglik
[1] 4271.951

Also try the same except an ARIMA(1,0,1) model.

>  mymod2 <- arima.mle(trinity, model=list( list(ndiff=1, period=12),

+ list(order=c(1,0,1))))
 > mymod2$loglik
[1] 4256.779

Just for kicks we’ll try the raw AR(1) and raw MA(1) and see how they do.

> 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

Either of these is demonstrably works (likelihood ratio test) than the combined models. Lets quickly look at the diagnostics for the two candidate models.

> arima.diag(mymod)
> arima.diag(mymod2)

I could not not figure out how to automatically get the pacfs, so here is the manual approach to that

> 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")

The plots are very similar, but I (again) prefer the second one. The summary of our fit is

$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

254 255