Did the five-factor model survive the years Fama and French never saw?

# Frozen before any data was loaded. Every cell below reads from PREREG, so no
# value is typed twice and the plan cannot drift from the code.

PREREG = {
    "hypothesis":
        "The five-factor model of Fama and French (2015) holds up in the 150 "
        "months after that paper's sample ends. Specifically: RMW and CMA "
        "still earn positive premia; HML is still redundant once RMW and CMA "
        "are present; the five-factor model still beats the three-factor "
        "model on the paper's own test portfolios; and the failure the paper "
        "reports on small stocks with low profitability and high investment "
        "is no worse than it was.",

    "universe":
        "Ken French's US 2x3 five-factor series and his 25 Size-B/M, "
        "25 Size-OP, 25 Size-Inv and 32 Size-OP-Inv portfolios. The authors' "
        "own published returns, not a universe rebuilt here.",

    "start": "1963-07-01",
    "end":   "2026-06-30",

    "replication_window":
        "1963-07 to 2013-12, 606 months. The paper's own sample. Used once, "
        "to confirm this pull reproduces its published Tables 4 and 5.",

    "holdout":
        "2014-01 to 2026-06, 150 months. Every month the paper could not "
        "have seen. Touched once, in the last cell of Results.",

    "rebalance":
        "None. French's factors and portfolios are already monthly returns.",

    "params": {
        "as_of": "2026-09-01",
        "factor_file": "F-F_Research_Data_5_Factors_2x3",
        "test_portfolios": {
            "25 Size-B/M":     "25_Portfolios_5x5",
            "25 Size-OP":      "25_Portfolios_ME_OP_5x5",
            "25 Size-Inv":     "25_Portfolios_ME_INV_5x5",
            "32 Size-OP-Inv":  "32_Portfolios_ME_OP_INV_2x4x4",
        },
        "models": {
            "HML":         ["Mkt-RF", "SMB", "HML"],
            "HML RMW":     ["Mkt-RF", "SMB", "HML", "RMW"],
            "HML CMA":     ["Mkt-RF", "SMB", "HML", "CMA"],
            "RMW CMA":     ["Mkt-RF", "SMB", "RMW", "CMA"],
            "HML RMW CMA": ["Mkt-RF", "SMB", "HML", "RMW", "CMA"],
        },
        "t_method":
            "Plain OLS t-statistic on the mean, matching the paper. Table 4's "
            "t-statistics reproduce as mean/(sd/sqrt(n)) to rounding.",
        "nw_lags_robustness": 6,      # author's choice, not the paper's
        "alpha_t_threshold": 2.0,
        "subperiods": {              # robustness only, badly underpowered
            "2014-2019": ("2014-01-01", "2019-12-31"),
            "2020-2026": ("2020-01-01", "2026-06-30"),
        },
    },

    # Table 4, 2x3 factors, 1963-07 to 2013-12, 606 months.
    "paper_table4_2x3": {
        "Mkt-RF": {"mean": 0.50, "sd": 4.49, "t": 2.74},
        "SMB":    {"mean": 0.29, "sd": 3.07, "t": 2.31},
        "HML":    {"mean": 0.37, "sd": 2.88, "t": 3.20},
        "RMW":    {"mean": 0.25, "sd": 2.14, "t": 2.92},
        "CMA":    {"mean": 0.33, "sd": 2.01, "t": 4.07},
    },

    # Table 5, 2x3 factors, 25 Size-B/M portfolios (Panel A). Models are named
    # by the factors beyond Mkt-RF and SMB that they include.
    "paper_table5_2x3_size_bm": {
        "HML":          {"grs": 3.62, "avg_abs_alpha": 0.102},
        "HML RMW":      {"grs": 3.13, "avg_abs_alpha": 0.095},
        "HML CMA":      {"grs": 3.52, "avg_abs_alpha": 0.101},
        "RMW CMA":      {"grs": 2.84, "avg_abs_alpha": 0.100},
        "HML RMW CMA":  {"grs": 2.84, "avg_abs_alpha": 0.094},
    },

    "not_tested":
        "The abstract's claim that performance is insensitive to factor "
        "definition. French publishes only the 2x3 US construction; the 2x2 "
        "and 2x2x2x2 versions in Table 4 are not available, and rebuilding "
        "them would substitute this author's choices for the paper's.",

    "success":
        "Decided per test, before the holdout is opened. "
        "PREMIA: CONTINUED if the holdout 95% confidence interval for a "
        "factor's mean contains the paper's published mean; CHANGED if it "
        "excludes it. "
        "REDUNDANCY: a factor is redundant if the intercept from regressing "
        "it on the other four has |t| below 2.0. "
        "GRS: the five-factor model beats the three-factor model if its GRS "
        "statistic is lower on at least three of the four portfolio sets. "
        "HML is still redundant if dropping it changes GRS by less than the "
        "gap the paper reports between its three- and five-factor models. "
        "SMALL GROWTH: the named failure has worsened if the intercept on "
        "the small, low-profitability, high-investment portfolio is more "
        "negative out of sample than in the replication window. "
        "Significance of any single premium is not the criterion. Minimum "
        "detectable means and implied power are computed from the paper's "
        "own standard deviations and printed before any holdout result.",
}

