nb5.2 — FF93 Factor Regressions: Loadings, R^2, \alpha, and the GRS Test

Companion to Ch 5.2, the reader’s guide to Fama & French (1993), “Common Risk Factors in the Returns on Stocks and Bonds.”

The guide drew a sharp line between two papers that students forever confuse. FF92 ran a cross-section and learned that size and book-to-market are priced. FF93 ran a time series and asked the structurally different question: is there a small set of factors — market, size (SMB), value (HML) — such that regressing any test portfolio’s excess return on those three has (a) a high R^2, so the factors capture the common month-to-month movement, and (b) an intercept \alpha indistinguishable from zero, so no priced risk is left on the table? High R^2 + zero \alpha = a working factor model. That two-column reading is the literacy this whole week teaches.

This notebook turns that reading into running code. We will:

Get the data two ways. First a real, free pull of Ken French’s three factors and his 25 size×BE/ME portfolios via pandas_datareader — wrapped in try/except so the notebook never hard-fails offline. If the network (or the library) is unavailable, we fall back to a
seeded synthetic data-generating process built from a known three-factor model with tiny alphas, so the regressions can recover the truth and we can grade ourselves.
Run the 25 time-series regressions, tabulate loadings \hat b,\hat s,\hat h, the intercept \hat\alpha with its t, and R^2, and read the two-column story.
Compare to CAPM (market-only) to see R^2 rise and intercepts shrink when SMB/HML are added.
Build the GRS F-statistic by hand for the joint null \alpha_1=\cdots=\alpha_{25}=0, and interpret it through Week 1’s multiple-testing lens.

Setup

One seeded generator for the whole notebook (the reproducibility rule from CONVENTIONS.md), and matplotlib forced onto its headless Agg backend so the notebook runs anywhere, screen or no screen.

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

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

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

pd.set_option("display.float_format", lambda v: f"{v:8.4f}")
np.set_printoptions(precision=6, suppress=True)
print("numpy", np.__version__, "| pandas", pd.__version__, "| statsmodels", statsmodels.__version__)

Path A — the real, free Ken French data (attempt the pull)

Ken French’s Data Library is free, public data — not the licensed CRSP/Compustat we keep on GMU infrastructure. Anyone can download it. The convenient door is pandas_datareader.data.DataReader(name, "famafrench"), which fetches a named dataset and hands back a dictionary of monthly tables.

Two names we want:

"F-F_Research_Data_Factors" — the monthly Mkt-RF, SMB, HML factors plus the risk-free rate RF (all in percent per month).
"25_Portfolios_5x5" — the 25 test portfolios on a 5×5 size×BE/ME sort (also percent per month, as raw returns, so we subtract RF ourselves to get excess returns).

We wrap the whole thing in try/except. Networks fail, libraries go missing, French occasionally reorganizes a file — none of that should crash a teaching notebook. If anything goes wrong we set USED_REAL = False and drop to the synthetic generator in Path B. Nothing downstream cares which path produced the numbers, which is exactly the discipline you want: write the analysis once, feed it whichever data you have.

# Attempt the free, real Ken French pull. Falls through to synthetic on ANY failure.
USED_REAL = False
ff_factors = ff_port = None
START, END = "1963-07", "1991-12"  # the classic FF93 stock window

try:
    import pandas_datareader.data as web
    fac = web.DataReader("F-F_Research_Data_Factors", "famafrench",
                         start=START, end=END)[0]
    prt = web.DataReader("25_Portfolios_5x5", "famafrench",
                         start=START, end=END)[0]
    # Align on common months, then form excess returns for the 25 portfolios.
    common = fac.index.intersection(prt.index)
    fac = fac.loc[common].copy()
    prt = prt.loc[common].copy()
    rf = fac["RF"]
    factors = fac[["Mkt-RF", "SMB", "HML"]].copy()
    excess = prt.sub(rf, axis=0).copy()          # raw return minus RF, month by month
    port_names = list(excess.columns)
    USED_REAL = True
    print(f"Real Ken French data loaded: {len(factors)} months, "
          f"{excess.shape[1]} test portfolios.")
