nb7.4 — Building the Analysis Dataset: Merge, Lag, and Validate

In Chapter 7.2, Sam pulled two raw files off WRDS: monthly stock returns from CRSP and annual accounting numbers from Compustat. Each file is clean on its own. But his question — do high book-to-market (“value”) firms earn higher returns? — lives in neither file: the returns are in CRSP, the book equity is in Compustat. To answer it he has to join them, and the moment he types merge, every quiet decision becomes a loud one with consequences.

The trick of this notebook, stated up front: a dataset is not data; it is a sequence of decisions, and the decisions are the dataset. We will build Sam’s CRSP–Compustat value panel end to end — link, lag, screen, winsorize — and at every step do the two things amateurs skip: count the rows and check the match rate.

Everything here runs offline on synthetic data we generate ourselves, seeded so it is identical every run. The numbers are fake; the machinery and the bugs it catches are exactly what you will use on real CRSP/Compustat. We never touch licensed data — and per the CONVENTIONS, real CRSP/Compustat extracts stay read-only on GMU infrastructure with a pinned snapshot date.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")  # non-interactive backend: render to figures, never pop a window
import matplotlib.pyplot as plt
import tempfile, os, logging

rng = np.random.default_rng(42)  # one seed -> identical synthetic data every run

print("pandas", pd.__version__, "| numpy", np.__version__)

# All notebook outputs (plots, the final panel) go to a temp dir, never the repo.
WORKDIR = tempfile.mkdtemp(prefix="nb74_")
print("scratch dir:", WORKDIR)

1. Generate three synthetic raw tables

Real CRSP and Compustat were built by different people, decades apart, and do not use the same identifier for a firm. CRSP tracks securities and keys them by PERMNO. Compustat tracks companies and keys them by GVKEY. A PERMNO is a security; a GVKEY is a company; one company can have several share classes (several PERMNOs). You cannot merge on a ticker — tickers are recycled nicknames, not identifiers.

We synthesize three tables that mimic this world:

  1. CRSP-like — monthly returns keyed by (permno, date), with a dlret delisting return for firms that die.
  2. Compustat-like — annual book equity keyed by (gvkey, datadate) where datadate is the fiscal-year-end.
  3. CCM-like link — the crosswalk mapping permno <-> gvkey, but time-varying: each link is valid only within a [linkdt, linkenddt] window. That date window is the whole game — it is why a crosswalk is not just a {permno: gvkey} dictionary.

We deliberately build in entries (IPOs), exits (delistings), one company with two share classes (two PERMNOs), heavy-tailed book equity, some negative book equity, and some missing returns — so every downstream lesson has something real to catch.

# ---- Universe of firms ---------------------------------------------------
N_FIRMS = 60
YEARS   = range(2015, 2021)          # 2015..2020 fiscal years
MONTHS  = pd.date_range("2015-01-31", "2021-12-31", freq="ME")  # month-end dates

gvkeys  = np.arange(1001, 1001 + N_FIRMS)            # company keys
permnos = np.arange(80001, 80001 + N_FIRMS)          # one base security per firm

# Each firm enters in a random year (IPO) and some exit early (delisting),
# so the per-year firm count CHURNS rather than only growing.
entry_year = rng.integers(2015, 2018, size=N_FIRMS)            # IPO year
exit_year  = entry_year + rng.integers(2, 8, size=N_FIRMS)     # death/last year
exit_year  = np.minimum(exit_year, 2022)
firms = pd.DataFrame({"gvkey": gvkeys, "permno": permnos,
                      "entry_year": entry_year, "exit_year": exit_year})
firms["delists"] = firms["exit_year"] <= 2021     # firms that die inside the window
print(firms.head())
print(f"\n{firms['delists'].sum()} of {N_FIRMS} firms delist before 2022 (these are the survivorship lesson)")

1a. The CRSP-like monthly returns table

We emit one row per (permno, month) for every month the firm is alive. A firm with book-to-market skill in real life would have returns correlated with its characteristics; here we just need a plausible spread. The firms that delist get a final, often catastrophic, delisting return (dlret) in their last month — CRSP retains delisted securities precisely so you can incorporate this exit return instead of silently dropping the firm.

