4. Beyond Three Factors: Profitability, Investment and Evidence

# ============================================================================
# Setup: Import Libraries
# ============================================================================

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


try:
    plt.style.use('seaborn-v0_8-darkgrid')
except OSError:
    plt.style.use('default')
sns.set_palette("husl")
plt.rcParams["font.family"] = "DejaVu Sans"
np.random.seed(42)

print("Libraries imported successfully!")
from analysis_utils import AS_OF, END_EXCLUSIVE, load_french, load_prices, monthly_returns, coverage, fit_hac, mean_hac, nw_lags, load_beta_portfolios

download_ff_csv = load_french

# --- 3-Factor data ---
print("Downloading Fama-French 3-Factor data...")
ff3 = download_ff_csv("F-F_Research_Data_Factors_CSV.zip")
ff3.columns = ['Mkt_RF', 'SMB', 'HML', 'RF']

# --- 5-Factor data ---
print("Downloading Fama-French 5-Factor data...")
ff5 = download_ff_csv("F-F_Research_Data_5_Factors_2x3_CSV.zip")
ff5.columns = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'RF']

# --- Momentum factor ---
print("Downloading Momentum factor (UMD)...")
mom = download_ff_csv("F-F_Momentum_Factor_CSV.zip")
mom.columns = ['UMD']

# --- Short-Term Reversal factor ---
print("Downloading Short-Term Reversal factor (ST_Rev)...")
st_rev = download_ff_csv("F-F_ST_Reversal_Factor_CSV.zip")
st_rev.columns = ['ST_Rev']

# --- Long-Term Reversal factor ---
print("Downloading Long-Term Reversal factor (LT_Rev)...")
lt_rev = download_ff_csv("F-F_LT_Reversal_Factor_CSV.zip")
lt_rev.columns = ['LT_Rev']

# Merge everything on the common date index
factors = ff5.copy()
factors = factors.join(mom, how='inner')
factors = factors.join(st_rev, how='inner')
factors = factors.join(lt_rev, how='inner')

# Filter to analysis period
factors = factors.loc['2000-01-01':].join(ff3['SMB'].rename('SMB_FF3'), how='inner')
coverage('US combined factor sample', factors)

print(f"Combined factor dataset: {factors.shape[0]} months, {factors.shape[1]} columns")
print(f"Date range: {factors.index[0].date()} to {factors.index[-1].date()}")
print(f"Columns: {list(factors.columns)}")
print(f"\nFirst 5 rows:")
print(factors.head())

def run_ff_regression(y, X_cols, data, model_name="Model", stock_name="Stock"):
    columns = list(X_cols)
    if columns in [['Mkt_RF', 'SMB', 'HML'], ['Mkt_RF', 'SMB', 'HML', 'UMD']]:
        columns = ['SMB_FF3' if c == 'SMB' else c for c in columns]
    target = data['Stock_Excess_Return'] if y is None else y
    return fit_hac(target, data[columns])

etf_tickers = {'SPY': 'SPY (S&P 500)', 'IWN': 'IWN (Small Cap Value)',
               'VTV': 'VTV (Large Cap Value)'}
etf_dataframes = {}
for ticker in etf_tickers:
    returns = monthly_returns(load_prices(ticker)).rename('ETF_Return')
    edf = factors.join(returns, how='inner').dropna()
    edf['Stock_Excess_Return'] = edf['ETF_Return'] - edf['RF']
    etf_dataframes[ticker] = edf
common_idx = etf_dataframes['SPY'].index
for edf in etf_dataframes.values():
    common_idx = common_idx.intersection(edf.index)
etf_dataframes = {t: d.loc[common_idx].copy() for t, d in etf_dataframes.items()}
df = etf_dataframes['SPY'].copy()
coverage('Common ETF regression sample (limited by VTV inception)', df)

# ============================================================================
# DDM Illustration: How Profitability, Investment, and B/M Map to Expected Returns
# ============================================================================
# This cell provides a numerical illustration of the Fama-French (2015) DDM
# comparative statics. We solve for the implied discount rate (expected return)
# under different assumptions about profitability and investment growth.
#
# The simplified DDM (Gordon Growth variant with clean surplus):
#   M = (Y - dB) / (r - g)   =>   r = (Y - dB)/M + g
# where Y = earnings, dB = net investment, M = market value, g = growth rate.
# Equivalently, dividing by book equity B:
#   M/B = (Y/B - dB/B) / (r - g)
# so:
#   r = (Y/B - dB/B) / (M/B) + g

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker

# --- Parameters ---
growth_rate = 0.02          # long-run nominal growth rate (g)
profitability_range = np.linspace(0.06, 0.20, 100)  # Y/B: ROE from 6% to 20%
investment_range = np.linspace(0.01, 0.11, 100)      # dB/B: reinvestment rate

# Fix one dimension and vary the other to show comparative statics

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# --- Panel A: Expected return vs. Profitability (Y/B) ---
# Fix: M/B = 1.5,  dB/B = 0.05
mb_fixed = 1.5
inv_fixed = 0.05
r_vs_prof = (profitability_range - inv_fixed) / mb_fixed + growth_rate

ax = axes[0]
ax.plot(profitability_range * 100, r_vs_prof * 100, color='#2166ac', linewidth=2.5)
ax.set_xlabel('Profitability (Y/B, %)', fontsize=11)
ax.set_ylabel('Implied Expected Return r (%)', fontsize=11)
ax.set_title('Panel A: Higher Profitability → Higher r\n(Motivates RMW)', fontsize=11, fontweight='bold')
ax.axhline(y=(0.10 - inv_fixed) / mb_fixed * 100 + growth_rate * 100, color='gray',
           linestyle='--', alpha=0.4)
ax.annotate('Holding M/B and investment fixed,\nhigher profitability requires\na higher discount rate',
            xy=(0.12 * 100, ((0.12 - inv_fixed) / mb_fixed + growth_rate) * 100),
            xytext=(0.08 * 100, ((0.17 - inv_fixed) / mb_fixed + growth_rate) * 100),
            fontsize=9, fontstyle='italic',
            arrowprops=dict(arrowstyle='->', color='gray'))
