CAPM — Single-Factor Model, Beta, and the Market Risk Premium

Companion notebook to the CAPM — Single-Factor Model article (Factor Models, Module 1). Uses numpy throughout; no other external dependencies.

§0 — Imports and simulation setup: define true parameters, seed RNG
§1 — Simulate market and asset returns; estimate β via OLS on excess returns
§2 — Variance decomposition: systematic (β²σₑ²) vs. idiosyncratic (σε²); verify R²
§3 — Jensen’s alpha: compute α̂, t-statistic, 95% confidence interval
§4 — Portfolio beta linearity: verify β̣ = weighted average of individual betas
§5 — Validation summary

Press Shift+Enter or click ▶ Run on each cell. Run cells in order — §2–§5 depend on §0 and §1.

§0 — Imports and Simulation Parameters

We use numpy for matrix operations and OLS. All true parameters are declared as named constants so results are reproducible and the ground truth is always accessible for validation. The simulation uses daily returns over T = 500 trading periods.

Convention: all returns are expressed as decimals (not percentages). Annualisation uses √252 for volatility and ×252 for expected returns. The risk-free rate is held constant at 3% per annum throughout.

import numpy as np

rng = np.random.default_rng(seed=42)

# ── True parameters (annualised) ──────────────────────────────────────────────
BETA_TRUE        = 1.35    # true systematic sensitivity
SIGMA_M_ANNUAL   = 0.15    # market annualised vol (15%)
SIGMA_EPS_ANNUAL = 0.20    # idiosyncratic annualised vol (20%)
MU_M_ANNUAL      = 0.08    # market annualised expected excess return (8%)
RF_ANNUAL        = 0.03    # risk-free rate per annum (3%)

# ── Convert to daily ──────────────────────────────────────────────────────────
TRADING_DAYS = 252
T            = 500          # number of simulated periods

sigma_M   = SIGMA_M_ANNUAL   / np.sqrt(TRADING_DAYS)
sigma_eps = SIGMA_EPS_ANNUAL / np.sqrt(TRADING_DAYS)
mu_M_daily = MU_M_ANNUAL / TRADING_DAYS
rf_daily   = RF_ANNUAL   / TRADING_DAYS

# ── Expected return of asset i under CAPM ────────────────────────────────────
mu_i_annual = RF_ANNUAL + BETA_TRUE * MU_M_ANNUAL  # SML equation
print(f"True parameters:")
print(f"  beta_true        = {BETA_TRUE:.2f}")
print(f"  sigma_M (annual) = {SIGMA_M_ANNUAL:.1%}")
print(f"  sigma_eps (ann.) = {SIGMA_EPS_ANNUAL:.1%}")
print(f"  E[R_M] (annual)  = {RF_ANNUAL + MU_M_ANNUAL:.1%}  (rf={RF_ANNUAL:.1%}, MRP={MU_M_ANNUAL:.1%})")
print(f"  E[R_i] (CAPM)    = {mu_i_annual:.2%}  (= rf + beta * MRP = {RF_ANNUAL:.2%} + {BETA_TRUE}*{MU_M_ANNUAL:.2%})")
print()
print(f"Simulation: T = {T} daily periods  (approx {T/TRADING_DAYS:.1f} years of trading data)")
print("✓ Parameters loaded.")

§1 — Simulate Returns and Estimate β via OLS

We generate market excess returns from a normal distribution with known parameters, then generate asset excess returns via the market model: r_i^e = alpha + beta * r_M^e + eps, with alpha = 0 (CAPM equilibrium).

OLS estimator: β̂ = Cov(r_i^e, r_M^e) / Var(r_M^e). We verify that the estimate is close to the true value of 1.35.

Key formula: β̂ = (r_M^e · r_i^e - T·μ̂_M·μ̂_i) / (|r_M^e|² - T·μ̂_M²) — numerically equivalent to np.cov with ddof=1.

