nb5.1 — FF92 Portfolio Sorts: the Size and Book-to-Market Cross-Section

Companion to Ch 5.1, the Reader’s Guide to Fama & French (1992), “The Cross-Section of Expected Stock Returns.” This is replication, not lecture: Ch 5.1 told you what to look for; here you build it.

Week 5 flips the verb. For four weeks you were a producer of empirical work; now you read a frontier paper the way a referee reads it and decide for yourself whether to believe it. Fama and French (1992) is the perfect text because it kills a reigning theory (CAPM) using exactly the machine you built in Lab 2 — the Fama–MacBeth cross-sectional regression. You already own the tools. This notebook is where you watch them do the killing, on data you control.

Sam’s stake. Sam trades on momentum and walked in believing “riskier stocks pay more” — the CAPM claim that a stock’s reward is proportional to its market beta. This notebook is the autopsy on that belief. We will show, on a panel built to be honest, that once you know a firm’s size and book-to-market, beta has almost nothing left to say.

The three replication exercises from Ch 5.1 §7, in order:

  1. Univariate sorts. Rank firms into size deciles, then into BE/ME deciles; read the average-return spread off each. Small beats big; value beats growth.
  2. The 5\times5 double sort. Form a size × BE/ME grid to disentangle the two effects — each cell holds firms alike on one dimension, spread on the other. This is “controlling for a confounder,” done nonparametrically with buckets.
  3. The Fama–MacBeth horse race. Three specifications — beta alone; size & BE/ME alone; all three together — and watch beta’s slope collapse to zero while size and BE/ME survive. This is Lab 2, verbatim.

A word on the data, because it is the whole reason this notebook is built the way it is.

0. Two data paths — and why we run on synthetic

CRSP (returns, prices, shares) and Compustat (book equity) are the workhorse pairing of empirical finance and the exact databases Fama and French used. They are also licensed. Per CONVENTIONS.md §5, licensed data stays read-only on GMU infrastructure (Hopper / WRDS Cloud); we never fabricate CRSP numbers and never export rows to a laptop.

So this notebook ships two paths:

  • Path A (real data, GMU only). The next cell shows the genuine wrds connection and the CRSP-monthly + Compustat-book-equity query pattern, with a pinned snapshot date. It is illustrative — it will not run off-campus, and that is intentional. Read it for the schema: permno, date, ret, market cap (prc × shrout), and book equity from Compustat (ceq + balance-sheet adjustments).
  • Path B (synthetic fallback, runs everywhere). A seeded, deterministic firm-month panel with the same schema and known premia planted into the data-generating process: small > big, value > growth, and beta priced into nothing once size and BE/ME are present. Because we built the truth, we can check the sorts and the Fama–MacBeth regression actually recover it.

Everything below Path A runs on Path B. The point of replication is not to memorize Fama and French’s exact 1963–1990 numbers; it is to build the machinery and confirm it finds a premium that is really there.

# =====================================================================
# PATH A: REAL CRSP / COMPUSTAT via WRDS — RUN ONLY ON GMU INFRASTRUCTURE
# Licensed data; read-only on Hopper / WRDS Cloud. Do NOT export rows.
# This cell is illustrative: it will not execute off-campus. Read it for the SCHEMA.
# =====================================================================
RUN_PATH_A = False   # flip to True only inside GMU's WRDS environment

