import warnings
warnings.filterwarnings("ignore") # statsmodels emits harmless sqrt warnings on rare degenerate draws
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg") # non-interactive backend: render to buffers, never pop a window
import matplotlib.pyplot as plt
import statsmodels
import statsmodels.formula.api as smf
rng = np.random.default_rng(42) # ONE seed for the whole notebook -> reproducible Monte Carlo
print("numpy", np.__version__, "| pandas", pd.__version__, "| statsmodels", statsmodels.__version__)Notebook 5.4 — Replicating Petersen (2009): The Clustering Taxonomy
The one-sentence result. When the errors in a panel are correlated within firm, cluster by firm; when they are correlated within period, cluster by time; when both, cluster two ways — and any method whose assumed correlation structure does not match the truth reports a standard error that is too small, so you walk around overconfident.
In Chapter 5.4 you read Petersen (2009) as a paper and noticed its cleverest move: it is a Monte Carlo study, so the author identifies the truth by building it himself. He writes down a data-generating process where he sets the true slope \beta and the true error correlation, generates many artificial panels, and then asks of each standard-error method one blunt question — does its reported standard error match the actual sampling variability of \hat\beta? No appeal to a stylized fact, no argument from authority. Just: build a world where you know the answer, and watch which estimator recovers it.
This notebook hands you that machine. We will reproduce Petersen’s headline taxonomy with our own hands — Devon’s hands, really, since he keeps insisting his 100,000 on-chain transactions are “100,000 independent data points” and we are about to show him they are not. We build three worlds — a pure firm effect, a pure time effect, and both — run a small Monte Carlo in each, and tabulate the empirical 95% confidence-interval coverage of five methods: OLS, White (heteroskedasticity-robust), cluster-by-firm, cluster-by-time, and two-way clustering. The correct method should trap the true \beta about 95% of the time; the wrong ones, far less.
Petersen, M. A. (2009). Estimating Standard Errors in Finance Panel Data Sets: Comparing Approaches. Review of Financial Studies, 22(1), 435–480.
1. The data-generating process we control
Here is the world. We have N_{\text{firm}} firms observed over T years — a balanced panel of N = N_{\text{firm}} \times T firm-year rows. The truth is a clean line, y_{it} = \beta\, x_{it} + \varepsilon_{it}, \qquad \beta = 0.30 , with no intercept term in the truth (we still estimate one). The whole experiment lives in two design choices.
First, the regressor carries persistent structure. We build x_{it} from three pieces — a firm-level piece a_i (the same value for firm i in every year), a time-level piece b_t (shared by all firms in year t), and pure noise: x_{it} = a_i + b_t + \nu_{it}. This matters because of the Moulton logic from Ch 2.4: ignored within-cluster error correlation only inflates the danger when the regressor is also correlated within that cluster. A regressor that is pure noise within firm cannot be hurt by a firm effect, no matter how strong. By giving x both a firm and a time component, we let the damage bite along whichever error dimension we switch on.
Second, the error is a switchboard. We assemble \varepsilon_{it} from a firm component \gamma_i, a time component \delta_t, and idiosyncratic noise u_{it}:
\varepsilon_{it} = \gamma_i + \delta_t + u_{it}.
The structure argument decides which switches are on. "firm" keeps only \gamma_i (a firm mispriced this year tends to be mispriced next year). "time" keeps only \delta_t (a common shock — the whole market drops in 2008). "both" keeps both. That single switchboard is the entire taxonomy.
def simulate_panel(structure, r, n_firms=60, n_years=30, beta_true=0.30,
firm_sd=1.0, time_sd=1.0, idio_sd=1.0):
"""Generate one balanced panel from a KNOWN process.
structure : 'firm' (only gamma_i), 'time' (only delta_t), or 'both'.
r : a numpy Generator (passed in so the Monte Carlo controls the stream).
Returns a tidy DataFrame with columns y, x, firm, year.
"""
N = n_firms * n_years
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,...,T-1, 0,1,...,T-1, ...
# --- regressor: a persistent firm piece + a common time piece + noise ---
x_firm = np.repeat(r.normal(size=n_firms), n_years) # one level per firm
x_time = np.tile(r.normal(size=n_years), n_firms) # one level per year
x = x_firm + x_time + r.normal(size=N)
# --- error: switch components on/off by structure ---
use_firm = structure in ("firm", "both")
use_time = structure in ("time", "both")
gamma = np.repeat(r.normal(size=n_firms), n_years) if use_firm else np.zeros(N)
delta = np.tile(r.normal(size=n_years), n_firms) if use_time else np.zeros(N)
eps = firm_sd * gamma + time_sd * delta + idio_sd * r.normal(size=N)
y = beta_true * x + eps
return pd.DataFrame({"y": y, "x": x, "firm": firm, "year": year})
BETA_TRUE = 0.30
demo = simulate_panel("firm", rng, n_firms=60, n_years=30)
print(demo.head())
print(f"\nPanel shape: {demo['firm'].nunique()} firms x {demo['year'].nunique()} years "
f"= {len(demo)} firm-years")2. One fit, five standard errors
For any single simulated panel we run one pooled OLS regression — so there is exactly one \hat\beta — and then ask statsmodels for the standard error of x five different ways through get_robustcov_results. Each “flavor” keeps the point estimate identical and only changes its assumption about the error covariance matrix \boldsymbol\Omega:
- OLS (classical): \boldsymbol\Omega = \sigma^2 \mathbf I — equal variances, zero correlation.
- White (HC1): lets the diagonal of \boldsymbol\Omega vary (heteroskedasticity), but still assumes zero off-diagonal correlation. The paper’s big “gotcha”: robust does not mean robust to everything.
- Cluster by firm: allows arbitrary correlation among a firm’s own rows across years.
- Cluster by time: allows arbitrary correlation among all firms within a year.
- Two-way: allows both at once (
groupsis a two-column frame).
A trap worth flagging again.
get_robustcov_results(...)returns results whose.paramsand.bseare plain numpy arrays, not pandas Series indexed by name. Sores.bse["x"]fails. We grab the integer column position ofxonce and index by position everywhere.
def fit_ses(df):
"""Fit ONE pooled OLS and return beta_hat plus the SE of x under five methods."""
fit = smf.ols("y ~ x", data=df).fit()
j = list(fit.params.index).index("x") # integer position of the x coefficient
beta = float(np.asarray(fit.params)[j])
def se(res): # pull bse[x] by POSITION, never by name
return float(np.asarray(res.bse)[j])
return {
"beta": beta,
"OLS": se(fit),
"White": se(fit.get_robustcov_results("HC1")),
"Cluster firm": se(fit.get_robustcov_results("cluster", groups=df["firm"])),
"Cluster time": se(fit.get_robustcov_results("cluster", groups=df["year"])),
"Two-way": se(fit.get_robustcov_results("cluster", groups=df[["firm", "year"]])),
}
# Sanity check on a single firm-effect panel: cluster-firm SE should dwarf OLS/White.
one = fit_ses(simulate_panel("firm", rng))
print(f"beta_hat = {one['beta']:.3f} (true = {BETA_TRUE})")
for m in ["OLS", "White", "Cluster firm", "Cluster time", "Two-way"]:
print(f" {m:13s} SE = {one[m]:.4f}")3. The Monte Carlo engine
Now the move that makes the truth knowable. We draw S independent panels from the same known process. On each we record \hat\beta and the five standard errors. Then we compute two things Petersen puts side by side in his tables:
- the true standard error — the actual standard deviation of our S estimates \hat\beta, which we can compute because we have many draws; this is the ruler;
- the average estimated SE — the mean of each method’s reported \widehat{\text{se}} across draws; a correct method’s average sits on the ruler, a biased method’s sits below it.
And we compute coverage directly: across the S panels, what fraction of a method’s \hat\beta \pm 1.96\,\widehat{\text{se}} intervals actually trapped the true \beta = 0.30? A correct method covers about 95%; an overconfident one covers far less. We keep S = 400 so the whole notebook runs in well under a minute — modest, but plenty to see the taxonomy.
METHODS = ["OLS", "White", "Cluster firm", "Cluster time", "Two-way"]
def monte_carlo(structure, reps=500, seed=42, **panel_kw):
"""Run `reps` panels of a given structure; return a coverage/SE summary table."""
r = np.random.default_rng(seed) # fresh stream per call -> reproducible
betas = []
ses = {m: [] for m in METHODS}
hits = {m: 0 for m in METHODS}
for _ in range(reps):
res = fit_ses(simulate_panel(structure, r, **panel_kw))
betas.append(res["beta"])
for m in METHODS:
s = res[m]
ses[m].append(s)
lo, hi = res["beta"] - 1.96 * s, res["beta"] + 1.96 * s
if lo <= BETA_TRUE <= hi:
hits[m] += 1
true_se = float(np.std(betas, ddof=1)) # the ruler
rows = [{
"method": m,
"mean_SE": float(np.nanmean(ses[m])),
"true_SE": true_se,
"ratio": float(np.nanmean(ses[m])) / true_se, # < 1 means overconfident
"coverage": hits[m] / reps,
} for m in METHODS]
return pd.DataFrame(rows).set_index("method")
REPS = 400
print(f"Monte Carlo with S = {REPS} reps per world. Running...")4. World 1 — a pure firm effect
Exercise 1 from the reader’s guide. Errors are persistent within firm, no time component, and the regressor carries a firm component. This is the canonical corporate-finance case. Read the table column by column: anchor on true_SE (the ruler), then see how far each method’s mean_SE falls below it, and watch coverage collapse for the methods blind to the firm correlation.
tab_firm = monte_carlo("firm", reps=REPS).round(4)
print("WORLD 1: pure FIRM effect\n")
print(tab_firm.to_string())
print("\nReading: OLS and White SEs sit far BELOW the true SE (ratio << 1) -> they undercover.")
print("Cluster-by-FIRM matches the ruler and recovers ~95% coverage. Cluster-by-time fails.")5. World 2 — a pure time effect
Exercise 2: flip the structure. Drop the firm component, switch on a common shock \delta_t that hits all firms in a period. Now the symmetry of Petersen’s argument should appear — the method that was the hero in World 1, cluster-by-firm, becomes one of the failures, and cluster-by-time is the one that matches the ruler. Same code, same five methods; only the truth changed.
tab_time = monte_carlo("time", reps=REPS).round(4)
print("WORLD 2: pure TIME effect\n")
print(tab_time.to_string())
print("\nReading: now Cluster-by-FIRM undercovers and Cluster-by-TIME is the one near ~95%.")
print("The method whose assumed correlation matches the truth is the method that works.")5.1 An aside: time fixed effects largely absorb the time effect
The reader’s guide notes that fixed effects and clustering do different jobs. A time effect that lives entirely in an additive, mean-shifting \delta_t can be swept out by adding a dummy per year (a time fixed effect) — there is then little within-period correlation left for the standard error to worry about. We add year dummies via C(year) and check that even plain OLS standard errors behave far better than they did without the dummies. (This works precisely because our \delta_t is a pure level shift; richer time dependence would still need clustering.)
def coverage_with_time_fe(reps=400, seed=11, **panel_kw):
"""Time-effect world, but estimate y ~ x + C(year); report OLS-SE coverage of x."""
r = np.random.default_rng(seed)
hits, betas = 0, []
for _ in range(reps):
df = simulate_panel("time", r, **panel_kw)
fit = smf.ols("y ~ x + C(year)", data=df).fit() # year dummies sweep out delta_t
b = fit.params["x"]; s = fit.bse["x"]
betas.append(b)
if b - 1.96 * s <= BETA_TRUE <= b + 1.96 * s:
hits += 1
return hits / reps, float(np.std(betas, ddof=1))
cov_fe, tse_fe = coverage_with_time_fe()
print("Time-effect world, model y ~ x + C(year) (time fixed effects):")
print(f" plain OLS-SE coverage of x = {cov_fe:.3f} (vs ~{tab_time.loc['OLS','coverage']:.3f} "
f"with NO year dummies)")
print(" -> the year dummies absorb the common shock, so even classical SEs behave.")6. World 3 — both effects at once
Exercise 3: set the firm and time components to comparable sizes. Now neither one-way clustering is enough — clustering by firm ignores the within-period correlation and clustering by time ignores the within-firm correlation, so each one-way method still undercovers somewhat. Only two-way clustering, which allows both blocks of \boldsymbol\Omega to be dense at once, lands near 95%.
tab_both = monte_carlo("both", reps=REPS).round(4)
print("WORLD 3: BOTH effects\n")
print(tab_both.to_string())
print("\nReading: each ONE-way clustering still undercovers (it fixes only half the problem).")
print("TWO-WAY clustering is the one near ~95%. This is Petersen's 'both -> cluster two ways'.")7. The coverage taxonomy in one picture
Stack the three worlds into a single method × true-structure coverage table — this is Petersen’s taxonomy, condensed. Read down each column (a world): the cell nearest 0.95 names the correct clustering for that world, and it lies exactly where the method’s assumed structure matches the truth. Off that match, coverage falls — the further the assumption from the truth, the more overconfident you are.
cov_grid = pd.DataFrame({
"Firm effect": tab_firm["coverage"],
"Time effect": tab_time["coverage"],
"Both": tab_both["coverage"],
})
print("Empirical 95%-CI coverage (rows = SE method, cols = TRUE structure):\n")
print(cov_grid.round(3).to_string())
fig, ax = plt.subplots(figsize=(6.6, 4.4))
data = cov_grid.values
im = ax.imshow(data, cmap="RdYlGn", vmin=0.40, vmax=1.0, aspect="auto")
ax.set_xticks(range(cov_grid.shape[1])); ax.set_xticklabels(cov_grid.columns)
ax.set_yticks(range(cov_grid.shape[0])); ax.set_yticklabels(cov_grid.index)
ax.set_xlabel("TRUE correlation structure"); ax.set_ylabel("SE method")
ax.set_title("95%-CI coverage: green ≈ 0.95 (correct), red = overconfident")
for i in range(data.shape[0]):
for k in range(data.shape[1]):
ax.text(k, i, f"{data[i, k]:.2f}", ha="center", va="center",
color="black", fontsize=10)
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label="coverage")
fig.tight_layout()
plt.show()
print("\nThe green diagonal of 'right clustering' is the whole paper in one image.")8. Fama–MacBeth vs. cluster-by-time
The reader’s guide stresses one comparison among the historical fixes: Fama–MacBeth corrects for a time effect but not a firm effect. Fama–MacBeth runs a separate cross-sectional regression in each period, collects the T period slopes \hat\beta_t, and reports their average with a standard error \operatorname{sd}(\hat\beta_t)/\sqrt{T} and a t-distribution with T-1 degrees of freedom. Because it lets each period speak once, a common per-period shock is differenced away — so under a time effect it should cover correctly, neck-and-neck with cluster-by-time. But under a firm effect, a firm’s persistence leaks across the period regressions and Fama–MacBeth undercovers. We check both claims in the Monte Carlo.
from scipy import stats
def _ols_slope(yv, xv):
"""Tiny 2-column OLS (intercept + x); return the slope. No patsy -> fast."""
X = np.column_stack([np.ones_like(xv), xv])
beta, *_ = np.linalg.lstsq(X, yv, rcond=None)
return beta[1]
def fama_macbeth(df):
"""Return (avg slope, FM standard error, T) from period-by-period cross-sections."""
yv, xv, yr = df["y"].to_numpy(), df["x"].to_numpy(), df["year"].to_numpy()
slopes = [_ols_slope(yv[yr == t], xv[yr == t]) for t in np.unique(yr)]
slopes = np.asarray(slopes)
T = slopes.size
return slopes.mean(), slopes.std(ddof=1) / np.sqrt(T), T
def compare_fm(structure, reps=250, seed=7, **panel_kw):
"""Coverage of Fama-MacBeth vs cluster-by-time in a given world."""
r = np.random.default_rng(seed)
fm_hits = ct_hits = 0
for _ in range(reps):
df = simulate_panel(structure, r, **panel_kw)
b_fm, se_fm, T = fama_macbeth(df)
crit = stats.t.ppf(0.975, T - 1) # FM uses t with T-1 df, not 1.96
if b_fm - crit * se_fm <= BETA_TRUE <= b_fm + crit * se_fm:
fm_hits += 1
res = fit_ses(df)
s = res["Cluster time"]
if res["beta"] - 1.96 * s <= BETA_TRUE <= res["beta"] + 1.96 * s:
ct_hits += 1
return fm_hits / reps, ct_hits / reps
fm_time, ct_time = compare_fm("time")
fm_firm, ct_firm = compare_fm("firm")
fm_tab = pd.DataFrame({
"Fama-MacBeth": [fm_time, fm_firm],
"Cluster time": [ct_time, ct_firm],
}, index=["Time effect (FM should WORK)", "Firm effect (FM should FAIL)"]).round(3)
print("95%-CI coverage:\n")
print(fm_tab.to_string())
print("\nTime world: Fama-MacBeth ~= cluster-by-time, both near 0.95.")
print("Firm world: Fama-MacBeth undercovers -> it is the WRONG tool when firms persist.")9. The few-clusters caveat — break it on purpose
Petersen’s own caveat (Ch 5.4 §6): the cluster-robust estimator is consistent as the number of clusters grows, not as the number of observations grows. With only a handful of clusters it is itself noisy and biased downward, quietly reintroducing the overconfidence it was meant to cure. We stay in the time-effect world — where cluster-by-time is the correct method — and shrink the number of years (the number of time clusters) from 30 down to 5. Watch coverage fall back below nominal even though we are “doing it right.” This is the failure you only truly believe after you have made it happen on your own screen.
rows = []
for ny in [30, 20, 12, 8, 5]:
tab = monte_carlo("time", reps=250, n_years=ny)
rows.append({"n_years (= #time clusters)": ny,
"cluster-by-time coverage": round(tab.loc["Cluster time", "coverage"], 3)})
few = pd.DataFrame(rows).set_index("n_years (= #time clusters)")
print("Time-effect world, CORRECT method (cluster-by-time), shrinking the cluster count:\n")
print(few.to_string())
fig, ax = plt.subplots(figsize=(6.6, 4.0))
ax.plot(few.index, few["cluster-by-time coverage"], "o-", color="#3b6ea5")
ax.axhline(0.95, color="green", ls="--", lw=1, label="nominal 0.95")
ax.set_xlabel("number of time clusters (years)")
ax.set_ylabel("empirical coverage")
ax.set_title("Few-clusters caveat: 'right' clustering still undercovers when G is small")
ax.invert_xaxis(); ax.legend()
fig.tight_layout()
plt.show()
print("\nEven the correct clustering drifts below 0.95 as the cluster count shrinks.")
print("Below ~30-50 clusters, reach for the t(G-1) critical value or a wild cluster bootstrap.")Your Turn — push the taxonomy
You now own the author’s machine. Edit the cells above and re-run; each experiment is one keyword away.
- Tune the regressor’s persistence. Right now x carries equal firm and time pieces. Replace the regressor in
simulate_panelwithx = 0.1 * x_firm + x_time + r.normal(size=N)so x is nearly pure-noise within firm. Re-run World 1. Predict first: with a weak firm component in x, does the firm effect still wreck the OLS standard error? (Moulton says the bias scales with the regressor’s within-cluster correlation — so it should shrink.) - Crank the firm component. Pass
firm_sd=3.0tomonte_carlo("firm", ...). How low does OLS coverage go? Tabulate the OLSratio(mean_SE / true_SE) asfirm_sdclimbs 1, 2, 3, 4. - Two-way with few clusters in one dimension. Run World 3 with
n_firms=8butn_years=30. Two-way clustering needs enough clusters in both dimensions — does coverage hold, or does the thin firm dimension drag it down? - A referee question to carry forward. If a published paper clusters by industry with twelve industries, what would you ask the authors to report, and what alternative inference (from §9) would you demand?
When these stop surprising you, the recipe from Ch 2.4 has become something you know.