# ── Simulate market excess returns ───────────────────────────────────────────
r_M_excess = rng.normal(loc=mu_M_daily, scale=sigma_M, size=T)

# ── Simulate idiosyncratic shocks ─────────────────────────────────────────────
eps = rng.normal(loc=0.0, scale=sigma_eps, size=T)

# ── Asset excess returns: market model with true alpha=0 ──────────────────────
alpha_true = 0.0
r_i_excess = alpha_true + BETA_TRUE * r_M_excess + eps

# ── OLS estimation ────────────────────────────────────────────────────────────
# beta_hat = Cov(r_i^e, r_M^e) / Var(r_M^e)
cov_matrix = np.cov(r_i_excess, r_M_excess, ddof=1)   # shape (2, 2)
beta_hat   = cov_matrix[0, 1] / cov_matrix[1, 1]       # Cov(i,M) / Var(M)
alpha_hat  = np.mean(r_i_excess) - beta_hat * np.mean(r_M_excess)

# ── OLS residuals ─────────────────────────────────────────────────────────────
eps_hat    = r_i_excess - (alpha_hat + beta_hat * r_M_excess)
sigma_eps_hat = np.std(eps_hat, ddof=2)                # ddof=2: two parameters estimated

print(f"OLS Beta Estimation Results (T = {T}):")
print(f"  beta_true  = {BETA_TRUE:.4f}")
print(f"  beta_hat   = {beta_hat:.4f}  (error = {beta_hat - BETA_TRUE:+.4f})")
print(f"  alpha_hat  = {alpha_hat:.6f}  (daily;  true = {alpha_true:.6f})")
print()

# Asymptotic standard error: SE(beta) = sigma_eps / (sigma_M * sqrt(T))
sigma_M_hat  = np.std(r_M_excess, ddof=1)
se_beta = sigma_eps_hat / (sigma_M_hat * np.sqrt(T))
print(f"  SE(beta_hat) = {se_beta:.4f}  (formula: sigma_eps / (sigma_M * sqrt(T)))")
print(f"  95% CI beta: [{beta_hat - 1.96*se_beta:.3f}, {beta_hat + 1.96*se_beta:.3f}]")
print(f"  True beta in 95% CI? {'YES' if abs(beta_hat - BETA_TRUE) < 1.96*se_beta else 'NO'}")
print()
print(f"Annualised sigma_eps_hat = {sigma_eps_hat * np.sqrt(TRADING_DAYS):.2%}  (true = {SIGMA_EPS_ANNUAL:.2%})")

§2 — Variance Decomposition: Systematic vs. Idiosyncratic

Total variance of excess returns decomposes as: Var(r_i^e) = β̂² × Var(r_M^e) + Var(ε̂). The R² of the regression equals the fraction of total variance explained by the market factor.

R² = systematic variance / total variance. For large-cap US equities, R² typically ranges from 0.20 to 0.80. Low R² indicates a stock with predominantly idiosyncratic risk — it is only weakly correlated with the market.

# ── Variance decomposition ────────────────────────────────────────────────────
var_total      = np.var(r_i_excess, ddof=1)
var_systematic = beta_hat**2 * np.var(r_M_excess, ddof=1)
var_idio       = np.var(eps_hat, ddof=2)

R2_decomp    = var_systematic / var_total          # from decomposition formula
R2_corr      = np.corrcoef(r_i_excess, r_M_excess)[0, 1]**2  # from correlation

# True values (population)
var_total_true      = BETA_TRUE**2 * sigma_M**2 + sigma_eps**2
var_systematic_true = BETA_TRUE**2 * sigma_M**2
var_idio_true       = sigma_eps**2
R2_true             = var_systematic_true / var_total_true