if RUN_PATH_A:
    import wrds

    CRSP_SNAPSHOT = "2025-12-31"   # PINNED snapshot date — record this in your write-up
    db = wrds.Connection(wrds_username="your_gmu_wrds_user")  # creds via ~/.pgpass, never in code

    # --- CRSP monthly: returns + the pieces of market equity (price x shares outstanding) ---
    # Common stocks (shrcd 10,11) on NYSE/AMEX/Nasdaq (exchcd 1,2,3); nonfinancial filter applied later.
    crsp = db.raw_sql(
        """
        select a.permno, a.date, a.ret,
               abs(a.prc) * a.shrout as mktcap          -- market equity in $thousands
        from crsp.msf as a
        left join crsp.msenames as b
            on a.permno = b.permno
           and b.namedt <= a.date
           and a.date <= b.nameendt
        where b.shrcd in (10, 11)
          and b.exchcd in (1, 2, 3)
          and a.date between '1963-07-01' and '1990-12-31'   -- FF92 return window [CHECK exact]
        """,
        date_cols=["date"],
    )

    # --- Compustat annual: BOOK EQUITY (book value of common equity) ---
    # FF book equity = stockholders' equity + deferred taxes - preferred stock (with fallbacks).
    comp = db.raw_sql(
        """
        select gvkey, datadate, fyear,
               coalesce(seq, ceq + pstk, at - lt) as she,   -- stockholders' equity, with fallbacks
               coalesce(txditc, 0)                as deftax, -- deferred taxes + ITC
               coalesce(pstkrv, pstkl, pstk, 0)   as ps      -- preferred stock (redemption/liq/par)
        from comp.funda
        where indfmt = 'INDL' and datafmt = 'STD'
          and popsrc = 'D'    and consol  = 'C'
          and datadate between '1961-01-01' and '1989-12-31'
        """,
        date_cols=["datadate"],
    )
    db.close()

    # Book equity, then the CRITICAL timing lag (Ch 5.1 section 3, guards look-ahead bias):
    #   returns from July of year t .. June of t+1   <-  book equity from fiscal year ending t-1.
    # Link gvkey<->permno via crsp.ccmxpf_lnkhist, compute BE/ME = book_equity / market_equity_at_Dec(t-1),
    # then proceed with the SAME sort/FM code used on the synthetic panel below.
    # (Linking + lagging omitted here for brevity; the synthetic path exercises the identical logic.)

1. Setup — one seeded generator, headless plotting

We fix one seeded NumPy generator so the notebook reproduces bit-for-bit (the CONVENTIONS.md §5 rule), and force matplotlib’s non-interactive Agg backend so it runs headless — laptop, server, or CI — with no display attached. The seed is 42.

import matplotlib
matplotlib.use("Agg")   # headless, non-interactive backend

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)   # one seeded generator for the whole notebook

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

2. Path B — a synthetic firm-month panel with KNOWN premia

We simulate N = 1500 firms over T = 240 months (20 years), the same schema CRSP/Compustat hand you: a permno per firm, a month, a return ret, a mktcap (market equity), and a bm (book-to-market, BE/ME). The data-generating process is engineered so the three facts of FF92 are literally true in the data, with sizes you can read off the printout.

The firm characteristics (fixed per firm, the cross-sectional spread we sort on).

  • Size. Each firm has a log market cap logme \sim \mathcal N(0,1). Small firms have low logme. Real market caps are enormously right-skewed, so mktcap = exp(logme) gives the familiar few-giants-many-minnows shape.
  • Book-to-market. Each firm has a log BE/ME logbm \sim \mathcal N(0,1), drawn with a mild negative correlation to size (big firms lean a touch growth-y, as in the data), so the sorts have to work to disentangle them.
  • Beta. Each firm has a market beta beta centered at 1. Crucially, beta is correlated with size — small firms are higher-beta — exactly the tangle Ch 5.1 §5 warns about: a naive sort on beta is partly a disguised sort on smallness.

The return equation — where the premia are planted. Each month a firm’s return is

r_{it} = \underbrace{\beta_i\, r_{m,t}}_{\text{market exposure}} \;+\; \underbrace{\lambda_{\text{size}}\,(-z^{\text{size}}_i)}_{\text{small} > \text{big}} \;+\; \underbrace{\lambda_{\text{value}}\,(+z^{\text{bm}}_i)}_{\text{value} > \text{growth}} \;+\; \varepsilon_{it},

where z^{\text{size}} and z^{\text{bm}} are the standardized size and BE/ME characteristics. The planted monthly premia are LAM_SIZE = 0.0025 (0.25%/month for being one s.d. smaller) and LAM_VALUE = 0.0035 (0.35%/month for being one s.d. more value). Beta enters only through the market term \beta_i r_{m,t} — and because the market’s average excess return is set to zero (MKT_MEAN = 0), beta earns no average premium. That is the trapdoor: beta moves returns month to month, but over the long run it pays nothing once you hold size and BE/ME fixed. We built the DGP so “beta adds little” is true by construction; the FM regression in §6 then has to find that.

