TeX 62.4%
Python 37.5%
1"""2WP 6 -- Revision of the connectedness analysis (editorial-board review, 2026-10-17).34Starts from the daily realized-volatility panel produced by 01_analysis.py5(results/daily_rv_all.csv, annualized realized volatility per instrument and calendar day)6and re-estimates every VAR/GFEVD object on a correctly aligned panel:78 * trading-day panel: a day enters only if all ten instruments have an own realized9 volatility (no forward filling; the original script forward-filled weekends and10 holidays, so that 1,137 of its 3,651 rows were stale copies);11 * log realized volatility (Andersen et al., 2003) instead of levels;12 * VAR(4) -- HQ criterion -- with H = 10, as in the original paper;13 * residual-bootstrap confidence intervals for TSI, TO, FROM, NET (Kilian, 1998);14 * robustness grid (transformation, p, H, sample, weekend-inclusive crypto RV, subsets);15 * rolling windows of 200 / 250 / 500 trading days; calendar-year sub-samples;16 * regime sub-samples (pre-COVID, COVID crash, 2022 hiking cycle, spot-ETF era);17 * Granger-causality tests with conventional F and HAC-robust Wald statistics;18 * replication of the original (flawed) construction for the audit trail.1920Outputs: results/rev_*.csv, paper/tables/*.tex (generated, not typed), figures/fig*.png.21Run time: about 3 minutes on a laptop (bootstrap included).22"""23import os24import time25import warnings26import numpy as np27import pandas as pd28import matplotlib29matplotlib.use('Agg')30import matplotlib.pyplot as plt31import matplotlib.patches as mpatches32import networkx as nx33from statsmodels.tsa.api import VAR34from statsmodels.tsa.stattools import adfuller, kpss, acf35import statsmodels.api as sm3637warnings.filterwarnings('ignore')38T0 = time.time()3940ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))41RESULTS = os.path.join(ROOT, 'results')42FIGURES = os.path.join(ROOT, 'figures')43TABLES = os.path.join(ROOT, 'paper', 'tables')44PAPER_FIG = os.path.join(ROOT, 'paper', 'figures')45for d in (RESULTS, FIGURES, TABLES, PAPER_FIG):46 os.makedirs(d, exist_ok=True)4748ASSETS = ['ES', 'NQ', 'EURUSD', 'USDJPY', 'GBPUSD', 'GC', 'CL', 'BTC', 'ETH', 'VIX']49ASSET_CLASS = {'ES': 'Equity', 'NQ': 'Equity', 'EURUSD': 'FX', 'USDJPY': 'FX', 'GBPUSD': 'FX',50 'GC': 'Commodity', 'CL': 'Commodity', 'BTC': 'Crypto', 'ETH': 'Crypto', 'VIX': 'Volatility'}51CLASS_COLOR = {'Equity': '#1f77b4', 'FX': '#2ca02c', 'Commodity': '#ff7f0e', 'Crypto': '#9467bd',52 'Volatility': '#8c564b'}53P_BASE, H_BASE, W_BASE = 4, 10, 20054B_BOOT = 50055RNG = np.random.default_rng(20261017)56SAMPLE_START, SAMPLE_END = '2016-03-11', '2026-03-11'5758plt.rcParams.update({'font.size': 10, 'axes.titlesize': 11, 'axes.labelsize': 10, 'figure.dpi': 150,59 'savefig.dpi': 150, 'font.family': 'DejaVu Sans',60 'axes.spines.top': False, 'axes.spines.right': False})616263def log(msg):64 print(f"[{time.time() - T0:6.1f}s] {msg}", flush=True)656667def savefig(fig, name):68 for d in (FIGURES, PAPER_FIG):69 fig.savefig(os.path.join(d, name), dpi=150, bbox_inches='tight')70 plt.close(fig)717273# ─── 1. Panels ────────────────────────────────────────────────────────────────74rv_all = pd.read_csv(os.path.join(RESULTS, 'daily_rv_all.csv'), parse_dates=['date']).set_index('date')75rv_all = rv_all[ASSETS]7677# (a) trading-day panel: all ten instruments observed (VIX trades only on CBOE days)78panel = rv_all.loc[SAMPLE_START:SAMPLE_END].dropna()79log(f"trading-day panel: {panel.index[0].date()} -> {panel.index[-1].date()}, T = {len(panel)}; "80 f"weekday counts = {panel.index.dayofweek.value_counts().sort_index().to_dict()}")8182# (b) replication of the original construction (calendar days, forward fill up to 5 days, levels)83orig_panel = rv_all.ffill(limit=5).dropna(thresh=len(ASSETS) - 2).dropna()84n_weekend = int((orig_panel.index.dayofweek >= 5).sum())85n_ffilled_rows = int((orig_panel.index.isin(rv_all.dropna().index) == False).sum())86log(f"original construction: T = {len(orig_panel)}, weekend rows = {n_weekend}, "87 f"rows with at least one forward-filled value = {n_ffilled_rows}")8889# (c) weekend-inclusive crypto RV: average daily variance over the calendar days since the90# previous trading day (so Monday carries Saturday, Sunday and Monday variance)91def weekend_inclusive(series, trading_days):92 var = series.dropna() ** 293 out = {}94 prev = None95 for d in trading_days:96 if prev is None:97 out[d] = var.get(d, np.nan)98 else:99 win = var.loc[prev + pd.Timedelta(days=1):d]100 out[d] = win.mean() if len(win) else np.nan101 prev = d102 return pd.Series(out) ** 0.5103104105panel_wk = panel.copy()106for c in ('BTC', 'ETH'):107 panel_wk[c] = weekend_inclusive(rv_all[c], panel.index)108panel_wk = panel_wk.dropna()109110111# ─── 2. Descriptive statistics of daily RV ────────────────────────────────────112rows = []113for a in ASSETS:114 s = panel[a]115 ls = np.log(s)116 adf_l = adfuller(ls, regression='c', autolag='AIC')117 kp_l = kpss(ls, regression='c', nlags='auto')118 rows.append({'Asset': a, 'Class': ASSET_CLASS[a], 'Mean_RV_pct': 100 * s.mean(), 'Median_RV_pct': 100 * s.median(),119 'Max_RV_pct': 100 * s.max(), 'Skew_level': s.skew(), 'Kurt_level': s.kurt(),120 'Skew_log': ls.skew(), 'Kurt_log': ls.kurt(), 'AC1_log': acf(ls, nlags=22)[1],121 'AC22_log': acf(ls, nlags=22)[22], 'ADF_log': adf_l[0], 'ADF_p': adf_l[1],122 'KPSS_log': kp_l[0], 'KPSS_p': kp_l[1]})123desc = pd.DataFrame(rows).set_index('Asset')124desc.to_csv(os.path.join(RESULTS, 'rev_table_rv_descriptives.csv'), float_format='%.4f')125log("descriptives done")126127128# ─── 3. GFEVD machinery ───────────────────────────────────────────────────────129def gfevd_from_var(coefs, sigma, H):130 """Normalized generalized FEVD (Koop et al., 1996; Pesaran and Shin, 1998)."""131 n = sigma.shape[0]132 p = coefs.shape[0]133 psi = [np.eye(n)]134 for h in range(1, H):135 acc = np.zeros((n, n))136 for j in range(1, min(h, p) + 1):137 acc += coefs[j - 1] @ psi[h - j]138 psi.append(acc)139 num = np.zeros((n, n))140 den = np.zeros(n)141 for h in range(H):142 ps = psi[h] @ sigma143 num += ps ** 2144 den += np.einsum('ij,ij->i', ps, psi[h])145 theta = num / (np.diag(sigma)[None, :] * den[:, None])146 return theta / theta.sum(axis=1, keepdims=True)147148149def connectedness(y, p=P_BASE, H=H_BASE):150 """Fit VAR(p) with constant, return dict of connectedness objects (in percent)."""151 res = VAR(y.values).fit(p, trend='c')152 th = gfevd_from_var(res.coefs, res.sigma_u, H) * 100153 off = th.copy()154 np.fill_diagonal(off, 0)155 out = {'theta': th, 'FROM': off.sum(axis=1), 'TO': off.sum(axis=0), 'TSI': off.sum() / th.shape[0],156 'res': res}157 out['NET'] = out['TO'] - out['FROM']158 out['NETPAIR'] = th.T - th # [i,j] = theta_ji - theta_ij : net from i to j159 return out160161162def bootstrap_connectedness(y, p=P_BASE, H=H_BASE, B=B_BOOT, rng=RNG):163 """Residual (i.i.d.) bootstrap of the VAR; returns arrays of TSI, TO, FROM, NET."""164 base = connectedness(y, p, H)165 res = base['res']166 T, n = y.shape167 u = res.resid - res.resid.mean(axis=0)168 A = res.coefs169 c = res.intercept170 Y = y.values171 tsi_b, to_b, from_b, net_b = [], [], [], []172 for b in range(B):173 idx = rng.integers(0, len(u), len(u))174 ub = u[idx]175 yb = np.empty_like(Y)176 yb[:p] = Y[:p]177 for t in range(p, T):178 acc = c.copy()179 for l in range(p):180 acc += A[l] @ yb[t - l - 1]181 yb[t] = acc + ub[t - p]182 rb = VAR(yb).fit(p, trend='c')183 th = gfevd_from_var(rb.coefs, rb.sigma_u, H) * 100184 off = th.copy()185 np.fill_diagonal(off, 0)186 fr, to = off.sum(axis=1), off.sum(axis=0)187 tsi_b.append(off.sum() / n)188 to_b.append(to)189 from_b.append(fr)190 net_b.append(to - fr)191 return base, np.array(tsi_b), np.array(to_b), np.array(from_b), np.array(net_b)192193194def ci(arr, q=(5, 95)):195 return np.percentile(arr, q, axis=0)196197198# ─── 4. Full-sample connectedness (log RV, trading days) with bootstrap ───────199ylog = np.log(panel)200log("full-sample bootstrap ...")201base, tsi_b, to_b, from_b, net_b = bootstrap_connectedness(ylog)202tsi_ci, to_ci, from_ci, net_ci = ci(tsi_b), ci(to_b), ci(from_b), ci(net_b)203log(f"full-sample TSI = {base['TSI']:.2f} [{tsi_ci[0]:.2f}, {tsi_ci[1]:.2f}]")204205full_tbl = pd.DataFrame(base['theta'], index=ASSETS, columns=ASSETS)206full_tbl['FROM'] = base['FROM']207full_tbl.loc['TO'] = list(base['TO']) + [base['TSI']]208full_tbl.loc['NET'] = list(base['NET']) + [np.nan]209full_tbl.to_csv(os.path.join(RESULTS, 'rev_table_connectedness_full.csv'), float_format='%.2f')210ci_tbl = pd.DataFrame({'TO': base['TO'], 'TO_lo': to_ci[0], 'TO_hi': to_ci[1], 'FROM': base['FROM'],211 'FROM_lo': from_ci[0], 'FROM_hi': from_ci[1], 'NET': base['NET'], 'NET_lo': net_ci[0],212 'NET_hi': net_ci[1]}, index=ASSETS)213ci_tbl.loc['TSI', ['NET', 'NET_lo', 'NET_hi']] = [base['TSI'], tsi_ci[0], tsi_ci[1]]214ci_tbl.to_csv(os.path.join(RESULTS, 'rev_table_connectedness_full_ci.csv'), float_format='%.2f')215216# original construction replicated (audit) and intermediate variants217orig = connectedness(orig_panel, P_BASE, H_BASE)218lvl_td = connectedness(panel, P_BASE, H_BASE)219audit = pd.DataFrame([220 {'Construction': 'Original: calendar days, forward-filled, RV levels', 'T': len(orig_panel), 'TSI': orig['TSI'],221 'NET_ES': orig['NET'][0], 'NET_NQ': orig['NET'][1], 'NET_VIX': orig['NET'][9], 'NET_BTC': orig['NET'][7],222 'NET_USDJPY': orig['NET'][3]},223 {'Construction': 'Trading days, RV levels', 'T': len(panel), 'TSI': lvl_td['TSI'], 'NET_ES': lvl_td['NET'][0],224 'NET_NQ': lvl_td['NET'][1], 'NET_VIX': lvl_td['NET'][9], 'NET_BTC': lvl_td['NET'][7],225 'NET_USDJPY': lvl_td['NET'][3]},226 {'Construction': 'Trading days, log RV (revised baseline)', 'T': len(panel), 'TSI': base['TSI'],227 'NET_ES': base['NET'][0], 'NET_NQ': base['NET'][1], 'NET_VIX': base['NET'][9], 'NET_BTC': base['NET'][7],228 'NET_USDJPY': base['NET'][3]},229])230audit.to_csv(os.path.join(RESULTS, 'rev_audit_construction.csv'), index=False, float_format='%.2f')231log("audit table done")232233# ─── 5. Robustness grid ───────────────────────────────────────────────────────234rob_rows = []235236237def rob(label, y, p=P_BASE, H=H_BASE):238 r = connectedness(y, p, H)239 d = {'Specification': label, 'T': len(y), 'p': p, 'H': H, 'TSI': r['TSI']}240 for a in ['ES', 'NQ', 'USDJPY', 'GC', 'CL', 'BTC', 'ETH', 'VIX']:241 if a in y.columns:242 d[f'NET_{a}'] = r['NET'][list(y.columns).index(a)]243 else:244 d[f'NET_{a}'] = np.nan245 if 'BTC' in y.columns:246 d['Own_BTC'] = r['theta'][list(y.columns).index('BTC'), list(y.columns).index('BTC')]247 rob_rows.append(d)248 return r249250251rob('Baseline: log RV, trading days, VAR(4), H = 10', ylog)252rob('RV in levels (trading days)', panel)253for p in (1, 2, 3, 6):254 rob(f'VAR({p})', ylog, p=p)255for H in (5, 20):256 rob(f'H = {H}', ylog, H=H)257rob('Weekend-inclusive crypto RV', np.log(panel_wk))258rob('Sample 2018-01 to 2026-03', ylog.loc['2018-01-01':])259rob('Sample 2016-03 to 2019-12', ylog.loc[:'2019-12-31'])260rob('Sample 2020-01 to 2026-03', ylog.loc['2020-01-01':])261rob('Excluding VIX', ylog.drop(columns=['VIX']))262rob('Excluding BTC and ETH', ylog.drop(columns=['BTC', 'ETH']))263rob('Excluding ETH (BTC only)', ylog.drop(columns=['ETH']))264rob_df = pd.DataFrame(rob_rows)265rob_df.to_csv(os.path.join(RESULTS, 'rev_table_robustness.csv'), index=False, float_format='%.2f')266log("robustness grid done")267268# ─── 6. Rolling windows ───────────────────────────────────────────────────────269def rolling(y, w, p=P_BASE, H=H_BASE):270 dates, tsi, net = [], [], []271 Y = y.values272 for i in range(len(y) - w + 1):273 r = VAR(Y[i:i + w]).fit(p, trend='c')274 th = gfevd_from_var(r.coefs, r.sigma_u, H) * 100275 off = th.copy()276 np.fill_diagonal(off, 0)277 dates.append(y.index[i + w - 1])278 tsi.append(off.sum() / th.shape[0])279 net.append(off.sum(axis=0) - off.sum(axis=1))280 out = pd.DataFrame(net, index=dates, columns=[f'NET_{a}' for a in y.columns])281 out.insert(0, 'TSI', tsi)282 out.index.name = 'date'283 return out284285286roll = {}287for w in (200, 250, 500):288 roll[w] = rolling(ylog, w)289 roll[w].to_csv(os.path.join(RESULTS, f'rev_rolling_w{w}.csv'), float_format='%.3f')290 log(f"rolling w={w}: {len(roll[w])} windows, TSI mean {roll[w]['TSI'].mean():.2f}, "291 f"min {roll[w]['TSI'].min():.2f} ({roll[w]['TSI'].idxmin().date()}), "292 f"max {roll[w]['TSI'].max():.2f} ({roll[w]['TSI'].idxmax().date()})")293roll_orig = pd.read_csv(os.path.join(RESULTS, 'rolling_total_spillover.csv'), parse_dates=['date']).set_index('date')294roll_summary = pd.DataFrame({295 'window': ['original (200 calendar days, levels, forward-filled)', '200 trading days', '250 trading days',296 '500 trading days'],297 'n_windows': [len(roll_orig), len(roll[200]), len(roll[250]), len(roll[500])],298 'mean': [roll_orig['total_spillover'].mean()] + [roll[w]['TSI'].mean() for w in (200, 250, 500)],299 'min': [roll_orig['total_spillover'].min()] + [roll[w]['TSI'].min() for w in (200, 250, 500)],300 'min_date': [roll_orig['total_spillover'].idxmin().date()] + [roll[w]['TSI'].idxmin().date() for w in (200, 250, 500)],301 'max': [roll_orig['total_spillover'].max()] + [roll[w]['TSI'].max() for w in (200, 250, 500)],302 'max_date': [roll_orig['total_spillover'].idxmax().date()] + [roll[w]['TSI'].idxmax().date() for w in (200, 250, 500)],303})304roll_summary.to_csv(os.path.join(RESULTS, 'rev_rolling_summary.csv'), index=False, float_format='%.2f')305306# ─── 7. Regimes (with bootstrap CIs on NET) ───────────────────────────────────307REGIMES = {308 'Pre-COVID (2019-01 to 2020-02)': ('2019-01-01', '2020-02-28'),309 'COVID crash (2020-03 to 2020-06)': ('2020-03-01', '2020-06-30'),310 'Rate hikes (2022-01 to 2022-12)': ('2022-01-01', '2022-12-31'),311 'Spot-ETF era (2024-01 to 2026-03)': ('2024-01-11', '2026-03-11'),312}313reg_rows, reg_full = [], {}314for name, (s, e) in REGIMES.items():315 y = ylog.loc[s:e]316 p = P_BASE if len(y) >= 150 else 1317 b, tb, tob, frb, nb = bootstrap_connectedness(y, p=p, B=B_BOOT)318 tci, nci = ci(tb), ci(nb)319 row = {'Regime': name, 'T': len(y), 'p': p, 'TSI': b['TSI'], 'TSI_lo': tci[0], 'TSI_hi': tci[1]}320 for k, a in enumerate(ASSETS):321 row[f'NET_{a}'] = b['NET'][k]322 row[f'NET_{a}_lo'] = nci[0][k]323 row[f'NET_{a}_hi'] = nci[1][k]324 reg_rows.append(row)325 tbl = pd.DataFrame(b['theta'], index=ASSETS, columns=ASSETS)326 tbl['FROM'] = b['FROM']327 tbl.loc['TO'] = list(b['TO']) + [b['TSI']]328 tbl.loc['NET'] = list(b['NET']) + [np.nan]329 reg_full[name] = tbl330 slug = name.split(' (')[0].lower().replace(' ', '_').replace('-', '')331 tbl.to_csv(os.path.join(RESULTS, f'rev_table_regime_{slug}.csv'), float_format='%.2f')332 log(f"regime {name}: T={len(y)}, p={p}, TSI={b['TSI']:.1f} [{tci[0]:.1f},{tci[1]:.1f}], "333 f"NET BTC={b['NET'][7]:+.1f} [{nci[0][7]:+.1f},{nci[1][7]:+.1f}], NET USDJPY={b['NET'][3]:+.1f}")334 # same regime with p = 4 for comparison when p = 1 was used335 if p != P_BASE:336 alt = connectedness(y, P_BASE, H_BASE)337 log(f" (VAR(4) on the same window: TSI={alt['TSI']:.1f}, NET USDJPY={alt['NET'][3]:+.1f}, "338 f"NET BTC={alt['NET'][7]:+.1f})")339 reg_rows[-1]['TSI_var4'] = alt['TSI']340 reg_rows[-1]['NET_USDJPY_var4'] = alt['NET'][3]341 reg_rows[-1]['NET_BTC_var4'] = alt['NET'][7]342reg_df = pd.DataFrame(reg_rows)343reg_df.to_csv(os.path.join(RESULTS, 'rev_table_regimes.csv'), index=False, float_format='%.2f')344345# ─── 8. Calendar-year sub-samples (crypto integration) ────────────────────────346yr_rows = []347for yr in range(2016, 2027):348 y = ylog.loc[str(yr)]349 if len(y) < 60:350 continue351 b, tb, tob, frb, nb = bootstrap_connectedness(y, p=1 if len(y) < 150 else 2, B=300)352 tci, nci = ci(tb), ci(nb)353 i_btc, i_es, i_eth = ASSETS.index('BTC'), ASSETS.index('ES'), ASSETS.index('ETH')354 th = b['theta']355 yr_rows.append({'Year': yr, 'T': len(y), 'TSI': b['TSI'], 'TSI_lo': tci[0], 'TSI_hi': tci[1],356 'NET_BTC': b['NET'][i_btc], 'NET_BTC_lo': nci[0][i_btc], 'NET_BTC_hi': nci[1][i_btc],357 'NET_ETH': b['NET'][i_eth], 'NET_ETH_lo': nci[0][i_eth], 'NET_ETH_hi': nci[1][i_eth],358 'FROM_BTC': b['FROM'][i_btc], 'FROM_BTC_lo': ci(frb)[0][i_btc], 'FROM_BTC_hi': ci(frb)[1][i_btc],359 'FROM_BTC_excl_ETH': b['FROM'][i_btc] - th[i_btc, i_eth],360 'BTC_from_ES': th[i_btc, i_es], 'ES_from_BTC': th[i_es, i_btc],361 'NET_ES': b['NET'][i_es], 'NET_ES_lo': nci[0][i_es], 'NET_ES_hi': nci[1][i_es]})362yr_df = pd.DataFrame(yr_rows)363yr_df.to_csv(os.path.join(RESULTS, 'rev_table_yearly.csv'), index=False, float_format='%.2f')364log("yearly sub-samples done")365366# ─── 9. Granger causality (log RV, p lags, F and HAC Wald) ────────────────────367PAIRS = [('ES', 'VIX'), ('NQ', 'VIX'), ('VIX', 'ES'), ('VIX', 'NQ'), ('ES', 'NQ'), ('NQ', 'ES'),368 ('ES', 'BTC'), ('NQ', 'BTC'), ('VIX', 'BTC'), ('BTC', 'ES'), ('BTC', 'NQ'), ('BTC', 'VIX'),369 ('USDJPY', 'ES'), ('ES', 'USDJPY'), ('GC', 'BTC'), ('BTC', 'GC')]370371372def granger(y, cause, effect, p=P_BASE, hac_lags=10):373 df = pd.DataFrame({'y': y[effect], 'x': y[cause]})374 X = pd.concat([df['y'].shift(l).rename(f'y_l{l}') for l in range(1, p + 1)] +375 [df['x'].shift(l).rename(f'x_l{l}') for l in range(1, p + 1)], axis=1)376 X = sm.add_constant(X)377 data = pd.concat([df['y'], X], axis=1).dropna()378 yy, XX = data['y'], data.drop(columns='y')379 m_u = sm.OLS(yy, XX).fit()380 m_r = sm.OLS(yy, XX[[c for c in XX.columns if not c.startswith('x_')]]).fit()381 T = len(yy)382 F = ((m_r.ssr - m_u.ssr) / p) / (m_u.ssr / (T - 2 * p - 1))383 from scipy import stats384 pF = 1 - stats.f.cdf(F, p, T - 2 * p - 1)385 m_hac = sm.OLS(yy, XX).fit(cov_type='HAC', cov_kwds={'maxlags': hac_lags})386 R = np.zeros((p, XX.shape[1]))387 for k, c in enumerate([c for c in XX.columns if c.startswith('x_')]):388 R[k, list(XX.columns).index(c)] = 1389 w = m_hac.wald_test(R, use_f=False, scalar=True)390 return {'Cause': cause, 'Effect': effect, 'T': T, 'F': F, 'p_F': pF, 'Wald_HAC': float(w.statistic),391 'p_HAC': float(w.pvalue), 'sum_coef': float(m_u.params[[c for c in XX.columns if c.startswith('x_')]].sum())}392393394gr = pd.DataFrame([granger(ylog, c, e) for c, e in PAIRS])395gr.to_csv(os.path.join(RESULTS, 'rev_table_granger.csv'), index=False, float_format='%.4f')396log("Granger tests done")397398# ─── 10. Lag-order criteria ───────────────────────────────────────────────────399sel = VAR(ylog.values).select_order(maxlags=10)400lag_tbl = pd.DataFrame({'p': range(0, 11), 'AIC': sel.ics['aic'], 'BIC': sel.ics['bic'], 'HQ': sel.ics['hqic']})401lag_tbl.to_csv(os.path.join(RESULTS, 'rev_table_lag_selection.csv'), index=False, float_format='%.3f')402log(f"lag selection: {sel.selected_orders}")403404405# ─── 11. LaTeX tables (generated) ─────────────────────────────────────────────406def fmt(x, d=1, sign=False):407 if pd.isna(x):408 return ''409 s = f"{x:+.{d}f}" if sign else f"{x:.{d}f}"410 return s.replace('-', '$-$').replace('+', '$+$') if sign else s411412413def write(name, text):414 with open(os.path.join(TABLES, name), 'w') as f:415 f.write(text)416417418# 11a. RV descriptives419lines = [r"\begin{tabular}{@{}llrrrrrrrrr@{}}", r"\toprule",420 r"Asset & Class & Mean & Median & Max & Skew & Skew & Kurt & $\rho_1$ & ADF & KPSS \\",421 r" & & \multicolumn{3}{c}{annualized RV (\%)} & level & log & log & log & log & log \\",422 r"\cmidrule(lr){3-5}\cmidrule(lr){6-6}\cmidrule(lr){7-9}\cmidrule(lr){10-11}"]423lines[-1] = r"\midrule"424for a in ASSETS:425 r = desc.loc[a]426 kp = fmt(r['KPSS_log'], 2) + ('' if r['KPSS_p'] > 0.05 else r'$^{*}$')427 adf_st = r'$^{***}$' if r['ADF_p'] < 0.01 else r'$^{**}$' if r['ADF_p'] < 0.05 else r'$^{*}$' if r['ADF_p'] < 0.10 else ''428 lines.append(f"{a} & {r['Class']} & {fmt(r['Mean_RV_pct'])} & {fmt(r['Median_RV_pct'])} & {fmt(r['Max_RV_pct'], 0)} & "429 f"{fmt(r['Skew_level'], 2)} & {fmt(r['Skew_log'], 2)} & {fmt(r['Kurt_log'], 2)} & {fmt(r['AC1_log'], 2)} & "430 f"{fmt(r['ADF_log'], 1)}{adf_st} & {kp} \\\\")431lines += [r"\bottomrule", r"\end{tabular}"]432write('tab_rv_descriptives.tex', "\n".join(lines))433434# 11b. full connectedness table (with CI rows)435def conn_table_tex(tbl, to_ci_=None, from_ci_=None, net_ci_=None, tsi_ci_=None):436 cols = ASSETS437 L = [r"\begin{tabular}{@{}l*{10}{r}r@{}}", r"\toprule",438 r"& \multicolumn{10}{c}{Shock origin $j$} & \\", r"\cmidrule(lr){2-11}",439 "Recipient $i$ & " + " & ".join(cols) + r" & FROM \\", r"\midrule"]440 for i, a in enumerate(cols):441 vals = [fmt(tbl.loc[a, c], 1) for c in cols]442 vals = [f"\\textbf{{{v}}}" if c == a else v for v, c in zip(vals, cols)]443 L.append(f"{a} & " + " & ".join(vals) + f" & {fmt(tbl.loc[a, 'FROM'], 1)} \\\\")444 if from_ci_ is not None:445 L.append(" & " + " & ".join([''] * 10) + f" & {{\\scriptsize[{fmt(from_ci_[0][i], 1)}, {fmt(from_ci_[1][i], 1)}]}} \\\\")446 L.append(r"\midrule")447 L.append("TO & " + " & ".join(fmt(tbl.loc['TO', c], 1) for c in cols) + f" & \\textbf{{{fmt(tbl.loc['TO', 'FROM'], 1)}}} \\\\")448 if to_ci_ is not None:449 L.append(" & " + " & ".join(f"{{\\scriptsize[{fmt(to_ci_[0][k], 1)}, {fmt(to_ci_[1][k], 1)}]}}" for k in range(10)) +450 f" & {{\\scriptsize[{fmt(tsi_ci_[0], 1)}, {fmt(tsi_ci_[1], 1)}]}} \\\\")451 L.append("NET & " + " & ".join(fmt(tbl.loc['NET', c], 1, sign=True) for c in cols) + r" & \\")452 if net_ci_ is not None:453 L.append(" & " + " & ".join(f"{{\\scriptsize[{fmt(net_ci_[0][k], 1, True)}, {fmt(net_ci_[1][k], 1, True)}]}}" for k in range(10)) + r" & \\")454 L += [r"\bottomrule", r"\end{tabular}"]455 return "\n".join(L)456457458write('tab_connectedness_full.tex', conn_table_tex(full_tbl, to_ci, from_ci, net_ci, tsi_ci))459for name, tbl in reg_full.items():460 slug = name.split(' (')[0].lower().replace(' ', '_').replace('-', '')461 write(f'tab_regime_{slug}.tex', conn_table_tex(tbl))462463# 11c. audit / construction table464L = [r"\begin{tabular}{@{}lrrrrrrr@{}}", r"\toprule",465 r"Construction of the panel & $T$ & TSI & ES & NQ & VIX & BTC & USDJPY \\", r"\midrule"]466for _, r in audit.iterrows():467 L.append(f"{r['Construction']} & {int(r['T'])} & {fmt(r['TSI'])} & {fmt(r['NET_ES'], 1, True)} & {fmt(r['NET_NQ'], 1, True)} & "468 f"{fmt(r['NET_VIX'], 1, True)} & {fmt(r['NET_BTC'], 1, True)} & {fmt(r['NET_USDJPY'], 1, True)} \\\\")469L += [r"\bottomrule", r"\end{tabular}"]470write('tab_audit_construction.tex', "\n".join(L))471472# 11d. robustness grid473L = [r"\begin{tabular}{@{}lrrrrrrrrr@{}}", r"\toprule",474 r"Specification & $T$ & TSI & ES & NQ & USDJPY & CL & BTC & VIX & BTC own \\", r"\midrule"]475for _, r in rob_df.iterrows():476 L.append(f"{r['Specification']} & {int(r['T'])} & {fmt(r['TSI'])} & {fmt(r['NET_ES'], 1, True)} & {fmt(r['NET_NQ'], 1, True)} & "477 f"{fmt(r['NET_USDJPY'], 1, True)} & {fmt(r['NET_CL'], 1, True)} & {fmt(r['NET_BTC'], 1, True)} & "478 f"{fmt(r['NET_VIX'], 1, True)} & {fmt(r.get('Own_BTC', np.nan))} \\\\")479L += [r"\bottomrule", r"\end{tabular}"]480write('tab_robustness.tex', "\n".join(L))481482# 11e. regimes summary with CIs483L = [r"\begin{tabular}{@{}l" + "r" * len(REGIMES) + "@{}}", r"\toprule",484 " & " + " & ".join(n.split(' (')[0] for n in REGIMES) + r" \\",485 " & " + " & ".join("{\\scriptsize " + n.split(' (')[1].rstrip(')') + "}" for n in REGIMES) + r" \\", r"\midrule",486 "$T$ (trading days) & " + " & ".join(str(int(r['T'])) for _, r in reg_df.iterrows()) + r" \\",487 "VAR lag $p$ & " + " & ".join(str(int(r['p'])) for _, r in reg_df.iterrows()) + r" \\",488 "TSI (\\%) & " + " & ".join(fmt(r['TSI']) for _, r in reg_df.iterrows()) + r" \\",489 " & " + " & ".join(f"{{\\scriptsize[{fmt(r['TSI_lo'])}, {fmt(r['TSI_hi'])}]}}" for _, r in reg_df.iterrows()) + r" \\",490 r"\midrule", r"\multicolumn{" + str(1 + len(REGIMES)) + r"}{@{}l}{\textit{Net directional spillovers (pp)}} \\"]491for a in ASSETS:492 L.append(f"{a} & " + " & ".join(fmt(r[f'NET_{a}'], 1, True) for _, r in reg_df.iterrows()) + r" \\")493 L.append(" & " + " & ".join(f"{{\\scriptsize[{fmt(r[f'NET_{a}_lo'], 1, True)}, {fmt(r[f'NET_{a}_hi'], 1, True)}]}}"494 for _, r in reg_df.iterrows()) + r" \\")495L += [r"\bottomrule", r"\end{tabular}"]496write('tab_regimes.tex', "\n".join(L))497498# 11f. yearly table499L = [r"\begin{tabular}{@{}lrrrrrrr@{}}", r"\toprule",500 r"Year & $T$ & TSI & NET BTC & NET ETH & FROM BTC & BTC $\leftarrow$ ES & ES $\leftarrow$ BTC \\",501 r" & & (\%) & (pp) & (pp) & excl.\ ETH (\%) & (\%) & (\%) \\", r"\midrule"]502for _, r in yr_df.iterrows():503 L.append(f"{int(r['Year'])} & {int(r['T'])} & {fmt(r['TSI'])} {{\\scriptsize[{fmt(r['TSI_lo'])}, {fmt(r['TSI_hi'])}]}} & "504 f"{fmt(r['NET_BTC'], 1, True)} {{\\scriptsize[{fmt(r['NET_BTC_lo'], 1, True)}, {fmt(r['NET_BTC_hi'], 1, True)}]}} & "505 f"{fmt(r['NET_ETH'], 1, True)} & {fmt(r['FROM_BTC_excl_ETH'])} & {fmt(r['BTC_from_ES'])} & {fmt(r['ES_from_BTC'])} \\\\")506L += [r"\bottomrule", r"\end{tabular}"]507write('tab_yearly.tex', "\n".join(L))508509# 11g. Granger510def stars(p):511 return r'$^{***}$' if p < 0.01 else r'$^{**}$' if p < 0.05 else r'$^{*}$' if p < 0.10 else ''512513514L = [r"\begin{tabular}{@{}llrrrrr@{}}", r"\toprule",515 r"Cause & Effect & $T$ & $F$ & $p$ & Wald (HAC) & $p$ \\", r"\midrule"]516for _, r in gr.iterrows():517 L.append(f"{r['Cause']} & {r['Effect']} & {int(r['T'])} & {fmt(r['F'], 2)}{stars(r['p_F'])} & {fmt(r['p_F'], 3)} & "518 f"{fmt(r['Wald_HAC'], 2)}{stars(r['p_HAC'])} & {fmt(r['p_HAC'], 3)} \\\\")519L += [r"\bottomrule", r"\end{tabular}"]520write('tab_granger.tex', "\n".join(L))521522# 11h. lag selection523L = [r"\begin{tabular}{@{}lrrrrrrrrrr@{}}", r"\toprule", "Criterion & " + " & ".join(f"$p={p}$" for p in range(1, 11)) + r" \\", r"\midrule"]524for crit in ('AIC', 'BIC', 'HQ'):525 vals = lag_tbl[crit].values[1:]526 best = np.argmin(vals)527 L.append(f"{crit} & " + " & ".join((f"\\textbf{{{v:.2f}}}" if k == best else f"{v:.2f}") for k, v in enumerate(vals)) + r" \\")528L += [r"\bottomrule", r"\end{tabular}"]529write('tab_lag_selection.tex', "\n".join(L))530531# 11i. rolling summary532L = [r"\begin{tabular}{@{}lrrrrrr@{}}", r"\toprule", r"Window & Windows & Mean & Min & (date) & Max & (date) \\", r"\midrule"]533for _, r in roll_summary.iterrows():534 L.append(f"{r['window']} & {int(r['n_windows'])} & {fmt(r['mean'])} & {fmt(r['min'])} & {r['min_date']} & {fmt(r['max'])} & {r['max_date']} \\\\")535L += [r"\bottomrule", r"\end{tabular}"]536write('tab_rolling_summary.tex', "\n".join(L))537log("LaTeX tables written")538539# ─── 12. Figures ──────────────────────────────────────────────────────────────540# Fig 2: heatmap541fig, ax = plt.subplots(figsize=(10, 7.5))542mat = base['theta']543im = ax.imshow(mat, cmap='Blues', vmin=0, vmax=mat.max())544for i in range(10):545 for j in range(10):546 ax.text(j, i, f"{mat[i, j]:.1f}", ha='center', va='center', fontsize=8,547 color='white' if mat[i, j] > 0.55 * mat.max() else 'black')548ax.set_xticks(range(10)); ax.set_yticks(range(10))549ax.set_xticklabels(ASSETS, rotation=45, ha='right'); ax.set_yticklabels(ASSETS)550ax.set_xlabel('Shock origin $j$'); ax.set_ylabel('Recipient $i$')551ax.set_title(f'Generalized FEVD shares (%), log RV, VAR(4), H = 10 — TSI = {base["TSI"]:.1f}%')552plt.colorbar(im, ax=ax, label='% of 10-day forecast-error variance')553plt.tight_layout(); savefig(fig, 'fig2_connectedness_heatmap.png')554555# Fig 3: rolling TSI (three windows) + original for reference556fig, ax = plt.subplots(figsize=(13, 5))557ax.plot(roll[200].index, roll[200]['TSI'], color='#1f77b4', lw=1.0, label='200 trading days')558ax.plot(roll[250].index, roll[250]['TSI'], color='#ff7f0e', lw=0.9, alpha=0.9, label='250 trading days')559ax.plot(roll[500].index, roll[500]['TSI'], color='#2ca02c', lw=1.3, label='500 trading days')560ax.axhline(base['TSI'], color='gray', ls='--', lw=0.8, label=f'full sample ({base["TSI"]:.1f}%)')561events = {'COVID crash\n(Mar 2020)': '2020-03-16', 'First Fed hike\n(Mar 2022)': '2022-03-16',562 'FTX\n(Nov 2022)': '2022-11-11', 'Spot BTC ETFs\n(Jan 2024)': '2024-01-11', 'Volatility spike\n(Aug 2024)': '2024-08-05'}563for k, (lab, d) in enumerate(events.items()):564 d = pd.Timestamp(d)565 ax.axvline(d, color='k', ls=':', lw=0.8, alpha=0.6)566 ax.text(d, ax.get_ylim()[1] - (k % 2) * 3, lab, ha='center', va='top', fontsize=7)567ax.set_ylabel('Total spillover index (%)'); ax.set_xlabel('End of window')568ax.set_title('Rolling total spillover index — log RV, VAR(4), H = 10, trading-day panel')569ax.legend(loc='lower left', ncol=4, fontsize=8)570plt.tight_layout(); savefig(fig, 'fig3_rolling_spillover.png')571572# Fig 4: NET by regime with CI bars573fig, ax = plt.subplots(figsize=(12, 5.5))574x = np.arange(10); nreg = len(REGIMES); width = 0.8 / nreg575pal = ['#1f77b4', '#d62728', '#ff7f0e', '#9467bd']576for k, (_, r) in enumerate(reg_df.iterrows()):577 vals = np.array([r[f'NET_{a}'] for a in ASSETS])578 lo = np.array([r[f'NET_{a}_lo'] for a in ASSETS]); hi = np.array([r[f'NET_{a}_hi'] for a in ASSETS])579 off = (k - (nreg - 1) / 2) * width580 ax.bar(x + off, vals, width * 0.9, color=pal[k], alpha=0.85, label=f"{r['Regime']} (TSI {r['TSI']:.1f}%)")581 ax.errorbar(x + off, vals, yerr=[vals - lo, hi - vals], fmt='none', ecolor='k', elinewidth=0.6, capsize=1.5)582ax.axhline(0, color='k', lw=0.8); ax.set_xticks(x); ax.set_xticklabels(ASSETS)583ax.set_ylabel('Net directional spillover (pp)'); ax.legend(fontsize=8)584ax.set_title('Net directional spillovers by regime, with 90% bootstrap intervals')585plt.tight_layout(); savefig(fig, 'fig4_net_spillovers_regimes.png')586587# Fig 6: rolling NET BTC vs ES (w = 200 and 500)588fig, ax = plt.subplots(figsize=(13, 5))589ax.plot(roll[200].index, roll[200]['NET_BTC'], color='#9467bd', lw=0.8, alpha=0.6, label='BTC, 200 days')590ax.plot(roll[500].index, roll[500]['NET_BTC'], color='#9467bd', lw=1.6, label='BTC, 500 days')591ax.plot(roll[200].index, roll[200]['NET_ES'], color='#1f77b4', lw=0.8, alpha=0.6, label='ES, 200 days')592ax.plot(roll[500].index, roll[500]['NET_ES'], color='#1f77b4', lw=1.6, label='ES, 500 days')593ax.axhline(0, color='k', lw=0.8, ls='--')594for lab, d in events.items():595 ax.axvline(pd.Timestamp(d), color='k', ls=':', lw=0.8, alpha=0.6)596ax.set_ylabel('Net directional spillover (pp)'); ax.set_xlabel('End of window')597ax.set_title('Rolling net directional spillovers: Bitcoin and E-mini S&P 500')598ax.legend(ncol=4, fontsize=8, loc='lower left')599plt.tight_layout(); savefig(fig, 'fig6_net_directional_btc_es.png')600601# Fig 7: network of net pairwise directional spillovers602netpair = base['NETPAIR']603thr = np.percentile(np.abs(netpair[np.triu_indices(10, 1)]), 60)604G = nx.DiGraph()605for a in ASSETS:606 G.add_node(a)607for i, a in enumerate(ASSETS):608 for j, b in enumerate(ASSETS):609 if i < j:610 v = netpair[i, j]611 if abs(v) > thr:612 if v > 0:613 G.add_edge(a, b, weight=v)614 else:615 G.add_edge(b, a, weight=-v)616order = ['ES', 'NQ', 'VIX', 'GC', 'CL', 'EURUSD', 'USDJPY', 'GBPUSD', 'BTC', 'ETH']617pos = {a: (np.cos(2 * np.pi * k / 10 + np.pi / 2), np.sin(2 * np.pi * k / 10 + np.pi / 2)) for k, a in enumerate(order)}618fig, ax = plt.subplots(figsize=(10.5, 8))619sizes = [300 + 25 * base['TO'][ASSETS.index(a)] for a in G.nodes()]620colors = [CLASS_COLOR[ASSET_CLASS[a]] for a in G.nodes()]621wts = np.array([G[u][v]['weight'] for u, v in G.edges()])622nx.draw_networkx_edges(G, pos, ax=ax, width=0.6 + 3.5 * wts / wts.max(), alpha=0.7, edge_color='#555555',623 arrows=True, arrowsize=14, connectionstyle='arc3,rad=0.12', node_size=sizes)624nx.draw_networkx_nodes(G, pos, ax=ax, node_size=sizes, node_color=colors, alpha=0.92, edgecolors='white', linewidths=1.5)625nx.draw_networkx_labels(G, pos, ax=ax, font_size=9, font_weight='bold')626ax.legend(handles=[mpatches.Patch(color=c, label=k) for k, c in CLASS_COLOR.items()], loc='upper left', fontsize=8)627ax.set_title(f'Net pairwise directional spillovers above the 60th percentile (|net| > {thr:.1f} pp)\n'628 'arrow from net transmitter to net receiver; node size proportional to TO spillover')629ax.axis('off'); plt.tight_layout(); savefig(fig, 'fig7_network_graph.png')630631# Fig 8: yearly BTC net and FROM with CIs632fig, axes = plt.subplots(1, 2, figsize=(12, 4.2))633ax = axes[0]634ax.errorbar(yr_df['Year'], yr_df['NET_BTC'], yerr=[yr_df['NET_BTC'] - yr_df['NET_BTC_lo'], yr_df['NET_BTC_hi'] - yr_df['NET_BTC']],635 fmt='o-', color='#9467bd', capsize=3, label='BTC')636ax.errorbar(yr_df['Year'], yr_df['NET_ES'], yerr=[yr_df['NET_ES'] - yr_df['NET_ES_lo'], yr_df['NET_ES_hi'] - yr_df['NET_ES']],637 fmt='s-', color='#1f77b4', capsize=3, label='ES')638ax.axhline(0, color='k', lw=0.8, ls='--'); ax.set_ylabel('Net directional spillover (pp)'); ax.set_title('(a) Net spillovers by calendar year')639ax.legend()640ax = axes[1]641ax.plot(yr_df['Year'], yr_df['FROM_BTC_excl_ETH'], 'o-', color='#9467bd', label='BTC: FROM traditional assets (excl. ETH)')642ax.plot(yr_df['Year'], yr_df['BTC_from_ES'], 's--', color='#1f77b4', label='BTC $\\leftarrow$ ES')643ax.plot(yr_df['Year'], yr_df['ES_from_BTC'], '^--', color='#ff7f0e', label='ES $\\leftarrow$ BTC')644ax.set_ylabel('% of forecast-error variance'); ax.set_title('(b) Bitcoin–traditional-asset shares by year'); ax.legend(fontsize=8)645plt.tight_layout(); savefig(fig, 'fig8_yearly_crypto_integration.png')646log("figures written")647648# ─── 13. Key numbers for the paper ────────────────────────────────────────────649key = {650 'T_panel': len(panel), 'T_original_panel': len(orig_panel), 'weekend_rows_original': n_weekend,651 'ffilled_rows_original': n_ffilled_rows, 'TSI': base['TSI'], 'TSI_lo': tsi_ci[0], 'TSI_hi': tsi_ci[1],652 'TSI_original_replicated': orig['TSI'], 'TSI_levels_trading_days': lvl_td['TSI'],653}654for k, a in enumerate(ASSETS):655 key[f'NET_{a}'] = base['NET'][k]; key[f'NET_{a}_lo'] = net_ci[0][k]; key[f'NET_{a}_hi'] = net_ci[1][k]656 key[f'TO_{a}'] = base['TO'][k]; key[f'FROM_{a}'] = base['FROM'][k]; key[f'OWN_{a}'] = base['theta'][k, k]657key['BTC_ETH'] = base['theta'][7, 8]; key['ETH_BTC'] = base['theta'][8, 7]658key['VIX_from_ES'] = base['theta'][9, 0]; key['VIX_from_NQ'] = base['theta'][9, 1]659key['BTC_from_traditional'] = base['FROM'][7] - base['theta'][7, 8]660key['ETH_from_traditional'] = base['FROM'][8] - base['theta'][8, 7]661for w in (200, 250, 500):662 key[f'roll{w}_mean'] = roll[w]['TSI'].mean(); key[f'roll{w}_min'] = roll[w]['TSI'].min(); key[f'roll{w}_max'] = roll[w]['TSI'].max()663 key[f'roll{w}_min_date'] = str(roll[w]['TSI'].idxmin().date()); key[f'roll{w}_max_date'] = str(roll[w]['TSI'].idxmax().date())664 key[f'roll{w}_last_NET_BTC'] = roll[w]['NET_BTC'].iloc[-1]; key[f'roll{w}_last_NET_ES'] = roll[w]['NET_ES'].iloc[-1]665 key[f'roll{w}_NET_BTC_2025_mean'] = roll[w].loc['2025-01-01':'2025-12-31', 'NET_BTC'].mean()666 r17 = roll[w].loc['2017-01-01':'2017-12-31', 'NET_BTC']667 key[f'roll{w}_NET_BTC_2017_mean'] = r17.mean() if len(r17) else np.nan668 key[f'roll{w}_share_BTC_positive'] = (roll[w]['NET_BTC'] > 0).mean()669# rolling TSI around Feb 2018 and Mar 2020670key['roll200_2018-02-02'] = roll[200]['TSI'].get(pd.Timestamp('2018-02-02'), np.nan)671key['roll200_2018-02-06'] = roll[200]['TSI'].get(pd.Timestamp('2018-02-06'), np.nan)672key['roll200_2020-02-20'] = roll[200]['TSI'].get(pd.Timestamp('2020-02-20'), np.nan)673key['roll200_2020-03-31'] = roll[200]['TSI'].get(pd.Timestamp('2020-03-31'), np.nan)674key['roll200_2020-06-30'] = roll[200]['TSI'].get(pd.Timestamp('2020-06-30'), np.nan)675pd.Series(key).to_csv(os.path.join(RESULTS, 'rev_key_numbers.csv'), header=['value'])676log("key numbers saved; done")677print(pd.Series(key).to_string())678