#!/usr/bin/env python3 # ============================================================================= # Author: Simon-Pierre Boucher # Contact: contact@spboucher.ai # ============================================================================= """Step 13 — Extended robustness analyses (paper revision 1.1). Six additional checks, all computable from the processed panel alone: A. Winsorization sensitivity — the 5-day return regression under no winsorization and cutoffs of 0.5%, 1% (baseline), 2.5%, 5%. B. Newey-West lag sensitivity — HAC t-statistics at 5 (baseline), 10 and 22 lags. C. Rank-based information coefficients — daily cross-sectional Spearman correlations between each predictor and forward returns (stocks), with Fama-MacBeth-style time-series t-statistics. D. Decile sorts — long-short D10−D1 portfolios for the three strongest sort variables, versus the baseline quintile sorts. E. Leave-one-year-out — panel R² for the 5-day return regression and the HAR+IV RV model, excluding one calendar year at a time. F. Placebo test — features permuted within ticker (fixed seed), which should (and does) destroy predictability. Inputs : data/processed/merged_options_rv.parquet Outputs: results/extended_winsorization_sensitivity.csv, results/extended_newey_west_lags.csv, results/extended_information_coefficients.csv, results/extended_decile_sorts.csv, results/extended_leave_one_year_out.csv, results/extended_placebo.csv """ import warnings import numpy as np import pandas as pd import _bootstrap # noqa: F401 from wp7 import config from wp7.data_io import load_merged from wp7.econometrics import (add_constant, adjusted_r_squared, hc1_tstats, ols, r_squared, standardize, winsorize) from wp7.portfolio import portfolio_sort warnings.filterwarnings('ignore') FEATURES = config.RQ1_FEATURES RV_FEATURES = config.HAR_FEATURES + config.IV_FEATURES KEY_VARS = ['implied_kurtosis_proxy', 'pc_volume_ratio'] def fit_panel(sub: pd.DataFrame, features: list, target: str): """Standardized pooled OLS with HC1 t-stats on an already-prepared frame.""" X = add_constant(standardize(sub[features].values)) y = sub[target].values coefs = ols(X, y) r2 = r_squared(y, X @ coefs) _, t = hc1_tstats(X, y, coefs, len(features)) return coefs, t, r2, len(sub) # -------------------------------------------------------------------------- # A. Winsorization sensitivity # -------------------------------------------------------------------------- def winsorization_sensitivity(df: pd.DataFrame) -> None: print("\n--- A. WINSORIZATION SENSITIVITY (5-day returns) ---") rows = [] for q, label in [(None, 'None'), (0.005, '0.5%'), (0.01, '1% (baseline)'), (0.025, '2.5%'), (0.05, '5%')]: sub = df[FEATURES + ['ret_5d']].dropna() if q is not None: for f in FEATURES: sub[f] = winsorize(sub[f], q) sub['ret_5d'] = winsorize(sub['ret_5d'], q) _, t, r2, n = fit_panel(sub, FEATURES, 'ret_5d') t_by_var = dict(zip(['const'] + FEATURES, t)) rows.append({ 'winsorization': label, 'r2': r2, 'adj_r2': adjusted_r_squared(r2, n, len(FEATURES)), 'n_obs': n, 't_implied_kurtosis': t_by_var['implied_kurtosis_proxy'], 't_pc_volume_ratio': t_by_var['pc_volume_ratio'], 't_implied_skewness': t_by_var['implied_skewness'], 'n_significant_5pct': int(np.sum(np.abs(t[1:]) > 1.96)), }) print(f" {label:15s}: R²={r2:.4f}, t_kurt={t_by_var['implied_kurtosis_proxy']:6.2f}, " f"t_pc={t_by_var['pc_volume_ratio']:6.2f}, sig={rows[-1]['n_significant_5pct']}/10") pd.DataFrame(rows).to_csv( config.RESULTS_DIR / "extended_winsorization_sensitivity.csv", index=False) # -------------------------------------------------------------------------- # B. Newey-West lag sensitivity (vectorized HAC) # -------------------------------------------------------------------------- def nw_tstats_vectorized(X: np.ndarray, y: np.ndarray, coefs: np.ndarray, n_lags: int) -> np.ndarray: """Bartlett-kernel HAC t-stats; einsum form of the lag sums.""" e = y - X @ coefs Xe = X * e[:, None] S = Xe.T @ Xe for lag in range(1, n_lags + 1): w = 1 - lag / (n_lags + 1) cross = Xe[lag:].T @ Xe[:-lag] S += w * (cross + cross.T) try: XtX_inv = np.linalg.inv(X.T @ X) except np.linalg.LinAlgError: XtX_inv = np.linalg.pinv(X.T @ X) se = np.sqrt(np.abs(np.diag(XtX_inv @ S @ XtX_inv))) return coefs / np.where(se > 0, se, 1) def newey_west_lags(df: pd.DataFrame) -> None: print("\n--- B. NEWEY-WEST LAG SENSITIVITY ---") rows = [] for target, horizon in [('ret_1d', '1-Day'), ('ret_5d', '5-Day')]: sub = df[FEATURES + [target]].dropna() for f in FEATURES: sub[f] = winsorize(sub[f]) sub[target] = winsorize(sub[target]) X = add_constant(standardize(sub[FEATURES].values)) y = sub[target].values coefs = ols(X, y) for lag in [5, 10, 22]: t = nw_tstats_vectorized(X, y, coefs, lag) for i, name in enumerate(['const'] + FEATURES): rows.append({'target': horizon, 'n_lags': lag, 'variable': name, 'coefficient': coefs[i], 'nw_t_stat': t[i], 'sig_5pct': abs(t[i]) > 1.96}) n_sig = int(np.sum(np.abs(t[1:]) > 1.96)) print(f" {horizon} | NW({lag:2d}): {n_sig}/10 significant at 5%") pd.DataFrame(rows).to_csv( config.RESULTS_DIR / "extended_newey_west_lags.csv", index=False) # -------------------------------------------------------------------------- # C. Daily cross-sectional Spearman information coefficients (stocks) # -------------------------------------------------------------------------- def information_coefficients(stocks: pd.DataFrame) -> None: print("\n--- C. RANK-BASED INFORMATION COEFFICIENTS (stocks) ---") rows = [] for target, horizon in [('ret_1d', '1-Day'), ('ret_5d', '5-Day')]: cols = FEATURES + [target] daily_ics = {f: [] for f in FEATURES} for _, g in stocks[cols + ['trade_date']].groupby('trade_date'): g = g[cols].dropna() if len(g) < 20: continue ranks = g.rank().values ranks = (ranks - ranks.mean(axis=0)) / ranks.std(axis=0) y = ranks[:, -1] for j, f in enumerate(FEATURES): daily_ics[f].append(np.mean(ranks[:, j] * y)) for f in FEATURES: ics = np.array(daily_ics[f]) t_stat = ics.mean() / (ics.std() / np.sqrt(len(ics))) rows.append({ 'target': horizon, 'variable': f, 'mean_ic': ics.mean(), 'std_ic': ics.std(), 't_stat': t_stat, 'pct_positive': (ics > 0).mean(), 'n_days': len(ics), }) print(f" {horizon} | {f:25s}: IC={ics.mean():+.4f}, t={t_stat:7.2f}, " f"%>0={(ics > 0).mean()*100:5.1f}") pd.DataFrame(rows).to_csv( config.RESULTS_DIR / "extended_information_coefficients.csv", index=False) # -------------------------------------------------------------------------- # D. Decile sorts # -------------------------------------------------------------------------- def decile_sorts(stocks: pd.DataFrame) -> None: print("\n--- D. DECILE SORTS (D10 - D1, 5-day returns) ---") rows = [] for sort_var in ['implied_kurtosis_proxy', 'pc_volume_ratio', 'iv_skew_25d']: for n_q, scheme in [(5, 'Quintile (baseline)'), (10, 'Decile')]: res = portfolio_sort(stocks, sort_var, 'ret_5d', n_quantiles=n_q) ls_key = f'LS_{n_q}_1' if res is None or ls_key not in res: continue r = res[ls_key] rows.append({ 'sort_variable': config.SORT_VARIABLES[sort_var], 'scheme': scheme, 'mean_daily_bps': r['mean_daily'] * 10000, 'annualized_return_pct': r['annualized_return'] * 100, 'annualized_vol_pct': r['annualized_vol'] * 100, 'sharpe_ratio': r['sharpe'], 't_statistic': r['t_stat'], 'n_days': r['n_days'], }) print(f" {config.SORT_VARIABLES[sort_var]:25s} | {scheme:20s}: " f"ann.ret={r['annualized_return']*100:7.2f}%, " f"Sharpe={r['sharpe']:6.3f}, t={r['t_stat']:7.2f}") pd.DataFrame(rows).to_csv( config.RESULTS_DIR / "extended_decile_sorts.csv", index=False) # -------------------------------------------------------------------------- # E. Leave-one-year-out # -------------------------------------------------------------------------- def leave_one_year_out(df: pd.DataFrame) -> None: print("\n--- E. LEAVE-ONE-YEAR-OUT PANEL R² ---") years = sorted(df['trade_date'].dt.year.unique()) rows = [] for year in years: sub_df = df[df['trade_date'].dt.year != year] sub = sub_df[FEATURES + ['ret_5d']].dropna() for f in FEATURES: sub[f] = winsorize(sub[f]) sub['ret_5d'] = winsorize(sub['ret_5d']) _, _, r2_ret, n_ret = fit_panel(sub, FEATURES, 'ret_5d') sub_rv = sub_df[RV_FEATURES + ['rv_fwd_1d']].dropna() for f in RV_FEATURES: sub_rv[f] = winsorize(sub_rv[f]) sub_rv['rv_fwd_1d'] = winsorize(sub_rv['rv_fwd_1d']) _, _, r2_rv, _ = fit_panel(sub_rv, RV_FEATURES, 'rv_fwd_1d') rows.append({'excluded_year': year, 'r2_ret_5d': r2_ret, 'r2_rv_har_iv': r2_rv, 'n_obs_ret': n_ret}) print(f" excl. {year}: 5D-ret R²={r2_ret:.4f}, HAR+IV RV R²={r2_rv:.4f}") pd.DataFrame(rows).to_csv( config.RESULTS_DIR / "extended_leave_one_year_out.csv", index=False) # -------------------------------------------------------------------------- # F. Placebo: within-ticker permutation of the feature block # -------------------------------------------------------------------------- def placebo_test(df: pd.DataFrame, n_draws: int = 10) -> None: print("\n--- F. PLACEBO TEST (features permuted within ticker) ---") base = df[FEATURES + ['ret_5d', 'ticker']].dropna().reset_index(drop=True) for f in FEATURES: base[f] = winsorize(base[f]) base['ret_5d'] = winsorize(base['ret_5d']) _, _, r2_actual, n = fit_panel(base, FEATURES, 'ret_5d') print(f" Actual R²: {r2_actual:.4f} (N={n:,})") rng = np.random.default_rng(42) ticker_codes = base['ticker'].astype('category').cat.codes.values feat_matrix = base[FEATURES].values r2_placebos = [] for draw in range(n_draws): perm_matrix = np.empty_like(feat_matrix) for code in np.unique(ticker_codes): idx = np.flatnonzero(ticker_codes == code) perm_matrix[idx] = feat_matrix[idx[rng.permutation(len(idx))]] X = add_constant(standardize(perm_matrix)) y = base['ret_5d'].values coefs = ols(X, y) r2_placebos.append(r_squared(y, X @ coefs)) r2_placebos = np.array(r2_placebos) print(f" Placebo R² over {n_draws} draws: mean={r2_placebos.mean():.6f}, " f"max={r2_placebos.max():.6f}") pd.DataFrame({ 'draw': ['actual'] + [str(i + 1) for i in range(n_draws)], 'r2': [r2_actual] + list(r2_placebos), }).to_csv(config.RESULTS_DIR / "extended_placebo.csv", index=False) def main(): print("=" * 70) print("EXTENDED ROBUSTNESS ANALYSES") print("=" * 70) config.ensure_output_dirs() df = load_merged() stocks = df[~df['ticker'].isin(config.NON_STOCK_TICKERS)].copy() winsorization_sensitivity(df) newey_west_lags(df) information_coefficients(stocks) decile_sorts(stocks) leave_one_year_out(df) placebo_test(df) print("\nEXTENDED ROBUSTNESS COMPLETE.") if __name__ == "__main__": main()