Worked Solutions

# ============================================================================
# NB01 Exercise 1 — Hand Computation (verified with NumPy)
# ============================================================================

x = np.array([2, 5, 3, 8, 7], dtype=float)
y = np.array([4, 7, 5, 10, 9], dtype=float)
n = len(x)

# ── Part 1: Mean, Variance, Standard Deviation ──
mean_x = x.sum() / n
mean_y = y.sum() / n

var_x = ((x - mean_x)**2).sum() / (n - 1)   # sample variance (Bessel's correction)
var_y = ((y - mean_y)**2).sum() / (n - 1)

std_x = np.sqrt(var_x)
std_y = np.sqrt(var_y)

print("PART 1: Descriptive Statistics")
print("=" * 50)
print(f"  x = {x}")
print(f"  y = {y}")
print(f"\n  mean(x) = {x.sum()}/{n} = {mean_x:.2f}")
print(f"  mean(y) = {y.sum()}/{n} = {mean_y:.2f}")
print(f"\n  var(x)  = Σ(xᵢ − x̄)² / (n−1) = {((x - mean_x)**2).sum():.2f} / {n-1} = {var_x:.4f}")
print(f"  var(y)  = {((y - mean_y)**2).sum():.2f} / {n-1} = {var_y:.4f}")
print(f"\n  std(x)  = √{var_x:.4f} = {std_x:.4f}")
print(f"  std(y)  = √{var_y:.4f} = {std_y:.4f}")

# Verify
print(f"\n  Verify with NumPy:  std(x, ddof=1) = {x.std(ddof=1):.4f}  ✓")

# ── Part 2: Covariance and Correlation ──
cov_xy = ((x - mean_x) * (y - mean_y)).sum() / (n - 1)
corr_xy = cov_xy / (std_x * std_y)

print(f"\nPART 2: Covariance and Correlation")
print("=" * 50)
print(f"  Cov(x,y) = Σ(xᵢ-x̄)(yᵢ-ȳ)/(n-1) = {((x-mean_x)*(y-mean_y)).sum():.2f}/{n-1} = {cov_xy:.4f}")
print(f"  Corr(x,y) = Cov(x,y)/(s_x · s_y) = {cov_xy:.4f}/({std_x:.4f}·{std_y:.4f}) = {corr_xy:.4f}")
print(f"\n  Verify: np.corrcoef = {np.corrcoef(x, y)[0, 1]:.4f}  ✓")

# ── Part 3: Normal Equations ──
X = np.column_stack([np.ones(n), x])   # design matrix with intercept
XtX     = X.T @ X
XtX_inv = np.linalg.inv(XtX)
Xty     = X.T @ y
beta_hat = XtX_inv @ Xty

print(f"\nPART 3: OLS via Normal Equations")
print("=" * 50)
print(f"  Design matrix X (n×2):\n{X}")
print(f"\n  X'X =\n{XtX}")
print(f"\n  (X'X)⁻¹ =\n{XtX_inv}")
print(f"\n  X'y = {Xty}")
print(f"\n  β̂ = (X'X)⁻¹ X'y = [{beta_hat[0]:.6f}, {beta_hat[1]:.6f}]")
print(f"\n  Interpretation:  ŷ = {beta_hat[0]:.4f} + {beta_hat[1]:.4f} · x")
print(f"    → intercept ≈ {beta_hat[0]:.2f}, slope ≈ {beta_hat[1]:.2f}")

# Quick sanity check: slope should be close to Cov(x,y)/Var(x)
print(f"\n  Sanity: Cov(x,y)/Var(x) = {cov_xy/var_x:.6f}  (should match slope)  ✓")

# ============================================================================
# NB01 Exercise 2 — CLT with Beta(0.5, 0.5)
# ============================================================================

np.random.seed(42)
n_sims = 5000
a, b = 0.5, 0.5
mu_beta  = a / (a + b)                          # = 0.5
var_beta = (a * b) / ((a + b)**2 * (a + b + 1)) # = 0.125