We also blank out a small fraction of returns at random, so the missing-data probe later has something to find.

rows = []
for _, f in firms.iterrows():
    # months this firm is alive: from June of entry_year through Dec of exit_year
    start = pd.Timestamp(f.entry_year, 6, 30)
    end   = pd.Timestamp(min(f.exit_year, 2021), 12, 31)
    alive = MONTHS[(MONTHS >= start) & (MONTHS <= end)]
    if len(alive) == 0:
        continue
    # value-ish firms get a slightly higher mean monthly return (pure synthetic)
    mu = rng.normal(0.008, 0.004)
    ret = rng.normal(mu, 0.06, size=len(alive))
    dlret = np.full(len(alive), np.nan)
    if f.delists and f.exit_year <= 2021:
        dlret[-1] = rng.normal(-0.35, 0.10)   # catastrophic exit return in last month
    rows.append(pd.DataFrame({"permno": f.permno, "date": alive,
                              "ret": ret, "dlret": dlret}))

crsp_raw = pd.concat(rows, ignore_index=True)   # ONE concat, never df.append in a loop

# Blank ~3% of returns at random to exercise the missing-data probe later.
miss_idx = rng.choice(crsp_raw.index, size=int(0.03 * len(crsp_raw)), replace=False)
crsp_raw.loc[miss_idx, "ret"] = np.nan

# Add the second share class for ONE company (same gvkey, a NEW permno) ----
twin_gvkey = 1001                                  # firm 0 has two share classes
twin = crsp_raw[crsp_raw["permno"] == 80001].copy()
twin["permno"] = 90001                             # class-C permno
twin["ret"] = twin["ret"] + rng.normal(0, 0.005, size=len(twin))
crsp_raw = pd.concat([crsp_raw, twin], ignore_index=True)

crsp_raw["year"] = crsp_raw["date"].dt.year
print(f"crsp_raw: {len(crsp_raw):,} rows, {crsp_raw['permno'].nunique()} permnos")
crsp_raw.head()

1b. The Compustat-like fundamentals table

One row per (gvkey, fiscal-year). The key field is datadate — the fiscal-year-end date, not the date the number became public. A firm with a December fiscal-year-end does not file its 10-K on December 31; it files weeks or months later. We will respect that lag in Section 4.

Book equity is heavy-tailed by construction (most firms modest, a few enormous) and a handful of firm-years have negative book equity — a real Compustat phenomenon that makes book-to-market explode or flip sign, and a standard Fama–French screen drops them.

frames = []
for _, f in firms.iterrows():
    fy = [y for y in YEARS if f.entry_year <= y <= f.exit_year]
    if not fy:
        continue
    # heavy-tailed book equity: lognormal base scale per firm
    base = rng.lognormal(mean=5.0, sigma=1.2)
    be = base * (1 + rng.normal(0, 0.15, size=len(fy)))
    frames.append(pd.DataFrame({"gvkey": f.gvkey,
                                "datadate": [pd.Timestamp(y, 12, 31) for y in fy],
                                "book_equity": be}))
comp_raw = pd.concat(frames, ignore_index=True)

# Inject a few NEGATIVE book-equity firm-years (distressed firms) ----------
neg_idx = rng.choice(comp_raw.index, size=8, replace=False)
comp_raw.loc[neg_idx, "book_equity"] = -np.abs(rng.normal(20, 10, size=8))

# Inject a few EXTREME book-equity values (data-error / huge-firm tail) -----
big_idx = rng.choice(comp_raw.index, size=5, replace=False)
comp_raw.loc[big_idx, "book_equity"] = comp_raw.loc[big_idx, "book_equity"] * 60

comp_raw["fyear"] = comp_raw["datadate"].dt.year
print(f"comp_raw: {len(comp_raw):,} firm-years, {comp_raw['gvkey'].nunique()} gvkeys")
print(f"negative book-equity firm-years: {(comp_raw['book_equity'] < 0).sum()}")
comp_raw.head()

2. Merge discipline: cardinality is a claim you must verify