except Exception as e:
    # ImportError (no pandas_datareader), URLError (offline), HTTP/parse errors -> fallback.
    print(f"Real pull unavailable ({type(e).__name__}: {e}).")
    print("Falling back to the seeded synthetic DGP in Path B.")

Path B — the seeded synthetic three-factor world (the offline guarantee)

If Path A did not load, we manufacture a world where the three-factor model is true by construction, so we know the right answer and can check that our regressions recover it. This is the same Monte-Carlo move from Week 2: write down the data-generating process (DGP), simulate from it, and confirm the estimator finds the parameters you planted.

The DGP mirrors FF93’s economics:

R_{it}-R_{ft} \;=\; \alpha_i \;+\; b_i\,(R_{Mt}-R_{ft}) \;+\; s_i\,\text{SMB}_t \;+\; h_i\,\text{HML}_t \;+\; \varepsilon_{it}.

We lay the 25 portfolios on a 5×5 grid of size (rows S1→S5, small→big) and book-to-market (columns BM1→BM5, growth→value), and we plant loadings that march with the sorts, exactly the pattern the guide told you to look for:

SMB loading s_i runs from strongly positive in the small row to negative in the big row (small stocks co-move with the size factor).
HML loading h_i runs from negative in the growth column to strongly positive in the value column (value stocks co-move with the value factor).
Market loading b_i sits near 1 everywhere, drifting a touch above 1 for small stocks.

The planted \alpha_i are tiny — a few basis points per month of mean-zero noise — so the model is almost exact. That mimics reality: FF93’s intercepts are small but not all perfectly zero, with small-growth the stubborn corner. We deliberately give the small-growth cell (S1,BM1) a visibly larger alpha so it misbehaves the way the real data do.

if not USED_REAL:
    T = 342  # months, matching the classic July 1963 - Dec 1991 window length
    months = pd.period_range("1963-07", periods=T, freq="M")

    # --- factor returns (percent/month), loosely calibrated to FF magnitudes ---
    rf = pd.Series(0.40, index=months, name="RF")  # ~0.4%/month flat risk-free
    mkt_rf = rng.normal(0.50, 4.5, T)   # equity premium ~0.5%/mo, vol ~4.5
    smb    = rng.normal(0.30, 3.0, T)   # size premium, smaller and less volatile
    hml    = rng.normal(0.45, 3.0, T)   # value premium (FF historical ~0.3-0.5%/mo)
    factors = pd.DataFrame({"Mkt-RF": mkt_rf, "SMB": smb, "HML": hml}, index=months)

    # --- planted loadings on the 5x5 grid (size rows S1..S5, BM cols BM1..BM5) ---
    size_levels = ["S1", "S2", "S3", "S4", "S5"]   # small -> big
    bm_levels   = ["BM1", "BM2", "BM3", "BM4", "BM5"]  # growth -> value
    s_row = np.array([ 1.30,  0.95,  0.60,  0.25, -0.20])   # SMB loading by size (small high)
    h_col = np.array([-0.45, -0.10,  0.20,  0.55,  0.90])   # HML loading by BM (value high)

    port_names, true_b, true_s, true_h, true_a = [], [], [], [], []
    cols = {}
    for r, sz in enumerate(size_levels):
        for c, bm in enumerate(bm_levels):
            name = f"{sz} {bm}"
            b = 1.0 + 0.06 * (4 - r)            # market beta ~1, a touch higher for small
            s = s_row[r]
            h = h_col[c]
            # tiny alphas (bp/month); make small-growth (S1,BM1) the stubborn corner
            a = rng.normal(0.0, 0.03)
            if sz == "S1" and bm == "BM1":
                a = -0.28  # the notorious small-growth pricing error
            eps = rng.normal(0.0, 1.5, T)       # idiosyncratic, diversifiable noise
            ret = a + b * mkt_rf + s * smb + h * hml + eps
            cols[name] = ret
            port_names.append(name)
            true_b.append(b); true_s.append(s); true_h.append(h); true_a.append(a)

    excess = pd.DataFrame(cols, index=months)   # these ARE excess returns by construction
    truth = pd.DataFrame({"true_alpha": true_a, "true_b": true_b,
                          "true_s": true_s, "true_h": true_h}, index=port_names)
    print(f"Synthetic data built: {T} months, {len(port_names)} test portfolios.")
    print("\nPlanted loadings (first 6 portfolios):")
    print(truth.head(6))