sample_sizes = [5, 30, 100]

fig, axes = plt.subplots(1, len(sample_sizes) + 1, figsize=(18, 4))

# Panel 0: Show the raw Beta(0.5, 0.5) density (bimodal)
x_grid = np.linspace(0.001, 0.999, 500)
axes[0].plot(x_grid, stats.beta.pdf(x_grid, a, b), 'k-', linewidth=2)
axes[0].set_title('Beta(0.5, 0.5) density\n(bimodal!)', fontweight='bold')
axes[0].set_xlabel('x')
axes[0].set_ylabel('Density')
axes[0].grid(True, alpha=0.3)

# Panels 1–3: Histograms of standardised sample means
for i, n in enumerate(sample_sizes):
    means = np.array([np.random.beta(a, b, size=n).mean() for _ in range(n_sims)])
    z = (means - mu_beta) / np.sqrt(var_beta / n)   # standardise to N(0,1) if CLT holds
    
    ax = axes[i + 1]
    ax.hist(z, bins=50, density=True, alpha=0.6, color='steelblue', edgecolor='white')
    
    # Overlay N(0,1) density
    t_grid = np.linspace(-4, 4, 200)
    ax.plot(t_grid, stats.norm.pdf(t_grid), 'r-', linewidth=2, label='N(0,1)')
    
    ax.set_title(f'n = {n}', fontweight='bold', fontsize=13)
    ax.set_xlabel('Standardised z')
    ax.legend(fontsize=9)
    ax.set_xlim(-4, 4)
    ax.grid(True, alpha=0.3)

fig.suptitle('CLT Demonstration: Beta(0.5, 0.5) — Sample Means Become Normal',
             fontweight='bold', fontsize=13, y=1.03)
plt.tight_layout()
plt.show()

# Symmetry is distinct from the absolute third moment in Berry–Esseen.
print("The bounded Beta distribution has finite moments and satisfies the iid CLT.")
print("Its signed third central moment is zero, but E[|X-mu|^3] is positive.")
print("Symmetry removes a skewness term in an Edgeworth expansion; it does not make the Berry–Esseen bound zero.")
print("Inspect the histograms: visual similarity at a given n is not a universal error bound.")

x_orig = np.array([2., -1., 3., 0., 1.])
y_orig = np.array([3., -2., 4., 1., 0.])
y_changed = y_orig * np.array([1, 1, 1, 3, 3])
res = sm.OLS(y_changed, sm.add_constant(x_orig)).fit()
# No finite-sample correction here, to match the hand-derived HC0/HAC formula.
hc0 = res.get_robustcov_results(cov_type='HC0', use_t=False)
hac = res.get_robustcov_results(cov_type='HAC', maxlags=1, use_correction=False, use_t=False)
comparison = pd.DataFrame({'Classic': res.bse, 'HC0': hc0.bse, 'HAC(1), uncorrected': hac.bse}, index=['Intercept', 'Slope']).T
comparison['Slope SE change %'] = 100*(comparison.Slope/res.bse[1]-1)
display(comparison.round(5))
print('Largest absolute slope-SE change:', comparison['Slope SE change %'].abs().idxmax())
print('Five observations illustrate mechanics; robust asymptotic inference is unreliable here.')

# ============================================================================
# NB01 Exercise 4 — Monte Carlo: Rejection Rate vs Autocorrelation
# ============================================================================

np.random.seed(0)
num_mc = 2000
n_mc = 150
rho_values = [0.0, 0.3, 0.6, 0.9]

results_mc = []  # list of (rho, classic_rate, nw_rate)

