spb/wp7_uqo Public
UQO Working Paper No. 7 — Options-implied information for cross-asset return and volatility prediction: evidence from 3.8B option contracts.
Python 66.5%
TeX 32.7%
Makefile 0.8%
1#!/usr/bin/env python32# =============================================================================3# Author: Simon-Pierre Boucher4# Contact: contact@spboucher.ai5# =============================================================================6"""Step 06 — RQ5: ML on the full SPX options surface vs VIX for realized-7variance forecasting.89Builds 21 daily SPX surface features and SPX 5-minute realized variance from10the raw stores, then compares — out-of-sample, temporal split 2020-01-01 —11the VIX benchmark, OLS models (HAR-RV, IV surface, combined), Random Forest12and Gradient Boosting. Entirely raw-dependent: exits gracefully when the raw13stores are absent (the shipped result CSVs remain authoritative).1415Inputs : options.duckdb, index_5min.duckdb16Outputs: results/rq5_model_comparison.csv, results/rq5_feature_importance.csv17"""1819import sys20import warnings2122import numpy as np23import pandas as pd24from numpy.linalg import lstsq2526import _bootstrap # noqa: F40127from wp7 import config28from wp7.data_io import RawDataUnavailableError, open_raw_db2930warnings.filterwarnings('ignore')3132SURFACE_FEATURES = [33 'iv_atm_1w', 'iv_atm_2w', 'iv_atm_1m', 'iv_atm_2m', 'iv_atm_3m', 'iv_atm_6m',34 'skew_25d_1m', 'skew_25d_3m', 'skew_10d_1m', 'butterfly_1m',35 'ts_slope_3m_1m', 'ts_slope_6m_1m',36 'total_gamma_oi', 'net_gamma', 'total_vega_oi', 'avg_theta',37 'pc_vol_ratio', 'pc_oi_ratio', 'total_volume', 'total_oi', 'avg_spread_pct',38]39HAR_FEATURES = ['rv_lag1', 'rv_w', 'rv_m']40TARGETS = {'1-Day RV': 'rv_fwd_1d', '5-Day RV': 'rv_fwd_5d', '22-Day RV': 'rv_fwd_22d'}41SPLIT_DATE = pd.Timestamp('2020-01-01')424344def build_dataset() -> pd.DataFrame:45 """SPX surface features + SPX 5-min RV + VIX-implied daily variance."""46 opt_con = open_raw_db("options")47 print("\n--- BUILDING SPX OPTIONS SURFACE FEATURES ---")48 spx_surface = opt_con.execute("""49 WITH base AS (50 SELECT trade_date, strike, expiry_date, call_put,51 (expiry_date - trade_date) AS dte,52 bid_price, ask_price, bid_iv, ask_iv,53 open_interest, volume, delta, gamma, vega, theta, rho,54 CASE WHEN bid_iv > 0 AND ask_iv > 0 THEN (bid_iv + ask_iv)/2.055 WHEN ask_iv > 0 THEN ask_iv ELSE bid_iv END AS mid_iv,56 (bid_price + ask_price) / 2.0 AS mid_price57 FROM option_chain58 WHERE ticker = 'SPX'59 AND ask_price > 060 AND (expiry_date - trade_date) BETWEEN 1 AND 36561 )62 SELECT63 trade_date,6465 -- ATM IV by tenor66 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,67 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,68 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,69 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,70 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,71 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,7273 -- Skew by tenor (25-delta put minus 25-delta call)74 AVG(CASE WHEN call_put='p' AND ABS(delta+0.25)<0.07 AND dte BETWEEN 25 AND 35 THEN mid_iv END) -75 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,76 AVG(CASE WHEN call_put='p' AND ABS(delta+0.25)<0.07 AND dte BETWEEN 85 AND 95 THEN mid_iv END) -77 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,7879 -- Deep OTM skew (10-delta)80 AVG(CASE WHEN call_put='p' AND ABS(delta+0.10)<0.05 AND dte BETWEEN 25 AND 35 THEN mid_iv END) -81 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,8283 -- Butterfly (wings / ATM)84 (AVG(CASE WHEN ABS(delta)<0.15 AND ABS(delta)>0.03 AND dte BETWEEN 25 AND 35 THEN mid_iv END) /85 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)86 ) AS butterfly_1m,8788 -- Term structure slopes89 AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 85 AND 95 THEN mid_iv END) -90 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,91 AVG(CASE WHEN call_put='c' AND ABS(delta-0.5)<0.1 AND dte BETWEEN 170 AND 200 THEN mid_iv END) -92 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,9394 -- Aggregate Greeks95 SUM(gamma * open_interest) AS total_gamma_oi,96 SUM(CASE WHEN call_put='c' THEN gamma * open_interest ELSE 0 END) -97 SUM(CASE WHEN call_put='p' THEN gamma * open_interest ELSE 0 END) AS net_gamma,98 SUM(vega * open_interest) AS total_vega_oi,99 AVG(theta) AS avg_theta,100101 -- Volume/OI ratios102 SUM(CASE WHEN call_put='p' THEN volume ELSE 0 END)::DOUBLE /103 NULLIF(SUM(CASE WHEN call_put='c' THEN volume ELSE 0 END),0) AS pc_vol_ratio,104 SUM(CASE WHEN call_put='p' THEN open_interest ELSE 0 END)::DOUBLE /105 NULLIF(SUM(CASE WHEN call_put='c' THEN open_interest ELSE 0 END),0) AS pc_oi_ratio,106 SUM(volume) AS total_volume,107 SUM(open_interest) AS total_oi,108109 -- Liquidity proxy110 AVG(CASE WHEN mid_price > 0 THEN (ask_price - bid_price) / mid_price END) AS avg_spread_pct,111112 COUNT(*) AS n_contracts113114 FROM base115 GROUP BY trade_date116 HAVING COUNT(*) >= 50117 ORDER BY trade_date118 """).fetchdf()119 opt_con.close()120 print(f" SPX surface features: {len(spx_surface)} days, {spx_surface.shape[1]} columns")121122 idx_con = open_raw_db("indices_5min")123 spx_rv = idx_con.execute("""124 WITH bars AS (125 SELECT datetime, close,126 CAST(datetime AS DATE) AS trade_date,127 LAG(close) OVER (ORDER BY datetime) AS prev_close128 FROM ohlcv WHERE symbol = 'SPX'129 ),130 rets AS (131 SELECT trade_date, LN(close / NULLIF(prev_close, 0)) AS log_ret132 FROM bars WHERE prev_close > 0 AND close > 0133 )134 SELECT trade_date,135 SUM(log_ret * log_ret) AS rv_5min,136 COUNT(*) AS n_obs,137 SUM(log_ret) AS daily_ret138 FROM rets139 GROUP BY trade_date140 HAVING COUNT(*) >= 20141 ORDER BY trade_date142 """).fetchdf()143 vix_daily = idx_con.execute("""144 SELECT CAST(datetime AS DATE) AS trade_date,145 LAST(close) AS vix_close146 FROM ohlcv WHERE symbol = 'VIX'147 GROUP BY CAST(datetime AS DATE)148 ORDER BY trade_date149 """).fetchdf()150 idx_con.close()151152 for frame in (spx_rv, vix_daily, spx_surface):153 frame['trade_date'] = pd.to_datetime(frame['trade_date'])154155 # Forward RV targets and HAR components156 spx_rv = spx_rv.sort_values('trade_date')157 spx_rv['rv_fwd_1d'] = spx_rv['rv_5min'].shift(-1)158 spx_rv['rv_fwd_5d'] = spx_rv['rv_5min'].shift(-1).rolling(5, min_periods=3).sum()159 spx_rv['rv_fwd_22d'] = spx_rv['rv_5min'].shift(-1).rolling(22, min_periods=10).sum()160 spx_rv['rv_lag1'] = spx_rv['rv_5min'].shift(1)161 spx_rv['rv_w'] = spx_rv['rv_5min'].rolling(5, min_periods=3).mean()162 spx_rv['rv_m'] = spx_rv['rv_5min'].rolling(22, min_periods=10).mean()163164 # VIX-implied daily variance: (VIX/100)² / 252165 vix_daily['vix_implied_var'] = (vix_daily['vix_close'] / 100) ** 2 / 252166167 ml_data = spx_rv.merge(vix_daily, on='trade_date', how='inner')168 ml_data = ml_data.merge(spx_surface, on='trade_date', how='inner')169 ml_data = ml_data.sort_values('trade_date').reset_index(drop=True)170 ml_data = ml_data.dropna(subset=['rv_fwd_1d'])171172 # Fill and winsorize surface features173 for col in SURFACE_FEATURES:174 if col in ml_data.columns:175 ml_data[col] = ml_data[col].ffill().bfill()176 ml_data[col] = ml_data[col].fillna(ml_data[col].median())177 lo, hi = ml_data[col].quantile(0.01), ml_data[col].quantile(0.99)178 ml_data[col] = ml_data[col].clip(lo, hi)179180 ml_data = ml_data.replace([np.inf, -np.inf], np.nan)181 ml_data = ml_data.dropna(subset=['rv_fwd_1d'] + HAR_FEATURES)182183 print(f"\n ML dataset: {len(ml_data)} days, {ml_data.shape[1]} features")184 print(f" Date range: {ml_data['trade_date'].min()} to {ml_data['trade_date'].max()}")185 return ml_data186187188def oos_metrics(y_true, y_pred) -> dict:189 mse = np.mean((y_true - y_pred) ** 2)190 mae = np.mean(np.abs(y_true - y_pred))191 ss_res = np.sum((y_true - y_pred) ** 2)192 ss_tot = np.sum((y_true - y_true.mean()) ** 2)193 return {'mse': mse, 'mae': mae, 'r2_oos': 1 - ss_res / ss_tot}194195196def main():197 print("=" * 70)198 print("RQ5: ML OPTIONS SURFACE vs VIX FOR SPX RV FORECASTING")199 print("=" * 70)200 config.ensure_output_dirs()201202 ml_data = build_dataset()203 all_features = HAR_FEATURES + SURFACE_FEATURES204205 train = ml_data[ml_data['trade_date'] < SPLIT_DATE].copy()206 test = ml_data[ml_data['trade_date'] >= SPLIT_DATE].copy()207 print(f"\n Train: {len(train)} days (up to {SPLIT_DATE.date()})")208 print(f" Test: {len(test)} days (from {SPLIT_DATE.date()})")209210 # ── VIX benchmark ──211 print("\n--- BENCHMARK: VIX ---")212 vix_results = {}213 horizon_scale = {'rv_fwd_1d': 1, 'rv_fwd_5d': 5, 'rv_fwd_22d': 22}214 for target_label, target_col in TARGETS.items():215 sub_test = test[['vix_implied_var', target_col]].dropna()216 y_pred = sub_test['vix_implied_var'] * horizon_scale[target_col]217 res = oos_metrics(sub_test[target_col], y_pred)218 vix_results[target_label] = res219 print(f" {target_label}: MSE={res['mse']:.10f}, MAE={res['mae']:.8f}, "220 f"R²_oos={res['r2_oos']:.4f}")221222 # ── OLS models ──223 print("\n--- OLS MODELS ---")224 ols_results = {}225 models_ols = {'HAR-RV': HAR_FEATURES, 'IV Surface': SURFACE_FEATURES,226 'HAR + IV Surface': all_features}227 for model_name, features in models_ols.items():228 for target_label, target_col in TARGETS.items():229 feats = [c for c in features if c in train.columns]230 cols_needed = feats + [target_col]231 tr = train[cols_needed].replace([np.inf, -np.inf], np.nan).dropna()232 te = test[cols_needed].replace([np.inf, -np.inf], np.nan).dropna()233 if len(tr) < 50 or len(te) < 20:234 continue235236 X_tr = tr[feats].values237 m_tr, s_tr = X_tr.mean(axis=0), X_tr.std(axis=0)238 s_tr[s_tr == 0] = 1239 X_tr = np.column_stack([np.ones(len(X_tr)), (X_tr - m_tr) / s_tr])240 coefs, _, _, _ = lstsq(X_tr, tr[target_col].values, rcond=None)241242 X_te = np.column_stack([np.ones(len(te)),243 (te[feats].values - m_tr) / s_tr])244 y_pred = np.maximum(X_te @ coefs, 0)245 res = oos_metrics(te[target_col].values, y_pred)246 ols_results[f"{model_name}|{target_label}"] = \247 {'model': model_name, 'target': target_label, **res}248 print(f" {model_name} → {target_label}: MSE={res['mse']:.10f}, "249 f"R²_oos={res['r2_oos']:.4f}")250251 # ── Random Forest and Gradient Boosting ──252 print("\n--- RANDOM FOREST / GRADIENT BOOSTING ---")253 from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor254 from sklearn.metrics import (mean_absolute_error, mean_squared_error, r2_score)255 from sklearn.preprocessing import StandardScaler256257 rf_results, gbm_results = {}, {}258 features = [f for f in all_features if f in train.columns]259 for target_label, target_col in TARGETS.items():260 tr = train[features + [target_col]].replace([np.inf, -np.inf], np.nan).dropna()261 te = test[features + [target_col]].replace([np.inf, -np.inf], np.nan).dropna()262 if len(tr) < 50 or len(te) < 20:263 continue264265 scaler = StandardScaler()266 X_tr = scaler.fit_transform(tr[features].values)267 y_tr = tr[target_col].values268 X_te = scaler.transform(te[features].values)269 y_te = te[target_col].values270271 rf = RandomForestRegressor(n_estimators=200, max_depth=10,272 min_samples_leaf=20, random_state=42, n_jobs=-1)273 rf.fit(X_tr, y_tr)274 y_pred_rf = np.maximum(rf.predict(X_te), 0)275 rf_results[target_label] = {'mse': mean_squared_error(y_te, y_pred_rf),276 'mae': mean_absolute_error(y_te, y_pred_rf),277 'r2_oos': r2_score(y_te, y_pred_rf)}278 print(f" RF → {target_label}: MSE={rf_results[target_label]['mse']:.10f}, "279 f"R²_oos={rf_results[target_label]['r2_oos']:.4f}")280 importances = pd.Series(rf.feature_importances_, index=features).sort_values(ascending=False)281 print(f" Top 5 features: {dict(importances.head(5).round(4))}")282283 gbm = GradientBoostingRegressor(n_estimators=200, max_depth=5,284 learning_rate=0.05, subsample=0.8,285 min_samples_leaf=20, random_state=42)286 gbm.fit(X_tr, y_tr)287 y_pred_gbm = np.maximum(gbm.predict(X_te), 0)288 gbm_results[target_label] = {'mse': mean_squared_error(y_te, y_pred_gbm),289 'mae': mean_absolute_error(y_te, y_pred_gbm),290 'r2_oos': r2_score(y_te, y_pred_gbm)}291 print(f" GBM → {target_label}: MSE={gbm_results[target_label]['mse']:.10f}, "292 f"R²_oos={gbm_results[target_label]['r2_oos']:.4f}")293 imp_gbm = pd.Series(gbm.feature_importances_, index=features).sort_values(ascending=False)294 print(f" Top 5 features: {dict(imp_gbm.head(5).round(4))}")295296 # ── Comparison table ──297 print("\n" + "=" * 70)298 print("COMPREHENSIVE MODEL COMPARISON (OUT-OF-SAMPLE)")299 print("=" * 70)300 comparison = []301 for target_label in TARGETS:302 if target_label in vix_results:303 v = vix_results[target_label]304 comparison.append({'Model': 'VIX (benchmark)', 'Target': target_label,305 'MSE': v['mse'], 'MAE': v['mae'], 'R²_OOS': v['r2_oos']})306 for mn in ['HAR-RV', 'IV Surface', 'HAR + IV Surface']:307 key = f"{mn}|{target_label}"308 if key in ols_results:309 r = ols_results[key]310 comparison.append({'Model': f'OLS: {mn}', 'Target': target_label,311 'MSE': r['mse'], 'MAE': r['mae'], 'R²_OOS': r['r2_oos']})312 for label, res_dict in [('Random Forest', rf_results),313 ('Gradient Boosting', gbm_results)]:314 if target_label in res_dict:315 r = res_dict[target_label]316 comparison.append({'Model': label, 'Target': target_label,317 'MSE': r['mse'], 'MAE': r['mae'], 'R²_OOS': r['r2_oos']})318319 comp_df = pd.DataFrame(comparison)320 comp_df.to_csv(config.RESULTS_DIR / "rq5_model_comparison.csv", index=False)321322 for target_label in TARGETS:323 sub = comp_df[comp_df['Target'] == target_label].copy()324 sub['MSE'] = sub['MSE'].map(lambda x: f"{x:.10f}")325 sub['MAE'] = sub['MAE'].map(lambda x: f"{x:.8f}")326 sub['R²_OOS'] = sub['R²_OOS'].map(lambda x: f"{x:.4f}")327 print(f"\n {target_label}:")328 print(sub[['Model', 'MSE', 'MAE', 'R²_OOS']].to_string(index=False))329330 # ── Feature importance (RF, full training sample, 1-day RV) ──331 print("\n--- FEATURE IMPORTANCE (RF, 1-Day RV target) ---")332 tr_final = train[features + ['rv_fwd_1d']].replace([np.inf, -np.inf], np.nan).dropna()333 scaler = StandardScaler()334 X = scaler.fit_transform(tr_final[features].values)335 y = tr_final['rv_fwd_1d'].values336 rf_final = RandomForestRegressor(n_estimators=300, max_depth=10,337 min_samples_leaf=20, random_state=42, n_jobs=-1)338 rf_final.fit(X, y)339 imp = pd.DataFrame({'feature': features,340 'importance': rf_final.feature_importances_}) \341 .sort_values('importance', ascending=False)342 imp.to_csv(config.RESULTS_DIR / "rq5_feature_importance.csv", index=False)343 print(imp.to_string(index=False))344345 print("\nRQ5 COMPLETE.")346347348if __name__ == "__main__":349 try:350 main()351 except RawDataUnavailableError as exc:352 print(f"\n[SKIPPED] {exc}", file=sys.stderr)353 print("The shipped results/rq5_*.csv files remain authoritative.", file=sys.stderr)354 sys.exit(2)355