Practical Factor Investing: ETF Regression Examples

# ============================================================================
# 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
from statsmodels.stats.diagnostic import acorr_ljungbox
import yfinance as yf
import urllib.request
import zipfile
import tempfile
import os
import warnings
warnings.filterwarnings('ignore', category=FutureWarning)

try:
    plt.style.use('seaborn-v0_8-darkgrid')
except:
    plt.style.use('default')
sns.set_palette("husl")
np.random.seed(42)

print("Libraries imported successfully!")

# ============================================================================
# Step 1: Download Fama-French Factors (DM Daily + DM Monthly + EM Monthly)
# ============================================================================

def load_ff_csv_from_zip(url, freq='auto'):
    """Download and parse a Fama-French CSV from a zip file.
    
    Parameters
    ----------
    url : str
        URL to the zip file on Kenneth French's data library.
    freq : str
        'daily', 'monthly', or 'auto' (tries daily first, then monthly).
    """
    with tempfile.TemporaryDirectory() as tmpdir:
        zip_path = os.path.join(tmpdir, 'data.zip')
        urllib.request.urlretrieve(url, zip_path)
        with zipfile.ZipFile(zip_path, 'r') as z:
            csv_name = [f for f in z.namelist() if f.endswith('.CSV') or f.endswith('.csv')][0]
            with z.open(csv_name) as f:
                lines = f.read().decode('utf-8').splitlines()
    
    # Find the header row (contains 'Mkt-RF' or 'WML' or 'Mom')
    header_idx = None
    for i, line in enumerate(lines):
        if 'Mkt-RF' in line or 'WML' in line or 'Mom' in line:
            header_idx = i
            break
    
    if header_idx is None:
        raise ValueError("Could not find header row in CSV")
    
    # Parse from header row
    from io import StringIO
    data_text = '\n'.join(lines[header_idx:])
    df = pd.read_csv(StringIO(data_text), index_col=0)
    
    # Clean: remove non-numeric rows and annual data section
    df.index = df.index.astype(str).str.strip()
    
    # Detect frequency: YYYYMMDD (daily) or YYYYMM (monthly)
    daily_rows = df[df.index.str.match(r'^\d{8}$')]
    monthly_rows = df[df.index.str.match(r'^\d{6}$')]
    
    if freq == 'daily' and len(daily_rows) > 0:
        df = daily_rows
        df.index = pd.to_datetime(df.index, format='%Y%m%d')
    elif freq == 'monthly' and len(monthly_rows) > 0:
        df = monthly_rows
        df.index = pd.to_datetime(df.index, format='%Y%m')
    elif len(daily_rows) > 0:
        df = daily_rows
        df.index = pd.to_datetime(df.index, format='%Y%m%d')
    elif len(monthly_rows) > 0:
        df = monthly_rows
        df.index = pd.to_datetime(df.index, format='%Y%m')
    else:
        raise ValueError("No valid date rows found")
    
    df = df.apply(pd.to_numeric, errors='coerce')
    df = df.dropna()
    
    return df

# ---- Developed Markets Factors (Daily — kept for descriptive stats) ----
print("Loading Developed Markets Fama-French 5 Factors (Daily)...")
ff5_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/Developed_5_Factors_Daily_CSV.zip'
ff5 = load_ff_csv_from_zip(ff5_url, freq='daily')
ff5.columns = ff5.columns.str.strip()
ff5 = ff5.rename(columns={'Mkt-RF': 'Mkt_RF'})

print(f"  Shape: {ff5.shape}")
print(f"  Date range: {ff5.index[0].date()} to {ff5.index[-1].date()}")

print("\nLoading Developed Markets Momentum Factor (Daily)...")
mom_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/Developed_Mom_Factor_Daily_CSV.zip'
mom = load_ff_csv_from_zip(mom_url, freq='daily')
mom.columns = mom.columns.str.strip()
mom_col = mom.columns[0]  # Usually 'WML'
mom = mom.rename(columns={mom_col: 'UMD'})

# Merge Developed FF5 + Momentum (daily — used for descriptive stats & cumulative charts)
factors = ff5.join(mom[['UMD']], how='inner')
print(f"  Combined DM daily factors: {factors.shape}")

# ---- Developed Markets Factors (Monthly — for AVWS/AVWC/IWMO regressions) ----
print("\n" + "-"*70)
print("Loading Developed Markets Fama-French 5 Factors (Monthly)...")
dm_ff5_m_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/Developed_5_Factors_CSV.zip'
dm_ff5_m = load_ff_csv_from_zip(dm_ff5_m_url, freq='monthly')
dm_ff5_m.columns = dm_ff5_m.columns.str.strip()
dm_ff5_m = dm_ff5_m.rename(columns={'Mkt-RF': 'Mkt_RF'})

print(f"  Shape: {dm_ff5_m.shape}")
print(f"  Date range: {dm_ff5_m.index[0].date()} to {dm_ff5_m.index[-1].date()}")
print(f"  Columns: {list(dm_ff5_m.columns)}")

print("\nLoading Developed Markets Momentum Factor (Monthly)...")
dm_mom_m_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/Developed_Mom_Factor_CSV.zip'
dm_mom_m = load_ff_csv_from_zip(dm_mom_m_url, freq='monthly')
dm_mom_m.columns = dm_mom_m.columns.str.strip()
dm_mom_m_col = dm_mom_m.columns[0]
dm_mom_m = dm_mom_m.rename(columns={dm_mom_m_col: 'UMD'})

# Merge Developed FF5 + Momentum (monthly)
dm_factors_monthly = dm_ff5_m.join(dm_mom_m[['UMD']], how='inner')
print(f"  Combined DM monthly factors (FF5+Mom): {dm_factors_monthly.shape}")
print(f"  Columns: {list(dm_factors_monthly.columns)}")

# ---- Emerging Markets Factors (Monthly) — for AVEM ----
print("\n" + "-"*70)
print("Loading Emerging Markets Fama-French 5 Factors (Monthly)...")
em_ff5_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/Emerging_5_Factors_CSV.zip'
em_ff5 = load_ff_csv_from_zip(em_ff5_url, freq='monthly')
em_ff5.columns = em_ff5.columns.str.strip()
em_ff5 = em_ff5.rename(columns={'Mkt-RF': 'Mkt_RF'})

print(f"  Shape: {em_ff5.shape}")
print(f"  Date range: {em_ff5.index[0].date()} to {em_ff5.index[-1].date()}")
print(f"  Columns: {list(em_ff5.columns)}")

# Note: EM Momentum is available separately but we use FF5 only for AVEM
# to keep the model consistent (the EM WML series can have gaps)
em_factors = em_ff5.copy()

print(f"\nSample DM monthly factors (last 5 rows):")
print(dm_factors_monthly.tail())

# ============================================================================
# Step 2: Download ETF Prices from Xetra + EUR/USD FX Rate
# ============================================================================

etf_tickers = {
    'AVWS.DE': 'Global Small Cap Value',
    'AVWC.DE': 'Global Equity',
    'AVEM.DE': 'Emerging Markets Equity',
    'IS3R.DE': 'MSCI World Momentum (IWMO)',
}

print("="*70)
print("DOWNLOADING ETF DATA (XETRA)")
print("="*70)

# Download ETF prices
etf_prices = {}
for ticker, name in etf_tickers.items():
    data = yf.download(ticker, period='max', progress=False)
    # Handle multi-level columns from yfinance
    if isinstance(data.columns, pd.MultiIndex):
        data.columns = data.columns.get_level_values(0)
    etf_prices[ticker] = data['Close']
    print(f"\n{ticker} ({name}):")
    print(f"  Observations: {len(data)}")
    print(f"  Date range: {data.index[0].date()} to {data.index[-1].date()}")
    print(f"  First price: {data['Close'].iloc[0]:.2f} EUR")
    print(f"  Last price: {data['Close'].iloc[-1]:.2f} EUR")

# ---- Data cleaning: IS3R.DE (IWMO) ----
# Yahoo Finance reports erratic prices for IS3R.DE in its early trading period
# (Oct 2014 – Jun 2017): isolated single-day price spikes of ±20% that are
# clearly not real market moves for an MSCI World Momentum ETF. These appear
# to be stale/incorrect prices from periods of very low liquidity on Xetra.
#
# Two-step cleaning:
#  1. Drop the first calendar month (Oct 2014). The ETF's initial trading days
#     produce a spurious −19.5% monthly return that is absent from the
#     underlying MSCI World Momentum index.
#  2. Remove any remaining daily observation with |price change| > 10%.

