Notebook 2.4 — Standard-Error Flavors on a Clustered Panel

Maya has a clean-looking result and a problem she cannot name yet. Studying fair lending, she regressed loan interest rates on a borrower-tier variable and got a coefficient of \hat\beta_1 = 0.42 with a standard error of 0.07. The t-statistic is 0.42 / 0.07 = 6.0 — not “marginally significant,” but a coefficient screaming through a megaphone. Then her mentor asked one question: “Are those 9,000 loans really 9,000 independent pieces of evidence?”

They are not. Loans cluster by lender; lenders set rates in correlated ways. Once Maya accounts for that, the same \hat\beta_1 = 0.42 comes back with a standard error near 0.19 and the t-statistic collapses from 6.0 to about 2.2. Her point estimate never moved. What was wrong was the standard error — the number that tells you how much to trust the estimate.

The one-sentence result: when errors are heteroskedastic or correlated, \hat{\boldsymbol\beta} is still fine, but the classical standard-error formula lies, and you must replace it with a robust, clustered, or HAC formula.

In this notebook we build a finance panel — many firms over many years — with two diseases deliberately baked in: (a) heteroskedasticity (error scatter that varies across firms) and (b) within-firm correlation (each firm carries a persistent shock that makes its yearly errors hold hands). Then we compute, for the one and the same \hat\beta, seven different standard errors — classical, HC1, HC2, HC3, plain White (HC0), one-way cluster, two-way cluster, and Newey–West/HAC — and watch the t-statistic deflate as we stop pretending the errors are well behaved.

import numpy as npimport pandas as pdimport matplotlibmatplotlib.use("Agg")  # non-interactive backend: render to buffers, never pop a windowimport matplotlib.pyplot as pltimport statsmodelsimport statsmodels.api as smimport statsmodels.formula.api as smffrom linearmodels.panel import PooledOLSrng = np.random.default_rng(42)  # one seed for the whole notebook -> reproducibleprint("pandas", pd.__version__, "| statsmodels", statsmodels.__version__)

1. The sandwich, in one box

Everything below is the same variance formula, the sandwich:

