#!/usr/bin/env python3 # ============================================================================= # Author: Simon-Pierre Boucher # Contact: contact@spboucher.ai # ============================================================================= """Step 04 — RQ3: Implied vs realized correlation divergence as a stress predictor. Stage 1 (raw-dependent): build the daily implied-correlation index from SPX and constituent ATM IVs, the 22-day realized correlation from daily returns, merge with VIX, define stress events, and save the combined series. Stage 2 (parquet-only): predictive stress regressions and crisis-window analysis. When the raw stores are absent, stage 1 is skipped and stage 2 runs from the shipped ``correlation_divergence.parquet``. Inputs : options.duckdb + index_5min.duckdb (stage 1, optional), data/processed/realized_vol.parquet (stage 1), data/processed/correlation_divergence.parquet (stage 2 fallback) Outputs: data/processed/correlation_divergence.parquet (stage 1), results/rq3_stress_prediction.csv, results/rq3_crisis_analysis.csv """ import warnings import numpy as np import pandas as pd import _bootstrap # noqa: F401 from wp7 import config from wp7.data_io import (RawDataUnavailableError, load_correlation_divergence, load_realized_vol, open_raw_db) from wp7.econometrics import add_constant, hc1_tstats, ols, r_squared, standardize warnings.filterwarnings('ignore') STRESS_HORIZONS = [5, 10, 20] # -------------------------------------------------------------------------- # Stage 1 — build the correlation-divergence dataset (needs raw stores) # -------------------------------------------------------------------------- def build_correlation_dataset() -> pd.DataFrame: """Implied correlation (CBOE methodology), realized correlation, VIX.""" constituents = config.SPX_CONSTITUENTS const_str = ",".join([f"'{t}'" for t in constituents]) # 3A. SPX and constituent 30-day ATM implied volatilities opt_con = open_raw_db("options") print(" Extracting SPX implied volatility...") spx_iv = opt_con.execute(""" SELECT trade_date, AVG(CASE WHEN bid_iv > 0 AND ask_iv > 0 THEN (bid_iv + ask_iv)/2.0 WHEN ask_iv > 0 THEN ask_iv ELSE bid_iv END) AS spx_iv_atm FROM option_chain WHERE ticker = 'SPX' AND call_put = 'c' AND ABS(delta - 0.5) < 0.1 AND (expiry_date - trade_date) BETWEEN 20 AND 40 AND ask_price > 0 GROUP BY trade_date ORDER BY trade_date """).fetchdf() print(f" SPX IV: {len(spx_iv)} days") print(" Extracting constituent implied volatilities...") const_iv = opt_con.execute(f""" SELECT ticker, trade_date, AVG(CASE WHEN bid_iv > 0 AND ask_iv > 0 THEN (bid_iv + ask_iv)/2.0 WHEN ask_iv > 0 THEN ask_iv ELSE bid_iv END) AS iv_atm FROM option_chain WHERE ticker IN ({const_str}) AND call_put = 'c' AND ABS(delta - 0.5) < 0.1 AND (expiry_date - trade_date) BETWEEN 20 AND 40 AND ask_price > 0 GROUP BY ticker, trade_date ORDER BY trade_date """).fetchdf() print(f" Constituent IVs: {len(const_iv)} rows") opt_con.close() # Implied correlation: rho = (sigma²_idx − n⁻¹·avg(sigma²_i)) / # ((1 − n⁻¹)·avg(sigma_i)²) const_iv['iv_sq'] = const_iv['iv_atm'] ** 2 avg_const_iv = const_iv.groupby('trade_date').agg( avg_const_iv_sq=('iv_sq', 'mean'), avg_const_iv=('iv_atm', 'mean'), n_constituents=('ticker', 'nunique'), ).reset_index() ic_df = pd.merge(spx_iv, avg_const_iv, on='trade_date', how='inner') ic_df['spx_iv_sq'] = ic_df['spx_iv_atm'] ** 2 n = ic_df['n_constituents'] ic_df['implied_corr'] = (ic_df['spx_iv_sq'] - (1.0 / n) * ic_df['avg_const_iv_sq']) / \ ((1.0 - 1.0 / n) * ic_df['avg_const_iv'] ** 2) ic_df['implied_corr'] = ic_df['implied_corr'].clip(0, 1) print(f" Implied correlation computed: {len(ic_df)} days") # 3B. Realized correlation: 22-day rolling average pairwise correlation # of constituent daily returns. print("\n Computing realized correlations from daily returns...") rv = load_realized_vol() rv = rv[rv['ticker'].isin(constituents)] ret_pivot = rv.pivot_table(index='trade_date', columns='ticker', values='daily_return', aggfunc='first') ret_pivot = ret_pivot.dropna(axis=0, thresh=10) rolling_corr = [] dates = ret_pivot.index.tolist() for i in range(21, len(dates)): window = ret_pivot.iloc[i - 21:i + 1] corr_matrix = window.corr() mask = np.triu(np.ones(corr_matrix.shape, dtype=bool), k=1) rolling_corr.append({ 'trade_date': dates[i], 'realized_corr': corr_matrix.values[mask].mean(), 'n_pairs': mask.sum(), }) realized_corr_df = pd.DataFrame(rolling_corr) print(f" Realized correlation computed: {len(realized_corr_df)} days") # 3C. Merge implied and realized correlation ic_df['trade_date'] = pd.to_datetime(ic_df['trade_date']) corr_merged = pd.merge(ic_df, realized_corr_df, on='trade_date', how='inner') corr_merged['corr_divergence'] = corr_merged['implied_corr'] - corr_merged['realized_corr'] corr_merged['corr_ratio'] = corr_merged['implied_corr'] / \ corr_merged['realized_corr'].clip(0.01) print(f" Merged correlation data: {len(corr_merged)} days") print(f" Avg implied corr: {corr_merged['implied_corr'].mean():.4f}") print(f" Avg realized corr: {corr_merged['realized_corr'].mean():.4f}") print(f" Avg divergence: {corr_merged['corr_divergence'].mean():.4f}") # 3D. Stress events from VIX levels and 5-day VIX spikes con = open_raw_db("indices_5min") vix_daily = con.execute(""" SELECT CAST(datetime AS DATE) AS trade_date, LAST(close) AS vix_close, MAX(close) AS vix_high FROM ohlcv WHERE symbol = 'VIX' GROUP BY CAST(datetime AS DATE) ORDER BY trade_date """).fetchdf() con.close() vix_daily['trade_date'] = pd.to_datetime(vix_daily['trade_date']) corr_merged = pd.merge(corr_merged, vix_daily, on='trade_date', how='left') corr_merged['vix_5d_change'] = corr_merged['vix_close'].pct_change(5) corr_merged['stress_event'] = ((corr_merged['vix_close'] > 25) | (corr_merged['vix_5d_change'] > 0.20)).astype(int) for horizon in STRESS_HORIZONS: corr_merged[f'stress_fwd_{horizon}d'] = corr_merged['stress_event'].rolling( horizon, min_periods=1).max().shift(-horizon) corr_merged.to_parquet(config.CORRELATION_DIVERGENCE_PARQUET, index=False) return corr_merged # -------------------------------------------------------------------------- # Stage 2 — stress regressions and crisis windows (parquet-only) # -------------------------------------------------------------------------- def stress_regressions(corr_merged: pd.DataFrame) -> None: """Linear probability models: divergence measures → forward stress.""" print("\n--- STRESS PREDICTION REGRESSIONS ---") features = ['corr_divergence', 'corr_ratio', 'implied_corr', 'realized_corr', 'spx_iv_atm', 'vix_close'] rows = [] for horizon in STRESS_HORIZONS: target = f'stress_fwd_{horizon}d' sub = corr_merged[features + [target]].dropna() if len(sub) < 100: continue X = add_constant(standardize(sub[features].values)) y = sub[target].values coefs = ols(X, y) r2 = r_squared(y, X @ coefs) _, t_stats = hc1_tstats(X, y, coefs, len(features)) for i, fname in enumerate(['const'] + features): rows.append({ 'horizon': f'{horizon}d', 'variable': fname, 'coefficient': coefs[i], 't_stat': t_stats[i], 'significant': abs(t_stats[i]) > 1.96, }) print(f"\n {horizon}-day stress prediction: R²={r2:.4f}, N={len(sub)}") for i, fname in enumerate(['const'] + features): sig = "*" if abs(t_stats[i]) > 1.96 else "" sig2 = "*" if abs(t_stats[i]) > 2.576 else "" print(f" {fname:20s}: β={coefs[i]:8.4f}, t={t_stats[i]:7.3f} {sig}{sig2}") pd.DataFrame(rows).to_csv(config.RESULTS_DIR / "rq3_stress_prediction.csv", index=False) def crisis_analysis(corr_merged: pd.DataFrame) -> None: """Correlation divergence in 30-day pre / 10-day post windows around crises.""" print("\n--- CORRELATION DIVERGENCE AROUND CRISES ---") rows = [] for name, date_str in config.CRISES.items(): date = pd.Timestamp(date_str) pre = corr_merged[(corr_merged['trade_date'] >= date - pd.Timedelta(days=30)) & (corr_merged['trade_date'] < date)] post = corr_merged[(corr_merged['trade_date'] >= date) & (corr_merged['trade_date'] < date + pd.Timedelta(days=10))] if len(pre) > 0 and len(post) > 0: rows.append({ 'crisis': name, 'pre_impl_corr': pre['implied_corr'].mean(), 'pre_real_corr': pre['realized_corr'].mean(), 'pre_divergence': pre['corr_divergence'].mean(), 'post_impl_corr': post['implied_corr'].mean(), 'post_real_corr': post['realized_corr'].mean(), 'post_divergence': post['corr_divergence'].mean(), 'divergence_change': post['corr_divergence'].mean() - pre['corr_divergence'].mean(), 'pre_vix': pre['vix_close'].mean() if 'vix_close' in pre.columns else np.nan, 'post_vix': post['vix_close'].mean() if 'vix_close' in post.columns else np.nan, }) print(f" {name}:") print(f" Pre: IC={pre['implied_corr'].mean():.4f}, " f"RC={pre['realized_corr'].mean():.4f}, " f"Div={pre['corr_divergence'].mean():.4f}") print(f" Post: IC={post['implied_corr'].mean():.4f}, " f"RC={post['realized_corr'].mean():.4f}, " f"Div={post['corr_divergence'].mean():.4f}") pd.DataFrame(rows).to_csv(config.RESULTS_DIR / "rq3_crisis_analysis.csv", index=False) def main(): print("=" * 70) print("RQ3: IMPLIED vs REALIZED CORRELATION DIVERGENCE") print("=" * 70) config.ensure_output_dirs() try: corr_merged = build_correlation_dataset() except RawDataUnavailableError as exc: print(f"\n[Stage 1 skipped — raw stores unavailable]\n{exc}\n") print("Loading shipped correlation_divergence.parquet instead.") corr_merged = load_correlation_divergence() stress_regressions(corr_merged) crisis_analysis(corr_merged) print("\nRQ3 COMPLETE.") if __name__ == "__main__": main()