import sys
import warnings

# Library import notices are not findings, and would otherwise be
# published as output. Silenced before the imports that raise them.
warnings.filterwarnings("ignore")


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

sys.path.append("..")
from lib import style
from lib.data import famafrench, to_month_end

style.use()
SEED = 0
np.random.seed(SEED)

AS_OF = PREREG["params"]["as_of"]

factors = to_month_end(famafrench(PREREG["params"]["factor_file"], AS_OF, 0))
portfolios = {
    label: to_month_end(famafrench(name, AS_OF, 0))
    for label, name in PREREG["params"]["test_portfolios"].items()
}

REP = slice(PREREG["start"], "2013-12-31")     # the paper's own months
OOS = slice("2014-01-01", PREREG["end"])       # months it could not have seen

print(f"factors      {factors.index.min():%Y-%m} to {factors.index.max():%Y-%m}"
      f"  n={len(factors)}  cols={list(factors.columns)}")
for label, df in portfolios.items():
    print(f"{label:16} {df.index.min():%Y-%m} to {df.index.max():%Y-%m}"
          f"  n={len(df)}  portfolios={df.shape[1]}")
print()
print(f"replication  n={len(factors.loc[REP])} months")
print(f"holdout      n={len(factors.loc[OOS])} months")

def t_stat(x):
    """Plain t-statistic on a mean, matching the paper's Table 4."""
    x = np.asarray(x, dtype=float)
    return x.mean() / (x.std(ddof=1) / np.sqrt(len(x)))


def newey_west_t(x, lags):
    """Autocorrelation-robust t on a mean. Robustness only, not the primary test."""
    x = np.asarray(x, dtype=float)
    n = len(x)
    e = x - x.mean()
    gamma0 = (e @ e) / n
    var = gamma0
    for k in range(1, lags + 1):
        gk = (e[k:] @ e[:-k]) / n
        var += 2 * (1 - k / (lags + 1)) * gk
    return x.mean() / np.sqrt(var / n)


def grs(excess, F):
    """Gibbons, Ross & Shanken (1989).

    `excess` is T x N test-asset excess returns, `F` is T x K factor returns.
    Returns the statistic, its p-value, and the average absolute intercept.
    """
    R = np.asarray(excess, dtype=float)
    Fv = np.asarray(F, dtype=float)
    T, N = R.shape
    K = Fv.shape[1]
    X = np.column_stack([np.ones(T), Fv])
    B = np.linalg.lstsq(X, R, rcond=None)[0]
    alpha = B[0]
    resid = R - X @ B
    Sigma = (resid.T @ resid) / (T - K - 1)
    mu = Fv.mean(axis=0)
    Omega = np.cov(Fv, rowvar=False, ddof=1)
    sharpe_sq = mu @ np.linalg.solve(np.atleast_2d(Omega), mu)
    stat = ((T - N - K) / N) * (alpha @ np.linalg.solve(Sigma, alpha)) / (1 + sharpe_sq)
    p = stats.f.sf(stat, N, T - N - K)
    return stat, p, np.abs(alpha).mean()


def detectable(sd, n, threshold):
    """Smallest mean this many months could call significant, and the power
    to detect a premium that simply continued at the paper's published size."""
    se = sd / np.sqrt(n)
    return se * threshold, se