ax.yaxis.set_major_formatter(mticker.FormatStrFormatter('%.1f'))
ax.grid(True, alpha=0.3)

# --- Panel B: Expected return vs. Investment (dB/B) ---
# Fix: M/B = 1.5,  Y/B = 0.12
prof_fixed = 0.12
r_vs_inv = (prof_fixed - investment_range) / mb_fixed + growth_rate

ax = axes[1]
ax.plot(investment_range * 100, r_vs_inv * 100, color='#b2182b', linewidth=2.5)
ax.set_xlabel('Investment (ΔB/B, %)', fontsize=11)
ax.set_ylabel('Implied Expected Return r (%)', fontsize=11)
ax.set_title('Panel B: Higher Investment → Lower r\n(Motivates CMA)', fontsize=11, fontweight='bold')
ax.annotate('Holding M/B and profitability fixed,\nhigher investment (aggressive)\nrequires a lower discount rate',
            xy=(0.08 * 100, ((prof_fixed - 0.08) / mb_fixed + growth_rate) * 100),
            xytext=(0.04 * 100, ((prof_fixed - 0.10) / mb_fixed + growth_rate) * 100),
            fontsize=9, fontstyle='italic',
            arrowprops=dict(arrowstyle='->', color='gray'))
ax.yaxis.set_major_formatter(mticker.FormatStrFormatter('%.1f'))
ax.grid(True, alpha=0.3)

# --- Panel C: Expected return vs. B/M (inverse of M/B) ---
# Fix: Y/B = 0.12,  dB/B = 0.05
bm_range = np.linspace(0.3, 1.5, 100)   # B/M from 0.3 (growth) to 1.5 (deep value)
mb_range = 1.0 / bm_range
r_vs_bm = (prof_fixed - inv_fixed) / mb_range + growth_rate

ax = axes[2]
ax.plot(bm_range, r_vs_bm * 100, color='#1a9850', linewidth=2.5)
ax.set_xlabel('Book-to-Market (B/M)', fontsize=11)
ax.set_ylabel('Implied Expected Return r (%)', fontsize=11)
ax.set_title('Panel C: Higher B/M → Higher r\n(Motivates HML)', fontsize=11, fontweight='bold')
ax.annotate('Holding profitability and investment\nfixed, cheaper stocks (high B/M)\nhave higher expected returns',
            xy=(1.0, ((prof_fixed - inv_fixed) * 1.0 + growth_rate) * 100),
            xytext=(0.5, ((prof_fixed - inv_fixed) * 1.2 + growth_rate) * 100),
            fontsize=9, fontstyle='italic',
            arrowprops=dict(arrowstyle='->', color='gray'))
ax.yaxis.set_major_formatter(mticker.FormatStrFormatter('%.1f'))
ax.grid(True, alpha=0.3)

plt.suptitle('Dividend Discount Model: The Theoretical Foundation of the FF5 Factors\n'
             r'$r = \frac{Y/B\ -\ \Delta B/B}{M/B} + g$' + '    (simplified Gordon Growth DDM with clean surplus)',
             fontsize=12, y=1.08)
plt.tight_layout()
plt.show()

# --- Summary table ---
print("=" * 70)
print("DDM COMPARATIVE STATICS — NUMERICAL EXAMPLE")
print("=" * 70)
print(f"\nBaseline: Y/B = 12%,  ΔB/B = 5%,  M/B = 1.5,  g = 2%")
baseline_r = (0.12 - 0.05) / 1.5 + 0.02
print(f"  → Implied expected return: {baseline_r:.2%}\n")

scenarios = [
    ("High profitability (Y/B = 18%)",  0.18, 0.05, 1.5),
    ("Low profitability  (Y/B = 6%)",   0.06, 0.05, 1.5),
    ("Conservative inv.  (ΔB/B = 2%)",  0.12, 0.02, 1.5),
    ("Aggressive inv.    (ΔB/B = 10%)", 0.12, 0.10, 1.5),
    ("Deep value stock   (M/B = 0.8)",  0.12, 0.05, 0.8),
    ("Growth stock       (M/B = 3.0)",  0.12, 0.05, 3.0),
]

print(f"{'Scenario':<38} {'Y/B':>6} {'ΔB/B':>6} {'M/B':>6} {'r':>8}")
print("-" * 70)
for label, yb, db, mb in scenarios:
    r = (yb - db) / mb + growth_rate
    delta = r - baseline_r
    arrow = "↑" if delta > 0 else "↓"
    print(f"{label:<38} {yb:>5.0%} {db:>5.0%} {mb:>6.1f} {r:>7.2%} {arrow} ({delta:+.2%})")

print("\n→ These comparative statics motivate the factors under the imposed assumptions.")
print("  The DDM shows they are not ad hoc — they follow from present-value logic.")
print("Illustration only: fixed dividend-growth assumptions are imposed, not estimated.")
print("Treat Y and investment as next-period cash-flow components; sustained changes in reinvestment can also change growth.")

# ============================================================================
# Visualize: Factor risk-return profile
# ============================================================================

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

factor_cols = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA']
ann_means = factors[factor_cols].mean() * 12
ann_stds = factors[factor_cols].std() * np.sqrt(12)
sharpes = ann_means / ann_stds

colors = {'Mkt_RF': '#2c3e50', 'SMB': '#e74c3c', 'HML': '#2ecc71',
          'RMW': '#3498db', 'CMA': '#9b59b6'}

# --- Left: Risk-return scatter ---
ax = axes[0]
for col in factor_cols:
    ax.scatter(ann_stds[col], ann_means[col], s=120, color=colors[col],
               edgecolors='white', linewidth=1.5, zorder=3)
    ax.annotate(col, (ann_stds[col], ann_means[col]),
                textcoords="offset points", xytext=(8, 5), fontsize=10, fontweight='bold',
                color=colors[col])
ax.axhline(y=0, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)
ax.set_xlabel('Annualized Volatility', fontsize=11)
ax.set_ylabel('Annualized Mean Return', fontsize=11)
ax.set_title('Factor Risk-Return Profile', fontweight='bold', fontsize=13)
ax.grid(True, alpha=0.3)