print("Variance Decomposition (daily units):")
print(f"  {'Component':<30} {'Estimated':>12} {'True (pop.)':>14} {'Ann. vol (est.)':>16}")
print(f"  {'-'*74}")
print(f"  {'Total variance':.<30} {var_total:>12.8f} {var_total_true:>14.8f} {np.sqrt(var_total * TRADING_DAYS):>15.2%}")
print(f"  {'Systematic (beta^2 * sigma_M^2)':.<30} {var_systematic:>12.8f} {var_systematic_true:>14.8f} {np.sqrt(var_systematic * TRADING_DAYS):>15.2%}")
print(f"  {'Idiosyncratic (Var(eps))':.<30} {var_idio:>12.8f} {var_idio_true:>14.8f} {np.sqrt(var_idio * TRADING_DAYS):>15.2%}")
print()
print(f"  R^2 (from decomposition):  {R2_decomp:.4f}")
print(f"  R^2 (from correlation^2):  {R2_corr:.4f}  (should match)")
print(f"  R^2 (true population):     {R2_true:.4f}")
print()

# Verify: systematic + idiosyncratic = total (up to rounding)
ok = abs((var_systematic + var_idio) - var_total) < 1e-12
print(f"  Check: systematic + idiosyncratic = total?  {'✓' if ok else '✗'}")
print(f"  Residual = {(var_systematic + var_idio) - var_total:.2e}")

§3 — Jensen’s Alpha: Estimation and Hypothesis Test

Jensen’s alpha is the OLS intercept in the market model regression. CAPM predicts α̂ = 0 in equilibrium. We test H₀: α = 0 using a standard t-test.

SE(α̂) = σ̂_ε × √(1/T + μ̂₂_M / Σ(ṙ_M)²) — the standard OLS formula for the intercept standard error.

Important: Because we simulated with true α = 0, we expect the t-statistic to be small. Rejecting H₀ here would indicate a simulation error, not a genuine alpha signal.

# ── Standard error of the OLS intercept ──────────────────────────────────────
# SE(alpha) = sigma_eps * sqrt(1/T + mean(r_M^e)^2 / sum((r_M^e - mean(r_M^e))^2))
r_M_demeaned_sq_sum = np.sum((r_M_excess - np.mean(r_M_excess))**2)
se_alpha = sigma_eps_hat * np.sqrt(1.0/T + np.mean(r_M_excess)**2 / r_M_demeaned_sq_sum)

# t-statistic under H0: alpha = 0
t_stat = alpha_hat / se_alpha

# t critical value at 5% two-tailed, df = T - 2
# t_{0.975, 498} ~ 1.965 (large sample: use 1.96)
t_critical = 1.96

# 95% confidence interval for alpha (daily)
ci_lo = alpha_hat - t_critical * se_alpha
ci_hi = alpha_hat + t_critical * se_alpha

# Annualise alpha for interpretability
alpha_annual = alpha_hat * TRADING_DAYS
se_alpha_annual = se_alpha * TRADING_DAYS

print("Jensen's Alpha Test (H0: alpha = 0):")
print(f"  alpha_hat (daily)    = {alpha_hat:.6f}  (annualised: {alpha_annual:.4%})")
print(f"  SE(alpha) (daily)    = {se_alpha:.6f}  (annualised: {se_alpha_annual:.4%})")
print(f"  t-statistic          = {t_stat:.4f}  (critical: ±{t_critical:.3f}  at 5% two-tailed)")
print(f"  95% CI (daily):        [{ci_lo:.6f}, {ci_hi:.6f}]")
print(f"  95% CI (annual):       [{ci_lo*TRADING_DAYS:.4%}, {ci_hi*TRADING_DAYS:.4%}]")
print()
reject = abs(t_stat) > t_critical
print(f"  Reject H0 (alpha=0) at 5%?  {'YES — unexpected; check simulation' if reject else 'NO — consistent with true alpha = 0'}")
print()
print("Note: since true alpha = 0, we expect ~95% of simulations to NOT reject H0.")
print(f"  True alpha = {alpha_true:.6f}  |  Estimated = {alpha_hat:.6f}")

§4 — Portfolio Beta: Linearity Verification