raw_prices = etf_prices['IS3R.DE']

# Step 1: trim first calendar month (launch-period artefact)
first_full_month = raw_prices.index[0] + pd.offsets.MonthBegin(1)
raw_prices = raw_prices[raw_prices.index >= first_full_month]
print(f"\n  IS3R.DE DATA CLEANING:")
print(f"    Trimmed first month (before {first_full_month.date()}) — launch-period artefact")

# Step 2: remove remaining daily spikes > 10%
daily_pct = raw_prices.pct_change().abs()
bad_days = daily_pct[daily_pct > 0.10].index
if len(bad_days) > 0:
    bad_and_next = set(bad_days)
    for d in bad_days:
        next_idx = raw_prices.index.get_indexer([d], method='pad')[0] + 1
        if next_idx < len(raw_prices.index):
            bad_and_next.add(raw_prices.index[next_idx])
    raw_prices = raw_prices.drop(index=list(bad_and_next), errors='ignore')
    print(f"    Removed {len(bad_and_next)} days with |daily Δ| > 10%")
    print(f"    Spike dates span: {min(bad_days).date()} to {max(bad_days).date()}")
else:
    print("    No remaining daily spikes > 10%")
print(f"    Clean observations: {len(raw_prices)} (trimmed from {len(etf_prices['IS3R.DE'])})")
etf_prices['IS3R.DE'] = raw_prices

# Download EUR/USD exchange rate
print("\n" + "-"*70)
print("Downloading EUR/USD exchange rate...")
fx_data = yf.download('EURUSD=X', start='2014-09-01', progress=False)
if isinstance(fx_data.columns, pd.MultiIndex):
    fx_data.columns = fx_data.columns.get_level_values(0)
eurusd = fx_data['Close']
print(f"  Observations: {len(eurusd)}")
print(f"  Date range: {eurusd.index[0].date()} to {eurusd.index[-1].date()}")

# Compute daily returns in EUR, then convert to USD
print("\n" + "="*70)
print("COMPUTING DAILY RETURNS (USD-CONVERTED)")
print("="*70)

etf_returns_usd = {}
for ticker, name in etf_tickers.items():
    prices_eur = etf_prices[ticker]
    
    # Daily EUR returns
    ret_eur = prices_eur.pct_change().dropna()
    
    # EUR/USD rate changes
    fx_aligned = eurusd.reindex(ret_eur.index, method='ffill')
    fx_ret = fx_aligned.pct_change()
    
    # Convert to USD: R_usd = (1 + R_eur)(1 + R_fx) - 1
    ret_usd = (1 + ret_eur) * (1 + fx_ret) - 1
    ret_usd = ret_usd.dropna()
    
    # Convert to percentage points to match FF factor data
    etf_returns_usd[ticker] = ret_usd * 100
    
    print(f"\n{ticker} ({name}):")
    print(f"  Daily return obs: {len(ret_usd)}")
    print(f"  Mean daily return: {ret_usd.mean()*100:.4f}%")
    print(f"  Daily volatility: {ret_usd.std()*100:.4f}%")
    print(f"  Ann. return (approx): {ret_usd.mean()*252*100:.2f}%")
    print(f"  Ann. volatility (approx): {ret_usd.std()*np.sqrt(252)*100:.2f}%")

# ============================================================================
# Step 3: Merge ETF Returns with Factor Data
# ============================================================================

print("="*70)
print("MERGING ETF RETURNS WITH FAMA-FRENCH FACTORS")
print("="*70)

# ---- A. Daily: merge all ETFs with DM daily factors (descriptive stats & cumulative charts) ----
etf_dataframes = {}
for ticker, name in etf_tickers.items():
    ret = etf_returns_usd[ticker]
    merged = pd.DataFrame({'ETF_Return': ret})
    merged = merged.join(factors, how='inner')
    merged['Excess_Return'] = merged['ETF_Return'] - merged['RF']
    etf_dataframes[ticker] = merged
    
    print(f"\n{ticker} ({name}) — DM daily (descriptive stats only):")
    print(f"  Matched observations: {len(merged)}")
    if len(merged) > 0:
        print(f"  Date range: {merged.index[0].date()} to {merged.index[-1].date()}")
        print(f"  Mean excess return: {merged['Excess_Return'].mean():.4f}%/day")

# ---- B. Monthly: build monthly excess returns for ALL ETFs (for regressions) ----
print("\n" + "="*70)
print("MONTHLY RETURN DATA (FOR REGRESSIONS)")
print("="*70)

etf_monthly_dfs = {}

def build_monthly_df(ticker, monthly_factors, factor_label):
    """Aggregate daily→monthly and merge with the appropriate monthly factor set."""
    daily_ret = etf_returns_usd[ticker] / 100          # decimal
    monthly_ret = (1 + daily_ret).resample('ME').prod() - 1
    monthly_ret = monthly_ret * 100                     # back to % to match FF
    
    df = pd.DataFrame({'ETF_Return': monthly_ret})
    df.index = df.index.to_period('M').to_timestamp()   # align to month-start
    df = df.join(monthly_factors, how='inner')
    df['Excess_Return'] = df['ETF_Return'] - df['RF']
    
    print(f"\n  {ticker} — {factor_label}:")
    print(f"    Monthly obs: {len(df)}")
    if len(df) > 0:
        print(f"    Date range: {df.index[0].date()} to {df.index[-1].date()}")
        print(f"    Mean monthly excess: {df['Excess_Return'].mean():.2f}%")
    return df

# AVWS, AVWC, and IWMO → DM monthly FF5+Mom
for ticker in ['AVWS.DE', 'AVWC.DE', 'IS3R.DE']:
    etf_monthly_dfs[ticker] = build_monthly_df(ticker, dm_factors_monthly,
                                                'DM monthly FF5+Mom')

# AVEM → EM monthly FF5
etf_monthly_dfs['AVEM.DE'] = build_monthly_df('AVEM.DE', em_factors,
                                               'EM monthly FF5')

# Backward-compatibility alias
avem_monthly_df = etf_monthly_dfs['AVEM.DE']

# ---- Summary table ----
print("\n" + "="*70)
print("DATA SUMMARY")
print("="*70)
print(f"{'ETF':<12} {'Daily':>6} {'Monthly':>8} {'Start (m)':>12} {'End (m)':>12} {'Ann.Ret(m)':>11} {'Factors'}")
print("-"*85)
for ticker, name in etf_tickers.items():
    df_d = etf_dataframes[ticker]
    df_m = etf_monthly_dfs[ticker]
    ann_ret = df_m['Excess_Return'].mean() * 12
    flabel = 'DM monthly' if ticker != 'AVEM.DE' else 'EM monthly'
    print(f"{ticker:<12} {len(df_d):>6} {len(df_m):>8} "
          f"{df_m.index[0].strftime('%Y-%m-%d'):>12} "
          f"{df_m.index[-1].strftime('%Y-%m-%d'):>12} "
          f"{ann_ret:>10.2f}%  {flabel}")

# ============================================================================
# Descriptive Statistics and Return Distributions
# ============================================================================

print("="*70)
print("DESCRIPTIVE STATISTICS — DAILY EXCESS RETURNS (% per day)")
print("="*70)

desc_stats = []
for ticker, name in etf_tickers.items():
    r = etf_dataframes[ticker]['Excess_Return']
    jb_stat, jb_p, skew_val, kurt_val = jarque_bera(r)
    
    # Ljung-Box test for autocorrelation (lag 5)
    lb_result = acorr_ljungbox(r, lags=[5], return_df=True)
    lb_p = lb_result['lb_pvalue'].values[0]
    
    # Ljung-Box on squared returns (volatility clustering)
    lb_sq = acorr_ljungbox(r**2, lags=[5], return_df=True)
    lb_sq_p = lb_sq['lb_pvalue'].values[0]
    
    short = ticker.replace('.DE', '')
    # Use original ticker name for IS3R
    display_name = 'IWMO' if ticker == 'IS3R.DE' else short
    
    desc_stats.append({
        'ETF': display_name,
        'N': len(r),
        'Mean (%)': f"{r.mean():.4f}",
        'Std (%)': f"{r.std():.4f}",
        'Skewness': f"{skew_val:.3f}",
        'Kurtosis': f"{kurt_val:.3f}",
        'JB p-value': f"{jb_p:.4f}",
        'LB(5) p': f"{lb_p:.4f}",
        'LB²(5) p': f"{lb_sq_p:.4f}",
        'Min (%)': f"{r.min():.2f}",
        'Max (%)': f"{r.max():.2f}",
    })

