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

\parskip=0.5em
\parindent=0pt

\frenchspacing

\long\def\comment#1{}

\let\bar\overline

\def\tp{^T}

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

\def\Var{\mbox{Var}\,}

\def\Varsigma{{\mathbf{\not{\hspace*{-0.33em}\Sigma}}}}
\def\tr{\mbox{tr\,}}

\def\diag{\mbox{\it diag\,}}

\let\vc\Varsigma
\let\phi\varphi

\def\Corr{\mbox{Corr\,}}

\let\ul\mb
\let\by\times

\def\ul#1{\underline{#1}}
\def\ulem#1{\ul{\em #1}}

% \def\mb#1{{\mathbf #1}}
\let\mb\ul

\def\Re{{\cal R}}

\newenvironment{boxed}{\begin{center}\begin{tabular}{|p{5.5in}|}
\hline}{\\ \hline\end{tabular}\end{center}}

\newtheorem{definition}{Definition}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\newtheorem{lemma}{Lemma}

\begin{document}

\begin{center}\Large\bf 
        Generalized Least Squares and Time Series Regression
\end{center}

% This material follows fairly closely Fox Section 14.3. 



\section{Generalized Least Squares}

A simple and natural extension of weighted least squares (WLS) is {\em
  generalized least-squares (GLS) regression}.  Formally, GLS looks
  exactly like WLS:

Suppose that the model is 
\[
        \mb y = X \mb \beta + \mb \epsilon
\]
where
\[
        \mb \epsilon_{n\times 1} \sim N_n(\mb 0,\vc)
\]
i.e. now the random vector $\mb \epsilon$ is multivariate normal with
mean $\mb 0$ and variance-covariance matrix $\Var(\mb \epsilon) =
\vc$.  

If we know the exact value of $\vc_{n\times n}$, we can compute a
square root $\Gamma_{n\times n}$ of $\vc^{-1}$, as we have seen many
times before:
\[
        \vc^{-1} = \Gamma\tp \Gamma
\]
For example we can compute $\Gamma$ using the eigenvector-eigenvalue
decomposition of $\vc$.  For computational reasons it is usually
better to choose $\Gamma$ to be upper-triangular or lower-triangular;
this can be done in Splus using the function {\tt chol()}.  Once we
have $\Gamma$, we can immediately reduce to the standard model by
considering instead
\[
        \mb y^* = \Gamma\mb y = \Gamma X\mb \beta + \Gamma\epsilon
        = X^* \mb \beta + \epsilon^*
\]
Note that in our discussion of WLS we simply calculated $\Gamma$
directly as
\[
        \Gamma = W^{1/2} = \left[\begin{array}{cccc}
            w_i & 0 & \cdots & 0 \\
            0 & w_2 & \cdots & 0 \\
            \vdots & \vdots & \ddots & \vdots \\
            0 & 0       & \cdots & w_n
          \end{array}\right]
\]
We can immediately apply our usual least-squares or normal-theory MLE
machinery to show that the MLE's are
\[
\begin{array}{rcccccl}
        \hat{\mb \beta}_{GLS} & = & ({X^*}\tp X^*)^{-1}X^* \mb y^* 
        & = & (X\tp \Gamma\tp\Gamma X)^{-1}X\tp \Gamma\tp\Gamma \mb y 
        & = & (X\tp \vc^{-1} X)^{-1}X\tp \vc^{-1} \mb y \\
\end{array}
\]
Once again, we could also obtain these results by directly
differentiating the log of the likelihood
\begin{eqnarray}
         L(\ul y; \ul\beta,\vc) 
        & = &
         \frac{1}{(2\pi)^{n/2}|\vc|^{1/2}}
                \exp \left[ - \frac12 (\ul y - X\ul \beta)\tp
                \vc^{-1} (\ul y - X\ul\beta) \right] 
        \label{gls-lik} \\
        & = &
         \frac{1}{(2\pi)^{n/2}|\vc|^{1/2}}
                \exp \left[ - \frac12 (\Gamma[\ul y - X\ul \beta])\tp
                (\Gamma[\ul y - X\ul\beta]) \right]
        \nonumber \\
        & = &
         \frac{1}{(2\pi)^{n/2}|\vc|^{1/2}}
                \exp \left[ - \frac12 (\ul y^* - X^*\ul \beta)\tp
                (\ul y^* - X^*\ul\beta) \right]
        \nonumber
\end{eqnarray}
As before, the solutions $\hat{\mb \beta}_{GLS}$ are obtained by
minimizing the ``weighted'' squared raw residuals
\[
        Q ~ = ~ \sum_i (y_i^* - \hat y_i^*)^2 
        ~ = ~ (\ul y^* - X^*\ul \beta)\tp
                (\ul y^* - X^*\ul\beta)
        ~ = ~ (\ul y - X\ul \beta)\tp
                \vc^{-1} (\ul y - X\ul\beta)
\]
If we really have the right $\vc$, it is easy to show that
\begin{equation}
\left.
\begin{array}{rcl}
        E[\hat{\mb\beta}_{GLS}|X] & = & \mb\beta , \\
%
        \hat{\rm Var}\,(\hat{\mb\beta}_{GLS}) & = & 
        (X\tp \vc^{-1} X)^{-1} ,
\end{array}
\right\}
\label{GLS-stuff}
\end{equation}
that the $\hat{\mb
\beta}_{GLS}$ are BLUE (Gauss-Markov theorem), etc.

