Software 10.5.1 (Least squares with conventional and robust standard errors). Base R computes the sandwich in a few lines; the sandwich package does the same with vcovHC(fit, type = "HC3"):

> st <- data.frame(state.x77)
> st$TI  <- st$Income * st$Population / 1000     # total income, $ million
> st$Pop <- st$Population / 1000                  # population, millions
> fit <- lm(TI ~ Pop, data = st)
> summary(fit)$coef
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -1431.807  466.14137 -3.071615 3.501713e-03
Pop          4904.810   76.05455 64.490683 2.727591e-48
> X <- model.matrix(fit); e <- resid(fit); n <- nrow(X); k <- ncol(X)
> XtXi <- solve(crossprod(X))
> h <- rowSums((X %*% XtXi) * X)                  # leverages
> hc <- function(w) sqrt(diag(XtXi %*% crossprod(X * w) %*% XtXi))
> round(rbind(conventional = sqrt(diag(vcov(fit))),
+             HC0 = hc(e), HC1 = hc(e) * sqrt(n / (n - k)),
+             HC3 = hc(e / (1 - h))), 2)
             (Intercept)    Pop
conventional      466.14  76.05
HC0               446.63 129.69
HC1               455.84 132.36
HC3               566.20 172.61

Reading the output. Each additional million people adds about $4 905 million of income, about $4 900 a head, close to average per-capita income, as one would expect. crossprod(X * w) computes \(\sum w_i^2x_ix_i'\), the meat of the sandwich, with \(w_i = e_i\) for HC0 and \(e_i/(1 - p_{ii})\) for HC3. The robust standard error of the slope is \(130\) to \(173\), against a conventional \(76\): the conventional formula understates the uncertainty by a factor of two or more. That is the direction Theorem 10.1.1 predicts, since the large residuals occur at the populous states, where \(x_i^2\) is largest. The robust \(t\) ratio is still enormous, \(4905/173 = 28\), so the conclusion survives; with a weaker relationship it might not.

Stata: regress TI Pop, vce(robust) gives HC1, and vce(hc3) gives HC3. EViews: in the estimation dialog, Options \(\to \) Coefficient covariance \(\to \) Huber-White.

Software 10.5.2 (Breusch–Pagan, White and Goldfeld–Quandt tests). The three tests, each in a few lines:

> e2 <- e^2; aux <- lm(e2 ~ Pop, data = st)
> c(koenker = n * summary(aux)$r.squared,
+   p = pchisq(n * summary(aux)$r.squared, 1, lower.tail = FALSE))
     koenker            p
1.734921e+01 3.110270e-05
> g <- e2 / mean(e2); a2 <- lm(g ~ Pop, data = st)
> c(BP = sum((fitted(a2) - mean(g))^2) / 2)
      BP
25.92774
> w <- lm(e2 ~ Pop + I(Pop^2), data = st)
> c(white = n * summary(w)$r.squared,
+   p = pchisq(n * summary(w)$r.squared, 2, lower.tail = FALSE))
       white            p
1.735017e+01 1.707886e-04
> o <- order(st$Pop); lo <- st[o[1:20], ]; hi <- st[o[31:50], ]
> r1 <- deviance(lm(TI ~ Pop, lo)); r2 <- deviance(lm(TI ~ Pop, hi))
> c(GQ = r2 / r1, crit = qf(0.95, 18, 18),
+   p = pf(r2 / r1, 18, 18, lower.tail = FALSE))
          GQ         crit            p
2.170859e+01 2.217197e+00 1.095935e-08

Reading the output. Every test rejects homoscedasticity decisively. Koenker’s \(nR^2 = 17.3\) on \(1\) degree of freedom. The original Breusch–Pagan statistic, \(\frac 12\text {ESS}\) from the regression of \(e_i^2/\hat \sigma ^2\), is larger at \(25.9\), because it assumes normal kurtosis and these residuals have much fatter tails; the gap illustrates Koenker’s point. White’s test, with \(\text {Pop}\) and \(\text {Pop}^2\), gives \(17.4\) on \(2\) degrees of freedom. The Goldfeld–Quandt test orders the states by population, drops the middle ten, and finds that the residual variance of the twenty largest states is \(21.7\) times that of the twenty smallest. (Both subsets have \(20 - 2 = 18\) degrees of freedom, so the ratio of sums of squares is the ratio of variances.)

Stata: estat hettest (Breusch–Pagan; add , iid for Koenker) and estat imtest, white. EViews: View \(\to \) Residual Diagnostics \(\to \) Heteroskedasticity Tests, choosing Breusch-Pagan-Godfrey, White or Harvey.

Software 10.5.3 (Weighted and feasible generalised least squares). Three choices of weights: variance proportional to population, to population squared (the per-capita model of Example 10.4.2), and estimated from the data:

> summary(lm(TI ~ Pop, data = st, weights = 1 / Pop))$coef
            Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -327.110  237.11681 -1.379531 1.741244e-01
Pop         4644.662   91.94615 50.515024 2.846440e-43
> summary(lm(TI ~ Pop, data = st, weights = 1 / Pop^2))$coef
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept)  170.1277   129.9784  1.308892 1.968028e-01
Pop         4327.1723   119.7114 36.146715 1.790659e-36
> v <- lm(log(e2) ~ log(Pop), data = st); coef(v)
(Intercept)    log(Pop)
 13.7423055   0.4956839
> fg <- lm(TI ~ Pop, data = st, weights = 1 / exp(fitted(v)))
> summary(fg)$coef
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -790.1068   337.9007 -2.338281 2.358648e-02
Pop         4784.5722    82.4403 58.036811 4.041469e-46

Reading the output. The slope moves from \(4\,905\) (least squares) to between \(4\,327\) and \(4\,785\) depending on the weights. These differences are larger than one would like if all the estimators were consistent for the same thing. They hint that the straight line is not the whole story, and that the largest states have a slightly different income per head; heteroscedasticity tests can be signalling misspecification. The estimated variance function says \(\sigma _i^2 \propto \text {Pop}^{0.50}\), much less steep than either assumed weighting, but \(\log e_i^2\) is a noisy dependent variable, and its slope has a large standard error. In practice one would report the feasible GLS estimates with robust standard errors, and check whether the conclusions depend on the choice of weights.

Stata: regress TI Pop [aweight = 1/Pop]; for feasible GLS, predict the fitted log variance from the auxiliary regression and weight by its exponential. EViews: Options \(\to \) Weights in the estimation dialog, or the ARCH/heteroskedasticity dialog for the multiplicative model.