Research Design, Monte Carlo, Replication, and Reporting

1. From question to estimand

A good empirical paper answers four questions, in order.

1. Estimand. What population quantity do you want? "The effect of class size on test scores" is not yet an estimand. "The average effect, on third-grade reading scores, of reducing class size from 22 to 15 for students in Tennessee public schools in 1985–89" is one.
2. Identification. Under which assumptions does that quantity equal a feature of the observable distribution (a regression coefficient, a Wald ratio, a discontinuity)? State them in words and in notation, and say which can be partially tested.
3. Estimation and inference. Which estimator, which variance estimator (matching the dependence structure, Lecture 4), and which finite-sample corrections?
4. Validation. Placebo outcomes, placebo treatments, pre-trends, balance, density tests, sensitivity analysis.

Pre-analysis plans. For experiments, and increasingly for observational work using new data, register the primary outcomes, specifications and hypotheses before seeing the outcome data (the AEA RCT Registry, OSF). A plan does not forbid exploration. It makes the distinction between confirmatory and exploratory results visible (Olken, 2015).

2. Monte Carlo experiments

Purpose. Simulation answers questions that asymptotic theory leaves open: Is the bias of my estimator material at \(n=200\)? Do my 95% intervals cover? How many clusters do I need?

Design checklist.

• Data-generating process (DGP). Make it realistic: calibrate it to your data (sample size, number of clusters, regressor distributions, \(R^2\), error dependence). Vary the features that theory says matter (instrument strength, persistence, heteroskedasticity).
• Estimators and procedures compared on the same simulated datasets (common random numbers reduce the noise in comparisons).
• Performance measures for an estimator \(\widehat\theta_r\), \(r=1,\dots,R\):

Report median bias and interquartile ranges when the estimator may lack moments (just-identified IV has no finite mean).

• Monte Carlo standard errors. Simulation results are estimates. For a rejection rate \(\widehat p\), the MC standard error is \(\sqrt{\widehat p(1-\widehat p)/R}\). With \(R=1{,}000\) and true size 5%, it is 0.0069, so an observed 6.2% is not evidence of size distortion. To detect a 1-point distortion reliably, you need \(R\approx10{,}000\).
• Reproducibility. Fix and report seeds, save simulated results, and parallelize with independent random streams (for example, numpy.random.SeedSequence.spawn).

3. The bootstrap

The bootstrap approximates the sampling distribution of a statistic by resampling from the data (Efron, 1979). The resampling scheme must mimic the dependence structure of the data.

Asymptotic refinement. Bootstrapping an asymptotically pivotal statistic, such as a \(t\)-statistic, gives more accurate tests and intervals than normal approximations (the percentile-\(t\) method). Bootstrapping the estimator itself (the percentile method) does not provide this refinement.

When the bootstrap fails. Examples: parameters on the boundary of the parameter space, the maximum of a sample, non-differentiable functionals, nearest-neighbor matching estimators (Abadie and Imbens, 2008), weak instruments, and unit-root autoregressions (\(\phi=1\)). Subsampling or the \(m\)-out-of-\(n\) bootstrap can help in some of these cases.

4. Multiple testing and specification robustness

Multiple hypotheses. Testing \(m\) hypotheses at 5% each gives a probability of at least one false rejection of up to \(1-0.95^m\): 40% for \(m=10\).

• Family-wise error rate (FWER). Bonferroni tests each hypothesis at \(\alpha/m\). Holm's step-down procedure is uniformly more powerful: sort p-values \(p_{(1)}\le\dots\le p_{(m)}\) and reject \(H_{(j)}\) while \(p_{(j)}\le\alpha/(m-j+1)\). Romano–Wolf (2005) uses the bootstrap to exploit the correlation between test statistics, which gives more power when outcomes are correlated.
• False discovery rate (FDR). Benjamini–Hochberg controls the expected share of false rejections among rejections: reject \(H_{(1)},\dots,H_{(k)}\) for the largest \(k\) with \(p_{(k)}\le k\alpha/m\). It is appropriate for exploratory analysis with many outcomes. Anderson (2008) applies these methods to early-childhood interventions.
• Index outcomes. Combine related outcomes into a standardized summary index (Kling, Liebman and Katz, 2007), which reduces the number of tests and increases power.

Specification robustness. Report how the main estimate changes across defensible choices of controls, samples, functional forms and estimators: a "specification curve" or multiverse analysis (Simonsohn, Simmons and Nelson, 2020). Robustness tables that only show specifications that work are not robustness checks.

5. Reporting standards

• Lead with effect sizes and confidence intervals in meaningful units (percent, standard deviations, dollars per person). Interpret magnitude against a benchmark.
• Avoid significance stars as the main message. "Not significant" is not "no effect": report the interval and what it rules out.
• Report the first-stage \(F\) (effective \(F\)) for IV, the number of clusters, the variance estimator, and the number of treated units.
• Show the raw data where possible: binned scatter plots, event-study plots with pre-periods, RD plots with binned means.
• Be explicit about the estimand. LATE, ATT and ATE differ, and so do weighted and unweighted estimates.

6. Replication packages

Modern journals require a replication package that reproduces every table and figure from raw data. The AEA Data and Code Availability Policy and the Social Science Data Editors' template README set the standard. A package should contain:

1. a README with data sources, access conditions (including for confidential data), the software and versions used, expected run time, and a table mapping every table and figure to the program that produces it;
2. raw data (or instructions to obtain them) kept separate from derived data, with the cleaning code;
3. a master script that runs everything in order with a single command;
4. a fixed environment: requirements.txt or environment.yml (Python), renv.lock (R), the Stata version plus ado-files;
5. fixed seeds for anything random (bootstrap, simulation, sample splits);
6. relative paths only, and no manual steps.

