#!/usr/bin/env python3 # ============================================================================= # Author: Simon-Pierre Boucher # Contact: contact@spboucher.ai # ============================================================================= """Step 06 — RQ5: ML on the full SPX options surface vs VIX for realized- variance forecasting. Builds 21 daily SPX surface features and SPX 5-minute realized variance from the raw stores, then compares — out-of-sample, temporal split 2020-01-01 — the VIX benchmark, OLS models (HAR-RV, IV surface, combined), Random Forest and Gradient Boosting. Entirely raw-dependent: exits gracefully when the raw stores are absent (the shipped result CSVs remain authoritative). Inputs : options.duckdb, index_5min.duckdb Outputs: results/rq5_model_comparison.csv, results/rq5_feature_importance.csv """ import sys import warnings import numpy as np import pandas as pd from numpy.linalg import lstsq import _bootstrap # noqa: F401 from wp7 import config from wp7.data_io import RawDataUnavailableError, open_raw_db warnings.filterwarnings('ignore') SURFACE_FEATURES = [ 'iv_atm_1w', 'iv_atm_2w', 'iv_atm_1m', 'iv_atm_2m', 'iv_atm_3m', 'iv_atm_6m', 'skew_25d_1m', 'skew_25d_3m', 'skew_10d_1m', 'butterfly_1m', 'ts_slope_3m_1m', 'ts_slope_6m_1m', 'total_gamma_oi', 'net_gamma', 'total_vega_oi', 'avg_theta', 'pc_vol_ratio', 'pc_oi_ratio', 'total_volume', 'total_oi', 'avg_spread_pct', ] HAR_FEATURES = ['rv_lag1', 'rv_w', 'rv_m'] TARGETS = {'1-Day RV': 'rv_fwd_1d', '5-Day RV': 'rv_fwd_5d', '22-Day RV': 'rv_fwd_22d'} SPLIT_DATE = pd.Timestamp('2020-01-01') def build_dataset() -> pd.DataFrame: """SPX surface features + SPX 5-min RV + VIX-implied daily variance.""" opt_con = open_raw_db("options") print("\n--- BUILDING SPX OPTIONS SURFACE FEATURES ---") spx_surface = opt_con.execute(""" WITH base AS ( SELECT trade_date, strike, expiry_date, call_put, (expiry_date - trade_date) AS dte, bid_price, ask_price, bid_iv, ask_iv, open_interest, volume, delta, gamma, vega, theta, rho, 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 mid_iv, (bid_price + ask_price) / 2.0 AS mid_price FROM option_chain WHERE ticker = 'SPX' AND ask_price > 0 AND (expiry_date - trade_date) BETWEEN 1 AND 365 ) SELECT trade_date, -- ATM IV by tenor AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 5 AND 10 THEN mid_iv END) AS iv_atm_1w, AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 13 AND 17 THEN mid_iv END) AS iv_atm_2w, AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 25 AND 35 THEN mid_iv END) AS iv_atm_1m, AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 55 AND 65 THEN mid_iv END) AS iv_atm_2m, AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 85 AND 95 THEN mid_iv END) AS iv_atm_3m, AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 170 AND 200 THEN mid_iv END) AS iv_atm_6m, -- Skew by tenor (25-delta put minus 25-delta call) AVG(CASE WHEN call_put='p' AND ABS(delta+0.25)<0.07 AND dte BETWEEN 25 AND 35 THEN mid_iv END) - AVG(CASE WHEN call_put='c' AND ABS(delta-0.25)<0.07 AND dte BETWEEN 25 AND 35 THEN mid_iv END) AS skew_25d_1m, AVG(CASE WHEN call_put='p' AND ABS(delta+0.25)<0.07 AND dte BETWEEN 85 AND 95 THEN mid_iv END) - AVG(CASE WHEN call_put='c' AND ABS(delta-0.25)<0.07 AND dte BETWEEN 85 AND 95 THEN mid_iv END) AS skew_25d_3m, -- Deep OTM skew (10-delta) AVG(CASE WHEN call_put='p' AND ABS(delta+0.10)<0.05 AND dte BETWEEN 25 AND 35 THEN mid_iv END) - AVG(CASE WHEN call_put='c' AND ABS(delta-0.10)<0.05 AND dte BETWEEN 25 AND 35 THEN mid_iv END) AS skew_10d_1m, -- Butterfly (wings / ATM) (AVG(CASE WHEN ABS(delta)<0.15 AND ABS(delta)>0.03 AND dte BETWEEN 25 AND 35 THEN mid_iv END) / NULLIF(AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 25 AND 35 THEN mid_iv END),0) ) AS butterfly_1m, -- Term structure slopes AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 85 AND 95 THEN mid_iv END) - AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 25 AND 35 THEN mid_iv END) AS ts_slope_3m_1m, AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 170 AND 200 THEN mid_iv END) - AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 25 AND 35 THEN mid_iv END) AS ts_slope_6m_1m, -- Aggregate Greeks SUM(gamma * open_interest) AS total_gamma_oi, SUM(CASE WHEN call_put='c' THEN gamma * open_interest ELSE 0 END) - SUM(CASE WHEN call_put='p' THEN gamma * open_interest ELSE 0 END) AS net_gamma, SUM(vega * open_interest) AS total_vega_oi, AVG(theta) AS avg_theta, -- Volume/OI ratios SUM(CASE WHEN call_put='p' THEN volume ELSE 0 END)::DOUBLE / NULLIF(SUM(CASE WHEN call_put='c' THEN volume ELSE 0 END),0) AS pc_vol_ratio, SUM(CASE WHEN call_put='p' THEN open_interest ELSE 0 END)::DOUBLE / NULLIF(SUM(CASE WHEN call_put='c' THEN open_interest ELSE 0 END),0) AS pc_oi_ratio, SUM(volume) AS total_volume, SUM(open_interest) AS total_oi, -- Liquidity proxy AVG(CASE WHEN mid_price > 0 THEN (ask_price - bid_price) / mid_price END) AS avg_spread_pct, COUNT(*) AS n_contracts FROM base GROUP BY trade_date HAVING COUNT(*) >= 50 ORDER BY trade_date """).fetchdf() opt_con.close() print(f" SPX surface features: {len(spx_surface)} days, {spx_surface.shape[1]} columns") idx_con = open_raw_db("indices_5min") spx_rv = idx_con.execute(""" WITH bars AS ( SELECT datetime, close, CAST(datetime AS DATE) AS trade_date, LAG(close) OVER (ORDER BY datetime) AS prev_close FROM ohlcv WHERE symbol = 'SPX' ), rets AS ( SELECT trade_date, LN(close / NULLIF(prev_close, 0)) AS log_ret FROM bars WHERE prev_close > 0 AND close > 0 ) SELECT trade_date, SUM(log_ret * log_ret) AS rv_5min, COUNT(*) AS n_obs, SUM(log_ret) AS daily_ret FROM rets GROUP BY trade_date HAVING COUNT(*) >= 20 ORDER BY trade_date """).fetchdf() vix_daily = idx_con.execute(""" SELECT CAST(datetime AS DATE) AS trade_date, LAST(close) AS vix_close FROM ohlcv WHERE symbol = 'VIX' GROUP BY CAST(datetime AS DATE) ORDER BY trade_date """).fetchdf() idx_con.close() for frame in (spx_rv, vix_daily, spx_surface): frame['trade_date'] = pd.to_datetime(frame['trade_date']) # Forward RV targets and HAR components spx_rv = spx_rv.sort_values('trade_date') spx_rv['rv_fwd_1d'] = spx_rv['rv_5min'].shift(-1) spx_rv['rv_fwd_5d'] = spx_rv['rv_5min'].shift(-1).rolling(5, min_periods=3).sum() spx_rv['rv_fwd_22d'] = spx_rv['rv_5min'].shift(-1).rolling(22, min_periods=10).sum() spx_rv['rv_lag1'] = spx_rv['rv_5min'].shift(1) spx_rv['rv_w'] = spx_rv['rv_5min'].rolling(5, min_periods=3).mean() spx_rv['rv_m'] = spx_rv['rv_5min'].rolling(22, min_periods=10).mean() # VIX-implied daily variance: (VIX/100)² / 252 vix_daily['vix_implied_var'] = (vix_daily['vix_close'] / 100) ** 2 / 252 ml_data = spx_rv.merge(vix_daily, on='trade_date', how='inner') ml_data = ml_data.merge(spx_surface, on='trade_date', how='inner') ml_data = ml_data.sort_values('trade_date').reset_index(drop=True) ml_data = ml_data.dropna(subset=['rv_fwd_1d']) # Fill and winsorize surface features for col in SURFACE_FEATURES: if col in ml_data.columns: ml_data[col] = ml_data[col].ffill().bfill() ml_data[col] = ml_data[col].fillna(ml_data[col].median()) lo, hi = ml_data[col].quantile(0.01), ml_data[col].quantile(0.99) ml_data[col] = ml_data[col].clip(lo, hi) ml_data = ml_data.replace([np.inf, -np.inf], np.nan) ml_data = ml_data.dropna(subset=['rv_fwd_1d'] + HAR_FEATURES) print(f"\n ML dataset: {len(ml_data)} days, {ml_data.shape[1]} features") print(f" Date range: {ml_data['trade_date'].min()} to {ml_data['trade_date'].max()}") return ml_data def oos_metrics(y_true, y_pred) -> dict: mse = np.mean((y_true - y_pred) ** 2) mae = np.mean(np.abs(y_true - y_pred)) ss_res = np.sum((y_true - y_pred) ** 2) ss_tot = np.sum((y_true - y_true.mean()) ** 2) return {'mse': mse, 'mae': mae, 'r2_oos': 1 - ss_res / ss_tot} def main(): print("=" * 70) print("RQ5: ML OPTIONS SURFACE vs VIX FOR SPX RV FORECASTING") print("=" * 70) config.ensure_output_dirs() ml_data = build_dataset() all_features = HAR_FEATURES + SURFACE_FEATURES train = ml_data[ml_data['trade_date'] < SPLIT_DATE].copy() test = ml_data[ml_data['trade_date'] >= SPLIT_DATE].copy() print(f"\n Train: {len(train)} days (up to {SPLIT_DATE.date()})") print(f" Test: {len(test)} days (from {SPLIT_DATE.date()})") # ── VIX benchmark ── print("\n--- BENCHMARK: VIX ---") vix_results = {} horizon_scale = {'rv_fwd_1d': 1, 'rv_fwd_5d': 5, 'rv_fwd_22d': 22} for target_label, target_col in TARGETS.items(): sub_test = test[['vix_implied_var', target_col]].dropna() y_pred = sub_test['vix_implied_var'] * horizon_scale[target_col] res = oos_metrics(sub_test[target_col], y_pred) vix_results[target_label] = res print(f" {target_label}: MSE={res['mse']:.10f}, MAE={res['mae']:.8f}, " f"R²_oos={res['r2_oos']:.4f}") # ── OLS models ── print("\n--- OLS MODELS ---") ols_results = {} models_ols = {'HAR-RV': HAR_FEATURES, 'IV Surface': SURFACE_FEATURES, 'HAR + IV Surface': all_features} for model_name, features in models_ols.items(): for target_label, target_col in TARGETS.items(): feats = [c for c in features if c in train.columns] cols_needed = feats + [target_col] tr = train[cols_needed].replace([np.inf, -np.inf], np.nan).dropna() te = test[cols_needed].replace([np.inf, -np.inf], np.nan).dropna() if len(tr) < 50 or len(te) < 20: continue X_tr = tr[feats].values m_tr, s_tr = X_tr.mean(axis=0), X_tr.std(axis=0) s_tr[s_tr == 0] = 1 X_tr = np.column_stack([np.ones(len(X_tr)), (X_tr - m_tr) / s_tr]) coefs, _, _, _ = lstsq(X_tr, tr[target_col].values, rcond=None) X_te = np.column_stack([np.ones(len(te)), (te[feats].values - m_tr) / s_tr]) y_pred = np.maximum(X_te @ coefs, 0) res = oos_metrics(te[target_col].values, y_pred) ols_results[f"{model_name}|{target_label}"] = \ {'model': model_name, 'target': target_label, **res} print(f" {model_name} → {target_label}: MSE={res['mse']:.10f}, " f"R²_oos={res['r2_oos']:.4f}") # ── Random Forest and Gradient Boosting ── print("\n--- RANDOM FOREST / GRADIENT BOOSTING ---") from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor from sklearn.metrics import (mean_absolute_error, mean_squared_error, r2_score) from sklearn.preprocessing import StandardScaler rf_results, gbm_results = {}, {} features = [f for f in all_features if f in train.columns] for target_label, target_col in TARGETS.items(): tr = train[features + [target_col]].replace([np.inf, -np.inf], np.nan).dropna() te = test[features + [target_col]].replace([np.inf, -np.inf], np.nan).dropna() if len(tr) < 50 or len(te) < 20: continue scaler = StandardScaler() X_tr = scaler.fit_transform(tr[features].values) y_tr = tr[target_col].values X_te = scaler.transform(te[features].values) y_te = te[target_col].values rf = RandomForestRegressor(n_estimators=200, max_depth=10, min_samples_leaf=20, random_state=42, n_jobs=-1) rf.fit(X_tr, y_tr) y_pred_rf = np.maximum(rf.predict(X_te), 0) rf_results[target_label] = {'mse': mean_squared_error(y_te, y_pred_rf), 'mae': mean_absolute_error(y_te, y_pred_rf), 'r2_oos': r2_score(y_te, y_pred_rf)} print(f" RF → {target_label}: MSE={rf_results[target_label]['mse']:.10f}, " f"R²_oos={rf_results[target_label]['r2_oos']:.4f}") importances = pd.Series(rf.feature_importances_, index=features).sort_values(ascending=False) print(f" Top 5 features: {dict(importances.head(5).round(4))}") gbm = GradientBoostingRegressor(n_estimators=200, max_depth=5, learning_rate=0.05, subsample=0.8, min_samples_leaf=20, random_state=42) gbm.fit(X_tr, y_tr) y_pred_gbm = np.maximum(gbm.predict(X_te), 0) gbm_results[target_label] = {'mse': mean_squared_error(y_te, y_pred_gbm), 'mae': mean_absolute_error(y_te, y_pred_gbm), 'r2_oos': r2_score(y_te, y_pred_gbm)} print(f" GBM → {target_label}: MSE={gbm_results[target_label]['mse']:.10f}, " f"R²_oos={gbm_results[target_label]['r2_oos']:.4f}") imp_gbm = pd.Series(gbm.feature_importances_, index=features).sort_values(ascending=False) print(f" Top 5 features: {dict(imp_gbm.head(5).round(4))}") # ── Comparison table ── print("\n" + "=" * 70) print("COMPREHENSIVE MODEL COMPARISON (OUT-OF-SAMPLE)") print("=" * 70) comparison = [] for target_label in TARGETS: if target_label in vix_results: v = vix_results[target_label] comparison.append({'Model': 'VIX (benchmark)', 'Target': target_label, 'MSE': v['mse'], 'MAE': v['mae'], 'R²_OOS': v['r2_oos']}) for mn in ['HAR-RV', 'IV Surface', 'HAR + IV Surface']: key = f"{mn}|{target_label}" if key in ols_results: r = ols_results[key] comparison.append({'Model': f'OLS: {mn}', 'Target': target_label, 'MSE': r['mse'], 'MAE': r['mae'], 'R²_OOS': r['r2_oos']}) for label, res_dict in [('Random Forest', rf_results), ('Gradient Boosting', gbm_results)]: if target_label in res_dict: r = res_dict[target_label] comparison.append({'Model': label, 'Target': target_label, 'MSE': r['mse'], 'MAE': r['mae'], 'R²_OOS': r['r2_oos']}) comp_df = pd.DataFrame(comparison) comp_df.to_csv(config.RESULTS_DIR / "rq5_model_comparison.csv", index=False) for target_label in TARGETS: sub = comp_df[comp_df['Target'] == target_label].copy() sub['MSE'] = sub['MSE'].map(lambda x: f"{x:.10f}") sub['MAE'] = sub['MAE'].map(lambda x: f"{x:.8f}") sub['R²_OOS'] = sub['R²_OOS'].map(lambda x: f"{x:.4f}") print(f"\n {target_label}:") print(sub[['Model', 'MSE', 'MAE', 'R²_OOS']].to_string(index=False)) # ── Feature importance (RF, full training sample, 1-day RV) ── print("\n--- FEATURE IMPORTANCE (RF, 1-Day RV target) ---") tr_final = train[features + ['rv_fwd_1d']].replace([np.inf, -np.inf], np.nan).dropna() scaler = StandardScaler() X = scaler.fit_transform(tr_final[features].values) y = tr_final['rv_fwd_1d'].values rf_final = RandomForestRegressor(n_estimators=300, max_depth=10, min_samples_leaf=20, random_state=42, n_jobs=-1) rf_final.fit(X, y) imp = pd.DataFrame({'feature': features, 'importance': rf_final.feature_importances_}) \ .sort_values('importance', ascending=False) imp.to_csv(config.RESULTS_DIR / "rq5_feature_importance.csv", index=False) print(imp.to_string(index=False)) print("\nRQ5 COMPLETE.") if __name__ == "__main__": try: main() except RawDataUnavailableError as exc: print(f"\n[SKIPPED] {exc}", file=sys.stderr) print("The shipped results/rq5_*.csv files remain authoritative.", file=sys.stderr) sys.exit(2)