nb8.1 — The Specification-Curve Generator

Week 8, Chapter 8.1 companion. Maya runs her pre-registered confirmatory result once, then makes the entire garden of forking paths visible.

In nb7.5 Maya did the hardest thing in the course: she filed a pre-analysis plan, flipped the PAP_FILED gate, and ran her one pre-registered specification — exactly as written, reported whatever it showed. That single number, \widehat{\text{ATT}} \approx -1.89 pp, is her confirmatory result. It is the number her paper defends.

But a clean confirmatory run is necessary, not sufficient. The one regression she ran is one branch of a tree she could defensibly have climbed differently: cluster by county instead of state, add state-specific trends, drop the COVID year, define the outcome a touch differently. A referee is entitled to ask: if you had walked those forks, would the -1.89 still be there, or would it dissolve?

This notebook answers that question with the instrument from Simonsohn, Simmons & Nelson (2020) — the specification curve. We will:

Generate the same synthetic fair-lending panel as nb7.5, with a known planted true ATT, so we can grade every estimate against the truth.
Run the pre-registered primary specification once — the privileged, labeled point in the multiverse.
Enumerate every defensible analytic choice (controls, fixed effects, sample window, clustering, outcome definition) and take the cross-product into a grid of 4\times2\times3\times3\times2 = 144 specifications. Run them all.
Build the canonical two-panel spec curve: sorted estimates with CIs on top, a choice-indicator matrix below — so you can read downward from any region of the curve to the choices that produced it.
Summarize with the §8.1.3 checklist: share significant, sign stability, median-vs-primary — and read whether the result is robust.
Contrast with a fragile multiverse (Sam’s momentum anomaly) where the sign flips across specs — so you see both faces of the method.

Simonsohn, U., Simmons, J. P., & Nelson, L. D. (2020). Specification Curve Analysis. Nature Human Behaviour, 4, 1208–1214.

A boundary, restated from the chapter: the spec curve tests analytic-choice multiplicity (“did my choices drive the result?”). It is silent on identification (“did my design identify the effect?”) — that is nb8.2’s job, with placebos and Oster’s \delta. A robust spec curve on an invalid design is a confident wrong answer.

0. Setup — one pinned generator, headless plotting

Same conventions as every notebook this week: a single np.random.default_rng(42) drives the whole notebook so results are reproducible, and matplotlib.use("Agg") renders to buffers (no screen needed).

import matplotlib
matplotlib.use("Agg")  # headless: render to buffers, never to a screen

import itertools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import pyfixest as pf

rng = np.random.default_rng(42)  # one generator, pinned, drives the whole notebook

plt.rcParams["figure.dpi"] = 110
plt.rcParams["axes.grid"] = False

pd.set_option("display.width", 110)
pd.set_option("display.max_columns", 30)

print("numpy   ", np.__version__)
print("pandas  ", pd.__version__)
print("pyfixest", pf.__version__)

1. The analysis panel — Maya’s HMDA-style fair-lending DiD (with a known truth)

We regenerate the exact panel from nb7.4/nb7.5 so this notebook stands on the same ground. One row per county–year; some states adopt a fair-lending examination program in a staggered set of cohort years; never-adopters are clean controls. The outcome is the county-level minority–white mortgage-denial gap in percentage points, and the program narrows the gap — so the planted effect is negative.

The data-generating process, written forward in CONVENTIONS §3 notation, is the two-way fixed-effects model

\texttt{gap}_{it} = \alpha_i + \lambda_t + \tau\cdot\texttt{post\_treat}_{it} + \mathbf{x}_{it}'\boldsymbol{\gamma} + \varepsilon_{it},

with planted true ATT \tau = -1.80. Because we build the world, we know the right answer and can grade every one of the 144 specifications against it.

Two details make this notebook’s panel a superset of nb7.5’s, so the multiverse has something to vary:

We carry an extra applicant-composition control, loanamt (standardized loan amount), drawn in a separate pass after the gap so the gap draws match nb7.5 byte-for-byte — the primary spec will recover the same -1.89. loanamt is an ancillary control: a reasonable researcher might include it or not.
We carry a second, equally defensible outcome operationalization, gap_alt — the denial gap measured with a slightly different composition adjustment. Same percentage-point scale, defensibly different number. This is the “outcome definition” fork from §8.1.2.

# --- Planted truth and panel dimensions (identical to nb7.5) -----------------
TAU_TRUE   = -1.80          # programs NARROW the denial gap by 1.8 pp (the effect we must recover)
YEARS      = list(range(2014, 2024))
N_STATES   = 30             # 30 states -> 30 clusters (treatment varies at the STATE level)
COUNTIES_PER_STATE = 8      # 240 counties, ~2,400 county-year rows
GAMMA_INCOME = -0.9
GAMMA_LTV    =  1.2
SIGMA        =  0.8