for rho in rho_values:
    reject_classic = 0
    reject_nw = 0
    
    for sim in range(num_mc):
        x_mc = np.random.randn(n_mc) * 0.04
        
        # AR(1) errors with heteroscedasticity and fat tails
        innov = stats.t.rvs(df=5, size=n_mc + 500)
        eps = np.zeros(n_mc + 500)
        for t in range(n_mc + 500):
            eps[t] = rho * eps[t-1] + innov[t] if t > 0 else innov[t]
        eps = eps[-n_mc:] * (0.01 + 0.4 * np.abs(x_mc))
        
        y_mc = 0.0 + 0.0 * x_mc + eps   # true β = 0
        
        X_mc = sm.add_constant(x_mc)
        res_mc = sm.OLS(y_mc, X_mc).fit()
        lag_mc = int(np.floor(4 * (n_mc / 100) ** (2/9)))
        res_nw_mc = res_mc.get_robustcov_results(cov_type='HAC', maxlags=lag_mc, use_correction=True, use_t=False)
        
        if res_mc.pvalues[1] < 0.05:
            reject_classic += 1
        if res_nw_mc.pvalues[1] < 0.05:
            reject_nw += 1
    
    rate_c = reject_classic / num_mc * 100
    rate_nw = reject_nw / num_mc * 100
    results_mc.append((rho, rate_c, rate_nw))
    print(f"  ρ = {rho:.1f}:  Classic = {rate_c:.1f}%,  Newey-West = {rate_nw:.1f}%")

# ── Plot ──
rhos    = [r[0] for r in results_mc]
classic = [r[1] for r in results_mc]
nw      = [r[2] for r in results_mc]

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(rhos, classic, 'o-', color='#d9534f', linewidth=2, markersize=8, label='Classic OLS')
ax.plot(rhos, nw,      's-', color='#5cb85c', linewidth=2, markersize=8, label='Newey-West HAC')
ax.axhline(5, color='black', linestyle='--', linewidth=1.5, label='Nominal 5% level')
ax.set_xlabel('Autocorrelation coefficient ρ', fontsize=12)
ax.set_ylabel('Rejection rate (%)', fontsize=12)
ax.set_title('False Rejection Rate Under H₀ vs Autocorrelation Strength\n'
             f'({num_mc} simulations, n={n_mc}, heteroscedastic t(5) errors)',
             fontweight='bold', fontsize=12)
ax.set_xticks(rho_values)
ax.set_ylim(0, max(classic + nw + [5]) * 1.25)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print('Monte Carlo standard errors, percentage points:')
for rho, classic_rate, hac_rate in results_mc:
    print(rho, [round(100*np.sqrt((r/100)*(1-r/100)/num_mc), 2) for r in [classic_rate, hac_rate]])
print('HAC is asymptotic and need not attain the nominal 5% size in this design.')
print('The slope score is x_t*epsilon_t: iid mean-zero x can weaken score dependence even when errors are persistent.')

# ============================================================================
# NB02 Exercise 1 — Stationarity: AAPL Prices vs Log Returns
# ============================================================================

prices = load_prices('AAPL').loc['2015-01-01':]
coverage('AAPL daily prices', prices)

log_returns = np.log(prices / prices.shift(1)).dropna()

# ── Part 1: ADF tests ──
print("ADF TEST RESULTS")
print("=" * 60)

adf_price = adfuller(prices.values, maxlag=20, autolag='AIC')
print(f"\nAAPL PRICES:")
print(f"  ADF statistic:  {adf_price[0]:.4f}")
print(f"  p-value:        {adf_price[1]:.4f}")
print(f"  Conclusion:     {'Reject unit root' if adf_price[1] < 0.05 else 'Do not reject unit root'}")

adf_ret = adfuller(log_returns.values, maxlag=20, autolag='AIC')
print(f"\nAAPL LOG RETURNS:")
print(f"  ADF statistic:  {adf_ret[0]:.4f}")
print(f"  p-value:        {adf_ret[1]:.6f}")
print(f"  Conclusion:     {'Reject unit root' if adf_ret[1] < 0.05 else 'Do not reject unit root'}")

print('ADF uses a constant, AIC lag selection and a unit-root null; it does not test full stationarity.')

