Software 11.7.1 (Durbin–Watson statistic and its exact distribution). The statistic, the residual correlogram, and the exact null distribution of \(d\) for these regressors, simulated from part (ii) of Theorem 11.3.2 as \(MU\) with \(U\) standard normal:

> sb <- data.frame(Seatbelts); sb$t <- 1:192
> sb$month <- factor(cycle(Seatbelts[, "drivers"]))
> fit <- lm(log(drivers) ~ log(kms) + log(PetrolPrice) + law + month + t,
+           data = sb)
> e <- resid(fit); n <- length(e); X <- model.matrix(fit)
> d <- sum(diff(e)^2) / sum(e^2)
> c(d = d, r1 = sum(e[-1] * e[-n]) / sum(e^2))
        d        r1
1.3359231 0.3302132
> round(acf(e, lag.max = 6, plot = FALSE)$acf[2:7], 3)
[1] 0.330 0.295 0.147 0.026 0.129 0.003
> M <- diag(n) - X %*% solve(crossprod(X), t(X))
> A <- crossprod(diff(diag(n)))                  # d = e'Ae / e'e
> c(Ed = sum(diag(M %*% A)) / (n - ncol(X)))
      Ed
2.028209
> set.seed(1)
> dsim <- replicate(20000, {u <- M %*% rnorm(n); sum(diff(u)^2) / sum(u^2)})
> c(mean = mean(dsim), quantile(dsim, c(0.01, 0.05)), min = min(dsim),
+   as.small = sum(dsim <= d))
    mean       1%       5%      min as.small
2.027099 1.687188 1.784067 1.390448 0.000000

