\documentstyle[11pt,fullpage,psfig]{article}

\input info.tex

\newcommand{\Var}{\mbox{Var}\,}

\newcommand{\eq}[1]{(\ref{#1})}

\newcommand{\combos}[2]{\left(\begin{array}{c}#1\\#2\end{array}\right)}

%\slides

\begin{document}

\begin{center}\bf
Notes for Tuesday, April 14, 1998
\end{center}



\section*{Cyclic Behavior of Time Series: the Spectrum and the Periodogram}

Many times series show cyclic fluctuations (hours, days, weeks, months,
seasons, years), but the underlying sine and cosine functions to use
to model this are usually not obvious.

One way to approach this problem is to look at the {\em discrete
  Fourier transform (DFT)} of the time series.  The sequence of
squared Fourier coefficients in the DFT is called the {\em spectrum}
of the time series.  The {\em periodogram} is an estimate of the
spectrum, based on the DFT of the observed time series, that can help
us to see what the important cycles in a time series are.  

This sort of analysis is called ``spectral analysis'', or ``looking at
the time series in the frequency domain.''  We will first look at the
spectrum and periodogram and spectrum from the point of view of
Fourier analysis, and then we will look at an interesting connection
with linear regression.  Finally we will look at some examples.

\subsection*{Short Digression: Fourier Analysis of {\boldmath
$f(t),~t\in[0,1]$}} 
 
A standard theorem from differential equations, advanced calculus, or
even a sophomore/junior linear algebra class, says that any reasonable
real function $f(t)$ with domain $[0,1]$ (or periodic function with
period 1) has a Fourier series representation,
\[
	f(t) = \sum_{k=1}^\infty a_k \cos(2\pi k t) + b_k \sin(2\pi k t)
	= \sum_{k=-\infty}^\infty c_k e^{2\pi i k t} ,
\]
where $a_k$ and $b_k$ are real coefficients, $e^{i\theta} =
\cos(\theta) + i\sin(\theta)$, $c_k = (a_k - ib_k)/2$, and $c_{-k} =
\overline{c_k} = (a_k + ib_k)/2$.  The coefficients are obtained by
$$
	a_k = 2\int_0^1 f(t) \cos(2\pi k t) dt ~~~~ \mbox{and} ~~~~
	b_k = 2\int_0^1 f(t) \sin(2\pi k t) dt
	\eqno{(*)}
$$
or equivalently $c_k = \int_0^1 f(t) \exp(2\pi i k t) dt$.  The set of
coefficients $\{a_k, b_k:~k=1,2,\ldots\}$ or equivalently
$\{c_k:~k=1,2,\ldots\}$ is called the {\em Discrete Fourier Transform
(DFT) of $f(t)$.}

Recalling that for any complex number $z=a+ib$, $z\overline z =
(a+bi)(a-bi) = a^2 +b^2 = |z|^2$, we also know that
\[
	\int_0^1 f^2(t) dt = \int_0^1 f(t)\overline{f(t)}dt
	= \int_0^1 f(t) \sum_{k=-\infty}^\infty c_k e^{-2\pi i k t}\,\,dt
	= \sum_{k=-\infty}^\infty c_k\overline{c_k}
	= \sum_{k=-\infty}^\infty |c_k|^2
	= 2\sum_{k=1}^\infty (a_k^2 + b_k^2)/4
\]
this is called {\em Parseval's identity.}  The function
\[
	S(k) = |c_k|^2 = (a_k^2 + b_k^2)/4
\]
is called the {\em spectrum} of $f(t)$; a particular estimate of
$S(k)$ from time series data, described below, is called a {\em
periodogram}.

{\bf The point:} {\em If $S(k)\approx 0$ for all $k$ beyond some $K$,
we might as well truncate the Fourier series for $f(t)$ to just $K$
terms.}


\subsection*{Fourier Analysis of a Time Series as
Linear Regression} 

This is not quite the way spectral analysis of time series is usually
developed, but it will give you the right intuition for the problem.

Consider a very simple linear regression model for a time series $y_t$
exhibiting cyclic variation with known period $p$, and overall mean
$\sum_{t=1}^n y_t/n \approx 0$:
\[
        y_t = \alpha \cos(\omega t) + \beta \sin(\omega t) + z_t ,
\] 
where $z_t$ is a white noise sequence and $\omega=2\pi/p$ is the known
frequency, i.e. [after some trig identities; try fooling around with
$e^{i(a+b)}$] we are fitting a phase-shifted sinusoid
\[
	\sqrt{\alpha^2 + \beta^2}\,\sin(\omega t + t_0)
\]
to $y_t$.  This may be rewritten in matrix form as
\[
	\left[\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n\end{array}\right]
	~=~
	\left[\begin{array}{cc}
	\cos(\omega 1) & \sin(\omega 1) \\
	\cos(\omega 2) & \sin(\omega 2) \\
	\vdots & \vdots \\
	\cos(\omega n) & \sin(\omega n)
	\end{array}\right]
	\,\,\,
	\left[\begin{array}{c} \alpha \\ \beta \end{array}\right]
	~+~
	\left[\begin{array}{c} z_1 \\ z_2 \\ \vdots \\ z_n\end{array}\right]
	~,~~~\mbox{or}~~~~
	\underline y = X\underline\beta + \underline z
\]
We may use the usual matrix formula for the regression coefficients
$\hat{\underline\beta} = (X'X)^{-1}X'\underline y$ together with some trig
integrals from Calc II, to show that
\begin{eqnarray*} 
\hat{\alpha} & \approx & \frac{2}{n}  \sum_{t=1}^n y_t\,\cos(\omega t) \\ 
\hat{\beta} & \approx & \frac{2}{n} \sum_{t=1}^n y_t\,\sin(\omega t) 
\end{eqnarray*} 
Clearly, these are discrete approximations to the integrals in $(*)$
defining a one-term DFT for $f(t)=y_{t/p}$.
% , and it may be shown that they are reasonable
% estimators of the $A_k$ and $B_k$ in the DFT of $y_t$.

Moreover, the regression sum of squares is
\[
	\sum_{t=1}^n \hat y_t^2 
	= \sum_{t=1}^n [\hat\alpha \cos(\omega t) + 
		\hat\beta \sin(\omega t)]^2 
	\approx \frac{n}{2} (\hat\alpha^2 + \hat\beta^2)
\]
(a version of Parseval's identity), so that $R^2 = \sum_t \hat y_t^2 /
\sum_t y_t^2$ is determined by the squared amplitude of the fitted
sinusoid.

This is a quite simplistic model.  More realistic would be if we not
only did not know the amplitude but also didn't know the frequency of
the series.

So consider a model of the form
\[
y_t = \sum_{k=1}^K \left \{ \alpha_k \cos(\omega_k t) + \beta_k
\sin(\omega_k t) \right \} + z_t ;
\]
we are now fitting $K$ shifted sinusoids to $y_t$.  This is again an
easy model to estimate, provided we pick the right set of frequencies.
If we take $\omega_k = 2 \pi k / n $, with $k < n/2$ then some nice
things happen:
\begin{itemize}
\item For starters, {\em this is the only set of frequencies that
we really can estimate}: smaller periods would wiggle more often than
the data, longer periods would not see an entire cycle in the data!

\item By fooling around with $x=e^{2\pi i k t/n}$ in the identity
$(1-x^n)=(1-x)(1+x+x^2+\cdots+x^{n-1})$, one can show that {\em
$e^{2\pi i k t/n}$, $k=1,2,\ldots$ represents a set of orthogonal
functions.}  That is the correlation between any pair of these
functions (measured at $t=1,2,\ldots,n$) is zero.

\item The nicest thing that happens (because of the orthogonality) is
that the estimates of $\hat{\alpha}$ and $\hat{\beta}$ are the same as
before, namely,
\begin{eqnarray*}
	\hat{\alpha}_k & = & \frac{2}{n} 
	\left \{ \sum_{t=1}^n y_t\,\cos(\omega_k t) \right \}  \\ 
	\hat{\beta}_k  & = & \frac{2}{n}
	\left \{ \sum_{t=1}^n y_t\,\sin(\omega_k t) \right \} 
\end{eqnarray*}

\item Some further manipulations show that
\[
	\sum_{t=1}^n 
		[y_t \cos(\omega_k t) + y_t\sin(\omega_k t)]^2 =
%  	\left\{ \sum{y_t \cos(\omega_k t)} \right\}^2 
%  	+ \left\{ \sum{y_t \sin(\omega_k t)} \right\}^2  =
	\frac{n}{2}(\hat{\alpha}^2_k + \hat{\beta}^2_k)
\]
so that $\hat{\alpha}_k^2 + \hat{\beta}_k^2$ again gives the sums of
squares associated with each of the $K$ fitted sinusoids.
% 
% This is quite intuitive, the terms are just the inner product of the
% data with the (basis) function.  This is what regression (and
% particularly ANOVA) was all about.  

\item Finally, (via some math, or common sense), if all the
frequencies are fit, then the regression model fits exactly.

\end{itemize}


\subsection*{The Periodogram}

Now we define the {\em periodogram function}
$I(\omega)$ for $ 0 \leq \omega \leq \pi $, by
\[
	I(\omega) = 
	\frac{1}{n}\sum_{t=1}^n [ y_t \cos(\omega t) + y_t\sin(\omega t) ]^2
% 	= \frac{1}{n}\left[ \left\{ \sum y_t \cos(\omega t) \right\}^2 +
% 		          \left\{ \sum y_t \sin(\omega t) \right\}^2 
% 		   \right ]
	.
\]
Of course $I(2\pi k/n) = (\hat{\alpha}^2_k + \hat{\beta}^2_k )/2$ is
an estimate of (twice) the spectrum $S(k) = (\hat{\alpha}^2_k +
\hat{\beta}^2_k )/4$ of the time series, and is 
% 
% The function $I(2\pi k/n) = (\alpha_k^2 + \beta_k^2)/2$ is
proportional to the squared amplitude of the sinusoid
\[
	\sqrt{\alpha_k^2 + \beta_k^2}\,\sin(\omega_k t + t_k)
\]
fitted at frequency $\omega_k = 2\pi k/n$, and represents that power
(amount of energy) that can be estimated at each of these fundamental
frequencies.  

Notice that we can ``parition'' the total sum of squares for the data
$y_t$ as
\[
	\sum_{t=1}^n y_t^2 = nI(0) + 2n \sum_{k=1}^{K_{max}} I(2 \pi k/n)
\]
where $I(0)$ corresponds to the intercept $\sum_{t=1}^n y_t$ if it is
nonzero, and $K_{max}=\lfloor n/2\rfloor$ is the largest integer less than
or equal to $n/2$.  On the one hand, this is just Parseval's identity
again.  On the other, it suggests the following ANOVA table, for
$1\leq K\leq K_{max}$, when $\sum y_t = 0$:

\begin{center}
\begin{tabular}{cr@{$~=~$}c@{$~=~$}l}
\em $k$ & \multicolumn{2}{c}{\em SS} 
	& \multicolumn{1}{c}{\em Periodogram} \\
\hline
% \\
% Intercept & $SS_0$ & $n\overline y^2$ & $nI(0)$ \\
\\
$1$	& $SS_1$ &
	$\sum_{t=1}^n [y_t \cos(\omega_1 t) +y_t \sin(\omega_1 t)]^2$ 
	& $2n\,I(2\pi/n)$ \\
%        & $\frac{n}{2}(\hat{\alpha}^2_1 + \hat{\beta}^2_1)$ \\
\\
$2$	& $SS_2$ &
	$\sum_{t=1}^n [y_t \cos(\omega_2 t) +y_t \sin(\omega_2 t)]^2$ 
	& $2n\,I(2\pi 2/n)$ \\
%        & $\frac{n}{2}(\hat{\alpha}^2_2 + \hat{\beta}^2_2)$ \\
\\
$\vdots$ & \multicolumn{1}{c}{$\vdots$} & \multicolumn{1}{l}{~}
	& $~~~~~\vdots$ \\
\\
$K$	& $SS_K$ &
	$\sum_{t=1}^n [y_t \cos(\omega_K t) +y_t \sin(\omega_K t)]^2$ 
	& $2n\,I(2\pi K/n)$ \\
%        & $\frac{n}{2}(\hat{\alpha}^2_K + \hat{\beta}^2_K)$ \\
\\
\hline
\\
Regression & $SSR$ &
%	$\sum_{t=1}^n (\hat y_t - \overline y)^2$ $~=~$ $\sum_{k=1}^K SS_k$
	$\sum_{t=1}^n \hat y_t^2$ $~=~$ $\sum_{k=1}^K SS_k$
	& $2n\,\sum_{k=1}^K I(2\pi k/n)$ \\
%	& $\frac{n}{2}\sum_{k=1}^K (\hat\alpha_k^2 + \hat\beta_k^2)$ \\
\\
\hline
\\
Error   & $SSE$ &
	$\sum_1^n (y_t - \hat y_t)^2$ $~=~$ $\sum_{k=K+1}^{n/2} SS_k$
	& $2n\sum_{k=K+1}^{n/2} I(2\pi k/n)$ \\
%	& $\frac{n}{2}\sum_{k=K+1}^{K_{max}} (\hat\alpha_k^2 + \hat\beta_k^2)$ \\
\\
\hline
\hline
\\
Total   & $SST$ &
%	$\sum_1^n (y_t - \overline y)^2$ $~=~$ $\sum_1^n y_t^2 - n(\overline y)^2$
	$\sum_1^n y_t^2$ $~=~$ $\sum_1^n \hat y_t^2 + \sum_1^n(y_t-\hat y_t)^2$
	& $2n\sum_{k=1}^K I(2\pi k/n)$ \\
%	& $\frac{n}{2}\sum_{k=1}^{K_{max}} (\hat\alpha_k^2 + \hat\beta_k^2)$ \\
\\
\hline
\end{tabular}
\end{center}
If you go back to your old linear regression notes, and examine the
matrix form of regression $\underline y = X\underline \beta +
\underline \varepsilon$ when the columns of $X$ are orthogonal, you
will see that this is the same ANOVA table: each SS is the dot product
of $\underline y$ with the corresponding column(s) of $X$, and
(because of orthogonality) it doesn't matter what order the
independent variables are entered in the regression.
 

% \subsection*{Aside: The Periodogram and the Autocovariance Function}
% 
% A more formal development of the regression analysis above says that
% if $y_t$, $t\in[0,T]$, is a random process (time series) with mean
% $E[y_t]=0$ and constant, finite variance $\mbox{Var}(y_t) = E[y_t^2]=
% \sigma^2$, then $y_t$ can be well-approximated by a truncated Fourier
% series
% \[
% 	y_t = \sum_{k=1}^{\infty} A_k \cos(\omega_k t) + B_k
% 	\sin(\omega_k t)
% 	= \sum_{k=-\infty}^\infty C_k e^{ i\omega_k t}
% \]
% where $\omega_k = 2\pi k/T$ as before and now $A_k$, $B_k$ and $C_k =
% \overline C_{-k} = (A_k - iB_k)/2$ are {\em random} coefficients,
% obtained in the same way as before, 
% \[
% 	C_k = \frac{1}{T}\int_0^T y_t e^{\omega_k t} dt .
% \]
% As random variables themselves, the $C_k$ satisfy $E[C_k]=0$ and
% \[
% 	E[C_j\overline C_k] = 
% 	\left\{
% 	\begin{array}{cl}
% 	0, & \mbox{for $k\neq j$;} \\
% 	E[|C_k|^2], & \mbox{for $k=j$}
% 	\end{array}
% 	\right.
% \]
% We begin by defining the spectrum as 
% \[
% 	S(\omega_k) = E[|C_k|^2]
% \]
% as before.  For a stationary time series $y_t$ the autocovariance
% function is
% \[
% 	\gamma(t) = E[y_s y_{s+t}]
% \]
% and is the same for any $s$.  Using the Fourier series representation
% above, it follows that
% \begin{eqnarray*}
% 	\gamma(t) & = & E[y_s y_{s+t}] ~ = ~ E[y_s\overline{y_{s+t}}] \\
% 	& = & E\left[ 
% 		\sum_{k=-\infty}^\infty C_k e^{ i\omega_k s}
% 		\sum_{k=-\infty}^\infty \overline{C_k} e^{ - i\omega_k (s+t)}
% 		\right] \\
% 	& = & E\left[
% 		\sum_{k=-\infty}^\infty |C_k|^2 e^{ - i\omega_k t}
% 		\right] \\
% 	& = & \sum_{k=-\infty}^\infty S(\omega_k) e^{ - i\omega_k t} \\
% 	& = & \sum_{k=-\infty}^\infty S(\omega_k) e^{ i\omega_k t} 
% \end{eqnarray*}
% [since $\overline{\gamma(t)}=\gamma(t)$].  It follows that
% $S(\omega_k)$ is the DFT of $\gamma(t)$.  In practice, this means that
% we can go back and forth between estimates of $\gamma(t)$ and
% estimates of $S(\omega_k)$:
% \begin{itemize}
% \item If we have estimated $\hat S(\omega_k)$, say, from a regression
% analysis like the one developed above, then we can use an ``inverse
% DFT'' to obtain an estimate $\hat\gamma(t)$ of the autocovariance.
% \item If we have estimated $\hat\gamma(t)$, then its DFT is a
% periodogram, i.e. it is an estimate $\hat S(\omega_k)$ of the
% spectrum.
% \end{itemize}
% 

\subsection*{Some Examples}


{\bf Male and Female Respiratory Deaths in the UK.} First, we re-plot
the male and female deaths in the UK from various respiratory
illnesses.

\centerline{\psfig{file=ukplots.eps,angle=-90,width=4in}}

In the male sequence, there appears to be a clear 12-month cycle, and
a hint of a cycle at about five or six months.  Not much else is
evident from these plots.

The DFT of a discrete sequence $y_t$ can be conveniently computed in
SPLUS with the so-called ``fast Fourier transform'', {\tt fft()}.
Thus, a {\em really} raw periodogram may be calculated as follows (the
{\tt Mod} function computes the modulus of a complex number):
\begin{verbatim}
402 > par(mfrow=c(4,1))
402 > length(uk.male)
[1] 72
402 > plot((1:35)/6,10/log(10)*log(Mod(fft(uk.male)))[2:36],type="l")
402 > title("Really raw periodogram")
\end{verbatim}
The values of the periodogram $I(\omega_k)$ have been plotted on a
decibel scale, $10\cdot \log_{10} S(\omega_k)$, since that is how the
periodogram is usually plotted.  

The overall mean is not zero here, so we should take it out before we
compute the periodogram.  Often, not only a mean but a linear trend is
taken out of the data before calculating the periodogram, for example:
\begin{verbatim}
402 > uk.resids _ residuals(lm(resp ~ time, data.frame(time=1:72,resp=c(uk.male))))
402 > plot((1:35)/6,10/log(10)*log(Mod(fft(uk.resids)))[2:36],type="l")
402 > title("Really raw periodogram, detrended")
\end{verbatim}
The function {\tt spectrum} in SPLUS also computes a periodogram from
the detrended time series, but in addition it smooths the observations
at the beginning and end of the series; and allows us to run one or
more moving-average smooths over the final periodogram.
\begin{verbatim}
402 > spectrum(uk.male)
402 > spectrum(uk.male,spans=c(3,3)) 
402 > # smooth with running means of length 3, twice 
\end{verbatim}
Graphs of these various periodogram estimates follow.  By convention,
all results are plotted on a decibel ($10\log_{10}$) scale.


Notice that the main power corresponds to a 12 months cycle, but there
is also significant power at some other frequencies (a half- or
one-third-year cycle and a two-year cycle).  Remember, this is like
looking at an ANOVA table to decide which terms to keep in the
regression.

Raw and smoothed periodograms of the {\tt uk.female} series are
essentially indistinguishable from these.

\centerline{\psfig{figure=periodogram.eps,width=6in}} 

{\bf Lutenizing hormone.}  We will look at the {\tt late1} series from
these data.
\begin{verbatim}
402 > ts.plot(late1)
402 > title("Lutenizing Hormone, Late1 Series")
402 > length(late1)
[1] 48
402 > l1.resids _ residuals(lm(resp~time,data.frame(resp=c(late1),time=1:48)))
402 > plot((1:23)/48,10/log(10)*log(Mod(fft(l1.resids)))[2:24],type="l")  
402 > title("Really raw periodogram, detrended")  
402 > spectrum(late1)                  
402 > spectrum(late1,span=3) # follow with a single running mean of length 3
\end{verbatim}

The plots below suggest there is a lot of power at three wavelengths,
8 = .17*48, 14 = .29*48 and 19 = .40*48, although there appears to
be ``junk'' at shorter wavelengths as well.

\centerline{\psfig{figure=periodogram2.eps,width=6in}}


\section*{Serial Dependence in Time Series}

The second popular way to analyze time series data is ``state space
analysis'' or ``looking at the time series in the space domain''.
This method concentrates on the serial dependence rather than the
cyclic behavior of the series.  The most widely used models are ARMA
and ARIMA processes.  We introduce the basic AR and MA processes here.

\subsection*{Autocovariance and autocorrelation function}

Recall that the {\em autocovariance function\/} of a stationary random
function $\{Y(t)\}$ is
\[
          \gamma(k) = \mbox{Cov}\{Y(t), Y(t-k)\} 
\]
This does not not depend on $t$ because of the the assumption of
stationarity.  Since $\gamma(0)$ is the variance, the
{\em autocorrelation function\/} is
\[
          \rho(k) = \gamma(k)/\gamma(0)  
\]


\subsection*{Basic time series models: AR and MA}

The two basic building-blocks for the ``noise'' part of a time series
model are the autoregressive (AR) and moving average (MA) processes.
MA processes are easier to describe first.

\begin{itemize}

\item A \underline{\em moving average} process of order $q$, denotes $MA(q)$,
  has the form
\[
        y_t = \mu + \beta_0 z_t + \beta_1 z_{t-1} + \ldots + \beta_q z_{t-q}
\]
or, in terms of the {\em backshift operator} $B z_t = z_{t-1}$,
\begin{eqnarray*}
        y_t-\mu  & = & (\beta_0 + \beta_1 B + \ldots \beta_q B^q) z_t \\
        & \equiv & \Theta(B) z_t,
\end{eqnarray*}
where $z_t$ is a mean-zero white
noise process.  It is easy to show that $y_t$ is always stationary,
and
\[
        \gamma(k) = \left\{
        \begin{array}{cl}
        \sum_{t=0}^{q-|k|} \beta_t \beta_{t+|k|} & |k|\leq q \\
        0 & |k|>q 
        \end{array}
        \right.
\]

Here are plots of an $MA(1)$ and an $MA(2)$ process.
% \begin{verbatim}
% 402 > ma.1 _ arima.sim(model=list(ma=.5))
% 402 > ma.2 _ arima.sim(model=list(ma=c(.2,.3)))
% 402 > ts.plot(ma.1); title("MA(1) with beta1=0.5")              
% 402 >  ts.plot(ma.2); title("MA(2) with beta1=0.2, beta2=0.3")  
% 402 > acf(ma.1)
% 402 > acf(ma.2)
% \end{verbatim}

\centerline{\psfig{figure=ma.eps,angle=-90,width=5in}}

The autocorrelation plot for a $MA(q)$ process drops to zero (except
for noise) after lag $q$, which make the {\tt acf} plot a good first
step in guessing the order $q$ of the process.

\item An \underline{\em autoregressive} process of order $p$, denoted $AR(p)$
  has the form
\[
        y_t  = \alpha_1 y_{t-1} + \alpha_2 y_{t-2} + \cdots + \alpha_p
        y_{t-p} + z_t 
\]
where again $z_t$ is a mean-zero white noise process.  In terms of the {\em
  backshift operator}, $B y_t = y_{t-1}$, we may write
\[
         (1 - \alpha_1 B - \alpha_2 B^2 - \cdots - \alpha_p B^p) y_t = z_t .
\]
If $E[y_t]=\mu\neq0$, then we would write
\[
        (1 - \alpha_1 B - \alpha_2 B^2 - \cdots - \alpha_p B^p)(y_t-\mu)=z_t .
\]
Not all values of the coefficients $\alpha_k$ lead to a stationary
series.  For example we saw last time that the random walk $y_t =
y_{t-1} + z_t$ was {\em not} a stationary process.  More generally the
complex roots of the polynomial
\[
        \Phi(x) = 1 - \alpha_1 x - \alpha_2 x^2 - \cdots - \alpha_p x^p
\]
must have all $p$ of its roots outside the unit circle in the complex
plane.  In the case of an AR(1) process, this is easy to show!

One can show that the autocovariance function for a stationary $AR(p)$
process satisfies the ``Yule-Walker'' relations
\[
        \sum_{k=1}^p \gamma(k-i)\alpha_k  =  \gamma(k) 
                \mbox{~~~ $i=1,2,\ldots,p$} ,
\]
We can use the Yule-Walker equations to solve for $\gamma(k)$ in terms
of the $\alpha$'s, if they are known, or solve for the $\alpha$'s in
terms of the $\gamma(k)$'s.  This suggests a simple way to estimate
the coefficients of an $AR$ process.  Also,
\[
        \sigma^2_y  =  \sum_{k=1}^p \gamma(k)\alpha_k + \sigma^2_z .
\]
As an exercise, show that the autocorrelation function of an $AR(1)$
process $y_t = \alpha y_{t-1} + z_t$ satisfies $\rho_k = \alpha
\rho_{k-1}$, and so, inserting $\rho_0 = 1 $, we get $\rho_k =
\alpha^k$.  What does $\rho_k$ look like for an $AR(2)$ process?
Here are plots of simulated $AR(1)$ and $AR(2)$ processes.

% \begin{verbatim}
% 402 > par(mfrow=c(2,2))
% 402 > ar.1 _ arima.sim(model=list(ar=.5))
% 402 > ts.plot(ar.1); title("AR(1) with alpha1=0.5")
% 402 > ar.2 _ arima.sim(model=list(ar=c(.2,.3)))
% 402 > ts.plot(ar.2); title("AR(2) with alpha1=0.2, alpha2=0.3")
% 402 > acf(ar.1)
% 402 > acf(ar.2)
% \end{verbatim}

\centerline{\psfig{figure=ar.eps,angle=-90,height=3.5in}}

The {\tt acf} plot of an $AR(1)$ decays exponentially to zero,
possibly with some oscillation in sign.  The $AR(2)$ {\tt acf} plot
also decays to zero (this gradual decay in the {\tt acf} plot
distinguishes pure-AR from pure-MA processes) but it isn't easy to
tell from the plot what the order $p$ of the process should be.

\end{itemize}

\subsection*{Partial Autocorrelation Function}

Here are some more examples of AR processes and their acf's:

\begin{itemize}

\item The following graph shows plots of two AR processes,
\[
y_t = .7 y_{t-1} + .2 y_{t-2} + z_t 
\]
and
\[
y_t = .7 y_{t-1} + z_t 
\]
and their autocorrelation functions.  Notice that the ACF's have long
lags and it is not easy to diagnose anything.  
\begin{verbatim}
402 > tsplot(dd <- arima.sim(200, model=list(ar=c(.7, .2))))
402 > acf(dd)
402 > tsplot(dd <- arima.sim(200, model=list(ar=c(.7))))    
402 > acf(dd)                                           
\end{verbatim}
\centerline{\psfig{figure=4acf.eps,angle=-90,width=6in}}

\item Next we plot two more MA processes and their acf's:
\begin{verbatim}
402 > tsplot(dd <- arima.sim(200, model=list(ma=c(.6))))
402 > acf(dd)
402 > tsplot(dd <- arima.sim(200, model=list(ma=c(.6, .4))))
402 > acf(dd)
\end{verbatim}
\centerline{\psfig{figure=4acfma.eps,angle=-90,width=6in}}

The two models that we fit are
\[
y_t = .6 z_{t-1} + z_t 
\]
and
\[
y_t = .6 z_{t-1} + .4 z_{t-2} + z_t 
\]

Just as before, that the ACF tails off after lag 1 and 2.  That
is the ACF provides a diagnostic of the order of the MA process.
\end{itemize}
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.

The form of an autoregressive model,
\[
        y_t  = \alpha_1 y_{t-1} + \alpha_2 y_{t-2} + \cdots + \alpha_p
        y_{t-p} + z_t 
\]
suggests that instead of plotting the covariances between lagged
observations, we should be plotting the coefficients $\alpha_k$
themselves: essentially, we want to fit AR models of orders 1, 2, 3,
and so on.  Then at lag 1 we plot the coefficient $\alpha_1$ of
the fitted AR(1) model.  At lag 2 we plot the coefficient $\alpha_2$ of
the AR(2) model, and so on.  When $p$ is too high, the plotted
coefficients should be essentially zero.  This type of plot is called
a partial autocorrelation function (PACF), and it works really well.  

To actually calculate the PACF (and plot it in SPLUS) requires the use
of the acf function with an extra argument.
\begin{verbatim}
acf(dd,type="p")
\end{verbatim}
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. 
\begin{verbatim}
402 > acf(arima.sim(200, model=list(ar=c(.7, .2))), type="p") 
402 > acf(arima.sim(200, model=list(ar=c(.7, .2))), type="p")
402 > acf(arima.sim(200, model=list(ar=c(.7))), type="p")    
402 > acf(arima.sim(200, model=list(ar=c(.7))), type="p")
\end{verbatim}
\centerline{\psfig{figure=4pacfar.eps,angle=-90,width=6in}}

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
\begin{verbatim}
402 > acf( arima.sim(200, model=list(ma=c(-.9))), type="p")
402 > acf( arima.sim(200, model=list(ma=c(.9))), type="p") 
402 > acf( arima.sim(200, model=list(ma=c(.9, .6))), type="p")
402 > acf( arima.sim(200, model=list(ma=c(.9, .6))), type="p")
\end{verbatim}
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.

\centerline{\psfig{figure=4pacfma.eps,angle=-90,width=6in}}

Thus,
\begin{itemize}
\item To get an idea what the order $q$ of a pure $MA(q)$ process
should be, use an ACF plot.
\item To get an idea what the order $p$ of a pure $AR(p)$ process
should be, use a PACF plot.
\end{itemize}
There are better ways to choose $p$ and $q$, but these make useful
illustrations and can help you diagnose whether you are ``way off''.

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