# ── Part 2: Rolling 60-day mean and std dev ──
fig, axes = plt.subplots(2, 2, figsize=(15, 8))

axes[0, 0].plot(prices.index, prices, color='darkblue', linewidth=0.8)
axes[0, 0].set_title('AAPL Daily Close Price', fontweight='bold')
axes[0, 0].set_ylabel('Price ($)')
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].plot(log_returns.index, log_returns * 100, color='steelblue', linewidth=0.5, alpha=0.7)
axes[0, 1].set_title('AAPL Daily Log Returns', fontweight='bold')
axes[0, 1].set_ylabel('Return (%)')
axes[0, 1].grid(True, alpha=0.3)

rolling_mean = log_returns.rolling(60).mean() * 100
axes[1, 0].plot(rolling_mean.index, rolling_mean, color='darkred', linewidth=1)
axes[1, 0].axhline(log_returns.mean() * 100, color='black', linestyle='--', alpha=0.5,
                     label=f'Full-sample mean = {log_returns.mean()*100:.3f}%')
axes[1, 0].set_title('Rolling 60-Day Mean of Log Returns', fontweight='bold')
axes[1, 0].set_ylabel('Mean return (%)')
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)

rolling_std = log_returns.rolling(60).std() * 100
axes[1, 1].plot(rolling_std.index, rolling_std, color='darkorange', linewidth=1)
axes[1, 1].axhline(log_returns.std() * 100, color='black', linestyle='--', alpha=0.5,
                     label=f'Full-sample std = {log_returns.std()*100:.3f}%')
axes[1, 1].set_title('Rolling 60-Day Std Dev (Volatility)', fontweight='bold')
axes[1, 1].set_ylabel('Std dev (%)')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)

plt.suptitle('AAPL: Stationarity Analysis', fontweight='bold', fontsize=14, y=1.01)
plt.tight_layout()
plt.show()

print('Rolling moments are diagnostics. Changing conditional volatility is compatible with stationary GARCH, but neither the plot nor ADF proves that model.')

fig, axes = plt.subplots(1, 2, figsize=(13, 4))
for ax, (label, series) in zip(axes, [('2000–2011', ff3_full.Mkt_RF.loc[:'2011-12']), ('2012–latest', ff3_full.Mkt_RF.loc['2012-01':])]):
    lm, p, f, fp = het_arch((series-series.mean()).to_numpy(), nlags=6, ddof=1)
    print(f'{label}: N={len(series)}, monthly SD={series.std():.2%}, ARCH LM={lm:.3f}, p={p:.4g}')
    ax.plot(series.index, 100*series)
    ax.set(title=label, ylabel='Market excess return (%)')
plt.tight_layout()
plt.show()
print('Differences reflect both the sample and test power. Non-rejection does not establish constant conditional variance.')

tickers = ['SPY', 'IWN', 'VTV']
monthly_etfs = pd.concat({tic: monthly_returns(load_prices(tic)) for tic in tickers}, axis=1)
comparison_data = ff3_short.join(monthly_etfs, how='inner').dropna()
coverage('Common ETF comparison', comparison_data)
results_dict = {tic: fit_hac(comparison_data[tic]-comparison_data.RF, comparison_data[['Mkt_RF','SMB','HML']]) for tic in tickers}
rows = []
for tic, res in results_dict.items():
    rows.append({'ETF': tic, 'Alpha %/yr': 1200*res.params[0], 'Alpha p': res.pvalues[0], 'Market beta': res.params[1], 'SMB': res.params[2], 'HML': res.params[3], 'R²': res.rsquared})
comparison_table = pd.DataFrame(rows).set_index('ETF')
display(comparison_table.round(4))
print('Highest sample R²:', comparison_table['R²'].idxmax())
print('IWN minus SPY SMB:', results_dict['IWN'].params[2]-results_dict['SPY'].params[2])
print('High R² measures within-sample fit, not superior investment performance.')

