spb/wp7_uqo Public

1#!/usr/bin/env python32# =============================================================================3# Author: Simon-Pierre Boucher4# Contact: contact@spboucher.ai5# =============================================================================6"""Step 07 — Extended descriptive statistics and data-quality analysis.78Panels A–E and G run from the processed panel alone. Panel F (VIX regimes)9and Panel H (raw options data quality) additionally need the raw stores and10are skipped gracefully when those are absent.1112Inputs : data/processed/merged_options_rv.parquet13         (+ index_5min.duckdb and options.duckdb for panels F/H)14Outputs: results/descriptive_*.csv (8 tables)15"""1617import warnings1819import pandas as pd2021import _bootstrap  # noqa: F40122from wp7 import config23from wp7.data_io import RawDataUnavailableError, load_merged, load_vix_daily, open_raw_db2425warnings.filterwarnings('ignore')2627SUMMARY_VARS = [28    'iv_atm_30d', 'iv_atm_90d', 'iv_term_slope', 'iv_skew_25d',29    'implied_skewness', 'implied_kurtosis_proxy',30    'pc_volume_ratio', 'pc_oi_ratio', 'net_gamma_exposure',31    'avg_vega_30d', 'avg_theta_30d', 'total_option_volume', 'total_oi',32    'rv_daily', 'rvol_daily', 'rv_weekly', 'daily_return',33    'realized_skew', 'realized_kurt',34    'ret_1d', 'ret_5d', 'rv_fwd_1d', 'rv_fwd_5d',35]3637CORR_VARS = ['iv_atm_30d', 'iv_term_slope', 'iv_skew_25d', 'implied_skewness',38             'implied_kurtosis_proxy', 'pc_volume_ratio', 'pc_oi_ratio',39             'rv_daily', 'rv_weekly', 'ret_1d', 'ret_5d']4041AUTOCORR_VARS = ['iv_atm_30d', 'iv_skew_25d', 'rv_daily', 'daily_return', 'pc_volume_ratio']424344def group_stats(data: pd.DataFrame, label: str) -> dict:45    """One summary row per asset group."""46    return {47        'Group': label,48        'N_obs': len(data),49        'N_tickers': data['ticker'].nunique(),50        'Date_min': str(data['trade_date'].min().date()),51        'Date_max': str(data['trade_date'].max().date()),52        'Mean_IV_ATM': data['iv_atm_30d'].mean(),53        'Std_IV_ATM': data['iv_atm_30d'].std(),54        'Mean_RV': data['rv_daily'].mean(),55        'Mean_Skew': data['iv_skew_25d'].mean(),56        'Mean_Ret_1d': data['ret_1d'].mean(),57        'Std_Ret_1d': data['ret_1d'].std(),58        'Mean_PC_ratio': data['pc_volume_ratio'].mean(),59    }606162def main():63    print("=" * 70)64    print("EXTENDED DESCRIPTIVE STATISTICS")65    print("=" * 70)66    config.ensure_output_dirs()67    res_dir = config.RESULTS_DIR6869    df = load_merged()70    df['year'] = df['trade_date'].dt.year7172    # ── Panel A: summary statistics ──73    print("\n--- PANEL A: SUMMARY STATISTICS ---")74    summary = df[SUMMARY_VARS].describe(75        percentiles=[0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99]).T76    summary['skewness'] = df[SUMMARY_VARS].skew()77    summary['kurtosis'] = df[SUMMARY_VARS].kurtosis()78    summary['pct_missing'] = df[SUMMARY_VARS].isnull().mean() * 10079    summary.to_csv(res_dir / "descriptive_summary_stats.csv")80    print(summary[['count', 'mean', 'std', '1%', '50%', '99%',81                   'skewness', 'kurtosis', 'pct_missing']].round(4).to_string())8283    # ── Panel B: coverage by year ──84    print("\n--- PANEL B: COVERAGE BY YEAR ---")85    coverage = df.groupby('year').agg(86        n_obs=('ticker', 'count'),87        n_tickers=('ticker', 'nunique'),88        avg_iv_atm=('iv_atm_30d', 'mean'),89        avg_rv=('rv_daily', 'mean'),90        avg_skew=('iv_skew_25d', 'mean'),91        avg_ret=('daily_return', 'mean'),92        std_ret=('daily_return', 'std'),93    ).round(6)94    coverage.to_csv(res_dir / "descriptive_coverage_by_year.csv")95    print(coverage.to_string())9697    # ── Panel C: coverage by asset group ──98    print("\n--- PANEL C: BY ASSET GROUP ---")99    stocks_list = [t for t in df['ticker'].unique() if t not in config.NON_STOCK_TICKERS]100    groups = pd.DataFrame([101        group_stats(df[df['ticker'].isin(stocks_list)], 'Stocks'),102        group_stats(df[df['ticker'].isin(config.ETF_TICKERS)], 'ETFs'),103        group_stats(df[df['ticker'].isin(config.INDEX_OPTION_TICKERS)], 'Indices'),104        group_stats(df, 'All'),105    ])106    groups.to_csv(res_dir / "descriptive_by_group.csv", index=False)107    print(groups.to_string(index=False))108109    # ── Panel D: correlation matrix ──110    print("\n--- PANEL D: CORRELATION MATRIX ---")111    corr_matrix = df[CORR_VARS].corr().round(3)112    corr_matrix.to_csv(res_dir / "descriptive_correlation_matrix.csv")113    print(corr_matrix.to_string())114115    # ── Panel E: autocorrelation structure ──116    print("\n--- PANEL E: AUTOCORRELATION STRUCTURE ---")117    autocorr_results = []118    for var in AUTOCORR_VARS:119        for lag in [1, 5, 10, 22]:120            ac = df.groupby('ticker')[var].apply(lambda x: x.autocorr(lag=lag)).mean()121            autocorr_results.append({'variable': var, 'lag': lag, 'avg_autocorr': ac})122    autocorr_df = pd.DataFrame(autocorr_results)123    autocorr_df.to_csv(res_dir / "descriptive_autocorrelations.csv", index=False)124    print(autocorr_df.pivot(index='variable', columns='lag', values='avg_autocorr')125          .round(4).to_string())126127    # ── Panel G: cross-sectional dispersion by year ──128    print("\n--- PANEL G: CROSS-SECTIONAL DISPERSION ---")129    cs_disp = df.groupby('year').agg(130        iv_atm_cs_std=('iv_atm_30d', 'std'),131        skew_cs_std=('iv_skew_25d', 'std'),132        rv_cs_std=('rv_daily', 'std'),133        ret_cs_std=('daily_return', 'std'),134        n_tickers=('ticker', 'nunique'),135    ).round(6)136    cs_disp.to_csv(res_dir / "descriptive_cross_sectional_dispersion.csv")137    print(cs_disp.to_string())138139    # ── Panel F: statistics by VIX regime (needs raw VIX series) ──140    print("\n--- PANEL F: STATISTICS BY VIX REGIME ---")141    try:142        vix = load_vix_daily()143        dfv = df.merge(vix, on='trade_date', how='left')144        dfv['vix_regime'] = pd.cut(dfv['vix_close'], bins=config.VIX_REGIME_BINS,145                                   labels=config.VIX_REGIME_LABELS_VERBOSE)146        regime_stats = dfv.groupby('vix_regime', observed=True).agg(147            n_obs=('ticker', 'count'),148            pct_obs=('ticker', lambda x: len(x) / len(dfv) * 100),149            mean_iv_atm=('iv_atm_30d', 'mean'),150            mean_rv=('rv_daily', 'mean'),151            mean_skew=('iv_skew_25d', 'mean'),152            mean_ret_1d=('ret_1d', 'mean'),153            std_ret_1d=('ret_1d', 'std'),154            mean_pc_ratio=('pc_volume_ratio', 'mean'),155            mean_impl_skew=('implied_skewness', 'mean'),156        ).round(6)157        regime_stats.to_csv(res_dir / "descriptive_vix_regimes.csv")158        print(regime_stats.to_string())159    except RawDataUnavailableError as exc:160        print(f"  [Panel F skipped — raw stores unavailable]\n  {exc}")161162    # ── Panel H: raw options data quality (needs options.duckdb) ──163    print("\n--- PANEL H: OPTIONS DATA QUALITY ---")164    try:165        opt_con = open_raw_db("options")166        quality = opt_con.execute("""167            SELECT168                EXTRACT(YEAR FROM trade_date) AS year,169                COUNT(*) AS n_records,170                COUNT(DISTINCT ticker) AS n_tickers,171                AVG(CASE WHEN bid_iv > 0 AND ask_iv > 0 THEN 1.0 ELSE 0.0 END) AS pct_valid_iv,172                AVG(CASE WHEN volume > 0 THEN 1.0 ELSE 0.0 END) AS pct_with_volume,173                AVG(CASE WHEN open_interest > 0 THEN 1.0 ELSE 0.0 END) AS pct_with_oi,174                AVG(CASE WHEN delta IS NOT NULL AND delta != 0 THEN 1.0 ELSE 0.0 END) AS pct_valid_greeks,175                AVG(ask_price - bid_price) AS avg_spread176            FROM option_chain177            GROUP BY year178            ORDER BY year179        """).fetchdf()180        opt_con.close()181        quality.to_csv(res_dir / "descriptive_options_quality.csv", index=False)182        print(quality.round(4).to_string(index=False))183    except RawDataUnavailableError as exc:184        print(f"  [Panel H skipped — raw stores unavailable]\n  {exc}")185186    print("\nDESCRIPTIVE STATISTICS COMPLETE.")187188189if __name__ == "__main__":190    main()191