# Staggered adoption; never-adopters (NaN) are clean controls.
adopt_cohort = {}
for s in range(N_STATES):
    if   s < 8:   adopt_cohort[s] = 2018
    elif s < 15:  adopt_cohort[s] = 2020
    else:         adopt_cohort[s] = np.nan

state_level = {s: rng.normal(0.0, 1.0) for s in range(N_STATES)}          # state baseline shifts
year_shock  = {t: 0.10 * (t - 2014) + rng.normal(0.0, 0.3) for t in YEARS}  # national year shocks

rows = []
for s in range(N_STATES):
    G = adopt_cohort[s]
    for c in range(COUNTIES_PER_STATE):
        county_id = s * COUNTIES_PER_STATE + c
        alpha_i = state_level[s] + rng.normal(0.0, 0.7)   # county permanent level
        for t in YEARS:
            income = rng.normal(0.0, 1.0)                 # standardized log applicant income
            ltv    = rng.normal(0.0, 1.0)                 # standardized loan-to-value
            post_treat = 1 if (not np.isnan(G) and t >= G) else 0
            # NOTE: same draw order as nb7.5 so the gap matches byte-for-byte.
            gap = (alpha_i + year_shock[t]
                   + TAU_TRUE * post_treat
                   + GAMMA_INCOME * income + GAMMA_LTV * ltv
                   + rng.normal(0.0, SIGMA))
            rows.append({"county": county_id, "state": s, "year": t, "cohort_G": G,
                         "post_treat": post_treat, "income": income, "ltv": ltv, "gap": gap})

panel = pd.DataFrame(rows)
# Ancillary control, drawn in a SEPARATE pass so the gap above is unchanged from nb7.5.
panel["loanamt"] = rng.normal(0.0, 1.0, len(panel))          # standardized loan amount
panel["adopter"] = panel["cohort_G"].notna().astype(int)
# Second defensible outcome operationalization: same pp scale, slightly different composition adjustment.
panel["gap_alt"] = panel["gap"] + 0.15 * panel["income"] - 0.10 * panel["ltv"]

print(f"county-year rows : {len(panel):,}")
print(f"states (clusters): {panel['state'].nunique()}   counties: {panel['county'].nunique()}")
print(f"adopting states  : {panel.loc[panel['adopter']==1,'state'].nunique()}"
      f"   never-adopters: {panel.loc[panel['adopter']==0,'state'].nunique()}")
print(f"PLANTED true ATT : {TAU_TRUE:+.2f} pp  (used only to grade the specs)")
panel.head()

2. The confirmatory run — execute the pre-registered primary spec ONCE

Before any multiverse, we honor the freeze. Maya’s pap-filed commit named, in CONVENTIONS §4 form, a single staggered difference-in-differences:

Outcome: gap (county–year minority–white denial gap, pp). Treatment: post_treat (= 1 once the county’s state has an active exam program). Controls: income, ltv. Fixed effects: county and year. Clustering: by state (the level treatment varies). Sample: the full balanced 2014–2023 panel. Identifying assumption: absent the programs, adopting and never-adopting states’ denial gaps would have moved in parallel.

This is one of the 144 specifications we are about to enumerate — but it is the privileged, labeled one. The confirmatory run comes first and is reported whatever it shows; the curve is context around it, never a way to choose it (§8.1.4, the “fishing-license” trap). We run it once and read off the number.

PRIMARY_SPEC = {
    "controls":     ["income", "ltv"],
    "fe":           "county+year",
    "sample":       "full",
    "cluster":      "state",
    "outcome":      "gap",
}

primary = pf.feols("gap ~ post_treat + income + ltv | county + year",
                   data=panel, vcov={"CRV1": "state"})

b_primary  = float(primary.coef()["post_treat"])
se_primary = float(primary.se()["post_treat"])
p_primary  = float(primary.pvalue()["post_treat"])
ci_primary = primary.confint().loc["post_treat"]
lo_primary, hi_primary = float(ci_primary.iloc[0]), float(ci_primary.iloc[1])

print("=" * 62)
print("  PRE-REGISTERED PRIMARY SPECIFICATION (confirmatory result)")
print("=" * 62)
print(f"  ATT on post_treat        : {b_primary:+.4f} pp")
print(f"  clustered SE (by state)  : {se_primary:.4f}")
print(f"  p-value (two-sided)      : {p_primary:.4g}")
print(f"  95% CI (clustered)       : [{lo_primary:+.4f}, {hi_primary:+.4f}]")
print("-" * 62)
print(f"  planted true ATT         : {TAU_TRUE:+.4f} pp")
print(f"  CI covers the truth?     : {lo_primary <= TAU_TRUE <= hi_primary}")
print("=" * 62)

