Is the security market line still flat? Testing Frazzini and Pedersen out of sample

# Frozen before any data is loaded. Every cell below reads from PREREG,
# so no value is typed twice and the plan cannot drift from the code.
#
# Provenance of every parameter, because a number attributed to a paper
# nobody opened is a fabricated citation wearing a number.
#
# From Frazzini & Pedersen (2014), JFE 111(1), 1-25, read out of the
# author copy at w4.stern.nyu.edu:
#   - the prediction tested: "the alphas decline almost monotonically from
#     the low-beta to high-beta portfolios" (their Proposition 1, Table III)
#   - the BAB construction: each leg rescaled to a beta of one at formation,
#     long low-beta, short high-beta, self-financing at the risk-free rate
#     (their Eq. 17)
#   - their US sample ends March 2012, which is where the holdout starts
#
# Mine, not theirs, and named as mine:
#   - post_rank_beta_window_months. Frazzini and Pedersen estimate stock
#     betas from a 1-year volatility window and a 5-year correlation window
#     and shrink 0.6 of the way toward 1. That is a stock-level device and
#     the sort here is already done by French, so the rescaling beta is a
#     plain 60-month rolling regression of the decile on the market. Their
#     shrinkage is deliberately not applied.
#   - cost_bps and hac_lags.

PREREG = {
    "hypothesis": (
        "CAPM alphas fall as pre-ranking beta rises across Ken French's "
        "beta-sorted deciles; a beta-neutral portfolio long the low-beta "
        "decile and short the high-beta decile, each rescaled to a beta of "
        "one, earns positive three-factor alpha net of costs; and that "
        "portfolio has a positive Sharpe ratio in the fourteen years after "
        "Frazzini and Pedersen's sample ends."
    ),
    "universe": (
        "Ken French Portfolios_Formed_on_BETA, value-weighted deciles built "
        "from NYSE, AMEX and NASDAQ common stocks on CRSP, sorted each June "
        "on Scholes-Williams betas from the prior 60 months. French builds "
        "the portfolios from CRSP himself, so no survivorship bias enters "
        "through universe selection. The sort is not Frazzini and Pedersen's "
        "beta estimator, which is the main gap between this test and theirs."
    ),
    "start":     "1963-07-31",
    "end":       "2026-06-30",
    "holdout":   "2012-04-30 onwards, touched once in the last cell of Results",
    "rebalance": "monthly, from French's monthly portfolio returns",
    "params": {
        "n_deciles":                    10,
        "post_rank_beta_window_months": 60,
        "post_rank_beta_min_months":    36,
        "cost_bps":                     10,
        "hac_lags":                     6,
    },
    "success": (
        "All three must hold. (a) In sample, the Spearman rank correlation "
        "between decile number and CAPM alpha is negative with p < 0.05. "
        "(b) In sample, the rescaled long-short portfolio has a positive "
        "Fama-French three-factor alpha with |t| > 2 under Newey-West "
        "standard errors, net of 10bp round-trip costs. (c) In the holdout, "
        "that portfolio's net Sharpe ratio is positive."
    ),
}

import sys
import warnings

# Library import notices are not findings, and they 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 statsmodels.api as sm
from scipy import stats
import matplotlib.pyplot as plt

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

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

AS_OF = "2026-08-28"          # as-of date of every cached pull in this note
P = PREREG["params"]
DECILES = ["Lo 10"] + [f"Dec {i}" for i in range(2, 10)] + ["Hi 10"]

dec = data.to_month_end(data.famafrench("Portfolios_Formed_on_BETA", AS_OF, table=0))
ff  = data.to_month_end(data.famafrench("F-F_Research_Data_Factors", AS_OF))

m = dec[DECILES].join(ff[["Mkt-RF", "RF", "SMB", "HML"]], how="inner")
m = m.loc[PREREG["start"]:PREREG["end"]]
ex = m[DECILES].sub(m["RF"], axis=0)          # excess returns, percent per month

ins = ex.loc[:"2012-03-31"]                    # paper's sample ends March 2012
oos = ex.loc["2012-04-30":]                    # holdout, touched once at the end

