Repository Public Research preview
1#!/usr/bin/env python32# =============================================================================3# Author: Simon-Pierre Boucher4# Contact: contact@spboucher.ai5# =============================================================================6"""Step 14 — Generate every LaTeX table fragment of the paper from results/.78Added in the 2026-10-17 revision so that no number in the manuscript is typed9by hand: each fragment in paper/tables/ is \\input by the sections. Tables10that depend on raw-only outputs (RQ4 Greeks decay, RQ5, feature importance)11are built from the frozen CSVs shipped in results/ and are labelled as such12in the paper.1314Inputs : results/*.csv15Outputs: paper/tables/*.tex16"""1718import numpy as np19import pandas as pd2021import _bootstrap # noqa: F40122from wp7 import config2324OUT = config.PROJECT_ROOT / "paper" / "tables"25R = config.RESULTS_DIR2627VAR_LABEL = {28 'iv_atm_30d': r'$IV_{ATM,30d}$', 'iv_term_slope': 'IV term slope',29 'iv_skew_25d': r'Skew$_{25\delta}$', 'implied_skewness': 'Implied skewness',30 'implied_kurtosis_proxy': 'Implied kurtosis', 'pc_volume_ratio': 'PC volume ratio',31 'pc_oi_ratio': 'PC OI ratio', 'net_gamma_exposure': 'Net gamma exposure',32 'rv_daily': r'$RV_t$', 'rv_w': r'$RV_t^{(w)}$', 'rv_lag1': r'$RV_{t-1}$',33 'log_volume': 'log option volume', 'log_oi': 'log open interest',34 'abs_return': r'$|r_t|$', 'ret_lag1': r'$r_{t}$ (lagged return)', 'ret_lag5': r'$\sum_{j=1}^{5} r_{t-j}$',35 'const': 'Constant',36}37SORT_LABEL = {38 'ATM IV (30d)': r'$IV_{ATM,30d}$', 'Volatility Skew (25d)': r'Skew$_{25\delta}$',39 'Implied Skewness': 'Implied skewness', 'Implied Kurtosis': 'Implied kurtosis',40 'Put-Call Volume Ratio': 'PC volume ratio', 'Put-Call OI Ratio': 'PC OI ratio',41 'IV Term Structure Slope': 'IV term slope',42}43MODEL_LABEL = {44 'RV + r2 (ARCH proxy)': r'$RV_t + r_t^2$ (ARCH-type)', 'HAR-RV': 'HAR-RV',45 'IV only': 'IV only (30d, 90d)', 'IV surface': 'IV surface',46 'HAR-RV + IV surface': 'HAR-RV $+$ IV surface',47}48SUBPERIOD_LABEL = {49 'Pre-GFC Recovery (2010-2012)': 'Post-GFC (2010--12)', 'Bull Market (2013-2016)': 'Bull (2013--16)',50 'Low Vol Era (2017-2018)': 'Low vol (2017--18)', 'Pre-COVID (2019)': 'Pre-COVID (2019)',51 'COVID Period (2020)': 'COVID (2020)', 'Post-COVID Bull (2021)': 'Post-COVID (2021)',52 'Rate Hiking (2022)': 'Rate hikes (2022)', 'Recovery (2023-2024)': 'Recovery (2023--24)',53 'Recent (2025)': 'Recent (2025)',54}55REGIME_LABEL = {'VeryLow': 'Very low ($<$15)', 'Low': 'Low (15--20)', 'Medium': 'Medium (20--25)',56 'High': 'High (25--35)', 'Crisis': 'Crisis ($>$35)'}575859def f(x, d=2, pct=False, plus=False):60 if x is None or (isinstance(x, float) and np.isnan(x)):61 return '--'62 s = f"{x * (100 if pct else 1):{'+' if plus else ''}.{d}f}"63 return s.replace('-', r'$-$')646566def stars(t):67 a = abs(t)68 return r'\sym{***}' if a > 2.576 else (r'\sym{**}' if a > 1.96 else (r'\sym{*}' if a > 1.645 else ''))697071def num(n):72 return f"{int(n):,}".replace(',', '{,}')737475def write(name, body):76 (OUT / f"{name}.tex").write_text(body)77 print(f" wrote paper/tables/{name}.tex")787980def table(cols, header_rows, body_rows, notes, caption, label, size=r'\small', colsep=None):81 lines = [r'\begin{table}[H]', r'\centering', rf'\caption{{{caption}}}', rf'\label{{{label}}}',82 r'\begin{threeparttable}', size]83 if colsep:84 lines.append(rf'\setlength{{\tabcolsep}}{{{colsep}}}')85 lines += [rf'\begin{{tabular}}{{{cols}}}', r'\toprule']86 lines += header_rows + [r'\midrule'] + body_rows + [r'\bottomrule', r'\end{tabular}']87 if notes:88 lines += [r'\begin{tablenotes}[flushleft]\footnotesize', rf'\item \textit{{Notes.}} {notes}', r'\end{tablenotes}']89 lines += [r'\end{threeparttable}', r'\end{table}']90 return '\n'.join(lines) + '\n'919293# --------------------------------------------------------------------------94def tab_rq1_panel():95 m = pd.read_csv(R / "rq1_meta.csv")96 m = m[m['r_squared'].notna()]97 rows = []98 for grp in ['Stocks', 'ETFs', 'Indices', 'All']:99 r1 = m[m['label'] == f'{grp} | 1-Day'].iloc[0]100 r5 = m[m['label'] == f'{grp} | 5-Day'].iloc[0]101 rows.append(f"{grp} & {num(r1['n_obs'])} & {int(r1['n_tickers'])} & "102 f"{f(r1['r_squared'],2,True)} & {f(r1['adj_r_squared'],2,True)} & {int(r1['n_sig_5pct'])} & {int(r1['n_sig_bonf'])} & "103 f"{f(r5['r_squared'],2,True)} & {f(r5['adj_r_squared'],2,True)} & {int(r5['n_sig_5pct'])} & {int(r5['n_sig_bonf'])} \\\\")104 hdr = [r'& & & \multicolumn{4}{c}{1-day return} & \multicolumn{4}{c}{5-day return} \\',105 r'\cmidrule(lr){4-7}\cmidrule(lr){8-11}',106 r'Group & $N$ & Tickers & $R^2$ (\%) & Adj.\ $R^2$ (\%) & \#sig$_{5\%}$ & \#sig$_{B}$ & $R^2$ (\%) & Adj.\ $R^2$ (\%) & \#sig$_{5\%}$ & \#sig$_{B}$ \\']107 notes = (r"Pooled OLS of the $h$-day forward log return on the ten standardized predictors of Table~\ref{tab:variables} "108 r"(plus $RV_t$ and $RV_t^{(w)}$), HC1 standard errors; all variables winsorized at the 1st and 99th percentiles. "109 r"The 5-day return is $\sum_{j=1}^{5} r_{t+j}$. \#sig$_{5\%}$ counts predictors with $|t|>1.96$; \#sig$_{B}$ counts those "110 r"clearing the Bonferroni threshold for the 20 tests of each horizon pair ($|t|>3.02$). $N$ is the number of ticker-days with all ten predictors observed.")111 write("tab_rq1_panel", table('@{}lrrrrrrrrrr@{}', hdr, rows, notes,112 'Return predictability: pooled panel regressions by asset group.', 'tab:rq1',113 size=r'\footnotesize', colsep='4pt'))114115116def tab_rq1_coef():117 r = pd.read_csv(R / "rq1_regression_results.csv")118 r = r[r['group'] == 'Stocks']119 rows = []120 for v in config.RQ1_FEATURES:121 cells = [VAR_LABEL[v]]122 for hz in ['1-Day', '5-Day']:123 po = r[(r['method'] == 'Pooled OLS') & (r['horizon'] == hz) & (r['variable'] == v)].iloc[0]124 fm = r[(r['method'] == 'Fama-MacBeth') & (r['horizon'] == hz) & (r['variable'] == v)].iloc[0]125 cells += [f"{f(po['coefficient_bps'],1)}{stars(po['t_stat'])}", f"({f(po['t_stat'],2)})",126 f"{f(fm['fm_coefficient']*1e4,1)}{stars(fm['fm_t_stat'])}", f"({f(fm['fm_t_stat'],2)})"]127 rows.append(' & '.join(cells) + r' \\')128 meta = pd.read_csv(R / "rq1_meta.csv")129 n1 = meta[(meta['label'] == 'Stocks | 1-Day') & meta['n_periods'].notna()]['n_periods'].iloc[0]130 n5 = meta[(meta['label'] == 'Stocks | 5-Day') & meta['n_periods'].notna()]['n_periods'].iloc[0]131 hdr = [r'& \multicolumn{4}{c}{1-day return} & \multicolumn{4}{c}{5-day return} \\', r'\cmidrule(lr){2-5}\cmidrule(lr){6-9}',132 r'& \multicolumn{2}{c}{Pooled OLS} & \multicolumn{2}{c}{Fama--MacBeth} & \multicolumn{2}{c}{Pooled OLS} & \multicolumn{2}{c}{Fama--MacBeth} \\',133 r'\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}\cmidrule(lr){8-9}',134 r'Predictor & $b$ & $(t)$ & $b$ & $(t)$ & $b$ & $(t)$ & $b$ & $(t)$ \\']135 notes = (rf"Individual stocks (49 tickers, 79{{,}}943 ticker-days). Coefficients are in basis points of $h$-day return per one-standard-deviation "136 rf"change in the predictor. Pooled OLS: HC1 $t$-statistics. Fama--MacBeth: daily cross-sectional regressions "137 rf"({int(n1):,} days at 1 day, {int(n5):,} at 5 days; at least 10 stocks per day); $t$-statistics from the time series of slopes with "138 rf"Newey--West standard errors (0 lags at 1 day, 4 lags at 5 days). \sym{{***}}, \sym{{**}}, \sym{{*}}: $|t|>2.576$, $1.96$, $1.645$.").replace(',', '{,}').replace('{,}}', ',}')139 notes = notes.replace('{,} ', ', ').replace('{,}$', ',$')140 write("tab_rq1_coef", table('@{}lrrrrrrrr@{}', hdr, rows, notes,141 'Return predictability for individual stocks: coefficients and $t$-statistics.', 'tab:rq1_coef', size=r'\footnotesize'))142143144def tab_sorts():145 d = pd.read_csv(R / "portfolio_sort_results.csv")146 ls = d[d['quintile'] == 'L/S(5-1)']147 rows = []148 for hz in ['5-Day', '1-Day']:149 rows.append(rf'\multicolumn{{7}}{{@{{}}l}}{{\textit{{Panel {"A" if hz=="5-Day" else "B"}: {hz.lower()} forward returns}}}} \\')150 sub = ls[ls['return_horizon'] == hz].sort_values('sharpe_ratio', ascending=False)151 for _, r in sub.iterrows():152 rows.append(f"{SORT_LABEL[r['sort_variable']]} & {f(r['annualized_return_pct'],2)} & {f(r['annualized_vol_pct'],2)} & "153 f"{f(r['sharpe_ratio'],3)} & {f(r['t_statistic'],2)}{stars(r['t_statistic'])} & {f(r['t_statistic_iid'],2)} & {num(r['n_days'])} \\\\")154 if hz == '5-Day':155 rows.append(r'\addlinespace')156 hdr = [r'Sorting variable & Ann.\ ret.\ (\%) & Ann.\ vol.\ (\%) & Sharpe & $t_{NW}$ & $t_{iid}$ & $N_{days}$ \\']157 notes = (r"Equal-weighted quintile portfolios of individual stocks formed on each trading day; the long-short portfolio is Q5$-$Q1. "158 r"Panel A: 5-day forward returns $\sum_{j=1}^{5} r_{t+j}$, annualized with 52 periods; the daily series of 5-day returns overlaps, so "159 r"$t_{NW}$ uses Newey--West standard errors with 4 lags while $t_{iid}$ ignores the overlap (the statistic reported in version~1.1). "160 r"Panel B: next-day returns, annualized with 252 periods. \sym{***}, \sym{**}, \sym{*}: $|t_{NW}|>2.576$, $1.96$, $1.645$.")161 write("tab_sorts", table('@{}lrrrrrr@{}', hdr, rows, notes,162 'Long-short quintile portfolios sorted on option-implied variables (individual stocks).', 'tab:sorts'))163164165def tab_dsort():166 d = pd.read_csv(R / "double_sort_iv_skew.csv")167 def cell(a, b):168 sel = d[(d['q1'].astype(str) == str(a)) & (d['q2'].astype(str) == str(b))]169 if sel.empty:170 return '', ''171 r = sel.iloc[0]172 return f"{f(r['mean_bps'],1)}{stars(r['t_stat'])}", f"({f(r['t_stat'],2)})"173 rows = []174 for a, lab in [(1, 'Low IV'), (2, 'Medium IV'), (3, 'High IV'), ('H-L', r'\textit{High$-$Low IV}')]:175 m, t = zip(*[cell(a, b) for b in (1, 2, 3, 'H-L')])176 rows.append(f"{lab} & " + ' & '.join(m) + r' \\')177 rows.append(' & ' + ' & '.join(t) + r' \\')178 if a == 3:179 rows.append(r'\midrule')180 hdr = [r'& Low skew & Medium & High skew & \textit{High$-$Low skew} \\']181 notes = (r"Independent daily tercile sorts of individual stocks on $IV_{ATM,30d}$ and implied skewness; cells report the time-series mean of the daily "182 r"equal-weighted cell portfolio's 5-day forward return in basis points per week, with Newey--West (4 lags) $t$-statistics in parentheses. "183 r"The last row and column are differences between the daily High and Low portfolio returns. "184 r"\sym{***}, \sym{**}, \sym{*}: $|t|>2.576$, $1.96$, $1.645$.")185 write("tab_dsort", table('@{}lrrrr@{}', hdr, rows, notes,186 r'Double sort: ATM\,IV $\times$ implied skewness $\to$ 5-day returns (bps/week).', 'tab:dsort'))187188189def tab_calendar():190 ct = pd.read_csv(R / "calendar_time_strategies.csv")191 tc = pd.read_csv(R / "transaction_cost_analysis.csv")192 rows = []193 for h, pan in [(5, 'A'), (1, 'B')]:194 rows.append(rf'\multicolumn{{9}}{{@{{}}l}}{{\textit{{Panel {pan}: {h}-day holding period}}}} \\')195 sub = ct[ct['holding_days'] == h].sort_values('sharpe', ascending=False)196 for _, r in sub.iterrows():197 net = tc[(tc['holding_days'] == h) & (tc['sort_variable'] == r['sort_variable'])].set_index('tc_bps')198 rows.append(f"{SORT_LABEL[r['sort_variable']]} & {f(r['ann_return_pct'],2)} & {f(r['ann_vol_pct'],2)} & {f(r['sharpe'],2)} & "199 f"{f(r['t_nw'],2)}{stars(r['t_nw'])} & {f(r['avg_daily_turnover'],2)} & {f(r['breakeven_cost_bps'],1)} & "200 f"{f(net.loc[10,'sharpe'],2)} & {f(net.loc[25,'sharpe'],2)} \\\\")201 if h == 5:202 rows.append(r'\addlinespace')203 hdr = [r'& \multicolumn{4}{c}{Gross} & & & \multicolumn{2}{c}{Net Sharpe} \\', r'\cmidrule(lr){2-5}\cmidrule(lr){8-9}',204 r'Sorting variable & Ret.\ (\%) & Vol.\ (\%) & Sharpe & $t_{NW}$ & Turnover & Break-even (bps) & 10 bps & 25 bps \\']205 notes = (r"Calendar-time implementation of the long-short quintile strategies (49 stocks): each day a new Q5$-$Q1 cohort is formed and held for $h$ "206 r"trading days, so the strategy holds $h$ overlapping cohorts with weight $1/h$ each and its daily return series does not overlap "207 r"\citep{jegadeesh1993returns}. Returns and volatilities are annualized with 252 days; $t_{NW}$ uses 5 Newey--West lags. "208 r"Turnover is the average daily two-sided traded value $\sum_i |W_{t,i}-W_{t-1,i}|$ per unit of long leg; the break-even cost is the one-way "209 r"proportional cost (basis points of traded value) at which the mean net return is zero; net Sharpe ratios deduct one-way costs of 10 and 25 bps "210 r"\citep{novymarx2016taxonomy}. \sym{***}, \sym{**}, \sym{*}: $|t_{NW}|>2.576$, $1.96$, $1.645$.")211 write("tab_calendar", table('@{}lrrrrrrrr@{}', hdr, rows, notes,212 'Calendar-time long-short strategies: gross performance, turnover and net-of-cost Sharpe ratios.', 'tab:calendar', size=r'\footnotesize'))213214215def tab_rq2_insample():216 d = pd.read_csv(R / "rq2_insample.csv")217 rows = []218 for m in config.RQ2_MODELS:219 r1 = d[(d['target'] == '1-Day RV') & (d['model'] == m)].iloc[0]220 r5 = d[(d['target'] == '5-Day RV') & (d['model'] == m)].iloc[0]221 rows.append(f"{MODEL_LABEL[m]} & {int(r1['n_features'])} & {f(r1['r2'],3)} & {f(r1['delta_r2_vs_har_pct'],1,plus=True)} & "222 f"{f(r5['r2'],3)} & {f(r5['delta_r2_vs_har_pct'],1,plus=True)} \\\\")223 n = int(d['n_obs'].iloc[0])224 hdr = [r'& & \multicolumn{2}{c}{$RV_{t+1}$} & \multicolumn{2}{c}{$\sum_{j=1}^{5} RV_{t+j}$} \\', r'\cmidrule(lr){3-4}\cmidrule(lr){5-6}',225 r'Model & $k$ & $R^2$ & $\Delta R^2$ (\%) & $R^2$ & $\Delta R^2$ (\%) \\']226 notes = (rf"Pooled OLS on the common sample of {num(n)} ticker-days (69 tickers) for which every regressor of the largest model is observed; "227 r"regressors and targets winsorized at the 1st and 99th percentiles within ticker and standardized. $k$ is the number of regressors. "228 r"$\Delta R^2$ is the change relative to HAR-RV in percent of the HAR-RV $R^2$. Model definitions in Section~\ref{sec:method}.")229 write("tab_rq2_insample", table('@{}lrrrrr@{}', hdr, rows, notes,230 'In-sample fit of realized-variance forecasting models (common sample).', 'tab:rq2'))231232233def tab_rq2_oos():234 oos = pd.read_csv(R / "rq2_oos_results.csv")235 dm = pd.read_csv(R / "rq2_diebold_mariano.csv")236 rows = []237 for tgt, pan, lab in [('1-Day RV', 'A', r'$RV_{t+1}$'), ('5-Day RV', 'B', r'$\sum_{j=1}^{5} RV_{t+j}$')]:238 rows.append(rf'\multicolumn{{6}}{{@{{}}l}}{{\textit{{Panel {pan}: target {lab}}}}} \\')239 sub = oos[oos['target'] == tgt]240 n_t = sub['ticker'].nunique()241 for m in config.RQ2_MODELS:242 s = sub[sub['model'] == m]243 rows.append(f"{MODEL_LABEL[m]} & {f(s['r2_oos'].median(),3)} & {f(s['r2_oos'].mean(),3)} & "244 f"{f(s['qlike'].median(),3)} & {f(s['qlike'].mean(),3)} & {int((s['r2_oos'] > 0).sum())}/{n_t} \\\\")245 d = dm[dm['target'] == tgt]246 rows.append(rf"\multicolumn{{6}}{{@{{}}l}}{{\quad HAR-RV $+$ IV surface vs.\ HAR-RV: median MSE ratio {f(d['mse_ratio'].median(),3)}, "247 rf"DM$_{{MSE}}>1.96$ for {int((d['dm_mse']>1.96).sum())}/{len(d)} tickers, $<-1.96$ for {int((d['dm_mse']<-1.96).sum())}; }} \\")248 rows.append(rf"\multicolumn{{6}}{{@{{}}l}}{{\quad median QLIKE ratio {f(d['qlike_ratio'].median(),3)}, "249 rf"DM$_{{QLIKE}}>1.96$ for {int((d['dm_qlike']>1.96).sum())}/{len(d)}, $<-1.96$ for {int((d['dm_qlike']<-1.96).sum())}.}} \\")250 if tgt == '1-Day RV':251 rows.append(r'\addlinespace')252 hdr = [r'& \multicolumn{2}{c}{$R^2_{OOS}$} & \multicolumn{2}{c}{QLIKE} & \\', r'\cmidrule(lr){2-3}\cmidrule(lr){4-5}',253 r'Model & Median & Mean & Median & Mean & $R^2_{OOS}>0$ \\']254 notes = (r"Rolling-window forecasts for the 30 tickers with the longest histories: each model is re-estimated by OLS on 500 trading days and used for the "255 r"next 250 days, the window then advancing by 250 days; the stacked forecasts yield per-ticker $R^2_{OOS}=1-\sum(RV-\widehat{RV})^2/\sum(RV-\overline{RV})^2$ "256 r"and QLIKE $=\frac{1}{T}\sum[RV/\widehat{RV}-\ln(RV/\widehat{RV})-1]$ \citep{patton2011volatility}. Columns report the median and mean across tickers "257 r"and the number of tickers with positive $R^2_{OOS}$. Diebold--Mariano statistics \citep{diebold1995comparing} are computed per ticker on the loss "258 r"differential with Newey--West variance ($h$ lags); positive values favour the IV-augmented model. Per-ticker statistics are in Table~\ref{tab:app_dm}.")259 write("tab_rq2_oos", table('@{}lrrrrr@{}', hdr, rows, notes,260 'Out-of-sample realized-variance forecasting: HAR-RV against IV-surface models.', 'tab:rq2_oos', size=r'\footnotesize'))261262263def tab_app_dm():264 dm = pd.read_csv(R / "rq2_diebold_mariano.csv")265 t1 = dm[dm['target'] == '1-Day RV'].set_index('ticker')266 t5 = dm[dm['target'] == '5-Day RV'].set_index('ticker')267 rows = []268 for tk in sorted(t1.index):269 a, b = t1.loc[tk], t5.loc[tk]270 rows.append(f"{tk} & {num(a['n_forecasts'])} & {f(a['mse_ratio'],3)} & {f(a['dm_mse'],2)} & {f(a['qlike_ratio'],3)} & {f(a['dm_qlike'],2)} & "271 f"{f(b['mse_ratio'],3)} & {f(b['dm_mse'],2)} & {f(b['qlike_ratio'],3)} & {f(b['dm_qlike'],2)} \\\\")272 hdr = [r'& & \multicolumn{4}{c}{$RV_{t+1}$} & \multicolumn{4}{c}{$\sum_{j=1}^{5} RV_{t+j}$} \\', r'\cmidrule(lr){3-6}\cmidrule(lr){7-10}',273 r'Ticker & $T$ & MSE ratio & DM & QLIKE ratio & DM & MSE ratio & DM & QLIKE ratio & DM \\']274 notes = (r"Ratios are losses of HAR-RV $+$ IV surface over HAR-RV (values below one favour the IV-augmented model); DM statistics are positive when the "275 r"IV-augmented model has the lower loss. $T$ is the number of out-of-sample forecasts. Rolling scheme as in Table~\ref{tab:rq2_oos}.")276 write("tab_app_dm", table('@{}lrrrrrrrrr@{}', hdr, rows, notes,277 'Per-ticker Diebold--Mariano tests: HAR-RV $+$ IV surface versus HAR-RV.', 'tab:app_dm', size=r'\scriptsize'))278279280def tab_rq3():281 d = pd.read_csv(R / "rq3_stress_prediction.csv")282 specs = [('A: ratio', 'A', r'$\hat\rho_{impl}/\hat\rho_{real}$'), ('B: divergence', 'B', r'$\hat\rho_{impl}-\hat\rho_{real}$'),283 ('C: components', 'C', 'implied and realized correlation entered separately')]284 lab = {'corr_ratio': r'Corr.\ ratio', 'corr_divergence': r'Corr.\ divergence', 'implied_corr': r'$\hat\rho_{impl}$',285 'realized_corr': r'$\hat\rho_{real}$', 'vix_close': 'VIX level', 'spx_iv_atm': 'SPX ATM IV'}286 rows = []287 for spec, pan, desc in specs:288 rows.append(rf'\multicolumn{{7}}{{@{{}}l}}{{\textit{{Panel {pan}: {desc}}}}} \\')289 s = d[d['spec'] == spec]290 for v in [x for x in s['variable'].unique() if x != 'const']:291 cells = [lab[v]]292 for hz in ['5d', '10d', '20d']:293 r = s[(s['horizon'] == hz) & (s['variable'] == v)].iloc[0]294 cells += [f"{f(r['coefficient']*100,2)}{stars(r['t_nw'])}", f"({f(r['t_nw'],2)})"]295 rows.append(' & '.join(cells) + r' \\')296 r2 = [s[(s['horizon'] == hz)]['r2'].iloc[0] for hz in ['5d', '10d', '20d']]297 rows.append(r'$R^2$ & ' + ' & '.join([rf'\multicolumn{{2}}{{c}}{{{f(x,3)}}}' for x in r2]) + r' \\')298 rows.append(r'\addlinespace')299 ns = [d[d['horizon'] == hz]['n_obs'].iloc[0] for hz in ['5d', '10d', '20d']]300 br = [d[d['horizon'] == hz]['base_rate'].iloc[0] for hz in ['5d', '10d', '20d']]301 rows.append(r'$N$ (days) & ' + ' & '.join([rf'\multicolumn{{2}}{{c}}{{{num(x)}}}' for x in ns]) + r' \\')302 rows.append(r'Base rate & ' + ' & '.join([rf'\multicolumn{{2}}{{c}}{{{f(x,3)}}}' for x in br]) + r' \\')303 hdr = [r'& \multicolumn{2}{c}{5-day} & \multicolumn{2}{c}{10-day} & \multicolumn{2}{c}{20-day} \\', r'\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}',304 r'Predictor & $b$ & $(t)$ & $b$ & $(t)$ & $b$ & $(t)$ \\']305 notes = (r"Linear probability models of the indicator that a stress event (VIX $>25$ or a 5-day VIX increase above 20\%) occurs within the next $h$ trading days, "306 r"on standardized predictors observed at $t$; February 2016--December 2025 (the period over which the realized-correlation series is available). "307 r"Coefficients are in percentage points of probability per standard deviation; $t$-statistics in parentheses use Newey--West standard errors with $h$ lags "308 r"because the forward indicator overlaps. The base rate is the unconditional frequency of the indicator. "309 r"\sym{***}, \sym{**}, \sym{*}: $|t|>2.576$, $1.96$, $1.645$.")310 write("tab_rq3", table('@{}lrrrrrr@{}', hdr, rows, notes,311 'Implied versus realized correlation as a predictor of market stress.', 'tab:rq3', size=r'\footnotesize'))312313314def tab_granger():315 g = pd.read_csv(R / "granger_causality.csv")316 lab = {'ATM_IV → RV': r'ATM IV $\to$ RV', 'RV → ATM_IV': r'RV $\to$ ATM IV', 'Skew → RV': r'Skew$_{25\delta}$ $\to$ RV',317 'RV → Skew': r'RV $\to$ Skew$_{25\delta}$', 'ATM_IV → Return': r'ATM IV $\to$ Return', 'Return → ATM_IV': r'Return $\to$ ATM IV',318 'Skew → Return': r'Skew$_{25\delta}$ $\to$ Return', 'Return → Skew': r'Return $\to$ Skew$_{25\delta}$',319 'PC_Ratio → Return': r'PC volume ratio $\to$ Return', 'Return → PC_Ratio': r'Return $\to$ PC volume ratio',320 'Impl_Skew → RV': r'Implied skewness $\to$ RV', 'Impl_Kurt → RV': r'Implied kurtosis $\to$ RV'}321 rows = [f"{lab[r['test']]} & {f(r['avg_f_stat'],2)} & {f(r['avg_p_value'],3)} & {f(r['pct_significant_5pct']*100,1)} & {int(r['n_tickers'])} \\\\"322 for _, r in g.iterrows()]323 hdr = [r'Null hypothesis rejected: $X \to Y$ & Mean $F$ & Mean $p$ & \% tickers sig.\ (5\%) & Tickers \\']324 notes = (r"Per-ticker bivariate Granger causality $F$-tests with 5 lags \citep{granger1969investigating} on winsorized daily series; the table reports the "325 r"cross-ticker mean of the $F$-statistic and of the $p$-value, and the share of tickers rejecting the null of no Granger causality at the 5\% level.")326 write("tab_granger", table('@{}lrrrr@{}', hdr, rows, notes, 'Granger causality tests (5 lags).', 'tab:granger'))327328329def tab_fevd():330 d = pd.read_csv(R / "fevd_results.csv").set_index('horizon')331 v = pd.read_csv(R / "var_results.csv")332 rows = [f"{h} & {f(d.loc[h,'fevd_rv_by_iv_rvfirst'],1)} & {f(d.loc[h,'fevd_rv_by_iv_ivfirst'],1)} & "333 f"{f(d.loc[h,'fevd_iv_by_rv_rvfirst'],1)} & {f(d.loc[h,'fevd_iv_by_rv_ivfirst'],1)} \\\\" for h in [1, 2, 5, 10, 15, 20]]334 hdr = [r'& \multicolumn{2}{c}{RV variance due to IV shock} & \multicolumn{2}{c}{IV variance due to RV shock} \\', r'\cmidrule(lr){2-3}\cmidrule(lr){4-5}',335 r'Horizon (days) & RV first & IV first & RV first & IV first \\']336 notes = (rf"Cholesky-orthogonalized forecast-error variance decomposition of the bivariate VAR(5) in standardized ATM IV and RV, averaged over the "337 rf"{len(v)} tickers with the longest histories (mean $N={num(v['n_obs'].mean())}$ days per ticker). ``RV first'' orders RV before IV (RV shocks may move IV "338 rf"contemporaneously, not the reverse); ``IV first'' is the opposite ordering \citep{{sims1980macroeconomics,kilian2017structural}}. "339 rf"The mean contemporaneous residual correlation is {f(v['resid_corr'].mean(),2)}; mean equation $R^2$: {f(v['r2_iv_eq'].mean(),2)} (IV) and {f(v['r2_rv_eq'].mean(),2)} (RV).")340 write("tab_fevd", table('@{}rrrrr@{}', hdr, rows, notes,341 r'Forecast-error variance decomposition (\%), bivariate VAR(5) in ATM IV and RV.', 'tab:fevd'))342343344def tab_robust_se():345 nw = pd.read_csv(R / "robustness_newey_west.csv")346 dc = pd.read_csv(R / "robustness_double_clustered.csv")347 rows = []348 for hz, pan in [('5-Day', 'A'), ('1-Day', 'B')]:349 a = nw[nw['target'] == hz].set_index('variable'); b = dc[dc['target'] == hz].set_index('variable')350 rows.append(rf"\multicolumn{{6}}{{@{{}}l}}{{\textit{{Panel {pan}: {hz.lower()} return}} ($R^2 = {f(a['r2'].iloc[0],2,True)}\%$, $N = {num(a['n_obs'].iloc[0])}$)}} \\")351 for v in config.RQ1_FEATURES:352 rows.append(f"{VAR_LABEL[v]} & {f(a.loc[v,'coefficient_bps'],2)} & {f(a.loc[v,'hc1_t_stat'],2)} & {f(a.loc[v,'nw_t_stat'],2)} & "353 f"{f(a.loc[v,'nw22_t_stat'],2)} & {f(b.loc[v,'dc_t_stat'],2)}{stars(b.loc[v,'dc_t_stat'])} \\\\")354 if hz == '5-Day':355 rows.append(r'\addlinespace')356 hdr = [r'& & \multicolumn{4}{c}{$t$-statistic} \\', r'\cmidrule(lr){3-6}', r'Predictor & $b$ (bps) & HC1 & NW(5) & NW(22) & DC \\']357 notes = (r"Pooled regressions on the full 69-ticker panel. $b$ is in basis points per standard deviation. HC1: heteroskedasticity-robust \citep{white1980heteroskedasticity}; "358 r"NW($L$): Newey--West with $L$ lags \citep{newey1987simple}; DC: two-way clustering by ticker and date \citep{cameron2011robust,petersen2009estimating}. "359 r"Stars refer to the DC statistic: \sym{***}, \sym{**}, \sym{*}: $|t|>2.576$, $1.96$, $1.645$.")360 write("tab_robust_se", table('@{}lrrrrr@{}', hdr, rows, notes,361 'Robustness of return predictability to the standard-error specification.', 'tab:robust_se', size=r'\footnotesize'))362363364def tab_subperiod():365 d = pd.read_csv(R / "subperiod_results.csv")366 rows = []367 for sp in config.SUBPERIODS:368 s = d[d['subperiod'] == sp]369 r5 = s[s['horizon'] == '5D'].iloc[0]; r1 = s[s['horizon'] == '1D'].iloc[0]370 har = s[s['model'] == 'HAR-RV'].iloc[0]; hiv = s[s['model'] == 'HAR+IV'].iloc[0]371 rows.append(f"{SUBPERIOD_LABEL[sp]} & {num(r5['n_obs'])} & {f(r1['r2'],2,True)} & {f(r5['r2'],2,True)} & {num(har['n_obs'])} & {f(har['r2'],3)} & {f(hiv['r2'],3)} & "372 f"{f((hiv['r2']-har['r2'])/har['r2']*100,1,plus=True)} \\\\")373 hdr = [r'& \multicolumn{3}{c}{Return predictability} & \multicolumn{4}{c}{$RV_{t+1}$ forecasting} \\', r'\cmidrule(lr){2-4}\cmidrule(lr){5-8}',374 r'Subperiod & $N$ & 1-day $R^2$ (\%) & 5-day $R^2$ (\%) & $N$ & HAR-RV & HAR $+$ IV & $\Delta R^2$ (\%) \\']375 notes = (r"Pooled regressions re-estimated within each subperiod (full panel). Return regressions: ten predictors, HC1, pooled winsorization. RV regressions: "376 r"HAR-RV and HAR-RV $+$ IV surface estimated on the same observations (complete IV surface), per-ticker winsorization; $\Delta R^2$ is relative to HAR-RV.")377 write("tab_subperiod", table('@{}lrrrrrrr@{}', hdr, rows, notes, 'Subperiod stability.', 'tab:subperiod', size=r'\footnotesize'))378379380def tab_vix_regime():381 d = pd.read_csv(R / "regime_results.csv")382 rows = []383 for rg in config.VIX_REGIME_LABELS:384 s = d[d['regime'] == rg]385 r1 = s[s['target'] == '1D'].iloc[0]; r5 = s[s['target'] == '5D'].iloc[0]; rv = s[s['target'] == 'RV_1D'].iloc[0]386 rows.append(f"{REGIME_LABEL[rg]} & {f(r1['pct_panel_obs'],1)} & {num(r5['n_obs'])} & {f(r1['r2'],2,True)} & {f(r5['r2'],2,True)} & {f(rv['r2'],3)} \\\\")387 hdr = [r'VIX regime & \% of panel & $N$ & 1-day $R^2$ (\%) & 5-day $R^2$ (\%) & HAR $+$ IV $R^2$ \\']388 notes = (r"Observations are assigned to regimes by the VIX close on day $t$. Return regressions as in Table~\ref{tab:rq1} (full panel); the last column is the "389 r"in-sample $R^2$ of the HAR-RV $+$ IV surface model for $RV_{t+1}$ within the regime.")390 write("tab_vix_regime", table('@{}lrrrrr@{}', hdr, rows, notes, 'Predictability by VIX regime.', 'tab:vix_regime'))391392393def tab_winsor():394 w = pd.read_csv(R / "extended_winsorization_sensitivity.csv")395 e = pd.read_csv(R / "extended_newey_west_lags.csv")396 e = e[e['target'] == '5-Day']397 rows = [rf'\multicolumn{{6}}{{@{{}}l}}{{\textit{{Panel A: winsorization cutoff (5-day return, HC1 $t$-statistics)}}}} \\']398 for _, r in w.iterrows():399 lab = 'None' if pd.isna(r['winsorization']) else r['winsorization'].replace('%', r'\%')400 rows.append(f"{lab} & {f(r['r2'],2,True)} & {f(r['t_implied_skewness'],2)} & {f(r['t_implied_kurtosis'],2)} & {f(r['t_pc_volume_ratio'],2)} & {int(r['n_significant_5pct'])}/10 \\\\")401 rows.append(r'\addlinespace')402 rows.append(rf'\multicolumn{{6}}{{@{{}}l}}{{\textit{{Panel B: Newey--West lag length (5-day return, baseline winsorization)}}}} \\')403 for L in [5, 10, 22]:404 s = e[e['n_lags'] == L].set_index('variable')405 rows.append(f"NW({L}) & & {f(s.loc['implied_skewness','nw_t_stat'],2)} & {f(s.loc['implied_kurtosis_proxy','nw_t_stat'],2)} & "406 f"{f(s.loc['pc_volume_ratio','nw_t_stat'],2)} & {int((s.drop('const')['nw_t_stat'].abs() > 1.96).sum())}/10 \\\\")407 hdr = [r'& $R^2$ (\%) & $t_{\text{impl.\ skew}}$ & $t_{\text{impl.\ kurt}}$ & $t_{\text{PC vol}}$ & \#sig.\ (5\%) \\']408 notes = (r"Pooled 5-day return regression on the full panel (ten predictors). Panel A varies the symmetric winsorization cutoff applied to predictors and target. "409 r"Panel B varies the Newey--West lag length. The last column counts predictors with $|t|>1.96$.")410 write("tab_winsor", table('@{}lrrrrr@{}', hdr, rows, notes, 'Winsorization and HAC-lag sensitivity (5-day returns).', 'tab:winsor'))411412413def tab_ic():414 d = pd.read_csv(R / "extended_information_coefficients.csv")415 d5 = d[d['target'] == '5-Day'].set_index('variable'); d1 = d[d['target'] == '1-Day'].set_index('variable')416 order = d5.sort_values('mean_ic', ascending=False).index417 rows = [f"{VAR_LABEL[v]} & {f(d5.loc[v,'mean_ic'],4)} & {f(d5.loc[v,'t_stat'],2)}{stars(d5.loc[v,'t_stat'])} & {f(d5.loc[v,'pct_positive']*100,1)} & "418 f"{f(d1.loc[v,'mean_ic'],4)} & {f(d1.loc[v,'t_stat'],2)}{stars(d1.loc[v,'t_stat'])} & {f(d1.loc[v,'pct_positive']*100,1)} \\\\" for v in order]419 hdr = [r'& \multicolumn{3}{c}{5-day return} & \multicolumn{3}{c}{1-day return} \\', r'\cmidrule(lr){2-4}\cmidrule(lr){5-7}',420 r'Predictor & Mean IC & $t$ & \% days $>0$ & Mean IC & $t$ & \% days $>0$ \\']421 n = int(d5['n_days'].max())422 notes = (rf"For each trading day with at least 20 stocks, the Spearman rank correlation between the predictor and the forward return is computed across stocks "423 rf"({n:,} days); the table reports the time-series mean, its $t$-statistic (i.i.d.\ standard error, as in version~1.1; with 4 Newey--West lags the 5-day "424 rf"statistics would be roughly halved) and the share of days with a positive IC. \sym{{***}}, \sym{{**}}, \sym{{*}}: $|t|>2.576$, $1.96$, $1.645$.").replace(f"{n:,}", num(n))425 write("tab_ic", table('@{}lrrrrrr@{}', hdr, rows, notes,426 'Daily cross-sectional Spearman information coefficients (individual stocks).', 'tab:ic'))427428429def tab_decile():430 d = pd.read_csv(R / "extended_decile_sorts.csv")431 rows = []432 for sv in d['sort_variable'].unique():433 s = d[d['sort_variable'] == sv]434 for i, (_, r) in enumerate(s.iterrows()):435 scheme = 'Q5$-$Q1' if 'Quintile' in r['scheme'] else 'D10$-$D1'436 rows.append(f"{SORT_LABEL[sv] if i == 0 else ''} & {scheme} & {f(r['annualized_return_pct'],2)} & {f(r['sharpe_ratio'],3)} & {f(r['t_statistic'],2)}{stars(r['t_statistic'])} \\\\")437 rows.append(r'\addlinespace')438 rows = rows[:-1]439 hdr = [r'Sorting variable & Scheme & Ann.\ ret.\ (\%) & Sharpe & $t_{NW}$ \\']440 notes = (r"Equal-weighted daily-rebalanced portfolios of individual stocks, 5-day forward returns, annualization with 52 periods; long-short is the top minus the bottom "441 r"portfolio. $t_{NW}$: Newey--West with 4 lags. \sym{***}, \sym{**}, \sym{*}: $|t|>2.576$, $1.96$, $1.645$.")442 write("tab_decile", table('@{}llrrr@{}', hdr, rows, notes, 'Quintile versus decile long-short sorts (5-day returns).', 'tab:decile'))443444445def tab_loyo():446 d = pd.read_csv(R / "extended_leave_one_year_out.csv")447 rows = [f"{int(r['excluded_year'])} & {f(r['r2_ret_5d'],2,True)} & {f(r['r2_rv_har'],3)} & {f(r['r2_rv_har_iv'],3)} & {f(r['delta_r2_pct'],1,plus=True)} \\\\"448 for _, r in d.iterrows()]449 hdr = [r'Excluded year & 5-day return $R^2$ (\%) & HAR-RV $R^2$ & HAR $+$ IV $R^2$ & $\Delta R^2$ (\%) \\']450 notes = (r"Each row re-estimates the pooled 5-day return regression (ten predictors, HC1) and the two realized-variance models for $RV_{t+1}$ (common sample, "451 r"per-ticker winsorization) on the full panel excluding the indicated calendar year.")452 write("tab_loyo", table('@{}lrrrr@{}', hdr, rows, notes, 'Leave-one-year-out panel $R^2$.', 'tab:loyo'))453454455def tab_quantile():456 q = pd.read_csv(R / "robustness_quantile_regression.csv")457 piv = q.pivot(index='variable', columns='tau', values='coefficient') * 1e4458 rows = [f"{VAR_LABEL[v]} & " + ' & '.join(f(piv.loc[v, t], 1) for t in [0.10, 0.25, 0.50, 0.75, 0.90]) + r' \\'459 for v in config.RQ1_FEATURES]460 hdr = [r'Predictor & $\tau=0.10$ & $\tau=0.25$ & $\tau=0.50$ & $\tau=0.75$ & $\tau=0.90$ \\']461 notes = (r"Quantile regression coefficients \citep{koenker1978regression} of the 5-day return on the standardized predictors (basis points per standard deviation), "462 r"full panel, estimated by iteratively reweighted least squares. No standard errors are computed; the table documents the sign and magnitude of the "463 r"coefficients across the conditional distribution only.")464 write("tab_quantile", table('@{}lrrrrr@{}', hdr, rows, notes, 'Quantile regression coefficients (5-day returns, bps).', 'tab:quantile'))465466467def tab_controls():468 c = pd.read_csv(R / "robustness_with_controls.csv")469 c = c[c['target'] == '5-Day'].set_index('variable')470 base = pd.read_csv(R / "rq1_meta.csv"); base = base[base['label'] == 'All | 5-Day']['r_squared'].dropna().iloc[0]471 rows = [f"{VAR_LABEL[v]} & {f(c.loc[v,'coefficient']*1e4,2)} & {f(c.loc[v,'t_stat'],2)}{stars(c.loc[v,'t_stat'])} \\\\"472 for v in config.RQ1_FEATURES + ['log_volume', 'log_oi', 'abs_return', 'ret_lag1', 'ret_lag5']]473 rows.insert(len(config.RQ1_FEATURES), r'\midrule')474 hdr = [r'Regressor & $b$ (bps) & $t_{HC1}$ \\']475 notes = (rf"Pooled 5-day return regression on the full panel with five additional controls; $R^2 = {f(c['r2'].iloc[0],2,True)}\%$ against {f(base,2,True)}\% without controls. "476 r"\sym{***}, \sym{**}, \sym{*}: $|t|>2.576$, $1.96$, $1.645$.")477 write("tab_controls", table('@{}lrr@{}', hdr, rows, notes, 'Return predictability with additional controls (5-day returns).', 'tab:controls'))478479480def tab_desc_groups():481 g = pd.read_csv(R / "descriptive_by_group.csv").set_index('Group')482 def col(k, d=4, mult=1):483 return ' & '.join(f(g.loc[x, k] * mult, d) for x in ['Stocks', 'ETFs', 'Indices', 'All'])484 rows = [rf"$N$ (ticker-days) & " + ' & '.join(num(g.loc[x, 'N_obs']) for x in ['Stocks', 'ETFs', 'Indices', 'All']) + r' \\',485 rf"$N$ tickers & " + ' & '.join(str(int(g.loc[x, 'N_tickers'])) for x in ['Stocks', 'ETFs', 'Indices', 'All']) + r' \\',486 rf"Mean $IV_{{ATM,30d}}$ & {col('Mean_IV_ATM')} \\", rf"Mean skew$_{{25\delta}}$ & {col('Mean_Skew')} \\",487 rf"Mean PC volume ratio & {col('Mean_PC_ratio', 3)} \\", rf"Mean daily RV ($\times 10^{{4}}$) & {col('Mean_RV', 2, 1e4)} \\",488 rf"Mean daily return (\%) & {col('Mean_Ret_1d', 3, 100)} \\", rf"Std.\ daily return (\%) & {col('Std_Ret_1d', 2, 100)} \\"]489 hdr = [r'& Stocks & ETFs & Indices & All \\']490 notes = (r"Sample period: January 2010--December 2025. $IV_{ATM,30d}$ is the 30-day at-the-money call implied volatility. Skew is the 25-delta put minus 25-delta call IV. "491 r"PC volume ratio is total put volume over total call volume. RV is the realized variance from 5-minute returns; the daily return is the sum of 5-minute log returns.")492 write("tab_desc_groups", table('@{}lrrrr@{}', hdr, rows, notes, 'Descriptive statistics by asset group.', 'tab:desc_groups'))493494495def tab_autocorr():496 a = pd.read_csv(R / "descriptive_autocorrelations.csv").pivot(index='variable', columns='lag', values='avg_autocorr')497 lab = {'iv_atm_30d': 'ATM IV (30d)', 'iv_skew_25d': r'Skew (25$\delta$)', 'rv_daily': 'Realized variance',498 'daily_return': 'Daily return', 'pc_volume_ratio': 'PC volume ratio'}499 rows = [f"{lab[v]} & " + ' & '.join(f(a.loc[v, L], 3) for L in [1, 5, 10, 22]) + r' \\' for v in lab]500 hdr = [r'Variable & Lag 1 & Lag 5 & Lag 10 & Lag 22 \\']501 write("tab_autocorr", table('@{}lrrrr@{}', hdr, rows, 'Cross-sectional average of within-ticker autocorrelations.',502 'Average autocorrelation structure.', 'tab:autocorr'))503504505def tab_app_crises():506 c = pd.read_csv(R / "rq3_crisis_analysis.csv")507 names = {'Flash Crash (2010-05-06)': 'Flash Crash (May 2010)', 'Euro Crisis (2011-08-05)': 'Euro crisis (Aug.\\ 2011)',508 'China Deval (2015-08-24)': 'China devaluation (Aug.\\ 2015)', 'Volmageddon (2018-02-05)': 'Volmageddon (Feb.\\ 2018)',509 'COVID Crash (2020-03-16)': 'COVID-19 (Mar.\\ 2020)', 'Meme Stocks (2021-01-27)': 'Meme stocks (Jan.\\ 2021)',510 'Rate Shock (2022-06-13)': 'Rate shock (June 2022)', 'SVB Crisis (2023-03-10)': 'SVB (Mar.\\ 2023)', 'Aug VIX Spike (2024-08-05)': 'VIX spike (Aug.\\ 2024)'}511 rows = [f"{names[r['crisis']]} & {f(r['pre_impl_corr'],3)} & {f(r['post_impl_corr'],3)} & {f(r['post_impl_corr']-r['pre_impl_corr'],3,plus=True)} & "512 f"{f(r['pre_real_corr'],3)} & {f(r['post_real_corr'],3)} & {f(r['pre_vix'],1)} & {f(r['post_vix'],1)} \\\\" for _, r in c.iterrows()]513 hdr = [r'& \multicolumn{3}{c}{Implied correlation} & \multicolumn{2}{c}{Realized corr.} & \multicolumn{2}{c}{VIX} \\', r'\cmidrule(lr){2-4}\cmidrule(lr){5-6}\cmidrule(lr){7-8}',514 r'Event & Pre & Post & $\Delta$ & Pre & Post & Pre & Post \\']515 notes = (r"Pre: mean over the 30 calendar days before the event date; Post: mean over the 10 calendar days from the event date. The realized-correlation series starts in "516 r"February 2016, hence the missing entries for earlier events.")517 write("tab_app_crises", table('@{}lrrrrrrr@{}', hdr, rows, notes, 'Implied and realized correlation around major market events.', 'tab:app_crises', size=r'\footnotesize'))518519520def tab_rq5():521 d = pd.read_csv(R / "rq5_model_comparison.csv")522 d = d[(d['Target'] == '1-Day RV') & (~d['Model'].str.contains('IV Surface'))]523 lab = {'VIX (benchmark)': r'VIX benchmark, $(VIX_t/100)^2/252$', 'OLS: HAR-RV': 'OLS: HAR-RV', 'Random Forest': 'Random forest (surface $+$ HAR)',524 'Gradient Boosting': 'Gradient boosting (surface $+$ HAR)'}525 rows = [f"{lab[r['Model']]} & {f(r['MSE']*1e7,3)} & {f(r['MAE']*1e4,3)} & {f(r['R²_OOS'],3)} \\\\" for _, r in d.iterrows()]526 hdr = [r'Model & MSE ($\times10^{7}$) & MAE ($\times10^{4}$) & $R^2_{OOS}$ \\']527 notes = (r"Target: next-day 5-minute realized variance of the S\&P~500. Training 2010--2019 (2,515 days), test 2020--2025 (1,507 days). Frozen output of "528 r"\texttt{scripts/06} (raw-dependent; not re-run in the 2026-10-17 revision). The OLS surface models of the original run are omitted because their "529 r"design matrix contained exact linear dependencies (term-structure slopes defined as differences of included tenors) and their forecasts were numerically degenerate; "530 r"the 5- and 22-day targets of the original run were mis-aligned (Section~\ref{sec:data}) and are not reported.")531 write("tab_rq5", table('@{}lrrr@{}', hdr, rows, notes, 'Out-of-sample SPX realized-variance forecasting, 2020--2025 (next-day target).', 'tab:rq5'))532533534def tab_app_importance():535 d = pd.read_csv(R / "rq5_feature_importance.csv")536 lab = {'iv_atm_1w': r'$IV_{ATM,1w}$', 'iv_atm_2w': r'$IV_{ATM,2w}$', 'iv_atm_1m': r'$IV_{ATM,1m}$', 'iv_atm_2m': r'$IV_{ATM,2m}$',537 'iv_atm_3m': r'$IV_{ATM,3m}$', 'iv_atm_6m': r'$IV_{ATM,6m}$', 'total_volume': 'Total option volume', 'pc_vol_ratio': 'PC volume ratio',538 'net_gamma': 'Net gamma exposure', 'rv_w': r'$RV^{(w)}$', 'skew_25d_3m': r'Skew (25$\delta$, 3m)', 'avg_theta': r'Average $\Theta$'}539 top = d.head(10)540 rows = [f"{lab.get(r['feature'], r['feature'].replace('_', r'\_'))} & {f(r['importance'],3)} \\\\" for _, r in top.iterrows()]541 rows += [r'\midrule', rf"\textit{{All other features ({len(d)-10})}} & \textit{{{f(d['importance'].iloc[10:].sum(),3)}}} \\"]542 hdr = [r'Feature & Importance \\']543 write("tab_app_importance", table('@{}lr@{}', hdr, rows, r"Mean-decrease-in-impurity importances of the random forest fitted on 2010--2019 (frozen output of \texttt{scripts/06}).",544 'Random forest feature importance (next-day SPX RV, 2010--2019 training).', 'tab:app_importance'))545546547def tab_rq4():548 d = pd.read_csv(R / "rq4_greeks_decay.csv")549 lab = {'01_1w': '1 week (1--7)', '02_2w': '2 weeks (8--14)', '03_1m': '1 month (15--30)', '04_2m': '2 months (31--60)',550 '05_3m': '3 months (61--90)', '06_6m': '6 months (91--180)'}551 rows = []552 for b in sorted(d['dte_bucket'].unique()):553 rv = d[(d['dte_bucket'] == b) & (d['target'] == 'RV')].iloc[0]; rt = d[(d['dte_bucket'] == b) & (d['target'] == 'Return')].iloc[0]554 rows.append(f"{lab[b]} & {num(rv['n_obs'])} & {f(rv['r_squared'],2,True)} & {f(rv['iv_t'],1)} & {f(rt['r_squared'],2,True)} & {f(rt['iv_t'],2)} \\\\")555 hdr = [r'& & \multicolumn{2}{c}{$RV_{t+1}$} & \multicolumn{2}{c}{$r_{t+1}$} \\', r'\cmidrule(lr){3-4}\cmidrule(lr){5-6}',556 r'DTE bucket (days) & $N$ & $R^2$ (\%) & $t_{IV}$ & $R^2$ (\%) & $t_{IV}$ \\']557 notes = (r"Ten liquid underlyings (AAPL, MSFT, AMZN, GOOG, TSLA, NVDA, META, JPM, SPY, QQQ). Within each days-to-expiry bucket the next-day RV or return is regressed on the "558 r"bucket's average gamma, vega, $|\Theta|$, IV and call/put gamma-weighted open interest (HC1). Only the $R^2$ and the IV $t$-statistic are reported: the Greeks are "559 r"nearly collinear within a bucket and their individual $t$-statistics are not interpretable. Frozen output of \texttt{scripts/05} part~A (raw-dependent).")560 write("tab_rq4", table('@{}lrrrrr@{}', hdr, rows, notes, 'Predictive content of bucket-level Greeks by days to expiry.', 'tab:rq4'))561562563def tab_magnet():564 d = pd.read_csv(R / "rq4_price_magnet.csv")565 lab = {'Q1_Low': 'Q1 (lowest)', 'Q2': 'Q2', 'Q3': 'Q3', 'Q4': 'Q4', 'Q5_High': 'Q5 (highest)'}566 rows = [f"{lab[r['oi_conc_quintile']]} & {num(r['n_obs'])} & {f(r['pct_moved_toward']*100,1)} & {f(r['avg_dist_open']*100,2)} & {f(r['avg_dist_close']*100,2)} \\\\"567 for _, r in d.iterrows()]568 tot = d['n_obs'].sum(); p = (d['pct_moved_toward'] * d['n_obs']).sum() / tot569 se = np.sqrt(0.25 / tot)570 rows += [r'\midrule', rf"All & {num(tot)} & {f(p*100,1)} & & \\"]571 hdr = [r'OI-concentration quintile & $N$ & Moved toward strike (\%) & Open distance (\%) & Close distance (\%) \\']572 notes = (rf"Ticker-days (ten liquid underlyings, 2010--2025) on which the open was within 0.1\%--10\% of the maximum-open-interest strike among options expiring within 30 days. "573 rf"``Moved toward strike'' is the share of days on which the close was relatively closer to that strike than the open; under the null of no attraction the "574 rf"expected share is 50\% (binomial standard error {f(se*100,2)} percentage points for the pooled sample, ignoring dependence). OI concentration is the share of "575 rf"the maximum strike in total open interest. Built from the shipped \texttt{{price\_magnet\_data.parquet}}.")576 write("tab_magnet", table('@{}lrrrr@{}', hdr, rows, notes, 'Price-magnet test: does the underlying drift toward the maximum-open-interest strike?', 'tab:magnet'))577578579def app_tickers():580 tex = lambda xs: ', '.join(x.replace('_', r'\_') for x in xs)581 body = (rf"The 49 stocks are {tex(config.MAJOR_TICKERS)}. "582 rf"The 17 ETFs are {tex(config.ETF_TICKERS)}. The three indices are {tex(config.INDEX_OPTION_TICKERS)}. "583 rf"The 30 constituents used for the implied-correlation index (RQ3) are {tex(config.SPX_CONSTITUENTS)}; "584 rf"the ten underlyings of the Greeks and price-magnet analyses (RQ4) are {tex(config.RQ4_TICKERS)}. "585 r"Lists are taken from \texttt{src/wp7/config.py}; BRK\_B denotes Berkshire Hathaway class B.")586 write("app_tickers", body + "\n")587588589def main():590 OUT.mkdir(parents=True, exist_ok=True)591 app_tickers()592 for fn in [tab_rq1_panel, tab_rq1_coef, tab_sorts, tab_dsort, tab_calendar, tab_rq2_insample, tab_rq2_oos, tab_app_dm,593 tab_rq3, tab_granger, tab_fevd, tab_robust_se, tab_subperiod, tab_vix_regime, tab_winsor, tab_ic, tab_decile,594 tab_loyo, tab_quantile, tab_controls, tab_desc_groups, tab_autocorr, tab_app_crises, tab_rq5, tab_app_importance,595 tab_rq4, tab_magnet]:596 fn()597 print("TABLES COMPLETE.")598599600if __name__ == "__main__":601 main()602