# The honest-machinery check: the pre-registered spec recovers the planted effect.
assert lo_primary <= TAU_TRUE <= hi_primary, "primary CI failed to cover the planted ATT"
assert b_primary < 0, "expected a NEGATIVE effect (programs narrow the gap)"
assert abs(b_primary - (-1.89)) < 0.05, "primary should land near the chapter's -1.89"
print("\nConfirmatory result recovered ~-1.89 pp -- consistent with nb7.5 and Chapter 8.1.1.")

3. Enumerate the defensible analytic choices (the garden, written down)

Step 1 of Simonsohn–Simmons–Nelson is the act of writing the garden down: for each decision, list every option a competent researcher could defend — and only the defensible ones. Padding the multiverse with specs you would never defend in a seminar is p-hacking with extra steps (§8.1.4). The choices for Maya’s DiD partition into the CONVENTIONS §4 slots:

Controls (4): none / income only / income+ltv / all three (add loanamt).
Fixed effects (2): county+year baseline / county+year plus state-specific linear trends (a common parallel-trends hedge).
Sample window (3): full / drop 2020 (COVID) / drop the first-adopting cohort (thinly identified).
Clustering / inference (3): cluster by state (pre-registered) / cluster by county / HC1 heteroskedasticity-robust.
Outcome definition (2): the denial gap gap / the re-operationalized gap_alt.

The cross-product is 4 \times 2 \times 3 \times 3 \times 2 = 144 specifications. Notice that one of the 144 is the pre-registered primary — the confirmatory run lives inside the multiverse, as a labeled point.

Here the Week 1 arithmetic comes back to bite. With 144 specifications, even if the true effect were exactly zero, the chance that at least one lands at p<0.05 by luck is overwhelming (1-(1-0.05)^{144}\approx 1). The garden of forking paths is multiple testing wearing the disguise of “I only ran one regression.” The spec curve does not eliminate that problem — it shows you the whole distribution so a single starred specification can never masquerade as the only one you considered.

# Each slot maps a human-readable LABEL -> the machinery that implements that choice.

CONTROL_SETS = {
    "ctrl:none":        [],
    "ctrl:income":      ["income"],
    "ctrl:income+ltv":  ["income", "ltv"],     # <- the pre-registered control set
    "ctrl:all3":        ["income", "ltv", "loanamt"],
}
FE_SETS = {
    "fe:county+year":          False,          # baseline (pre-registered): plain county + year FE
    "fe:county+year+trend":    True,           # add state-specific linear time trends
}
SAMPLE_RULES = {
    "samp:full":        lambda d: d,                                          # pre-registered
    "samp:drop2020":    lambda d: d[d["year"] != 2020].copy(),
    "samp:drop-cohort": lambda d: d[d["cohort_G"] != 2018].copy(),            # drop 1st-adopting cohort
}
CLUSTER_RULES = {
    "vcov:state":   {"CRV1": "state"},          # pre-registered (level treatment varies)
    "vcov:county":  {"CRV1": "county"},
    "vcov:hc1":     "hetero",
}
OUTCOME_DEFS = {
    "out:gap":      "gap",                       # pre-registered
    "out:gap_alt":  "gap_alt",
}

n_specs = (len(CONTROL_SETS) * len(FE_SETS) * len(SAMPLE_RULES)
           * len(CLUSTER_RULES) * len(OUTCOME_DEFS))
print(f"Defensible choices per slot: controls={len(CONTROL_SETS)}, FE={len(FE_SETS)}, "
      f"sample={len(SAMPLE_RULES)}, cluster={len(CLUSTER_RULES)}, outcome={len(OUTCOME_DEFS)}")
print(f"Cross-product = {len(CONTROL_SETS)} x {len(FE_SETS)} x {len(SAMPLE_RULES)} x "
      f"{len(CLUSTER_RULES)} x {len(OUTCOME_DEFS)} = {n_specs} specifications")
assert n_specs == 144, "expected 144 specifications in the multiverse"

# The pre-registered choice labels -- so we can MARK the primary inside the curve later.
PRIMARY_LABELS = {"ctrl:income+ltv", "fe:county+year", "samp:full", "vcov:state", "out:gap"}

4. Run all 144 specifications and collect the choice-record for each

Step 2 is mechanical and embarrassingly fast: estimate the ATT on post_treat under every combination, and store the estimate, its 95% CI, its p-value, and which choice was made in each slot. That choice-record is the load-bearing part — it is what lets the bottom panel of the curve become a diagnostic.

A note on the state-trend fixed-effects option: a “state-specific linear time trend” is a separate slope per state on calendar time. We implement it by adding i(state, year_c) interaction terms (state dummies interacted with centered year) to the regressor list — pyfixest’s i() syntax builds exactly those. Everything else flows straight from the dictionaries above; nothing is hand-edited per spec.