desc_df = pd.DataFrame(desc_stats).set_index('ETF')
print(desc_df.to_string())

print("""
INTERPRETATION:
━━━━━━━━━━━━━━
• JB p-value < 0.05 → Returns are NOT normally distributed
  
• LB(5) p-value < 0.05 → significant autocorrelation detected

• LB²(5) p-value < 0.05 → ARCH effects present

→ These diagnostics confirm that Newey-West standard errors are appropriate
  for our factor regressions (they account for both autocorrelation and 
  heteroskedasticity).
""")

# ============================================================================
# Visualisation: Return Distributions and Cumulative Returns
# ============================================================================

n_etfs = len(etf_tickers)
fig, axes = plt.subplots(2, n_etfs, figsize=(5 * n_etfs, 10))

for i, (ticker, name) in enumerate(etf_tickers.items()):
    r = etf_dataframes[ticker]['Excess_Return']
    short_name = 'IWMO' if ticker == 'IS3R.DE' else ticker.replace('.DE', '')
    
    # Top row: Return histograms with normal overlay
    ax = axes[0, i]
    ax.hist(r, bins=50, density=True, alpha=0.7, color='steelblue', 
            edgecolor='white', linewidth=0.5)
    x = np.linspace(r.min(), r.max(), 100)
    ax.plot(x, stats.norm.pdf(x, r.mean(), r.std()), 'r-', linewidth=2,
            label='Normal')
    ax.set_title(f'{short_name}: Daily Excess Returns', fontweight='bold')
    ax.set_xlabel('Return (%)')
    ax.set_ylabel('Density')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    
    # Bottom row: Cumulative returns
    ax = axes[1, i]
    cum_ret = (1 + r / 100).cumprod()
    ax.plot(cum_ret.index, cum_ret, color='darkblue', linewidth=1.5)
    ax.axhline(1, color='red', linestyle='--', linewidth=0.8, alpha=0.5)
    ax.set_title(f'{short_name}: Cumulative Growth of $1', fontweight='bold')
    ax.set_xlabel('Date')
    ax.set_ylabel('Value ($)')
    ax.grid(True, alpha=0.3)
    ax.tick_params(axis='x', rotation=30)

plt.suptitle('ETF Return Distributions and Cumulative Performance',
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('etf_descriptive.png', dpi=100, bbox_inches='tight')
plt.show()

print("""
TOP ROW: Return distributions. Compare the histogram to the red normal curve.
The histograms are LEPTOKURTIC. The most visually striking feature is the SHARP, 
TALL PEAK around zero — the histogram towers above the red normal curve at the 
centre, showing that returns cluster more tightly around zero on most days than a 
normal distribution would predict.

The other side of leptokurtosis — FATTER TAILS — is harder to spot visually in a 
histogram because extreme returns are rare by definition (just a few sparse bars 
far from zero). However, the fat tails are confirmed by:
  • The kurtosis statistics in the table above (excess kurtosis > 3)
  • The Jarque-Bera test rejecting normality (JB p-value < 0.05)
  • The QQ plots in the residual diagnostics (Section 4), where deviations from 
    the straight line at the extremes make fat tails clearly visible

Together, "peaked centre + fat tails" is the hallmark of excess kurtosis in 
financial returns. It means that a normal distribution simultaneously overestimates 
the frequency of modest moves and underestimates the frequency of extreme moves. 
This is exactly why we use Newey-West standard errors throughout: standard OLS 
errors assume normality, which would understate the true uncertainty.

BOTTOM ROW: Cumulative returns since inception. The Avantis ETFs have a very short 
period (~17 months), while IWMO has ~11 years of history. These cumulative paths 
reflect market movements and noise; the short-sample ETFs do not yet provide 
reliable estimates of long-run factor premia.
""")

# ============================================================================
# Monthly Descriptive Statistics and Cross-ETF Correlations
# ============================================================================
# Since all regressions use monthly data, we show the monthly-level statistics
# and the correlation structure that drives portfolio diversification.

print("="*70)
print("MONTHLY EXCESS RETURN STATISTICS (% per month)")
print("="*70)

stats_table = []
for ticker, name in etf_tickers.items():
    r = etf_monthly_dfs[ticker]['Excess_Return']
    short = 'IWMO' if ticker == 'IS3R.DE' else ticker.replace('.DE', '')
    ann_factor = 12
    stats_table.append({
        'ETF': short,
        'N': len(r),
        'Mean (%)': f"{r.mean():.2f}",
        'Std (%)': f"{r.std():.2f}",
        'Ann. Ret (%)': f"{r.mean()*ann_factor:.1f}",
        'Ann. Vol (%)': f"{r.std()*np.sqrt(ann_factor):.1f}",
        'Sharpe': f"{r.mean()/r.std()*np.sqrt(ann_factor):.2f}" if r.std() > 0 else 'n/a',
        'Min (%)': f"{r.min():.2f}",
        'Max (%)': f"{r.max():.2f}",
    })
print(pd.DataFrame(stats_table).set_index('ETF').to_string())

# ---- Cross-ETF Correlation Matrix (monthly excess returns, common sample) ----
print("\n" + "="*70)
print("CROSS-ETF CORRELATION MATRIX (monthly excess returns, common sample)")
print("="*70)

etf_list = list(etf_tickers.keys())
etf_short = ['IWMO' if t == 'IS3R.DE' else t.replace('.DE', '') for t in etf_list]
etf_order = ['AVWS', 'AVWC', 'AVEM', 'IWMO']

# Find common monthly dates across all ETFs
common_months = etf_monthly_dfs[etf_list[0]].index
for t in etf_list[1:]:
    common_months = common_months.intersection(etf_monthly_dfs[t].index)

corr_df = pd.DataFrame({
    short: etf_monthly_dfs[t].loc[common_months, 'Excess_Return']
    for t, short in zip(etf_list, etf_short)
})[etf_order]

corr = corr_df.corr()
print(f"\nCommon sample: {len(common_months)} months "
      f"({common_months[0].strftime('%Y-%m')} to {common_months[-1].strftime('%Y-%m')})\n")
print(corr.round(3).to_string())

# Heatmap
fig, ax = plt.subplots(figsize=(6, 5))
sns.heatmap(corr, annot=True, fmt='.2f', cmap='RdYlBu_r', center=0,
            vmin=-1, vmax=1, square=True, linewidths=0.5, ax=ax,
            xticklabels=etf_order, yticklabels=etf_order)
ax.set_title('Cross-ETF Correlation (Monthly Excess Returns, Common Sample)',
             fontweight='bold', fontsize=11)
plt.tight_layout()
plt.show()

print("""
These correlations drive portfolio diversification. High correlations between
AVWS and AVWC are expected — both hold developed-market equities. IWMO's
distinct correlation pattern (it loads on momentum, not value) makes it a
useful diversifier, which motivates the portfolio extension in Section 8.
""")

# ============================================================================
# Factor Regression Helper Function
# ============================================================================

def run_factor_regression(etf_data, factor_cols, ticker_name, nw_lags=None):
    """
    Run OLS regression with Newey-West standard errors.
    
    Parameters
    ----------
    etf_data : DataFrame with 'Excess_Return' and factor columns
    factor_cols : list of factor column names
    ticker_name : string for display
    nw_lags : int or None (auto = int(4*(T/100)^(2/9)))
    
    Returns
    -------
    statsmodels RegressionResults
    """
    y = etf_data['Excess_Return'].values
    X = sm.add_constant(etf_data[factor_cols].values)
    
    # Automatic Newey-West lag selection (Andrews, 1991 rule of thumb)
    T = len(y)
    if nw_lags is None:
        nw_lags = int(np.ceil(4 * (T / 100) ** (2/9)))
    
    model = sm.OLS(y, X)
    results = model.fit(cov_type='HAC', cov_kwds={'maxlags': nw_lags})
    
    return results, nw_lags

def print_regression_results(results, factor_cols, ticker_name, nw_lags, 
                              freq='daily'):
    """Pretty-print regression results.
    
    Parameters
    ----------
    freq : 'daily' or 'monthly' — controls how alpha is annualised
    """
    param_names = ['Alpha (α)'] + factor_cols
    ann_mult = 252 if freq == 'daily' else 12
    freq_label = 'daily' if freq == 'daily' else 'monthly'
    
    print(f"\n{'='*65}")
    print(f"  {ticker_name}")
    print(f"  Newey-West lags: {nw_lags} | N = {int(results.nobs)} | R² = {results.rsquared:.4f}")
    print(f"{'='*65}")
    print(f"{'Factor':<14} {'Coeff':>10} {'Std Err':>10} {'t-stat':>10} {'p-value':>10} {'Sig':>5}")
    print(f"{'-'*65}")
    
    for i, name in enumerate(param_names):
        coef = results.params[i]
        se = results.bse[i]
        t = results.tvalues[i]
        p = results.pvalues[i]
        sig = '***' if p < 0.01 else '**' if p < 0.05 else '*' if p < 0.10 else ''
        print(f"{name:<14} {coef:>10.4f} {se:>10.4f} {t:>10.3f} {p:>10.4f} {sig:>5}")
    
    # Annualise alpha
    alpha = results.params[0]
    alpha_annual = alpha * ann_mult
    print(f"\n  Annualised alpha ≈ {alpha_annual:.2f}% (= {alpha:.4f}% × {ann_mult})")
    print(f"  Residual std ({freq_label}) = {np.sqrt(results.mse_resid):.4f}%")

print("Regression helper functions defined.")

# ============================================================================
# Run FF3, Full Factor Model, and 5-Factor (no CMA) Regressions — ALL Monthly
# ============================================================================

ff3_cols    = ['Mkt_RF', 'SMB', 'HML']
ff5_cols    = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA']
ff6_cols    = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD']
noCMA_cols  = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'UMD']  # 5-factor without CMA