else:
    truth = None  # no ground truth when we used the real data
print("\nDATA SOURCE USED:", "REAL Ken French" if USED_REAL else "SYNTHETIC fallback")

A first look at the factors

Whichever path ran, we now hold a factors frame (Mkt-RF, SMB, HML) and an excess frame of 25 test-portfolio excess returns. Before regressing, do the cheap sanity checks the guide’s first table demands: are the average factor returns positive (did small and value actually pay?), and are the factors roughly orthogonal (low correlation, so the loadings stay interpretable and standard errors stay small)? In real data SMB/HML correlate only weakly with the market and with each other — by design, from the 2×3 averaging.

summary = pd.DataFrame({
    "mean (%/mo)": factors.mean(),
    "std (%/mo)":  factors.std(),
    "annualized mean (%)": factors.mean() * 12,
})
print("Factor summary statistics")
print(summary, "\n")

print("Factor correlation matrix (should be near-diagonal)")
print(factors.corr())

The pattern to be explained: average excess returns on the 25 portfolios

This is the explanandum — the grid the model has to reproduce. Read it as a 5×5 table. Moving across BE/ME (growth → value, BM1 → BM5) average returns should rise: the value effect. Moving down size (small → big, S1 → S5) they should fall: the size effect. The highest-returning corner is small-value (top-right); small-growth (top-left) is the famous problem child. Burn this grid in — every later table is asking “did the three factors reproduce this?”

# Reshape the 25 mean excess returns into the 5x5 size x BM grid.
mean_exc = excess.mean()

def to_grid(series, rows=("S1","S2","S3","S4","S5"),
            cols=("BM1","BM2","BM3","BM4","BM5")):
    "Lay a 25-vector of 'Sx BMy' values onto the 5x5 grid; tolerant of name variants."
    g = pd.DataFrame(index=list(rows), columns=list(cols), dtype=float)
    for r in rows:
        for c in cols:
            for nm in series.index:
                key = str(nm).upper().replace("-", " ")
                if (r in key.split()) and (c in key.split()):
                    g.loc[r, c] = series[nm]
                    break
    return g

grid = to_grid(mean_exc)
if grid.isna().all().all():        # real-data column labels differ -> fall back to raw order
    grid = pd.DataFrame(mean_exc.values.reshape(5, 5),
                        index=["S1","S2","S3","S4","S5"],
                        columns=["BM1","BM2","BM3","BM4","BM5"])
print("Average excess return (%/mo), rows = size (small->big), cols = BM (growth->value)")
print(grid)

The heart of the paper: 25 time-series regressions

For each test portfolio i we run plain multiple OLS over all T months (Ch 2.1):

R_{it}-R_{ft} \;=\; \alpha_i + b_i\,(R_{Mt}-R_{ft}) + s_i\,\text{SMB}_t + h_i\,\text{HML}_t + \varepsilon_{it}.

statsmodels does the algebra; our job is to harvest the right four numbers per regression — the intercept \hat\alpha_i and its t, the three loadings, and the R^2 — and stack the 25 residual series into a T\times 25 matrix \hat{\mathbf E}. We need that residual matrix later: the GRS test depends on how the portfolios’ residuals co-move, not just on each regression in isolation.

