TeX 62.4%
Python 37.5%
1"""2WP 6 -- Intraday periodicity and trading-session analysis (revision, 2026-10-17).34The original 01_analysis.py labelled FirstRate time stamps as UTC, whereas they are in5U.S. Eastern time (the ES series has its daily Globex maintenance break at 17:00 and the6VIX series starts at 03:15 / 09:30 Eastern). This script rebuilds the intraday analysis7with an explicit America/New_York -> UTC conversion from one-minute FirstRate bars held8locally (sibling repository WP 12, same vendor; path in WP6_MINUTE_DIR), aggregated to the9paper's 5-minute grid over the common sample 2016-03-11 to 2026-03-11:1011 * hourly profile of average absolute 5-minute returns, in UTC, for ES, NQ, VIX, BTC;12 * session realized variances (Asian 00:00-08:00, European 08:00-14:30,13 U.S. 14:30-21:00 UTC) and their shares of the daily total;14 * session-ordered predictive regressions: does Bitcoin's Asian/European-session15 volatility help predict the U.S.-session volatility of ES (and vice versa), and16 does weekend Bitcoin volatility predict Monday equity volatility?1718Outputs: results/rev_intraday_*.csv, paper/tables/tab_sessions*.tex, figures/fig5_hourly_volatility.png.19Run time: about 2 minutes (loading 1-minute parquet files).20"""21import os22import time23import warnings24import numpy as np25import pandas as pd26import matplotlib27matplotlib.use('Agg')28import matplotlib.pyplot as plt29import statsmodels.api as sm3031warnings.filterwarnings('ignore')32T0 = time.time()33ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))34RESULTS = os.path.join(ROOT, 'results')35FIGURES = os.path.join(ROOT, 'figures')36TABLES = os.path.join(ROOT, 'paper', 'tables')37PAPER_FIG = os.path.join(ROOT, 'paper', 'figures')38MINUTE_DIR = os.environ.get('WP6_MINUTE_DIR', os.path.join(ROOT, '..', 'realized-volatility-horse-race', 'data', 'raw'))39FILES = {'ES': 'futures_ES_1min.parquet', 'NQ': 'futures_NQ_1min.parquet', 'VIX': 'index_VIX_1min.parquet',40 'BTC': 'crypto_BTC_1min.parquet'}41START, END = '2016-03-11', '2026-03-11'42COLORS = {'ES': '#1f77b4', 'NQ': '#aec7e8', 'VIX': '#8c564b', 'BTC': '#9467bd'}43SESSIONS = {'Asian': (0.0, 8.0), 'European': (8.0, 14.5), 'US': (14.5, 21.0)}444546def log(m):47 print(f"[{time.time() - T0:6.1f}s] {m}", flush=True)484950def load_5min_utc(sym):51 df = pd.read_parquet(os.path.join(MINUTE_DIR, FILES[sym]), columns=['datetime', 'close'])52 s = df.set_index('datetime')['close'].sort_index()53 s = s.loc[START:END]54 # FirstRate stamps are naive U.S. Eastern (bar-start convention); convert to UTC55 idx = s.index.tz_localize('America/New_York', ambiguous='NaT', nonexistent='NaT')56 s = s[~idx.isna()]57 s.index = idx[~idx.isna()].tz_convert('UTC')58 s5 = s.resample('5min').last().dropna()59 r = np.log(s5).diff().dropna()60 # drop returns spanning a gap longer than one hour (session breaks, weekends)61 gaps = s5.index.to_series().diff().dt.total_seconds().reindex(r.index)62 r = r[gaps <= 3600]63 return r646566rets = {}67for sym in FILES:68 rets[sym] = load_5min_utc(sym)69 log(f"{sym}: {len(rets[sym]):,} five-minute returns, {rets[sym].index.min().date()} -> {rets[sym].index.max().date()}")7071# ─── 1. Hourly profile (UTC) ──────────────────────────────────────────────────72hourly = pd.DataFrame({sym: r.abs().groupby(r.index.hour).mean() * 100 for sym, r in rets.items()})73hourly.index.name = 'hour_utc'74hourly.to_csv(os.path.join(RESULTS, 'rev_intraday_hourly_absret_utc.csv'), float_format='%.5f')75hourly_et = pd.DataFrame({sym: r.abs().groupby(r.index.tz_convert('America/New_York').hour).mean() * 10076 for sym, r in rets.items()})77hourly_et.index.name = 'hour_eastern'78hourly_et.to_csv(os.path.join(RESULTS, 'rev_intraday_hourly_absret_eastern.csv'), float_format='%.5f')7980fig, axes = plt.subplots(2, 2, figsize=(13, 7.5))81for ax, sym in zip(axes.flatten(), ['VIX', 'ES', 'NQ', 'BTC']):82 h = hourly[sym]83 ax.bar(h.index, h.values, color=COLORS[sym], alpha=0.85, edgecolor='white', linewidth=0.5)84 ax.axvspan(-0.5, 7.5, alpha=0.07, color='blue', label='Asian 00:00–08:00 UTC')85 ax.axvspan(7.5, 14.0, alpha=0.07, color='green', label='European 08:00–14:30 UTC')86 ax.axvspan(14.0, 20.5, alpha=0.07, color='orange', label='U.S. 14:30–21:00 UTC')87 ax.set_title(f'{sym} — mean |5-min return| by hour (UTC)')88 ax.set_xlabel('Hour of day (UTC)'); ax.set_ylabel('Mean |return| (%)')89 ax.set_xticks(range(0, 24, 2)); ax.set_xlim(-0.6, 23.6)90 if sym == 'VIX':91 ax.legend(fontsize=7.5, loc='upper right')92plt.suptitle('Intraday periodicity of 5-minute absolute returns, 2016-03 to 2026-03 (time stamps converted from U.S. Eastern)',93 fontsize=11)94plt.tight_layout()95for d in (FIGURES, PAPER_FIG):96 fig.savefig(os.path.join(d, 'fig5_hourly_volatility.png'), dpi=150, bbox_inches='tight')97plt.close(fig)98log("hourly figure written")99100# ─── 2. Session realized variances per UTC day ────────────────────────────────101def session_rv(r):102 hour = r.index.hour + r.index.minute / 60103 day = r.index.tz_convert('UTC').normalize().tz_localize(None)104 out = {}105 for name, (a, b) in SESSIONS.items():106 m = (hour >= a) & (hour < b)107 rv = (r[m] ** 2).groupby(day[m]).sum()108 n = r[m].groupby(day[m]).size()109 rv = rv[n >= 12] # at least one hour of bars110 out[name] = rv111 out['Day'] = (r ** 2).groupby(day).sum()112 return pd.DataFrame(out)113114115sess = {sym: session_rv(r) for sym, r in rets.items()}116rows = []117for sym, df in sess.items():118 wd = df[df.index.dayofweek < 5]119 for name in SESSIONS:120 s = wd[name].dropna()121 rows.append({'Asset': sym, 'Session': name, 'Days': len(s), 'Ann_vol_pct': 100 * np.sqrt(252 * s).mean(),122 'Mean_absret_pct': np.nan,123 'Share_of_day_pct': 100 * (wd[name] / wd[['Asian', 'European', 'US']].sum(axis=1)).mean()})124sess_tbl = pd.DataFrame(rows)125# mean absolute return per session126for sym, r in rets.items():127 hour = r.index.hour + r.index.minute / 60128 for name, (a, b) in SESSIONS.items():129 m = (hour >= a) & (hour < b) & (r.index.dayofweek < 5)130 sess_tbl.loc[(sess_tbl.Asset == sym) & (sess_tbl.Session == name), 'Mean_absret_pct'] = 100 * r[m].abs().mean()131sess_tbl.to_csv(os.path.join(RESULTS, 'rev_intraday_sessions.csv'), index=False, float_format='%.4f')132log("session table done")133134# weekend share for BTC: average over ISO weeks of (Saturday + Sunday variance) / (weekly variance),135# which is robust to the extreme 2017 days; day-of-week profile as mean annualized daily RV136btc_day = sess['BTC']['Day']137wk_var = btc_day.groupby(btc_day.index.to_period('W')).sum()138wk_end = btc_day[btc_day.index.dayofweek >= 5].groupby(btc_day[btc_day.index.dayofweek >= 5].index.to_period('W')).sum()139ratio = (wk_end / wk_var).dropna()140ratio = ratio[btc_day.groupby(btc_day.index.to_period('W')).size() == 7]141wk_share = ratio.mean()142btc_rv_day = 100 * np.sqrt(252 * btc_day)143btc_dow = btc_rv_day.groupby(btc_rv_day.index.dayofweek).mean()144btc_dow.to_csv(os.path.join(RESULTS, 'rev_intraday_btc_dayofweek.csv'), header=['ann_vol_pct'])145log(f"BTC weekend share of weekly variance: {100 * wk_share:.1f}%")146147# ─── 3. Session-ordered predictive regressions ───────────────────────────────148def predictive(dep, pred_same_day=None, controls_lag=4, label=''):149 """OLS of log dep RV on its own lags (+ optional same-day predictor and its lags); HAC Wald on predictor block."""150 df = pd.DataFrame({'y': np.log(dep.where(dep > 0))})151 X = [df['y'].shift(l).rename(f'y_l{l}') for l in range(1, controls_lag + 1)]152 pred_cols = []153 if pred_same_day is not None:154 x = np.log(pred_same_day.where(pred_same_day > 0)).reindex(df.index)155 for l in range(0, controls_lag):156 X.append(x.shift(l).rename(f'x_l{l}'))157 pred_cols.append(f'x_l{l}')158 X = sm.add_constant(pd.concat(X, axis=1))159 data = pd.concat([df['y'], X], axis=1).dropna()160 m = sm.OLS(data['y'], data.drop(columns='y')).fit(cov_type='HAC', cov_kwds={'maxlags': 10})161 R = np.zeros((len(pred_cols), m.params.shape[0]))162 for k, c in enumerate(pred_cols):163 R[k, list(m.params.index).index(c)] = 1164 w = m.wald_test(R, use_f=False, scalar=True)165 base = sm.OLS(data['y'], data.drop(columns=['y'] + pred_cols)).fit()166 return {'Test': label, 'T': int(m.nobs), 'coef_x_l0': float(m.params['x_l0']), 'se_x_l0': float(m.bse['x_l0']),167 'Wald_HAC': float(w.statistic), 'p': float(w.pvalue), 'R2_with': m.rsquared, 'R2_without': base.rsquared}168169170wd = {sym: df[df.index.dayofweek < 5] for sym, df in sess.items()}171tests = []172# (i) Bitcoin earlier sessions -> ES / NQ / VIX U.S. session (same day)173for tgt in ['ES', 'NQ', 'VIX']:174 for s in ['Asian', 'European']:175 tests.append(predictive(wd[tgt]['US'].dropna(), wd['BTC'][s], label=f'BTC {s} (t) $\\rightarrow$ {tgt} U.S. (t)'))176# (ii) ES / VIX U.S. session (t-1) -> Bitcoin Asian session (t): shift predictor by one day177for src in ['ES', 'NQ', 'VIX']:178 prev = wd[src]['US'].dropna().copy()179 prev.index = prev.index + pd.offsets.BDay(1)180 tests.append(predictive(wd['BTC']['Asian'].dropna(), prev, label=f'{src} U.S. (t$-$1) $\\rightarrow$ BTC Asian (t)'))181# (iii) ES European (t) -> BTC U.S. (t) and BTC U.S. (t) -> ES U.S.(t) skipped: contemporaneous overlap182# (iv) ES Asian (Globex overnight, t) -> BTC U.S. (t) for symmetry183tests.append(predictive(wd['BTC']['US'].dropna(), wd['ES']['Asian'], label='ES Asian (t) $\\rightarrow$ BTC U.S. (t)'))184tests.append(predictive(wd['ES']['US'].dropna(), wd['ES']['Asian'], label='ES Asian (t) $\\rightarrow$ ES U.S. (t) [own overnight]'))185pred_df = pd.DataFrame(tests)186pred_df.to_csv(os.path.join(RESULTS, 'rev_intraday_session_predictive.csv'), index=False, float_format='%.4f')187log("session predictive regressions done")188189# (v) weekend Bitcoin variance -> Monday ES U.S.-session variance, controlling for Friday ES and Friday BTC190es_us = wd['ES']['US'].dropna()191mon = es_us[es_us.index.dayofweek == 0]192rows = []193for d in mon.index:194 fri = d - pd.Timedelta(days=3)195 sat, sun = d - pd.Timedelta(days=2), d - pd.Timedelta(days=1)196 bwk = btc_day.get(sat, np.nan) + btc_day.get(sun, np.nan)197 rows.append({'date': d, 'es_mon': mon[d], 'es_fri': es_us.get(fri, np.nan), 'btc_wk': bwk,198 'btc_fri': btc_day.get(fri, np.nan), 'btc_mon_asia': wd['BTC']['Asian'].get(d, np.nan)})199wk = pd.DataFrame(rows).set_index('date')200wk = np.log(wk.where(wk > 0)).dropna()201m1 = sm.OLS(wk['es_mon'], sm.add_constant(wk[['es_fri', 'btc_fri']])).fit(cov_type='HAC', cov_kwds={'maxlags': 4})202m2 = sm.OLS(wk['es_mon'], sm.add_constant(wk[['es_fri', 'btc_fri', 'btc_wk']])).fit(cov_type='HAC', cov_kwds={'maxlags': 4})203m3 = sm.OLS(wk['es_mon'], sm.add_constant(wk[['es_fri', 'btc_fri', 'btc_wk', 'btc_mon_asia']])).fit(cov_type='HAC', cov_kwds={'maxlags': 4})204wk_res = pd.DataFrame({205 'model': ['Friday controls only', '+ weekend BTC', '+ weekend BTC + Monday Asian BTC'],206 'N': [int(m1.nobs), int(m2.nobs), int(m3.nobs)],207 'coef_btc_wk': [np.nan, m2.params['btc_wk'], m3.params['btc_wk']],208 'se_btc_wk': [np.nan, m2.bse['btc_wk'], m3.bse['btc_wk']],209 'p_btc_wk': [np.nan, m2.pvalues['btc_wk'], m3.pvalues['btc_wk']],210 'coef_btc_mon_asia': [np.nan, np.nan, m3.params['btc_mon_asia']],211 'p_btc_mon_asia': [np.nan, np.nan, m3.pvalues['btc_mon_asia']],212 'coef_es_fri': [m1.params['es_fri'], m2.params['es_fri'], m3.params['es_fri']],213 'R2': [m1.rsquared, m2.rsquared, m3.rsquared]})214wk_res.to_csv(os.path.join(RESULTS, 'rev_intraday_weekend_btc_monday_es.csv'), index=False, float_format='%.4f')215log("weekend regression done")216217# ─── 4. LaTeX tables ──────────────────────────────────────────────────────────218def fmt(x, d=2, sign=False):219 if pd.isna(x):220 return '--'221 s = f"{x:+.{d}f}" if sign else f"{x:.{d}f}"222 return s.replace('-', '$-$').replace('+', '$+$') if sign else s223224225def stars(p):226 return r'$^{***}$' if p < 0.01 else r'$^{**}$' if p < 0.05 else r'$^{*}$' if p < 0.10 else ''227228229L = [r"\begin{tabular}{@{}llrrrr@{}}", r"\toprule",230 r"Asset & Session (UTC) & Days & Mean $|r|$ (\%) & Annualized vol.\ (\%) & Share of day (\%) \\", r"\midrule"]231for sym in ['ES', 'NQ', 'VIX', 'BTC']:232 sub = sess_tbl[sess_tbl.Asset == sym]233 for k, (_, r) in enumerate(sub.iterrows()):234 lab = {'Asian': 'Asian 00:00--08:00', 'European': 'European 08:00--14:30', 'US': 'U.S. 14:30--21:00'}[r['Session']]235 L.append(f"{sym if k == 0 else ''} & {lab} & {int(r['Days'])} & {fmt(r['Mean_absret_pct'], 3)} & "236 f"{fmt(r['Ann_vol_pct'], 1)} & {fmt(r['Share_of_day_pct'], 1)} \\\\")237 if sym != 'BTC':238 L.append(r"\addlinespace")239L += [r"\bottomrule", r"\end{tabular}"]240with open(os.path.join(TABLES, 'tab_sessions.tex'), 'w') as f:241 f.write("\n".join(L))242243L = [r"\begin{tabular}{@{}lrrrrr@{}}", r"\toprule",244 r"Predictive relation & $T$ & Coef.\ (same day) & Wald (HAC) & $p$ & $\Delta R^2$ \\", r"\midrule"]245for _, r in pred_df.iterrows():246 L.append(f"{r['Test']} & {int(r['T'])} & {fmt(r['coef_x_l0'], 3, True)} ({fmt(r['se_x_l0'], 3)}) & "247 f"{fmt(r['Wald_HAC'], 2)}{stars(r['p'])} & {fmt(r['p'], 3)} & {fmt(r['R2_with'] - r['R2_without'], 3)} \\\\")248L += [r"\bottomrule", r"\end{tabular}"]249with open(os.path.join(TABLES, 'tab_sessions_predictive.tex'), 'w') as f:250 f.write("\n".join(L))251252L = [r"\begin{tabular}{@{}lrrrrr@{}}", r"\toprule",253 r"Specification & $N$ & $\ln \mathrm{RV}^{\mathrm{ES}}_{\mathrm{Fri}}$ & $\ln \mathrm{RV}^{\mathrm{BTC}}_{\mathrm{weekend}}$ & $\ln \mathrm{RV}^{\mathrm{BTC}}_{\mathrm{Mon,\,Asian}}$ & $R^2$ \\",254 r"\midrule"]255for _, r in wk_res.iterrows():256 bw = '--' if pd.isna(r['coef_btc_wk']) else f"{fmt(r['coef_btc_wk'], 3, True)}{stars(r['p_btc_wk'])} ({fmt(r['se_btc_wk'], 3)})"257 ba = '--' if pd.isna(r['coef_btc_mon_asia']) else f"{fmt(r['coef_btc_mon_asia'], 3, True)}{stars(r['p_btc_mon_asia'])}"258 L.append(f"{r['model']} & {int(r['N'])} & {fmt(r['coef_es_fri'], 3, True)} & {bw} & {ba} & {fmt(r['R2'], 3)} \\\\")259L += [r"\bottomrule", r"\end{tabular}"]260with open(os.path.join(TABLES, 'tab_weekend_btc.tex'), 'w') as f:261 f.write("\n".join(L))262263key = {'btc_weekend_share_pct': 100 * wk_share}264for sym in FILES:265 h = hourly[sym]266 key[f'{sym}_peak_hour_utc'] = int(h.idxmax()); key[f'{sym}_peak_absret'] = h.max()267 key[f'{sym}_min_hour_utc'] = int(h.idxmin()); key[f'{sym}_min_absret'] = h.min()268 key[f'{sym}_peak_to_trough'] = h.max() / h.min()269for d, v in btc_dow.items():270 key[f'BTC_dow{d}_annvol'] = v271pd.Series(key).to_csv(os.path.join(RESULTS, 'rev_intraday_key_numbers.csv'), header=['value'])272print(pd.Series(key).to_string())273print(sess_tbl.to_string())274print(pred_df.to_string())275print(wk_res.to_string())276log("done")277