\documentstyle[12pt,fullpage,psfig,times]{article}

\parskip=0.5em
\parindent=0pt

\frenchspacing

\long\def\comment#1{}

\def\tp{^T}

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

\def\vc#1{{\underline{#1}}}

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

\def\Varsigma{\mb{\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}}

\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 
        Odds and Ends from the Theory of Linear Regression
\end{center}



\section{MLE Properties of Least-Squares}

\subsection{The Multivariate Normal Density}

The univariate normal density ($n=1$) we all remember is 
\[
        \phi(y) = \frac{1}{\sqrt{2\pi}\sigma} \exp \left[-\frac12 
        \left( \frac{y-\mu}{\sigma}\right)^2 \right]
\]
for $y \in \Re$.  We often write $Y \sim N(\mu,\sigma^2)$ to convey
that the random variable $Y$ has this density.

The multivariate normal density, for $\ul y \in \Re^n$, is
\begin{equation}
        \phi(\ul y) = \frac{1}{(2\pi)^{n/2}|\vc|^{1/2}}
        \exp [ -\frac12 (\ul y - \ul \mu)\tp \vc^{-1} (\ul y - \ul \mu) ]
\label{mvn}
\end{equation}
The mean vector $\ul \mu$ replaces the univariate mean $\mu$ and the
variance-covariance matrix $\vc$ replaces the univariate variance
$\sigma^2$. We often write $\ul Y \sim N_n(\ul\mu,\vc)$, or just $\ul
Y \sim N(\ul\mu,\vc)$ if the dimension $n$ is clear from context.
Sometimes we also say that the components of $\ul Y$ are ``jointly
normal'' or ``multivariate normal''.

One way to build up the multivariate normal is to start with a
$n$-dimensional vector of independent identically distributed $N(0,1)$
variables $\ul Z$ $=$ $(Z_1, \ldots, Z_n)$.  Since densities for
independent random variables multiply, the density for $\ul Z$ must be
\begin{eqnarray*}
        \phi(\ul z) & = & \prod _1^n \phi(z_i)  \\
        & = & \frac{1}{(2\pi)^{n/2}} 
        \exp [ - \frac12 \sum_1^n z_i^2 ] \\
        & = & \frac{1}{(2\pi)^{n/2}} 
        \exp [ - \frac12 \ul z\tp I^{-1} \ul z ]
\end{eqnarray*}
where $I$ is the $n\by n$ identity matrix.  

Then if $\ul Y = A\ul Z + \ul \mu$, it is clear that $E[\ul Y]$ $=$
$\ul \mu$ and $\Var(\ul Y)$ $=$ $E[(\ul Y - \ul \mu)(\ul Y - \ul
\mu)\tp]$ $=$ $AA\tp$ $=$ $\vc$ (so for example it is immediately
obvious that $\vc$ must by symmtric and positive definite).  

The multivariate change-of-variables formula from calculus will ensure
that the density for $\ul Y$ looks like \eq{mvn}.  In fact, it is
generally true that if $\ul Y = A\ul Y + \ul b$, where $A$ is $m\by n$
and $\ul b \in \Re^m$, and $\ul Y \sim N_n(\ul \mu, \vc)$, then $\ul Y
\sim N_m(A\ul \mu + \ul b, A\vc A\tp)$.

\comment{

NOTHING WRONG WITH THIS, JUST SAVING SPACE

\clearpage

As a more explicit example, in the case $k=2$ we often write
\[
        \ul \mu = \left[ \begin{array}{c}\mu_1 \\ \mu_2 \end{array} \right] ;
        \vc = \left[ \begin{array}{cc} \sigma_1^2 & \rho\sigma_1\sigma_2 \\
                                    \rho\sigma_1\sigma_2 & \sigma_2^2
                \end{array} \right]
\]
in which case you should be able to figure out from \eq{mvn} that
\small
\begin{eqnarray*}
        \lefteqn{\phi(y_1, y_2)} \\
         & = & \frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}
        \exp \left\{ -\frac1{2(1-\rho^2)} \left[ 
        \left(\frac{y_1 - \mu_1}{\sigma_1}\right)^2
        - 2\rho \left(\frac{y_1 - \mu_1}{\sigma_1}\right)\left(\frac{y_2 -
                \mu_2}{\sigma_2}\right)
        + \left(\frac{y_2 - \mu_2}{\sigma_2}\right) \right] \right\} \\
        & = & 
        \frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}
        \exp \left\{ -\frac1{2(1-\rho^2)} \left[ z_1^2 - 2z_1z_2 + z_2^2 
                \right] \right\}