Xf = sm.add_constant(factors[["Mkt-RF", "SMB", "HML"]])  # T x 4 design (const + 3 factors)

rows = []
resid_mat = np.empty((len(excess), len(port_names)))  # T x N residuals for GRS
for j, nm in enumerate(port_names):
    y = excess[nm]
    res = sm.OLS(y.values, Xf.values).fit()
    alpha, b, s, h = res.params      # order follows add_constant: [const, Mkt-RF, SMB, HML]
    rows.append({
        "portfolio": nm,
        "alpha": alpha, "t(alpha)": res.tvalues[0],
        "b_mkt": b, "s_smb": s, "h_hml": h,
        "R2": res.rsquared,
    })
    resid_mat[:, j] = res.resid

ff3 = pd.DataFrame(rows).set_index("portfolio")
print("Three-factor time-series regressions (first 8 of 25 portfolios)")
print(ff3.head(8))
print(f"\nMean R^2 across the 25 portfolios: {ff3['R2'].mean():.4f}")

Pass 1 — the R^2 column (common variation)

Scan only the R^2 column. They should be high across the board: the three factors explain the large majority of each portfolio’s month-to-month variance. This is the “common risk factors” of the title — the shared movement that cannot be diversified away. (In the synthetic world the ceiling is set by how much idiosyncratic noise we planted; in real FF data the 25 R^2s typically sit in the 0.90s.)

r2_grid = to_grid(ff3["R2"])
if r2_grid.isna().all().all():
    r2_grid = pd.DataFrame(ff3["R2"].values.reshape(5, 5),
                           index=["S1","S2","S3","S4","S5"],
                           columns=["BM1","BM2","BM3","BM4","BM5"])
print("R^2 grid (rows = size, cols = BM)")
print(r2_grid)
print(f"\nMin R^2 = {ff3['R2'].min():.4f}   Max R^2 = {ff3['R2'].max():.4f}")

Pass 2 — the loadings (do they march with the sorts?)

Now read the loadings along the dimensions the portfolios were sorted on. The SMB loading \hat s_i should run down the size grid from strongly positive (small) toward zero/negative (big). The HML loading \hat h_i should run across the BM grid from negative (growth) to strongly positive (value). That monotone march is the model “recognizing” the very characteristics used to form the portfolios — what people mean when they say the loadings “make sense.”

def show_grid(col, label):
    g = to_grid(ff3[col])
    if g.isna().all().all():
        g = pd.DataFrame(ff3[col].values.reshape(5, 5),
                         index=["S1","S2","S3","S4","S5"],
                         columns=["BM1","BM2","BM3","BM4","BM5"])
    print(f"{label} (rows = size small->big, cols = BM growth->value)")
    print(g, "\n")
    return g

g_smb = show_grid("s_smb", "SMB loading  s_hat  (should fall down the size rows)")
g_hml = show_grid("h_hml", "HML loading  h_hat  (should rise across the BM cols)")
print("Row-mean SMB loading by size:", g_smb.mean(axis=1).round(3).to_dict())
print("Col-mean HML loading by BM:  ", g_hml.mean(axis=0).round(3).to_dict())

Pass 3 — the intercepts (pricing errors)

Finally the \hat\alpha_i column and its t-statistics. The claim is that they are economically small and mostly statistically indistinguishable from zero — tiny relative to the average returns in the explanandum grid. Watch the small-growth corner (S1, BM1): it is the model’s most stubborn pricing error in real data, and we planted exactly that pathology into the synthetic world so you can see it bite.

A word of Week-1 caution before you celebrate or panic over any single cell: with 25 separate t-tests at the 5% level, you expect about one or two to clear |t|>2 by chance alone. That is the multiple-comparisons trap (Ch 1.5). Eyeballing 25 t-stats is not a test. The honest instrument is the joint test — coming up next.