N, T = 1500, 240                 # 1500 firms, 240 months (20 years)
LAM_SIZE  = 0.0025               # planted SIZE premium: +0.25%/month per s.d. SMALLER
LAM_VALUE = 0.0035               # planted VALUE premium: +0.35%/month per s.d. higher BE/ME
MKT_MEAN  = 0.0                  # market average excess return = 0  => beta earns NO premium
MKT_SD    = 0.045                # monthly market volatility (~4.5%)

# --- Firm characteristics: fixed per firm (this is the cross-section we sort on) ---
logme = rng.normal(0.0, 1.0, size=N)                      # log market cap; small firms = low logme
# BE/ME mildly NEGATIVELY correlated with size: big firms lean a touch growth-y.
logbm = -0.30 * logme + np.sqrt(1 - 0.30**2) * rng.normal(0.0, 1.0, size=N)
# Beta correlated with SIZE: small firms are higher-beta (the FF92 section-5 tangle).
beta  = 1.0 - 0.35 * logme + 0.30 * rng.normal(0.0, 1.0, size=N)

mktcap = np.exp(logme)           # market equity: right-skewed, few giants & many minnows
bm     = np.exp(logbm)           # book-to-market ratio (BE/ME), strictly positive

# Standardize characteristics so the planted premia are "per one standard deviation"
z_size = (logme - logme.mean()) / logme.std()            # high => BIG firm
z_bm   = (logbm - logbm.mean()) / logbm.std()            # high => VALUE firm

# --- The market return path (one draw per month, shared by all firms) ---
mkt = rng.normal(MKT_MEAN, MKT_SD, size=T)               # market excess return r_{m,t}

# --- Build the firm-month panel ---
firm_idx  = np.repeat(np.arange(N), T)                   # 0,0,...,0,1,1,...  (each firm T times)
month_idx = np.tile(np.arange(T), N)

ret = (beta[firm_idx] * mkt[month_idx]                   # market exposure (beta priced into NOTHING)
       + LAM_SIZE  * (-z_size[firm_idx])                 # small > big
       + LAM_VALUE * ( z_bm[firm_idx])                   # value > growth
       + rng.normal(0.0, 0.06, size=N * T))              # idiosyncratic noise (~6%/month)

panel = pd.DataFrame({
    "permno": firm_idx,                                  # firm id (mimics CRSP permno)
    "month":  month_idx,
    "ret":    ret,
    "mktcap": mktcap[firm_idx],                           # market equity (size)
    "bm":     bm[firm_idx],                               # book-to-market (BE/ME)
    "beta":   beta[firm_idx],                             # pre-formation beta (known here; estimated in section 5)
})

print(f"panel: {len(panel):,} firm-months  ({N} firms x {T} months)")
print(f"planted premia (monthly):  size = {LAM_SIZE:.4f}   value = {LAM_VALUE:.4f}   beta-premium = 0 by design")
print(panel.head())

Looking at the world before we estimate

A referee never trusts a regression before eyeballing the raw data. Three sanity checks: market equity should be heavily right-skewed; size and beta should be negatively correlated (the §5 tangle); and the planted correlations should match what we typed.

fig, axes = plt.subplots(1, 2, figsize=(12, 4.0))

axes[0].hist(np.log(panel.groupby("permno")["mktcap"].first()), bins=40, color="#4C72B0", edgecolor="white")
axes[0].set_title("Cross-section of firm size (log market equity)")
axes[0].set_xlabel("log market cap"); axes[0].set_ylabel("number of firms")

axes[1].scatter(z_size, beta, s=8, alpha=0.35, color="#C44E52")
axes[1].set_title("Beta vs size: small firms are higher-beta (the FF92 tangle)")
axes[1].set_xlabel("standardized size  (low = small firm)"); axes[1].set_ylabel("beta")
fig.tight_layout(); fig.savefig("_fig_world.png", dpi=90); plt.close(fig)

print(f"corr(size, beta)   = {np.corrcoef(z_size, beta)[0,1]:+.3f}   <- small firms high-beta (the tangle to break)")
print(f"corr(size, BE/ME)  = {np.corrcoef(z_size, z_bm)[0,1]:+.3f}   <- big firms lean growth-y")
print("Saved beta-vs-size scatter to _fig_world.png")

3. Exercise 1 — univariate sorts (the gentlest replication)

