spb/wp7_uqo Public

1#!/usr/bin/env python32# =============================================================================3# Author: Simon-Pierre Boucher4# Contact: contact@spboucher.ai5# =============================================================================6"""Step 11 — Granger causality and VAR analysis.78Lead-lag structure between option-implied measures and returns / realized9volatility: per-ticker bivariate Granger F-tests (5 lags), bivariate VAR(5)10with impulse-response functions, and forecast-error variance decomposition.1112Inputs : data/processed/merged_options_rv.parquet13Outputs: results/granger_causality.csv, results/var_results.csv,14         results/irf_results.csv, results/fevd_results.csv15"""1617import warnings1819import numpy as np20import pandas as pd2122import _bootstrap  # noqa: F40123from wp7 import config24from wp7.data_io import load_merged25from wp7.econometrics import estimate_var, granger_test, winsorize2627warnings.filterwarnings('ignore')2829TEST_PAIRS = [30    ('iv_atm_30d', 'rv_daily', 'ATM_IV → RV'),31    ('rv_daily', 'iv_atm_30d', 'RV → ATM_IV'),32    ('iv_skew_25d', 'rv_daily', 'Skew → RV'),33    ('rv_daily', 'iv_skew_25d', 'RV → Skew'),34    ('iv_atm_30d', 'daily_return', 'ATM_IV → Return'),35    ('daily_return', 'iv_atm_30d', 'Return → ATM_IV'),36    ('iv_skew_25d', 'daily_return', 'Skew → Return'),37    ('daily_return', 'iv_skew_25d', 'Return → Skew'),38    ('pc_volume_ratio', 'daily_return', 'PC_Ratio → Return'),39    ('daily_return', 'pc_volume_ratio', 'Return → PC_Ratio'),40    ('implied_skewness', 'rv_daily', 'Impl_Skew → RV'),41    ('implied_kurtosis_proxy', 'rv_daily', 'Impl_Kurt → RV'),42]4344IRF_PERIODS = 20454647def main():48    print("=" * 70)49    print("GRANGER CAUSALITY & VAR ANALYSIS")50    print("=" * 70)51    config.ensure_output_dirs()5253    df = load_merged()54    df = df.sort_values(['ticker', 'trade_date'])5556    # ── A. Granger causality tests (per ticker, averaged) ──57    print("\n--- A. GRANGER CAUSALITY TESTS ---")58    granger_results = []59    for x_var, y_var, label in TEST_PAIRS:60        ticker_results = []61        for ticker in df['ticker'].unique():62            td = df[df['ticker'] == ticker][[x_var, y_var]].dropna()63            if len(td) < 100:64                continue65            x = winsorize(td[x_var]).values66            y = winsorize(td[y_var]).values67            res = granger_test(y, x, max_lag=5)68            if res:69                ticker_results.append(res)7071        if ticker_results:72            avg_f = np.mean([r['f_stat'] for r in ticker_results])73            avg_p = np.mean([r['p_value'] for r in ticker_results])74            pct_sig = np.mean([1 if r['p_value'] < 0.05 else 0 for r in ticker_results])75            granger_results.append({76                'test': label,77                'x_causes_y': f"{x_var} → {y_var}",78                'avg_f_stat': avg_f,79                'avg_p_value': avg_p,80                'pct_significant_5pct': pct_sig,81                'n_tickers': len(ticker_results),82            })83            sig_str = "***" if avg_p < 0.01 else \84                ("**" if avg_p < 0.05 else ("*" if avg_p < 0.10 else ""))85            print(f"  {label:25s}: F={avg_f:8.3f}, p={avg_p:.4f}, "86                  f"{pct_sig*100:.1f}% sig  {sig_str}")8788    pd.DataFrame(granger_results).to_csv(89        config.RESULTS_DIR / "granger_causality.csv", index=False)9091    # ── B. Bivariate VAR(5): ATM IV ↔ RV (20 largest tickers) ──92    print("\n--- B. BIVARIATE VAR: IV_ATM ↔ RV (pooled) ---")93    var_results, irf_all = [], []94    top_tickers = df.groupby('ticker').size().nlargest(20).index.tolist()9596    for ticker in top_tickers:97        td = df[df['ticker'] == ticker][['iv_atm_30d', 'rv_daily']].dropna()98        if len(td) < 200:99            continue100        y1 = winsorize(td['iv_atm_30d']).values101        y2 = winsorize(td['rv_daily']).values102        y1 = (y1 - y1.mean()) / y1.std()103        y2 = (y2 - y2.mean()) / y2.std()104105        res = estimate_var(y1, y2, lags=5, irf_periods=IRF_PERIODS)106        if res:107            var_results.append({108                'ticker': ticker,109                'r2_iv_eq': res['r2_eq1'],110                'r2_rv_eq': res['r2_eq2'],111                'n_obs': res['n_obs'],112            })113            for h in range(IRF_PERIODS):114                irf_all.append({115                    'ticker': ticker, 'horizon': h,116                    'iv_to_iv': res['irf'][h, 0, 0],117                    'rv_to_iv': res['irf'][h, 0, 1],118                    'iv_to_rv': res['irf'][h, 1, 0],119                    'rv_to_rv': res['irf'][h, 1, 1],120                })121122    var_df = pd.DataFrame(var_results)123    var_df.to_csv(config.RESULTS_DIR / "var_results.csv", index=False)124    print(f"\n  VAR estimated for {len(var_results)} tickers")125    print(f"  Avg R² (IV equation): {var_df['r2_iv_eq'].mean():.4f}")126    print(f"  Avg R² (RV equation): {var_df['r2_rv_eq'].mean():.4f}")127128    irf_df = pd.DataFrame(irf_all)129    irf_df.to_csv(config.RESULTS_DIR / "irf_results.csv", index=False)130131    avg_irf = irf_df.groupby('horizon')[['iv_to_iv', 'rv_to_iv',132                                         'iv_to_rv', 'rv_to_rv']].mean()133    print("\n  Average Impulse Response Function:")134    print(f"  {'h':>3} {'IV→IV':>8} {'RV→IV':>8} {'IV→RV':>8} {'RV→RV':>8}")135    for h in range(0, IRF_PERIODS, 2):136        row = avg_irf.loc[h]137        print(f"  {h:>3} {row['iv_to_iv']:>8.4f} {row['rv_to_iv']:>8.4f} "138              f"{row['iv_to_rv']:>8.4f} {row['rv_to_rv']:>8.4f}")139140    # ── C. Forecast-error variance decomposition ──141    print("\n--- C. FORECAST ERROR VARIANCE DECOMPOSITION ---")142    fevd_results = []143    for h in range(1, IRF_PERIODS + 1):144        sub_irf = irf_df[irf_df['horizon'] < h]145        avg = sub_irf.groupby('horizon')[['iv_to_rv', 'rv_to_rv']].mean()146        total_var = (avg['iv_to_rv'] ** 2 + avg['rv_to_rv'] ** 2).sum()147        iv_share = (avg['iv_to_rv'] ** 2).sum() / total_var if total_var > 0 else 0148        fevd_results.append({149            'horizon': h,150            'pct_rv_explained_by_iv': iv_share * 100,151            'pct_rv_explained_by_rv': (1 - iv_share) * 100,152        })153154    fevd_df = pd.DataFrame(fevd_results)155    fevd_df.to_csv(config.RESULTS_DIR / "fevd_results.csv", index=False)156    print(f"  {'Horizon':>8} {'% RV by IV':>12} {'% RV by RV':>12}")157    for _, row in fevd_df.iterrows():158        print(f"  {int(row['horizon']):>8} {row['pct_rv_explained_by_iv']:>12.2f} "159              f"{row['pct_rv_explained_by_rv']:>12.2f}")160161    print("\nGRANGER CAUSALITY & VAR ANALYSIS COMPLETE.")162163164if __name__ == "__main__":165    main()166