# The Method of Maximum Likelihood
If we make the assumption that the error terms are normally
distributed, the maximum likelihood estimators (MLE) coincide with
the least square estimators. But MLE can also be applied
to a wide variety of models other than regression models, and it
generally generates estimators with excellent asymptotic
properties. The major disadvantage of MLE is that it requires
stronger distributional assumptions than does the method of
moments.
## Basic Concepts
Maximum likelihood estimation is to find the set of parameters which makes the current sample most likely given the statistical model.
Suppose we have a sample of $n$ independent and identically distributed (IID) observations. Each of them comes from some distribution function $f(\cdot, {\symbf{\theta)}}$, where parameters ${\symbf{\theta}}$ is fixed. The joint distribution of all the data points would be
$$ f(x_1, x_2, \dots, x_n | \theta) = f(x_1 | \theta) \times f(x_2 | \theta) \times \dots \times f(x_n | \theta). $$ {#eq-mle-1}
It is referred to as the likelihood function of the model for the
given data set.
Then it's natural to estimate $\theta$ by picking the $\theta$ that maximize this joint distribution function. A parameter vector $\hat {\symbf{\theta}}$ at which the
likelihood takes on its maximum value is called a maximum likelihood
estimate, or MLE, of the parameters.
Usually the log form of the likelihood function is preferred for computational purposes.
## MLE for a linear model
Consider the classical normal linear model:
$$
{\symbf{y=X\beta+u, \quad u \sim N(0, \sigma^2 I)
}}$$ {#eq-mle-2}
In this case, we have two parameters to estimate, $\sigma$ and ${\symbf{\beta}}$. The distribution of ${\symbf{y | X}}$ is
$$ f({{\symbf{y}}} \mid {{\symbf{X}}}, \beta, \sigma) = (2 \pi \sigma^2)^{-n/2} \exp\left(- \frac{({{\symbf{y- X \beta}}})'({{\symbf{y- X \beta}}})}{2 \sigma^2}\right)$$ {#eq-mle-3}
The log-likelihood function is
$$ l({{\symbf{y, \beta, \sigma}}}) = -\frac{n}{2} \log(2 \pi) - \frac{n}{2} \log(\sigma^2) - \frac{1}{2 \sigma^2} {\symbf{(y-X \beta)'(y-X \beta) }}$$ {#eq-mle-4}
Maximizing this function gives exactly the OLS estimator for ${\symbf{\beta}}$. That is algebra, and it holds whatever the true error distribution happens to be: the Gaussian log-likelihood is a decreasing function of the sum of squared residuals, so whatever minimizes one maximizes the other. For $\sigma^2$ the MLE divides by $n$ rather than $n-k$, so it differs in finite samples and agrees asymptotically. What normality buys is not the formula but its standing. Under normality the OLS estimator *is* the maximum likelihood estimator and the model-based variance is valid. Without normality the same formula is a quasi-MLE: still consistent for ${\symbf{\beta}}$ under the usual conditions, but its variance needs the sandwich correction discussed below.
## Asymptotic Properties of MLE
1. Consistency
$$
{{\mathrm{plim}}} ({\symbf{\hat \theta)=\theta.
}}$$ {#eq-mle-5}
This is to say that MLE estimator can get arbitrarily close to the true parameter if the sample size gets to infinity.
2. Asymptotic Normality
$$
{{\symbf{\sqrt{n}(\hat \theta - \theta)}}} \xrightarrow{d} N({\symbf{0, I^{-1}(\theta)),
}}$$ {#eq-mle-6}
where ${\symbf{I(\theta)}}$ below is the *per-observation* information.
Equivalently, in finite-sample terms,
$$
{{\symbf{\hat \theta}}} \approx N({\symbf{\theta, I_n^{-1}(\theta)), \qquad I_n(\theta) = n\, I(\theta),
}}$$ {#eq-mle-7}
with ${\symbf{I_n(\theta)}}$ the *full-sample* information. The $\sqrt{n}$
scaling (equivalently, dividing the full-sample information by $n$) is
what makes the statement well-defined; writing
$\hat\theta \sim N(\theta, I^{-1}(\theta))$ without specifying which
information is meant is ambiguous.
This states that the asymptotic distribution of ${\symbf{\hat \theta}}$
is normal with mean ${\symbf{\theta}}$ and variance given by the inverse
of the information, where the per-observation information matrix
${\symbf{I(\theta)}}$ is defined using the single-observation log-likelihood
contribution $l_t$ (so that the full-sample log-likelihood is
$l=\sum_{t=1}^n l_t$):
$$
{{\symbf{I(\theta)}}}=E[(\frac{\partial l_t}{\partial
\theta})(\frac{\partial l_t}{\partial
\theta})']=-E[(\frac{\partial^2 l_t}{\partial \theta \partial
\theta'})]
$$ {#eq-mle-8}
The information matrix is the expected value of the outer product of the gradient of a single observation's log-likelihood contribution, or the negative value of the expected value of its Hessian (second derivative). Using the *full-sample* $l$ here instead of $l_t$ would give the full-sample information $I_n(\theta) = n\,I(\theta)$, not $I(\theta)$ itself.
3. Asymptotic efficiency.
If ${\symbf{\hat \theta}}$ is the MLE estimator of $\theta$ and ${\symbf{V}}$ denotes its
*asymptotic* variance matrix -- the variance of the limiting distribution of
$\sqrt n (\hat \theta - \theta)$, not of ${\symbf{\hat \theta}}$ itself, the two
differing by a factor of $n$ exactly as with the information matrix above --
then
$$
\sqrt n (\hat \theta -\theta) \rightarrow^d N({\symbf{0, V), \qquad V = I^{-1}(\theta).
}}$$ {#eq-mle-9}
If ${\symbf{\tilde V}}$ denotes the asymptotic variance matrix of any other *regular*
consistent, asymptotically normal estimator, then ${\symbf{\tilde V- V}}$ is a positive
semidefinite matrix. That is, the MLE is the best estimator, in terms of
asymptotic variance, in that class. The regularity qualifier is not pedantry:
without it the statement is false, since superefficient estimators such as
Hodges' can beat the bound at a single point of the parameter space, at the cost
of much worse behaviour in a neighbourhood of it.
4. Invariance.
If $\hat \theta$ is the MLE of $\theta$ and $g(\theta)$ is a
continuous function of $\theta$, then $g(\hat \theta)$ is the MLE of
$g(\theta)$.
## When the Model Is Wrong: Quasi-MLE and M-Estimation
### What MLE Actually Minimizes
The asymptotic properties above assume that the true density belongs
to the parametric family $f(y \mid \theta)$. In practice this is
rarely literally true. What does MLE estimate when the model is
wrong?
Write $g(y)$ for the true density and $f(y \mid \theta)$ for the
parametric model. Maximizing the log-likelihood is equivalent to
minimizing the Kullback-Leibler divergence from the true density to
the model:
$$
\mathrm{KL}\bigl(g \,\|\, f(\cdot \mid \theta)\bigr)
= \int g(y) \log \frac{g(y)}{f(y \mid \theta)}\, dy
= E_g\!\left[\log g(Y)\right]
- E_g\!\left[\log f(Y \mid \theta)\right].
$$ {#eq-kl-div}
The first term does not depend on $\theta$, so minimizing KL over
$\theta$ is the same as maximizing $E_g[\log f(Y \mid \theta)]$,
which is the population version of the sample log-likelihood. The
MLE converges to the value $\theta^*$ that makes $f(\cdot \mid
\theta)$ closest to $g$ in the KL sense, whether or not $g$ belongs
to the family. This $\theta^*$ is called the pseudo-true value.
### The Information Matrix Inequality
Under correct specification, the two forms of the information matrix
are equal:
$$
J(\theta) = E\!\left[
\frac{\partial l_t}{\partial \theta}
\frac{\partial l_t}{\partial \theta}'
\right]
= -E\!\left[
\frac{\partial^2 l_t}{\partial \theta \,\partial \theta'}
\right]
= H(\theta).
$$ {#eq-info-equality}
This is the information matrix equality. Under misspecification it
fails: $J(\theta^*) \neq H(\theta^*)$ in general. The asymptotic
variance of the MLE is no longer $I^{-1}(\theta)$ but the sandwich
$$
V = H(\theta^*)^{-1}\, J(\theta^*)\, H(\theta^*)^{-1}.
$$ {#eq-sandwich-mle}
When the model is correctly specified, $J = H$ and the sandwich
collapses to $H^{-1} = I^{-1}$, recovering the standard result.
When it is not, $H^{-1}$ is simply the wrong variance -- inconsistent, rather than merely
too small. There is no general ordering between $H^{-1}$ and the sandwich: it can be too
small in some parameter directions and too large in others.
### Quasi-Maximum Likelihood Estimation
The sandwich variance makes MLE usable even under misspecification,
provided we are clear about what parameter is being estimated. A
quasi-maximum likelihood estimator (QMLE) is an MLE computed from a
likelihood that may be wrong, together with a sandwich variance that
accounts for the misspecification.
The leading case in econometrics is Poisson QMLE for count or
nonnegative data. The Poisson log-likelihood for observation $t$ is
$$
l_t = y_t \log \mu_t - \mu_t - \log(y_t!), \quad
\mu_t = \exp({\symbf{x}}_t'\beta).
$$ {#eq-poisson-qmle}
The score equation $\sum_t (y_t - \mu_t)\, {\symbf{x}}_t = 0$ depends
only on the conditional mean $E[Y_t \mid {\symbf{x}}_t] = \mu_t$,
not on whether $Y_t$ is actually Poisson. If the conditional mean
is correctly specified, the QMLE $\hat\beta$ is consistent for
$\beta$ regardless of the true distribution of $Y_t$ --- the data
need not even be counts. The variance, however, must be estimated
with the sandwich, since the Poisson variance assumption
$\mathrm{Var}(Y_t \mid {\symbf{x}}_t) = \mu_t$ will generally be wrong.
A fuller treatment is in the [Count Data chapter](count-data.qmd).
The same logic applies to logit and probit for binary data: as long
as the conditional mean $P(Y=1 \mid X) = F(X'\beta)$ is correctly
specified, the MLE is consistent for $\beta$ even if the error
distribution is not exactly logistic or normal. In this sense,
binary response MLE is already quasi-MLE.
### M-Estimators
The QMLE perspective generalizes further. An M-estimator is any
estimator defined as the solution to
$$
\sum_{t=1}^n \psi(y_t, {\symbf{x}}_t, \theta) = 0,
$$ {#eq-m-est}
where $\psi$ is a known vector-valued function. MLE is the special
case where $\psi = \partial l_t / \partial \theta$, the score. OLS
is the case $\psi(y_t, {\symbf{x}}_t, \beta) = {\symbf{x}}_t(y_t -
{\symbf{x}}_t'\beta)$, the population orthogonality condition from the
[OLS chapter](ols.qmd). GMM generalizes further by allowing more
equations than parameters; the [GMM chapter](gmm.qmd) develops this.
All M-estimators share the same asymptotic structure under
regularity conditions:
$$
\sqrt{n}\,(\hat\theta - \theta^*)
\xrightarrow{d}
N\!\left(0,\;
A^{-1}\, B\, A^{-1\prime}
\right),
$$ {#eq-m-est-avar}
where $A = E[\partial \psi / \partial \theta']$ and
$B = E[\psi \psi']$. This is the sandwich again. When $\psi$ is a
score from a correctly specified likelihood, $A = -H$, $B = J = H$,
and the sandwich reduces to $H^{-1}$.
The M-estimator framework unifies OLS, MLE, QMLE, and GMM as
special cases of the same asymptotic theory. The sandwich variance
is the general formula; the simpler $I^{-1}$ or $\sigma^2(X'X)^{-1}$
are special cases that hold under correct specification.
## Example
Let's look at the example of regression of car's mpg on disp, hp and wt. This time we use MLE. (R has a mle() function which I cannot make it work on this example. So I am using optim().)
```{r}
ols.lf <- function(theta, y, X) {
beta <- theta[-1]
sigma2 <- theta[1]
# Return a large penalty rather than NA: L-BFGS-B's line search can
# transiently evaluate points that violate the box constraint, and
# `fn` returning NA there causes a fatal "L-BFGS-B needs finite values
# of 'fn'" error. A finite penalty keeps the optimizer on track instead.
if (sigma2 <= 0) return(1e10)
n <- nrow(X)
e <- y - X%*%beta
logl <- ((-n/2)*log(2*pi)) - ((n/2)*log(sigma2)) - ((t(e)%*%e)/(2*sigma2))
return(-logl)
}
X <- cbind(1,mtcars[,c('disp','hp','wt')])
y <- mtcars$mpg
# Deliberately naive starting values (sigma2 = 1, slopes at +/-1): the point
# of the exercise is that the optimizer finds the OLS solution without being
# handed it. For harder likelihoods than this one, good starting values matter
# a great deal -- they reduce the risk of converging to a spurious local
# maximum, or of failing to converge at all.
optim(c(1,1,-1,-1,-1), method="L-BFGS-B", fn=ols.lf, lower=c(1e-6,-Inf,-Inf,-Inf,-Inf), upper=rep(Inf,5), y=y, X=as.matrix(X))
# it should match the OLS estimator
lm1 <- lm(mpg~disp+hp+wt, data=mtcars)
summary(lm1)
```
## Hypotheses Testing
### Linear Hypotheses Testing
Linear hypothesis can be framed as
$$
{\symbf{H_0: \ \ R\beta -r=0,
}}$$ {#eq-mle-10}
where ${\symbf{R}}$ is a $q \times k$ matrix of known constants with $q \le k$ and
$\mathrm{rank}({\symbf{R}}) = q$, and ${\symbf{r}}$ is a $q$-vector of known constants.
The restrictions may number as many as the parameters -- testing that every coefficient
is zero is the case $q = k$ -- but they must be linearly independent.
We assume that $u_i$ are normally distributed.
$$
{\symbf{u \sim N(0, \sigma^2 I)
}}$$ {#eq-mle-11}
Even if this assumption does not hold, we still have asymptotic normality for ${\symbf{\hat \beta}}$, so the tests below remain valid asymptotically.
Since linear combination of normal variables are also normally
distributed,
$$
\begin{aligned}
{\symbf{\hat \beta }}&{\symbf{\sim N(\beta, \sigma^2 (X'X)^{-1})}}\\{\symbf{
R \hat \beta }}&{\symbf{\sim N(R \beta, \sigma^2 R(X'X)^{-1}R')}}\\{\symbf{
R (\hat \beta - \beta) }}&{\symbf{\sim N(0, \sigma^2 R(X'X)^{-1}R')
}}\end{aligned}
$$ {#eq-mle-12}
Under null hypothesis,
$$
{\symbf{(R \hat \beta - r) \sim N(0, \sigma^2 R(X'X)^{-1}R')
}}$$ {#eq-mle-13}
Therefore,
$$
{{\symbf{(R \hat \beta - r)'[\sigma^2 R(X'X)^{-1}R']^{-1} (R \hat
\beta - r)}}} \sim \chi^2(q)
$$ {#eq-mle-14}
Here we don't know $\sigma^2$. However, it can be shown that
$$
{{\symbf{\frac{\hat u' \hat u}{\sigma^2}}}} \sim \chi^2(n-k)
$$ {#eq-mle-15}
Then we have
$$
{{\symbf{\frac{(R \hat \beta - r)'[ R(X'X)^{-1}R']^{-1} (R \hat \beta
- r)/q}{\hat u' \hat u / (n-k)}}}} \sim F(q, n-k)
$$ {#eq-mle-16}
This is an example of the $F$ test for any linear hypotheses testing. $t$ test will be a special case of this.
### The three classical tests: LR, Wald and LM
The $F$ test above leans on normal errors. The general machinery for testing
${\symbf{H_0: R\theta - r = 0}}$ after maximum likelihood needs only the asymptotic
results already established, and it comes in three forms. Write $\hat\theta_u$ for
the unrestricted MLE, $\hat\theta_r$ for the MLE subject to the restriction,
$s(\theta) = \partial \ell / \partial \theta$ for the score, and
${\symbf{I_n(\theta) = n\,I(\theta)}}$ for the full-sample information.
**Likelihood ratio.** Compare the two maximized log-likelihoods:
$$
LR = 2\bigl[\ell(\hat\theta_u) - \ell(\hat\theta_r)\bigr] \xrightarrow{d} \chi^2(q).
$$ {#eq-mle-17}
**Wald.** Ask whether the *unrestricted* estimate is far from satisfying the
restriction, measuring distance in units of its own estimated variance:
$$
W = ({{\symbf{R}}}\hat\theta_u - {{\symbf{r}}})' \bigl[{{\symbf{R\, I}}}_n^{-1}(\hat\theta_u)\,{{\symbf{R}}}'\bigr]^{-1} ({{\symbf{R}}}\hat\theta_u - {{\symbf{r}}}) \xrightarrow{d} \chi^2(q).
$$ {#eq-mle-18}
**Lagrange multiplier (score).** Ask whether the score, evaluated at the
*restricted* estimate, is far from zero -- if the restriction is harmless, the
restricted point should already be nearly a maximum:
$$
LM = s(\hat\theta_r)'\, {{\symbf{I}}}_n^{-1}(\hat\theta_r)\, s(\hat\theta_r) \xrightarrow{d} \chi^2(q).
$$ {#eq-mle-19}
The three are asymptotically equivalent under the null and have the same local
power, so the choice among them is one of convenience: **Wald needs only the
unrestricted fit, LM needs only the restricted fit, and LR needs both.** That is
why LM is the natural choice for specification testing, where the restricted model
is the easy one to estimate (White's heteroskedasticity test and the
Breusch-Godfrey serial-correlation test are both LM tests run as auxiliary
regressions), while Wald is what software reports by default, since the
unrestricted fit is what was just computed.
In the classical normal linear model the three are monotone transformations of the
same $F$ statistic above, so nothing is lost by using $F$ there. Two cautions
otherwise. In finite samples the three can disagree, sometimes enough to straddle
a conventional critical value, and there is a known ordering
$W \ge LR \ge LM$ in the linear model under normality. And the Wald statistic is
**not invariant to nonlinear reparameterization**: testing $\gamma = 1$ and
testing $\log\gamma = 0$ are the same hypothesis but give different Wald
statistics, whereas $LR$ is invariant because it compares likelihood *values*,
which reparameterization does not change. When a hypothesis can be written several
algebraically equivalent ways, prefer $LR$.
## Test of Structural Change (Chow test)
Chow test is a test between two groups of observations. For example,
we have data on consumption function of prewar period and postwar
period. It is natural to think of a test on whether consumption
function parameters differ between prewar and postwar periods. The
null hypothesis is that there is no difference.
Let ${\symbf{y_i}}$, ${\symbf{X_i}}$ ($i=1,2)$ indicate the appropriate
partitioning of the data. The unrestricted model may be written
$${\symbf{}}\begin{bmatrix}{\symbf{y_1 }}\\{\symbf{ y_2
}}\end{bmatrix}=
\begin{bmatrix}{\symbf{X_1 }}&{\symbf{ 0 }}\\{\symbf{ 0 }}&{\symbf{ X_2
}}\end{bmatrix}
\begin{bmatrix}{\symbf{\beta_1 }}\\{\symbf{ \beta_2
}}\end{bmatrix}
+u \quad u \sim N(0, \sigma^2I)
$$ {#eq-mle-20}
The null hypothesis of no structural break is
$$
{\symbf{H_0: \ \ \beta_1 = \beta_2.
}}$$ {#eq-mle-21}
Running OLS on two equations separately, we have
$${\symbf{}}\begin{bmatrix}{\symbf{\hat \beta_1 }}\\{\symbf{ \hat \beta_2
}}\end{bmatrix}=
\begin{bmatrix}{\symbf{X_1' X_1 }}&{\symbf{ 0 }}\\{\symbf{ 0 }}&{\symbf{ X_2' X_2
}}\end{bmatrix}^{-1}
\begin{bmatrix}{\symbf{X_1' y_1 }}\\{\symbf{ X_2' y_2
}}\end{bmatrix}
=\begin{bmatrix}{\symbf{(X_1' X_1)^{-1}X_1'y_1 }}\\{\symbf{ (X_2' X_2)^{-1}
X_2' y_2
}}\end{bmatrix}
$$ {#eq-mle-22}
Restricted model under null hypothesis
$${\symbf{}}\begin{bmatrix}{\symbf{y_1 }}\\{\symbf{ y_2
}}\end{bmatrix}=
\begin{bmatrix}{\symbf{X_1 }}\\{\symbf{ X_2
}}\end{bmatrix}
{\symbf{\beta +u
}}$$ {#eq-mle-23}
The test of the null hypothesis is given by
$$
{{\symbf{F= \frac{(\hat u_*' \hat u_*- \hat u' \hat u)/k}{\hat u' \hat
u/ (n-2k)}}}} \sim F(k , n-2k)
$$ {#eq-mle-24}
where $\hat u_*$ is the residual of the restricted model, $\hat u$ is the stacked residual of two unrestricted models.
For example, for prewar period (suppose there is $n_1$ observations):
$$
y_1 = \beta_0 + \beta_1 x_1 + \beta_2 x_2
$$ {#eq-mle-25}
For postwar period (suppose there is $n_2$ observations):
$$
y_2 = \gamma_0 + \gamma_1 x_1 + \gamma_2 x_2
$$ {#eq-mle-26}
The null hypothesis of no structural break is
$$
{\symbf{H_0: \ \ \beta_0 = \gamma_0, \beta_1 = \gamma_1, \beta_2 = \gamma_2.
}}$$ {#eq-mle-27}
Running OLS on two equations separately, we have two sums of squares
of error ($SSR_1$ and $SSR_2$). The unrestricted error sum of squares
is
$$
SSR_U = SSR_1 + SSR_2
$$ {#eq-mle-28}
Then run the regression on the stacked sample of $n_1+n_2$ observations
and get $SSR_R$.
Then
$$
\frac{(SSR_R-SSR_U)/k}{SSR_U/(n_1+n_2-2k)} \sim F_{k, n_1+n_2-2k}.
$$ {#eq-mle-29}
$k$ is number of restrictions; in this example, it is 3.
The other way to do the same test is easier: simply include
interaction terms between group dummy and explanatory variables in the
regression.
For our earlier example, we would estimate the model:
$$
y=\alpha_0 + \alpha_1 x_1 + \alpha_2 x_2 + \eta_0 d + \eta_1 x_1 d + \eta_2
x_2 d
$$ {#eq-mle-30}
Here $d$ is a dummy variable for postwar period. The Chow test is
simply a test on the hypothesis that all the coefficients involving
$d$ are zero. A regression of the above model with $n+m$ observations
and a $F$ test are easy to do.
$$
{\symbf{H_0: \ \ \eta_0 = \eta_1 = \eta_2 = 0.
}}$$ {#eq-mle-31}
### An example of Chow test in R
Here we use an example from the "strucchange" library.
```{r}
## Example 7.4 from Greene (1993), "Econometric Analysis"
## Chow test on Longley data
data("longley")
library(strucchange)
## use structural change test from the library.
sctest(Employed ~ Year + GNP.deflator + GNP + Armed.Forces, data = longley,
type = "Chow", point = 7)
## which is equivalent to segmenting the regression via
fac <- factor(c(rep(1, 7), rep(2, 9)))
## here fac is the factor for segmenting it at point 7.
fm0 <- lm(Employed ~ Year + GNP.deflator + GNP + Armed.Forces, data = longley)
fm1 <- lm(Employed ~ fac/(Year + GNP.deflator + GNP + Armed.Forces), data = longley)
## Here anova is returning the F test (Chow test). Equivalent to sctest results.
anova(fm0, fm1)
## estimates from Table 7.5 in Greene (1993)
summary(fm0)
summary(fm1)
```
## Heteroskedasticity and Autocorrelation Consistent Standard Errors
### Variance estimator under homoscedasticity
In OLS regression, the variance-covariance matrix of $\hat \beta$ is $${{\mathrm{var}}} {\symbf{(\hat \beta) = {E}[(\hat \beta-\beta)(\hat \beta -\beta)']
= {E}[(X'X)^{-1}X'uu'X(X'X)^{-1}] = \sigma^2 (X'X)^{-1} }}$$ {#eq-mle-32}
if $${{\mathrm{var}}} {\symbf{(u)=\sigma^2 I}}$$ {#eq-mle-33}.
That is, the error term has mean zero and constant variance
(homoscedasticity). The pairwise correlation between error terms
is always zero (no serial correlation).
The estimation of $\sigma^2$:
$$
s^2=\frac{ {\symbf{\hat u' \hat u}}}{n-k}
$$ {#eq-mle-34}
is an unbiased estimator of $\sigma^2$ (proof omitted).
The standard estimate of the variance-covariance matrix of the OLS
parameter estimates under the assumption of IID errors is
$$
{{\mathrm{\hat {Var}}}}(\hat \beta)=s^2 {\symbf{(X'X)^{-1}
}}$$ {#eq-mle-35}
### White’s estimator
We made strong assumption that the error terms of the regression
model are IID when we estimate the variance-covariance matrix of
OLS estimator. Under this assumption, the usual estimator of
variance-covariance matrix of ${\symbf{\hat \beta}}$ is consistent. Now
let's relax this assumption to only independent but not
identically distributed. Again, the linear regression models is
$$
{\symbf{y=X\beta+u, \quad {E} (u)=0, \quad {E} (u u')=\Omega,
}}$$ {#eq-mle-36}
where ${\symbf{\Omega}}$ is the error variance-covariance matrix with
diagonal elements being $\sigma_t^2$ for $t^{th}$ element,
off-diagonal elements being zero. In other words, the error terms
are heteroscedastic.
The variance-covariance matrix of the OLS estimator ${\symbf{\hat
\beta}}$ is equal to
$$
\begin{aligned}
{{\mathrm{Var}}} {\symbf{(\hat \beta) }}&{\symbf{= {E}[(\hat \beta
-\beta)(\hat
\beta -\beta)'] }}\\{\symbf{
}}&{\symbf{= [(X'X)^{-1}X'({E}(uu'))X(X'X)^{-1}] }}\\{\symbf{
}}&{\symbf{= (X'X)^{-1} X' \Omega X (X'X)^{-1}
}}\end{aligned}
$$ {#eq-mle-37}
If we know $\sigma_t^2$, then we would be able to estimate this
"sandwich covariance matrix". But we don't.
$$
{{\mathrm{Var}}} {\symbf{(\hat \beta) = \frac{1}{n}
[\frac{1}{n}(X'X)]^{-1}[\frac{1}{n}X'\Omega
X][\frac{1}{n}(X'X)]^{-1}
}}$$ {#eq-mle-38}
Let $y_t$ denote the $t$th observation on the dependent variable,
and $x_t'=[1 x_{2t} \cdots x_{kt}]$ denote the $t$th row of the
${\symbf{X}}$ matrix. Then
$$
{\symbf{X' \Omega X=\sum_{t=1}^n \sigma_t^2 x_t x_t'
}}$$ {#eq-mle-39}
The White estimator replaces the unknown $\sigma_t^2$ by $\hat u_t^2$,
the estimated OLS residuals. This provides a consistent estimator of
the variance matrix for the OLS coefficient vector and is particularly
useful since it does not require any specific assumptions about the
form of the heteroscedasticity.
Therefore,
$$
\begin{aligned}
\hat {{\mathrm{Var}}} {\symbf{(\hat \beta) }}&{\symbf{= \frac{1}{n}
[\frac{1}{n}(X'X)]^{-1}[\frac{1}{n}X' \hat \Omega
X][\frac{1}{n}(X'X)]^{-1}}}\\{\symbf{
\hat \Omega }}&{\symbf{= {diag} ({\hat u_1^2, \hat u_2^2, \cdots,
\hat u_n^2})
}}\end{aligned}
$$ {#eq-mle-40}
### Newey-West estimator
White's estimator deals with the situation that we have
heteroskedasticity (a diagonal $\Omega$) of unknown form. When we
have serial correlation of unknown form (a non-diagonal $\Omega$), we
can estimate the variance-covariance matrix by a heteroskedasticity
and autocorrelation consistent, or HAC, estimator. Newey-West
estimator is the most popular HAC estimator.
Given a time series data set, suppose we are interested in estimating
the mean vector (suppose we have more than one variable) and its
variance. We know that given IID data, we can apply central limit
theorem: sample mean is a consistent estimator of the population mean
and it's variance can be calculated since asymptotically the sample
mean conforms to a normal distribution and the variance can be
estimated, relatively easily. However, in the case of time series
data, autocorrelation usually exists. We may be concerned the CLT may
not work in this case.
Fortunately, autocorrelation does not by itself defeat inference, though the
conditions need care. Suppose ${\symbf{y_t}}$ is covariance-stationary (the
covariance is not a function of time) *and* its autocovariances are absolutely
summable, $\sum_{v=-\infty}^{\infty} \lVert \Gamma_v \rVert < \infty$. Then, as in
Hamilton (1994), the sample mean satisfies:
$${\symbf{\bar y_t \to \mu , }}$$ {#eq-mle-41}
$$ {{\symbf{S}}} = \lim_{T \to \infty} {\symbf{ {T \cdot E[(\bar y_T - \mu)(\bar y_T -\mu)' ]} = \sum_{v=-\infty}^{\infty} \Gamma_v . }}$$ {#eq-mle-42}
where ${\symbf{\Gamma_v}}$ is the variance-covariance matrix for ${\symbf{y_t}}$ and ${\symbf{y_{t-v}}}$.
The first one says the law of large numbers still holds. The second gives the
long-run variance used for the standard error, and absolute summability is what
makes that sum exist.
Covariance stationarity on its own is not enough for either statement. It can
fail even the law of large numbers. Let $y_t = z + \varepsilon_t$, where $z$ is
drawn once and then held fixed and $\varepsilon_t$ is white noise. Every
autocovariance equals $\mathrm{Var}(z)$, so the process is covariance-stationary,
the autocovariances are not summable, and $\bar y_T$ converges to $z$ rather than
to $\mu$. Asymptotic normality asks for more again: a mixing, near-epoch
dependence, or martingale-difference condition, together with moment
restrictions. Covariance stationarity buys the variance formula, not the
central limit theorem.
If the data were generated by a vector MA(q) process, then
$$ {{\symbf{S}}} = \sum_{v=-q}^{q} \Gamma_v . $$ {#eq-mle-43}
A natural estimate is
$$ {\symbf{\hat S = \hat \Gamma_0 + \sum_{v=1}^{q} (\hat \Gamma_v + \hat \Gamma_v' ),}}$$ {#eq-mle-44}
where $$ {\symbf{\hat \Gamma_v = (1/T) \sum_{t=v+1}^{T} (y_t- \bar y)(y_{t-v} - \bar y)'. }}$$ {#eq-mle-45}
This gives a consistent estimate of ${\symbf{S}}$; however, it sometimes is
not positive semidefinite.
Newey-West (1987) suggested putting in a weight:
$$ {{\symbf{\hat S = \hat \Gamma_0}}} + \sum_{v=1}^{q} {(1- \frac{v}{q+1})} ({\symbf{\hat \Gamma_v + \hat \Gamma_v' ), }}$$ {#eq-mle-46}
Here $q$ is a truncation bandwidth, not a true MA order. For an MA($q$) process the two
coincide, but the Newey-West estimator is used precisely when the dependence has no known
finite order: the Bartlett weights $1 - v/(q+1)$ downweight higher lags and keep
$\hat S$ positive semidefinite, and consistency comes from letting $q$ grow with $T$,
slowly enough that the poorly estimated high-order autocovariances do not dominate.
Consider a linear regression model:
$$ y_t={ {\symbf{x_t' \beta}}} + u_t $$ {#eq-mle-47}
Suppose we have the OLS estimator ${\symbf{b_T}}$, then
$$ \sqrt{T} ({{\symbf{b_T - \beta}}}) = [(1/T) {\sum_{t=1}^T {\symbf{x_t
x_t'}}}]^{-1} [(1/{\sqrt{T}}) {\sum_{t=1}^{T} {\symbf{x_t u_t}}}] $$ {#eq-mle-48}
The first term converges in probability to some constant. The second
term is $\sqrt{T}$ times the sample mean of the vector ${\symbf{x_t u_t}}$,
which is exactly the form to which a central limit theorem applies.
Under general conditions,
$$ \sqrt{T} ({{\symbf{b_T - \beta}}}) \rightarrow^L N(0, Q^{-1}SQ^{-1}) $$ {#eq-mle-49}
where $Q = {{\mathrm{plim}}} \, (1/T)\sum_{t=1}^T {{\symbf{x_t x_t'}}} = {{\mathrm{E}}}[{{\symbf{x_t x_t'}}}]$, and $S$ can be estimated by
$$ {{\symbf{\hat S_T = \hat \Gamma_{0T}}}} + \sum_{v=1}^{q} {(1- \frac{v}{q+1})} ({\symbf{\hat \Gamma_{v,T} + \hat \Gamma_{v,T}' ), }}$$ {#eq-mle-50}
where
$$ \hat \Gamma_{v,T}= (1/T) \sum_{t=v+1}^T (x_t \hat u_{t, T} \hat u_{t-v, T} x_{t-v}'), $$ {#eq-mle-51}
where $\hat u_{t,T}$ is the OLS residual for data $t$ in a sample of size $T$.
Overall, the variance of ${\symbf{b_T}}$ is approximated by
$$\hat \Sigma_{NW} = [\sum_{t=1}^T x_t x_t']^{-1} [ \sum_{t=1}^T \hat u_t^2 x_t x_t' + \sum_{v=1}^{q} {(1- \frac{v}{q+1})} \sum_{t=v+1}^T (x_t \hat u_{t, T} \hat u_{t-v, T} x_{t-v}' + x_{t-v} \hat u_{t-v, T} \hat u_{t, T} x_{t}') ] [\sum_{t=1}^T x_t x_t']^{-1} $$ {#eq-mle-52}
This estimation obviously depends on the selection of $q$, the lag
length beyond which the autocorrelation of $x_t u_t$ and $x_{t-v} u_{t-v}$ is
treated as negligible. One common rule of thumb sets
$q = 0.75 \cdot T^{\frac{1}{3}}$; it is a convention rather than a theorem, and
other growth rates are in use.
Newey-West (1994) has suggested a way to automatically select the
bandwidth $q$. Here we omitted the discussion. Both Stata and R now
also implement Newey-West (1994) estimator, with no need to specify
$q$.