Ch 5.1 §4 says: read the sort tables first, because they assume almost nothing. A univariate sort is the simplest possible nonparametric estimator of “does this characteristic line up with returns?” Rank every firm into ten size buckets (deciles), form the ten portfolios, and average each decile’s monthly return over the whole sample. Plot average return against the decile.

CAPM-or-not, the FF92 fact is: the line slopes down — small firms (decile 1) earn more than big firms (decile 10). We use pd.qcut to cut into deciles, and we average equal-weighted within each decile, each month, then across months — the honest order (a portfolio return is the cross-sectional mean that month; the sample mean of those is the bucket’s average).

def decile_sort(df, sort_col, ascending_label):
    """Sort firms into 10 buckets by sort_col; return each decile's avg monthly return.

    Equal-weighted: average within (decile, month) first, then across months. No chained indexing.
    """
    d = df.copy()
    d["decile"] = pd.qcut(d[sort_col], 10, labels=False) + 1          # deciles 1..10
    by_month = d.groupby(["decile", "month"], observed=True)["ret"].mean()   # portfolio return each month
    avg = by_month.groupby("decile", observed=True).mean()                   # average over months
    avg.name = "avg_ret"
    return avg

size_sort = decile_sort(panel, "mktcap", "small->big")
print("SIZE deciles (1 = smallest, 10 = biggest):  average monthly return")
print((size_sort * 100).round(3).to_string())

spread = (size_sort.loc[1] - size_sort.loc[10]) * 100
print(f"\nSmall-minus-Big spread (decile 1 - decile 10) = {spread:+.3f} %/month")
print(f"  -> POSITIVE means small beats big (the FF92 size effect).  Planted: ~{LAM_SIZE*100:.2f}%/sd.")
assert size_sort.loc[1] > size_sort.loc[10], "size effect should be present: small > big"

Now the same machine on book-to-market. Slide from decile 1 (low BE/ME = growth) to decile 10 (high BE/ME = value) and watch average return rise — value beats growth. Two live patterns, built by the same three lines of code.

bm_sort = decile_sort(panel, "bm", "growth->value")
print("BE/ME deciles (1 = growth, 10 = value):  average monthly return")
print((bm_sort * 100).round(3).to_string())

spread_v = (bm_sort.loc[10] - bm_sort.loc[1]) * 100
print(f"\nValue-minus-Growth spread (decile 10 - decile 1) = {spread_v:+.3f} %/month")
print(f"  -> POSITIVE means value beats growth (the FF92 value effect).  Planted: ~{LAM_VALUE*100:.2f}%/sd.")
assert bm_sort.loc[10] > bm_sort.loc[1], "value effect should be present: value > growth"

# Plot both sorts side by side
fig, axes = plt.subplots(1, 2, figsize=(12, 4.0))
axes[0].plot(size_sort.index, size_sort.values * 100, "o-", color="#4C72B0")
axes[0].set_title("Size sort: small beats big (line slopes DOWN)")
axes[0].set_xlabel("size decile (1=small, 10=big)"); axes[0].set_ylabel("avg return %/month")
axes[1].plot(bm_sort.index, bm_sort.values * 100, "o-", color="#55A868")
axes[1].set_title("BE/ME sort: value beats growth (line slopes UP)")
axes[1].set_xlabel("BE/ME decile (1=growth, 10=value)"); axes[1].set_ylabel("avg return %/month")
fig.tight_layout(); fig.savefig("_fig_univariate.png", dpi=90); plt.close(fig)
print("Saved univariate sorts to _fig_univariate.png")

Break it on purpose: the look-ahead trap

Ch 5.1 §3 and §7 hammer one point: returns from July of year t are matched to accounting data known before July of t, with a deliberate timing gap, to block look-ahead bias — using information no investor could have had yet. To feel why the gap matters, do the wrong thing on purpose. Suppose the size you sort on is contaminated with a peek at the firm’s own future return (as it would be if you let this month’s price, which already moved with this month’s return, define “size” for explaining this month’s return). The sort then manufactures a spread that is not a premium at all — it is the outcome leaking into the regressor.

We add a small slug of the contemporaneous return into a fake “size” signal and re-sort. The spread inflates relative to the honest sort: the picture changes, and a referee who didn’t know the timing rule would over-credit the size effect.