import warnings

def run_spec(controls_lbl, fe_lbl, sample_lbl, cluster_lbl, outcome_lbl, data):
    """Estimate the ATT for one cell of the multiverse; return a tidy record."""
    ctrls   = CONTROL_SETS[controls_lbl]
    add_tr  = FE_SETS[fe_lbl]
    samp_fn = SAMPLE_RULES[sample_lbl]
    vcov    = CLUSTER_RULES[cluster_lbl]
    ycol    = OUTCOME_DEFS[outcome_lbl]

    d = samp_fn(data)
    rhs = ["post_treat"] + list(ctrls)
    if add_tr:
        d = d.copy()
        d["year_c"] = d["year"] - 2014                 # centered calendar time
        rhs = rhs + ["i(state, year_c)"]               # state-specific linear trends

    formula = f"{ycol} ~ " + " + ".join(rhs) + " | county + year"

    # With county + year FE plus a full set of state trends, ONE state-trend column is
    # mechanically collinear and pyfixest drops it -- expected and harmless. We silence
    # pyfixest's UserWarnings inside the loop to keep the 144-spec output readable.
    with warnings.catch_warnings():
        warnings.filterwarnings("ignore", category=UserWarning)
        m = pf.feols(formula, data=d, vcov=vcov)
    beta = float(m.coef()["post_treat"])
    se   = float(m.se()["post_treat"])
    ci   = m.confint().loc["post_treat"]
    return {
        "controls": controls_lbl, "fe": fe_lbl, "sample": sample_lbl,
        "cluster": cluster_lbl, "outcome": outcome_lbl,
        "beta": beta, "se": se,
        "lo": float(ci.iloc[0]), "hi": float(ci.iloc[1]),
        "pval": float(m.pvalue()["post_treat"]),
    }

records = []
for cl, fl, sl, vl, ol in itertools.product(
        CONTROL_SETS, FE_SETS, SAMPLE_RULES, CLUSTER_RULES, OUTCOME_DEFS):
    records.append(run_spec(cl, fl, sl, vl, ol, panel))

specs = pd.DataFrame(records)
specs["sig"] = specs["pval"] < 0.05
# Flag the single pre-registered specification inside the multiverse.
specs["is_primary"] = (
    (specs["controls"] == "ctrl:income+ltv") & (specs["fe"] == "fe:county+year")
    & (specs["sample"] == "samp:full") & (specs["cluster"] == "vcov:state")
    & (specs["outcome"] == "out:gap")
)

print(f"ran {len(specs)} specifications")
print(f"primary row present in multiverse? {specs['is_primary'].sum() == 1}")
assert specs["is_primary"].sum() == 1, "exactly one row must be the pre-registered primary"
# Sanity: the flagged primary row reproduces the confirmatory coefficient.
assert abs(float(specs.loc[specs['is_primary'], 'beta'].iloc[0]) - b_primary) < 1e-9
specs.head()

5. Summarize the multiverse — the §8.1.3 reading checklist

Before drawing the picture, read the numbers. The four-point checklist from §8.1.3:

Where does the sign live? If estimates span both signs, the direction is choice-dependent — the most serious fragility. If they share a sign, direction is robust even if magnitude wanders.
Where does the pre-registered spec sit? Near the center (representative) or out in a tail?
What fraction is significant, and where do the exceptions cluster?
Which choices drive the spread? We approximate this by the mean estimate within each choice level — the slot whose levels most separate tells you the fork your result is most sensitive to.

share_sig  = specs["sig"].mean()
share_neg  = (specs["beta"] < 0).mean()
share_pos  = (specs["beta"] > 0).mean()
median_est = specs["beta"].median()

print("MULTIVERSE SUMMARY (Maya's HMDA DiD)")
print("-" * 50)
print(f"  specifications              : {len(specs)}")
print(f"  median estimate             : {median_est:+.3f} pp")
print(f"  pre-registered primary      : {b_primary:+.3f} pp")
print(f"  estimate range              : [{specs['beta'].min():+.3f}, {specs['beta'].max():+.3f}]")
print(f"  share significant (p<.05)   : {share_sig:.1%}")
print(f"  share negative (correct sign): {share_neg:.1%}")
print(f"  share positive (wrong sign)  : {share_pos:.1%}")
print("-" * 50)

# Where does the primary sit in the sorted distribution? (percentile rank)
pct = (specs["beta"] < b_primary).mean()
print(f"  primary's percentile in the sorted estimates: {pct:.0%} "
      f"({'center-ish' if 0.25 <= pct <= 0.75 else 'tail'})")