In practice we have to either
\begin{itemize}
\item Know $ \vc$ in advance; or
\item Estimate $\vc$ from the data.
\end{itemize}
Of course, knowing $\vc$ in advance is unlikely for any realistic
problem.  To estimate $\vc$ from the data, we can precompute an
estimate of $\vc$ from the OLS residuals (which has the usual
multiple-looks-at-the-data problems) or we can try to maximize the
likelihood \eq{gls-lik} simultaneously in $\mb\beta$ and $\vc$.

For the latter approach to work, $\vc$ cannot be completely general,
because $\vc$ would then have $n(n-1)/2$ parameters to estimate, which
grows like the square of the sample size!  Instead we need to be able
to assert something about the form of $\vc$.  One very common and
useful form is obtained by assuming that the errors $\mb\epsilon$ form
a simple {\em time series}, and so we now digress to introduce some
simple time series ideas.




\section{Serially Correlated Errors --- Time Series}

A {\em time series} is simply a sequence of observations indexed by
time.  The errors $\epsilon_i$ in a linear regression model may or may
not be indexed by time (depends on the application of course!) but the
simple kinds of correlation structures that arise in time series
models are useful as models of correlation structure in linear
regression as well.  

The basic idea of time series correlation structure is that
$\epsilon_i$ and $\epsilon_j$ will be highly correlated if $i$ and $j$
are close together, and will have little correlation if $i$ and $j$
are far apart.  Clearly this kind of structure works for other
interpretations of $i$ and $j$ than points in time.  For example, if
$i$ and $j$ represent points in space, then we are dealing with an
area of statistics called {\em spatial statistics}; see for example
Venables and Ripley (1997, Chapter 16).

In this section we introduce a few simple building blocks of time
series models that are useful for building correlation models in
linear regression.  In all cases we are dealing with the time series
$
        \epsilon_1, \epsilon_2, \ldots
$
that arises as errors in linear regression.  In particular we will
assume
\begin{itemize}
\item $E[\epsilon_i] = 0$ for all $i$
\item $\Var(\epsilon_i) = \sigma^2_\epsilon$ for all $i$
\end{itemize}
and we are mostly interested in learning about the correlations
\[
        \mbox{Corr}\,(\epsilon_i,\epsilon_j)
        = 
        \mbox{Cov}\,(\epsilon_i, \epsilon_j)/\sigma^2_\epsilon
\]

\subsection{Stationary Errors and the Autocovariance and
Autocorrelation Functions} 