# --- Right: Sharpe ratio bar chart ---
ax2 = axes[1]
bar_colors = [colors[col] for col in factor_cols]
bars = ax2.bar(factor_cols, sharpes, color=bar_colors, edgecolor='white', width=0.55)
ax2.axhline(y=0, color='gray', linestyle='--', linewidth=0.8)
for bar, sr in zip(bars, sharpes):
    ax2.text(bar.get_x() + bar.get_width()/2,
             bar.get_height() + (0.01 if sr >= 0 else -0.04),
             f'{sr:.2f}', ha='center', fontsize=10, fontweight='bold')
ax2.set_ylabel('Sharpe Ratio', fontsize=11)
ax2.set_title('Factor Sharpe Ratios (Annualized)', fontweight='bold', fontsize=13)
ax2.grid(True, alpha=0.3, axis='y')

fig.suptitle('Fama-French 5-Factor Performance Overview', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

# ============================================================================
# 5-Factor Model: Factor Correlation Matrix
# ============================================================================

print("FACTOR CORRELATION MATRIX")
print("="*70)

corr = factors[factor_cols].corr()
print(corr.round(3))

# Visualize
fig, ax = plt.subplots(figsize=(8, 6))
mask = np.triu(np.ones_like(corr, dtype=bool), k=1)
sns.heatmap(corr, annot=True, fmt='.2f', cmap='RdBu_r', center=0,
            mask=mask, square=True, ax=ax, vmin=-1, vmax=1,
            linewidths=0.5)
ax.set_title('Factor Correlation Matrix (2000–latest released month)', fontweight='bold', fontsize=13)
plt.tight_layout()
plt.show()

print(f"HML–CMA correlation: {corr.loc['HML', 'CMA']:.3f}")
print("Style factors are not constructed to be independent or market-neutral.")

# ============================================================================
# Visualize: 3-Factor vs 5-Factor model comparison
# ============================================================================

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

tickers = list(etf_tickers.keys())
labels = [etf_tickers[t] for t in tickers]
short_labels = tickers

# Re-extract coefficients for plotting
coef_data_3f = {}
coef_data_5f = {}
for ticker in tickers:
    edf = etf_dataframes[ticker]
    res3 = run_ff_regression(None, ['Mkt_RF', 'SMB', 'HML'], edf)
    res5 = run_ff_regression(None, ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA'], edf)
    coef_data_3f[ticker] = {'Alpha': res3.params[0], 'Mkt_RF': res3.params[1],
                             'SMB': res3.params[2], 'HML': res3.params[3],
                             'R²': res3.rsquared}
    coef_data_5f[ticker] = {'Alpha': res5.params[0], 'Mkt_RF': res5.params[1],
                             'SMB': res5.params[2], 'HML': res5.params[3],
                             'RMW': res5.params[4], 'CMA': res5.params[5],
                             'R²': res5.rsquared}

# --- Left: R² comparison ---
ax = axes[0]
x_pos = np.arange(len(tickers))
width = 0.3
bars1 = ax.bar(x_pos - width/2, [coef_data_3f[t]['R²'] for t in tickers],
               width, label='3-Factor', color='#3498db', edgecolor='white')
bars2 = ax.bar(x_pos + width/2, [coef_data_5f[t]['R²'] for t in tickers],
               width, label='5-Factor', color='#e67e22', edgecolor='white')
for bar in list(bars1) + list(bars2):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.003,
            f'{bar.get_height():.3f}', ha='center', fontsize=9, fontweight='bold')
ax.set_xticks(x_pos)
ax.set_xticklabels(short_labels, fontsize=11)
ax.set_ylabel('R²', fontsize=11)
ax.set_title('Model Fit: 3-Factor vs 5-Factor', fontweight='bold', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis='y')
ax.set_ylim(0, 1.08)

# --- Right: Coefficient comparison for shared factors ---
ax2 = axes[1]
shared_factors = ['Mkt_RF', 'SMB', 'HML']
x_pos2 = np.arange(len(shared_factors))
group_width = 0.8
bar_w = group_width / len(tickers) / 2  # 2 models per ticker

color_3f = ['#3498db', '#2980b9', '#1a5276']
color_5f = ['#e67e22', '#d35400', '#a04000']

for i, ticker in enumerate(tickers):
    offset_3f = -group_width/2 + i * (group_width/len(tickers)) + bar_w * 0.1
    offset_5f = offset_3f + bar_w
    vals_3f = [coef_data_3f[ticker][f] for f in shared_factors]
    vals_5f = [coef_data_5f[ticker][f] for f in shared_factors]
    ax2.bar(x_pos2 + offset_3f, vals_3f, bar_w * 0.9, color=color_3f[i], alpha=0.7,
            label=f'{ticker} (3F)' if i == 0 else f'{ticker} (3F)')
    ax2.bar(x_pos2 + offset_5f, vals_5f, bar_w * 0.9, color=color_5f[i], alpha=0.7,
            label=f'{ticker} (5F)' if i == 0 else f'{ticker} (5F)')

ax2.set_xticks(x_pos2)
ax2.set_xticklabels(shared_factors, fontsize=11)
ax2.set_ylabel('Coefficient', fontsize=11)
ax2.set_title('How Betas Change: 3F → 5F', fontweight='bold', fontsize=13)
ax2.axhline(y=0, color='gray', linestyle='--', linewidth=0.8)
ax2.legend(fontsize=7, ncol=3, loc='upper right')
ax2.grid(True, alpha=0.3, axis='y')

fig.suptitle('3-Factor vs 5-Factor Model Comparison', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

print("Changes in asset HML coefficients do not test whether HML itself is spanned by the other factors.")

# ============================================================================
# Momentum Factor Analysis
# ============================================================================

print("="*70)
print("MOMENTUM FACTOR (UMD) ANALYSIS")
print("="*70)

# Cumulative returns of all factors
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Cumulative factor returns
ax = axes[0, 0]
for col, color, lbl in [('Mkt_RF', 'darkblue', 'Market'),
                          ('SMB', 'teal', 'SMB (Size)'),
                          ('HML', 'purple', 'HML (Value)'),
                          ('RMW', 'darkgreen', 'RMW (Profit.)'),
                          ('CMA', 'brown', 'CMA (Invest.)'),
                          ('UMD', 'red', 'UMD (Momentum)')]:
    cum = (1 + factors[col]).cumprod()
    ax.plot(factors.index, cum, label=lbl, linewidth=1.5, color=color)
ax.set_title('Cumulative Factor Returns (2000–latest released month)', fontweight='bold')
ax.set_ylabel('Hypothetical return index (start=1)')
ax.legend(fontsize=8, loc='upper left')
ax.set_yscale('log')
ax.grid(True, alpha=0.3)

# Plot 2: Momentum monthly returns
ax = axes[0, 1]
ax.bar(factors.index, factors['UMD'] * 100, color='red', alpha=0.6, width=25)
ax.set_title('Monthly Momentum Returns (UMD)', fontweight='bold')
ax.set_ylabel('Return (%)')
ax.grid(True, alpha=0.3)

# Plot 3: Rolling 12-month momentum return
ax = axes[1, 0]
rolling_mom = ((1 + factors['UMD']).rolling(12).apply(np.prod, raw=True) - 1) * 100
ax.plot(factors.index, rolling_mom, color='red', linewidth=1.5)
ax.axhline(0, color='black', linewidth=0.5, linestyle='--')
ax.fill_between(factors.index, rolling_mom, 0,
                where=rolling_mom < 0, color='red', alpha=0.3, label='Negative trailing returns')
ax.set_title('Rolling 12-Month Momentum Return', fontweight='bold')
ax.set_ylabel('Cumulative 12M Return (%)')
ax.legend()
ax.grid(True, alpha=0.3)

# Plot 4: Distribution of momentum returns vs market
ax = axes[1, 1]
ax.hist(factors['UMD'] * 100, bins=40, alpha=0.6, color='red', label='UMD', density=True)
ax.hist(factors['Mkt_RF'] * 100, bins=40, alpha=0.4, color='blue', label='Mkt-RF', density=True)
ax.set_title('Distribution: Momentum vs Market', fontweight='bold')
ax.set_xlabel('Monthly Return (%)')
ax.set_ylabel('Density')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Momentum statistics
print(f"\nMomentum (UMD) Summary Statistics:")
print(f"  Mean monthly return:    {factors['UMD'].mean()*100:.2f}%")
print(f"  Monthly std dev:        {factors['UMD'].std()*100:.2f}%")
print(f"  Annualized return:      {factors['UMD'].mean()*1200:.2f}%")
print(f"  Annualized volatility:  {factors['UMD'].std()*np.sqrt(12)*100:.2f}%")
print(f"  Sharpe ratio:           {factors['UMD'].mean()/factors['UMD'].std()*np.sqrt(12):.2f}")
print(f"  Skewness:               {factors['UMD'].skew():.2f}")
print(f"  Kurtosis:               {factors['UMD'].kurtosis():.2f}")
print(f"  Worst month:            {factors['UMD'].min()*100:.2f}%")
print(f"  Best month:             {factors['UMD'].max()*100:.2f}%")

print('Inspect the reported skewness, tail outcomes and sample mean. Historical crash exposure does not establish a future premium.')

# ============================================================================
# Profitability Factor Deep Dive
# ============================================================================

print("="*70)
print("PROFITABILITY FACTOR (RMW) DEEP DIVE")
print("="*70)

# RMW statistics
rmw = factors['RMW']
print(f"\nRMW (Robust Minus Weak) Summary:")
print(f"  Mean monthly:       {rmw.mean()*100:.2f}%")
print(f"  Annualized return:  {rmw.mean()*1200:.2f}%")
print(f"  Ann. volatility:    {rmw.std()*np.sqrt(12)*100:.2f}%")
print(f"  Sharpe ratio:       {rmw.mean()/rmw.std()*np.sqrt(12):.2f}")
print(f"  HAC z (mean=0):    {mean_hac(rmw).tvalues[0]:.2f}")
print(f"  Skewness:           {rmw.skew():.2f}")
print(f"  Correlation w/ HML: {factors['RMW'].corr(factors['HML']):.3f}")
print(f"  Correlation w/ UMD: {factors['RMW'].corr(factors['UMD']):.3f}")

# Interaction: Value × Profitability
# Compare the measured HML/RMW correlation and historical returns
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Plot 1: RMW vs HML scatter
ax = axes[0]
ax.scatter(factors['HML'] * 100, factors['RMW'] * 100, alpha=0.4, s=30)
z = np.polyfit(factors['HML'], factors['RMW'], 1)
p = np.poly1d(z)
x_range = np.linspace(factors['HML'].min(), factors['HML'].max(), 100)
ax.plot(x_range * 100, p(x_range) * 100, 'r-', linewidth=2)
ax.set_xlabel('HML Return (%)')
ax.set_ylabel('RMW Return (%)')
ax.set_title('Value vs Profitability\n(Sample correlation, not a diversification guarantee)', fontweight='bold')
ax.grid(True, alpha=0.3)
corr_val = factors['HML'].corr(factors['RMW'])
ax.annotate(f'r = {corr_val:.2f}', xy=(0.05, 0.95), xycoords='axes fraction',
            fontsize=12, fontweight='bold', va='top')

# Plot 2: Cumulative returns of HML, RMW, and combined
ax = axes[1]
combined = (factors['HML'] + factors['RMW']) / 2  # Equal-weight combination
cum_hml = (1 + factors['HML']).cumprod()
cum_rmw = (1 + factors['RMW']).cumprod()
cum_combined = (1 + combined).cumprod()
ax.plot(factors.index, cum_hml, label='HML (Value)', color='purple', linewidth=1.5)
ax.plot(factors.index, cum_rmw, label='RMW (Profitability)', color='darkgreen', linewidth=1.5)
ax.plot(factors.index, cum_combined, label='50/50 Blend', color='black', linewidth=2, linestyle='--')
ax.set_title('Cumulative Returns:\nValue, Profitability, and Blend', fontweight='bold')
ax.set_ylabel('Hypothetical return index (start=1)')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Plot 3: Rolling correlation between HML and RMW
ax = axes[2]
rolling_corr = factors['HML'].rolling(36).corr(factors['RMW'])
ax.plot(factors.index, rolling_corr, color='darkred', linewidth=1.5)
ax.axhline(0, color='black', linewidth=0.5, linestyle='--')
ax.set_title('Rolling 36-Month Correlation\nHML vs RMW', fontweight='bold')
ax.set_ylabel('Correlation')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print('Read the measured full-sample and rolling correlations. Correlations of characteristics, factor spreads and fund returns are different quantities.')

# ============================================================================
# Reversal Factors Analysis
# ============================================================================

print("="*70)
print("SHORT-TERM AND LONG-TERM REVERSAL FACTORS")
print("="*70)

for col, name in [('ST_Rev', 'Short-Term Reversal'), ('LT_Rev', 'Long-Term Reversal')]:
    r = factors[col]
    tstat = mean_hac(r).tvalues[0]
    print(f"\n{name} ({col}):")
    print(f"  Mean monthly:       {r.mean()*100:.2f}%")
    print(f"  Annualized return:  {r.mean()*1200:.2f}%")
    print(f"  Ann. volatility:    {r.std()*np.sqrt(12)*100:.2f}%")
    print(f"  Sharpe ratio:       {r.mean()/r.std()*np.sqrt(12):.2f}")
    print(f"  t-stat:             {tstat:.2f}")

# Correlation with other factors
print("\n\nCORRELATION OF REVERSAL FACTORS WITH STANDARD FACTORS:")
rev_corr_cols = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD', 'ST_Rev', 'LT_Rev']
print(factors[rev_corr_cols].corr()[['ST_Rev', 'LT_Rev']].round(3))

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Cumulative reversal factor returns
ax = axes[0, 0]
cum_st = (1 + factors['ST_Rev']).cumprod()
cum_lt = (1 + factors['LT_Rev']).cumprod()
cum_mom = (1 + factors['UMD']).cumprod()
ax.plot(factors.index, cum_st, label='ST Reversal', color='orange', linewidth=1.5)
ax.plot(factors.index, cum_lt, label='LT Reversal', color='darkblue', linewidth=1.5)
ax.plot(factors.index, cum_mom, label='Momentum (UMD)', color='red', linewidth=1.5, linestyle='--')
ax.set_title('Reversal vs Momentum: Cumulative Returns', fontweight='bold')
ax.set_ylabel('Hypothetical return index (start=1)')
ax.legend()
ax.grid(True, alpha=0.3)

# Plot 2: Monthly return comparison
ax = axes[0, 1]
ax.scatter(factors['UMD'] * 100, factors['ST_Rev'] * 100, alpha=0.4, s=30, color='orange')
ax.set_xlabel('Momentum (UMD) Return (%)')
ax.set_ylabel('ST Reversal Return (%)')
corr_sm = factors['UMD'].corr(factors['ST_Rev'])
ax.set_title(f'Momentum vs ST Reversal (r = {corr_sm:.2f})', fontweight='bold')
ax.grid(True, alpha=0.3)
ax.axhline(0, color='gray', linewidth=0.5)
ax.axvline(0, color='gray', linewidth=0.5)

# Plot 3: LT Reversal vs HML
ax = axes[1, 0]
ax.scatter(factors['HML'] * 100, factors['LT_Rev'] * 100, alpha=0.4, s=30, color='purple')
ax.set_xlabel('HML (Value) Return (%)')
ax.set_ylabel('LT Reversal Return (%)')
corr_lh = factors['HML'].corr(factors['LT_Rev'])
ax.set_title(f'Value vs LT Reversal (r = {corr_lh:.2f})', fontweight='bold')
ax.grid(True, alpha=0.3)
ax.axhline(0, color='gray', linewidth=0.5)
ax.axvline(0, color='gray', linewidth=0.5)

# Plot 4: Autocorrelation of market returns at different lags
# Aggregate market autocorrelation is distinct from cross-sectional reversal sorts.
ax = axes[1, 1]
mkt_returns = factors['Mkt_RF']
lags = list(range(1, 61))  # 1 to 60 months
autocorrs = [mkt_returns.autocorr(lag=l) for l in lags]
colors_ac = ['orange' if l <= 1 else ('red' if l <= 12 else 'darkblue') for l in lags]
ax.bar(lags, autocorrs, color=colors_ac, alpha=0.7, width=0.8)
ax.axhline(0, color='black', linewidth=0.5)
# Add significance band (approximate 95% CI for white noise)
n_obs = len(mkt_returns)
sig_bound = 1.96 / np.sqrt(n_obs)
ax.axhline(sig_bound, color='gray', linewidth=0.8, linestyle='--', alpha=0.6)
ax.axhline(-sig_bound, color='gray', linewidth=0.8, linestyle='--', alpha=0.6)
ax.set_xlabel('Lag (months)')
ax.set_ylabel('Autocorrelation')
ax.set_title('Market Return Autocorrelation by Lag\n(orange=1mo, red=2-12mo, blue=13-60mo)', fontweight='bold')
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print('Read the reported reversal correlations. Aggregate market autocorrelation is not a test of cross-sectional momentum or reversal.')

# ============================================================================
# The Beta Anomaly: Empirical Security Market Line
# ============================================================================

print("="*70)
print("THE BETA ANOMALY — FLAT SECURITY MARKET LINE")
print("="*70)

# ---- Download beta-sorted decile portfolios from Ken French ----
# This file contains monthly returns for 10 portfolios sorted by pre-ranking
# market beta. We'll use value-weighted returns.

download_beta_portfolios = load_beta_portfolios

print("Downloading beta-sorted decile portfolios from Ken French...")
beta_ports = download_beta_portfolios()

# Filter to our analysis period
beta_ports = beta_ports.loc['2000-01-01':]
beta_ports = beta_ports.loc[beta_ports.index.intersection(factors.index)]
print(f"Beta-sorted decile portfolios: {beta_ports.shape[0]} months")
print(f"Deciles: D1 (lowest beta) through D10 (highest beta)\n")

# ---- Compute average excess returns and realized betas ----
rf_aligned = factors.loc[beta_ports.index, 'RF']
mkt_aligned = factors.loc[beta_ports.index, 'Mkt_RF']

avg_excess_returns = []
realized_betas = []
decile_labels = []

print(f"{'Decile':<10s} {'Avg Excess Ret (ann %)':>22s} {'Realized Beta':>14s} "
      f"{'Volatility (ann %)':>20s} {'Sharpe Ratio':>13s}")
print("-" * 82)

for col in beta_ports.columns:
    excess = beta_ports[col] - rf_aligned
    
    # Realized beta from regression on market
    X = sm.add_constant(mkt_aligned.values)
    res = sm.OLS(excess.values, X).fit()
    beta_hat = res.params[1]
    
    ann_ret = excess.mean() * 12 * 100
    ann_vol = excess.std() * np.sqrt(12) * 100
    sharpe = excess.mean() / excess.std() * np.sqrt(12)
    
    avg_excess_returns.append(ann_ret)
    realized_betas.append(beta_hat)
    decile_labels.append(col)
    
    print(f"{col:<10s} {ann_ret:22.2f} {beta_hat:14.2f} {ann_vol:20.2f} {sharpe:13.2f}")

avg_excess_returns = np.array(avg_excess_returns)
realized_betas = np.array(realized_betas)

# ---- The Punchline: Plot the Empirical SML ----
fig, axes = plt.subplots(1, 2, figsize=(16, 7))

# Left panel: Empirical SML vs CAPM prediction
ax = axes[0]
ax.scatter(realized_betas, avg_excess_returns, s=120, c='darkblue', zorder=5,
           edgecolors='white', linewidth=1.5)

for i, label in enumerate(decile_labels):
    ax.annotate(label, (realized_betas[i], avg_excess_returns[i]),
                textcoords="offset points", xytext=(8, 5), fontsize=9,
                fontweight='bold', color='darkblue')

# Theoretical CAPM line: passes through (0, 0) and (1, market premium)
mkt_premium_ann = mkt_aligned.mean() * 12 * 100
beta_range = np.linspace(0, max(realized_betas) * 1.1, 100)
capm_line = beta_range * mkt_premium_ann
ax.plot(beta_range, capm_line, 'r--', linewidth=2, label='CAPM Prediction',
        alpha=0.8)

# Best-fit line through the data
slope, intercept = np.polyfit(realized_betas, avg_excess_returns, 1)
fit_line = slope * beta_range + intercept
ax.plot(beta_range, fit_line, 'b-', linewidth=2.5,
        label=f'Empirical Fit (slope={slope:.1f})', alpha=0.9)

ax.set_xlabel('Realized Market Beta', fontsize=12, fontweight='bold')
ax.set_ylabel('Average Annualized Excess Return (%)', fontsize=12, fontweight='bold')
ax.set_title('The Empirical Security Market Line\n(Beta-Sorted Decile Portfolios, 2000–latest released month)',
             fontsize=13, fontweight='bold')
ax.legend(fontsize=11, loc='upper left')
ax.grid(True, alpha=0.3)

# Right panel: Sharpe ratio by decile
ax = axes[1]
sharpe_ratios = []
for col in beta_ports.columns:
    excess = beta_ports[col] - rf_aligned
    sharpe_ratios.append(excess.mean() / excess.std() * np.sqrt(12))

colors = plt.cm.RdYlGn_r(np.linspace(0.1, 0.9, 10))
bars = ax.bar(range(1, 11), sharpe_ratios, color=colors, edgecolor='black',
              linewidth=0.8)
ax.set_xlabel('Beta Decile (1=Lowest, 10=Highest)', fontsize=12, fontweight='bold')
ax.set_ylabel('Sharpe Ratio', fontsize=12, fontweight='bold')
ax.set_title('Risk-Adjusted Returns by Beta Decile\n(Higher is Better)',
             fontsize=13, fontweight='bold')
ax.set_xticks(range(1, 11))
ax.axhline(0, color='black', linewidth=0.5)
ax.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for bar, sr in zip(bars, sharpe_ratios):
    ax.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.01,
            f'{sr:.2f}', ha='center', va='bottom', fontsize=9, fontweight='bold')

plt.tight_layout()
plt.show()

# ---- Construct a simple BAB-like factor ----
print("\n" + "="*70)
print("SIMPLE BAB FACTOR: LONG LOW-BETA (D1) − SHORT HIGH-BETA (D10)")
print("="*70)

bab_simple = (beta_ports['D1'] - rf_aligned) - (beta_ports['D10'] - rf_aligned)
# More properly: BAB = (1/beta_low)(D1 - RF) - (1/beta_high)(D10 - RF)
# but for illustration the simple spread is instructive

# Estimate each beta from the prior 60 months, with a one-month information lag.
# This decile proxy is not the stock-level Frazzini–Pedersen BAB replication.
low_excess = beta_ports['D1'] - rf_aligned
high_excess = beta_ports['D10'] - rf_aligned
rolling_variance = mkt_aligned.rolling(60).var()
beta_low = (low_excess.rolling(60).cov(mkt_aligned) / rolling_variance).shift(1)
beta_high = (high_excess.rolling(60).cov(mkt_aligned) / rolling_variance).shift(1)
valid = (beta_low > 0.1) & (beta_high > 0.1)
bab_leveraged = (low_excess / beta_low - high_excess / beta_high).where(valid).dropna()
coverage('BAB-like proxy (prior 60-month beta scaling)', bab_leveraged)
bab_simple_comparable = bab_simple.loc[bab_leveraged.index]

for name, ret in [('Simple (D1 − D10)', bab_simple_comparable), 
                   ('Lagged-beta decile proxy', bab_leveraged)]:
    ann_ret = ret.mean() * 12 * 100
    ann_vol = ret.std() * np.sqrt(12) * 100
    sharpe = ret.mean() / ret.std() * np.sqrt(12)
    tstat = mean_hac(ret).tvalues[0]
    cum = (1 + ret).cumprod()
    maxdd = (cum / cum.cummax().clip(lower=1) - 1).min() * 100
    print(f"\n  {name}:")
    print(f"    Ann. Return:   {ann_ret:7.2f}%")
    print(f"    Ann. Vol:      {ann_vol:7.2f}%")
    print(f"    Sharpe Ratio:  {sharpe:7.2f}")
    print(f"    t-statistic:   {tstat:7.2f}")
    print(f"    Max Drawdown:  {maxdd:7.1f}%")

# ---- Alpha test: regress BAB on FF5 + Momentum ----
print("\n" + "="*70)
print("DOES BAB HAVE ALPHA AFTER FF5 + MOMENTUM?")
print("="*70)

bab_aligned = bab_leveraged.reindex(factors.index).dropna()
ff_aligned = factors.loc[bab_aligned.index, ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD']]

X_ff = sm.add_constant(ff_aligned.values)
res_bab = sm.OLS(bab_aligned.values, X_ff).fit(
    cov_type='HAC', cov_kwds={'maxlags': 6, 'use_correction': True}, use_t=False
)

coef_names = ['Alpha'] + list(ff_aligned.columns)
print(f"\n{'Factor':<10s} {'Coef':>8s} {'t-stat':>8s} {'p-value':>9s}")
print("-" * 38)
for i, name in enumerate(coef_names):
    print(f"{name:<10s} {res_bab.params[i]:8.4f} {res_bab.tvalues[i]:8.2f} "
          f"{res_bab.pvalues[i]:9.4f}")
print(f"\nR² = {res_bab.rsquared:.3f}")
print(f"Alpha (annualized) = {res_bab.params[0] * 12 * 100:.2f}%")

# ---- Cumulative BAB vs individual factors ----
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

ax = axes[0]
cum_bab = (1 + bab_leveraged).cumprod()
cum_mkt = (1 + factors.loc[bab_leveraged.index, 'Mkt_RF']).cumprod()
cum_hml = (1 + factors.loc[bab_leveraged.index, 'HML']).cumprod()
cum_rmw = (1 + factors.loc[bab_leveraged.index, 'RMW']).cumprod()

ax.plot(cum_bab.index, cum_bab, linewidth=2.5, label='BAB (Leveraged)', color='darkblue')
ax.plot(cum_mkt.index, cum_mkt, linewidth=1.5, label='Mkt-RF', color='red', alpha=0.7)
ax.plot(cum_hml.index, cum_hml, linewidth=1.5, label='HML', color='green', alpha=0.7)
ax.plot(cum_rmw.index, cum_rmw, linewidth=1.5, label='RMW', color='purple', alpha=0.7)
ax.set_title('Cumulative Returns: BAB vs Major Factors', fontsize=13, fontweight='bold')
ax.set_ylabel('Hypothetical return index (start=1)')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Drawdown of BAB
ax = axes[1]
dd_bab = (cum_bab / cum_bab.cummax().clip(lower=1) - 1) * 100
ax.fill_between(dd_bab.index, dd_bab, 0, alpha=0.4, color='darkblue')
ax.plot(dd_bab.index, dd_bab, linewidth=1, color='darkblue')
ax.set_title('BAB Factor Drawdowns', fontsize=13, fontweight='bold')
ax.set_ylabel('Drawdown (%)')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Empirical cross-sectional slope {slope:.2f}; sample CAPM slope {mkt_premium_ann:.2f}.")
print(f"Proxy HAC alpha p-value: {res_bab.pvalues[0]:.4f}.")
print("This descriptive decile proxy omits trading, stock-borrow and financing spreads.")
print("Predetermined beta scaling avoids full-sample lookahead, but estimated exposures")
print("can change and the spread need not be market neutral out of sample.")

# ============================================================================
# Research Deep Dive: Factor Performance Across Decades
# ============================================================================

print("="*70)
print("FACTOR PERFORMANCE BY DECADE")
print("="*70)

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

# Define sub-periods
periods = {
    '2000-2004': ('2000-01-01', '2004-12-31'),
    '2005-2009': ('2005-01-01', '2009-12-31'),
    '2010-2014': ('2010-01-01', '2014-12-31'),
    '2015-2019': ('2015-01-01', '2019-12-31'),
    '2020-2024': ('2020-01-01', '2024-12-31'),
    f'2025–{factors.index[-1]:%Y-%m} (partial)': ('2025-01-01', AS_OF.strftime('%Y-%m-%d')),
}

decade_returns = {}
for period_name, (start, end) in periods.items():
    sub = factors.loc[start:end, factor_cols_all]
    ann_ret = sub.mean() * 12 * 100  # Annualized, in %
    decade_returns[period_name] = ann_ret

decade_df = pd.DataFrame(decade_returns).T
print("\nAnnualized Returns by Sub-Period (%):")
print(decade_df.round(2))

# Visualization
fig, axes = plt.subplots(2, 1, figsize=(14, 10))

# Plot 1: Heatmap of factor returns by period
ax = axes[0]
sns.heatmap(decade_df, annot=True, fmt='.1f', cmap='RdYlGn', center=0,
            ax=ax, linewidths=0.5, cbar_kws={'label': 'Annualized Return (%)'})
ax.set_title('Factor Returns by Sub-Period (%)\n'
             'Green = positive return, Red = negative return', fontweight='bold', fontsize=12)
ax.set_ylabel('')

# Plot 2: HML return index; this is not a valuation spread
ax = axes[1]
# Show the value factor's "death and resurrection"
hml_cum = (1 + factors['HML']).cumprod()
ax.plot(factors.index, hml_cum, color='purple', linewidth=2, label='HML Cumulative')
ax.set_title('The Value Factor: Drawdown and Recovery (2000–latest released month)', fontweight='bold', fontsize=12)
ax.set_ylabel('Hypothetical return index (start=1)')
ax.grid(True, alpha=0.3)

# Annotate key periods
ax.axvspan('2018-01-01', '2020-09-01', alpha=0.15, color='red', label='Value drawdown')
ax.axvspan('2020-10-01', '2022-12-31', alpha=0.15, color='green', label='Value recovery')
ax.legend()

plt.tight_layout()
plt.show()

print("Compare realized means across periods; the latest block is incomplete.")
print("Subperiod differences are not by themselves a structural-break test or forecast.")

# ============================================================================
# The "Kitchen Sink" Model: 5 Factors + Momentum on Our ETFs
# ============================================================================

print("="*70)
print("COMPREHENSIVE MODEL: FF5 + MOMENTUM")
print("="*70)
print("\nAdding UMD to the 5-factor model gives us the broadest")
print("specification considered here. Let's see how each ETF maps onto all six factors.\n")

full_cols = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD']
full_names = ['Alpha', 'Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD']

all_params = {}
all_pvals = {}
all_r2 = {}

for ticker, label in etf_tickers.items():
    edf = etf_dataframes[ticker]
    res = run_ff_regression(None, full_cols, edf, stock_name=label)
    all_params[ticker] = res.params
    all_pvals[ticker] = res.pvalues
    all_r2[ticker] = res.rsquared

    print(f"\n{label}")
    print("-" * 60)
    for i, name in enumerate(full_names):
        sig = '***' if res.pvalues[i] < 0.01 else ('**' if res.pvalues[i] < 0.05 else ('*' if res.pvalues[i] < 0.10 else ''))
        print(f"  {name:10s} = {res.params[i]:9.5f}  (t={res.tvalues[i]:6.2f}, p={res.pvalues[i]:.4f}) {sig}")
    print(f"  R² = {res.rsquared:.4f}   Adj R² = {res.rsquared_adj:.4f}")

# Visualization: Factor loadings comparison
fig, ax = plt.subplots(figsize=(12, 6))
x = np.arange(len(full_names) - 1)  # Exclude alpha for visualization
width = 0.25

for i, (ticker, label) in enumerate(etf_tickers.items()):
    params_no_alpha = all_params[ticker][1:]  # Skip alpha
    bars = ax.bar(x + i * width, params_no_alpha, width, label=label, alpha=0.8)
    # Add significance markers
    for j, bar in enumerate(bars):
        p = all_pvals[ticker][j + 1]
        if p < 0.01:
            ax.text(bar.get_x() + bar.get_width()/2, bar.get_height(),
                    '***', ha='center', va='bottom', fontsize=8)

ax.set_xlabel('Factor')
ax.set_ylabel('Factor Loading (β)')
ax.set_title('Factor Loadings: FF5 + Momentum\n(*** = significant at 1% level, Newey-West SEs)',
             fontweight='bold', fontsize=12)
ax.set_xticks(x + width)
ax.set_xticklabels(full_names[1:])
ax.axhline(0, color='black', linewidth=0.5)
ax.legend()
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print("""
INTERPRETATION GUIDE:
- Bar height = factor sensitivity (positive = tilts toward that factor)
- *** = statistically significant at 1% level
- SPY should be dominated by Mkt_RF ≈ 1, everything else near zero
- IWN should show positive SMB, positive HML, and possibly negative UMD
- VTV should show positive HML, possibly positive RMW
""")

# Hypothetical continuous-time lognormal model, annual time steps.
rng = np.random.default_rng(42)
n_simulations, years = 20_000, 30
market_drift, market_vol, tracking_error = 0.08, 0.16, 0.04
market_shocks = rng.normal(size=(n_simulations, years))
tracking_shocks = rng.normal(size=(n_simulations, years))
market_log = (market_drift - market_vol**2 / 2) + market_vol * market_shocks
market_terminal = np.exp(market_log.sum(axis=1))
scenarios = [-0.01, 0.0, 0.01, 0.02]
fig, axes = plt.subplots(1, 4, figsize=(18, 4.5), sharey=True)
for ax, net_premium in zip(axes, scenarios):
    # Shared market path; independent tracking component only.
    relative_log = net_premium - tracking_error**2 / 2 + tracking_error * tracking_shocks
    tilted_terminal = np.exp((market_log + relative_log).sum(axis=1))
    probability = np.mean(tilted_terminal > market_terminal)
    analytical = stats.norm.cdf((net_premium - tracking_error**2 / 2) * np.sqrt(years) / tracking_error)
    assert abs(probability - analytical) < 0.02
    ax.hist(np.log10(tilted_terminal / market_terminal), bins=50, density=True, alpha=0.7)
    ax.axvline(0, color='black', linestyle='--')
    ax.set_title(f'Assumed net premium {net_premium:+.0%}\nP(outperform) = {probability:.1%}')
    ax.set_xlabel('Log10(relative terminal wealth)')
    print(f'Net drift {net_premium:+.1%}: simulated {probability:.2%}; analytical {analytical:.2%}')
axes[0].set_ylabel('Density')
fig.suptitle('Hypothetical 30-year relative performance — shared market paths')
plt.tight_layout()
plt.show()
print('These probabilities depend entirely on imposed parameters, not estimates from ETF data.')

# Rolling estimates use only the window ending at the labelled month.
window = 36
roll_cols = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD']
roll_names = ['Alpha monthly'] + roll_cols
for ticker, label in etf_tickers.items():
    edf = etf_dataframes[ticker]
    estimates, bounds, dates = [], [], []
    for i in range(window, len(edf)+1):
        d = edf.iloc[i-window:i]
        res = fit_hac(d.Stock_Excess_Return, d[roll_cols])
        estimates.append(res.params)
        bounds.append(res.conf_int())
        dates.append(edf.index[i-1])
    assert dates[-1] == edf.index[-1]
    estimates, bounds = np.asarray(estimates), np.asarray(bounds)
    fig, axes = plt.subplots(4, 2, figsize=(13, 11), sharex=True)
    for j, name in enumerate(roll_names):
        ax = axes.flat[j]
        ax.plot(dates, estimates[:,j], color='darkblue', linewidth=1)
        ax.fill_between(dates, bounds[:,j,0], bounds[:,j,1], alpha=.15, color='steelblue')
        ax.axhline(0, color='gray', linewidth=.6)
        ax.set_title(name)
    axes.flat[-1].set_visible(False)
    fig.suptitle(f'{label}: rolling {window}-month coefficients and pointwise 95% HAC intervals')
    fig.tight_layout()
    plt.show()
print('Overlapping windows are dependent; these exploratory intervals are not simultaneous bands or a structural-break test.')
print('Apparent beta shifts combine estimation noise and possible exposure changes. Passive funds can have nonzero model-relative alpha.')