spb/wp7_uqo Public

1#!/usr/bin/env python32# =============================================================================3# Author: Simon-Pierre Boucher4# Contact: contact@spboucher.ai5# =============================================================================6"""Step 13 — Extended robustness analyses (paper revision 1.1).78Six additional checks, all computable from the processed panel alone:910A. Winsorization sensitivity — the 5-day return regression under no11   winsorization and cutoffs of 0.5%, 1% (baseline), 2.5%, 5%.12B. Newey-West lag sensitivity — HAC t-statistics at 5 (baseline), 10 and13   22 lags.14C. Rank-based information coefficients — daily cross-sectional Spearman15   correlations between each predictor and forward returns (stocks),16   with Fama-MacBeth-style time-series t-statistics.17D. Decile sorts — long-short D10−D1 portfolios for the three strongest18   sort variables, versus the baseline quintile sorts.19E. Leave-one-year-out — panel R² for the 5-day return regression and the20   HAR+IV RV model, excluding one calendar year at a time.21F. Placebo test — features permuted within ticker (fixed seed), which22   should (and does) destroy predictability.2324Inputs : data/processed/merged_options_rv.parquet25Outputs: results/extended_winsorization_sensitivity.csv,26         results/extended_newey_west_lags.csv,27         results/extended_information_coefficients.csv,28         results/extended_decile_sorts.csv,29         results/extended_leave_one_year_out.csv,30         results/extended_placebo.csv31"""3233import warnings3435import numpy as np36import pandas as pd3738import _bootstrap  # noqa: F40139from wp7 import config40from wp7.data_io import load_merged41from wp7.econometrics import (add_constant, adjusted_r_squared, hc1_tstats, ols,42                              r_squared, standardize, winsorize)43from wp7.portfolio import portfolio_sort4445warnings.filterwarnings('ignore')4647FEATURES = config.RQ1_FEATURES48RV_FEATURES = config.HAR_FEATURES + config.IV_FEATURES49KEY_VARS = ['implied_kurtosis_proxy', 'pc_volume_ratio']505152def fit_panel(sub: pd.DataFrame, features: list, target: str):53    """Standardized pooled OLS with HC1 t-stats on an already-prepared frame."""54    X = add_constant(standardize(sub[features].values))55    y = sub[target].values56    coefs = ols(X, y)57    r2 = r_squared(y, X @ coefs)58    _, t = hc1_tstats(X, y, coefs, len(features))59    return coefs, t, r2, len(sub)606162# --------------------------------------------------------------------------63# A. Winsorization sensitivity64# --------------------------------------------------------------------------65def winsorization_sensitivity(df: pd.DataFrame) -> None:66    print("\n--- A. WINSORIZATION SENSITIVITY (5-day returns) ---")67    rows = []68    for q, label in [(None, 'None'), (0.005, '0.5%'), (0.01, '1% (baseline)'),69                     (0.025, '2.5%'), (0.05, '5%')]:70        sub = df[FEATURES + ['ret_5d']].dropna()71        if q is not None:72            for f in FEATURES:73                sub[f] = winsorize(sub[f], q)74            sub['ret_5d'] = winsorize(sub['ret_5d'], q)75        _, t, r2, n = fit_panel(sub, FEATURES, 'ret_5d')76        t_by_var = dict(zip(['const'] + FEATURES, t))77        rows.append({78            'winsorization': label,79            'r2': r2,80            'adj_r2': adjusted_r_squared(r2, n, len(FEATURES)),81            'n_obs': n,82            't_implied_kurtosis': t_by_var['implied_kurtosis_proxy'],83            't_pc_volume_ratio': t_by_var['pc_volume_ratio'],84            't_implied_skewness': t_by_var['implied_skewness'],85            'n_significant_5pct': int(np.sum(np.abs(t[1:]) > 1.96)),86        })87        print(f"  {label:15s}: R²={r2:.4f}, t_kurt={t_by_var['implied_kurtosis_proxy']:6.2f}, "88              f"t_pc={t_by_var['pc_volume_ratio']:6.2f}, sig={rows[-1]['n_significant_5pct']}/10")8990    pd.DataFrame(rows).to_csv(91        config.RESULTS_DIR / "extended_winsorization_sensitivity.csv", index=False)929394# --------------------------------------------------------------------------95# B. Newey-West lag sensitivity (vectorized HAC)96# --------------------------------------------------------------------------97def nw_tstats_vectorized(X: np.ndarray, y: np.ndarray, coefs: np.ndarray,98                         n_lags: int) -> np.ndarray:99    """Bartlett-kernel HAC t-stats; einsum form of the lag sums."""100    e = y - X @ coefs101    Xe = X * e[:, None]102    S = Xe.T @ Xe103    for lag in range(1, n_lags + 1):104        w = 1 - lag / (n_lags + 1)105        cross = Xe[lag:].T @ Xe[:-lag]106        S += w * (cross + cross.T)107    try:108        XtX_inv = np.linalg.inv(X.T @ X)109    except np.linalg.LinAlgError:110        XtX_inv = np.linalg.pinv(X.T @ X)111    se = np.sqrt(np.abs(np.diag(XtX_inv @ S @ XtX_inv)))112    return coefs / np.where(se > 0, se, 1)113114115def newey_west_lags(df: pd.DataFrame) -> None:116    print("\n--- B. NEWEY-WEST LAG SENSITIVITY ---")117    rows = []118    for target, horizon in [('ret_1d', '1-Day'), ('ret_5d', '5-Day')]:119        sub = df[FEATURES + [target]].dropna()120        for f in FEATURES:121            sub[f] = winsorize(sub[f])122        sub[target] = winsorize(sub[target])123        X = add_constant(standardize(sub[FEATURES].values))124        y = sub[target].values125        coefs = ols(X, y)126        for lag in [5, 10, 22]:127            t = nw_tstats_vectorized(X, y, coefs, lag)128            for i, name in enumerate(['const'] + FEATURES):129                rows.append({'target': horizon, 'n_lags': lag, 'variable': name,130                             'coefficient': coefs[i], 'nw_t_stat': t[i],131                             'sig_5pct': abs(t[i]) > 1.96})132            n_sig = int(np.sum(np.abs(t[1:]) > 1.96))133            print(f"  {horizon} | NW({lag:2d}): {n_sig}/10 significant at 5%")134135    pd.DataFrame(rows).to_csv(136        config.RESULTS_DIR / "extended_newey_west_lags.csv", index=False)137138139# --------------------------------------------------------------------------140# C. Daily cross-sectional Spearman information coefficients (stocks)141# --------------------------------------------------------------------------142def information_coefficients(stocks: pd.DataFrame) -> None:143    print("\n--- C. RANK-BASED INFORMATION COEFFICIENTS (stocks) ---")144    rows = []145    for target, horizon in [('ret_1d', '1-Day'), ('ret_5d', '5-Day')]:146        cols = FEATURES + [target]147        daily_ics = {f: [] for f in FEATURES}148        for _, g in stocks[cols + ['trade_date']].groupby('trade_date'):149            g = g[cols].dropna()150            if len(g) < 20:151                continue152            ranks = g.rank().values153            ranks = (ranks - ranks.mean(axis=0)) / ranks.std(axis=0)154            y = ranks[:, -1]155            for j, f in enumerate(FEATURES):156                daily_ics[f].append(np.mean(ranks[:, j] * y))157158        for f in FEATURES:159            ics = np.array(daily_ics[f])160            t_stat = ics.mean() / (ics.std() / np.sqrt(len(ics)))161            rows.append({162                'target': horizon, 'variable': f,163                'mean_ic': ics.mean(), 'std_ic': ics.std(),164                't_stat': t_stat, 'pct_positive': (ics > 0).mean(),165                'n_days': len(ics),166            })167            print(f"  {horizon} | {f:25s}: IC={ics.mean():+.4f}, t={t_stat:7.2f}, "168                  f"%>0={(ics > 0).mean()*100:5.1f}")169170    pd.DataFrame(rows).to_csv(171        config.RESULTS_DIR / "extended_information_coefficients.csv", index=False)172173174# --------------------------------------------------------------------------175# D. Decile sorts176# --------------------------------------------------------------------------177def decile_sorts(stocks: pd.DataFrame) -> None:178    print("\n--- D. DECILE SORTS (D10 - D1, 5-day returns) ---")179    rows = []180    for sort_var in ['implied_kurtosis_proxy', 'pc_volume_ratio', 'iv_skew_25d']:181        for n_q, scheme in [(5, 'Quintile (baseline)'), (10, 'Decile')]:182            res = portfolio_sort(stocks, sort_var, 'ret_5d', n_quantiles=n_q)183            ls_key = f'LS_{n_q}_1'184            if res is None or ls_key not in res:185                continue186            r = res[ls_key]187            rows.append({188                'sort_variable': config.SORT_VARIABLES[sort_var],189                'scheme': scheme,190                'mean_daily_bps': r['mean_daily'] * 10000,191                'annualized_return_pct': r['annualized_return'] * 100,192                'annualized_vol_pct': r['annualized_vol'] * 100,193                'sharpe_ratio': r['sharpe'],194                't_statistic': r['t_stat'],195                'n_days': r['n_days'],196            })197            print(f"  {config.SORT_VARIABLES[sort_var]:25s} | {scheme:20s}: "198                  f"ann.ret={r['annualized_return']*100:7.2f}%, "199                  f"Sharpe={r['sharpe']:6.3f}, t={r['t_stat']:7.2f}")200201    pd.DataFrame(rows).to_csv(202        config.RESULTS_DIR / "extended_decile_sorts.csv", index=False)203204205# --------------------------------------------------------------------------206# E. Leave-one-year-out207# --------------------------------------------------------------------------208def leave_one_year_out(df: pd.DataFrame) -> None:209    print("\n--- E. LEAVE-ONE-YEAR-OUT PANEL R² ---")210    years = sorted(df['trade_date'].dt.year.unique())211    rows = []212    for year in years:213        sub_df = df[df['trade_date'].dt.year != year]214215        sub = sub_df[FEATURES + ['ret_5d']].dropna()216        for f in FEATURES:217            sub[f] = winsorize(sub[f])218        sub['ret_5d'] = winsorize(sub['ret_5d'])219        _, _, r2_ret, n_ret = fit_panel(sub, FEATURES, 'ret_5d')220221        sub_rv = sub_df[RV_FEATURES + ['rv_fwd_1d']].dropna()222        for f in RV_FEATURES:223            sub_rv[f] = winsorize(sub_rv[f])224        sub_rv['rv_fwd_1d'] = winsorize(sub_rv['rv_fwd_1d'])225        _, _, r2_rv, _ = fit_panel(sub_rv, RV_FEATURES, 'rv_fwd_1d')226227        rows.append({'excluded_year': year, 'r2_ret_5d': r2_ret,228                     'r2_rv_har_iv': r2_rv, 'n_obs_ret': n_ret})229        print(f"  excl. {year}: 5D-ret R²={r2_ret:.4f}, HAR+IV RV R²={r2_rv:.4f}")230231    pd.DataFrame(rows).to_csv(232        config.RESULTS_DIR / "extended_leave_one_year_out.csv", index=False)233234235# --------------------------------------------------------------------------236# F. Placebo: within-ticker permutation of the feature block237# --------------------------------------------------------------------------238def placebo_test(df: pd.DataFrame, n_draws: int = 10) -> None:239    print("\n--- F. PLACEBO TEST (features permuted within ticker) ---")240    base = df[FEATURES + ['ret_5d', 'ticker']].dropna().reset_index(drop=True)241    for f in FEATURES:242        base[f] = winsorize(base[f])243    base['ret_5d'] = winsorize(base['ret_5d'])244245    _, _, r2_actual, n = fit_panel(base, FEATURES, 'ret_5d')246    print(f"  Actual R²: {r2_actual:.4f} (N={n:,})")247248    rng = np.random.default_rng(42)249    ticker_codes = base['ticker'].astype('category').cat.codes.values250    feat_matrix = base[FEATURES].values251    r2_placebos = []252    for draw in range(n_draws):253        perm_matrix = np.empty_like(feat_matrix)254        for code in np.unique(ticker_codes):255            idx = np.flatnonzero(ticker_codes == code)256            perm_matrix[idx] = feat_matrix[idx[rng.permutation(len(idx))]]257        X = add_constant(standardize(perm_matrix))258        y = base['ret_5d'].values259        coefs = ols(X, y)260        r2_placebos.append(r_squared(y, X @ coefs))261262    r2_placebos = np.array(r2_placebos)263    print(f"  Placebo R² over {n_draws} draws: mean={r2_placebos.mean():.6f}, "264          f"max={r2_placebos.max():.6f}")265266    pd.DataFrame({267        'draw': ['actual'] + [str(i + 1) for i in range(n_draws)],268        'r2': [r2_actual] + list(r2_placebos),269    }).to_csv(config.RESULTS_DIR / "extended_placebo.csv", index=False)270271272def main():273    print("=" * 70)274    print("EXTENDED ROBUSTNESS ANALYSES")275    print("=" * 70)276    config.ensure_output_dirs()277278    df = load_merged()279    stocks = df[~df['ticker'].isin(config.NON_STOCK_TICKERS)].copy()280281    winsorization_sensitivity(df)282    newey_west_lags(df)283    information_coefficients(stocks)284    decile_sorts(stocks)285    leave_one_year_out(df)286    placebo_test(df)287288    print("\nEXTENDED ROBUSTNESS COMPLETE.")289290291if __name__ == "__main__":292    main()293