\end{eqnarray*}
\normalsize
where $z_i = (y_i - \mu_i)/\sigma_i$; what is the matrix $A$ in the
transformation from $\ul Z$ to $\ul Y$ mentioned above?

}



\subsection{Maximum Likelihood Estimation}

Now let us derive maximum likelihood estimates for $\ul\beta$ and
$\sigma^2_\epsilon$ in the multiple regression model
\[
        \ul y_{n\by 1} = X_{n\by k'} \ul \beta_{k'\by 1} + \ul
        \epsilon_{n\by 1}    ~~~~ \ul\epsilon \sim
        N_n(\ul 0, \sigma^2_\epsilon I_{n\by n}) .
\]
Since this implies the $y_i$'s are independent
\[
        \ul y \sim N_n(X\beta, \sigma^2_\epsilon I_{n\by n})
\]
we can immediately write down the likelihood for $\ul\beta$ and
$\sigma^2_\epsilon$ by writing down the product of the individual
normals\footnote{{\em Warning:} In this handout, $\ul x_i\tp$ is the
  $i^{th}$ \ul{\em row} of the regression matrix $X$.},
\begin{eqnarray*}
        L(\ul y; \ul\beta,\sigma^2_\epsilon) 
        & = & \prod _1^n \frac{1}{\sqrt{2\pi\sigma^2_\epsilon}}
                \exp \left[ - \frac12  (y_i-\ul x_i\tp\ul\beta)^2 / 
                \sigma^2_\epsilon \right] \\
        & = & \frac{1}{(2\pi\sigma^2_\epsilon)^{n/2}} 
                \exp \left[ - \frac12 \sum_1^n (y_i-\ul x_i\tp\ul\beta)^2 /
                \sigma^2_\epsilon \right] \\ 
        & = & \frac{1}{(2\pi)^{n/2}|\sigma^2 I|^{1/2}} 
                \exp \left[ - \frac12 (\ul y - X\ul \beta)\tp 
                (\sigma^2_\epsilon I)^{-1} (\ul y - X\ul\beta) \right] ,