print(f"deciles   {len(dec):6,} rows  {dec.index.min().date()} -> {dec.index.max().date()}")
print(f"factors   {len(ff):6,} rows  {ff.index.min().date()} -> {ff.index.max().date()}")
print(f"merged    {len(m):6,} rows  {m.index.min().date()} -> {m.index.max().date()}")
print(f"in-sample {len(ins):6,} months     holdout {len(oos):,} months")
print(f"missing values in the merged panel: {int(m.isna().sum().sum())}")

def capm(y, x, lags=P["hac_lags"]):
    return sm.OLS(y, sm.add_constant(x)).fit(cov_type="HAC", cov_kwds={"maxlags": lags})

def rolling_beta(y, x, window=P["post_rank_beta_window_months"],
                 min_periods=P["post_rank_beta_min_months"]):
    """Beta of y on x from a trailing window, shifted so month t uses only
    data through t-1. Without the shift the portfolio would be using its own
    contemporaneous return to size itself."""
    cov = y.rolling(window, min_periods=min_periods).cov(x)
    var = x.rolling(window, min_periods=min_periods).var()
    return (cov / var).shift(1)

def bab_factor(ex_panel, mkt, lo="Lo 10", hi="Hi 10", cost_bps=P["cost_bps"]):
    """Frazzini and Pedersen Eq. (17) at portfolio level: each leg rescaled to
    a beta of one, long low, short high, self-financing at the risk-free rate."""
    w_lo = 1.0 / rolling_beta(ex_panel[lo], mkt)
    w_hi = 1.0 / rolling_beta(ex_panel[hi], mkt)
    gross = w_lo * ex_panel[lo] - w_hi * ex_panel[hi]
    turnover = w_lo.diff().abs().fillna(0) + w_hi.diff().abs().fillna(0)
    return gross - turnover * cost_bps / 100.0, w_lo, w_hi, turnover

# Criterion (a): does alpha decline with beta? In sample only.

mkt_ins = m.loc[ins.index, "Mkt-RF"]
rows = []
for i, col in enumerate(DECILES, start=1):
    fit = capm(ins[col], mkt_ins)
    rows.append({
        "decile": i,
        "beta":   fit.params["Mkt-RF"],
        "mean":   ins[col].mean(),
        "alpha":  fit.params["const"],
        "t":      fit.tvalues["const"],
        "sharpe": np.sqrt(12) * ins[col].mean() / ins[col].std(),
    })
tab = pd.DataFrame(rows)

rho,  pval  = stats.spearmanr(tab["decile"], tab["alpha"])
rho2, pval2 = stats.spearmanr(tab["decile"], tab["sharpe"])

print("1963-07 to 2012-03, value-weighted beta deciles, monthly percent\n")
print(f"{'decile':>6} {'beta':>7} {'mean':>7} {'alpha':>7} {'t':>7} {'sharpe':>7}")
for r in rows:
    print(f"{r['decile']:>6} {r['beta']:7.2f} {r['mean']:7.3f} "
          f"{r['alpha']:+7.3f} {r['t']:+7.2f} {r['sharpe']:7.2f}")
print(f"\nSpearman(decile, alpha)  rho {rho:+.3f}   p {pval:.4f}")
print(f"Spearman(decile, Sharpe) rho {rho2:+.3f}   p {pval2:.4f}")

fig, ax = plt.subplots(figsize=(8, 4.4))

b, r = tab["beta"].values, tab["mean"].values
mkt_beta, mkt_mean = 1.0, mkt_ins.mean()

# The CAPM line: what the deciles would earn if beta were priced as theory says.
xs = np.linspace(0.5, 1.75, 50)
ax.plot(xs, xs * mkt_mean, color=style.SERIES[1], lw=1.3, ls="--",
        label="CAPM: return proportional to beta")

slope, intercept = np.polyfit(b, r, 1)
ax.plot(xs, intercept + slope * xs, color=style.ACCENT, lw=1.8,
        label=f"Fitted through the deciles (slope {slope:.2f}%)")