Test it by running the package from a clean machine or container before submission. Version control with Git, with tagged releases for submitted versions, makes the history auditable.

import numpy as np
from statsmodels.stats.multitest import multipletests

ss = np.random.SeedSequence(2024)
rngs = [np.random.default_rng(s) for s in ss.spawn(4)]   # independent streams, e.g. one per design

def coverage(rng, n, R=4000):
    hits = 0
    for _ in range(R):
        x = rng.lognormal(size=n)                        # skewed data: is the t-interval OK?
        m, s = x.mean(), x.std(ddof=1)
        hits += abs(m - np.exp(0.5)) <= 1.96 * s / np.sqrt(n)
    p = hits / R
    return p, np.sqrt(p * (1 - p) / R)

for rng, n in zip(rngs, (10, 30, 100, 1000)):
    p, se = coverage(rng, n)
    print(f"n={n:5d}: coverage {p:.3f} (MC s.e. {se:.3f})")

pvals = np.array([0.001, 0.008, 0.012, 0.030, 0.040, 0.300])
for method in ("bonferroni", "holm", "fdr_bh"):
    print(method, multipletests(pvals, alpha=0.05, method=method)[0])

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
df = wd.data("wagepan")
df["occupation"] = df[[f"occ{k}" for k in range(1, 10)]].to_numpy().argmax(1) + 1
f = "lwage ~ union + educ + exper + expersq + C(year)"
# TODO: restricted model without union; y* = fitted_r + v_g * resid_r; refit unrestricted; cluster t each draw

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
df = wd.data("wagepan").reset_index(drop=True)
df["occupation"] = df[[f"occ{k}" for k in range(1, 10)]].to_numpy().argmax(1) + 1
f = "lwage ~ union + educ + exper + expersq + C(year)"
g = df["occupation"].to_numpy(); G = len(np.unique(g))
fit = lambda d: smf.ols(f, data=d).fit(cov_type="cluster", cov_kwds={"groups": g, "use_correction": True})
base = fit(df); t0 = base.tvalues["union"]
from scipy import stats
print(f"{G} clusters: union = {base.params['union']:.4f}, cluster t = {t0:.2f}, p (normal) = {2*stats.norm.sf(abs(t0)):.4f}, p (t_{G-1}) = {2*stats.t.sf(abs(t0), G-1):.4f}")

r = smf.ols(f.replace("union + ", ""), data=df).fit()
rng = np.random.default_rng(1); ts = []
groups = np.unique(g)
for _ in range(999):
    v = dict(zip(groups, rng.choice([-1.0, 1.0], size=G)))
    d = df.copy(); d["lwage"] = r.fittedvalues + r.resid * np.vectorize(v.get)(g)
    ts.append(fit(d).tvalues["union"])
p_wcb = np.mean(np.abs(ts) >= abs(t0))
print(f"wild cluster bootstrap-t p-value: {p_wcb:.3f}   (with G = {G}, Rademacher weights give only 2^{G} = {2**G} distinct draws)")

9 clusters: union = 0.1768, cluster t = 5.37, p (normal) = 0.0000, p (t_8) = 0.0007
wild cluster bootstrap-t p-value: 0.012   (with G = 9, Rademacher weights give only 2^9 = 512 distinct draws)

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
from statsmodels.stats.multitest import multipletests
d = wd.data("jtrain2")
d["pos78"] = (d.re78 > 0).astype(int); d["gt5"] = (d.re78 > 5).astype(int); d["gt10"] = (d.re78 > 10).astype(int); d["chg"] = d.re78 - d.re75
# TODO: loop over outcomes, collect robust p-values, apply multipletests

import numpy as np, pandas as pd, wooldridge as wd, statsmodels.formula.api as smf
from statsmodels.stats.multitest import multipletests
d = wd.data("jtrain2")
d["pos78"] = (d.re78 > 0).astype(int); d["gt5"] = (d.re78 > 5).astype(int)
d["gt10"] = (d.re78 > 10).astype(int); d["chg"] = d.re78 - d.re75
outs = ["re78", "unem78", "pos78", "gt5", "gt10", "chg"]
rows = []
for y in outs:
    r = smf.ols(f"{y} ~ train", data=d).fit(cov_type="HC1")
    rows.append((y, r.params["train"], r.pvalues["train"]))
tab = pd.DataFrame(rows, columns=["outcome", "effect", "p"]).set_index("outcome")
for meth, lab in (("bonferroni", "Bonferroni"), ("holm", "Holm"), ("fdr_bh", "BH (FDR)")):
    tab[lab] = multipletests(tab["p"], method=meth)[1]
print(tab.round(4).to_string())

c = d[d.train == 0]; z = []
for y in outs:
    s = (d[y] - c[y].mean()) / c[y].std()
    z.append(-s if y == "unem78" else s)
d["index"] = np.mean(z, axis=0)
ri = smf.ols("index ~ train", data=d).fit(cov_type="HC1")
print(f"\nsummary index: effect {ri.params['train']:.3f} control-SD units (robust se {ri.bse['train']:.3f}, p = {ri.pvalues['train']:.4f})")

effect       p  Bonferroni    Holm  BH (FDR)
outcome                                              
re78     1.7943  0.0075      0.0449  0.0449    0.0216
unem78  -0.1106  0.0108      0.0648  0.0540    0.0216
pos78    0.1106  0.0108      0.0648  0.0540    0.0216
gt5      0.0771  0.1042      0.6252  0.2084    0.1158
gt10     0.0601  0.1158      0.6951  0.2084    0.1158
chg      1.5292  0.0324      0.1947  0.0973    0.0487

summary index: effect 0.227 control-SD units (robust se 0.085, p = 0.0074)