\end{eqnarray*}
and hence the log-likelihood is
\[
        \ell(\ul\beta,\sigma^2_\epsilon) 
        = \log L((\ul y; \ul\beta,\sigma^2_\epsilon)
        = - \frac n2 \log 2\pi - \frac n2 \log \sigma^2_\epsilon 
                - \frac{1}{2\sigma^2_\epsilon} (\ul y - X\ul \beta)\tp
                (\ul y - X\ul\beta)
\]
Now we differentiate with respect to $\ul\beta$ and
$\sigma^2_\epsilon$ to obtain
\begin{eqnarray*}
        \frac{\partial\ell}{\partial\ul\beta} & = &
        - \frac{1}{2\sigma^2_\epsilon}[-2X\tp(y - X\ul\beta)] \\
        \frac{\partial\ell}{\partial\sigma^2_\epsilon} & = &
        - \frac n2 \frac{1}{\sigma^2_\epsilon} 
        + \frac{1}{\sigma^4_\epsilon} (\ul y - X\ul \beta)\tp
                (\ul y - X\ul\beta)
\end{eqnarray*}
Setting the first set of partial derivatives equal to zero produces the
normal equations again, and hence
\[
        \hat{\ul\beta} = (X\tp X)^{-1} X\tp \ul y ,
\]
the usual least-squares regression estimates.  Setting the second set
of partial derivatives equal to zero produces
\[
        \hat\sigma^2_\epsilon = \frac{ (\ul y - X\hat{\ul\beta})\tp
                (\ul y - X\hat{\ul\beta})}{n} = \frac{RSS}{n}
\]
which we now know is a biased estimate of $\sigma^2_\epsilon$ (the
correct denominator should be $n-k'$).


\subsection{Some Geometry; Eigensystems; Square Roots; Principal Components}

A useful geometric property of the multivariate normal distribution is
that it is ``log-quadratic''; $\log \phi(\ul y)$ is a simple linear
function of the symmetric positive-definite quadratic form $(\ul y -
\ul \mu)\tp \vc^{-1} (\ul y - \ul \mu)$ in $(\ul y - \ul \mu)$.  This
means that the level surfaces $\{\ul y: \phi(\ul y) = K\}$, or
equivalently the level surfaces $\{\ul y: \log \phi(\ul y)= c\}$ will
be ellipsoids.

Finding the principal axes of the ellipsoidal level sets (contours of
$\phi(\ul y)$) is easy with Lagrange multipliers.  To find the first
principal axis for example, you need a point $\ul y$ that is a maximum
distance from $\ul \mu$ [maximize the squared distance $(\ul y - \ul
\mu)\tp(\ul y - \ul \mu)$], but is still on the level set [subject to
the constraint $(\ul y - \ul \mu)\tp\vc^{-1}(\ul y - \ul \mu) = c$].
Differentiating the Lagrange-multiplier objective function
\begin{equation}        
        f(\ul y) = (\ul y - \ul \mu)\tp(\ul y - \ul \mu) - 
        \lambda [(\ul y - \ul \mu)\tp\vc^{-1}(\ul y - \ul \mu) - c]
\label{pc-lagrange}
\end{equation}
with respect to $\ul y$ and setting these derivatives equal to zero
leads to the eigenvalue/eigenvector problem
\[
        (\vc - \lambda I) (\ul y - \ul \mu) = 0
\]
So there must be multiple critical points for \eq{pc-lagrange}, one
for one corresponding to each eigenvalue of $\sigma$.  It now a matter
of calculation to see that the eigenvalues will satisfy
\[
        \lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p > 0
\]
with corresponding eigenvectors $\ul e_i$ $=$ $(\ul y_i - \ul \mu)$
lying along the principal axes of the ellipsoid $(\ul y - \ul
\mu)\tp\vc^{-1}(\ul y - \ul \mu) = c$; the half-lengths of the principal
axes will be $\sqrt{\lambda c}$.

Several interesting algebraic and geometric properties come out of
this analysis.  To every eigenvalue $\lambda_i$ there corresponds an
eigenvector $\ul e_i$, and it can be shown that all the $\ul e_i$ are
linearly independent.  If two eigenvectors $\ul e_i$ and $\ul e_j$
correspond to different eigenvalues $\lambda_i \neq \lambda_j$ then,
since $\ul e_j\tp \vc \ul e_i$ is a scalar (hence equal to its transpose),
\[
        \lambda_i \ul e_j\tp\ul e_i = \ul e_j\tp \vc \ul e_i = [\ul e_j\tp
        \vc \ul e_i]\tp = \ul e_i\tp \vc \ul e_j = \lambda_j \ul e_i\tp\ul e_j
\]
so it must be that $\ul e_i\tp\ul e_i = 0$: the eigenvectors from
different eigenvalues are orthogonal.  If $\lambda_i$ $=$ $\lambda_j$,
there will still be two eigenvectors, and using the Gram-Schmidt
procedure we can select two eigenvectors that are orthogonal to each
other and to all the other eigenvectors in the problem.  If the two
eigenvalues are equal, then the cross-section of the ellipse in the
plane containing the two linearly independent eigenvectors $\ul e_i$ and
$\ul y_j$ will be circular (it would always be an ellipse, and now the
major and minor principal axes will be of the same length
$2\sqrt{\lambda_i c}$); if three are equal the corresponding
cross-section will be spherical, etc.

We can assume all the $e_i$ are orthogonal; so let $\ul a_i$ = $\ul
e_i / \| \ul e_i \|$ and consider the matrix
\[
        A = \left[\begin{array}{cccc}
                \ul a_1 & \ul a_2 & \cdots & \ul a_p
        \end{array}\right]
\]
You should verify that $AA\tp = I$ (the matrix is orthogonal---it has
orthonormal columns); that $\vc A = A \, \diag(\lambda_1, \ldots,
\lambda_p)$ (the square $p\by p$ matrix with $\lambda$'s on the main
diagonal and zeros everywhere else); and therefore that $A\tp\vc A =
\diag(\lambda_1, \ldots, \lambda_p)$.  

Now consider the random vector $\ul W = A (\ul Y - \ul \mu)$.  This vector
$\ul W$ will be multivariate normal, and you should verify that it has
mean $\ul 0$ (the zero vector) and variance-covariance matrix
$\diag(\lambda_1, \ldots, \lambda_p)$.  Hence each $W_i$ is
uncorrelated with $W_j$ (in fact, independent); and $\Var(W_i) =
\lambda_i$.  The $W_i$'s are called the {\em principal components\/}
of $\ul Y$ and they play a large role in understanding the structure
of $\ul Y$.


\section{The Gauss-Markov Theorem}

In the multiple regression model
\[
        \ul y_{n\by 1} = X_{n\by k'} \ul \beta_{k'\by 1} + \ul
        \epsilon_{n\by 1} 
\]
we will only assume that the error terms $\epsilon_i$ are independent,
mean 0, constant variance $\sigma^2_\epsilon$.  

Let us define $M = (X\tp X)^{-1}X\tp$, so that $\hat{\ul\beta} = M\ul y$
are the least-squares estimators of $\ul\beta$, and let
$\tilde{\ul\beta} = U\ul y$ be any unbiased linear estimate of
$\beta$.   Finally, let $A=U-M$.  

\begin{lemma}
        For any unbiased linear estimate $U\ul y$ of $\ul\beta$,
        $A=U-M$ is orthogonal to the rows of $X$; hence $AM\tp = MA\tp = 0$.
\end{lemma}

{\em Proof:}  Since both $\tilde{\ul\beta}$ and $\hat{\ul\beta}$ are
unbiased, 
\[
        \ul\beta = E[\tilde{\ul\beta}] = E[(M+A)\ul y] 
        = E[M\ul y] + E[M\ul y] = E[\hat{\ul\beta}] + AE[\ul y]
        = \ul\beta + AX\ul\beta ,
\]
and since this must be true for all $\ul\beta$ it follows that
$AX=0$. It immediately follows that $AM\tp = AX(X\tp X)^{-1} = 0$. $\Box$

\begin{theorem} {\em Gauss-Markov Theorem.}

        $\hat{\ul\beta}$ is the Best (minimum-variance) Linear Unbiased
        Estimate (BLUE) of $\ul\beta$.
\end{theorem}

{\em Proof:} We will do this by minimizing the variance of the
arbitrary unbiased linear estimate $\tilde{\ul\beta} = U\ul y =
(M+A)\ul y$.  Since $AM\tp = MA\tp = 0$, it follows that
\begin{eqnarray*}
        \Var(\tilde{\ul\beta}) & = & \Var[(M+A)\ul y] \\
        & = & (M+A) \Var(\ul y) (M+A)\tp ~ = ~ (M+A) \sigma^2_\epsilon
        I (M+A)\tp \\
        & = & \sigma^2_\epsilon (MM\tp + MA\tp + AM\tp + AA\tp) \\
        & = & \sigma^2_\epsilon (MM\tp + AA\tp)
\end{eqnarray*}
and hence, reading off the diagonal elements of this
variance-covariance matrix,
\[
        \Var(\tilde\beta_j) = \sigma^2_\epsilon \left(
        \sum_{i=1}^n m_{ij}^2 + \sum_{i=1}^n a_{ij}^2 
        \right) .
\]
Since the $m_{ij}$'s are fixed, we choose different unbiased linear
estimates $\tilde{\ul\beta}$ by choosing different $a_{ij}$'s; from
the above expression it is clear that the only way to minimize the
variance of all the $\tilde\beta_j$'s is to take $a_{ij} \equiv 0$;
so that $\tilde\beta_j=\hat\beta_j$ for every $j$. $\Box$

{\em Remark:} Essentially the same technique can be used to prove the
following result:
\begin{corollary}
        Let $\phi=L\ul\beta$ be a linear function of the $\beta$'s.
        The BLUE of $\phi$ is $L\hat{\ul\beta}$ (i.e. plug in the
        least-squares estimates).
\end{corollary}




\end{document}





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