d = comparison_data
base = sm.OLS((d.SPY-d.RF).to_numpy(), sm.add_constant(d[['Mkt_RF','SMB','HML']].to_numpy())).fit()
rows = []
for lag in sorted({0, nw_lags(len(d)), 8, 12}):
    r = fit_hac(d.SPY-d.RF, d[['Mkt_RF','SMB','HML']], lags=lag)
    rows.append({'Lag': lag, 'Monthly alpha': r.params[0], 'Alpha SE': r.bse[0], 'Alpha p': r.pvalues[0], 'SMB SE': r.bse[2]})
display(pd.DataFrame(rows).round(6))
np.testing.assert_allclose(fit_hac(d.SPY-d.RF, d[['Mkt_RF','SMB','HML']], lags=0).cov_params(), base.get_robustcov_results(cov_type='HC1').cov_params())
print('Lag zero with n/(n-k) correction equals HC1; uncorrected lag zero equals HC0.')
print('Coefficients stay fixed. HAC standard errors need not increase monotonically with lag length.')

aapl_monthly = monthly_returns(load_prices('AAPL')).rename('AAPL')
d = ff3_short.join(aapl_monthly, how='inner').dropna()
coverage('AAPL FF3', d)
r = fit_hac(d.AAPL-d.RF, d[['Mkt_RF','SMB','HML']])
display(pd.DataFrame({'Estimate': r.params, 'SE': r.bse, 'p': r.pvalues, '95% low': r.conf_int()[:,0], '95% high': r.conf_int()[:,1]}, index=['Alpha monthly','Mkt_RF','SMB','HML']).round(5))
print(f'Annualised arithmetic alpha: {1200*r.params[0]:.2f}%; R²: {r.rsquared:.4f}')
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
sm.qqplot(r.resid, line='s', ax=axes[0])
axes[1].plot(d.index, 100*r.resid)
axes[1].set(title='AAPL FF3 residuals', ylabel='Monthly residual (%)')
plt.tight_layout()
plt.show()
display(acorr_ljungbox(r.resid, lags=[6,12], return_df=True))
print('Single-stock alpha mixes model omissions and stock-specific outcomes; it does not identify skill.')

small = pd.concat({tic: monthly_returns(load_prices(tic)) for tic in ['IWM','IWN']}, axis=1)
d = factors.join(small, how='inner').dropna()
coverage('Common small-cap comparison', d)
rows = []
for tic in ['IWM','IWN']:
    for label, cols in [('FF3', ['Mkt_RF','SMB_FF3','HML']), ('FF5', ['Mkt_RF','SMB','HML','RMW','CMA'])]:
        r = fit_hac(d[tic]-d.RF, d[cols])
        rows.append({'ETF': tic, 'Model': label, 'Alpha %/yr': 1200*r.params[0], 'Alpha p': r.pvalues[0], 'SMB': r.params[2], 'HML': r.params[3], 'HML p': r.pvalues[3], 'R²': r.rsquared})
display(pd.DataFrame(rows).round(4))
print('FF3 uses its own SMB. Published FF3 and FF5 are not exactly nested, so R² need not increase.')
print('Compare measured SMB and HML estimates; a blend mandate does not imply zero HML exposure.')

# ============================================================================
# NB04 Exercise 2 — Rolling Sharpe Ratios
# ============================================================================

# Build strategy returns (equal-weight combinations)
strats = pd.DataFrame({
    'HML only':        factors['HML'],
    'Value + Prof':    (factors['HML'] + factors['RMW']) / 2,
    'Val + Prof + Mom':(factors['HML'] + factors['RMW'] + factors['UMD']) / 3,
})

# Rolling 36-month annualised Sharpe ratio
window = 36
rolling_sharpe = (strats.rolling(window).mean() / strats.rolling(window).std()) * np.sqrt(12)
rolling_sharpe = rolling_sharpe.dropna()

