Heteroskedasticity, Dependence, and Robust Inference

1. Heteroskedasticity-consistent covariance matrices

From Lecture 2, \(\operatorname{avar}(\widehat\beta)=Q^{-1}\Omega Q^{-1}/n\) with \(\Omega=E[X_iX_i'u_i^2]\). Conditional on \(X\), with independent observations,

The White (1980) insight is that we do not need to estimate each \(\sigma_i^2\) consistently. We only need the \(k\times k\) average \(n^{-1}\sum_iX_iX_i'\sigma_i^2\), and replacing \(\sigma_i^2\) by \(\widehat u_i^2\) does this.

Why leverage matters. Under homoskedasticity, \(E[\widehat u_i^2]=\sigma^2(1-h_{ii})\). Residuals at high-leverage points are shrunk because the fit chases those points. HC2 undoes this exactly; HC3 over-corrects, which protects against the remaining downward bias under heteroskedasticity. With \(n\) in the thousands and no extreme leverage, the four estimators agree. With small samples or a few influential observations, they can differ substantially. Report HC3, or check that conclusions survive it (Long and Ervin, 2000; MacKinnon and White, 1985).

2. Testing for heteroskedasticity, and WLS

Breusch–Pagan / White tests. Regress \(\widehat u_i^2\) on a set of variables \(Z_i\) (the regressors for BP, regressors plus squares and cross-products for White) and use \(nR^2\to_d\chi^2_{\dim Z-1}\) under homoskedasticity. A rejection says the classical variance formula is wrong. Because robust standard errors are valid either way, the test matters mainly for deciding whether efficiency gains from reweighting are worthwhile.

Weighted least squares. If \(\operatorname{Var}(u_i\mid X_i)=\sigma^2h(X_i)\) with \(h\) known, OLS on the transformed model \(Y_i/\sqrt{h_i}=(X_i/\sqrt{h_i})'\beta+u_i/\sqrt{h_i}\) is BLUE (GLS). With \(h(X_i;\gamma)\) estimated in a first step, feasible GLS is asymptotically efficient if the variance model is correct. Practical advice (Romano and Wolf, 2017): use WLS for efficiency and robust standard errors for validity, because the variance model is almost never exactly right. Note that WLS and OLS estimate the same \(\beta\) only when the CEF is linear. Under misspecification, they converge to different weighted projections.

3. Clustered dependence

Suppose observations come in \(G\) groups (villages, schools, firms, states), indexed by \(g\), with \(n_g\) members. Errors are independent across clusters but arbitrarily correlated within them. Writing \(X_g\) and \(u_g\) for the cluster's stacked data,

The cluster-robust (CR) estimator (Liang and Zeger, 1986; Arellano, 1987) replaces \(E[u_gu_g']\) by \(\widehat u_g\widehat u_g'\):

It is consistent as \(G\to\infty\) with cluster sizes bounded. The effective sample size for inference is \(G\), not \(n\).

The Moulton factor. Why does clustering matter so much? Take equal cluster sizes \(m\), a regressor \(x\) with within-cluster correlation \(\rho_x\), and an error with within-cluster correlation \(\rho_u\). Then

When the regressor is constant within clusters (\(\rho_x=1\)), as with a state-level policy, and \(m=100\) individuals per state with \(\rho_u=0.05\), the variance ratio is about \(5.95\). Conventional standard errors are then too small by a factor of about \(\sqrt{5.95}=2.4\) (Moulton, 1990).

Few clusters. With \(G\) below roughly 30–50, or with very unbalanced clusters, CR \(t\)-tests over-reject even with the \(c\) correction. Remedies:

• use \(t_{G-1}\) critical values rather than normal ones;
• use the CR2 estimator (Bell–McCaffrey) with Imbens–Kolesár degrees of freedom;
• use the wild cluster bootstrap with the null imposed (Cameron, Gelbach and Miller, 2008; Roodman et al., 2019). With very few treated clusters, every method struggles; randomization inference may be the honest option (MacKinnon and Webb, 2017).

Where to cluster. Abadie, Athey, Imbens and Wooldridge (2023) show that clustering is justified by the design: sampling clusters from a population (sampling-based clustering) or assigning treatment at the cluster level (design-based clustering). Clustering "because the errors might be correlated" at the coarsest available level is not always conservative, and it can be very imprecise. Rule of thumb: cluster at the level at which treatment is assigned. In panels, cluster at least at the unit level to allow arbitrary serial correlation (Bertrand, Duflo and Mullainathan, 2004).

4. Serial correlation: HAC estimation

In a time-series regression \(y_t=x_t'\beta+u_t\), the score \(s_t=x_tu_t\) may be serially correlated. Then

the long-run variance: the spectral density of \(s_t\) at frequency zero, times \(2\pi\). The sample autocovariances \(\widehat\Gamma_j=T^{-1}\sum_t\widehat s_t\widehat s_{t-j}'\) cannot simply be summed over all lags: the result is not consistent and need not be positive semidefinite. Newey and West (1987) down-weight with a Bartlett kernel:

which is positive semidefinite by construction and consistent if \(L\to\infty\) with \(L/T\to0\).

Bandwidth. A common rule is \(L=\lfloor0.75\,T^{1/3}\rfloor\) (Stock and Watson). Andrews (1991) gives data-dependent choices. Small \(L\) with persistent scores gives size distortions; large \(L\) gives noisy estimates. Recent work recommends larger bandwidths combined with fixed-\(b\) critical values, or equal-weighted cosine (EWC) estimators with \(t\) critical values (Lazarus, Lewis, Stock and Watson, 2018).

HAC does not fix dynamic misspecification. If \(y_t\) depends on \(y_{t-1}\) and the lag is omitted, the problem is in the conditional mean, not just the variance. In particular, a regression with a lagged dependent variable and serially correlated errors is inconsistent (Lecture 9).

5. A unified view

All of the estimators above have the same sandwich form \(\widehat Q^{-1}\widehat\Omega\widehat Q^{-1}\). They differ only in the "meat":

import numpy as np
rng = np.random.default_rng(4)

def sim(G=40, m=50, rho_u=0.1):
    g = np.repeat(np.arange(G), m)
    d = (rng.random(G) < 0.5).astype(float)[g]               # treatment assigned by cluster
    u = np.sqrt(rho_u) * rng.normal(size=G)[g] + np.sqrt(1 - rho_u) * rng.normal(size=G * m)
    y = 0.0 * d + u                                         # true effect = 0
    X = np.column_stack([np.ones_like(d), d])
    XtXi = np.linalg.inv(X.T @ X)
    b = XtXi @ X.T @ y
    e = y - X @ b
    V_hc = XtXi @ (X * e[:, None]**2).T @ X @ XtXi
    S = np.zeros((2, 2))
    for c in range(G):
        sc = X[g == c].T @ e[g == c]
        S += np.outer(sc, sc)
    n, k = len(y), 2
    V_cr = XtXi @ S @ XtXi * G / (G - 1) * (n - 1) / (n - k)
    return abs(b[1] / np.sqrt(V_hc[1, 1])) > 1.96, abs(b[1] / np.sqrt(V_cr[1, 1])) > 2.02  # t_39 critical value

rej = np.array([sim() for _ in range(2000)]).mean(axis=0)
print("rejection rate under H0: HC =", rej[0], " cluster =", rej[1])

import numpy as np, pandas as pd, wooldridge as wd, statsmodels.formula.api as smf
from statsmodels.stats.diagnostic import het_breuschpagan, het_white
df = wd.data("hprice1")
lev = smf.ols("price ~ lotsize + sqrft + bdrms", data=df).fit()
# TODO: het_breuschpagan(lev.resid, lev.model.exog); compare lev.get_robustcov_results(cov_type=...)

import numpy as np, pandas as pd, wooldridge as wd, statsmodels.formula.api as smf
from statsmodels.stats.diagnostic import het_breuschpagan, het_white
df = wd.data("hprice1")
lev = smf.ols("price ~ lotsize + sqrft + bdrms", data=df).fit()
bp, wh = het_breuschpagan(lev.resid, lev.model.exog), het_white(lev.resid, lev.model.exog)
print(f"levels: Breusch-Pagan LM = {bp[0]:.2f} (p = {bp[1]:.4f});  White LM = {wh[0]:.2f} (p = {wh[1]:.4f})")

ses = {"classical": lev.bse}
for c in ("HC0", "HC1", "HC2", "HC3"):
    ses[c] = pd.Series(lev.get_robustcov_results(cov_type=c).bse, index=lev.params.index)
print(pd.DataFrame(ses).round(4))

lg = smf.ols("np.log(price) ~ np.log(lotsize) + np.log(sqrft) + bdrms", data=df).fit()
bp2 = het_breuschpagan(lg.resid, lg.model.exog)
print(f"\nlogs:   Breusch-Pagan LM = {bp2[0]:.2f} (p = {bp2[1]:.4f})")

w = smf.wls("price ~ lotsize + sqrft + bdrms", data=df, weights=1 / df["sqrft"]).fit(cov_type="HC1")
print("\nWLS (weights 1/sqrft), robust se:")
print(pd.DataFrame({"coef": w.params, "robust se": w.bse}).round(4))

levels: Breusch-Pagan LM = 14.09 (p = 0.0028);  White LM = 33.73 (p = 0.0001)
           classical      HC0      HC1      HC2      HC3
Intercept    29.4750  36.2843  37.1382  38.3813  41.0327
lotsize       0.0006   0.0012   0.0013   0.0029   0.0071
sqrft         0.0132   0.0173   0.0177   0.0226   0.0407
bdrms         9.0101   8.2837   8.4786   9.1866  11.5618

logs:   Breusch-Pagan LM = 4.22 (p = 0.2383)

WLS (weights 1/sqrft), robust se:
              coef  robust se
Intercept   4.1989    35.1199
lotsize     0.0016     0.0008
sqrft       0.1178     0.0173
bdrms      10.6073     9.1660

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
df = wd.data("wagepan")
f = "lwage ~ union + educ + black + hisp + exper + expersq + married + C(year)"
# TODO: fit(cov_type="cluster", cov_kwds={"groups": df["nr"]})

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
df = wd.data("wagepan")
f = "lwage ~ union + educ + black + hisp + exper + expersq + married + C(year)"
m = smf.ols(f, data=df)
se = {"classical": m.fit().bse["union"], "HC1": m.fit(cov_type="HC1").bse["union"],
      "cluster(nr)": m.fit(cov_type="cluster", cov_kwds={"groups": df["nr"]}).bse["union"]}
print({k: round(v, 4) for k, v in se.items()}, " coef:", round(m.fit().params["union"], 4))

def within_corr(v):
    w = v.to_numpy().reshape(-1, 8)                 # 545 workers x 8 years (sorted by nr, year)
    w = w - w.mean()
    num = sum(np.sum(np.outer(r, r)) - np.sum(r**2) for r in w)
    return num / (7 * np.sum(w**2))
df["res"] = m.fit().resid
rho_u, rho_x = within_corr(df["res"]), within_corr(df["union"].astype(float))
print(f"within-worker correlation: residual {rho_u:.3f}, union {rho_x:.3f}")
print(f"variance ratio cluster/HC1 = {(se['cluster(nr)']/se['HC1'])**2:.2f};  Moulton 1+(m-1)rho_x rho_u = {1 + 7*rho_x*rho_u:.2f}")

{'classical': np.float64(0.0172), 'HC1': np.float64(0.0163), 'cluster(nr)': np.float64(0.0274)}  coef: 0.1825
within-worker correlation: residual 0.458, union 0.528
variance ratio cluster/HC1 = 2.84;  Moulton 1+(m-1)rho_x rho_u = 2.69