# Define which full model each ETF uses
etf_model_spec = {
    'AVWS.DE': {'cols_full': ff6_cols, 'label': 'DM Monthly FF5+Mom', 'nw': None},
    'AVWC.DE': {'cols_full': ff6_cols, 'label': 'DM Monthly FF5+Mom', 'nw': None},
    'AVEM.DE': {'cols_full': ff5_cols, 'label': 'EM Monthly FF5',     'nw': 1},
    'IS3R.DE': {'cols_full': ff6_cols, 'label': 'DM Monthly FF5+Mom', 'nw': None},
}

results_3f = {}
results_6f = {}       # full model (FF6 for DM, FF5 for EM)
results_noCMA = {}    # 5-factor without CMA (DM ETFs only)

# ---- FF3 for all ETFs (monthly) ----
print("█"*70)
print("  FAMA-FRENCH 3-FACTOR REGRESSIONS — ALL MONTHLY")
print("█"*70)

for ticker, name in etf_tickers.items():
    spec = etf_model_spec[ticker]
    df = etf_monthly_dfs[ticker]
    nw = spec['nw']
    res, lags = run_factor_regression(df, ff3_cols, f"{ticker} — {name}", nw_lags=nw)
    results_3f[ticker] = res
    print_regression_results(res, ff3_cols,
                             f"{ticker} — {name} [3-Factor, {spec['label']}]",
                             lags, freq='monthly')

# ---- Full model for all ETFs (monthly) ----
print("\n\n" + "█"*70)
print("  FULL FACTOR MODEL — ALL MONTHLY")
print("  AVWS/AVWC/IWMO: FF5+Mom (6 factors)  |  AVEM: FF5 (5 factors)")
print("█"*70)

for ticker, name in etf_tickers.items():
    spec = etf_model_spec[ticker]
    df = etf_monthly_dfs[ticker]
    nw = spec['nw']
    full_cols = spec['cols_full']
    res, lags = run_factor_regression(df, full_cols, f"{ticker} — {name}", nw_lags=nw)
    results_6f[ticker] = res
    print_regression_results(res, full_cols,
                             f"{ticker} — {name} [{spec['label']}]",
                             lags, freq='monthly')

# ---- 5-Factor without CMA (MKT, SMB, HML, RMW, UMD) for DM ETFs ----
print("\n\n" + "█"*70)
print("  5-FACTOR WITHOUT CMA (MKT, SMB, HML, RMW, UMD)")
print("  AVWS, AVWC, and IWMO only (DM monthly)")
print("█"*70)
print("\n  CMA is dropped to check whether the remaining coefficients are")
print("  sensitive to the HML–CMA collinearity typical of short samples.\n")

for ticker, name in etf_tickers.items():
    if ticker == 'AVEM.DE':
        continue  # EM already uses FF5 without UMD
    df = etf_monthly_dfs[ticker]
    res, lags = run_factor_regression(df, noCMA_cols, f"{ticker} — {name}")
    results_noCMA[ticker] = res
    print_regression_results(res, noCMA_cols,
                             f"{ticker} — {name} [5-Factor no CMA, DM Monthly]",
                             lags, freq='monthly')

# Store references for backward compatibility
avem_em_3f = results_3f['AVEM.DE']
avem_em_6f = results_6f['AVEM.DE']
avem_em_lags_3f = 1
avem_em_lags_6f = 1

# ============================================================================
# Residual Diagnostics — Validating Our Regression Assumptions (All Monthly)
# ============================================================================
# Following the diagnostic approach from Notebook 3, Section 5

print("="*70)
print("RESIDUAL DIAGNOSTICS (FULL FACTOR MODEL — MONTHLY)")
print("="*70)
print("All regressions use monthly returns. AVWS/AVWC/IWMO: DM FF5+Mom; AVEM: EM FF5.\n")

n_etfs = len(etf_tickers)
fig, axes = plt.subplots(n_etfs, 3, figsize=(18, 4 * n_etfs + 1))