# Contaminate the sort variable with a peek at the SAME month's return (the look-ahead mistake).
contaminated = panel.copy()
contaminated["mktcap_peek"] = contaminated["mktcap"] * np.exp(-3.0 * contaminated["ret"])  # price-already-moved leak

honest = decile_sort(panel, "mktcap", "")                       # clean lagged-style sort
peeked = decile_sort(contaminated, "mktcap_peek", "")           # look-ahead-contaminated sort

honest_spread = (honest.loc[1] - honest.loc[10]) * 100
peek_spread   = (peeked.loc[1] - peeked.loc[10]) * 100
print(f"Honest size spread (no peek)        = {honest_spread:+.3f} %/month")
print(f"Look-ahead-contaminated size spread = {peek_spread:+.3f} %/month")
print(f"  -> contamination INFLATES the apparent premium by {peek_spread - honest_spread:+.3f} %/month.")
print("  This is why FF92's July-to-June timing gap is not bookkeeping fussiness but a guardrail.")

4. Exercise 2 — the 5×5 double sort (disentangle size from value)

A univariate sort cannot tell you whether the value premium is really about BE/ME or is size in disguise (recall size and BE/ME are correlated). The double sort breaks the tangle. Cut firms into 5 size quintiles and, independently, 5 BE/ME quintiles, and average returns in each of the 25 cells. Now each cell holds firms alike in size and alike in BE/ME, so:

  • Read across a row (fix size, vary BE/ME): the value effect, holding size roughly fixed.
  • Read down a column (fix BE/ME, vary size): the size effect, holding value roughly fixed.

This is exactly “controlling for a confounder,” done nonparametrically with buckets instead of a regression slope — the same logic Ch 5.1 §5 calls the paper’s cleverest move.

def quintile(x):
    return pd.qcut(x, 5, labels=False) + 1     # quintiles 1..5

dbl = panel.copy()
dbl["size_q"] = quintile(dbl["mktcap"])        # 1 = small ... 5 = big
dbl["bm_q"]   = quintile(dbl["bm"])            # 1 = growth ... 5 = value

# Portfolio return each (size_q, bm_q, month), then average across months -> 5x5 table.
by_m  = dbl.groupby(["size_q", "bm_q", "month"], observed=True)["ret"].mean()
table = by_m.groupby(["size_q", "bm_q"], observed=True).mean().unstack("bm_q") * 100  # %/month

table.index   = [f"S{q} (small)" if q == 1 else f"S{q} (big)" if q == 5 else f"S{q}" for q in table.index]
table.columns = [f"B{q} (grow)"  if q == 1 else f"B{q} (val)"  if q == 5 else f"B{q}" for q in table.columns]
print("5x5 double sort: average monthly return (%), rows = size, cols = BE/ME")
print(table.round(3).to_string())

Now read the table the way a referee does: compute the value spread within each size row and the size spread within each BE/ME column. If both survive across the board — not just on the diagonal — the two effects are genuinely distinct, not one masquerading as the other.

val_spread_by_size  = (table.iloc[:, -1] - table.iloc[:, 0])   # last BE/ME col minus first, per size row
size_spread_by_bm   = (table.iloc[0, :] - table.iloc[-1, :])   # smallest size row minus biggest, per BE/ME col

print("VALUE spread within each size row (value - growth), %/month:")
print(val_spread_by_size.round(3).to_string())
print("\nSIZE spread within each BE/ME column (small - big), %/month:")
print(size_spread_by_bm.round(3).to_string())

print(f"\nBoth spreads positive in EVERY bucket? value: {(val_spread_by_size > 0).all()};  "
      f"size: {(size_spread_by_bm > 0).all()}")
assert (val_spread_by_size > 0).all(), "value premium should hold within every size row"
assert (size_spread_by_bm  > 0).all(), "size premium should hold within every BE/ME column"

fig, ax = plt.subplots(figsize=(6.5, 5.0))
im = ax.imshow(table.values, cmap="viridis", aspect="auto")
ax.set_xticks(range(5)); ax.set_xticklabels(table.columns, rotation=30, ha="right")
ax.set_yticks(range(5)); ax.set_yticklabels(table.index)
ax.set_title("5x5 sort: avg return %/month\n(brighter = higher; bright corner = small-value)")
fig.colorbar(im, ax=ax, label="%/month"); fig.tight_layout()
fig.savefig("_fig_doublesort.png", dpi=90); plt.close(fig)
print("Saved double-sort heatmap to _fig_doublesort.png")