Every merge has a cardinality — a claim about how rows relate. One-to-one? Many-to-one? You think you know which, and you are frequently wrong because of a duplicate you did not know was there. A merge whose cardinality is wrong does not crash; it silently does the wrong thing — most destructively a many-to-many join that produces a Cartesian explosion and quietly reweights your regression.

The fix is one keyword. pandas.merge(..., validate=...) asserts the cardinality and raises MergeError if reality disagrees:

  • "one_to_one" / "1:1", "one_to_many", "many_to_one" / "m:1" — these check a uniqueness assumption.
  • "many_to_many" — checks nothing (the silent default).

Always pass validate=. It turns silent data corruption into a loud exception at the exact line that caused it. Below we link the primary CCM table onto CRSP returns asserting many_to_one (many monthly rows per security : one link).

# Assert: each CRSP (permno) row matches AT MOST ONE primary link row.
linked = crsp_raw.merge(
    ccm_primary[["permno", "gvkey", "linkdt", "linkenddt"]],
    on="permno",
    how="left",
    validate="many_to_one",   # many CRSP rows : one link -- ASSERTED, not hoped
)
print(f"linked OK with validate='many_to_one': {len(linked):,} rows")
linked.head(3)

2a. Demonstrate the explosion validate= catches

Now the cautionary tale. Suppose we had been sloppy and tried to link CRSP to the full CCM table (which contains the non-primary second-share-class row for gvkey 1001) and, worse, merged on gvkey against a Compustat table that has many rows per gvkey. That is a genuine many-to-many key, and an unguarded merge would explode it into a Cartesian product. With validate="one_to_one" asserted, pandas refuses.

We wrap it in try/except so the notebook keeps running — the point is to see the guard fire, not to crash.

# A deliberately WRONG merge: comp_raw has MANY rows per gvkey (one per fiscal
# year), and we ask pandas to assert one_to_one. Reality (many:1) disagrees,
# so validate= raises MergeError BEFORE any silent Cartesian explosion.
from pandas.errors import MergeError

# First, show what the unguarded (validate=None) merge would have produced:
exploded = linked.merge(comp_raw[["gvkey", "fyear", "book_equity"]],
                        on="gvkey", how="left")  # NO validate -> silent m:m
print(f"Unguarded merge exploded {len(linked):,} -> {len(exploded):,} rows "
      f"({len(exploded)/len(linked):.1f}x) -- every dup reweights the regression!\n")

# Now the guarded version catches it:
try:
    bad = linked.merge(comp_raw[["gvkey", "fyear", "book_equity"]],
                       on="gvkey", how="left",
                       validate="one_to_one")   # <-- the tripwire
    print("merge succeeded (unexpected)")
except MergeError as e:
    print("validate= CAUGHT the bad merge (as intended):")
    print("  MergeError:", e)

3. Survivorship bias: the firms that vanished are the ones that matter

Survivorship bias is the error of analyzing only the units that survived to the end of your sample, when survival is correlated with the outcome. If Sam downloads the list of firms trading today and pulls their history backward, every firm in that list survived — the bankruptcies and delistings (often the distressed, high-book-to-market firms whose prices collapsed) are gone. He keeps the value firms that recovered and discards the ones that died, inflating the “value premium” into an artifact.

The smell test is simple: count firms per year. A real universe churns — IPOs enter, delistings depart, the count wobbles. A survivorship-biased universe’s count only ever grows. Below we plot both: our honest synthetic universe (with exits) versus the same data after the naive “keep only firms alive in 2021” filter, and we show how the naive filter biases the mean return upward by discarding the catastrophic delisting returns.

# Honest universe: firms per year (should churn).
counts_all = crsp_raw.groupby("year")["permno"].nunique()

# Survivorship-biased universe: keep only permnos still alive in 2021.
survivors = crsp_raw.loc[crsp_raw["year"] == 2021, "permno"].unique()
crsp_surv = crsp_raw[crsp_raw["permno"].isin(survivors)].copy()
counts_surv = crsp_surv.groupby("year")["permno"].nunique()

