Software 9.5.1 (Maximum likelihood by numerical optimisation). Maximise the log-likelihood of Theorem 9.2.1 directly, as a check on the theory. Writing \(\sigma ^2 = e^{\phi }\) keeps it positive during the search:

> X <- model.matrix(fit); y <- LifeCycleSavings$sr; n <- length(y)
> negll <- function(th) { b <- th[1:5]; s2 <- exp(th[6])
+   r <- y - X %*% b
+   0.5 * n * log(2 * pi * s2) + sum(r^2) / (2 * s2) }
> o <- optim(c(coef(fit), log(10)), negll, method = "BFGS",
+            control = list(maxit = 1000, reltol = 1e-14))
> rbind(ML = o$par[1:5], OLS = coef(fit))
    (Intercept)      pop15     pop75           dpi      ddpi
ML     28.56609 -0.4611932 -1.691498 -0.0003369020 0.4096951
OLS    28.56609 -0.4611931 -1.691498 -0.0003369019 0.4096949
> c(sigma2_ML = exp(o$par[6]), RSS_over_n = deviance(fit) / n,
+   s2 = deviance(fit) / 45)
 sigma2_ML RSS_over_n         s2
  13.01426   13.01426   14.46029
> c(-o$value, logLik(fit))
[1] -135.0981 -135.0981

Reading the output. The numerical maximum reproduces least squares to the seventh significant figure, the remaining differences being optimiser tolerance. The MLE of \(\sigma ^2\) is \(\text {RSS}/n = 13.01\), smaller than \(s^2 = \text {RSS}/(n - k) = 14.46\) by the factor \(45/50\). The maximised log-likelihood \(-135.10\) is what R’s logLik reports; it is the quantity that AIC and BIC are built from. Here \(\text {AIC} = -2\ell + 2p\) with \(p = 6\) parameters, counting \(\sigma ^2\).

Stata: ml model lf myprog (sr = pop15 pop75 dpi ddpi) /lnsigma, then ml maximize. EViews: a LogL object with @logl specification.

Software 9.5.2 (The Wald, likelihood ratio and LM tests). The joint test of Software 7.7.1, that pop15 and pop75 have no effect, in its three asymptotic forms:

> fit0 <- lm(sr ~ dpi + ddpi, data = LifeCycleSavings)
> RU <- deviance(fit); RR <- deviance(fit0)
> c(W = n * (RR - RU) / RU, LR = n * log(RR / RU),
+   LM = n * (RR - RU) / RR, crit = qchisq(0.95, 2))
        W        LR        LM      crit
13.370338 11.848645 10.549366  5.991465
> e0 <- resid(fit0); n * summary(lm(e0 ~ X - 1))$r.squared
[1] 10.54937
> 2 * (logLik(fit) - logLik(fit0))
'log Lik.' 11.84864 (df=6)

Reading the output. All three exceed the \(5\%\) point \(5.99\) of \(\chi ^2_2\), and they are ordered \(W > \text {LR} > \text {LM}\), as Theorem 9.3.1 guarantees. \(W = \frac {50}{45} \times 2 \times 6.0167 = 13.37\), linking it to the exact \(F\) test. The auxiliary regression of the restricted residuals on all five regressors reproduces LM as \(nR^2\); X - 1 uses the columns of \(X\), which already include the constant. And twice the difference of the two logLik values is LR. In a sample of \(50\) the differences between the three are noticeable, \(10.5\) to \(13.4\). The exact \(F\) test, available here because the model is linear and normal, is preferable to all three.

Stata: lrtest after storing the two models with estimates store; test gives the Wald form. EViews: View \(\to \) Coefficient Diagnostics offers the Wald and Omitted Variables – Likelihood Ratio tests.

Software 9.5.3 (GLS against least squares: a Monte Carlo). Generate data with \(\Var (u_i) = 0.09X_i^2\), so the disturbance standard deviation grows with \(X\), and compare least squares with GLS, which here is weighted least squares with weights \(1/X_i^2\). R’s weights argument takes the reciprocals of the \(\omega _i\):

> set.seed(9)
> n <- 50; x <- runif(n, 1, 10)            # fixed regressor
> sim <- replicate(5000, {
+   y <- 2 + 0.5 * x + rnorm(n, sd = 0.3 * x)
+   c(ols = coef(lm(y ~ x))[2],
+     gls = coef(lm(y ~ x, weights = 1 / x^2))[2]) })
> rbind(mean = rowMeans(sim), sd = apply(sim, 1, sd))
         ols.x      gls.x
mean 0.5007509 0.49918251
sd   0.1029070 0.06762586

and the exact standard deviations from Theorems 9.4.1 and 9.4.2 are \(0.1036\) for least squares and \(0.0680\) for GLS.

Reading the output. Both estimators are unbiased: the means are \(0.501\) and \(0.499\) against the true \(0.5\). But GLS has about two-thirds of the standard deviation of least squares, equivalent to having more than twice as many observations (\(0.103^2/0.068^2 = 2.3\)). The simulated standard deviations agree with the exact formulas to two significant figures. Least squares wastes information by giving the noisy observations at large \(X\) as much weight as the precise ones at small \(X\).

Stata: regress y x [aweight = 1/x^2]. EViews: in the equation dialog, Options \(\to \) Weights, with the weight series \(1/x\) (EViews weights by standard deviation).