spb/wp7_uqo Public

1#!/usr/bin/env python32# =============================================================================3# Author: Simon-Pierre Boucher4# Contact: contact@spboucher.ai5# =============================================================================6"""Step 02 — RQ1: Do option-implied moments predict next-day/next-week returns?78Pooled OLS (HC1) and Fama-MacBeth panel regressions across asset groups9(stocks, ETFs, indices, all).1011Inputs : data/processed/merged_options_rv.parquet12Outputs: results/rq1_regression_results.csv, results/rq1_meta.csv13"""1415import warnings1617import numpy as np18import pandas as pd1920import _bootstrap  # noqa: F40121from wp7 import config22from wp7.data_io import load_merged23from wp7.econometrics import fama_macbeth, winsorize2425warnings.filterwarnings('ignore')2627TARGETS = ['ret_1d', 'ret_5d']282930def panel_ols(data, features, target, label=""):31    """Pooled OLS with HC1 t-statistics, as specified in the original study.3233    Historical quirks preserved on purpose (they define the published34    numbers): features are standardized with the *pandas* sample std35    (ddof=1) and missing standardized cells are zero-filled; the reported36    ``p_value`` column is an ad-hoc normal-tail approximation, not an exact37    two-sided p-value (significance flags use the usual 1.96/2.576 cutoffs).38    """39    sub = data[features + [target, 'ticker', 'trade_date']].dropna()40    if len(sub) < 100:41        return None4243    for f in features:44        sub[f] = winsorize(sub[f])45    sub[target] = winsorize(sub[target])4647    means = sub[features].mean()48    stds = sub[features].std()49    X = (sub[features] - means) / stds50    X = X.fillna(0)51    X.insert(0, 'const', 1.0)52    y = sub[target].values5354    coefs, _, _, _ = np.linalg.lstsq(X.values, y, rcond=None)5556    y_pred = X.values @ coefs57    ss_res = np.sum((y - y_pred) ** 2)58    ss_tot = np.sum((y - y.mean()) ** 2)59    r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 06061    n = len(y)62    k = len(features)63    adj_r2 = 1 - (1 - r2) * (n - 1) / (n - k - 1)6465    # HC1 sandwich (X' diag(e²) X via broadcasting)66    e = y - y_pred67    XtX_inv = np.linalg.inv(X.values.T @ X.values)68    S = (X.values * (e ** 2)[:, None]).T @ X.values * n / (n - k - 1)69    se = np.sqrt(np.diag(XtX_inv @ S @ XtX_inv))70    t_stats = coefs / se7172    results = pd.DataFrame({73        'variable': ['const'] + features,74        'coefficient': coefs,75        'std_error': se,76        't_stat': t_stats,77        'p_value': 2 * (1 - pd.Series(np.abs(t_stats)).apply(78            lambda x: min(1.0, 0.5 * np.exp(-0.5 * x**2) * np.sqrt(2 / np.pi) if x < 30 else 0)79        )).values,80        'significant_5pct': np.abs(t_stats) > 1.96,81        'significant_1pct': np.abs(t_stats) > 2.576,82    })8384    meta = {85        'label': label,86        'target': target,87        'n_obs': n,88        'n_tickers': sub['ticker'].nunique(),89        'r_squared': r2,90        'adj_r_squared': adj_r2,91    }92    return results, meta939495def main():96    print("=" * 70)97    print("RQ1: OPTION-IMPLIED MOMENTS AND RETURN PREDICTABILITY")98    print("=" * 70)99    config.ensure_output_dirs()100101    df = load_merged()102    print(f"Loaded {len(df):,} rows, {df['ticker'].nunique()} tickers")103104    stocks = df[~df['ticker'].isin(config.NON_STOCK_TICKERS)].copy()105    etfs = df[df['ticker'].isin(config.ETF_TICKERS)].copy()106    indices = df[df['ticker'].isin(config.INDEX_OPTION_TICKERS)].copy()107108    features = config.RQ1_FEATURES109    all_results, all_meta = [], []110111    for group_name, group_data in [("Stocks", stocks), ("ETFs", etfs),112                                   ("Indices", indices), ("All", df)]:113        for target in TARGETS:114            horizon = "1-Day" if target == 'ret_1d' else "5-Day"115            label = f"{group_name} | {horizon}"116117            # Pooled OLS118            res = panel_ols(group_data, features, target, label)119            if res:120                r, m = res121                r['group'] = group_name122                r['horizon'] = horizon123                r['method'] = 'Pooled OLS'124                all_results.append(r)125                all_meta.append(m)126                print(f"\n{label} (Pooled OLS): R²={m['r_squared']:.6f}, "127                      f"Adj-R²={m['adj_r_squared']:.6f}, N={m['n_obs']:,}")128                sig = r[r['significant_5pct'] & (r['variable'] != 'const')]129                if len(sig) > 0:130                    print(f"  Significant predictors: {', '.join(sig['variable'].tolist())}")131132            # Fama-MacBeth (stocks and full panel only)133            if group_name in ['Stocks', 'All']:134                res_fm = fama_macbeth(group_data, features, target)135                if res_fm:136                    r_fm, n_periods, n_tickers = res_fm137                    r_fm['group'] = group_name138                    r_fm['horizon'] = horizon139                    r_fm['method'] = 'Fama-MacBeth'140                    all_results.append(r_fm)141                    all_meta.append({'label': label, 'target': target,142                                     'n_periods': n_periods, 'n_tickers': n_tickers})143                    print(f"  Fama-MacBeth: N_periods={n_periods}")144                    sig_fm = r_fm[r_fm['fm_significant_5pct'] & (r_fm['variable'] != 'const')]145                    if len(sig_fm) > 0:146                        print(f"  FM Significant: {', '.join(sig_fm['variable'].tolist())}")147148    results_df = pd.concat(all_results, ignore_index=True)149    results_df.to_csv(config.RESULTS_DIR / "rq1_regression_results.csv", index=False)150151    meta_df = pd.DataFrame(all_meta)152    meta_df.to_csv(config.RESULTS_DIR / "rq1_meta.csv", index=False)153154    print("\n" + "=" * 70)155    print("RQ1 SUMMARY TABLE")156    print("=" * 70)157    summary_rows = [158        {'Group': m.get('label', ''),159         'N': m.get('n_obs', m.get('n_periods', '')),160         'R²': f"{m.get('r_squared', 0):.6f}",161         'Adj-R²': f"{m.get('adj_r_squared', 0):.6f}"}162        for _, m in meta_df.iterrows() if 'r_squared' in m163    ]164    if summary_rows:165        print(pd.DataFrame(summary_rows).to_string(index=False))166167    print("\nRQ1 COMPLETE.")168169170if __name__ == "__main__":171    main()172