alpha_grid = show_grid("alpha", "Intercept alpha_hat (%/mo) - pricing errors")
print("Number of |t(alpha)| > 2 across the 25 portfolios:",
      int((ff3["t(alpha)"].abs() > 2).sum()),
      "(expect ~1-2 by chance even if the model is true)")
print(f"Mean |alpha| = {ff3['alpha'].abs().mean():.4f} %/mo")

# Largest pricing error: should be the small-growth corner in both data paths.
worst = ff3["alpha"].abs().idxmax()
print(f"\nLargest |alpha|: {worst}  ->  alpha = {ff3.loc[worst,'alpha']:.4f}, "
      f"t = {ff3.loc[worst,'t(alpha)']:.2f}")

Did we recover the truth? (synthetic path only)

When we built the data ourselves we know every planted parameter, so we can grade the regressions: estimated loadings should sit on top of the true ones, and estimated alphas should hug the tiny planted alphas. On the real-data path there is no ground truth, so we skip this and trust the qualitative reading above.

if not USED_REAL:
    chk = truth.copy()
    chk["est_b"] = ff3["b_mkt"].values
    chk["est_s"] = ff3["s_smb"].values
    chk["est_h"] = ff3["h_hml"].values
    chk["est_alpha"] = ff3["alpha"].values
    err_b = (chk["est_b"] - chk["true_b"]).abs().max()
    err_s = (chk["est_s"] - chk["true_s"]).abs().max()
    err_h = (chk["est_h"] - chk["true_h"]).abs().max()
    print("Max abs recovery error across 25 portfolios:")
    print(f"  market b: {err_b:.4f}   SMB s: {err_s:.4f}   HML h: {err_h:.4f}")
    print(f"  mean |est_alpha|: {chk['est_alpha'].abs().mean():.4f} %/mo "
          f"(planted alphas were ~0, except small-growth)")
    assert err_b < 0.10 and err_s < 0.10 and err_h < 0.10, "loadings not recovered!"
    print("\nLoadings recovered to within 0.10 of truth -> the machinery works.")
else:
    print("Real data path: no planted ground truth to compare against (this is expected).")

The contrast: CAPM (market only)

The guide’s exercise 2 asks us to re-run the 25 regressions with only the market factor — the old CAPM — and watch two columns move. Dropping SMB and HML should (i) lower the R^2s, because size and value co-movement is now unexplained, and (ii) inflate the intercepts, because the average returns that SMB/HML were soaking up now leak into \alpha. That side-by-side is the cleanest single argument for why two extra factors earned their place.

Xc = sm.add_constant(factors[["Mkt-RF"]])  # CAPM design: const + market only

capm_rows = []
resid_capm = np.empty((len(excess), len(port_names)))
for j, nm in enumerate(port_names):
    res = sm.OLS(excess[nm].values, Xc.values).fit()
    capm_rows.append({"portfolio": nm, "alpha": res.params[0],
                      "t(alpha)": res.tvalues[0], "b_mkt": res.params[1],
                      "R2": res.rsquared})
    resid_capm[:, j] = res.resid
capm = pd.DataFrame(capm_rows).set_index("portfolio")

compare = pd.DataFrame({
    "CAPM mean R2":  [capm["R2"].mean()],
    "FF3  mean R2":  [ff3["R2"].mean()],
    "CAPM mean|alpha|": [capm["alpha"].abs().mean()],
    "FF3  mean|alpha|": [ff3["alpha"].abs().mean()],
    "CAPM #|t|>2": [int((capm["t(alpha)"].abs() > 2).sum())],
    "FF3  #|t|>2": [int((ff3["t(alpha)"].abs() > 2).sum())],
})
print("CAPM vs. three-factor model, averaged over the 25 portfolios")
print(compare.T)

The GRS test, by hand

Twenty-five separate t-tests cannot answer “are all the alphas jointly zero?” — that is the multiple-testing trap. Gibbons, Ross & Shanken (1989) give the right instrument: a single F-statistic for the joint null H_0:\alpha_1=\cdots=\alpha_N=0 that correctly accounts for the fact that the 25 residual series are correlated. We name it in the guide; here we build it, because seeing the formula demystifies it.