# Plot
fig, ax = plt.subplots(figsize=(14, 5.5))
colors = ['#1f77b4', '#ff7f0e', '#2ca02c']
for col, c in zip(rolling_sharpe.columns, colors):
    ax.plot(rolling_sharpe.index, rolling_sharpe[col], label=col, linewidth=1.2, color=c)

ax.axhline(0, color='black', linewidth=0.5)
ax.fill_between(rolling_sharpe.index, 0, rolling_sharpe.min(axis=1),
                where=rolling_sharpe.min(axis=1) < 0, alpha=0.1, color='red')
ax.set_title('36-Month Rolling Sharpe Ratios: Single vs Combined Factor Strategies',
             fontweight='bold', fontsize=13)
ax.set_ylabel('Annualised Sharpe Ratio')
ax.legend(loc='upper right', framealpha=0.9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Analysis
print("SUMMARY STATISTICS")
print("=" * 60)
print(f"{'Strategy':>25s}  {'Mean SR':>8s}  {'Std SR':>8s}  {'% < 0':>8s}")
print("-" * 60)
for col in rolling_sharpe.columns:
    sr = rolling_sharpe[col]
    print(f"{col:>25s}  {sr.mean():8.3f}  {sr.std():8.3f}  {(sr < 0).mean()*100:7.1f}%")

# Check if all three are simultaneously negative
all_neg = (rolling_sharpe < 0).all(axis=1)
if all_neg.any():
    periods = rolling_sharpe.index[all_neg]
    print(f"\nFirst/last months with all three negative (not necessarily consecutive):")
    print(f"  {periods[0].strftime('%Y-%m')} to {periods[-1].strftime('%Y-%m')} "
          f"({all_neg.sum()} months)")
else:
    print("\nNo period where ALL three strategies are simultaneously negative.")

print('The displayed sample determines whether diversification improved Sharpe or stability.')
print('These overlapping rolling windows are dependent and are not independent performance tests.')

rolling_3m = (1+factors.UMD).rolling(3).apply(np.prod, raw=True)-1
worst_end = rolling_3m.idxmin()
worst_start = worst_end-pd.DateOffset(months=2)
hedged = (factors.UMD+factors.HML)/2
crash = factors.loc[worst_start:worst_end, ['UMD','HML']].copy()
crash['Market total'] = (factors.Mkt_RF+factors.RF).reindex(crash.index)
crash['UMD + HML'] = hedged.reindex(crash.index)
print('Worst UMD window:', worst_start.date(), 'to', worst_end.date())
display((100*crash).round(2))
losses = (1+crash).prod()-1
print('Compounded three-month returns (%):')
display((100*losses).round(2))
print(f'Combination minus UMD return: {100*(losses["UMD + HML"]-losses.UMD):.2f} percentage points')
fig, axes = plt.subplots(1, 2, figsize=(13, 4))
near = pd.DataFrame({'UMD': factors.UMD, 'HML': factors.HML, 'UMD + HML': hedged, 'Market total': factors.Mkt_RF+factors.RF}).loc[worst_start-pd.DateOffset(months=6):worst_end+pd.DateOffset(months=6)]
(1+near).cumprod().plot(ax=axes[0], title='Hypothetical indices around the crash')
for name, series in [('UMD', factors.UMD), ('UMD + HML', hedged)]:
    wealth = (1+series).cumprod()
    dd = wealth/wealth.cummax().clip(lower=1)-1
    axes[1].plot(dd.index, dd*100, label=name)
    print(f'{name} full-sample max drawdown: {dd.min():.2%}')
axes[1].set(title='Drawdown over time', ylabel='Drawdown (%)')
axes[1].legend()
plt.tight_layout()
plt.show()
print('Compare the displayed market rebound with the momentum-crash mechanism; diversification is not a guaranteed hedge.')

beta_ports = load_beta_portfolios().loc['2000-01-01':]
common_idx = beta_ports.index.intersection(factors.index)
beta_ports = beta_ports.loc[common_idx]
rf_aligned = factors.RF.reindex(common_idx)
mkt_aligned = factors.Mkt_RF.reindex(common_idx)

# ── Step 2: Split into two subperiods ──
split_date = '2011-12-31'
periods = {
    '2000–2011': (beta_ports.loc[:split_date], rf_aligned.loc[:split_date], mkt_aligned.loc[:split_date]),
    '2012–latest': (beta_ports.loc['2012-01-01':], rf_aligned.loc['2012-01-01':], mkt_aligned.loc['2012-01-01':]),
}

fig, axes = plt.subplots(1, 2, figsize=(16, 6))

for ax, (period_name, (bp, rf, mkt)) in zip(axes, periods.items()):
    avg_excess = []
    realized_betas = []
    for col in bp.columns:
        excess = bp[col] - rf
        X = sm.add_constant(mkt.values)
        res = sm.OLS(excess.values, X).fit()
        realized_betas.append(res.params[1])
        avg_excess.append(excess.mean() * 12 * 100)
    
    avg_excess = np.array(avg_excess)
    realized_betas = np.array(realized_betas)
    
    ax.scatter(realized_betas, avg_excess, s=120, c='darkblue', zorder=5,
               edgecolors='white', linewidth=1.5)
    for i, label in enumerate(bp.columns):
        ax.annotate(label, (realized_betas[i], avg_excess[i]),
                    textcoords="offset points", xytext=(8, 5), fontsize=9)
    
    # CAPM line
    mkt_prem = mkt.mean() * 12 * 100
    beta_range = np.linspace(0, max(realized_betas) * 1.1, 100)
    ax.plot(beta_range, beta_range * mkt_prem, 'r--', linewidth=2, label='CAPM', alpha=0.7)
    
    # Empirical fit
    slope, intercept = np.polyfit(realized_betas, avg_excess, 1)
    ax.plot(beta_range, slope * beta_range + intercept, 'b-', linewidth=2.5,
            label=f'Empirical (slope={slope:.1f})', alpha=0.9)
    
    ax.set_xlabel('Realized Beta', fontsize=11)
    ax.set_ylabel('Ann. Excess Return (%)', fontsize=11)
    ax.set_title(f'Empirical SML: {period_name}\n(CAPM slope={mkt_prem:.1f})',
                 fontweight='bold', fontsize=12)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

plt.suptitle('Beta Anomaly Subsample Robustness', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

# ── Step 3: BAB alpha in each subperiod ──
print("\n" + "=" * 60)
print("UNSCALED D1–D10 SPREAD ALPHA BY SUBPERIOD (FF5 + Momentum controls)")
print("=" * 60)

factor_cols_full = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD']

for period_name, (bp, rf, mkt) in periods.items():
    # Simple BAB = D1 excess - D10 excess
    bab = (bp['D1'] - rf) - (bp['D10'] - rf)
    ff = factors.loc[bab.index, factor_cols_full].dropna()
    bab = bab.loc[ff.index]
    
    X_ff = sm.add_constant(ff.values)
    res_bab = fit_hac(bab, ff, lags=6)
    
    alpha_ann = res_bab.params[0] * 12 * 100
    print(f"\n  {period_name} ({len(bab)} months)")
    print(f"    Alpha (ann.): {alpha_ann:+.2f}%   t = {res_bab.tvalues[0]:.2f}   "
          f"p = {res_bab.pvalues[0]:.4f}")
    sig = "SIGNIFICANT" if res_bab.pvalues[0] < 0.05 else "NOT significant"
    print(f"    Conclusion:   {sig}")

print('Compare estimated slopes with the CAPM reference; do not assume their signs.')
print('A significant spread intercept is conditional on this model and sample; this is not published BAB.')
print('Changes in access to leverage could affect demand for high-beta stocks, but this two-period comparison cannot identify that causal mechanism.')
print('A fitted cross-sectional line through ten noisy estimated betas is descriptive, without an errors-in-variables correction or pricing test.')