By Theorem 8.1, the beta of a portfolio equals the weighted average of individual betas: β_P = ∑ w_iβ_i.

We simulate three assets with known betas and combine them into a portfolio. We then verify that (a) the weighted average beta and (b) the beta estimated directly from the portfolio return series agree, up to estimation noise.

Why this matters: Beta linearity is what makes portfolio hedging tractable. To reduce portfolio beta from 1.2 to 0.8, you need to short enough of the market index to reduce the weighted average by 0.4 — no repricing of individual names required.

# ── Three assets with distinct true betas ─────────────────────────────────────
beta_true_vec  = np.array([0.70, 1.20, 1.50])  # AAPL-like, SPY-like, high-beta tech
sigma_eps_vec  = np.array([0.18, 0.12, 0.25]) / np.sqrt(TRADING_DAYS)  # daily idio vols
portfolio_w    = np.array([0.50, 0.30, 0.20])  # long-only, sums to 1.0

print(f"Portfolio weights: {portfolio_w}  (sum = {portfolio_w.sum():.2f})")
print(f"True betas:        {beta_true_vec}")

# ── Portfolio beta by weighted average (exact under linearity) ─────────────
beta_portfolio_formula = portfolio_w @ beta_true_vec
print(f"\nPortfolio beta (weighted average): {beta_portfolio_formula:.4f}")
print(f"  = {portfolio_w[0]}*{beta_true_vec[0]} + {portfolio_w[1]}*{beta_true_vec[1]} + {portfolio_w[2]}*{beta_true_vec[2]}")

# ── Simulate return series for all three assets ────────────────────────────────
# Re-use the same market return series r_M_excess from §1
eps_multi = rng.normal(0.0, 1.0, size=(T, 3)) * sigma_eps_vec  # shape (T, 3)
r_assets  = beta_true_vec * r_M_excess[:, np.newaxis] + eps_multi  # shape (T, 3)

# ── Portfolio return series ────────────────────────────────────────────────────
r_portfolio = r_assets @ portfolio_w  # shape (T,)

# ── Estimate portfolio beta directly from portfolio return series ─────────────
cov_pM  = np.cov(r_portfolio, r_M_excess, ddof=1)
beta_portfolio_ols = cov_pM[0, 1] / cov_pM[1, 1]

# ── Estimate each asset's beta individually, then compute weighted average ────
beta_estimated_vec = np.array([
    np.cov(r_assets[:, i], r_M_excess, ddof=1)[0, 1] / np.var(r_M_excess, ddof=1)
    for i in range(3)
])
beta_portfolio_wavg = portfolio_w @ beta_estimated_vec

print(f"\nIndividual beta estimates (OLS):")
for i in range(3):
    print(f"  Asset {i+1}: beta_true={beta_true_vec[i]:.2f}  beta_hat={beta_estimated_vec[i]:.4f}")

print(f"\nPortfolio beta comparison:")
print(f"  True (w @ beta_true):           {beta_portfolio_formula:.4f}")
print(f"  OLS on portfolio series:        {beta_portfolio_ols:.4f}")
print(f"  Weighted avg of beta estimates: {beta_portfolio_wavg:.4f}")
print()

diff1 = abs(beta_portfolio_ols   - beta_portfolio_formula)
diff2 = abs(beta_portfolio_wavg  - beta_portfolio_formula)
print(f"  |OLS - formula| = {diff1:.4f}  (estimation noise, should be < 3*SE)")
print(f"  |WAvg - formula| = {diff2:.4f}  (same estimation noise)")
print()

# SE of portfolio beta estimate
r_portfolio_eps = r_portfolio - (np.mean(r_portfolio) - beta_portfolio_ols * np.mean(r_M_excess)) - beta_portfolio_ols * r_M_excess
sigma_p_eps = np.std(r_portfolio_eps, ddof=2)
se_beta_p   = sigma_p_eps / (np.std(r_M_excess, ddof=1) * np.sqrt(T))
print(f"  SE(portfolio beta) = {se_beta_p:.4f}")
print(f"  Deviations within 3*SE? OLS: {'✓' if diff1 < 3*se_beta_p else '✗'}  WAvg: {'✓' if diff2 < 3*se_beta_p else '✗'}")