for i, (ticker, name) in enumerate(etf_tickers.items()):
    res = results_6f[ticker]
    resid = res.resid
    short = 'IWMO' if ticker == 'IS3R.DE' else ticker.replace('.DE', '')
    dates = etf_monthly_dfs[ticker].index
    spec = etf_model_spec[ticker]
    max_acf_lags = min(8, len(resid) // 2 - 1)  # fewer lags for short monthly samples
    
    # Column 1: Residual time series (bar chart for monthly)
    ax = axes[i, 0]
    colours = ['steelblue' if r >= 0 else 'salmon' for r in resid]
    ax.bar(dates, resid, width=25, color=colours, alpha=0.8, edgecolor='none')
    ax.axhline(0, color='red', linestyle='--', linewidth=0.8)
    ax.set_title(f'{short}: Monthly Residuals ({spec["label"]})', fontweight='bold', fontsize=10)
    ax.set_ylabel('Residual (%)')
    ax.grid(True, alpha=0.3)
    ax.tick_params(axis='x', rotation=30)
    
    # Column 2: Residual ACF
    ax = axes[i, 1]
    from statsmodels.graphics.tsaplots import plot_acf
    plot_acf(resid, ax=ax, lags=max_acf_lags, alpha=0.05)
    ax.set_title(f'{short}: Residual ACF', fontweight='bold', fontsize=10)
    ax.grid(True, alpha=0.3)
    
    # Column 3: QQ Plot
    ax = axes[i, 2]
    stats.probplot(resid, dist='norm', plot=ax)
    ax.set_title(f'{short}: QQ Plot', fontweight='bold', fontsize=10)
    ax.grid(True, alpha=0.3)

plt.suptitle('Full Factor Model — Residual Diagnostics (All Monthly)',
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('etf_diagnostics.png', dpi=100, bbox_inches='tight')
plt.show()

# Formal tests for all ETFs
print("\nFORMAL DIAGNOSTIC TESTS (monthly residuals):")
print(f"{'ETF':<10} {'N':>4} {'Ljung-Box(4)':>14} {'p-value':>9} {'JB stat':>10} {'JB p':>9} {'Skew':>7} {'Kurt':>7}")
print("-"*80)

for ticker, name in etf_tickers.items():
    res = results_6f[ticker]
    resid = res.resid
    short = 'IWMO' if ticker == 'IS3R.DE' else ticker.replace('.DE', '')
    T = len(resid)
    max_lb = min(4, T // 3)
    
    # Ljung-Box
    lb_result = sm.stats.acorr_ljungbox(resid, lags=[max_lb], return_df=True)
    lb, lb_p = lb_result.iloc[0]['lb_stat'], lb_result.iloc[0]['lb_pvalue']
    
    # Jarque-Bera
    jb_stat, jb_p = stats.jarque_bera(resid)
    skew_val = stats.skew(resid)
    kurt_val = stats.kurtosis(resid)
    
    print(f"{short:<10} {T:>4} {lb:>14.2f} {lb_p:>9.4f} {jb_stat:>10.2f} {jb_p:>9.4f} {skew_val:>7.2f} {kurt_val:>7.2f}")

print("\nInterpretation:")
print("  • Ljung-Box tests for residual autocorrelation (H₀: no serial correlation)")
print("  • Jarque-Bera tests for normality (H₀: residuals are normally distributed)")
print("  • With ~15 monthly observations (Avantis ETFs), these tests have low power;")
print("    interpret cautiously. IWMO has ~130 months, giving much more reliable diagnostics.")
print("  • Newey-West SEs remain valid regardless of the JB test outcome")

# ============================================================================
# Comparative Factor Exposure Chart (All Monthly)
# ============================================================================

factor_names = ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD']
all_tickers = list(etf_tickers.keys())
etf_short_names = ['IWMO' if t == 'IS3R.DE' else t.replace('.DE', '') for t in all_tickers]
colors = ['#2196F3', '#4CAF50', '#FF9800', '#9C27B0']  # Blue, Green, Orange, Purple

# Collect coefficients from the full model (FF6 for DM, FF5 for EM)
coef_matrix = np.zeros((len(all_tickers), len(factor_names)))
sig_matrix = np.ones((len(all_tickers), len(factor_names)))
se_matrix = np.zeros((len(all_tickers), len(factor_names)))

for i, ticker in enumerate(all_tickers):
    res = results_6f[ticker]
    spec = etf_model_spec[ticker]
    model_factors = spec['cols_full']
    for j, factor in enumerate(factor_names):
        if factor in model_factors:
            k = model_factors.index(factor) + 1
            coef_matrix[i, j] = res.params[k]
            sig_matrix[i, j] = res.pvalues[k]
            se_matrix[i, j] = res.bse[k]

# Plot: Side-by-side bar chart of factor loadings + alpha comparison
fig, axes = plt.subplots(1, 2, figsize=(20, 7))

# Panel A: Factor Betas (full model)
ax = axes[0]
x = np.arange(len(factor_names))
n_bars = len(all_tickers)
width = 0.8 / n_bars
for i, (etf_name, color) in enumerate(zip(etf_short_names, colors)):
    ticker = all_tickers[i]
    freq_label = 'DM monthly' if ticker != 'AVEM.DE' else 'EM monthly'
    offset = (i - (n_bars - 1) / 2) * width
    bars = ax.bar(x + offset, coef_matrix[i], width, 
                  label=f'{etf_name} ({freq_label})', color=color, alpha=0.85, 
                  edgecolor='white')
    for j, bar in enumerate(bars):
        ypos = bar.get_height() + 0.01 if bar.get_height() >= 0 else bar.get_height() - 0.03
        if sig_matrix[i, j] < 0.01:
            ax.text(bar.get_x() + bar.get_width()/2, ypos, '***', 
                    ha='center', va='bottom' if bar.get_height() >= 0 else 'top',
                    fontsize=7, fontweight='bold')
        elif sig_matrix[i, j] < 0.05:
            ax.text(bar.get_x() + bar.get_width()/2, ypos, '**',
                    ha='center', va='bottom' if bar.get_height() >= 0 else 'top',
                    fontsize=7, fontweight='bold')
        elif sig_matrix[i, j] < 0.10:
            ax.text(bar.get_x() + bar.get_width()/2, ypos, '*',
                    ha='center', va='bottom' if bar.get_height() >= 0 else 'top',
                    fontsize=7)

ax.set_xticks(x)
ax.set_xticklabels(factor_names, fontsize=11)
ax.set_ylabel('Beta Coefficient', fontsize=12)
ax.set_title('Panel A: Factor Betas — Full Model (Monthly)', fontweight='bold', fontsize=13)
ax.axhline(0, color='black', linewidth=0.8)
ax.legend(loc='upper right', fontsize=9)
ax.grid(True, alpha=0.2, axis='y')

# Panel B: Alpha comparison
ax = axes[1]
alphas = []
alpha_ses = []
alpha_sigs = []
labels = []
r2_vals = []
bar_colours = []

for i, ticker in enumerate(all_tickers):
    res = results_6f[ticker]
    short = etf_short_names[i]
    spec = etf_model_spec[ticker]
    ann = 12
    a = res.params[0] * ann
    a_se = res.bse[0] * ann
    a_sig = res.pvalues[0]
    r2 = res.rsquared
    freq_label = 'DM monthly' if ticker != 'AVEM.DE' else 'EM monthly'
    label = f'{short}\n({freq_label})'
    alphas.append(a)
    alpha_ses.append(a_se)
    alpha_sigs.append(a_sig)
    labels.append(label)
    r2_vals.append(r2)
    bar_colours.append(colors[i])

x_alpha = np.arange(len(alphas))
bars = ax.bar(x_alpha, alphas, color=bar_colours, alpha=0.85, edgecolor='white', width=0.6)
ax.errorbar(x_alpha, alphas, yerr=[1.96*s for s in alpha_ses], 
            fmt='none', color='black', capsize=5, linewidth=1.5)

for j, bar in enumerate(bars):
    ypos = bar.get_height() + 1.96*alpha_ses[j] + 0.3 if bar.get_height() >= 0 \
           else bar.get_height() - 1.96*alpha_ses[j] - 0.3
    ax.text(bar.get_x() + bar.get_width()/2, ypos, f'R²={r2_vals[j]:.3f}',
            ha='center', fontsize=9, style='italic')
    sig_ypos = bar.get_height() + 0.1 if bar.get_height() >= 0 else bar.get_height() - 0.1
    if alpha_sigs[j] < 0.01:
        ax.text(bar.get_x() + bar.get_width()/2, sig_ypos, '***',
                ha='center', fontsize=10, fontweight='bold')
    elif alpha_sigs[j] < 0.05:
        ax.text(bar.get_x() + bar.get_width()/2, sig_ypos, '**',
                ha='center', fontsize=10, fontweight='bold')

ax.set_xticks(x_alpha)
ax.set_xticklabels(labels, fontsize=10)
ax.set_ylabel('Annualised Alpha (%)', fontsize=12)
ax.set_title('Panel B: Alpha Estimates (x12 annualised, ±1.96 SE)', fontweight='bold', fontsize=13)
ax.axhline(0, color='black', linewidth=0.8)
ax.grid(True, alpha=0.2, axis='y')

plt.suptitle('Comparative Factor Analysis (Monthly)',
             fontsize=15, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('etf_factor_comparison.png', dpi=100, bbox_inches='tight')
plt.show()

# ============================================================================
# Summary Table
# ============================================================================
print("\n" + "="*100)
print("FULL MODEL SUMMARY (FF5+Mom for AVWS/AVWC/IWMO, FF5 for AVEM)")
print("="*100)
print(f"{'ETF':<8} {'Model':<18} {'N':>4} {'R²':>7} {'α (ann)':>9} {'MKT':>7} {'SMB':>7} {'HML':>7} {'RMW':>7} {'CMA':>7} {'UMD':>7}")
print("-"*100)
for i, ticker in enumerate(all_tickers):
    res = results_6f[ticker]
    short = etf_short_names[i]
    spec = etf_model_spec[ticker]
    alpha_ann = res.params[0] * 12
    betas = []
    for factor in factor_names:
        if factor in spec['cols_full']:
            k = spec['cols_full'].index(factor) + 1
            coef = res.params[k]
            p = res.pvalues[k]
            sig = '***' if p < 0.01 else '**' if p < 0.05 else '*' if p < 0.10 else ''
            betas.append(f'{coef:>5.2f}{sig}')
        else:
            betas.append(f'{"n/a":>7}')
    print(f"{short:<8} {spec['label']:<18} {int(res.nobs):>4} {res.rsquared:>7.4f} {alpha_ann:>8.2f}% " 
          + ' '.join(f'{b:>7}' for b in betas))

# ============================================================================
# 6.1 — Build the Base Portfolio (50% AVWS / 40% AVWC / 10% AVEM)
# ============================================================================

# Fixed weights — no optimisation needed
base_weights = {
    'AVWS.DE': 0.50,
    'AVWC.DE': 0.40,
    'AVEM.DE': 0.10,
}
base_tickers = list(base_weights.keys())
base_labels  = [t.replace('.DE', '') for t in base_tickers]

# ---- A. Find common daily/monthly dates across ALL 4 ETFs ----
# (We need IWMO aligned too for comparison charts and Section 8)
common_dates = etf_dataframes['AVWS.DE'].index
for t in ['AVWC.DE', 'AVEM.DE', 'IS3R.DE']:
    common_dates = common_dates.intersection(etf_dataframes[t].index)

common_months = etf_monthly_dfs['AVWS.DE'].index
for t in ['AVWC.DE', 'AVEM.DE', 'IS3R.DE']:
    common_months = common_months.intersection(etf_monthly_dfs[t].index)

print(f"Common daily dates (all 4 ETFs): {len(common_dates)}")
print(f"Date range: {common_dates[0].date()} to {common_dates[-1].date()}")
print(f"\nCommon monthly dates: {len(common_months)}")
print(f"Date range: {common_months[0].date()} to {common_months[-1].date()}")

# ---- B. Build daily return matrix ----
port_data = pd.DataFrame(index=common_dates)
for t in base_tickers:
    short = t.replace('.DE', '')
    port_data[f'{short}_excess'] = etf_dataframes[t].loc[common_dates, 'Excess_Return']
    port_data[f'{short}_ret']    = etf_dataframes[t].loc[common_dates, 'ETF_Return']

# Include IWMO for comparison charts and Section 8
port_data['IWMO_excess'] = etf_dataframes['IS3R.DE'].loc[common_dates, 'Excess_Return']
port_data['IWMO_ret']    = etf_dataframes['IS3R.DE'].loc[common_dates, 'ETF_Return']

# Add DM factor columns
for col in ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD', 'RF']:
    port_data[col] = etf_dataframes['AVWS.DE'].loc[common_dates, col]

# ---- C. Build monthly return matrix ----
port_monthly = pd.DataFrame(index=common_months)
for t in base_tickers:
    short = t.replace('.DE', '')
    port_monthly[f'{short}_excess'] = etf_monthly_dfs[t].loc[common_months, 'Excess_Return']
    port_monthly[f'{short}_ret']    = etf_monthly_dfs[t].loc[common_months, 'ETF_Return']

port_monthly['IWMO_excess'] = etf_monthly_dfs['IS3R.DE'].loc[common_months, 'Excess_Return']
port_monthly['IWMO_ret']    = etf_monthly_dfs['IS3R.DE'].loc[common_months, 'ETF_Return']

for col in ['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA', 'UMD', 'RF']:
    port_monthly[col] = dm_factors_monthly.loc[common_months, col]

# ---- D. Compute base portfolio returns ----
port_data['Port_ret']        = sum(base_weights[t] * port_data[f'{t.replace(".DE","")}_ret']        for t in base_tickers)
port_data['Port_excess']     = sum(base_weights[t] * port_data[f'{t.replace(".DE","")}_excess']     for t in base_tickers)
port_monthly['Port_ret']     = sum(base_weights[t] * port_monthly[f'{t.replace(".DE","")}_ret']     for t in base_tickers)
port_monthly['Port_excess']  = sum(base_weights[t] * port_monthly[f'{t.replace(".DE","")}_excess']  for t in base_tickers)

# ---- E. Summary Statistics ----
print("\n" + "=" * 70)
print("BASE PORTFOLIO: 50% AVWS / 40% AVWC / 10% AVEM")
print("=" * 70)

port_ret = port_data['Port_ret']
ann = 252

summary = {
    'Mean daily return (%)':       port_ret.mean(),
    'Std. dev. daily (%)':         port_ret.std(),
    'Annualised return (%)':       port_ret.mean() * ann,
    'Annualised volatility (%)':   port_ret.std() * np.sqrt(ann),
    'Sharpe ratio (annualised)':   port_ret.mean() / port_ret.std() * np.sqrt(ann),
    'Skewness':                    stats.skew(port_ret),
    'Excess kurtosis':             stats.kurtosis(port_ret),
    'Max daily loss (%)':          port_ret.min(),
    'Max daily gain (%)':          port_ret.max(),
}
for k, v in summary.items():
    print(f"  {k:<32s} {v:>10.4f}")

print(f"\n  Monthly obs (for regression): {len(port_monthly)}")
print(f"  Mean monthly excess return:   {port_monthly['Port_excess'].mean():.2f}%")
print(f"  Ann. excess return (×12):     {port_monthly['Port_excess'].mean()*12:.2f}%")

# ============================================================================
# 6.2 — Cumulative Returns and Drawdown (Base Portfolio)
# ============================================================================

fig, axes = plt.subplots(2, 1, figsize=(14, 9), sharex=True,
                         gridspec_kw={'height_ratios': [2, 1]})

# Panel A: Cumulative returns
etf_colours = {'AVWS': '#1f77b4', 'AVWC': '#ff7f0e', 'AVEM': '#2ca02c', 'IWMO': '#9467bd'}

for short in base_labels + ['IWMO']:
    cum_etf = (1 + port_data[f'{short}_ret'] / 100).cumprod()
    style = dict(alpha=0.5, linewidth=1.2) if short != 'IWMO' else dict(alpha=0.4, linewidth=1, linestyle='--')
    axes[0].plot(cum_etf.index, cum_etf.values, label=short, color=etf_colours[short], **style)

cum_port = (1 + port_data['Port_ret'] / 100).cumprod()
axes[0].plot(cum_port.index, cum_port.values, label='Portfolio (50/40/10)',
             color='black', linewidth=2.5)
axes[0].set_ylabel('Cumulative Return ($1 invested)')
axes[0].legend(loc='upper left', fontsize=9)
axes[0].set_title('Panel A: Cumulative Returns')
axes[0].grid(True, alpha=0.3)
axes[0].axhline(1, color='grey', linestyle='--', alpha=0.5)

# Panel B: Drawdown
cum_max = cum_port.cummax()
drawdown = (cum_port - cum_max) / cum_max * 100
axes[1].fill_between(drawdown.index, drawdown.values, 0, color='red', alpha=0.3)
axes[1].plot(drawdown.index, drawdown.values, color='red', linewidth=0.8)
axes[1].set_ylabel('Drawdown (%)')
axes[1].set_title('Panel B: Portfolio Drawdown')
axes[1].grid(True, alpha=0.3)

max_dd = drawdown.min()
max_dd_date = drawdown.idxmin()
axes[1].annotate(f'Max DD: {max_dd:.1f}%',
                 xy=(max_dd_date, max_dd),
                 xytext=(max_dd_date + pd.Timedelta(days=30), max_dd + 2),
                 arrowprops=dict(arrowstyle='->', color='darkred'),
                 fontsize=10, color='darkred', fontweight='bold')

fig.suptitle('Base Portfolio: 50% AVWS / 40% AVWC / 10% AVEM',
             fontsize=14, fontweight='bold', y=1.01)
plt.tight_layout()
plt.show()

print(f"\n  Maximum drawdown: {max_dd:.2f}% on {max_dd_date.date()}")
print(f"  Days from peak to trough: "
      f"{(max_dd_date - cum_port[:max_dd_date].idxmax()).days}")

# ============================================================================
# 7.1 — Base Portfolio Factor Regressions (Monthly, DM Factors)
# ============================================================================

port_label = '50/40/10'

port_reg_data = port_monthly[['Port_excess'] + ff6_cols].copy()
port_reg_data = port_reg_data.rename(columns={'Port_excess': 'Excess_Return'})

# FF3
res_port_3f, lags_3f = run_factor_regression(port_reg_data, ff3_cols, f'Portfolio ({port_label})')
print("=" * 70)
print(f"BASE PORTFOLIO ({port_label}) — Fama-French 3-Factor (Monthly)")
print("=" * 70)
print_regression_results(res_port_3f, ff3_cols, f'Portfolio [3-Factor, DM Monthly]', lags_3f, freq='monthly')

print("\n")

# FF6
res_port_6f, lags_6f = run_factor_regression(port_reg_data, ff6_cols, f'Portfolio ({port_label})')
print("=" * 70)
print(f"BASE PORTFOLIO ({port_label}) — FF5 + Momentum (Monthly)")
print("=" * 70)
print_regression_results(res_port_6f, ff6_cols, f'Portfolio [6-Factor, DM Monthly]', lags_6f, freq='monthly')

print("\n")

# 5-Factor without CMA (consistent with individual ETF analysis in Section 4)
res_port_noCMA, lags_noCMA = run_factor_regression(port_reg_data, noCMA_cols, f'Portfolio ({port_label})')
print("=" * 70)
print(f"BASE PORTFOLIO ({port_label}) — 5-Factor without CMA (Monthly)")
print("=" * 70)
print_regression_results(res_port_noCMA, noCMA_cols, f'Portfolio [5-Factor no CMA, DM Monthly]', lags_noCMA, freq='monthly')

# Compare HML coefficient stability
hml_idx_6f = ff6_cols.index('HML') + 1
hml_idx_5f = noCMA_cols.index('HML') + 1
print(f"\n  HML coefficient:  6-Factor = {res_port_6f.params[hml_idx_6f]:.4f} (SE {res_port_6f.bse[hml_idx_6f]:.4f})")
print(f"                    5F-noCMA = {res_port_noCMA.params[hml_idx_5f]:.4f} (SE {res_port_noCMA.bse[hml_idx_5f]:.4f})")
print(f"  R² change:        {res_port_6f.rsquared:.4f} → {res_port_noCMA.rsquared:.4f} (Δ = {res_port_noCMA.rsquared - res_port_6f.rsquared:+.5f})")

# ============================================================================
# 7.2 — Verify Linearity: Weighted-Average Betas vs. Portfolio Betas
# ============================================================================

# Re-run individual regressions on the SAME DM factors & common sample
etf_results_common = {}
weight_map = base_weights.copy()
weight_map['IS3R.DE'] = 0.0  # IWMO not in base portfolio, but include for later

for ticker in etf_tickers:
    short = 'IWMO' if ticker == 'IS3R.DE' else ticker.replace('.DE', '')
    temp_df = port_monthly[[f'{short}_excess'] + ff6_cols].copy()
    temp_df = temp_df.rename(columns={f'{short}_excess': 'Excess_Return'})
    res_i, _ = run_factor_regression(temp_df, ff6_cols, ticker)
    etf_results_common[ticker] = res_i

param_names = ['Alpha'] + ff6_cols
print("=" * 75)
print("Linearity Verification: Portfolio β  vs.  Σ wᵢ · βᵢ  (FF6, DM Monthly)")
print("  Weights: 50% AVWS / 40% AVWC / 10% AVEM")
print("  All ETFs regressed on DM monthly factors, same common sample")
print("=" * 75)
print(f"  {'Factor':<10s}  {'Portfolio β':>12s}  {'Weighted Avg':>12s}  {'Difference':>12s}")
print("  " + "-" * 50)

for j, factor in enumerate(param_names):
    beta_port = res_port_6f.params[j]
    beta_wavg = sum(weight_map.get(t, 0) * etf_results_common[t].params[j]
                    for t in base_tickers)
    diff = beta_port - beta_wavg
    print(f"  {factor:<10s}  {beta_port:>12.6f}  {beta_wavg:>12.6f}  {diff:>12.2e}")

print()
print("  Differences are at machine precision — linearity of OLS confirmed.")
print("  When all ETFs use the same factors and sample, portfolio factor")
print("  exposure is fully determined by individual exposures and weights.")

# ============================================================================
# 7.3 — Portfolio Factor Exposure Decomposition (Base 50/40/10)
# ============================================================================

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

factor_labels = ff6_cols
colours = {'AVWS.DE': '#1f77b4', 'AVWC.DE': '#ff7f0e', 'AVEM.DE': '#2ca02c'}

# Panel A: Stacked contribution to portfolio betas
contributions = {}
for factor in factor_labels:
    j = ff6_cols.index(factor) + 1
    contributions[factor] = [weight_map.get(t, 0) * etf_results_common[t].params[j]
                             for t in base_tickers]

x = np.arange(len(factor_labels))
width = 0.5
bottom_pos = np.zeros(len(factor_labels))
bottom_neg = np.zeros(len(factor_labels))

for i, ticker in enumerate(base_tickers):
    short = ticker.replace('.DE', '')
    vals = np.array([contributions[f][i] for f in factor_labels])
    pos_vals = np.where(vals >= 0, vals, 0)
    neg_vals = np.where(vals < 0, vals, 0)
    axes[0].bar(x, pos_vals, width, bottom=bottom_pos,
                label=f'{short} ({weight_map[ticker]:.0%})',
                color=colours[ticker], alpha=0.8, edgecolor='white', linewidth=0.5)
    axes[0].bar(x, neg_vals, width, bottom=bottom_neg,
                color=colours[ticker], alpha=0.8, edgecolor='white', linewidth=0.5)
    bottom_pos += pos_vals
    bottom_neg += neg_vals

# Portfolio beta markers
port_betas = [res_port_6f.params[ff6_cols.index(f) + 1] for f in factor_labels]
axes[0].scatter(x, port_betas, color='black', zorder=5, s=60, marker='D',
                label='Portfolio β (total)')

axes[0].set_xticks(x)
axes[0].set_xticklabels(factor_labels, rotation=45, ha='right')
axes[0].set_ylabel('Beta Coefficient')
axes[0].set_title('Panel A: Decomposition of Portfolio Betas\n(DM Monthly FF6, common sample)',
                   fontweight='bold', fontsize=11)
axes[0].legend(loc='upper right', fontsize=8)
axes[0].axhline(0, color='black', linewidth=0.5)
axes[0].grid(True, alpha=0.2, axis='y')

# Panel B: Alpha decomposition
alpha_contributions = [weight_map.get(t, 0) * etf_results_common[t].params[0] * 12
                       for t in base_tickers]
port_alpha_ann = res_port_6f.params[0] * 12

bars = axes[1].bar(base_labels, alpha_contributions,
                   color=[colours[t] for t in base_tickers], alpha=0.8, edgecolor='white')
axes[1].axhline(port_alpha_ann, color='black', linewidth=2, linestyle='--',
                label=f'Portfolio α = {port_alpha_ann:.2f}% ann.')
axes[1].set_ylabel('Annualised Alpha Contribution (%)')
axes[1].set_title('Panel B: Alpha Contributions to Portfolio\n(weight × ETF alpha, DM Monthly FF6)',
                   fontweight='bold', fontsize=11)
axes[1].legend(fontsize=9)
axes[1].grid(True, alpha=0.2, axis='y')

for bar, val in zip(bars, alpha_contributions):
    ypos = bar.get_height() + 0.05 if bar.get_height() >= 0 else bar.get_height() - 0.15
    axes[1].text(bar.get_x() + bar.get_width()/2, ypos, f'{val:.2f}%',
                 ha='center', fontsize=10, fontweight='bold')

plt.suptitle('Portfolio Factor Exposure Decomposition (50% AVWS / 40% AVWC / 10% AVEM)',
             fontsize=13, fontweight='bold', y=1.03)
plt.tight_layout()
plt.show()

# ============================================================================
# 7.4 — Portfolio Residual Diagnostics (Base 50/40/10)
# ============================================================================
# For consistency with the individual ETF diagnostics in Section 4, we check
# the portfolio regression residuals for autocorrelation and normality.

resid = res_port_6f.resid
dates = port_monthly.loc[common_months].index[:len(resid)]
T = len(resid)
max_acf_lags = min(8, T // 2 - 1)

fig, axes = plt.subplots(1, 3, figsize=(18, 4))

# Residual time series
colours = ['steelblue' if r >= 0 else 'salmon' for r in resid]
axes[0].bar(dates, resid, width=25, color=colours, alpha=0.8, edgecolor='none')
axes[0].axhline(0, color='red', linestyle='--', linewidth=0.8)
axes[0].set_title('Monthly Residuals (FF6, DM)', fontweight='bold')
axes[0].set_ylabel('Residual (%)')
axes[0].grid(True, alpha=0.3)
axes[0].tick_params(axis='x', rotation=30)

# Residual ACF
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(resid, ax=axes[1], lags=max_acf_lags, alpha=0.05)
axes[1].set_title('Residual ACF', fontweight='bold')
axes[1].grid(True, alpha=0.3)

# QQ Plot
stats.probplot(resid, dist='norm', plot=axes[2])
axes[2].set_title('QQ Plot', fontweight='bold')
axes[2].grid(True, alpha=0.3)

fig.suptitle(f'Portfolio (50/40/10) — Residual Diagnostics (N={T})',
             fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

# Formal tests
max_lb = min(4, T // 3)
lb_result = sm.stats.acorr_ljungbox(resid, lags=[max_lb], return_df=True)
lb, lb_p = lb_result.iloc[0]['lb_stat'], lb_result.iloc[0]['lb_pvalue']
jb_stat, jb_p = stats.jarque_bera(resid)

print(f"  Ljung-Box({max_lb}): stat={lb:.2f}, p={lb_p:.4f}")
print(f"  Jarque-Bera:     stat={jb_stat:.2f}, p={jb_p:.4f}")
print(f"\n  With only {T} monthly observations, these tests have low power.")
print("  Newey-West SEs remain valid regardless.")

# ============================================================================
# 8.1 — Find the Optimal IWMO Allocation (preserving 5:4:1 base ratios)
# ============================================================================
from scipy.optimize import minimize_scalar

# Monthly excess returns for all 4 ETFs on the common sample
all_labels = ['AVWS', 'AVWC', 'AVEM', 'IWMO']
all_tickers = ['AVWS.DE', 'AVWC.DE', 'AVEM.DE', 'IS3R.DE']
excess_cols = [f'{l}_excess' for l in all_labels]
monthly_ret = port_monthly[excess_cols].values   # T × 4

mu  = monthly_ret.mean(axis=0)
cov = np.cov(monthly_ret, rowvar=False)

# Base ratios (Avantis only, before IWMO)
base_ratios = np.array([0.50, 0.40, 0.10])

def weights_for_iwmo(w_iwmo):
    """Given IWMO weight, return full 4-asset weight vector."""
    w_base = (1 - w_iwmo) * base_ratios
    return np.append(w_base, w_iwmo)

def neg_sharpe_iwmo(w_iwmo):
    """Negative annualised Sharpe as a function of IWMO weight alone."""
    w = weights_for_iwmo(w_iwmo)
    port_mu  = w @ mu
    port_vol = np.sqrt(w @ cov @ w)
    return -(port_mu / port_vol * np.sqrt(12)) if port_vol > 0 else 0

# Optimise over w_iwmo ∈ [0, 0.30]  (cap at 30% so the base trio dominates)
result = minimize_scalar(neg_sharpe_iwmo, bounds=(0, 0.30), method='bounded')
w_iwmo_opt = result.x

# Build comparison table
print("=" * 75)
print("OPTIMAL IWMO ALLOCATION (preserving 5:4:1 Avantis ratios)")
print("=" * 75)

print(f"\n{'IWMO %':>8s}  {'AVWS':>7s}  {'AVWC':>7s}  {'AVEM':>7s}  "
      f"{'Ann.Ret':>8s}  {'Ann.Vol':>8s}  {'Sharpe':>7s}")
print("-" * 65)

for w_test in [0.00, 0.05, 0.10, 0.15, 0.20, 0.25, w_iwmo_opt]:
    w = weights_for_iwmo(w_test)
    r = w @ mu * 12
    s = np.sqrt(w @ cov @ w) * np.sqrt(12)
    sr = r / s if s > 0 else 0
    marker = "  ← optimal" if abs(w_test - w_iwmo_opt) < 0.001 else ""
    print(f"  {w_test:>5.1%}   {w[0]:>6.1%}  {w[1]:>6.1%}  {w[2]:>6.1%}  "
          f"{r:>7.2f}%  {s:>7.2f}%  {sr:>6.3f}{marker}")

# Build the enhanced portfolio
w_enhanced = weights_for_iwmo(w_iwmo_opt)
enhanced_weights = dict(zip(all_tickers, w_enhanced))

print(f"\n→ Enhanced portfolio weights:")
for t, l in zip(all_tickers, all_labels):
    print(f"  {l}: {enhanced_weights[t]:.1%}")

# Compute enhanced portfolio daily & monthly returns
port_data['Enhanced_ret'] = sum(enhanced_weights[t] * port_data[f'{l}_ret']
                                for t, l in zip(all_tickers, all_labels))
port_data['Enhanced_excess'] = sum(enhanced_weights[t] * port_data[f'{l}_excess']
                                   for t, l in zip(all_tickers, all_labels))

port_monthly['Enhanced_ret'] = sum(enhanced_weights[t] * port_monthly[f'{l}_ret']
                                   for t, l in zip(all_tickers, all_labels))
port_monthly['Enhanced_excess'] = sum(enhanced_weights[t] * port_monthly[f'{l}_excess']
                                      for t, l in zip(all_tickers, all_labels))

# Summary comparison
print(f"\n{'':32s}  {'Base 50/40/10':>14s}  {'+ IWMO':>14s}")
print("-" * 65)
r_base = port_data['Port_ret']
r_enh  = port_data['Enhanced_ret']
print(f"  {'Ann. return (%)':32s}  {r_base.mean()*252:>13.2f}%  {r_enh.mean()*252:>13.2f}%")
print(f"  {'Ann. volatility (%)':32s}  {r_base.std()*np.sqrt(252):>13.2f}%  {r_enh.std()*np.sqrt(252):>13.2f}%")
sr_base = r_base.mean() / r_base.std() * np.sqrt(252)
sr_enh  = r_enh.mean() / r_enh.std() * np.sqrt(252)
print(f"  {'Sharpe ratio':32s}  {sr_base:>14.3f}  {sr_enh:>14.3f}")
print(f"  {'UMD exposure (monthly β)':32s}  {'~0 (no IWMO)':>14s}  {'see below':>14s}")

# ============================================================================
# 8.2 — Side-by-Side: Base vs. Enhanced Portfolio
# ============================================================================

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

# Panel A: Cumulative returns comparison
cum_base = (1 + port_data['Port_ret'] / 100).cumprod()
cum_enh  = (1 + port_data['Enhanced_ret'] / 100).cumprod()

axes[0].plot(cum_base.index, cum_base.values, label='Base (50/40/10)',
             color='#1f77b4', linewidth=2)
axes[0].plot(cum_enh.index, cum_enh.values, label=f'Enhanced (+{w_iwmo_opt:.0%} IWMO)',
             color='#d62728', linewidth=2)

# Also show individual IWMO for context
cum_iwmo = (1 + port_data['IWMO_ret'] / 100).cumprod()
axes[0].plot(cum_iwmo.index, cum_iwmo.values, label='IWMO (standalone)',
             color='#9467bd', linewidth=1, alpha=0.4, linestyle='--')

axes[0].axhline(1, color='grey', linestyle='--', alpha=0.5)
axes[0].set_ylabel('Cumulative Return ($1 invested)')
axes[0].set_title('Panel A: Cumulative Returns — Base vs. Enhanced', fontweight='bold')
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)

# Panel B: Factor exposure comparison (bar chart)
# Run regression on enhanced portfolio
enh_reg_data = port_monthly[['Enhanced_excess'] + ff6_cols].copy()
enh_reg_data = enh_reg_data.rename(columns={'Enhanced_excess': 'Excess_Return'})
res_enh_6f, _ = run_factor_regression(enh_reg_data, ff6_cols, 'Enhanced')

x = np.arange(len(ff6_cols))
width = 0.35
base_betas = [res_port_6f.params[ff6_cols.index(f) + 1] for f in ff6_cols]
enh_betas  = [res_enh_6f.params[ff6_cols.index(f) + 1] for f in ff6_cols]

axes[1].bar(x - width/2, base_betas, width, label='Base (50/40/10)',
            color='#1f77b4', alpha=0.8, edgecolor='white')
axes[1].bar(x + width/2, enh_betas, width, label=f'Enhanced (+{w_iwmo_opt:.0%} IWMO)',
            color='#d62728', alpha=0.8, edgecolor='white')

axes[1].set_xticks(x)
axes[1].set_xticklabels(ff6_cols, rotation=45, ha='right')
axes[1].set_ylabel('Beta Coefficient')
axes[1].set_title('Panel B: Factor Exposures — Base vs. Enhanced', fontweight='bold')
axes[1].legend(fontsize=9)
axes[1].axhline(0, color='black', linewidth=0.5)
axes[1].grid(True, alpha=0.2, axis='y')

plt.suptitle(f'Impact of Adding {w_iwmo_opt:.0%} IWMO to the 50/40/10 Portfolio',
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

# Print the enhanced portfolio regression
print()
print("=" * 70)
print(f"ENHANCED PORTFOLIO — FF5 + Momentum (Monthly)")
print("=" * 70)
print_regression_results(res_enh_6f, ff6_cols, 'Enhanced Portfolio [6-Factor, DM Monthly]', 
                         lags_6f, freq='monthly')

# Key comparison
print("\n" + "-" * 60)
print("KEY CHANGES from adding IWMO:")
print("-" * 60)
for f in ff6_cols:
    j = ff6_cols.index(f) + 1
    b_base = res_port_6f.params[j]
    b_enh  = res_enh_6f.params[j]
    delta  = b_enh - b_base
    if abs(delta) > 0.01:
        print(f"  {f:<8s}: {b_base:>+.3f} → {b_enh:>+.3f}  (Δ = {delta:>+.3f})")