# Effective return = ret compounded with delisting return in the exit month.
crsp_raw["ret_eff"]  = crsp_raw["ret"].add(crsp_raw["dlret"].fillna(0.0))
crsp_surv["ret_eff"] = crsp_surv["ret"].add(crsp_surv["dlret"].fillna(0.0))
mean_all  = crsp_raw["ret_eff"].mean()
mean_surv = crsp_surv["ret_eff"].mean()
print(f"Mean monthly return  -- full (with exits): {mean_all:+.4%}")
print(f"Mean monthly return  -- survivors only   : {mean_surv:+.4%}")
print(f"Survivorship bias (survivors - full)     : {mean_surv - mean_all:+.4%}  "
      f"(upward: we dropped the dying firms' crash returns)")
assert mean_surv > mean_all, "survivorship filter should bias the mean UP"

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(counts_all.index, counts_all.values, "o-", label="Full universe (churns)")
ax.plot(counts_surv.index, counts_surv.values, "s--",
        label="Survivors only (monotone -> biased)")
ax.set_xlabel("year"); ax.set_ylabel("# unique firms"); ax.legend()
ax.set_title("Survivorship smell test: does the firm count ever fall?")
fig.tight_layout(); fig.savefig(os.path.join(WORKDIR, "survivorship.png"), dpi=90)
plt.close(fig)
print("\nfull-universe per-year counts:\n", counts_all.to_string())

4. Look-ahead guard: the 6-month lag and merge_asof

Look-ahead bias is using information that was not yet knowable when the decision would have acted. Your code sees the whole dataset at once; an investor in March 2015 could not see June 2015’s data. The canonical source in finance is the reporting lag: Compustat stamps fundamentals with the fiscal-year-end (datadate), but the firm files its 10-K weeks-to-months later.

The conservative Fama–French convention is a 6-month lag: accounting data for a fiscal year ending in calendar year y is usable only from June 30 of year y+1. So fiscal-2014 book equity is matched to returns from July 2015 through June 2016. We compute an available_date and join with pd.merge_asof(direction="backward"), which structurally cannot attach a future record — it walks backward to the most recent available value. Then an assert makes the no-peeking guarantee a tripwire that runs on every build.

# Point-in-time availability date: 6-month lag (June 30 of the NEXT year).
fund = comp_raw[["gvkey", "datadate", "book_equity"]].copy()
fund["available_date"] = pd.to_datetime(
    (fund["datadate"].dt.year + 1).astype(str) + "-06-30"
)

# merge_asof requires BOTH sides sorted by the as-of key.
fund_s   = fund.sort_values("available_date").reset_index(drop=True)
linked_s = linked.dropna(subset=["gvkey"]).sort_values("date").reset_index(drop=True)
linked_s["gvkey"] = linked_s["gvkey"].astype(int)
fund_s["gvkey"]   = fund_s["gvkey"].astype(int)

panel = pd.merge_asof(
    linked_s, fund_s,
    left_on="date", right_on="available_date",
    by="gvkey",
    direction="backward",         # only look BACKWARD in time -- no peeking
)

# THE TRIPWIRE: a fundamental is never used before it was available.
used = panel.dropna(subset=["available_date"])
assert (used["date"] >= used["available_date"]).all(), \
    "LOOK-AHEAD: a return uses accounting data dated after the return!"
print(f"No-look-ahead assertion HOLDS over {len(used):,} matched rows.")

# Concrete check: a Dec-2015 fiscal number should only attach from Jul-2016 on.
ex = panel[(panel["gvkey"] == 1002) & panel["available_date"].notna()]
if len(ex):
    first = ex.sort_values("date").iloc[0]
    print(f"gvkey 1002: fiscal {first['datadate'].date()} -> "
          f"first usable return month {first['date'].date()} "
          f"(available_date {first['available_date'].date()})")

5. Winsorize the heavy-tailed variable, within the cross-section

Sam computes book-to-market (bm = book_equity / market_equity; we proxy market equity with a synthetic price level) and finds most values between 0.2 and 2, but a few are 400 and one is negative. Outliers dominate any mean and any OLS slope. The standard tool is winsorizing: instead of deleting extremes, cap them at a percentile (here 1% / 99%). The row stays in the sample — you do not lose N — but its leverage is bounded.