With T months, N test portfolios, and K factors, collect the N intercepts into a vector \hat{\boldsymbol\alpha}, the T\times N residuals into \hat{\mathbf E}, and the T\times K factor matrix into \mathbf F. Define the residual covariance \hat{\boldsymbol\Sigma}=\hat{\mathbf E}'\hat{\mathbf E}/T and the factor covariance \hat{\boldsymbol\Omega}=(\mathbf F-\bar{\mathbf F})'(\mathbf F-\bar{\mathbf F})/T with factor means \bar{\boldsymbol\mu}. The statistic is

\text{GRS} \;=\; \frac{T-N-K}{N}\, \Big(1+\bar{\boldsymbol\mu}'\,\hat{\boldsymbol\Omega}^{-1}\,\bar{\boldsymbol\mu}\Big)^{-1} \,\hat{\boldsymbol\alpha}'\,\hat{\boldsymbol\Sigma}^{-1}\,\hat{\boldsymbol\alpha} \;\sim\; F_{N,\;T-N-K}\quad\text{under } H_0.

Read the pieces. The core quadratic form \hat{\boldsymbol\alpha}'\hat{\boldsymbol\Sigma}^{-1} \hat{\boldsymbol\alpha} is the total pricing error, but weighted by the inverse residual covariance — an alpha in a portfolio whose residual barely moves counts for more than the same alpha in a noisy portfolio, and correlations between portfolios are netted out. The (1+\bar{\boldsymbol\mu}'\hat{\boldsymbol\Omega}^{-1}\bar{\boldsymbol\mu})^{-1} term is a Sharpe-ratio penalty for how large the factor premia themselves are. The leading (T-N-K)/N turns the whole thing into something that follows an F distribution with N and T-N-K degrees of freedom. A model passes the asset-pricing bar when GRS fails to reject — when the pricing errors are jointly indistinguishable from zero.

def grs_test(alphas, resid, factor_df):
    "GRS (1989) F-stat for H0: all intercepts = 0. resid is T x N, factor_df is T x K."
    alphas = np.asarray(alphas, float).reshape(-1)
    E = np.asarray(resid, float)
    F = np.asarray(factor_df, float)
    T, N = E.shape
    K = F.shape[1]
    Sigma = (E.T @ E) / T                      # N x N residual covariance (MLE divisor T)
    mu = F.mean(axis=0)                        # K factor means
    Fc = F - mu
    Omega = (Fc.T @ Fc) / T                    # K x K factor covariance
    sharpe_pen = 1.0 + mu @ np.linalg.solve(Omega, mu)
    quad = alphas @ np.linalg.solve(Sigma, alphas)
    stat = (T - N - K) / N * quad / sharpe_pen
    df1, df2 = N, T - N - K
    pval = stats.f.sf(stat, df1, df2)          # right tail
    return stat, pval, df1, df2

T_obs = len(excess)
# FF3: K = 3 factors
grs3, p3, d1_3, d2_3 = grs_test(ff3["alpha"].values, resid_mat,
                                factors[["Mkt-RF","SMB","HML"]].values)
# CAPM: K = 1 factor
grs1, p1, d1_1, d2_1 = grs_test(capm["alpha"].values, resid_capm,
                                factors[["Mkt-RF"]].values)

print(f"T = {T_obs} months, N = {len(port_names)} portfolios\n")
print(f"CAPM   : GRS F = {grs1:8.3f}  ~ F({d1_1},{d2_1})   p-value = {p1:.4g}")
print(f"FF 3-f : GRS F = {grs3:8.3f}  ~ F({d1_3},{d2_3})   p-value = {p3:.4g}")
print("\nReject H0 (all alphas = 0) at 5%?  ",
      f"CAPM: {'YES' if p1 < 0.05 else 'no'}   FF3: {'YES' if p3 < 0.05 else 'no'}")