print("\nMean estimate within each choice level (which fork moves the number?):")
for slot in ["controls", "fe", "sample", "cluster", "outcome"]:
    g = specs.groupby(slot)["beta"].mean().round(3)
    spread = g.max() - g.min()
    print(f"  [{slot:9s}] spread={spread:.3f} :: " +
          ", ".join(f"{k.split(':')[1]}={v:+.2f}" for k, v in g.items()))

6. The canonical specification curve (two stacked panels)

Step 3 — the figure that became the signature of the method. It has two panels sharing a horizontal axis that runs over all specifications sorted by the magnitude of the estimate (most negative on the left).

Top panel: the sorted point estimates as a curve, each with its 95%-CI whisker. Markers are filled when the spec is significant at 5%, hollow when not. The pre-registered primary is marked in a distinct color. Reference lines at zero and at the planted truth.
Bottom panel: a grid — one row per analytic choice, one column per specification, aligned to the top panel’s sort order. A tick at (\text{row}, \text{column}) means that spec used that choice. Read downward from any region of the curve to the choices that produced it.

We factor the drawing into a reusable function so we can apply it to both the robust and the fragile multiverses with one call.

def specification_curve(specs_df, choice_slots, truth=None, primary_mask=None,
                        title="", xlabel="estimate", figsize=(13, 9)):
    """Draw the canonical two-panel specification curve.

    specs_df      : DataFrame with columns beta, lo, hi, sig (bool) + the slot columns.
    choice_slots  : list of column names whose category-levels become the bottom-panel rows.
    truth         : optional float to draw as a reference line (the planted effect).
    primary_mask  : optional boolean Series marking the pre-registered specification.
    """
    df = specs_df.sort_values("beta").reset_index(drop=True).copy()
    x = np.arange(len(df))

    # Build the ordered list of (slot, level) rows for the bottom panel.
    row_labels, row_keys = [], []
    for slot in choice_slots:
        for lvl in sorted(df[slot].unique()):
            row_labels.append(lvl)
            row_keys.append((slot, lvl))

    fig, (ax_top, ax_bot) = plt.subplots(
        2, 1, figsize=figsize, sharex=True,
        gridspec_kw={"height_ratios": [2.4, max(2.0, 0.32 * len(row_keys))], "hspace": 0.05})

    # ---- top panel: sorted estimates + CIs -------------------------------
    sig = df["sig"].to_numpy()
    yerr = np.vstack([df["beta"] - df["lo"], df["hi"] - df["beta"]])
    # significant (filled) vs non-significant (hollow)
    ax_top.errorbar(x[sig], df["beta"][sig], yerr=yerr[:, sig], fmt="o", ms=4,
                    color="#3b6ea5", ecolor="#9bb8d6", elinewidth=0.8, capsize=0,
                    label="significant (p<.05)")
    ax_top.errorbar(x[~sig], df["beta"][~sig], yerr=yerr[:, ~sig], fmt="o", ms=4,
                    mfc="white", mec="#b5651d", color="#b5651d", ecolor="#e0c0a0",
                    elinewidth=0.8, capsize=0, label="not significant")
    ax_top.axhline(0.0, color="grey", ls=":", lw=1)
    if truth is not None:
        ax_top.axhline(truth, color="green", ls="--", lw=1.3, label=f"planted truth ({truth:+.2f})")
    # mark the primary by matching its beta position after the sort
    if primary_mask is not None and primary_mask.any():
        beta_primary = specs_df.loc[primary_mask, "beta"].iloc[0]
        # position(s) in sorted frame whose beta equals the primary's
        pos = np.where(np.isclose(df["beta"].to_numpy(), beta_primary))[0]
        if len(pos):
            ax_top.scatter(pos, [beta_primary] * len(pos), s=120, marker="D",
                           facecolor="none", edgecolor="red", linewidths=1.8, zorder=5,
                           label="pre-registered primary")
    ax_top.set_ylabel(xlabel)
    ax_top.set_title(title, fontsize=12)
    ax_top.legend(loc="upper left", fontsize=8, framealpha=0.9)
    ax_top.grid(axis="y", alpha=0.25)

    # ---- bottom panel: choice-indicator matrix ---------------------------
    M = np.zeros((len(row_keys), len(df)))
    for r, (slot, lvl) in enumerate(row_keys):
        M[r, :] = (df[slot].to_numpy() == lvl).astype(float)
    # draw ticks
    for r in range(len(row_keys)):
        cols = np.where(M[r] > 0)[0]
        ax_bot.scatter(cols, [r] * len(cols), s=6, marker="s", color="#333333")
    # faint slot separators
    boundary = 0
    for slot in choice_slots:
        n_lvl = df[slot].nunique()
        boundary += n_lvl
        ax_bot.axhline(boundary - 0.5, color="#cccccc", lw=0.6)
    ax_bot.set_yticks(range(len(row_labels)))
    ax_bot.set_yticklabels(row_labels, fontsize=8)
    ax_bot.set_ylim(-0.5, len(row_labels) - 0.5)
    ax_bot.invert_yaxis()
    ax_bot.set_xlabel("specifications, sorted by estimate (left = most negative)")
    ax_bot.set_xlim(-1, len(df))
    ax_bot.grid(False)
    fig.subplots_adjust(left=0.18, right=0.98, top=0.94, bottom=0.07)  # no tight_layout (shared-x panels)
    return fig