Three disciplines separate a defensible winsor from a fudge: (1) winsorize within the cross-section (by date), so “extreme” is judged against contemporaries, not a different decade; (2) state the cutoffs and apply them once — searching cutoffs until your t-stat clears 2 is p-hacking; (3) do not winsorize the variable that is the whole point (Priya’s catastrophe losses, Devon’s crypto tails are the signal). Here bm is a control-style ratio, so we cap it.

# Build book-to-market on the matched, screened panel.
panel = panel.query("book_equity > 0").copy()          # standard FF screen
# synthetic "market equity" so we have a denominator (offline proxy)
rng_me = np.random.default_rng(7)
panel["market_equity"] = rng_me.lognormal(5.0, 1.0, size=len(panel))
panel["bm"] = panel["book_equity"] / panel["market_equity"]

def winsorize(s: pd.Series, lower=0.01, upper=0.99) -> pd.Series:
    """Cap a series at its lower/upper quantiles. Returns a NEW series."""
    lo, hi = s.quantile([lower, upper])
    return s.clip(lower=lo, upper=hi)

# Winsorize WITHIN each cross-section (by date), not pooled across all years.
panel["bm_w"] = (panel.groupby("date")["bm"]
                      .transform(lambda s: winsorize(s, 0.01, 0.99)))

raw_max, w_max = panel["bm"].max(), panel["bm_w"].max()
print(f"bm        max={raw_max:8.2f}  mean={panel['bm'].mean():.3f}")
print(f"bm_w (1/99) max={w_max:8.2f}  mean={panel['bm_w'].mean():.3f}")
print(f"winsorizing pulled the max DOWN by {raw_max - w_max:.2f} "
      f"and moved the mean by {panel['bm_w'].mean() - panel['bm'].mean():+.4f}")
assert w_max <= raw_max, "winsorized max cannot exceed raw max"
assert panel["bm_w"].min() >= panel["bm"].min()

fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(panel["bm"], bins=40, alpha=0.5, label="raw bm")
ax.hist(panel["bm_w"], bins=40, alpha=0.5, label="winsorized bm_w (1/99)")
ax.set_xlabel("book-to-market"); ax.set_ylabel("count"); ax.legend()
ax.set_title("Winsorizing caps the tail, keeps the row")
fig.tight_layout(); fig.savefig(os.path.join(WORKDIR, "winsor.png"), dpi=90)
plt.close(fig)

6. Missing-data probe: why it is missing beats how much

The professional question is never “how much is missing?” but “why?” If returns are Missing Completely At Random (a server dropped records), dropping those rows loses precision but not unbiasedness. If missingness depends on the value itself (Missing Not At Random), dropping is selecting on the outcome — the survivorship error in a different hat. The rule is absolute: never silently drop. Measure it, probe whether it is informative with a groupby on a missingness flag, and log the cut.

panel["ret_missing"] = panel["ret"].isna()
pct_missing = panel["ret_missing"].mean()
print(f"returns missing: {pct_missing:.2%}")

# Probe: is missingness informative? Compare a stable observed characteristic --
# book equity (log scale to tame the heavy tail) -- across the missing vs.
# present groups. By construction we blanked returns COMPLETELY at random, so
# the two group means should be statistically indistinguishable (MCAR). On real
# HMDA data you would instead see denial rates DIFFER sharply across the flag,
# which would signal MNAR and warn you NOT to dropna() silently.
panel["log_be"] = np.log(panel["book_equity"])
probe = panel.groupby("ret_missing")["log_be"].agg(["count", "mean", "std"])
print("\nmissingness probe (groupby flag, log book equity):\n", probe.to_string())

m_miss, m_pres = probe.loc[True, "mean"], probe.loc[False, "mean"]
rel_gap = abs(m_miss - m_pres) / probe["std"].mean()   # gap in pooled-SD units
print(f"\n|mean(log_be | missing) - mean(log_be | present)| = {rel_gap:.2f} SD")
print("  small (<~0.3 SD) -> consistent with MCAR, listwise deletion is safe here;")
print("  a large gap would flag informative (MNAR) missingness -- do NOT drop silently.")