# What 150 months can and cannot answer, computed from the paper's own
# standard deviations before any holdout number is read.
n_oos = len(factors.loc[OOS])
thr = PREREG["params"]["alpha_t_threshold"]

rows = []
for name, pub in PREREG["paper_table4_2x3"].items():
    mde, se = detectable(pub["sd"], n_oos, thr)
    t_if_unchanged = pub["mean"] / se
    power = stats.norm.sf(thr - t_if_unchanged)
    rows.append({
        "factor": name,
        "paper mean": pub["mean"],
        "paper sd": pub["sd"],
        "min detectable mean": round(mde, 3),
        "t if unchanged": round(t_if_unchanged, 2),
        "power at |t|>2": round(power, 2),
    })

power_table = pd.DataFrame(rows).set_index("factor")
print(f"holdout n = {n_oos} months\n")
print(power_table.to_string())

# Does today's vintage of French's files still produce the paper's Table 4?
rep = factors.loc[REP]
rows = []
for name, pub in PREREG["paper_table4_2x3"].items():
    x = rep[name]
    rows.append({
        "factor": name,
        "mean": round(x.mean(), 2),      "paper mean": pub["mean"],
        "sd": round(x.std(ddof=1), 2),   "paper sd": pub["sd"],
        "t": round(t_stat(x), 2),        "paper t": pub["t"],
    })
replication = pd.DataFrame(rows).set_index("factor")
print(f"1963-07 to 2013-12, n = {len(rep)} (paper: 606)\n")
print(replication.to_string())

# The same check for Table 5: model comparisons on the 25 Size-B/M portfolios.
rep_f = factors.loc[REP]
rep_p = portfolios["25 Size-B/M"].loc[REP]
rep_ex = rep_p.sub(rep_f["RF"], axis=0)

rows = []
for label, cols in PREREG["params"]["models"].items():
    g, p, aa = grs(rep_ex, rep_f[cols])
    pub = PREREG["paper_table5_2x3_size_bm"][label]
    rows.append({
        "model": label,
        "GRS": round(g, 2),        "paper GRS": pub["grs"],
        "A|a|": round(aa, 3),      "paper A|a|": pub["avg_abs_alpha"],
    })
replication_grs = pd.DataFrame(rows).set_index("model")
print(f"25 Size-B/M, 1963-07 to 2013-12, n = {len(rep_ex)}\n")
print(replication_grs.to_string())

# The holdout is opened here and read once.
oos_f = factors.loc[OOS]

rows = []
for name, pub in PREREG["paper_table4_2x3"].items():
    x = oos_f[name]
    m, sd, n = x.mean(), x.std(ddof=1), len(x)
    se = sd / np.sqrt(n)
    lo, hi = m - 1.96 * se, m + 1.96 * se
    contains = lo <= pub["mean"] <= hi
    rows.append({
        "factor": name,
        "paper mean": pub["mean"],
        "holdout mean": round(m, 2),
        "95% low": round(lo, 2),
        "95% high": round(hi, 2),
        "t": round(t_stat(x), 2),
        "verdict": "CONTINUED" if contains else "CHANGED",
    })
premia = pd.DataFrame(rows).set_index("factor")
print(f"2014-01 to 2026-06, n = {len(oos_f)} months\n")
print(premia.to_string())

FACTORS = list(PREREG["paper_table4_2x3"])

def redundancy(window_f, label):
    rows = []
    for target in FACTORS:
        others = [c for c in FACTORS if c != target]
        y = window_f[target].values
        X = np.column_stack([np.ones(len(y)), window_f[others].values])
        beta, *_ = np.linalg.lstsq(X, y, rcond=None)
        resid = y - X @ beta
        dof = len(y) - X.shape[1]
        s2 = (resid @ resid) / dof
        se = np.sqrt(s2 * np.linalg.inv(X.T @ X)[0, 0])
        t = beta[0] / se
        rows.append({
            "factor": target,
            "intercept": round(beta[0], 3),
            "t": round(t, 2),
            "redundant": abs(t) < PREREG["params"]["alpha_t_threshold"],
        })
    return pd.DataFrame(rows).set_index("factor").add_suffix(f" [{label}]")