Reading the output. \(d = 1.34\), close to \(2(1 - \hat r) = 2(1 - 0.330) = 1.34\), as Theorem 11.3.2(i) says it should be with \(n = 192\). The residual autocorrelations decline slowly from \(0.33\) and \(0.30\). The exact mean of \(d\) under \(H_0\), \(\tr (MA)/(n - k) = 2.028\) from Theorem 11.3.2(iii), is confirmed by the simulation, which also gives the exact \(5\%\) critical value for these regressors, \(1.78\). Not one of the \(20\,000\) simulated values is as small as \(1.34\), and the smallest is \(1.39\): the evidence of positive autocorrelation is overwhelming. The matrix A is built as \(D'D\), with diff(diag(n)) the first-difference matrix \(D\).

Stata: after tsset, estat dwatson. EViews: the Durbin–Watson statistic is printed in every equation output; the correlogram is View \(\to \) Residual Diagnostics \(\to \) Correlogram – Q-statistics.

Software 11.7.2 (Breusch–Godfrey and Ljung–Box tests). The Breusch–Godfrey test at orders \(1\) and \(12\), the second because the data are monthly, and the Ljung–Box test with \(12\) lags:

> bg <- function(p) {
+   E <- embed(c(rep(0, p), e), p + 1)[, -1, drop = FALSE]
+   LM <- n * summary(lm(e ~ X[, -1] + E))$r.squared
+   c(LM = LM, df = p, p = pchisq(LM, p, lower.tail = FALSE))
+ }
> rbind(bg(1), bg(12))
           LM df            p
[1,] 21.45457  1 3.623109e-06
[2,] 39.18633 12 9.804716e-05
> Box.test(e, lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  e
X-squared = 53.781, df = 12, p-value = 2.99e-07

Reading the output. embed lines up \(e_t\) with its first \(p\) lags, after padding with \(p\) zeros as in Theorem 11.3.3. All three tests reject independence decisively: \(\text {LM} = 21.5\) on \(1\) degree of freedom, \(39.2\) on \(12\), and \(Q^* = 53.8\) on \(12\). The conventional standard errors and the Chow tests of Software 8.7.2 are therefore suspect.

Stata: estat bgodfrey, lags(1 12); wntestq on the saved residuals for the Ljung–Box \(Q\). EViews: View \(\to \) Residual Diagnostics \(\to \) Serial Correlation LM Test, entering the number of lags.

Software 11.7.3 (Feasible GLS and maximum likelihood). The Prais–Winsten transformation (11.2) is written as a function; dropping the first row gives Cochrane–Orcutt, which is then iterated. Maximum likelihood with AR(1) and AR(2) disturbances uses R’s arima, which accepts regressors through xreg:

> y <- log(sb$drivers); r <- sum(e[-1] * e[-n]) / sum(e^2)
> pw <- function(rho, first = TRUE) {      # first = FALSE: Cochrane-Orcutt
+   ys <- y[-1] - rho * y[-n]; Xs <- X[-1, ] - rho * X[-n, ]
+   if (first) {
+     ys <- c(sqrt(1 - rho^2) * y[1], ys)
+     Xs <- rbind(sqrt(1 - rho^2) * X[1, ], Xs)
+   }
+   lm(ys ~ Xs - 1)
+ }
> rho <- 0
> for (i in 1:50) {                         # iterated Cochrane-Orcutt
+   u <- y - X %*% coef(pw(rho, FALSE))
+   new <- sum(u[-1] * u[-n]) / sum(u[-n]^2)
+   if (abs(new - rho) < 1e-8) break
+   rho <- new
+ }
> c(r = r, rho.CO = rho, iterations = i)
         r     rho.CO iterations
 0.3302132  0.3719266 10.0000000
> ml1 <- arima(y, order = c(1, 0, 0), xreg = X[, -1])
> ml2 <- arima(y, order = c(2, 0, 0), xreg = X[, -1])
> round(c(AR1 = ml1$coef[1], AR2 = ml2$coef[1:2]), 3)
AR1.ar1 AR2.ar1 AR2.ar2
  0.370   0.294   0.221
> c(AIC1 = AIC(ml1), AIC2 = AIC(ml2),
+   Q1 = Box.test(resid(ml1), 12, "Ljung-Box", fitdf = 1)$p.value,
+   Q2 = Box.test(resid(ml2), 12, "Ljung-Box", fitdf = 2)$p.value)
         AIC1          AIC2            Q1            Q2
-440.27190797 -447.70936705    0.04433298    0.62001488

Reading the output. The residual estimate \(\hat r = 0.330\) rises to \(\hat \rho = 0.372\) when Cochrane–Orcutt is iterated to convergence, and maximum likelihood gives \(0.370\). The least squares residuals understate \(\rho \), because they are forced to be orthogonal to the regressors, which are themselves smooth. But the AR(1) model does not remove all the autocorrelation: the Ljung–Box test on its residuals, with \(12 - 1 = 11\) degrees of freedom, still gives \(p = 0.044\). The AR(2) model, with coefficients \(0.294\) and \(0.221\), leaves white-noise residuals (\(p = 0.62\)) and has the lower AIC. This is the pattern that the slowly declining correlogram of Software 11.7.1 suggested.

Stata: prais y x1 x2 (Prais–Winsten), prais y x1 x2, corc (Cochrane–Orcutt), and arima y x1 x2, ar(1/2) for maximum likelihood. EViews: add ar(1) ar(2) to the equation specification, ls y c x1 x2 ar(1) ar(2).

Software 11.7.4 (Newey–West standard errors and a comparison). Newey–West standard errors with the rule-of-thumb lag, and a table of the two coefficients of most interest across all the methods:

> nw <- function(L) {                      # Newey-West, Bartlett weights
+   g <- X * e; S <- crossprod(g)
+   for (j in 1:L) {
+     G <- crossprod(g[-(1:j), ], g[1:(n - j), ])
+     S <- S + (1 - j / (L + 1)) * (G + t(G))
+   }
+   B <- solve(crossprod(X)); sqrt(diag(B %*% S %*% B))
+ }
> floor(4 * (n / 100)^(2/9))               # rule-of-thumb lag
[1] 4
> se <- function(f) sqrt(diag(f$var.coef))
> v <- c("log(kms)", "law")
> tab <- rbind(
+   OLS     = c(coef(fit)[v], sqrt(diag(vcov(fit)))[v]),
+   NW4     = c(coef(fit)[v], nw(4)[v]),
+   PW      = c(summary(pw(r))$coef[c(2, 4), 1:2]),
+   CO      = c(summary(pw(rho, FALSE))$coef[c(2, 4), 1:2]),
+   ML.AR1  = c(ml1$coef[v], se(ml1)[v]),
+   ML.AR2  = c(ml2$coef[v], se(ml2)[v]))
> colnames(tab) <- c("b.kms", "b.law", "se.kms", "se.law")
> round(tab, 4)
        b.kms   b.law se.kms se.law
OLS    0.4321 -0.1504 0.1208 0.0213
NW4    0.4321 -0.1504 0.1180 0.0244
PW     0.2479 -0.1521 0.1282 0.0288
CO     0.2163 -0.1509 0.1294 0.0306
ML.AR1 0.2226 -0.1529 0.1318 0.0291
ML.AR2 0.2013 -0.1599 0.1272 0.0347

Reading the output. g holds the rows \(g_t' = e_tx_t'\), and G is \(\hat \Gamma _j\), exactly as in Theorem 11.4.3. The seat-belt law effect is remarkably stable: every method puts it between \(-0.150\) and \(-0.160\), a fall in deaths of \(14\) to \(15\%\). Its standard error, however, rises from \(0.021\) (least squares) to between \(0.024\) and \(0.035\) once the autocorrelation is allowed for. The \(t\) ratio falls from \(-7.1\) to about \(-4.6\) in the AR(2) model, still decisive. The kilometre elasticity is more fragile: it falls from \(0.43\) to between \(0.20\) and \(0.25\) under GLS, with standard errors of about \(0.13\), so it is no longer clearly different from zero. Least squares and GLS are both consistent if the model is right, so a shift of about one and a half standard errors between them is a hint that the dynamics of the relation between traffic and deaths are not fully captured by a static equation with AR disturbances, in the spirit of Theorem 11.4.2.

Stata: newey y x1 x2, lag(4). EViews: Options \(\to \) Coefficient covariance matrix \(\to \) HAC (Newey-West) in the estimation dialog.

Software 11.7.5 (A dynamic model: Durbin’s \(h\) fails). Add last month’s log deaths as a regressor, so that the model has a lagged dependent variable, and test its residuals:

> ylag <- y[-n]; dyn <- lm(y[-1] ~ ylag + X[-1, -1])
> ed <- resid(dyn); m <- length(ed)
> dd <- sum(diff(ed)^2) / sum(ed^2); V <- vcov(dyn)["ylag", "ylag"]
> c(b = coef(dyn)[["ylag"]], d = dd, r = 1 - dd/2, nV = m * V,
+   h = (1 - dd/2) * sqrt(m / (1 - m * V)))
          b           d           r          nV           h
 0.35315042  2.11649667 -0.05824834  0.95891278 -3.97142918
> Z <- model.matrix(dyn); e1 <- c(0, ed[-m])
> LM <- m * summary(lm(ed ~ Z[, -1] + e1))$r.squared
> c(LM = LM, p = pchisq(LM, 1, lower.tail = FALSE))
        LM          p
3.22045105 0.07272359

Reading the output. The lagged dependent variable has coefficient \(0.353\), and the Durbin–Watson statistic is now \(2.12\): by Theorem 11.2.2, \(d\) is biased towards \(2\) in a dynamic model, and it can no longer be used. Durbin’s \(h\) gives \(-3.97\), an apparently significant negative autocorrelation. But \(n\widehat \Var (b_1) = 0.959\) is close to one, so the factor \(\sqrt {n/(1 - n\widehat \Var (b_1))} = 68\) magnifies a tiny \(\hat r = -0.058\) into a large statistic, and \(h\) is not to be trusted. The Breusch–Godfrey test, valid in this setting, gives \(\text {LM} = 3.22\) with \(p = 0.073\): no significant first-order autocorrelation remains once the lag is included.

Stata: estat durbinalt (Durbin’s alternative, a form of the Breusch–Godfrey test) and estat bgodfrey. EViews: the Serial Correlation LM Test is valid with lagged dependent variables.

Software 11.7.6 (ARCH test and GARCH(1,1) estimation). R’s EuStockMarkets data give daily closing values of the DAX from 1991 to 1998. Percentage log returns are tested for ARCH and a GARCH(1,1) is fitted by maximising the likelihood with optim:

> r <- 100 * diff(log(EuStockMarkets[, "DAX"])); n <- length(r)
> u <- r - mean(r)
> c(n = n, sd = sd(r), kurtosis = mean(u^4) / mean(u^2)^2)
          n          sd    kurtosis
1859.000000    1.030084    9.279689
> Q <- function(x, ...) unname(Box.test(x, 10, "Ljung-Box", ...)$statistic)
> c(Q.r = Q(r), Q.r2 = Q(r^2))
       Q.r       Q.r2
  6.365577 110.746179
> E <- embed(u^2, 6)                      # u2[t], u2[t-1], ..., u2[t-5]
> LM <- (n - 5) * summary(lm(E[, 1] ~ E[, -1]))$r.squared
> c(LM = LM, p = pchisq(LM, 5, lower.tail = FALSE))
          LM            p
6.971090e+01 1.177043e-13
> hfun <- function(p) {                   # GARCH(1,1) conditional variances
+   u <- r - p[1]; h <- numeric(n); h[1] <- var(r)
+   for (t in 2:n) h[t] <- p[2] + p[3] * u[t - 1]^2 + p[4] * h[t - 1]
+   h
+ }
> nll <- function(p) {
+   if (p[2] <= 0 || p[3] < 0 || p[4] < 0 || p[3] + p[4] >= 1) return(1e10)
+   h <- hfun(p); 0.5 * sum(log(2 * pi * h) + (r - p[1])^2 / h)
+ }
> g <- optim(c(mu = 0, omega = 0.05, alpha = 0.1, beta = 0.8), nll,
+            hessian = TRUE)
> round(cbind(estimate = g$par, se = sqrt(diag(solve(g$hessian)))), 4)
      estimate     se
mu      0.0653 0.0216
omega   0.0476 0.0127
alpha   0.0684 0.0148
beta    0.8876 0.0236
> c(loglik = -g$value, persistence = sum(g$par[3:4]),
+   uncond.var = unname(g$par[2] / (1 - sum(g$par[3:4]))))
       loglik   persistence    uncond.var
-2594.7966981     0.9559975     1.0814973
> z <- (r - g$par[1]) / sqrt(hfun(g$par))
> c(Q.z2 = Q(z^2, fitdf = 2), kurt.z = mean(z^4), min.z = min(z),
+   kurt.z.drop.min = mean(z[-which.min(z)]^4))
           Q.z2          kurt.z           min.z kurt.z.drop.min
       0.893289       15.944858      -12.335787        3.490465

Reading the output. The kurtosis of the returns is \(9.3\), far above the normal \(3\). The returns are serially uncorrelated (Ljung–Box \(Q = 6.4\) on \(10\) degrees of freedom), but their squares are strongly autocorrelated (\(Q = 110.7\)), and Engle’s test with five lags gives \(\text {LM} = 69.7\), \(p \approx 10^{-13}\), exactly the pattern of Theorem 11.6.2. The fitted model is \[ h_t = 0.048 + 0.068\,u_{t-1}^2 + 0.888\,h_{t-1}, \] with persistence \(\hat \alpha + \hat \beta = 0.956\) and so a half-life of \(\log 0.5/\log 0.956 = 15\) trading days. The implied unconditional variance, \(1.08\), matches the sample variance \(1.03^2 = 1.06\). The standardised residuals \(z_t = \hat u_t/\sqrt {h_t}\) have no remaining ARCH (\(Q = 0.89\)). Their kurtosis is \(15.9\) only because of a single day, the August 1991 fall shown in Figure 11, which came in a calm period and is \(12.3\) conditional standard deviations; without it the kurtosis is \(3.5\). Such outliers are why practitioners often replace the normal by a Student \(t\) distribution for \(\varepsilon _t\).

Stata: estat archlm, lags(5) after regress; arch r, arch(1) garch(1). EViews: View \(\to \) Residual Diagnostics \(\to \) Heteroskedasticity Tests \(\to \) ARCH; estimate with Quick \(\to \) Estimate Equation, method ARCH.