spb/wp7_uqo Public

1#!/usr/bin/env python32# =============================================================================3# Author: Simon-Pierre Boucher4# Contact: contact@spboucher.ai5# =============================================================================6"""Step 10 — Robustness checks.78A. Newey-West HAC(5) standard errors9B. Double-clustered (ticker + date) standard errors, CGM (2011)10C. Additional control variables (volume, OI, lagged returns)11D. Quantile regressions (IRLS) at τ ∈ {0.10, 0.25, 0.50, 0.75, 0.90}12E. Excluding high-VIX periods (needs the raw VIX series; console-only)13F. Ticker-by-ticker R² distribution1415Note: the original implementation of section D materialized an n×n diagonal16weight matrix (~97 GB at this sample size) and could not complete; the IRLS17weighting is now applied by broadcasting — mathematically identical — so18sections D and F produce the two result files that were missing from the19original archive (see AUDIT.md §6.2).2021Inputs : data/processed/merged_options_rv.parquet22         (+ index_5min.duckdb for section E)23Outputs: results/robustness_newey_west.csv,24         results/robustness_double_clustered.csv,25         results/robustness_with_controls.csv,26         results/robustness_quantile_regression.csv,27         results/robustness_ticker_r2.csv28"""2930import warnings3132import numpy as np33import pandas as pd3435import _bootstrap  # noqa: F40136from wp7 import config37from wp7.data_io import RawDataUnavailableError, load_merged, load_vix_daily38from wp7.econometrics import (add_constant, double_clustered_tstats, hc1_tstats,39                              newey_west_tstats, ols, quantile_regression,40                              r_squared, standardize, winsorize)4142warnings.filterwarnings('ignore')4344FEATURES = config.RQ1_FEATURES45HORIZONS = [('ret_1d', '1-Day'), ('ret_5d', '5-Day')]464748def prepare(data: pd.DataFrame, features: list, target: str,49            keep: tuple = ()) -> pd.DataFrame:50    """Drop incomplete rows and winsorize features and target at 1%/99%."""51    sub = data[features + [target, *keep]].dropna()52    for f in features:53        sub[f] = winsorize(sub[f])54    sub[target] = winsorize(sub[target])55    return sub565758def main():59    print("=" * 70)60    print("ROBUSTNESS CHECKS")61    print("=" * 70)62    config.ensure_output_dirs()6364    df = load_merged()6566    # ── A. Newey-West HAC standard errors ──67    print("\n--- A. NEWEY-WEST HAC STANDARD ERRORS (lag=5) ---")68    nw_results = []69    for target, horizon in HORIZONS:70        sub = prepare(df, FEATURES, target)71        X = add_constant(standardize(sub[FEATURES].values))72        y = sub[target].values73        coefs = ols(X, y)74        se, t_stats = newey_west_tstats(X, y, coefs, n_lags=5)75        r2 = r_squared(y, X @ coefs)7677        res_df = pd.DataFrame({78            'variable': ['const'] + FEATURES,79            'coefficient': coefs,80            'nw_se': se,81            'nw_t_stat': t_stats,82            'nw_sig_5pct': np.abs(t_stats) > 1.96,83        })84        res_df['target'] = horizon85        res_df['method'] = 'Newey-West(5)'86        nw_results.append(res_df)8788        sig = res_df[(res_df['nw_sig_5pct']) & (res_df['variable'] != 'const')]89        print(f"\n  {horizon}: R²={r2:.6f}, N={len(sub):,}")90        print(f"  Significant (NW): {', '.join(sig['variable'].tolist())}")91        for _, row in res_df.iterrows():92            star = "**" if abs(row['nw_t_stat']) > 2.576 else \93                ("*" if abs(row['nw_t_stat']) > 1.96 else "")94            print(f"    {row['variable']:25s}: β={row['coefficient']:9.6f}  "95                  f"t_NW={row['nw_t_stat']:7.3f} {star}")9697    pd.concat(nw_results, ignore_index=True).to_csv(98        config.RESULTS_DIR / "robustness_newey_west.csv", index=False)99100    # ── B. Double-clustered standard errors (ticker + date) ──101    print("\n--- B. DOUBLE-CLUSTERED SE (ticker + date) ---")102    dc_results = []103    for target, horizon in HORIZONS:104        sub = prepare(df, FEATURES, target, keep=('ticker', 'trade_date'))105        X = add_constant(standardize(sub[FEATURES].values))106        y = sub[target].values107        coefs = ols(X, y)108        se, t_stats = double_clustered_tstats(109            X, y, coefs, sub['ticker'].values, sub['trade_date'].values)110        r2 = r_squared(y, X @ coefs)111112        res_df = pd.DataFrame({113            'variable': ['const'] + FEATURES,114            'coefficient': coefs,115            'dc_se': se,116            'dc_t_stat': t_stats,117            'dc_sig_5pct': np.abs(t_stats) > 1.96,118        })119        res_df['target'] = horizon120        dc_results.append(res_df)121122        sig = res_df[(res_df['dc_sig_5pct']) & (res_df['variable'] != 'const')]123        print(f"\n  {horizon}: R²={r2:.6f}, N={len(sub):,}")124        print(f"  Significant (DC): {', '.join(sig['variable'].tolist())}")125        for _, row in res_df.iterrows():126            star = "**" if abs(row['dc_t_stat']) > 2.576 else \127                ("*" if abs(row['dc_t_stat']) > 1.96 else "")128            print(f"    {row['variable']:25s}: β={row['coefficient']:9.6f}  "129                  f"t_DC={row['dc_t_stat']:7.3f} {star}")130131    pd.concat(dc_results, ignore_index=True).to_csv(132        config.RESULTS_DIR / "robustness_double_clustered.csv", index=False)133134    # ── C. Additional control variables ──135    print("\n--- C. WITH CONTROL VARIABLES (volume, log_oi, spread proxy) ---")136    df['log_volume'] = np.log1p(df['total_option_volume'])137    df['log_oi'] = np.log1p(df['total_oi'])138    df['abs_return'] = df['daily_return'].abs()139    df['ret_lag1'] = df.groupby('ticker')['daily_return'].shift(1)140    df['ret_lag5'] = df.groupby('ticker')['daily_return'].transform(141        lambda x: x.shift(1).rolling(5).sum())142143    controls = ['log_volume', 'log_oi', 'abs_return', 'ret_lag1', 'ret_lag5']144    features_ctrl = FEATURES + controls145146    ctrl_results = []147    for target, horizon in HORIZONS:148        sub = prepare(df, features_ctrl, target)149        X = add_constant(standardize(sub[features_ctrl].values))150        y = sub[target].values151        coefs = ols(X, y)152        r2 = r_squared(y, X @ coefs)153        _, t = hc1_tstats(X, y, coefs, len(features_ctrl))154155        print(f"\n  {horizon} with controls: R²={r2:.6f}, N={len(sub):,}")156        for i, name in enumerate(['const'] + features_ctrl):157            star = "**" if abs(t[i]) > 2.576 else ("*" if abs(t[i]) > 1.96 else "")158            print(f"    {name:25s}: β={coefs[i]:9.6f}  t={t[i]:7.3f} {star}")159            ctrl_results.append({'target': horizon, 'variable': name,160                                 'coefficient': coefs[i], 't_stat': t[i], 'r2': r2})161162    pd.DataFrame(ctrl_results).to_csv(163        config.RESULTS_DIR / "robustness_with_controls.csv", index=False)164165    # ── D. Quantile regressions ──166    print("\n--- D. QUANTILE REGRESSION (iterative reweighting) ---")167    qr_results = []168    for target, horizon in [('ret_5d', '5-Day')]:169        sub = prepare(df, FEATURES, target)170        X = add_constant(standardize(sub[FEATURES].values))171        y = sub[target].values172173        for tau in [0.10, 0.25, 0.50, 0.75, 0.90]:174            coefs = quantile_regression(X, y, tau)175            print(f"\n  {horizon} | τ={tau}:")176            for i, name in enumerate(['const'] + FEATURES):177                qr_results.append({'target': horizon, 'tau': tau,178                                   'variable': name, 'coefficient': coefs[i]})179                print(f"    {name:25s}: β={coefs[i]:9.6f}")180181    pd.DataFrame(qr_results).to_csv(182        config.RESULTS_DIR / "robustness_quantile_regression.csv", index=False)183184    # ── E. Excluding high-VIX periods (console-only, needs raw VIX) ──185    print("\n--- E. EXCLUDING HIGH-VIX PERIODS (VIX < 30) ---")186    try:187        vix = load_vix_daily()188        calm = df.merge(vix, on='trade_date', how='left')189        calm = calm[calm['vix_close'] < 30]190        for target, horizon in HORIZONS:191            sub = prepare(calm, FEATURES, target)192            X = add_constant(standardize(sub[FEATURES].values))193            y = sub[target].values194            coefs = ols(X, y)195            print(f"  {horizon} (VIX<30): R²={r_squared(y, X @ coefs):.6f}, N={len(sub):,}")196    except RawDataUnavailableError as exc:197        print(f"  [Section E skipped — raw stores unavailable]\n  {exc}")198199    # ── F. Ticker-by-ticker R² distribution (5-day returns) ──200    print("\n--- F. TICKER-BY-TICKER R² DISTRIBUTION (5-Day returns) ---")201    ticker_r2 = []202    for ticker in df['ticker'].unique():203        td = df[df['ticker'] == ticker]204        sub = td[FEATURES + ['ret_5d']].dropna()205        if len(sub) < 200:206            continue207        for f in FEATURES:208            sub[f] = winsorize(sub[f])209        sub['ret_5d'] = winsorize(sub['ret_5d'])210        X = add_constant(standardize(sub[FEATURES].values))211        y = sub['ret_5d'].values212        coefs = ols(X, y)213        ticker_r2.append({'ticker': ticker, 'r2_5d': r_squared(y, X @ coefs),214                          'n_obs': len(sub)})215216    ticker_r2_df = pd.DataFrame(ticker_r2)217    ticker_r2_df.to_csv(config.RESULTS_DIR / "robustness_ticker_r2.csv", index=False)218    print(f"  Distribution of R² across {len(ticker_r2_df)} tickers:")219    for stat, val in [('Mean', ticker_r2_df['r2_5d'].mean()),220                      ('Median', ticker_r2_df['r2_5d'].median()),221                      ('Std', ticker_r2_df['r2_5d'].std()),222                      ('Min', ticker_r2_df['r2_5d'].min()),223                      ('Max', ticker_r2_df['r2_5d'].max()),224                      ('Q25', ticker_r2_df['r2_5d'].quantile(0.25)),225                      ('Q75', ticker_r2_df['r2_5d'].quantile(0.75))]:226        print(f"    {stat}:   {val:.6f}")227228    print("\nROBUSTNESS CHECKS COMPLETE.")229230231if __name__ == "__main__":232    main()233