red_rep = redundancy(factors.loc[REP], "1963-2013")
red_oos = redundancy(oos_f, "2014-2026")
print("Each factor regressed on the other four.\n")
print(red_rep.join(red_oos).to_string())

def grs_table(window):
    out = {}
    for pf_label, df in portfolios.items():
        f_w = factors.loc[window]
        ex = df.loc[window].sub(f_w["RF"], axis=0)
        ex = ex.replace(-99.99, np.nan).dropna(axis=1, how="any")
        col = {}
        for m_label, cols in PREREG["params"]["models"].items():
            g, p, aa = grs(ex, f_w[cols])
            col[m_label] = round(g, 2)
        out[f"{pf_label} (N={ex.shape[1]})"] = col
    return pd.DataFrame(out)

grs_rep = grs_table(REP)
grs_oos = grs_table(OOS)
print("GRS statistic, lower is better. Replication window:\n")
print(grs_rep.to_string())
print("\nHoldout:\n")
print(grs_oos.to_string())

# The failure the paper names: small stocks, low profitability, high investment.
# In the 2x4x4 sort that is the small size group, lowest OP quartile, highest
# Inv quartile. French labels these "SMALL LoOP HiINV" or a close variant.
sm = portfolios["32 Size-OP-Inv"]
corner = [c for c in sm.columns if c.strip().upper().startswith("SMALL")
          and "LOOP" in c.replace(" ", "").upper()
          and "HIINV" in c.replace(" ", "").upper()]
print("candidate columns:", corner if corner else list(sm.columns[:8]))

five = PREREG["params"]["models"]["HML RMW CMA"]

def corner_alpha(window, col):
    f_w = factors.loc[window]
    y = (sm.loc[window, col] - f_w["RF"]).values
    X = np.column_stack([np.ones(len(y)), f_w[five].values])
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    resid = y - X @ beta
    dof = len(y) - X.shape[1]
    se = np.sqrt((resid @ resid) / dof * np.linalg.inv(X.T @ X)[0, 0])
    return beta[0], beta[0] / se

col = corner[0]
a_rep, t_rep = corner_alpha(REP, col)
a_oos, t_oos = corner_alpha(OOS, col)
print(f"\nportfolio: {col!r}")
print(f"five-factor intercept, 1963-2013 : {a_rep:+.3f} %/month  (t = {t_rep:+.2f})")
print(f"five-factor intercept, 2014-2026 : {a_oos:+.3f} %/month  (t = {t_oos:+.2f})")
print(f"worsened: {a_oos < a_rep}")

# The holdout only. Plotting from 1963 compresses the 150 months this note is
# about into a sliver, and the question is not what the factors did over sixty
# years, it is whether they kept paying after the paper went out. Dotted lines
# are what the paper's published means would have compounded to over the same
# months, so the gap between solid and dotted is the finding.
fig, ax = plt.subplots(figsize=(8.0, 4.2))

oos_cum = (1 + factors.loc[OOS, ["RMW", "CMA", "HML"]] / 100).cumprod()
months = np.arange(1, len(oos_cum) + 1)

# Colours are set per factor rather than through `style.emphasise`, which walks
# every line on the axes and would recolour each dotted reference away from the
# solid line it belongs to. A reader has to be able to pair them by eye.
PAIRS = [("RMW", style.ACCENT, 1.8), ("CMA", style.SERIES[1], 1.0),
         ("HML", style.SERIES[2], 1.0)]

for name, colour, width in PAIRS:
    ax.plot(oos_cum.index, oos_cum[name], color=colour, lw=width,
            zorder=3 if colour == style.ACCENT else 2, label=name)
    implied = (1 + PREREG["paper_table4_2x3"][name]["mean"] / 100) ** months
    ax.plot(oos_cum.index, implied, color=colour, lw=0.9, ls=(0, (2, 2)),
            alpha=0.55, zorder=1)

ax.axhline(1.0, color=style.MUTED, lw=0.6, alpha=0.5)
ax.set_ylabel("growth of 1, from January 2014")
ax.set_title("What the three factors paid after the paper was published")
ax.legend(loc="upper left", title="solid: actual   dotted: paper's mean")
style.check_points(ax)
plt.show()