Now draw Maya’s curve. Watch the shape: a near-flat band sitting entirely below zero, almost every marker filled, the red diamond (her pre-registered primary) sitting near the center. That shape is the verdict.

fig = specification_curve(
    specs,
    choice_slots=["controls", "fe", "sample", "cluster", "outcome"],
    truth=TAU_TRUE,
    primary_mask=specs["is_primary"],
    title="Specification curve -- Maya's fair-lending DiD (ROBUST: 144 defensible specs)",
    xlabel="ATT on minority-white denial gap (pp)",
)
fig.savefig("nb81_speccurve_robust.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("saved figure: nb81_speccurve_robust.png")

7. Read the curve — is the result robust?

Apply the checklist to what we just computed and write the one honest paragraph the chapter asks for.

sig_n = int(specs["sig"].sum())
print("READING MAYA'S SPECIFICATION CURVE")
print("=" * 60)
print(f"1. Sign:   {share_neg:.0%} of specs are negative, {share_pos:.0%} positive "
      f"-> direction is {'ROBUST (never flips)' if share_neg==1.0 else 'choice-dependent'}.")
print(f"2. Primary: {b_primary:+.2f} pp sits at the {pct:.0%} percentile of the sorted "
      f"estimates -> {'central / representative' if 0.25<=pct<=0.75 else 'in a tail'}.")
print(f"3. Significance: {sig_n}/{len(specs)} ({share_sig:.0%}) reach p<.05.")
print(f"4. Median {median_est:+.2f} pp vs primary {b_primary:+.2f} pp -> "
      f"{'close' if abs(median_est-b_primary)<0.3 else 'apart'}.")
print("=" * 60)

robust = (share_neg == 1.0) and (share_sig >= 0.90)
print(f"\nVERDICT: {'ROBUST' if robust else 'FRAGILE'} across the multiverse.")
assert robust, "Maya's multiverse should read as robust"

print("""
The honest paragraph for the paper:

  Across all 144 defensible specifications, the estimated effect of fair-lending
  examinations on the minority-white denial gap is negative, with a median of
  {med:+.1f} pp and a range of [{lo:+.1f}, {hi:+.1f}] pp. The pre-registered
  primary specification ({pri:+.1f} pp) sits near the center of this distribution.
  The effect is statistically significant in {sig}/{n} specifications ({pct:.0%});
  the sign never flips. The result is a robust feature of the data, not an artifact
  of analytic choice.
""".format(med=median_est, lo=specs['beta'].min(), hi=specs['beta'].max(),
           pri=b_primary, sig=sig_n, n=len(specs), pct=share_sig))

8. The cautionary contrast — a FRAGILE multiverse (Sam’s momentum anomaly)

Reveal-the-trick demands you see the picture that should terrify you. Sam asks whether past-12-month winners outperform losers next month — a portfolio sort with a long menu of forks: formation window, holding period, value- vs equal-weighting, NYSE vs all-stock breakpoints, including or excluding microcaps, and the sample period. Those forks multiply past 200 specifications.

We build a synthetic momentum multiverse with a planted structure that mirrors the real anomaly’s fragility: microcaps carry a spurious positive momentum tilt, while large/liquid stocks show mild reversal. So the sign of the long-short premium depends on how you weight and which universe you keep — exactly the fork the chapter warns about. There is no single “true ATT” here; the point is that the headline is an artifact of two choices.

import math

# ---- a synthetic stock universe -------------------------------------------
N_STOCKS = 600
size      = rng.normal(0.0, 1.0, N_STOCKS)        # log market cap (standardized)
is_micro  = size < -0.5                            # the smallest ~31% are "microcaps"
stocks    = pd.DataFrame({"stock": np.arange(N_STOCKS), "size": size, "is_micro": is_micro})

def momentum_longshort(formation, holding, weighting, universe, period, breakpoints, gen):
    """Return (estimate %/month, SE %/month, p-value) for one momentum specification."""
    df = stocks if universe == "uni:all" else stocks[~stocks["is_micro"]]
    # planted per-stock long-short contribution: microcaps +, large caps - (reversal)
    contrib = np.where(df["is_micro"].to_numpy(), 0.012, -0.004)
    contrib = contrib * {"form:12m": 1.0, "form:12m-skip": 0.9, "form:6m": 0.7}[formation]
    contrib = contrib * {"hold:1m": 1.0, "hold:3m": 0.85, "hold:6m": 0.7}[holding]
    contrib = contrib * {"per:full": 1.0, "per:post2000": 0.95, "per:ex0809": 1.05}[period]
    contrib = contrib * {"bp:nyse": 0.95, "bp:allstock": 1.05}[breakpoints]
    if weighting == "wt:vw":
        w = np.exp(df["size"].to_numpy()); w = w / w.sum()    # value-weighting buries microcaps
        ls = float(np.sum(w * contrib))
    else:
        ls = float(np.mean(contrib))                          # equal-weighting amplifies microcaps
    se  = 0.0020 + 0.0008 * gen.random()
    est = ls + gen.normal(0.0, se)
    t   = est / se
    pval = 2 * (1 - 0.5 * (1 + math.erf(abs(t) / math.sqrt(2))))
    return est * 100.0, se * 100.0, pval                       # to %/month

FORMATION = ["form:12m", "form:12m-skip", "form:6m"]
HOLDING   = ["hold:1m", "hold:3m", "hold:6m"]
WEIGHT    = ["wt:ew", "wt:vw"]
UNIVERSE  = ["uni:all", "uni:no-micro"]
PERIOD    = ["per:full", "per:post2000", "per:ex0809"]
BREAKPTS  = ["bp:nyse", "bp:allstock"]

mom_records = []
for f, h, w, u, p, bp in itertools.product(FORMATION, HOLDING, WEIGHT, UNIVERSE, PERIOD, BREAKPTS):
    est, se, pv = momentum_longshort(f, h, w, u, p, bp, rng)
    mom_records.append({"formation": f, "holding": h, "weighting": w, "universe": u,
                        "period": p, "breakpoints": bp,
                        "beta": est, "se": se, "lo": est - 1.96 * se, "hi": est + 1.96 * se,
                        "pval": pv})
mom = pd.DataFrame(mom_records)
mom["sig"] = mom["pval"] < 0.05

print(f"Sam's momentum multiverse: {len(mom)} specifications "
      f"({len(FORMATION)}x{len(HOLDING)}x{len(WEIGHT)}x{len(UNIVERSE)}x{len(PERIOD)}x{len(BREAKPTS)})")
print(f"  median estimate : {mom['beta'].median():+.3f} %/month")
print(f"  estimate range  : [{mom['beta'].min():+.3f}, {mom['beta'].max():+.3f}] %/month")
print(f"  share positive  : {(mom['beta']>0).mean():.0%}   share negative: {(mom['beta']<0).mean():.0%}")
print(f"  share significant: {mom['sig'].mean():.0%}")

Draw Sam’s curve with the same machinery. The shape is the opposite of Maya’s: it sweeps from significantly negative on the left, through a fat band of near-zero estimates, up to significantly positive on the right. The sign is not a feature of the world — it is a feature of the forks.

fig = specification_curve(
    mom,
    choice_slots=["formation", "holding", "weighting", "universe", "period", "breakpoints"],
    truth=None,
    primary_mask=None,
    title="Specification curve -- Sam's momentum anomaly (FRAGILE: sign flips across specs)",
    xlabel="long-short return (%/month)",
    figsize=(13, 11),
)
fig.savefig("nb81_speccurve_fragile.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("saved figure: nb81_speccurve_fragile.png")

Now the forensic read of the bottom panel: take every significantly positive specification and ask what they share. If they all use one or two specific choices, the bottom panel has found the fork that manufactures the headline.

sig_pos = mom[(mom["sig"]) & (mom["beta"] > 0)]
sig_neg = mom[(mom["sig"]) & (mom["beta"] < 0)]

print("READING SAM'S SPECIFICATION CURVE")
print("=" * 60)
print(f"1. Sign:   {(mom['beta']>0).mean():.0%} positive, {(mom['beta']<0).mean():.0%} negative "
      f"-> direction is CHOICE-DEPENDENT (the most serious fragility).")
print(f"2. Significant-and-positive specs : {len(sig_pos)}")
print(f"   Significant-and-negative specs : {len(sig_neg)}")
print("=" * 60)

if len(sig_pos):
    print("\nWhat do ALL significantly-POSITIVE specs share?")
    for slot in ["formation", "holding", "weighting", "universe", "period", "breakpoints"]:
        vals = set(sig_pos[slot])
        flag = "  <-- the fork that drives the headline" if len(vals) == 1 else ""
        print(f"  {slot:12s}: {sorted(vals)}{flag}")

all_ew   = (sig_pos["weighting"] == "wt:ew").all() if len(sig_pos) else False
all_micro = (sig_pos["universe"] == "uni:all").all() if len(sig_pos) else False
print(f"\nEvery significantly-positive spec uses equal-weighting?      {all_ew}")
print(f"Every significantly-positive spec includes microcaps?       {all_micro}")
assert all_ew and all_micro, "fragile case should pin the headline to ew + microcap-inclusive forks"

print("""
Sam's honest paragraph (the mirror image of Maya's):

  The sign of the momentum premium is NOT robust to the weighting and universe
  choices. It is positive and significant only under equal-weighting with microcaps
  included, and is statistically indistinguishable from zero -- or negative -- under
  the value-weighted, liquid-universe specifications that a tradable strategy would
  actually face. The headline "momentum earns positive returns" is a single branch:
  the equal-weighted-microcap branch. The curve drags it into the light.
""")

9. Connect to Week 1 — the garden of forking paths is multiple testing

The reason this whole exercise matters traces straight back to Chapter 1.5. If you test k independent true nulls at \alpha=0.05, the probability that at least one crosses the line by luck is 1-(1-0.05)^k. With k=20 that is already \approx 64\%; with k=144 it is essentially certain. The 144 specifications are not “one regression” — they measure how many regressions you could defensibly have run, which is the real size of your multiple-testing problem.

The p-hacker exploits this by running many and reporting the one that looks best. The specification-curve analyst runs many and reports all of them at once. Same multiplicity; opposite ethic. Below, the Week-1 arithmetic made literal.

def family_wise_at_least_one(k, alpha=0.05):
    """P(at least one false positive) among k independent true nulls, Chapter 1.5."""
    return 1.0 - (1.0 - alpha) ** k

for k in [1, 5, 20, 50, 144, 216]:
    print(f"  k={k:4d} independent tests -> P(>=1 spurious 'winner') = {family_wise_at_least_one(k):.4f}")

print(f"\nMaya's multiverse has {len(specs)} specs; Sam's has {len(mom)}.")
print("Reporting one starred spec out of that many, while hiding the rest, IS p-hacking.")
print("The spec curve's defense is exhaustive DISCLOSURE: show the whole forest, not the best tree.")

10. When the curve lies — four caveats you must state in the paper

A picture this persuasive is dangerous because it is persuasive (§8.1.4). State which caveats apply:

The curve only contains the forks you enumerated. It is exhaustive over choices you listed and silent on every one you did not. Maya’s curve holds the identification design fixed across all 144 specs — it can look gloriously robust while resting on one unexamined fork. Enumerate honestly and broadly, and say what you held fixed.
Not all specs are equally defensible, and the curve pretends they are. Clustering at the level treatment varies (state) is more defensible than by county. Padding the multiverse with weak-but-arguable specs to drag the median is p-hacking with a fancier chart.
The curve describes analytic-choice fragility, not identification. All 144 of Maya’s specs share the same parallel-trends assumption. If it fails, the curve robustly reports a beautiful artifact. That is nb8.2’s job (placebos, Oster’s \delta).
The curve can become a fishing license. Running the multiverse and then reporting “the median spec” is outcome-driven selection in a lab coat. The PAP fixes the primary in advance; the curve is reported around it, never used to choose it. That is why §2 (confirmatory) came before §3 (multiverse) in this notebook — the order is not negotiable.

Your Turn

Move the truth. Set TAU_TRUE = 0.0 and re-run §1–§7. The multiverse should now scatter around zero: with a true null, the share-significant collapses toward the false-positive rate, and the spec curve correctly reads “no robust effect.” Notice you do not get to fish for the one starred spec.
Find the fork that flips a sign. In Maya’s robust multiverse, no single choice flips the sign — that is the strong result. In Sam’s fragile one, two choices (weighting + universe) do. Modify the momentum DGP so that the holding period also matters (e.g., make the microcap tilt vanish at hold:6m). Re-run §8 and re-read the bottom panel: which rows now best separate positive from negative?
Add a fork you forgot. Caveat (1) says the curve is only as complete as your enumeration. Add a sixth slot to Maya’s grid — say, winsorizing the outcome at the 1st/99th percentile (on/off) — re-run, and check whether the verdict changes. If it does, your original enumeration was incomplete; if it does not, you have made the robustness claim stronger.
(Extension) The permutation test on the whole curve. Simonsohn–Simmons–Nelson attach a permutation test: under the sharp null that the effect is zero everywhere, re-shuffle the treatment assignment, recompute all the specifications on each shuffle, and ask whether the observed curve (e.g., the median estimate, or the count of significant-correctly-signed specs) is more extreme than the null curves produce. Implement it for Maya’s multiverse on a coarsened grid (say, fix the outcome and clustering to keep it fast) with ~200 shuffles, and report the whole-forest p-value. This tests the multiverse as a unit — the right unit once you admit “which specification” was itself a choice.