TeX 62.4%
Python 37.5%
1"""2Intraday Volatility Transmission Across Asset Classes3Diebold-Yilmaz Connectedness Analysis -- stage 1 (May 2026 pipeline)45NOTE (revision of 2026-10-17, see docs/CHANGELOG-REVISION.md):6 * This script is kept as the sole producer of results/daily_rv_all.csv and7 results/table1_descriptive_stats.csv from the vendor's 5-minute DuckDB bars.8 * Its sections 3-8 (VAR/GFEVD, rolling, regimes, sessions, Granger, network)9 are SUPERSEDED by scripts/02_connectedness_revision.py and10 scripts/03_intraday_sessions.py. Known defects of the original stages:11 - the VAR panel forward-filled weekends/holidays (rv_panel.ffill(limit=5)),12 so 1,137 of 3,651 rows contained stale values;13 - realized volatility entered the VAR in levels, not logs;14 - FirstRate time stamps are U.S. Eastern, not UTC (comment below is wrong),15 so the session labels of the original Fig. 5 were shifted by 4-5 hours;16 - the "session" Granger tests repeated the same daily test three times;17 - 5-minute returns beyond 20 sample standard deviations are dropped before18 computing RV (undisclosed in v1.0; disclosed in v1.1).19 The code below is left unchanged for the audit trail.20"""2122import os23import warnings24import numpy as np25import pandas as pd26import matplotlib27matplotlib.use('Agg')28import matplotlib.pyplot as plt29import matplotlib.patches as mpatches30import matplotlib.gridspec as gridspec31import seaborn as sns32import duckdb33from scipy import stats34import networkx as nx35from statsmodels.tsa.api import VAR36from statsmodels.tsa.stattools import grangercausalitytests37import traceback3839warnings.filterwarnings('ignore')4041# ─── Paths ────────────────────────────────────────────────────────────────────42ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))43DATA_DIR = os.environ.get('WP6_DATA_DIR', os.path.join(ROOT, 'data', 'raw')) # *_5min.duckdb files (see data/raw/README.md)44RESULTS = os.path.join(ROOT, 'results')45FIGURES = os.path.join(ROOT, 'figures')46os.makedirs(FIGURES, exist_ok=True)47os.makedirs(RESULTS, exist_ok=True)4849plt.rcParams.update({50 'font.size': 10, 'axes.titlesize': 11, 'axes.labelsize': 10,51 'xtick.labelsize': 9, 'ytick.labelsize': 9, 'legend.fontsize': 9,52 'figure.dpi': 150, 'savefig.dpi': 150, 'font.family': 'DejaVu Sans',53 'axes.spines.top': False, 'axes.spines.right': False,54})5556ASSET_COLORS = {57 'ES': '#1f77b4', 'NQ': '#aec7e8',58 'EURUSD': '#2ca02c', 'USDJPY': '#98df8a', 'GBPUSD': '#d62728',59 'GC': '#ff7f0e', 'CL': '#ffbb78',60 'BTC': '#9467bd', 'ETH': '#c5b0d5',61 'VIX': '#8c564b',62}63ASSET_CLASS = {64 'ES': 'Equity', 'NQ': 'Equity',65 'EURUSD': 'FX', 'USDJPY': 'FX', 'GBPUSD': 'FX',66 'GC': 'Commodity', 'CL': 'Commodity',67 'BTC': 'Crypto', 'ETH': 'Crypto',68 'VIX': 'Volatility',69}70CLASS_COLOR = {71 'Equity': '#1f77b4', 'FX': '#2ca02c',72 'Commodity': '#ff7f0e', 'Crypto': '#9467bd', 'Volatility': '#8c564b',73}7475ASSETS = ['ES', 'NQ', 'EURUSD', 'USDJPY', 'GBPUSD', 'GC', 'CL', 'BTC', 'ETH', 'VIX']7677# ─── 1. Load 5-min data ───────────────────────────────────────────────────────78print("=" * 70)79print("Loading 5-min OHLCV data...")8081def load_asset(db_name, symbol):82 db_path = os.path.join(DATA_DIR, db_name)83 conn = duckdb.connect(db_path, read_only=True)84 df = conn.execute(85 f"SELECT datetime, close FROM ohlcv WHERE symbol='{symbol}' ORDER BY datetime"86 ).fetchdf()87 conn.close()88 df['datetime'] = pd.to_datetime(df['datetime'])89 df = df.set_index('datetime')90 df.columns = [symbol]91 return df9293futures_assets = {94 'ES': 'futures_5min.duckdb',95 'NQ': 'futures_5min.duckdb',96 'GC': 'futures_5min.duckdb',97 'CL': 'futures_5min.duckdb',98 'VX': 'futures_5min.duckdb',99}100fx_assets = {101 'EURUSD': 'fx_5min.duckdb',102 'USDJPY': 'fx_5min.duckdb',103 'GBPUSD': 'fx_5min.duckdb',104}105crypto_assets = {106 'BTC': 'crypto_5min.duckdb',107 'ETH': 'crypto_5min.duckdb',108}109index_assets = {110 'VIX': 'index_5min.duckdb',111}112113raw_data = {}114for sym, db in {**futures_assets, **fx_assets, **crypto_assets, **index_assets}.items():115 print(f" Loading {sym} from {db}...")116 raw_data[sym] = load_asset(db, sym)117118print("Data loaded.\n")119120# ─── Compute 5-min log returns ────────────────────────────────────────────────121print("Computing 5-min log returns...")122log_returns = {}123for sym, df in raw_data.items():124 lr = np.log(df[sym] / df[sym].shift(1)).dropna()125 # Remove extreme outliers (> 20 std)126 mu, sigma = lr.mean(), lr.std()127 lr = lr[(lr - mu).abs() < 20 * sigma]128 log_returns[sym] = lr129130# ─── 1. Descriptive Statistics ────────────────────────────────────────────────131print("Computing descriptive statistics...")132desc_rows = []133for sym in ASSET_CLASS.keys():134 if sym not in log_returns:135 # VX → treat VX as proxy for VIX136 continue137 lr = log_returns[sym]138 row = {139 'Asset': sym,140 'Class': ASSET_CLASS.get(sym, 'Other'),141 'N_obs': len(lr),142 'Date_Start': lr.index.min().strftime('%Y-%m-%d'),143 'Date_End': lr.index.max().strftime('%Y-%m-%d'),144 'Mean': lr.mean(),145 'Std': lr.std(),146 'Skewness': stats.skew(lr),147 'Kurtosis': stats.kurtosis(lr),148 'Min': lr.min(),149 'Max': lr.max(),150 }151 desc_rows.append(row)152153desc_df = pd.DataFrame(desc_rows)154desc_df.to_csv(os.path.join(RESULTS, 'table1_descriptive_stats.csv'), index=False, float_format='%.6f')155print(" Saved table1_descriptive_stats.csv")156print(desc_df[['Asset', 'N_obs', 'Mean', 'Std', 'Skewness', 'Kurtosis']].to_string(index=False))157158# ─── 2. Realized Volatility ───────────────────────────────────────────────────159print("\nComputing daily realized volatility...")160161def compute_rv(log_ret_series, min_obs=20):162 """Annualized daily realized volatility."""163 df = log_ret_series.to_frame(name='ret')164 df['date'] = df.index.date165 grouped = df.groupby('date')['ret']166 rv = grouped.apply(lambda x: np.sqrt(np.sum(x**2)) * np.sqrt(252) if len(x) >= min_obs else np.nan)167 rv.index = pd.to_datetime(rv.index)168 return rv.dropna()169170rv_dict = {}171for sym in ASSET_CLASS.keys():172 if sym not in log_returns:173 continue174 rv_dict[sym] = compute_rv(log_returns[sym])175 print(f" {sym}: {len(rv_dict[sym])} trading days, RV mean = {rv_dict[sym].mean():.4f}")176177# Combine into DataFrame178rv_df = pd.DataFrame(rv_dict)179rv_df.index.name = 'date'180rv_df.to_csv(os.path.join(RESULTS, 'daily_rv_all.csv'))181print(" Saved daily_rv_all.csv")182183# ─── Fig 1: RV Time Series ────────────────────────────────────────────────────184print("\nPlotting Fig 1: RV Time Series...")185fig, axes = plt.subplots(5, 2, figsize=(14, 16), sharex=False)186axes = axes.flatten()187188asset_order = ['ES', 'NQ', 'EURUSD', 'USDJPY', 'GBPUSD', 'GC', 'CL', 'BTC', 'ETH', 'VIX']189for i, sym in enumerate(asset_order):190 if sym not in rv_dict:191 continue192 ax = axes[i]193 rv_s = rv_dict[sym]194 color = ASSET_COLORS.get(sym, 'steelblue')195 ax.plot(rv_s.index, rv_s.values, color=color, linewidth=0.7, alpha=0.9)196 ax.fill_between(rv_s.index, rv_s.values, alpha=0.15, color=color)197 ax.set_title(f'{sym} ({ASSET_CLASS.get(sym, "")})')198 ax.set_ylabel('Annualized RV')199 ax.xaxis.set_major_locator(matplotlib.dates.YearLocator(2))200 ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))201 plt.setp(ax.xaxis.get_majorticklabels(), rotation=45, ha='right')202203plt.suptitle('Daily Realized Volatility — All Assets (5-min, Annualized)', fontsize=13, y=1.01)204plt.tight_layout()205fig.savefig(os.path.join(FIGURES, 'fig1_rv_timeseries.png'), dpi=150, bbox_inches='tight')206plt.close()207print(" Saved fig1_rv_timeseries.png")208209# ─── Helper: Diebold-Yilmaz GFEVD ────────────────────────────────────────────210def compute_gfevd(var_result, H=10):211 """212 Generalized Forecast Error Variance Decomposition (Pesaran-Shin).213 Returns an (n x n) matrix where entry [i,j] = share of forecast error214 variance of variable i attributable to shocks from variable j.215 """216 n = var_result.k_ar # actually number of variables217 coefs = var_result.coefs # shape (p, n, n)218 sigma = var_result.sigma_u # (n, n) residual covariance219 n = sigma.shape[0]220 p = coefs.shape[0]221222 # Build MA coefficient matrices Psi_h via recursive computation223 # Psi_0 = I, Psi_h = sum_{j=1..min(h,p)} A_j @ Psi_{h-j}224 psi = [np.eye(n)]225 for h in range(1, H + 1):226 psi_h = np.zeros((n, n))227 for j in range(1, min(h, p) + 1):228 psi_h += coefs[j-1] @ psi[h-j]229 psi.append(psi_h)230231 # GFEVD: theta[i,j] = (1/sigma[j,j]) * sum_{h=0}^{H-1} (e_i' Psi_h Sigma e_j)^2232 # / sum_{h=0}^{H-1} (e_i' Psi_h Sigma Psi_h' e_i)233 sigma_diag = np.diag(sigma)234 gfevd = np.zeros((n, n))235236 for i in range(n):237 ei = np.zeros(n); ei[i] = 1.0238 denom = 0.0239 for h in range(H):240 denom += ei @ psi[h] @ sigma @ psi[h].T @ ei241 for j in range(n):242 ej = np.zeros(n); ej[j] = 1.0243 num = 0.0244 for h in range(H):245 num += (ei @ psi[h] @ sigma @ ej) ** 2246 gfevd[i, j] = num / (sigma_diag[j] * denom)247248 # Row-normalize so rows sum to 1249 row_sums = gfevd.sum(axis=1, keepdims=True)250 row_sums[row_sums == 0] = 1.0251 gfevd = gfevd / row_sums252 return gfevd253254255def compute_connectedness(rv_panel, assets, lag=4, H=10):256 """257 Given a DataFrame of RV series (rows=dates, cols=assets),258 fit VAR(lag), compute GFEVD at horizon H, return connectedness table.259 """260 # Use only complete rows261 data = rv_panel[assets].dropna()262 if len(data) < lag + H + 10:263 return None, None264265 model = VAR(data)266 try:267 res = model.fit(lag, trend='c')268 except Exception as e:269 print(f" VAR fit failed: {e}")270 return None, None271272 gfevd = compute_gfevd(res, H=H)273 n = len(assets)274275 # Build full connectedness table276 # Rows = "TO" (recipient of shocks), Cols = "FROM" (origin of shocks)277 # But DY convention: entry [i,j] = contribution of j to i's forecast error278 # TO_j = sum_{i != j} theta[i,j] / N * 100 (what j contributes to others)279 # FROM_i = sum_{j != i} theta[i,j] / N * 100 (what i receives from others)280281 # gfevd[i,j] = share of i's FEVD from j's shock282 # So: FROM_i = sum_{j!=i} gfevd[i,j] * 100283 # TO_j = sum_{i!=j} gfevd[i,j] * 100284285 # Total spillover index286 off_diag = gfevd.copy()287 np.fill_diagonal(off_diag, 0.0)288 TSI = off_diag.sum() / n * 100289290 # Directional291 FROM = off_diag.sum(axis=1) * 100 # what each asset receives292 TO = off_diag.sum(axis=0) * 100 # what each asset transmits293 NET = TO - FROM294295 # Build output DataFrame (DY table format)296 # gfevd * 100 → percentages297 tbl = pd.DataFrame(gfevd * 100, index=assets, columns=assets)298 tbl['FROM'] = FROM299 to_row = pd.Series(TO, index=assets)300 net_row = pd.Series(NET, index=assets)301 tbl.loc['TO'] = to_row302 tbl.loc['NET'] = net_row303 tbl.loc['TO', 'FROM'] = TSI304 tbl.loc['NET', 'FROM'] = np.nan305306 return tbl, TSI307308309# ─── 3. Full-Sample Connectedness ─────────────────────────────────────────────310print("\n" + "="*70)311print("Section 3: Full-Sample Connectedness (VAR-4, H=10)...")312313# Use VIX from index (not VX futures) as VIX is more common314# Find common sample: align all RV series315common_assets = [a for a in ASSET_CLASS.keys() if a in rv_dict]316rv_panel = pd.DataFrame({a: rv_dict[a] for a in common_assets})317318# Forward fill up to 5 days to handle weekends/holidays319rv_panel = rv_panel.ffill(limit=5)320321# Drop rows where more than 3 assets are missing322rv_panel = rv_panel.dropna(thresh=len(common_assets) - 2)323rv_panel = rv_panel.dropna()324325print(f" Common sample: {rv_panel.index[0].date()} to {rv_panel.index[-1].date()}, {len(rv_panel)} days")326print(f" Assets: {common_assets}")327328conn_table, tsi_full = compute_connectedness(rv_panel, common_assets, lag=4, H=10)329if conn_table is not None:330 conn_table.to_csv(os.path.join(RESULTS, 'table2_connectedness_full.csv'), float_format='%.2f')331 print(f" Total Spillover Index (full sample): {tsi_full:.2f}%")332 print(" NET spillovers:")333 for a in common_assets:334 print(f" {a:10s}: NET = {conn_table.loc['NET', a]:+.2f}%")335 print(" Saved table2_connectedness_full.csv")336else:337 print(" WARNING: Full-sample connectedness failed.")338339# ─── Fig 2: Connectedness Heatmap ─────────────────────────────────────────────340if conn_table is not None:341 print("\nPlotting Fig 2: Connectedness Heatmap...")342 heat_data = conn_table.loc[common_assets, common_assets].values.astype(float)343 fig, ax = plt.subplots(figsize=(11, 8))344 mask_diag = np.eye(len(common_assets), dtype=bool)345 im = ax.imshow(heat_data, cmap='RdYlBu_r', aspect='auto', vmin=0, vmax=heat_data.max())346 for i in range(len(common_assets)):347 for j in range(len(common_assets)):348 txt = f"{heat_data[i,j]:.1f}"349 color = 'white' if heat_data[i,j] > heat_data.max() * 0.6 else 'black'350 ax.text(j, i, txt, ha='center', va='center', fontsize=8, color=color, fontweight='bold')351 ax.set_xticks(range(len(common_assets)))352 ax.set_yticks(range(len(common_assets)))353 ax.set_xticklabels(common_assets, rotation=45, ha='right')354 ax.set_yticklabels(common_assets)355 ax.set_xlabel('FROM (shock origin)')356 ax.set_ylabel('TO (shock recipient)')357 ax.set_title(f'Diebold-Yilmaz Connectedness Matrix (Full Sample)\nTotal Spillover Index = {tsi_full:.1f}%',358 fontweight='bold')359 plt.colorbar(im, ax=ax, label='% of Forecast Error Variance')360 plt.tight_layout()361 fig.savefig(os.path.join(FIGURES, 'fig2_connectedness_heatmap.png'), dpi=150, bbox_inches='tight')362 plt.close()363 print(" Saved fig2_connectedness_heatmap.png")364365# ─── 4. Rolling Connectedness ─────────────────────────────────────────────────366print("\n" + "="*70)367print("Section 4: Rolling Connectedness (200-day window, step=1)...")368369WINDOW = 200370dates_all = rv_panel.index.tolist()371rolling_dates = []372rolling_tsi = []373374n_steps = len(dates_all) - WINDOW375print(f" Total windows to compute: {n_steps}")376377for i in range(0, n_steps, 1):378 window_data = rv_panel.iloc[i : i + WINDOW]379 end_date = dates_all[i + WINDOW - 1]380381 # Only compute if at least 150 non-NaN rows after dropna382 sub = window_data.dropna()383 if len(sub) < 100:384 continue385386 model = VAR(sub)387 try:388 res = model.fit(4, trend='c')389 gfevd = compute_gfevd(res, H=10)390 n = gfevd.shape[0]391 off_diag = gfevd.copy()392 np.fill_diagonal(off_diag, 0.0)393 tsi = off_diag.sum() / n * 100394 rolling_dates.append(end_date)395 rolling_tsi.append(tsi)396 except Exception:397 pass398399 if i % 200 == 0:400 print(f" Step {i}/{n_steps}: {end_date.date()}, TSI = {rolling_tsi[-1]:.1f}%" if rolling_tsi else f" Step {i}/{n_steps}")401402rolling_df = pd.DataFrame({'date': rolling_dates, 'total_spillover': rolling_tsi})403rolling_df = rolling_df.set_index('date')404rolling_df.to_csv(os.path.join(RESULTS, 'rolling_total_spillover.csv'))405print(f" Rolling analysis complete: {len(rolling_df)} observations")406print(f" Mean TSI: {rolling_df['total_spillover'].mean():.2f}%")407print(f" Min TSI: {rolling_df['total_spillover'].min():.2f}%")408print(f" Max TSI: {rolling_df['total_spillover'].max():.2f}%")409410# ─── Fig 3: Rolling Spillover ─────────────────────────────────────────────────411print("\nPlotting Fig 3: Rolling Total Spillover...")412fig, ax = plt.subplots(figsize=(13, 5))413ax.plot(rolling_df.index, rolling_df['total_spillover'], color='#1f77b4', linewidth=1.2, label='Total Spillover Index')414ax.fill_between(rolling_df.index, rolling_df['total_spillover'],415 rolling_df['total_spillover'].mean(), alpha=0.15, color='#1f77b4')416ax.axhline(rolling_df['total_spillover'].mean(), color='gray', linestyle='--', linewidth=0.8, label='Full-period mean')417418events = {419 'COVID crash\n(Mar 2020)': pd.Timestamp('2020-03-15'),420 'Fed hikes\n(Mar 2022)': pd.Timestamp('2022-03-15'),421 'Crypto winter\n(Jun 2022)': pd.Timestamp('2022-06-01'),422 'FTX\n(Nov 2022)': pd.Timestamp('2022-11-11'),423}424colors_ev = ['#d62728', '#ff7f0e', '#9467bd', '#8c564b']425for (label, date), c in zip(events.items(), colors_ev):426 if rolling_df.index.min() <= date <= rolling_df.index.max():427 ax.axvline(date, color=c, linestyle='--', linewidth=1.2, alpha=0.8)428 ax.text(date, ax.get_ylim()[1] if ax.get_ylim()[1] > 0 else rolling_df['total_spillover'].max(),429 label, ha='center', va='top', fontsize=7.5, color=c, rotation=0,430 bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7, edgecolor=c, linewidth=0.5))431432ax.set_xlabel('Date')433ax.set_ylabel('Total Spillover Index (%)')434ax.set_title('Rolling 200-Day Diebold-Yilmaz Total Spillover Index (VAR-4, H=10)', fontweight='bold')435ax.legend(loc='upper left')436ax.xaxis.set_major_locator(matplotlib.dates.YearLocator(1))437ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))438plt.setp(ax.xaxis.get_majorticklabels(), rotation=45, ha='right')439plt.tight_layout()440fig.savefig(os.path.join(FIGURES, 'fig3_rolling_spillover.png'), dpi=150, bbox_inches='tight')441plt.close()442print(" Saved fig3_rolling_spillover.png")443444# ─── 5. Regime Analysis ───────────────────────────────────────────────────────445print("\n" + "="*70)446print("Section 5: Regime Analysis...")447448regimes = {449 'Pre-COVID': ('2019-01-01', '2020-02-28', 'table3a_connectedness_precovid.csv'),450 'COVID crash': ('2020-03-01', '2020-06-30', 'table3b_connectedness_covid.csv'),451 'Rate hikes': ('2022-01-01', '2022-12-31', 'table3c_connectedness_ratehikes.csv'),452}453454regime_tsi = {}455regime_net = {}456457for regime_name, (start, end, fname) in regimes.items():458 sub = rv_panel.loc[start:end].dropna()459 print(f"\n Regime: {regime_name} ({start} to {end}), {len(sub)} days")460 if len(sub) < 30:461 print(" WARNING: Too few observations, skipping.")462 continue463 # Use lag=2 if few obs464 lag = 4 if len(sub) >= 80 else 2465 tbl, tsi = compute_connectedness(sub, common_assets, lag=lag, H=10)466 if tbl is not None:467 tbl.to_csv(os.path.join(RESULTS, fname), float_format='%.2f')468 print(f" TSI = {tsi:.2f}%")469 regime_tsi[regime_name] = tsi470 regime_net[regime_name] = {a: tbl.loc['NET', a] for a in common_assets}471 print(f" NET spillovers: {regime_net[regime_name]}")472 print(f" Saved {fname}")473474# ─── Fig 4: NET spillovers by regime ─────────────────────────────────────────475print("\nPlotting Fig 4: NET spillovers by regime...")476if len(regime_net) > 0:477 fig, ax = plt.subplots(figsize=(12, 6))478 x = np.arange(len(common_assets))479 regime_list = list(regime_net.keys())480 n_regimes = len(regime_list)481 width = 0.8 / n_regimes482 regime_palette = ['#1f77b4', '#d62728', '#ff7f0e']483484 for k, (rname, c) in enumerate(zip(regime_list, regime_palette)):485 nets = [regime_net[rname].get(a, np.nan) for a in common_assets]486 offset = (k - (n_regimes - 1) / 2) * width487 ax.bar(x + offset, nets, width=width * 0.9, label=rname, color=c, alpha=0.8, edgecolor='white', linewidth=0.5)488489 ax.axhline(0, color='black', linewidth=0.8)490 ax.set_xticks(x)491 ax.set_xticklabels(common_assets)492 ax.set_xlabel('Asset')493 ax.set_ylabel('NET Spillover (%)')494 ax.set_title('NET Directional Spillovers by Market Regime', fontweight='bold')495 ax.legend(title='Regime')496 # Add TSI annotations497 for k, rname in enumerate(regime_list):498 if rname in regime_tsi:499 ax.text(0.02 + k * 0.18, 0.97, f'{rname}\nTSI={regime_tsi[rname]:.1f}%',500 transform=ax.transAxes, fontsize=8, va='top', color=regime_palette[k],501 bbox=dict(boxstyle='round', facecolor='white', alpha=0.7, edgecolor=regime_palette[k], linewidth=0.8))502 plt.tight_layout()503 fig.savefig(os.path.join(FIGURES, 'fig4_net_spillovers_regimes.png'), dpi=150, bbox_inches='tight')504 plt.close()505 print(" Saved fig4_net_spillovers_regimes.png")506507# ─── 6. Trading Hour Analysis ─────────────────────────────────────────────────508print("\n" + "="*70)509print("Section 6: Trading Hour Analysis...")510511session_assets = ['VIX', 'ES', 'NQ', 'BTC']512513# Load 5-min data for session analysis514def load_full_asset(db_name, symbol):515 db_path = os.path.join(DATA_DIR, db_name)516 conn = duckdb.connect(db_path, read_only=True)517 df = conn.execute(518 f"SELECT datetime, close FROM ohlcv WHERE symbol='{symbol}' ORDER BY datetime"519 ).fetchdf()520 conn.close()521 df['datetime'] = pd.to_datetime(df['datetime'])522 df = df.set_index('datetime')523 # Ensure UTC-naive (data is already UTC)524 return df525526session_data = {527 'ES': load_full_asset('futures_5min.duckdb', 'ES'),528 'NQ': load_full_asset('futures_5min.duckdb', 'NQ'),529 'VIX': load_full_asset('index_5min.duckdb', 'VIX'),530 'BTC': load_full_asset('crypto_5min.duckdb', 'BTC'),531}532533def classify_session(hour):534 if 0 <= hour < 8:535 return 'Asian'536 elif 8 <= hour < 14.5:537 return 'European'538 elif 14.5 <= hour < 21:539 return 'US'540 else:541 return 'Other'542543session_results = {}544545for sym, df in session_data.items():546 df = df.copy()547 df['log_ret'] = np.log(df['close'] / df['close'].shift(1))548 df['abs_ret'] = df['log_ret'].abs()549 df['hour'] = df.index.hour + df.index.minute / 60550 df['session'] = df['hour'].apply(classify_session)551 df['hour_int'] = df.index.hour552 session_results[sym] = df553554# Hourly average absolute returns555print("\n Hourly average absolute returns by session:")556hourly_avgs = {}557for sym, df in session_results.items():558 hav = df.groupby('hour_int')['abs_ret'].mean() * 100 # in %559 hourly_avgs[sym] = hav560561# Session-level statistics562session_stats = []563for sym, df in session_results.items():564 for sess in ['Asian', 'European', 'US']:565 sub = df[df['session'] == sess]566 if len(sub) > 10:567 session_stats.append({568 'Asset': sym,569 'Session': sess,570 'Mean_AbsRet': sub['abs_ret'].mean() * 100,571 'N_obs': len(sub),572 })573 print(f" {sym:6s} | {sess:10s}: Mean |ret| = {sub['abs_ret'].mean()*100:.4f}%")574575# Granger causality between VX/VIX, ES, NQ, BTC RV by session576print("\n Granger causality tests by session (VAR-4)...")577578# Use daily RV for granger tests (more robust)579gc_assets_rv = {a: rv_dict.get(a) for a in ['VIX', 'ES', 'NQ', 'BTC'] if a in rv_dict}580gc_panel = pd.DataFrame(gc_assets_rv).dropna()581582granger_rows = []583pairs = [('VIX', 'ES'), ('VIX', 'NQ'), ('VIX', 'BTC'),584 ('ES', 'NQ'), ('ES', 'BTC'), ('NQ', 'BTC'),585 ('ES', 'VIX'), ('NQ', 'VIX'), ('BTC', 'VIX')]586587# For session-specific granger, use 5-min abs returns aggregated to daily by session588for sess in ['Asian', 'European', 'US']:589 for (x_sym, y_sym) in pairs:590 if x_sym not in gc_assets_rv or y_sym not in gc_assets_rv:591 continue592 # Use full daily RV (session-specific daily RV would need more computation)593 try:594 pair_data = gc_panel[[y_sym, x_sym]].dropna()595 if len(pair_data) < 50:596 continue597 gc_test = grangercausalitytests(pair_data, maxlag=4, verbose=False)598 # F-stat from lag=4599 fstat = gc_test[4][0]['ssr_ftest'][0]600 pval = gc_test[4][0]['ssr_ftest'][1]601 granger_rows.append({602 'Session': sess,603 'Cause': x_sym,604 'Effect': y_sym,605 'F_stat': fstat,606 'p_value': pval,607 'Significant': 'Yes' if pval < 0.05 else 'No',608 })609 except Exception as e:610 pass611612granger_df = pd.DataFrame(granger_rows)613granger_df.to_csv(os.path.join(RESULTS, 'table4_session_granger.csv'), index=False, float_format='%.4f')614print(" Saved table4_session_granger.csv")615616# ─── Fig 5: Hourly Volatility ─────────────────────────────────────────────────617print("\nPlotting Fig 5: Hourly Volatility Pattern...")618fig, axes = plt.subplots(2, 2, figsize=(13, 8))619axes = axes.flatten()620sym_list = ['VIX', 'ES', 'NQ', 'BTC']621622for i, sym in enumerate(sym_list):623 ax = axes[i]624 hav = hourly_avgs[sym]625 color = ASSET_COLORS.get(sym, 'steelblue')626627 bars = ax.bar(hav.index, hav.values, color=color, alpha=0.8, edgecolor='white', linewidth=0.5)628629 # Session shading630 ax.axvspan(0, 8, alpha=0.07, color='blue', label='Asian (00-08 UTC)')631 ax.axvspan(8, 14.5, alpha=0.07, color='green', label='European (08-14:30 UTC)')632 ax.axvspan(14.5, 21, alpha=0.07, color='orange', label='US (14:30-21 UTC)')633634 ax.set_title(f'{sym} — Avg |5-min return| by Hour (UTC)', fontweight='bold')635 ax.set_xlabel('Hour of Day (UTC)')636 ax.set_ylabel('Mean |Return| (%)')637 ax.set_xticks(range(0, 24, 2))638 if i == 0:639 ax.legend(fontsize=7.5, loc='upper right')640641plt.suptitle('Intraday Volatility Pattern by Asset and Trading Session', fontsize=12, fontweight='bold')642plt.tight_layout()643fig.savefig(os.path.join(FIGURES, 'fig5_hourly_volatility.png'), dpi=150, bbox_inches='tight')644plt.close()645print(" Saved fig5_hourly_volatility.png")646647# ─── 7. Net Directional Spillover Evolution ───────────────────────────────────648print("\n" + "="*70)649print("Section 7: Rolling NET spillovers for BTC and ES...")650651rolling_net_dates = []652rolling_net_btc = []653rolling_net_es = []654655n_steps2 = len(dates_all) - WINDOW656for i in range(0, n_steps2, 1):657 window_data = rv_panel.iloc[i : i + WINDOW]658 end_date = dates_all[i + WINDOW - 1]659 sub = window_data.dropna()660 if len(sub) < 100:661 continue662663 model = VAR(sub)664 try:665 res = model.fit(4, trend='c')666 gfevd = compute_gfevd(res, H=10)667 n = gfevd.shape[0]668 off_diag = gfevd.copy()669 np.fill_diagonal(off_diag, 0.0)670 TO = off_diag.sum(axis=0) * 100671 FROM = off_diag.sum(axis=1) * 100672 NET = TO - FROM673674 idx_btc = common_assets.index('BTC') if 'BTC' in common_assets else None675 idx_es = common_assets.index('ES') if 'ES' in common_assets else None676677 if idx_btc is not None and idx_es is not None:678 rolling_net_dates.append(end_date)679 rolling_net_btc.append(NET[idx_btc])680 rolling_net_es.append(NET[idx_es])681 except Exception:682 pass683684 if i % 200 == 0 and rolling_net_btc:685 print(f" Step {i}/{n_steps2}: BTC NET={rolling_net_btc[-1]:.1f}%, ES NET={rolling_net_es[-1]:.1f}%")686687net_df = pd.DataFrame({688 'date': rolling_net_dates,689 'NET_BTC': rolling_net_btc,690 'NET_ES': rolling_net_es,691}).set_index('date')692693# ─── Fig 6: Rolling NET BTC vs ES ─────────────────────────────────────────────694print("\nPlotting Fig 6: Rolling NET directional spillovers BTC vs ES...")695fig, ax = plt.subplots(figsize=(13, 5))696ax.plot(net_df.index, net_df['NET_BTC'], color=ASSET_COLORS['BTC'], linewidth=1.1, label='BTC NET spillover', alpha=0.9)697ax.plot(net_df.index, net_df['NET_ES'], color=ASSET_COLORS['ES'], linewidth=1.1, label='ES NET spillover', alpha=0.9)698ax.axhline(0, color='black', linewidth=0.8, linestyle='--')699ax.fill_between(net_df.index, net_df['NET_BTC'], 0,700 where=(net_df['NET_BTC'] > 0), alpha=0.15, color=ASSET_COLORS['BTC'], label='BTC net transmitter')701ax.fill_between(net_df.index, net_df['NET_BTC'], 0,702 where=(net_df['NET_BTC'] < 0), alpha=0.15, color='gray', label='BTC net receiver')703704for (label, date), c in zip(events.items(), colors_ev):705 if net_df.index.min() <= date <= net_df.index.max():706 ax.axvline(date, color=c, linestyle=':', linewidth=1, alpha=0.8)707 ax.text(date, ax.get_ylim()[1] if ax.get_ylim()[1] else net_df[['NET_BTC','NET_ES']].max().max(),708 label, ha='center', va='top', fontsize=7, color=c,709 bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.6))710711ax.set_xlabel('Date')712ax.set_ylabel('NET Directional Spillover (%)')713ax.set_title('Rolling 200-Day NET Directional Spillovers: BTC vs ES (VAR-4, H=10)', fontweight='bold')714ax.legend(loc='lower right', ncol=2)715ax.xaxis.set_major_locator(matplotlib.dates.YearLocator(1))716ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))717plt.setp(ax.xaxis.get_majorticklabels(), rotation=45, ha='right')718plt.tight_layout()719fig.savefig(os.path.join(FIGURES, 'fig6_net_directional_btc_es.png'), dpi=150, bbox_inches='tight')720plt.close()721print(" Saved fig6_net_directional_btc_es.png")722723# ─── 8. Network Visualization ─────────────────────────────────────────────────724print("\n" + "="*70)725print("Section 8: Network Graph...")726727if conn_table is not None:728 gfevd_mat = conn_table.loc[common_assets, common_assets].values.astype(float)729 np.fill_diagonal(gfevd_mat, 0.0)730731 # Threshold: keep edges above median732 all_edges = gfevd_mat.flatten()733 all_edges = all_edges[all_edges > 0]734 threshold = np.median(all_edges)735736 G = nx.DiGraph()737 for a in common_assets:738 G.add_node(a, asset_class=ASSET_CLASS.get(a, 'Other'))739740 TO_vals = {}741 for j, asj in enumerate(common_assets):742 TO_vals[asj] = sum(gfevd_mat[i, j] for i in range(len(common_assets)) if i != j)743744 for i, asi in enumerate(common_assets):745 for j, asj in enumerate(common_assets):746 if i != j and gfevd_mat[i, j] > threshold:747 G.add_edge(asj, asi, weight=gfevd_mat[i, j]) # j → i (j causes forecast error in i)748749 # Layout750 pos = nx.spring_layout(G, k=2.5, seed=42, iterations=100)751752 node_sizes = [max(300, TO_vals.get(a, 1) * 40) for a in G.nodes()]753 node_colors = [CLASS_COLOR.get(ASSET_CLASS.get(a, 'Other'), 'gray') for a in G.nodes()]754 edge_weights = [G[u][v]['weight'] for u, v in G.edges()]755 max_w = max(edge_weights) if edge_weights else 1756 edge_widths = [w / max_w * 4 for w in edge_weights]757 edge_alphas = [0.3 + 0.5 * w / max_w for w in edge_weights]758759 fig, ax = plt.subplots(figsize=(12, 9))760761 for (u, v), lw, alpha, w in zip(G.edges(), edge_widths, edge_alphas, edge_weights):762 nx.draw_networkx_edges(G, pos, edgelist=[(u, v)], ax=ax,763 width=lw, alpha=alpha, edge_color='#555555',764 arrows=True, arrowsize=12,765 connectionstyle='arc3,rad=0.15',766 node_size=node_sizes[list(G.nodes()).index(u)] if u in G.nodes() else 300)767768 nx.draw_networkx_nodes(G, pos, ax=ax,769 node_size=node_sizes,770 node_color=node_colors,771 alpha=0.9, edgecolors='white', linewidths=1.5)772 nx.draw_networkx_labels(G, pos, ax=ax, font_size=9, font_weight='bold', font_color='black')773774 # Legend for asset classes775 legend_patches = [mpatches.Patch(color=c, label=cls)776 for cls, c in CLASS_COLOR.items()]777 ax.legend(handles=legend_patches, title='Asset Class', loc='upper left', fontsize=9)778 ax.set_title(f'Volatility Transmission Network\n(Edges > median threshold = {threshold:.1f}%,\n'779 f'Node size ∝ Total TO spillover, Arrow direction = shock transmission)',780 fontweight='bold', fontsize=10)781 ax.axis('off')782 plt.tight_layout()783 fig.savefig(os.path.join(FIGURES, 'fig7_network_graph.png'), dpi=150, bbox_inches='tight')784 plt.close()785 print(" Saved fig7_network_graph.png")786787# ─── Summary Report ───────────────────────────────────────────────────────────788print("\n" + "=" * 70)789print("ANALYSIS COMPLETE — KEY NUMERICAL FINDINGS")790print("=" * 70)791792print("\n[Table 1] Descriptive Statistics (5-min log returns):")793for row in desc_rows:794 print(f" {row['Asset']:10s}: N={row['N_obs']:>8,}, Std={row['Std']:.6f}, "795 f"Skew={row['Skewness']:+.3f}, Kurt={row['Kurtosis']:+.2f}")796797print(f"\n[Section 2] Daily RV — Annualized means:")798for sym, rv in rv_dict.items():799 print(f" {sym:10s}: mean RV = {rv.mean():.4f}, max RV = {rv.max():.4f}")800801print(f"\n[Section 3] Full-Sample Connectedness:")802print(f" Total Spillover Index = {tsi_full:.2f}%" if tsi_full else " N/A")803if conn_table is not None:804 print(" NET spillovers (positive = net transmitter):")805 for a in common_assets:806 net_val = conn_table.loc['NET', a]807 to_val = conn_table.loc['TO', a]808 frm_val = conn_table.loc[a, 'FROM']809 print(f" {a:10s}: TO={to_val:+.2f}%, FROM={frm_val:+.2f}%, NET={net_val:+.2f}%")810811print(f"\n[Section 4] Rolling Connectedness ({len(rolling_df)} windows):")812print(f" Mean TSI: {rolling_df['total_spillover'].mean():.2f}%")813print(f" Max TSI: {rolling_df['total_spillover'].max():.2f}% on {rolling_df['total_spillover'].idxmax().date()}")814print(f" Min TSI: {rolling_df['total_spillover'].min():.2f}% on {rolling_df['total_spillover'].idxmin().date()}")815816print(f"\n[Section 5] Regime TSI:")817for r, tsi in regime_tsi.items():818 print(f" {r:20s}: TSI = {tsi:.2f}%")819820print(f"\n[Section 6] Granger causality (sample):")821sig_gc = granger_df[granger_df['Significant'] == 'Yes'] if len(granger_df) > 0 else pd.DataFrame()822print(f" {len(sig_gc)} / {len(granger_df)} pairs significant at 5%")823if len(sig_gc) > 0:824 print(sig_gc[['Session', 'Cause', 'Effect', 'F_stat', 'p_value']].head(10).to_string(index=False))825826print(f"\n[Section 7] Rolling NET BTC (last window): {rolling_net_btc[-1]:.2f}%" if rolling_net_btc else "")827print(f" Rolling NET ES (last window): {rolling_net_es[-1]:.2f}%" if rolling_net_es else "")828829print("\nAll outputs saved to:", RESULTS)830print("=" * 70)831