The error time series $\epsilon_1, \epsilon_2, \ldots$ is {\em
stationary}\footnote{This is ``second-order stationarity'', which is
the same as the time series notion of ``strong stationarity if the
$\epsilon_i$'s are jointly normal.} if the covariance between any two
observations
\[
        \mbox{Cov}\,(\epsilon_i, \epsilon_j)  = E[\epsilon_i\epsilon_j]
\]
depends only on the difference $k=|i-j|$ and not on the particular
locations $i$ and $j$.

For a stationary time series, the {\em autocovariance function\/} is
\[
          \gamma(k) = \mbox{Cov}\,(\epsilon_i, \epsilon_{i-k})
\]
Note that $\sigma^2_\epsilon=\Var(\epsilon_i) = \gamma(0)$.  The {\em
autocorrelation function\/} is defined to be
\[
          \rho(k) = \gamma(k)/\gamma(0)
          \equiv\mbox{Corr}\,(\epsilon_i, \epsilon_{i-k})
\]
A plot of $x=k$ vs $y=\rho(k)$ is called an {\em autocorrelation
function plot} (ACF plot) and is a basic tool in displaying the serial
dependence in time series data.  (We'll see examples below).



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

The two basic building-blocks for time series models 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$, denoted $MA(q)$,
  has the form
\[
        \epsilon_i = \mu + \beta_0 z_i + \beta_1 z_{i-1} + \ldots + \beta_q z_{i-q}
\]
or, in terms of the {\em backshift operator} $B z_i = z_{i-1}$,
\begin{eqnarray*}
        \epsilon_i  & = & (\beta_0 + \beta_1 B + \ldots \beta_q B^q) z_i \\
        & \equiv & \Theta(B) z_i,
\end{eqnarray*}
where the $z_i$ are i.i.d. $N(0,\sigma^2_z)$.  (``white noise'' in time
series parlance).  One can show that $\epsilon_i$ is always
stationary, and
\[
        \gamma(k) = \left\{
        \begin{array}{cl}
        \sum_{i=0}^{q-|k|} \beta_i \beta_{i+|k|} & |k|\leq q \\
        0 & |k|>q 
        \end{array}
        \right.
\]

Figure~\ref{ma} gives plots of a simulated $MA(1)$ and a simulated
$MA(2)$ process.
\begin{verbatim}
> ma.1 _ arima.sim(model=list(ma=.5))
> ma.2 _ arima.sim(model=list(ma=c(.2,.3)))
> ts.plot(ma.1); title("MA(1) with beta1=0.5")              
> ts.plot(ma.2); title("MA(2) with beta1=0.2, beta2=0.3")  
> acf(ma.1)
> acf(ma.2)
\end{verbatim}

\begin{figure}
\centerline{\psfig{figure=ma.eps,angle=-90,width=5in}}
\caption{An MA(1) process and an MA(2) process, with corresponding
autocorrelation functions.}\label{ma}
\end{figure}

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
\[
        \epsilon_i  = \alpha_1 \epsilon_{i-1} + \alpha_2 \epsilon_{i-2} + \cdots + \alpha_p
        \epsilon_{i-p} + z_i 
\]
where again $z_i$ are i.i.d. $N(0,\sigma^2_z)$.  In terms of the {\em
  backshift operator}, $B \epsilon_i = \epsilon_{i-1}$, we may write
\[
         (1 - \alpha_1 B - \alpha_2 B^2 - \cdots - \alpha_p B^p) \epsilon_i = z_i .
\]
Not all values of the coefficients $\alpha_k$ lead to a stationary
series.  For example the process $\epsilon_i = \epsilon_{i-1} + z_i$
{\em not} a stationary process, since this would imply
$\sigma^2_\epsilon = \sigma^2_\epsilon + \sigma^2_z$, which is
impossible unless $\sigma^2_z=0$. 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.  Also,
\[
        \sigma^2_\epsilon  =  \sum_{k=1}^p \gamma(k)\alpha_k + \sigma^2_z .
\]
As an exercise, show that the autocorrelation function of an $AR(1)$
process $\epsilon_i = \alpha \epsilon_{i-1} + z_i$ 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?
Figure~\ref{ar} contains plots of simulated $AR(1)$ and $AR(2)$
processes, and the corresponding autocorrelation functions.

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

\begin{figure}
\centerline{\psfig{figure=ar.eps,angle=-90,height=3.5in}}
\caption{An AR(1) process and an AR(2) process, with corresponding
autocorrelation functions.}\label{ar}
\end{figure}

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.  

There is another kind of time series plot called a {\em partial
autocorrelation (PACF) plot}, which is basically the graph of the
coefficients $\alpha_i$ in an $AR(p_{MAX})$ model where $p_{MAX}$ is
much larger than the true order of the process, that can be used to
graphically guess the order of an AR process.  In a partial
autocorrelation plot, the plotted values of $\alpha_i$ fall to zero
when $i>p$, the true order of the $AR(p)$ process.  You can get a PACF
plot in Splus by using {\tt acf(...,type="partial")}.

\item For completeness I also mention the \underline{\em ARMA} process of orders $p$ and $q$, denoted
$ARMA(p,q)$, which simply combines the structure of an $AR(p)$ and an $MA(q)$
model. This is easiest to express in terms of the backshift operator:
An $ARMA(p,q)$ process satisifes 
$
         \Phi(B) \epsilon_i = \Theta(B) z_i,
$
where again $z_i$ are i.i.d. $N(0,\sigma^2)$, and $\Phi(B)$ and $\Theta(B)$ are
polynomials of the same form as indicated above.  

\end{itemize}



\subsection{The Variance-Covariance Matrix of an AR(1) Process}

Errors $\epsilon_1, \ldots, \epsilon_n$ that follow an $ARMA(p,q)$
process generally have a fairly simple correlation structure, with
only about $p+q$ parameters in $\vc$, instead of the full $n(n-1)/2$
parameters in an unstructured $\vc$.  

As an example, we note that, from the general descriptions of AR and MA
processes above, it is easy to show that, for an $AR(1)$ process
\begin{itemize}
\item $\epsilon_i =  \alpha \epsilon_{i-1} +  z_i$;
\item $\sigma^2_\epsilon = \sigma^2_z / (1-\alpha^2)$;
\item $\mbox{Corr}\,(\epsilon_i,\epsilon_{i-k}) = \alpha^k$.  
\end{itemize}
and thus the form of the covariance matrix is
\[
        \vc = \mbox{Var}\,({\mb \epsilon}) 
        = 
        \sigma^2_z/(1-\alpha^2)
        \left[\begin{array}{ccccc}
        1 & \alpha^1 & \alpha^2 & \cdots & \alpha^{n-1} \\
        \alpha^1 & 1 & \alpha^1 & \cdots & \alpha^{n-2} \\
        \alpha^2 & \alpha^1 & 1 & \cdots & \alpha^{n-3} \\
        \vdots & \vdots & \vdots & \ddots & \vdots \\
        \alpha^{n-1} & \alpha^{n-2} & \alpha^{n-2} & \cdots & 1
        \end{array}\right]
        = 
        \sigma^2_\epsilon
        \left[\begin{array}{ccccc}
        1 & \rho & \rho^2 & \cdots & \rho^{n-1} \\
        \rho & 1 & \rho & \cdots & \rho^{n-2} \\
        \rho^2 & \rho & 1 & \cdots & \rho^{n-3} \\
        \vdots & \vdots & \vdots & \ddots & \vdots \\
        \rho^{n-1} & \rho^{n-2} & \rho^{n-2} & \cdots & 1
        \end{array}\right]
\]
where we are letting $\rho=\alpha=\mbox{Corr}\,(\epsilon_i,\epsilon_{i-1})$.


\section{Computation, for the AR(1) errors model} 

Now let us consider estimation for 
\[
        \mb y = X \mb \beta + \mb \epsilon, ~~~ \mb\epsilon \sim N(\mb0,\vc)
\]
where the errors have AR(1) structure
\[
        \vc ~ = ~ \sigma^2_z/(1-\rho^2)
        \left[\begin{array}{ccccc}
        1 & \rho & \rho^2 & \cdots & \rho^{n-1} \\
        \rho & 1 & \rho & \cdots & \rho^{n-2} \\
        \rho^2 & \rho & 1 & \cdots & \rho^{n-3} \\
        \vdots & \vdots & \vdots & \ddots & \vdots \\
        \rho^{n-1} & \rho^{n-2} & \rho^{n-2} & \cdots & 1
        \end{array}\right]
        ~ \equiv ~ \sigma^2_\epsilon ~ R
\]
In this case it is easy to verify that a lower-triangular square root
of $R^{-1}$ is 
\begin{equation}
        \Gamma = 
        \left[\begin{array}{ccccc}
        \sqrt{1-\rho^2} & 0 & 0 & \cdots & 0 \\
        -\rho & 1 & 0 & \cdots & 0 \\
        0  & -\rho & 1 & \cdots & 0 \\
        \vdots & \vdots & \vdots & \ddots & \vdots \\
        0 & 0 & 0 & \cdots & 1
        \end{array}\right]
        \label{rootR}
\end{equation}
and we can immediately apply OLS to the transformed equation
\[
        \mb y^* = X^*\mb\beta + \mb\epsilon
\]
and use the usual OLS estimates of $\mb\beta$, $\sigma^2_\epsilon$,
etc.  Thus all we need to know is $\rho$.  There are basically two
approaches:
\begin{itemize}
\item Fit the uncorrected OLS regression and estimate $\rho$ directly
from the equation 
\begin{equation}
        \hat\rho = \frac{\sum_{i=2}^n E_iE_{i-1}}{\sum_{i=1}^n E_i^2}
        \label{corr-est}
\end{equation}
and then plug this into the estimation framework just outlined.

\item Observe that the general likelihood \eq{gls-lik} simplifies, in
the case of AR(1) structured errors, to the log-likelihood
\begin{eqnarray*}
         \log L(\ul y; \ul\beta,\vc) 
        & = &
         \log\left\{\frac{1}{(2\pi)^{n/2}|\vc|^{1/2}}
                \exp \left[ - \frac12 (\ul y - X\ul \beta)\tp
                \vc^{-1} (\ul y - X\ul\beta) \right] \right\} \\
        & = & 
        -\frac{n}2 \log 2\pi - \frac{n}2 \log \sigma^2_z - \frac12
        \log (1-\rho^2) - \frac{1}{2\sigma^2_z} (\mb y^* - X^*\mb\beta)\tp
        (\mb y^* - X^*\mb\beta)
\end{eqnarray*}
(where $\mb y^* = \Gamma \mb y$ and $X^* = \Gamma X$ and $\Gamma$ is
from equation \ref{rootR});
and this can be maximized directly (we will give an example at the
end of the next section).  
\end{itemize}

\section{An example}

As an example we consider 
% Fox's (1997, Ch. 14) adjustment of the
data of Fox and Hartnagel (1979, {\em Can. Rev. of Sociology and
  Anthropology, 16,} 96--104) on yearly theft-conviction rates per
100,000 Canadian women aged 15 years and older, in the years
1935--1968.  The variables are
\begin{verbatim}
Year          1935--1968
Fert          Total Fertility Rate per 1000
Partic        Labour-Force Participation Rate per 1000
Coll          Post-Secondary Degree Rate per 10,000
FemAllConv    Female Indictable-Offense Conviction Rate per 100,000
FemTheftConv  Female Theft Conviction Rate per 100,000
MalAllConv    Male Indictable-Offense Conviction Rate per 100,000
MalTheftConv  Male Theft Conviction Rate per 100,000
\end{verbatim}
Note that there are some missing values in this data, and Fox's codes
``?'' have to be changed to ``NA'' in order to identify (and for the
purposes of illustration delete) these rows of the data matrix.
Once this is done, we read the data into Splus in the usual way:

\begin{verbatim}
> hart _ read.table("hartnagl.dat",na.strings="?")
> dim(hart)
[1] 38  8
> names(hart) _ scan(what="")
1: Year
2: Fert
3: Partic
4: Coll
5: FemAllConv
6: FemTheftConv
7: MalAllConv
8: MalTheftConv
9: 
> hart
   Year Fert Partic Coll FemAllConv FemTheftConv MalAllConv MalTheftConv 
 1 1931 3200    234 12.4       77.1           NA      778.7           NA
 2 1932 3084    234 12.9       92.9           NA      745.7           NA
 3 1933 2864    235 13.9       98.3           NA      768.3           NA
 4 1934 2803    237 13.6       88.1           NA      733.6           NA
 5 1935 2755    238 13.2       79.4         20.4      765.7        247.1
 6 1936 2696    240 13.2       91.0         22.1      816.5        254.9
 7 1937 2646    241 12.2      100.4         22.4      821.8        272.4
 8 1938 2701    242 12.6      108.9         21.8      956.8        285.8
 9 1939 2654    244 12.3      123.6         21.1     1035.7        292.2
10 1040 2766    245 12.0      157.3         21.4      951.6        256.0
11 1941 2832    246 11.7      154.3         25.3      850.9        205.8
12 1942 2964    268 11.2      143.9         27.1      769.7        188.0
13 1943 3041    333 11.5      147.5         29.0      811.2        205.8
14 1944 3010    335 11.1       97.3         24.2      865.4        207.9
15 1945 3018    331 12.5       76.6         24.7      866.8        197.8
16 1946 3374    253 15.0       72.8         20.5      968.0        195.9
17 1947 3545    243 17.6       68.9         20.7      894.0        198.9
18 1948 3441    241 21.2       66.7         22.8      830.8        198.6
19 1949 3456    242 22.7       55.0         18.5      811.2        216.1
20 1950 3455    237 21.4       55.0         19.3      775.4        212.0
21 1951 3503    242 20.7       55.0         20.0      739.8        208.0
22 1952 3641    240 19.5       54.7         19.4      740.3        199.3
23 1953 3721    233 36.6       47.9         16.1      725.1        176.4
24 1954 3828    232 18.1       49.0         16.5      739.9        166.8
25 1955 3831    236 17.7       44.4         15.2      663.8        179.7
26 1956 3858    245 18.7       40.2         15.3      633.7        198.8
27 1957 3925    256 20.8       43.8         17.1      716.5        262.5
28 1958 3880    261 22.7       49.5         21.2      761.0        263.4
29 1959 3935    265 24.6       50.4         22.3      665.3        253.6
30 1960 3895    278 28.5       59.6         28.5      742.6        280.9
31 1961 3840    291 31.3       72.4         34.4      789.8        290.9
32 1962 3756    290 38.0       70.6         36.0      776.8        274.7
33 1963 3669    290 42.0       82.0         43.8      842.2        296.9
34 1964 3502    307 48.5       89.5         49.7      799.1        281.3
35 1965 3145    313 59.2      101.3         58.6      763.0        264.2
36 1966 2812    325 69.7      116.7         71.4      804.6        286.6
37 1967 2586    339 80.4      115.2         70.6      781.1        272.0
38 1968 2441    338 90.4      122.9         73.0      849.7        274.7
> hart _ hart[!is.na(apply(hart,1,sum)),]
\end{verbatim}

The raw $y$ data (female theft conviction rates) is shown in
Figure~\ref{yts}.  One can see apparently-strong serial dependence in
the data, and also some change in the trend of the data after about 1960.
\begin{verbatim}
> attach(hart)
> ts.plot(FemTheftConv,axes=F,xlab="Year",ylab="Fem Theft Conv Rate")
Warning messages:
  Non time series data in: ts.plot(FemTheftConv, axes = F, xlab =
        "Year", ylab = "Fem Theft Conv Rate")
> axis(1,at=seq(1,1+40,by=5),label=seq(1935,1935+40,by=5))
> axis(2)
> box()
\end{verbatim}

\begin{figure}
\centerline{\psfig{figure=yts.ps,height=3.5in}}
\caption{Canadian women's theft-per-100,000 rate from the Fox and
Hartnagel data.}\label{yts}
\end{figure}

Now we will
\begin{itemize}
\item Do an OLS analysis and look for problems in the residuals
\item Do the plug-in estimate of $\rho$ for an AR(1) model of the
residuals
\item Do a joint maximization using an Splus function {\tt nlmin()} that
is designed to minimize fairly general nonlinear functions.
\end{itemize}

\subsection{The OLS analysis}

We regress female theft conviction rates on fertility rates,
participation in the labor force, and post-secondary degree rates, to
see if these things account for the rate at which women are convicted
for (and hence participate in) thefts.  In addition we regress on the
male theft conviction rate, as a stand-in for other environmental
factors that are not captured by the other three variables.

\begin{verbatim}
> summary(olsfit _ lm(FemTheftConv ~ Fert + Partic + Coll + MalTheftConv))

Call: lm(formula = FemTheftConv ~ Fert + Partic + Coll + MalTheftConv)
Residuals:
    Min     1Q  Median    3Q   Max 
 -8.974 -1.945 -0.4133 2.527 7.165

Coefficients:
                Value Std. Error  t value Pr(>|t|) 
 (Intercept)  -7.3341   9.4379    -0.7771   0.4434
        Fert  -0.0061   0.0014    -4.2017   0.0002
      Partic   0.1199   0.0234     5.1236   0.0000
        Coll   0.5516   0.0433    12.7525   0.0000
MalTheftConv   0.0393   0.0186     2.1193   0.0428

[some output omitted]

> par(mfrow=c(2,3))
> plot(olsfit)
> attach("../week10")
> plot(fitted(olsfit),residuals(olsfit))
> abline(h=0)
> qqenv(residuals(olsfit))
> plot(35:68,residuals(olsfit),type="l")
> abline(h=0)
> print(acf(residuals(olsfit)))
$acf:

, , 1
               [,1] 
 [1,]  1.0000000000
 [2,]  0.2441946566
 [3,] -0.1919368953
 [4,] -0.2654978931
 [5,] -0.0002950965
 [6,]  0.0247247219

[some output omitted]
\end{verbatim}

The ``usual'' residual diagnostics shown in the first two panels of
Figure~\ref{olsdiag} don't show much.  However, a time-series plot of
the residual does appear to reveal some serial dependence (though with
only 35 data points it is a bit hard to say).  The dependence is not
terribly strong, as the autocorrelation plot in Figure~\ref{olsdiag}
doesn't show any autocorrelations outside approximate 95\% testing
bands for $H_0$: correlation equals zero\footnote{These testing bands
are somewhat crude however and a better approach might be to use the
Durbin-Watson style test for nonzero autocorrelation mentioned in 
Rawlings, Pantula and Dickey (1998, pp. 354--355).}.
% Fox
% (1997, pp. 382ff.).}.

From the Splus code above, you can see that the {\tt acf()} function
also prints out the autocorrelations at lags 1, 2, \ldots.  In
particular we can deduce from the output that
\[
        \hat\rho = \widehat{\mbox{Cov}}(\epsilon_i,\epsilon_{i-1})
        = 0.2442
\]

\begin{figure}
\centerline{\psfig{figure=olsdiag.ps,height=3.5in}}
\caption{Some residual diagnostics for the OLS regression {\tt lm(formula
= FemTheftConv \~{} Fert + Partic + Coll + MalTheftConv)}.}\label{olsdiag} 
\end{figure}

It can be shown that the OLS $\hat{\mb \beta}$'s are still unbiased
and consistent estimates of $\mb\beta$, but they will no longer be
BLUE (the proper GLS estimates are BLUE), and their standard error
estimates will be biased in general; the correct standard errors of
the OLS estimates when $\mb\epsilon\sim N(\mb0,\vc)$ can be easily
shown to be
\[
        \mbox{Var}\,(\hat{\mb\beta}) ~=~
        (X\tp X)^{-1} X\tp \vc X (X\tp X)^{-1}
\]

Using the estimated $\hat\rho= 0.2442$, we can construct the matrix
$\Gamma$ in \eq{rootR} directly, form an estimate of the
variance-covariance matrix $\vc$, and apply this formula to get
corrected SE's for the OLS $\hat\beta$'s.
\begin{verbatim}
> n _ length(FemTheftConv)
> rho _ 0.2442
> Gam _ diag(c(sqrt(1-rho^2),rep(1,length=n-1))) - rho*rbind(0,cbind(diag(n-1),0))
> SE2 _ sum(residuals(olsfit)^2)/olsfit$df.residual
> VC _ SE2*solve(t(Gam)%*%Gam)
> X _ model.matrix(olsfit)
> XTXI _ solve(t(X)%*%X)
> VarBeta _ XTXI %*% t(X) %*% VC %*% X %*% XTXI
> sqrt(diag(VarBeta))

[1] 11.745943500  0.001859433  0.028712245  0.053070901  0.023239818

> # now compare with the OLS standard errors...

> rbind(summary(olsfit)$coefficients[,2],sqrt(diag(VarBeta)))

     (Intercept)        Fert     Partic       Coll MalTheftConv 
[1,]     9.43792 0.001449503 0.02340973 0.04325091   0.01855587
[2,]    11.74594 0.001859433 0.02871224 0.05307090   0.02323982
\end{verbatim}

Why don't we stop here, since we already seem to have a solution?
\begin{itemize}
\item The OLS estimates $\hat{\mb\beta}$ are not particularly
efficient; and less-so as $\rho$ grows;
\item The estimate $\sigma^2_\epsilon$ is also not particularly
efficient; and hence the corrected OLS standard error estimates are
not terribly efficient;
\item We have not incorporated uncertainty in estimating $\rho$ into
our SE estimates, and so the ``real'' SE estimates are likely to be
somewhat higher yet.
\end{itemize}

\subsection{The Plug-In GLS analysis}

Using the matrix $\Gamma$ in \eq{rootR} that we constructed above to
adjust the OLS SE's for $\hat\beta$'s, we can form the adjusted
regression $\mb y^* = X^*\mb\beta +\mb \epsilon$ and apply OLS to this
adjusted problem.
\begin{verbatim}
> ystar _ Gam %*% FemTheftConv
> Xstar _ Gam %*% X
> summary(glsfit1 _ lm(ystar ~ Xstar - 1))

Call: lm(formula = ystar ~ Xstar - 1)
Residuals:
    Min     1Q  Median    3Q   Max 
 -9.294 -1.738 -0.2283 2.157 6.915

Coefficients:
                     Value Std. Error  t value Pr(>|t|) 
 Xstar(Intercept)  -6.6417  11.1218    -0.5972   0.5550
        XstarFert  -0.0059   0.0018    -3.3042   0.0025
      XstarPartic   0.1156   0.0271     4.2672   0.0002
        XstarColl   0.5359   0.0501    10.7058   0.0000
XstarMalTheftConv   0.0399   0.0220     1.8168   0.0796

[some output omitted]

> # compare corrected-SE OLS estimates with plug-in GLS estimates

> cbind(cbind(summary(olsfit)$coefficients[,1],sqrt(diag(VarBeta))),
+ summary(glsfit1)$coefficients[,1:2])
                                                  Value   Std. Error 
 Xstar(Intercept) -7.33414345 11.745943500 -6.641719523 11.121808677
        XstarFert -0.00609037  0.001859433 -0.005879305  0.001779358
      XstarPartic  0.11994156  0.028712245  0.115624068  0.027095855
        XstarColl  0.55155755  0.053070901  0.535866449  0.050053675
XstarMalTheftConv  0.03932477  0.023239818  0.039933442  0.021980373
\end{verbatim}
The results are quite close to the OLS estimates with corrected SE's.
The corrected OLS standard errors are quite similar to the plug-in GLS
standard errors; and because $\rho$ is not large in this problem, OLS
$\hat\beta$'s are quite close to the GLS $\hat\beta$'s as well.

\begin{quote}
{\sl
Note that this is {\em not} expected to be a great analysis, since we
are not incorporating the uncertainty in estimating $\rho$ into the
estimates of the other parameters; in particular I would expect the
standard errors of the regression coefficients to be underestimated in
this approach.
}
\end{quote}

\subsection{The GLS MLE analysis}

A maximum likelihood analysis depends on maximizing the likelihood
\begin{eqnarray*}
         \log L(\ul y; \ul\beta,\vc) 
        & = &
         \log\left\{\frac{1}{(2\pi)^{n/2}|\vc|^{1/2}}
                \exp \left[ - \frac12 (\ul y - X\ul \beta)\tp
                \vc^{-1} (\ul y - X\ul\beta) \right] \right\} \\
        & = & 
        -\frac{n}2 \log 2\pi - \frac{n}2 \log \sigma^2_z - \frac12
        \log (1-\rho^2) - \frac{1}{2\sigma^2_z} (\mb y^* - X^*\mb\beta)\tp
        (\mb y^* - X^*\mb\beta)
\end{eqnarray*}
(where $\mb y^* = \Gamma \mb y$ and $X^* = \Gamma X$ and $\Gamma$ is
from equation \ref{rootR}) simultaneously in 
\[
\beta_0, ~~~~ \beta_1, ~~~~ \ldots, ~~~~
\beta_k, ~~~~ \sigma^2_z ~~~~ \mbox{and} ~~~~ \rho
\]

For a problem with relatively few parameters like our analysis of
female theft convictions in Canada, there is an Splus function {\tt
nlmin()} that can be used to maximize the likelihood above directly in
all these parameters.  The function {\tt nlmin()} minimizes nonlinear
functions in whatever parameters are appropriate.  In this case, if we
minimize something proportional to $-\log L(\ul y; \ul\beta,\vc)$ we
will have maximized the likelihood. 

First we need to set up some helper functions, and pre-compute some
quantities that are repeatedly used in an {\tt nlmin()} analysis.  The
following computations and function definition are in the file {\tt
gls-code.splus}.  Note that the code is fairly specific to our theft
conviction data.

\clearpage

\begin{quote}\small
\begin{verbatim}
# Set up for "ar1fun"

# attach(hart)  # (don't forget!)

olsfit _ lm(FemTheftConv ~ Fert + Partic + Coll + MalTheftConv)
X _ model.matrix(olsfit)
y _ FemTheftConv
n _ length(y)

# initial values for "nlmin"

rho _ acf(residuals(olsfit),plot=F)$acf[2,1,1]
sig2z _ sum(residuals(olsfit)^2)/n * (1-rho^2)
b.start _ c(coefficients(olsfit),rho,sig2z)

# "ar1fun" is proportional to -(log-likelihood) so by 
# minimizing "ar1fun" we'll maximize the log-likelihood for 
# regression with ar1 errors.

ar1fun _ function(b) {
#
# Here are the params:
# b[1] - beta0
# b[2] - beta1
# b[3] - beta2
# b[4] - beta3
# b[5] - beta4
# b[6] - rho
# b[7] - sig^2_z
#
# this assumes that "hart" is attached, and that
#
# olsfit _ lm(FemTheftConv ~ Fert + Partic + Coll + MalTheftConv))
# X _ model.matrix(olsfit)
# y _ FemTheftConv
# n _ length(y)
# 
# have already been computed.
#
# WARNING: The following code is likely very inefficient.  But I bet it
# works!
#
        rho _ b[6]
        s2z _ b[7] #
#
        Gam _ diag(c(sqrt(1-rho^2),rep(1,length=n-1))) -
                rho*rbind(0,cbind(diag(n-1),0)) #
#
#       note that minimizing Q below maximizes the log-likelihood 
#       for linear regression with AR(1) errors
#
        Q _ n/2 * log(s2z) + 1/2 * log(1-rho^2) + 
                        1/(2*s2z) * sum((Gam%*%(y - X%*%b[1:5]))^2)

        return(Q)
}

#
# now try, e.g.
#
# nlmin(ar1fun,b.start,print.level=1)
# and help(nlmin)
#
#
\end{verbatim}
\end{quote}
After a little fooling around, the following seemed to work well
\begin{verbatim}
> result _ nlmin(ar1fun,b.start,print.level=1,max.fcal=100,max.iter=100)
    it    nf      f       reldf      preldf     reldx
     0     1  5.863e1
     1     7  5.863e1    1.970e-5   3.944e-5   4.587e-7
     2    10  5.863e1    3.560e-7   2.549e-6   7.485e-8
     3    15  5.861e1    2.649e-4   4.922e-4   2.014e-5

[many iterations omitted]

    30    60  5.822e1    2.414e-8   2.225e-8   6.693e-4
    31    61  5.822e1    4.016e-10  3.786e-10  1.048e-4
    32    62  5.822e1    2.945e-12  2.864e-12  5.495e-6
 ***** relative function convergence *****
 function     5.821717e1    reldx        5.4952e-6
 function evaluations:  62   gradient evaluations: 303
 preldf       2.8637e-12  npreldf      2.8637e-12
 
> result
$x:
[1]  -5.079305725  -0.005753718   0.110158212   0.520073564   0.039469484
[6]   0.398738667  11.354708284

$converged:
[1] T

$conv.type:
[1] "relative function convergence"
\end{verbatim}

We can directly calculate standard errors using the formula
$\mbox{Var}\,(\hat{\mb\beta}_{GLS}) = (X\tp \vc^{-1}X)^{-1}$ from
Equation \eq{GLS-stuff} and compare with our earlier plug-in
estimates.
\begin{verbatim}
> betaGLS _ result$x[1:5]
> rhoGLS _ result$x[6]
> s2zGLS _ result$x[7]
> s2eGLS _ s2zGLS/(1-rhoGLS^2)
> s2eGLS
[1] 13.50132
> rho _ rhoGLS
> Gam _ diag(c(sqrt(1-rho^2),rep(1,length=n-1))) -
+       rho*rbind(0,cbind(diag(n-1),0)) 
> VCI _ t(Gam)%*%Gam/s2eGLS
> VarBetaGLS _ solve(t(X)%*%VCI%*%X)
> ests _ cbind(cbind(summary(olsfit)$coefficients[,1],sqrt(diag(VarBeta))),
+ summary(glsfit1)$coefficients[,1:2],betaGLS,sqrt(diag(VarBetaGLS)))
> dimnames(ests) _ list(c("Intercept","Fert","Partic","Coll","MalTheftConv"),
+ c("OLS","SE(corrected)","P-I GLS","P-I SE","MLE GLS","MLE SE"))
> round(ests,4)
                 OLS SE(corrected) P-I GLS  P-I SE MLE GLS  MLE SE 
   Intercept -7.3341       11.7459 -6.6417 11.1218 -5.0793 12.8435
        Fert -0.0061        0.0019 -0.0059  0.0018 -0.0058  0.0022
      Partic  0.1199        0.0287  0.1156  0.0271  0.1102  0.0304
        Coll  0.5516        0.0531  0.5359  0.0501  0.5201  0.0569
MalTheftConv  0.0393        0.0232  0.0399  0.0220  0.0395  0.0253
\end{verbatim}
In every case, as expected, the SE's from the maximum-likelihood
version of GLS are larger than either the ``corrected'' OLS SE's or
the plug-in GLS SE's. 

How do we assess the success of a GLS fit?  The easiest residuals to
look at are
\[
        \mb E = \mb y^* - \mb X^*\hat{\mb\beta}_{GLS} 
        = \Gamma (\mb y - X\hat{\mb\beta}_{GLS})
\]
since they are supposed to be estimating uncorrelated, equal-variance
$N(0,\sigma^2_\epsilon)$.  Note that these are the ``usual'' residuals
that we would ``expect'', but tilted by $\Gamma$ so as to be
estimating uncorrelated equal-variance quantities.
\begin{verbatim}
> yhat _ X%*%betaGLS[1:5]
> resids _ Gam %*% (y - yhat)
> par(mfrow=c(2,2))
> plot(yhat,resids)
> abline(h=0)
> plot(resids)
> abline(h=0)
> plot(yhat,y)
> abline(0,1)
> qqenv(resids)
\end{verbatim}
Some potentially useful plots are in Figure~\ref{gls-diag}.  Note that
the $y$ vs $\hat y$ plot in the lower left displays essentially $\mb y
- X\hat{\mb\beta}_{GLS}$ instead of $\Gamma(\mb y - X\hat{\mb\beta}_{GLS})$.
Quantities like standardized residuals and Cook's D's could also be
calculated, and might well be useful, but these plots illustrate the
general ideas.

\begin{figure}
\centerline{\psfig{figure=gls-diag.ps,height=3.5in}}
\caption{Some uncorrelated raw-residual diagnostic plots for the MLE
GLS fit {\tt FemTheftConv \~{} Fert + Partic + Coll +
MalTheftConv}.}\label{gls-diag}
\end{figure}


\section{Some Closing Remarks}

GLS is a useful general technique for solving linear regression
problems when the correlations structure of the regression errors is
not too complicated.  Maximum likelihood GLS is a reasonable approach
to fitting linear regression models with time series structured
errors.  

For problems with many more parameters, {\tt nlmin()} may need to be
replaced with a method you code yourself.  Two common methods are
\begin{itemize}
\item Set up a grid of proposed correlations, say, $\rho_1, \rho_2,
\ldots, \rho_M$, calculate $\hat{\mb\beta}_m$ for each $\rho_m$ and
choosing the $\hat{\mb\beta}_m$; or
\item Iterate back and forth between estimating $\mb\beta$ from the
current guess for $\rho$, and re-estimating $\rho$ by applying
formula~\eq{corr-est} to the new set of residuals obtained.
\end{itemize}
(Neither of the above are the best possible approach but they do
illustrate valid approaches to the joint-MLE method).  The second
method above is closely related to a general technique for maximizing
likelihoods (and posterior distributions) over different ``types'' of
parameters (like regression coefficients, vs. error correlations),
called the E-M (``expectation-maximization'') algorithm.

GLS and related methods also turn out to be useful in nonlinear
optimization problems, since many nonlinear problems can be converted
into linear ones by application of Newton's method and/or Taylor's
Theorem.   One version of this, which we may see later, is called {\em
iteratively reweighted least squares} (IRLS).

\end{document}






%  The section at
% the end introducing iteratively reweighted least squares and robust
% regression is like parts of Fox Section 14.3.



\section{An extension: Iteratively Reweighted Least Squares (IRLS)}

{\sf 

(or: application to robust regression: IRLS)

}



\section{Closing remarks}

GLS is a useful general technique for solving linear regression
problems when the correlations structure of the regression errors is
not too complicated.  Maximum likelihood GLS is a reasonable approach
to fitting linear regression models with time series structured
errors.  

However, GLS is much more generally useful.  If we need to maximize or
minimize a function
\[
        y = F(\mb\beta)
\]
over its parameters $\mb\beta$, and if $F$ is nicely differentiable,
it is usually equivalent to solve the equation
\begin{equation}        
        \mb 0 = \frac{\partial F}{\partial \mb\beta}
        \label{zero}
\end{equation}
(for example, if $F$ is a log-likelihood, then these are the
likelihood equations).  Note that $\partial F / \partial\mb\beta$ is a
vector function.  Newton's method generalizes to vector functions for
example and leads to solving iteratively using the recursion
\[
        \mb\beta_{m+1} = \mb\beta_m - 
        \left.\left[\frac{\partial^2 F}{\partial\mb\beta\partial\mb\beta}
        \right]^{-1} 
        \frac{\partial F}{\partial\mb\beta}\right|_{\mb\beta=\mb\beta_m}
\]
Rather than inverting the Hessian matrix $\partial^2 F /
\partial\mb\beta\partial\mb\beta$ directly, it is usually better to
solve
\[
        \left.\frac{\partial F}{\partial\mb\beta}\right|_{\mb\beta=\mb\beta_m}
        =
        \left.\frac{\partial^2 F}{\partial\mb\beta\partial\mb\beta}
        \right|_{\mb\beta=\mb\beta_m}\mb\beta^* 
\]
and then compute
\[
        \mb\beta_{m+1} = \mb\beta_m + \mb\beta^*
\]