7. The build log: raw -> intermediate -> analysis

The final discipline is to make the decisions reproducible and auditable with row-count logging at every step. Each transformation logs how many rows entered and left; that log is your audit trail. When a referee asks “where did the other 4,000 firm-months go?”, you point at it. Below, a tiny step() helper logs (rows, unique firms) after each stage, and we rerun the whole pipeline through it so it reads like the autobiography of the dataset. We end by saving the clean analysis panel to the temp dir.

logging.basicConfig(level=logging.INFO, format="%(message)s")
log = logging.getLogger("build")
for h in list(log.handlers):           # idempotent if the cell is re-run
    log.removeHandler(h)
_h = logging.StreamHandler(); _h.setFormatter(logging.Formatter("%(message)s"))
log.addHandler(_h); log.propagate = False

def step(df: pd.DataFrame, name: str) -> pd.DataFrame:
    """Log shape after a build step; return df unchanged for chaining."""
    log.info("[%-26s] rows=%7d  unique_firms=%5d",
             name, len(df), df["permno"].nunique())
    return df

# --- BUILD PIPELINE (each line is one auditable decision) ------------------
b = step(crsp_raw.copy(),                                  "00 raw CRSP")
b = step(b.merge(ccm_primary[["permno", "gvkey"]], on="permno",
                 how="inner", validate="many_to_one"),    "01 link CCM (m:1)")
b["gvkey"] = b["gvkey"].astype(int)
b = step(pd.merge_asof(b.sort_values("date"),
                       fund_s, left_on="date", right_on="available_date",
                       by="gvkey", direction="backward"),  "02 merge_asof (6mo lag)")
b = step(b.query("book_equity > 0").copy(),                "03 drop neg book eq")
b["bm"] = b["book_equity"] / rng_me.lognormal(5.0, 1.0, size=len(b))
b["bm_w"] = b.groupby("date")["bm"].transform(lambda s: winsorize(s, 0.01, 0.99))
analysis = step(b.dropna(subset=["ret", "bm_w"]).copy(),   "04 listwise (ret,bm)")

out_path = os.path.join(WORKDIR, "value_panel.parquet")
try:
    analysis.to_parquet(out_path)                          # needs pyarrow
    saved = out_path
except Exception:
    saved = os.path.join(WORKDIR, "value_panel.csv")
    analysis.to_csv(saved, index=False)                    # offline fallback
print(f"\nSaved clean analysis panel -> {saved}")
print(f"final shape: {analysis.shape[0]:,} rows x {analysis.shape[1]} cols, "
      f"{analysis['permno'].nunique()} firms")
assert analysis["ret"].notna().all() and analysis["bm_w"].notna().all()

Your Turn

You have built Sam’s value panel end to end — linked with validate= on every merge, lagged with merge_asof and proved no peeking with the assertion, run the survivorship smell test, winsorized by date, probed missingness, and emitted the full build log. Now extend it:

  1. Rerun the machinery on Maya’s HMDA-style data. Synthesize a mortgage-application table with applicant_income Missing Not At Random — make income missing more often for denied applications. Then redo the Section 6 probe: the groupby("income_missing")["denied"].mean() should now differ sharply between groups. What does that tell you about dropna()? Build the crosswalk yourself (there is no CCM for HMDA).

  2. Break the look-ahead guard on purpose. Change direction="backward" to direction="forward" in the merge_asof and confirm the no-look-ahead assert fires. This is the one assertion you should never ship without — feel how it catches the bug a backtest never would.

  3. Winsorize or not? Add Devon’s daily crypto returns (a single day can move 40%). Decide: cap at 1/99, model the fat tail, or leave it — and defend the choice against the principle that the variable which is the phenomenon should not be capped away. Re-winsorize bm at 5/95 instead of 1/99 and watch how much the mean moves; that sensitivity is exactly why you fix the cutoff before you look at the result.