Reading the GRS verdict

The number to compare is the p-value against 0.05. The three-factor model should land with a dramatically smaller GRS statistic than CAPM — its pricing errors are jointly far closer to zero. Whether FF3 fully clears the 5% bar depends on the sample and the construction; FF93’s honest story is that the model gets the intercepts down “close enough” to crush CAPM even when the joint test is not a spotless pass, with small-growth the usual culprit dragging the statistic up. In our synthetic world the planted alphas are essentially zero except the one small-growth cell, so you may see FF3 reject narrowly or not at all depending on how heavy that single corner sits — which is precisely the real lesson: one rebellious cell can move a joint test, and that is what the factor-zoo debates of the next decade fought over.

Let us put GRS into a picture: the inverse-covariance-weighted pricing error per model.

fig, ax = plt.subplots(figsize=(7, 4.5))
models = ["CAPM\n(market only)", "FF3\n(Mkt-RF, SMB, HML)"]
stats_vals = [grs1, grs3]
bars = ax.bar(models, stats_vals, color=["#b03a2e", "#1f618d"])
# critical value at 5% for the FF3 dof, drawn as the "pass/fail" line
crit = stats.f.ppf(0.95, d1_3, d2_3)
ax.axhline(crit, ls="--", color="gray",
           label=f"5% critical F ≈ {crit:.2f}")
for bar, v, p in zip(bars, stats_vals, [p1, p3]):
    ax.text(bar.get_x() + bar.get_width()/2, v, f"  F={v:.2f}\n  p={p:.3g}",
            ha="center", va="bottom", fontsize=9)
ax.set_ylabel("GRS F-statistic (joint test all α = 0)")
ax.set_title("Joint pricing-error test: lower is better")
ax.legend()
fig.tight_layout()
fig.savefig("grs_comparison.png", dpi=110)
plt.close(fig)
print("Figure saved to grs_comparison.png")
print(f"FF3 GRS = {grs3:.2f} vs CAPM GRS = {grs1:.2f}; 5% critical F ≈ {crit:.2f}")

What we built, and the two-column literacy to keep

We ran the FF93 machine end to end on whichever data loaded (real Ken French if the free pull worked, seeded synthetic otherwise) and read it the professional way:

R^2 near 1 across all 25 portfolios — the three factors capture the common variation, the part of returns many stocks share.
Loadings that march with the sorts — SMB falling down the size rows, HML rising across the value columns — the model recognizing the characteristics the portfolios were built on.
Intercepts small and mostly insignificant, with small-growth the stubborn corner, far better than CAPM’s bloated alphas.
A single GRS F-statistic that judges all 25 pricing errors jointly, dodging the multiple-testing trap, and that collapses for FF3 relative to CAPM.

High R^2 + zero \alpha = a working factor model. That is the ruler the rest of empirical asset pricing measures with: every Week-6 “anomaly” gets its three-factor alpha checked first.

Your Turn

Bolt on a momentum factor. If you used the real path, also pull "F-F_Momentum_Factor" (or synthesize a fourth factor WML and plant a small loading on it). Add it to Xf, re-run the 25 regressions, and re-compute GRS with K=4. Does the joint test improve? Does the small-growth alpha shrink — and what would it mean for the “risk vs. mispricing” debate if a fourth factor fixed the one cell FF3 cannot?
Stress-test by sample split. Cut the months in half (first vs. second period), re-run everything, and compare loadings and GRS across the two halves. Loadings should be stable; wild swings would warn that your “factor” is period-specific noise — the §6 robustness worry turned into a number.
Change one construction choice (synthetic path). In the DGP, shrink the planted alphas to exactly zero everywhere (including small-growth) and re-run. Confirm GRS now fails to reject for FF3 — the clean pass — so you can feel exactly how much that one rebellious corner was worth in the joint statistic.