5. Exercise 3 — the Fama–MacBeth horse race

The sorts showed the patterns; the Fama–MacBeth regression quantifies them and lets the characteristics fight each other in one specification. This is Lab 2, verbatim, so we reuse its machine exactly.

Pass 1 — pre-formation betas, estimated by portfolio. Ch 5.1 §2 and §6 stress that beta is estimated, not observed, and individual-stock betas are noisy. A noisy regressor suffers attenuation bias (the errors-in-variables problem from Week 2): its slope is pulled toward zero, so a skeptic could claim beta’s flat slope is just noise hiding a real effect. Fama and French’s fix — which we copy — is to assign betas via portfolios: sort firms into beta groups, estimate each group’s beta on its longer, smoother return series, and let every firm inherit its group’s beta. Portfolio betas are far less noisy, so attenuation cannot be the easy excuse.

# --- PASS 1: portfolio (pre-formation) betas to tame errors-in-variables ---
# Sort firms into 20 beta groups by their realized covariance with the market, then estimate
# each GROUP's beta on the group's average return series (longer, smoother => less noisy).
wide = panel.pivot_table(index="month", columns="permno", values="ret")  # T x N return matrix
firm_cov = wide.apply(lambda col: np.cov(col.values, mkt)[0, 1] / np.var(mkt))  # rough per-firm beta

beta_group = pd.qcut(firm_cov, 20, labels=False)                  # 20 beta portfolios
group_ret  = wide.T.groupby(beta_group.values).mean().T           # T x 20 group-average returns
# Each group's beta from a time-series regression of group return on the market.
X_mkt = sm.add_constant(mkt)
group_beta = group_ret.apply(lambda col: sm.OLS(col.values, X_mkt).fit().params[1])

# Map every firm to its group's (smooth) beta -> the pre-formation beta we use in Pass 2.
firm_beta_hat = beta_group.map(group_beta)                        # indexed by permno
firm_beta_hat.name = "beta_hat"

corr_true = np.corrcoef(firm_beta_hat.values, beta)[0, 1]
print(f"corr(portfolio beta_hat, TRUE beta) = {corr_true:+.3f}  (portfolio grouping recovers true beta well)")
print(firm_beta_hat.describe().round(3).to_string())

Pass 2 — one cross-sectional regression per month. Attach to every firm-month three regressors: the portfolio beta_hat, log size, and log BE/ME (standardized, so slopes are comparable and read in “per standard deviation” units, matching the planted premia). Then, separately in each month, regress that month’s returns on the chosen regressors. That yields a time series of monthly slopes \hat\gamma_{j,t} — one per regressor per month.

The Fama–MacBeth point estimate is the average of the monthly slopes, and — the move that makes the whole thing honest — its standard error is the ordinary standard error of that mean: the sample standard deviation of the slopes over \sqrt{T}.

\hat\gamma_j = \frac{1}{T}\sum_{t=1}^T \hat\gamma_{j,t}, \qquad \widehat{\operatorname{se}}(\hat\gamma_j) = \frac{\operatorname{sd}_t(\hat\gamma_{j,t})}{\sqrt{T}}.

Why this is the honest SE: in any month, every stock is shoved by the same market-wide shock, so within-month errors are violently correlated. Pooling all firm-months into one OLS would treat correlated observations as independent and report a standard error far too small (the “time effect” of Ch 2.4). Fama–MacBeth dissolves that — one month, one vote — exactly as you proved in Lab 2.

# Attach regressors to the panel (standardized, "per one s.d." units).
reg = panel.merge(firm_beta_hat, left_on="permno", right_index=True, how="left").copy()
reg["log_size"] = np.log(reg["mktcap"])
reg["log_bm"]   = np.log(reg["bm"])
for c in ["beta_hat", "log_size", "log_bm"]:
    reg[c] = (reg[c] - reg[c].mean()) / reg[c].std()        # standardize cross-section

