spb/wp7_uqo Public

1#!/usr/bin/env python32# =============================================================================3# Author: Simon-Pierre Boucher4# Contact: contact@spboucher.ai5# =============================================================================6"""Step 03 — RQ2: IV term structure + smile vs GARCH/HAR-RV for RV forecasting.78In-sample panel R² for five nested models, per-ticker rolling out-of-sample9evaluation, and Diebold-Mariano tests of HAR-RV vs HAR-RV + IV surface.1011Inputs : data/processed/merged_options_rv.parquet12Outputs: results/rq2_insample.csv, results/rq2_oos_results.csv,13         results/rq2_diebold_mariano.csv14"""1516import warnings1718import numpy as np19import pandas as pd2021import _bootstrap  # noqa: F40122from wp7 import config23from wp7.data_io import load_merged24from wp7.econometrics import add_constant, ols, standardize, winsorize25from wp7.forecasting import diebold_mariano, ols_forecast, rolling_evaluation2627warnings.filterwarnings('ignore')2829TARGETS = {'1-Day RV': 'rv_fwd_1d', '5-Day RV': 'rv_fwd_5d'}3031WINSORIZE_COLS = ['rv_fwd_1d', 'rv_fwd_5d', 'rv_lag1', 'rv_w', 'rv_m', 'sq_return',32                  'iv_atm_30d', 'iv_atm_90d', 'iv_term_slope', 'iv_skew_25d',33                  'implied_skewness', 'implied_kurtosis_proxy']343536def main():37    print("=" * 70)38    print("RQ2: IV SURFACE vs GARCH/HAR-RV FOR RV FORECASTING")39    print("=" * 70)40    config.ensure_output_dirs()4142    df = load_merged()43    df = df.sort_values(['ticker', 'trade_date']).reset_index(drop=True)44    print(f"Loaded {len(df):,} rows")4546    # GARCH(1,1) proxy regressor: yesterday's squared daily return47    df['sq_return'] = df['daily_return'] ** 24849    # Per-ticker winsorization of every model input and target50    for col in WINSORIZE_COLS:51        if col in df.columns:52            df[col] = df.groupby('ticker')[col].transform(lambda x: winsorize(x))5354    models = config.RQ2_MODELS5556    # ── In-sample panel R² ──57    print("\n--- IN-SAMPLE PANEL REGRESSIONS ---")58    insample_results = []59    for target_label, target_col in TARGETS.items():60        for model_name, features in models.items():61            sub = df[features + [target_col, 'ticker']].dropna()62            if len(sub) < 100:63                continue6465            X = add_constant(standardize(sub[features].values))66            y = sub[target_col].values67            coefs = ols(X, y)68            y_pred = X @ coefs69            ss_res = np.sum((y - y_pred) ** 2)70            ss_tot = np.sum((y - y.mean()) ** 2)71            r2 = 1 - ss_res / ss_tot72            adj_r2 = 1 - (1 - r2) * (len(y) - 1) / (len(y) - len(features) - 1)7374            insample_results.append({75                'target': target_label,76                'model': model_name,77                'r2': r2,78                'adj_r2': adj_r2,79                'n_obs': len(sub),80                'n_features': len(features),81            })82            print(f"  {target_label} | {model_name}: R²={r2:.6f}, "83                  f"Adj-R²={adj_r2:.6f}, N={len(sub):,}")8485    insample_df = pd.DataFrame(insample_results)86    insample_df.to_csv(config.RESULTS_DIR / "rq2_insample.csv", index=False)8788    # ── Out-of-sample rolling evaluation (30 tickers with most data) ──89    print("\n--- OUT-OF-SAMPLE ROLLING EVALUATION ---")90    top_tickers = df.groupby('ticker').size().nlargest(30).index.tolist()91    df_top = df[df['ticker'].isin(top_tickers)]9293    oos_results_all = []94    for target_label, target_col in TARGETS.items():95        print(f"\n  {target_label}:")96        oos = rolling_evaluation(df_top, models, target_col, target_label,97                                 window=500, step=250)98        oos_results_all.append(oos)99        if len(oos) > 0:100            summary = oos.groupby('model').agg({101                'avg_mse': 'mean', 'avg_mae': 'mean',102                'avg_r2_oos': 'mean', 'avg_qlike': 'mean',103                'n_windows': 'sum',104            }).round(6)105            print(summary.to_string())106107    oos_df = pd.concat(oos_results_all, ignore_index=True)108    oos_df.to_csv(config.RESULTS_DIR / "rq2_oos_results.csv", index=False)109110    # ── Diebold-Mariano: HAR-RV vs HAR-RV + IV surface (10 largest tickers) ──111    print("\n--- MODEL COMPARISON: DIEBOLD-MARIANO TESTS ---")112    f_har = ['rv_lag1', 'rv_w', 'rv_m']113    f_full = f_har + ['iv_atm_30d', 'iv_atm_90d', 'iv_term_slope',114                      'iv_skew_25d', 'implied_skewness', 'implied_kurtosis_proxy']115    dm_results = []116    for target_label, target_col in TARGETS.items():117        for ticker in top_tickers[:10]:118            td = df[df['ticker'] == ticker][f_full + [target_col]].dropna()119            if len(td) < 600:120                continue121122            train, test = td.iloc[:500], td.iloc[500:]123            r1 = ols_forecast(train, test, f_har, target_col)124            r2 = ols_forecast(train, test, f_full, target_col)125126            # Forecast errors on the test window (raw, before truncation at 0)127            X1 = np.column_stack([np.ones(len(test)), test[f_har].values])128            e1 = test[target_col].values - X1 @ r1['coefs']129            X2 = np.column_stack([np.ones(len(test)), test[f_full].values])130            e2 = test[target_col].values - X2 @ r2['coefs']131132            dm_stat = diebold_mariano(e1, e2)133            dm_results.append({134                'ticker': ticker,135                'target': target_label,136                'mse_har': r1['mse'],137                'mse_har_iv': r2['mse'],138                'mse_improvement_pct': (r1['mse'] - r2['mse']) / r1['mse'] * 100,139                'dm_statistic': dm_stat,140                'dm_significant_5pct': abs(dm_stat) > 1.96,141            })142143    dm_df = pd.DataFrame(dm_results)144    dm_df.to_csv(config.RESULTS_DIR / "rq2_diebold_mariano.csv", index=False)145    print(dm_df.to_string(index=False))146147    print("\n" + "=" * 70)148    print("RQ2 SUMMARY")149    print("=" * 70)150    print("\nIn-sample R² comparison:")151    print(insample_df.pivot(index='model', columns='target', values='r2')152          .round(6).to_string())153    print("\nRQ2 COMPLETE.")154155156if __name__ == "__main__":157    main()158