\operatorname{Var}(\hat{\boldsymbol\beta}\mid\mathbf{X}) = \underbrace{(\mathbf{X}'\mathbf{X})^{-1}}_{\text{bread}}\, \underbrace{\mathbf{X}'\boldsymbol\Omega\mathbf{X}}_{\text{filling}}\, \underbrace{(\mathbf{X}'\mathbf{X})^{-1}}_{\text{bread}} .

The two outer factors are the bread; the middle \mathbf{X}'\boldsymbol\Omega\mathbf{X} is the filling, where \boldsymbol\Omega = \operatorname{Var}(\boldsymbol\varepsilon\mid\mathbf{X}) is the N\times N variance–covariance matrix of the errors. The classical formula is just the sandwich for the special diet \boldsymbol\Omega = \sigma^2\mathbf{I} — same variance on every diagonal, zeros off-diagonal — which collapses to \sigma^2(\mathbf{X}'\mathbf{X})^{-1}. Every “flavor” of standard error in this notebook keeps the bread identical and only changes the recipe for the filling:

White / HC lets the diagonal of \boldsymbol\Omega vary (heteroskedasticity).
Clustered lets whole blocks of \boldsymbol\Omega be dense (within-group correlation).
HAC / Newey–West lets a fading band around the diagonal be nonzero (serial correlation).

We will not “fix” \hat{\boldsymbol\beta} anywhere — there is nothing wrong with it. The whole game is estimating that filling honestly.

2. A four-point sandwich you can check by hand

Before any simulation, let us collapse the sandwich to a single regressor with no intercept, y_i = \beta x_i + \varepsilon_i, so \mathbf{X}'\mathbf{X} = \sum_i x_i^2 is one number and the sandwich for a diagonal \boldsymbol\Omega is plain arithmetic:

\operatorname{Var}(\hat\beta) = \frac{\sum_i x_i^2\,\sigma_i^2}{\left(\sum_i x_i^2\right)^2}.

Take x = (1,2,3,4), so \sum_i x_i^2 = 30. Compare a homoskedastic world (every \sigma_i^2 = 1) with a heteroskedastic one where the scatter grows with x (\sigma_i^2 = x_i^2). The classical formula would report 1/30 in both worlds — it averages the variances and forgets that the big-variance points are also the high-leverage points. The honest sandwich does not.

x = np.array([1.0, 2.0, 3.0, 4.0])
sxx = np.sum(x**2)                         # X'X = 30

# World A: homoskedastic, every sigma_i^2 = 1
varA = np.sum(x**2 * 1.0) / sxx**2
# World B: heteroskedastic, sigma_i^2 = x_i^2  (scatter grows with x)
varB = np.sum(x**2 * x**2) / sxx**2

classical = 1.0 / sxx                      # what classical reports in BOTH worlds

print(f"X'X = sum(x^2)          = {sxx:.0f}")
print(f"classical Var(beta-hat) = {classical:.4f}   (reported in both worlds)")
print(f"sandwich  Var, World A  = {varA:.4f}   (homoskedastic -> matches classical)")
print(f"sandwich  Var, World B  = {varB:.4f}   (heteroskedastic -> {varB/classical:.1f}x larger)")
print(f"=> classical SE understates the truth by a factor of {np.sqrt(varB/classical):.1f} in World B")

So a reported t-stat of 3.5 in World B is really a t-stat of 1. That is the entire heteroskedasticity problem in four data points, with no calculus and no asymptotics — just plugging the right \sigma_i^2 into the middle of the sandwich. The rest of the notebook is about (a) what happens when off-diagonals are also nonzero, and (b) how to estimate the filling when nobody hands you the \sigma_i^2.

3. Build a finance panel with both diseases

Now the realistic case Maya and Priya actually live in: a panel of many firms i observed over many years t. We simulate 120 firms \times\,20 years (N = 2{,}400 firm-years) with a key regressor x (think: a leverage or risk measure) and an outcome y (think: a default-risk or pricing measure). The true slope is \beta = 0.30.

We deliberately engineer two failures of \boldsymbol\Omega = \sigma^2\mathbf{I}:

A persistent firm effect. Each firm draws one firm_shock that is added to every one of its yearly errors. That makes a firm’s errors correlated across all its years — the dense diagonal blocks of Section 4 in the chapter, the dominant problem in corporate-finance panels (Petersen 2009).
Heteroskedasticity. The idiosyncratic error scale is \exp(0.9\,|x_{\text{firm}}|), so some firms are far noisier than others — the varying diagonal of \boldsymbol\Omega.

We also add a milder common time effect (year_shock, shared by all firms in a year) so we have something for two-way clustering to catch. The regressor x itself carries a persistent firm-level piece (x_firm), which — per the Moulton logic — is exactly what makes ignoring clustering so costly.

One technical move: we residualize the error vector against x before forming y. This forces the in-sample correlation between x and the error to be exactly zero, so \hat\beta lands exactly on the true 0.30. That is not cheating — it just removes sampling noise in \hat\beta so the notebook’s whole point is unmistakable: as we change the SE flavor, the estimate sits perfectly still and only its uncertainty moves.

n_firms, n_years = 120, 20
N = n_firms * n_years
BETA_TRUE = 0.30

# panel index columns
firm = np.repeat(np.arange(n_firms), n_years)   # 0,0,...,0, 1,1,...,1, ...
year = np.tile(np.arange(n_years), n_firms)     # 0,1,...,19, 0,1,...,19, ...

# --- the key regressor: a persistent firm-level piece + idiosyncratic noise ---
x_firm = np.repeat(rng.normal(size=n_firms), n_years)   # one level per firm
x = x_firm + 0.6 * rng.normal(size=N)

# --- the error: persistent firm effect + common time effect + heteroskedastic idio ---
firm_shock = np.repeat(rng.normal(size=n_firms), n_years)  # persistent within-firm shock
year_shock = np.tile(rng.normal(size=n_years), n_firms)    # common shock per year
het_scale  = np.exp(0.9 * np.abs(x_firm))                  # noisier for some firms
idio       = 1.0 * het_scale * rng.normal(size=N)          # heteroskedastic idiosyncratic
eps        = 1.3 * firm_shock + 0.4 * year_shock + idio

# residualize eps on [1, x] so sample corr(x, eps) = 0  ->  beta-hat == BETA_TRUE exactly
Xmat = np.column_stack([np.ones(N), x])
eps = eps - Xmat @ np.linalg.lstsq(Xmat, eps, rcond=None)[0]

y = BETA_TRUE * x + eps

df = pd.DataFrame({"y": y, "x": x, "firm": firm, "year": year})
print(df.head())
print(f"\nPanel: {n_firms} firms x {n_years} years = {N} firm-years")

3.1 Confirm the diseases are really in the data

A simulation is only as good as the structure it actually contains. Two quick checks. First, the within-firm error correlation: if firm errors hold hands, then a firm’s average residual should be a big chunk of its individual residuals — formally, the between-firm variance of errors should be a large share of the total. Second, heteroskedasticity: the firm-level error standard deviation should vary a lot across firms.

tmp = df.copy()
tmp["eps"] = eps

# (1) within-firm correlation: share of error variance that lives BETWEEN firms.
firm_mean_eps = tmp.groupby("firm")["eps"].transform("mean")
between_share = np.var(firm_mean_eps, ddof=0) / np.var(tmp["eps"], ddof=0)

# (2) heteroskedasticity: spread of each firm's own error standard deviation.
firm_sd = tmp.groupby("firm")["eps"].std(ddof=1)

print(f"Between-firm share of error variance : {between_share:.2%}")
print(f"  (0% would mean no firm effect; we built in a big one)")
print(f"Firm-level error SD: min={firm_sd.min():.2f}, "
      f"median={firm_sd.median():.2f}, max={firm_sd.max():.2f}")
print(f"  (a wide range => heteroskedasticity across firms)")

3.2 The residual fan

The cheapest heteroskedasticity diagnostic is one plot: residuals against a regressor. Under homoskedasticity you see a structureless horizontal band of constant thickness; under heteroskedasticity you see a fan that widens as you move along the axis. Let us fit the pooled OLS once and look.

fit = smf.ols("y ~ x", data=df).fit()   # ONE OLS fit, reused for every SE flavor below
resid = fit.resid

fig, ax = plt.subplots(figsize=(7, 4.2))
ax.scatter(df["x"], resid, s=8, alpha=0.35, color="#3b6ea5")
ax.axhline(0, color="black", lw=0.8)
ax.set_xlabel("regressor  x")
ax.set_ylabel(r"OLS residual  $\hat\varepsilon_i$")
ax.set_title("Residual fan: scatter widens with |x|  (heteroskedasticity)")
fig.tight_layout()
plt.show()
print("Note the funnel shape: errors are far noisier at large |x|.")

4. One estimate, seven standard errors

Now the heart of the notebook. We have one fitted model, fit, and therefore one \hat\beta. We will ask statsmodels for the covariance matrix seven different ways via get_robustcov_results, and read the standard error of x from each.

A trap worth flagging. get_robustcov_results(...) returns results whose .params and .bse are plain numpy arrays, not pandas Series indexed by variable name. So res.bse["x"] will fail. We grab the column position of x once and index by integer everywhere.

xi = list(fit.params.index).index("x")   # integer column position of the x coefficient
print(f"'x' is coefficient number {xi} (0 is the intercept)")

def se_row(res, label):
    '''Pull (beta, se, t) for the x coefficient from a (possibly robust) result.

    get_robustcov_results returns numpy arrays, not pandas Series, so we index by
    POSITION (xi), never by name.
    '''
    b = float(np.asarray(res.params)[xi])
    s = float(np.asarray(res.bse)[xi])
    return {"flavor": label, "beta": b, "se": s, "t": b / s}

results = [
    se_row(fit,                                                   "Classical"),
    se_row(fit.get_robustcov_results("HC0"),                      "White HC0"),
    se_row(fit.get_robustcov_results("HC1"),                      "White HC1"),
    se_row(fit.get_robustcov_results("HC2"),                      "White HC2"),
    se_row(fit.get_robustcov_results("HC3"),                      "White HC3"),
    se_row(fit.get_robustcov_results("cluster", groups=df["firm"]),
           "Cluster: firm"),
    se_row(fit.get_robustcov_results("cluster", groups=df[["firm", "year"]]),
           "Two-way: firm+year"),
    se_row(fit.get_robustcov_results("HAC", maxlags=4),           "Newey-West HAC(4)"),
]

table = pd.DataFrame(results).set_index("flavor")
table = table.round({"beta": 3, "se": 4, "t": 2})
table

Read down the beta column first: it is 0.300 in every single row. The point estimate does not budge. Now read the t column. The classical t-stat is around 5.7 — a megaphone. The White (HC) flavors, which fix heteroskedasticity but still assume errors are uncorrelated across observations, knock it down to roughly 4.2. But the big move comes from clustering by firm: once we admit that a firm’s 20 yearly errors are not 20 independent facts, the t-stat collapses to about 2.2. Same \hat\beta, same data — the megaphone is now a normal speaking voice. That is Maya’s 6.0 \to 2.2, reproduced on simulated data you built yourself.

Notice the ordering of the standard errors: classical < White < HAC < clustered. The clustered SE is biggest because the firm effect is the dominant source of correlation, and only clustering accounts for it. The HC flavors barely differ from one another here because with N = 2{,}400 and no extreme leverage, the finite-sample corrections HC0–HC3 have little to correct.

4.1 The picture: standard error by method

order = ["Classical", "White HC0", "White HC1", "White HC2", "White HC3",
         "Newey-West HAC(4)", "Cluster: firm", "Two-way: firm+year"]
ses = table.loc[order, "se"].values
colors = ["#b0b0b0", "#7fb069", "#7fb069", "#7fb069", "#7fb069",
          "#e0a458", "#c1432f", "#8e2b1f"]

fig, ax = plt.subplots(figsize=(8.5, 4.6))
bars = ax.bar(range(len(order)), ses, color=colors)
ax.axhline(table.loc["Classical", "se"], color="#b0b0b0", ls="--", lw=1,
           label="classical (naive) level")
ax.set_xticks(range(len(order)))
ax.set_xticklabels(order, rotation=30, ha="right")
ax.set_ylabel(r"standard error of $\hat\beta_x$")
ax.set_title("Same estimate, very different uncertainty: SE by method")
for b_, v in zip(bars, ses):
    ax.text(b_.get_x() + b_.get_width()/2, v + 0.002, f"{v:.3f}",
            ha="center", va="bottom", fontsize=8)
ax.legend()
fig.tight_layout()
plt.show()

5. Cross-check the firm cluster with `linearmodels`

It is good hygiene to confirm the firm-clustered number with a second library. The linearmodels package was built for panels; its PooledOLS with cluster_entity=True clusters on the entity (here, firm) index. We set a MultiIndex of (firm, year) so it knows the panel structure, then compare to the statsmodels number above. They should agree to within a tiny finite-sample adjustment.

pdf = df.set_index(["firm", "year"])          # (entity, time) MultiIndex for linearmodels
lm = PooledOLS(pdf["y"], sm.add_constant(pdf[["x"]]))
lm_fit = lm.fit(cov_type="clustered", cluster_entity=True)

se_sm = table.loc["Cluster: firm", "se"]
se_lm = float(lm_fit.std_errors["x"])
print(f"statsmodels firm-cluster SE(x)   : {se_sm:.4f}")
print(f"linearmodels firm-cluster SE(x)  : {se_lm:.4f}")
print(f"difference                       : {abs(se_sm - se_lm):.4f}  (a small df adjustment)")

6. The few-clusters caveat

The cluster-robust estimator is consistent as the number of clusters G \to \infty, not as N \to \infty. With few clusters the clustered SE is itself noisy and biased downward — overconfident again, the very disease we are treating. The rule of thumb is to want at least 30–50 clusters before trusting the clustered t-test at face value.

Let us see the SE get unreliable. We re-estimate the firm-clustered SE using only the first G firms, for shrinking G, and watch it stop being stable.

def firm_cluster_se(sub):
    '''Refit OLS on a sub-panel and return the firm-clustered SE of x.'''
    f = smf.ols("y ~ x", data=sub).fit()
    j = list(f.params.index).index("x")
    rc = f.get_robustcov_results("cluster", groups=sub["firm"])
    return float(np.asarray(rc.bse)[j])

G_grid = [120, 80, 50, 30, 15, 8, 4]
rows = []
for G in G_grid:
    sub = df[df["firm"] < G].copy()        # keep first G firms; .copy() avoids chained assignment
    rows.append({"G_clusters": G, "n_obs": len(sub), "cluster_SE": firm_cluster_se(sub)})

few = pd.DataFrame(rows).set_index("G_clusters").round(4)
few

fig, ax = plt.subplots(figsize=(7, 4.2))
ax.plot(few.index, few["cluster_SE"], "o-", color="#8e2b1f")
ax.axvspan(0, 30, color="red", alpha=0.07)
ax.text(15, few["cluster_SE"].max()*0.9, "danger zone\n(G < 30)",
        ha="center", color="#8e2b1f", fontsize=9)
ax.set_xlabel("number of clusters  G  (firms kept)")
ax.set_ylabel("firm-clustered SE of  x")
ax.set_title("With few clusters the clustered SE gets jumpy and unreliable")
ax.invert_xaxis()  # read left-to-right as clusters shrink
fig.tight_layout()
plt.show()
print("As G shrinks below ~30 the SE wanders: the asymptotics live in G, not N.")

Your Turn — crank up the intra-cluster correlation

The naive (classical) SE understates the truth more the stronger the within-cluster correlation. We packaged the whole simulation into a function whose only knob is firm_shock_strength — the size of the persistent firm effect, i.e. how tightly each firm’s yearly errors hold hands. As you turn it up, the classical SE barely moves (it is blind to clustering) while the clustered SE climbs, so the ratio between them — the factor by which the naive number lies — grows.

Run the sweep below, read the rising ratio column, then try the experiments in the final cell.

def simulate_panel(firm_shock_strength, n_firms=120, n_years=20, seed=42):
    '''Rebuild the panel with a tunable persistent firm effect; return (classical_SE, cluster_SE).'''
    r = np.random.default_rng(seed)
    N = n_firms * n_years
    firm = np.repeat(np.arange(n_firms), n_years)
    x_firm = np.repeat(r.normal(size=n_firms), n_years)
    x = x_firm + 0.6 * r.normal(size=N)
    fshock = np.repeat(r.normal(size=n_firms), n_years)
    yshock = np.tile(r.normal(size=n_years), n_firms)
    idio = np.exp(0.9 * np.abs(x_firm)) * r.normal(size=N)
    e = firm_shock_strength * fshock + 0.4 * yshock + idio
    Xm = np.column_stack([np.ones(N), x])
    e = e - Xm @ np.linalg.lstsq(Xm, e, rcond=None)[0]   # clean beta-hat
    d = pd.DataFrame({"y": 0.30 * x + e, "x": x, "firm": firm})
    f = smf.ols("y ~ x", data=d).fit()
    j = list(f.params.index).index("x")
    s_cl = float(np.asarray(f.bse)[j])                                   # classical
    s_cluster = float(np.asarray(
        f.get_robustcov_results("cluster", groups=d["firm"]).bse)[j])    # firm-clustered
    return s_cl, s_cluster

sweep = []
for strength in [0.0, 0.5, 1.0, 1.5, 2.5, 4.0]:
    s_cl, s_cluster = simulate_panel(strength)
    sweep.append({"firm_shock_strength": strength,
                  "classical_SE": s_cl, "cluster_SE": s_cluster,
                  "ratio cluster/classical": s_cluster / s_cl})

sweep_df = pd.DataFrame(sweep).round(3)
print(sweep_df.to_string(index=False))
print("\nAs the firm effect grows, the classical SE barely moves but the clustered SE")
print("climbs -> the naive formula understates the truth more and more.")

Experiments to try (edit the cell above and re-run):

More observations, same clusters. In simulate_panel, raise n_years from 20 to 60 while keeping n_firms fixed. Does the clustered SE shrink the way the classical one does? (It barely does — piling more years onto the same firms does not buy independent evidence. Evidence comes from more firms.)
Kill the firm effect. Set firm_shock_strength = 0.0. The clustered and classical SEs should nearly coincide — when the errors really are well behaved, you lose almost nothing by defaulting to the honest formula.
Re-derive Maya’s collapse. Find a firm_shock_strength that makes ratio cluster/classical about 2.7. With a classical t-stat near 6, that ratio pushes the clustered t-stat down to roughly 6 / 2.7 \approx 2.2 — exactly the number Maya’s mentor forced her to confront.
Few clusters again. Call simulate_panel(2.5, n_firms=8) and compare its clustered SE to the n_firms=120 version. With only 8 clusters the number is no longer trustworthy — reach for the wild cluster bootstrap (Cameron, Gelbach & Miller 2008) in real work.

The lesson, one more time: getting \hat\beta right and getting its uncertainty right are two separate jobs. This whole notebook fixed only the second one — and it moved a t-stat from “publish” to “not yet” without the estimate ever budging.