ax.scatter(b, r, s=26, color=style.ACCENT, zorder=4)
ax.scatter([mkt_beta], [mkt_mean], s=42, marker="s", color=style.TEXT, zorder=5)
ax.annotate("market", (mkt_beta, mkt_mean), textcoords="offset points",
            xytext=(8, -11), fontsize=8.5, color=style.MUTED)
for i, (bi, ri) in enumerate(zip(b, r), start=1):
    if i in (1, 10):
        ax.annotate(f"decile {i}", (bi, ri), textcoords="offset points",
                    xytext=(6, 6), fontsize=8.5, color=style.MUTED)

ax.set_xlabel("Market beta, in sample")
ax.set_ylabel("Mean monthly excess return, %")
ax.set_title("The security market line is far flatter than the CAPM predicts")
ax.legend(loc="upper left")
plt.show()

# Criterion (b): the rescaled long-short, in sample.

bab, w_lo, w_hi, turnover = bab_factor(ex, m["Mkt-RF"])
bab_ins = bab.loc[:"2012-03-31"].dropna()

X = m.loc[bab_ins.index, ["Mkt-RF", "SMB", "HML"]]
fit = capm(bab_ins, X)
sharpe_ins = np.sqrt(12) * bab_ins.mean() / bab_ins.std()

print(f"leverage on the low-beta leg    mean {w_lo.mean():.2f}   max {w_lo.max():.2f}")
print(f"scaling on the high-beta leg    mean {w_hi.mean():.2f}   max {w_hi.max():.2f}")
print(f"monthly weight turnover         mean {turnover.mean():.3f}\n")
print(f"in sample, net of {P['cost_bps']}bp   n {len(bab_ins)} months")
print(f"  mean            {bab_ins.mean():+.3f}% per month")
print(f"  Sharpe          {sharpe_ins:+.3f}")
print(f"  three-factor alpha {fit.params['const']:+.3f}%  t {fit.tvalues['const']:+.2f}")

fig, ax = plt.subplots(figsize=(8, 4.2))

growth = (1 + bab.dropna() / 100).cumprod()
ax.plot(growth.index, growth.values, color=style.ACCENT, lw=1.6)
ax.axvspan(pd.Timestamp("2012-04-30"), growth.index.max(),
           color=style.SUBTLE, alpha=0.14, lw=0)
ax.annotate("holdout", (pd.Timestamp("2013-06-30"), growth.min() * 1.4),
            fontsize=8.5, color=style.MUTED)
ax.set_yscale("log")
ax.set_ylabel("Growth of $1, log scale")
ax.set_title("Beta-neutral low-minus-high, net of costs")
plt.show()

# Criterion (c): the holdout. Evaluated once.

bab_oos = bab.loc["2012-04-30":].dropna()
X_oos = m.loc[bab_oos.index, ["Mkt-RF", "SMB", "HML"]]
fit_oos = capm(bab_oos, X_oos)
sharpe_oos = np.sqrt(12) * bab_oos.mean() / bab_oos.std()

print(f"holdout 2012-04 to 2026-06, net of {P['cost_bps']}bp   n {len(bab_oos)} months")
print(f"  mean            {bab_oos.mean():+.3f}% per month")
print(f"  Sharpe          {sharpe_oos:+.3f}")
print(f"  three-factor alpha {fit_oos.params['const']:+.3f}%  t {fit_oos.tvalues['const']:+.2f}")
print()
print(f"  in sample:  mean {bab_ins.mean():+.3f}%  Sharpe {sharpe_ins:+.3f}  "
      f"alpha {fit.params['const']:+.3f}% (t {fit.tvalues['const']:+.2f})")
print(f"  holdout:    mean {bab_oos.mean():+.3f}%  Sharpe {sharpe_oos:+.3f}  "
      f"alpha {fit_oos.params['const']:+.3f}% (t {fit_oos.tvalues['const']:+.2f})")