§5 — Validation Summary

Final cell: verify all key claims from the article against the computed values. Each check is a precise numerical assertion.

print("=" * 65)
print("VALIDATION SUMMARY — CAPM Single-Factor Model")
print("=" * 65)
print()

checks = []

# 1. Beta estimate close to true value (within 2 SE)
c1 = abs(beta_hat - BETA_TRUE) < 2 * se_beta
checks.append(("beta_hat within 2 SE of beta_true (1.35)", c1))
print(f"1. beta_hat = {beta_hat:.4f}  true = {BETA_TRUE:.2f}  SE = {se_beta:.4f}  {'✓' if c1 else '✗'}")

# 2. Variance decomposition sums to total (up to floating-point noise)
c2 = abs((var_systematic + var_idio) - var_total) < 1e-12
checks.append(("systematic + idiosyncratic = total variance (identity)", c2))
print(f"2. var_sys + var_idio = {var_systematic + var_idio:.10f}  total = {var_total:.10f}  {'✓' if c2 else '✗'}")

# 3. Two R^2 computations agree
c3 = abs(R2_decomp - R2_corr) < 1e-10
checks.append(("R^2 from decomposition == R^2 from correlation^2", c3))
print(f"3. R2_decomp = {R2_decomp:.6f}  R2_corr = {R2_corr:.6f}  diff = {abs(R2_decomp-R2_corr):.2e}  {'✓' if c3 else '✗'}")

# 4. Alpha t-stat small (|t| < 1.96 since true alpha = 0)
c4 = abs(t_stat) < 1.96
checks.append(("Jensen's alpha t-stat < 1.96 (true alpha = 0)", c4))
print(f"4. |t_alpha| = {abs(t_stat):.4f}  < 1.96?  {'✓' if c4 else '✗'}")

# 5. Portfolio beta linearity: OLS on portfolio vs weighted average
c5 = abs(beta_portfolio_ols - beta_portfolio_wavg) < 3 * se_beta_p
checks.append(("Portfolio beta OLS ≈ weighted average of individual betas", c5))
print(f"5. beta_P_ols = {beta_portfolio_ols:.4f}  beta_P_wavg = {beta_portfolio_wavg:.4f}  diff = {abs(beta_portfolio_ols - beta_portfolio_wavg):.4f}  {'✓' if c5 else '✗'}")

# 6. CAPM SML: verify E[R_i] = rf + beta * MRP with true params
expected_return_capm = RF_ANNUAL + BETA_TRUE * MU_M_ANNUAL
c6 = abs(expected_return_capm - (RF_ANNUAL + BETA_TRUE * MU_M_ANNUAL)) < 1e-12
checks.append(("SML formula: E[R_i] = rf + beta * MRP (algebraic check)", c6))
print(f"6. E[R_i] = {RF_ANNUAL:.2%} + {BETA_TRUE} * {MU_M_ANNUAL:.2%} = {expected_return_capm:.2%}  {'✓' if c6 else '✗'}")

# 7. Beta linearity: w @ beta_true_vec == formula value exactly
c7 = abs(portfolio_w @ beta_true_vec - beta_portfolio_formula) < 1e-14
checks.append(("beta linearity: w @ beta_true_vec exact (floating-point)", c7))
print(f"7. w @ beta_true = {portfolio_w @ beta_true_vec:.6f}  formula = {beta_portfolio_formula:.6f}  {'✓' if c7 else '✗'}")

print()
passed = sum(1 for _, ok in checks if ok)
print(f"Results: {passed}/{len(checks)} claims verified {'✓' if passed == len(checks) else '— check failures above'}")
if passed == len(checks):
    print("All validations passed — article claims confirmed numerically.")