# One portfolio set, the paper's own Panel A, so the bars line up with a
# published number. The dashed line is where GRS stops being distinguishable
# from zero intercepts at 5%; a bar below it means nothing is left to explain.
fig, ax = plt.subplots(figsize=(8.0, 4.2))

col_rep = grs_rep["25 Size-B/M (N=25)"]
col_oos = grs_oos["25 Size-B/M (N=25)"]
labels = list(col_oos.index)
x = np.arange(len(labels))
w = 0.38

ax.bar(x - w / 2, col_rep.values, w, color=style.SERIES[2], label="1963-2013")
ax.bar(x + w / 2, col_oos.values, w, color=style.ACCENT, label="2014-2026")

crit = stats.f.ppf(0.95, 25, len(factors.loc[OOS]) - 25 - 5)
ax.axhline(crit, color=style.MUTED, lw=0.9, ls=(0, (4, 3)))
ax.annotate(f"5% critical value ({crit:.2f})", xy=(-0.45, crit),
            xytext=(0, 5), textcoords="offset points", ha="left",
            fontsize=8, color=style.MUTED)

ax.set_xticks(x)
ax.set_xticklabels(labels, rotation=20, ha="right")
ax.set_ylabel("GRS, 25 Size-B/M portfolios")
ax.set_title("Does the model leave less unexplained than it used to?")
# Headroom first, then a two-column legend above the bars. At the default
# limits the legend box sat on top of the tallest bar.
ax.set_ylim(0, 4.4)
ax.legend(loc="upper center", ncol=2)
plt.show()

rows = []
lags = PREREG["params"]["nw_lags_robustness"]
for name in FACTORS:
    x = oos_f[name]
    rows.append({
        "factor": name,
        "plain t": round(t_stat(x), 2),
        f"Newey-West t ({lags} lags)": round(newey_west_t(x, lags), 2),
    })
print("Autocorrelation correction. Not the primary test; the paper does not use one.\n")
print(pd.DataFrame(rows).set_index("factor").to_string())

# Splitting 150 months in half leaves each piece badly underpowered. This is
# here because it was preregistered, not because either half decides anything.
rows = []
for label, (a, b) in PREREG["params"]["subperiods"].items():
    w = factors.loc[a:b]
    for name in FACTORS:
        rows.append({
            "period": label, "n": len(w), "factor": name,
            "mean": round(w[name].mean(), 2),
            "t": round(t_stat(w[name]), 2),
        })
sub = pd.DataFrame(rows).pivot(index="factor", columns="period",
                               values=["mean", "t"])
print(sub.to_string())

# End-point sensitivity. The holdout stops in the middle of a sharp RMW
# drawdown: April, May and June 2026 ran -4.30, -8.42 and -4.64. A mean
# measured to the last available month inherits that, so here is the same mean
# with the tail trimmed. This was not preregistered; it is here because the
# data showed it.
rows = []
for cut in [0, 3, 6, 12]:
    w = oos_f.iloc[: len(oos_f) - cut] if cut else oos_f
    row = {"months dropped": cut, "n": len(w)}
    for name in FACTORS:
        row[name] = round(w[name].mean(), 2)
    rows.append(row)
print("Holdout mean with the final months removed.\n")
print(pd.DataFrame(rows).set_index("months dropped").to_string())


# Vintage drift: how far today's files have moved from the paper's 2014 pull.
drift = replication.copy()
drift["mean gap"] = (drift["mean"] - drift["paper mean"]).round(2)
drift["sd gap"] = (drift["sd"] - drift["paper sd"]).round(2)
drift["sd gap %"] = (100 * drift["sd gap"] / drift["paper sd"]).round(1)
print("Same months, same construction, twelve years of rebuilds apart.\n")
print(drift[["mean gap", "sd gap", "sd gap %"]].to_string())

import matplotlib, scipy
print(f"python      {sys.version.split()[0]}")
print(f"pandas      {pd.__version__}")
print(f"numpy       {np.__version__}")
print(f"scipy       {scipy.__version__}")
print(f"matplotlib  {matplotlib.__version__}")
print(f"seed        {SEED}")
print(f"data as-of  famafrench {AS_OF}")