def variant(label, table=0, cols=None, window=P["post_rank_beta_window_months"],
            cost=P["cost_bps"]):
    cols = cols or DECILES
    d = data.to_month_end(data.famafrench("Portfolios_Formed_on_BETA", AS_OF, table=table))
    mm = d[cols].join(ff[["Mkt-RF", "RF", "SMB", "HML"]], how="inner")
    mm = mm.loc[PREREG["start"]:PREREG["end"]]
    e = mm[cols].sub(mm["RF"], axis=0)
    w_lo = 1.0 / rolling_beta(e[cols[0]],  mm["Mkt-RF"], window=window)
    w_hi = 1.0 / rolling_beta(e[cols[-1]], mm["Mkt-RF"], window=window)
    g = w_lo * e[cols[0]] - w_hi * e[cols[-1]]
    turn = w_lo.diff().abs().fillna(0) + w_hi.diff().abs().fillna(0)
    r = (g - turn * cost / 100.0).loc[:"2012-03-31"].dropna()
    f = capm(r, mm.loc[r.index, ["Mkt-RF", "SMB", "HML"]])
    return (label, len(r), r.mean(), np.sqrt(12) * r.mean() / r.std(),
            f.params["const"], f.tvalues["const"])

QUINTILES = ["Lo 20", "Qnt 2", "Qnt 3", "Qnt 4", "Hi 20"]
variants = [
    variant("preregistered: VW deciles, 60m beta, 10bp"),
    variant("equal-weighted deciles", table=1),
    variant("beta window 36 months", window=36),
    variant("beta window 120 months", window=120),
    variant("quintiles instead of deciles", cols=QUINTILES),
    variant("cost 0bp", cost=0),
    variant("cost 50bp", cost=50),
    variant("cost 100bp", cost=100),
]

print(f"{'variant':44} {'n':>4} {'mean%':>7} {'Sharpe':>7} {'alpha%':>7} {'t':>6}")
for lab, n, mu, sr, a, t in variants:
    print(f"{lab:44} {n:>4} {mu:+7.3f} {sr:+7.3f} {a:+7.3f} {t:+6.2f}")

# Subperiod behaviour, preregistered specification, net of costs.

print("decade      n    mean%   Sharpe")
for start in range(1960, 2030, 10):
    s = bab.loc[f"{start}":f"{start + 9}"].dropna()
    if len(s) >= 24:
        print(f"{start}s   {len(s):>3}  {s.mean():+7.3f}  {np.sqrt(12) * s.mean() / s.std():+7.2f}")

fig, ax = plt.subplots(figsize=(8, 3.9))

labels, means, counts = [], [], []
for start in range(1960, 2030, 10):
    s = bab.loc[f"{start}":f"{start + 9}"].dropna()
    if len(s) >= 24:
        labels.append(f"{start}s")
        means.append(s.mean())
        counts.append(len(s))

# Only the two decades the title is about are loud. Painting every positive
# bar in the accent would make the chart say "mostly positive", which is the
# opposite of the point.
top2 = sorted(range(len(means)), key=lambda i: means[i], reverse=True)[:2]
colors = [style.ACCENT if i in top2 else style.SERIES[2] for i in range(len(means))]

ax.bar(labels, means, color=colors, width=0.62)
ax.axhline(0, color=style.MUTED, lw=0.7)

# Two of these decades are part-period stubs. Without the counts the chart
# would present 42 months and 120 months as equally solid.
for i, (v, n) in enumerate(zip(means, counts)):
    ax.annotate(f"n={n}", (i, 0), textcoords="offset points",
                xytext=(0, -16 if v > 0 else 8), ha="center",
                fontsize=7.5, color=style.MUTED)

ax.set_ylabel("Mean monthly return, %")
ax.set_title("Two decades carry the premium")
ax.tick_params(axis="x", pad=14)
plt.show()

print(f"python       {sys.version.split()[0]}")
print(f"pandas       {pd.__version__}")
print(f"numpy        {np.__version__}")
print(f"statsmodels  {sm.__version__}")
print(f"seed         {SEED}")
print(f"data as-of   {AS_OF}  (research/data/famafrench/{AS_OF}/)")