def fama_macbeth(df, regressors):
    """Lab-2 machine: one cross-sectional OLS per month, then average the slopes.

    Returns a tidy frame: coef (FM point estimate), se (FM standard error), tstat, for each regressor + const.
    """
    cols = ["const"] + regressors
    slopes = []
    for _, g in df.groupby("month", observed=True):
        X = sm.add_constant(g[regressors].to_numpy(), has_constant="add")
        b = sm.OLS(g["ret"].to_numpy(), X).fit().params         # length = 1 + len(regressors)
        slopes.append(b)
    S = pd.DataFrame(slopes, columns=cols)                      # T x (1 + k) table of monthly slopes
    Tn = len(S)
    out = pd.DataFrame({
        "coef": S.mean(),
        "se":   S.std(ddof=1) / np.sqrt(Tn),                    # FM SE = sd(slopes)/sqrt(T)
    })
    out["tstat"] = out["coef"] / out["se"]
    return out, S

The three columns: a sequence of horse races

Ch 5.1 §4 says to read the FM table as a sequence, not one block:

  1. Beta alone — the formal CAPM test. Its slope should be small and statistically indistinguishable from zero (we built the market mean to be zero).
  2. Size & BE/ME alone — both slopes large and significant (size negative: smaller firms, higher returns once we measure in raw log-size, where small = low; BE/ME positive: value wins).
  3. The kitchen sink — beta, size, BE/ME together. The decisive cell: size and BE/ME keep their slopes; beta’s stays near zero and insignificant.

Note on signs: our log_size is increasing in firm bigness, so a negative log_size slope is the size premium (small firms earn more). The planted per-s.d. magnitudes are size 0.25\% and value 0.35\% per month; standardizing changes the scale slightly but the signs and significance are what we read.

fm_beta,  _ = fama_macbeth(reg, ["beta_hat"])
fm_sv,    _ = fama_macbeth(reg, ["log_size", "log_bm"])
fm_all, S_all = fama_macbeth(reg, ["beta_hat", "log_size", "log_bm"])

def show(name, res):
    print(f"\n=== {name} ===")
    disp = (res[["coef", "se", "tstat"]] * [10000, 10000, 1]).round(2)   # coef/se in basis points
    disp.columns = ["coef (bps/mo)", "se (bps)", "t-stat"]
    print(disp.to_string())

show("Column 1: BETA alone (the CAPM test)", fm_beta)
show("Column 2: SIZE & BE/ME alone", fm_sv)
show("Column 3: KITCHEN SINK (beta + size + BE/ME)", fm_all)

The verdict, in our own t-statistics

The headline of the entire paper is the dog that didn’t bark: the number that should be big under CAPM — beta’s slope — is the number that isn’t. We check it numerically: beta’s t-stat is small in both the beta-alone column and the kitchen sink, while size and BE/ME stay significant (|t| well above 2) throughout.

t_beta_alone = fm_beta.loc["beta_hat", "tstat"]
t_beta_sink  = fm_all.loc["beta_hat",  "tstat"]
t_size_sink  = fm_all.loc["log_size",  "tstat"]
t_bm_sink    = fm_all.loc["log_bm",    "tstat"]

print(f"beta t-stat, ALONE        = {t_beta_alone:+.2f}   <- CAPM test: should be small (it is)")
print(f"beta t-stat, KITCHEN SINK = {t_beta_sink:+.2f}   <- stays near zero with size & BE/ME present")
print(f"size  t-stat, kitchen sink = {t_size_sink:+.2f}  (negative = small firms earn more)")
print(f"BE/ME t-stat, kitchen sink = {t_bm_sink:+.2f}  (positive = value beats growth)")

# The paper's thesis, asserted: beta loses, size & value survive together.
assert abs(t_beta_sink) < 2.0,  "beta should NOT be significant once size & BE/ME enter"
assert t_size_sink < -2.0,      "size premium should survive (negative, significant)"
assert t_bm_sink   >  2.0,      "value premium should survive (positive, significant)"
print("\nVERDICT: beta's slope collapses; size and BE/ME survive together. FF92 reproduced.")

Why the FM standard error is what it is

The FM SE is built from the month-to-month scatter of the slopes. Plot the monthly beta slope \hat\gamma_{\beta,t} from the kitchen-sink regression: it swings wildly around zero (driven by whether the market rose or fell that month), and its average is near zero — that average is the FM point estimate, and the roughness of the series is the FM standard error. Contrast the size and value slope series, which bounce around a clearly nonzero mean.

fig, ax = plt.subplots(figsize=(11, 4.0))
ax.plot(S_all.index, S_all["beta_hat"] * 100, color="#C44E52", lw=0.8, alpha=0.8, label="beta slope (mean ~0)")
ax.plot(S_all.index, S_all["log_bm"]   * 100, color="#55A868", lw=0.8, alpha=0.8, label="BE/ME slope (mean >0)")
ax.axhline(0, color="black", lw=0.8)
ax.axhline(S_all["beta_hat"].mean() * 100, color="#C44E52", ls="--", lw=1.4)
ax.axhline(S_all["log_bm"].mean()   * 100, color="#55A868", ls="--", lw=1.4)
ax.set_title("Monthly FM slopes (kitchen sink): beta scatters around 0; BE/ME around a positive mean")
ax.set_xlabel("month"); ax.set_ylabel("monthly slope (%/sd)"); ax.legend(loc="upper right")
fig.tight_layout(); fig.savefig("_fig_fm_slopes.png", dpi=90); plt.close(fig)

print(f"beta  slope: mean {S_all['beta_hat'].mean()*10000:+.2f} bps, sd {S_all['beta_hat'].std()*10000:.1f} bps "
      f"-> big sd, ~0 mean => insignificant")
print(f"BE/ME slope: mean {S_all['log_bm'].mean()*10000:+.2f} bps, sd {S_all['log_bm'].std()*10000:.1f} bps "
      f"-> clearly nonzero mean => significant")
print("Saved FM slope series to _fig_fm_slopes.png")

6. What to carry forward

You just rebuilt the spine of Fama & French (1992) on a panel where the truth was known:

  • Univariate sorts recovered the planted facts: small beats big, value beats growth — read straight off the decile means, no model imposed.
  • The 5×5 double sort showed both premia survive within every bucket of the other, so neither is a disguise for the other. That is nonparametric confounder control.
  • The Fama–MacBeth horse race put the characteristics in one ring: beta’s slope collapsed toward zero and lost significance, while size and BE/ME kept large, significant slopes — the dog that didn’t bark. The FM standard error came from the month-to-month scatter of slopes, one vote per month, exactly as in Lab 2.
  • We saw why beta’s flatness is not a measurement artifact: portfolio betas tracked the true betas tightly (corr above), so attenuation cannot rescue CAPM here. The premium that should have been large is genuinely near zero.

The one-line spec, in CONVENTIONS §4 discipline: outcome = a firm’s monthly return; key regressors = portfolio beta, log size, log BE/ME; controls = the other characteristics in the same regression; fixed effects = none (FM sense); standard errors = Fama–MacBeth (month-to-month scatter of slopes); sample = the synthetic panel (Path B) or 1963–1990 nonfinancial NYSE/AMEX/Nasdaq firms (Path A); identifying claim = none, causal — these are predictive associations in the cross-section.

Your Turn

Each of these is a small edit to the code above. Make the change, re-run, and predict the direction before you look.

  1. Resurrect beta. Set MKT_MEAN = 0.005 (a positive market risk premium) and re-run §5. Now beta should earn a premium. Does its FM slope turn positive and significant? You have just shown the machine is not rigged against beta — it reports a beta premium when one truly exists.
  2. Kill the value premium. Set LAM_VALUE = 0.0 and re-run the BE/ME sort (§3) and the kitchen sink (§5). Confirm the value spread collapses to noise and log_bm’s t-stat drops below 2. This is your falsification check: the method finds a premium only when one is planted.
  3. Individual betas (errors-in-variables, on purpose). In Pass 1, replace the portfolio firm_beta_hat with each firm’s own noisy firm_cov. Re-run the kitchen sink. The beta slope should attenuate even harder toward zero and the corr-with-truth should drop — a live demonstration of why Fama and French formed portfolios (Ch 5.1 §6).
  4. More months vs. more firms. Double T (to 480), then separately double N (to 3000), each time re-running §5. Which one shrinks the FM standard error more? (Hint: the FM SE lives in \sqrt{T}, the count of time clusters — the Ch 2.4 lesson wearing a new hat.)