TeX 62.4%
Python 37.5%
1"""2Intraday Volatility Transmission Across Asset Classes3Diebold-Yilmaz Connectedness Analysis4"""56import os7import warnings8import numpy as np9import pandas as pd10import matplotlib11matplotlib.use('Agg')12import matplotlib.pyplot as plt13import matplotlib.patches as mpatches14import matplotlib.gridspec as gridspec15import seaborn as sns16import duckdb17from scipy import stats18import networkx as nx19from statsmodels.tsa.api import VAR20from statsmodels.tsa.stattools import grangercausalitytests21import traceback2223warnings.filterwarnings('ignore')2425# ─── Paths ────────────────────────────────────────────────────────────────────26ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))27DATA_DIR = os.environ.get('WP6_DATA_DIR', os.path.join(ROOT, 'data', 'raw')) # *_5min.duckdb files (see data/raw/README.md)28RESULTS = os.path.join(ROOT, 'results')29FIGURES = os.path.join(ROOT, 'figures')30os.makedirs(FIGURES, exist_ok=True)31os.makedirs(RESULTS, exist_ok=True)3233plt.rcParams.update({34 'font.size': 10, 'axes.titlesize': 11, 'axes.labelsize': 10,35 'xtick.labelsize': 9, 'ytick.labelsize': 9, 'legend.fontsize': 9,36 'figure.dpi': 150, 'savefig.dpi': 150, 'font.family': 'DejaVu Sans',37 'axes.spines.top': False, 'axes.spines.right': False,38})3940ASSET_COLORS = {41 'ES': '#1f77b4', 'NQ': '#aec7e8',42 'EURUSD': '#2ca02c', 'USDJPY': '#98df8a', 'GBPUSD': '#d62728',43 'GC': '#ff7f0e', 'CL': '#ffbb78',44 'BTC': '#9467bd', 'ETH': '#c5b0d5',45 'VIX': '#8c564b',46}47ASSET_CLASS = {48 'ES': 'Equity', 'NQ': 'Equity',49 'EURUSD': 'FX', 'USDJPY': 'FX', 'GBPUSD': 'FX',50 'GC': 'Commodity', 'CL': 'Commodity',51 'BTC': 'Crypto', 'ETH': 'Crypto',52 'VIX': 'Volatility',53}54CLASS_COLOR = {55 'Equity': '#1f77b4', 'FX': '#2ca02c',56 'Commodity': '#ff7f0e', 'Crypto': '#9467bd', 'Volatility': '#8c564b',57}5859ASSETS = ['ES', 'NQ', 'EURUSD', 'USDJPY', 'GBPUSD', 'GC', 'CL', 'BTC', 'ETH', 'VIX']6061# ─── 1. Load 5-min data ───────────────────────────────────────────────────────62print("=" * 70)63print("Loading 5-min OHLCV data...")6465def load_asset(db_name, symbol):66 db_path = os.path.join(DATA_DIR, db_name)67 conn = duckdb.connect(db_path, read_only=True)68 df = conn.execute(69 f"SELECT datetime, close FROM ohlcv WHERE symbol='{symbol}' ORDER BY datetime"70 ).fetchdf()71 conn.close()72 df['datetime'] = pd.to_datetime(df['datetime'])73 df = df.set_index('datetime')74 df.columns = [symbol]75 return df7677futures_assets = {78 'ES': 'futures_5min.duckdb',79 'NQ': 'futures_5min.duckdb',80 'GC': 'futures_5min.duckdb',81 'CL': 'futures_5min.duckdb',82 'VX': 'futures_5min.duckdb',83}84fx_assets = {85 'EURUSD': 'fx_5min.duckdb',86 'USDJPY': 'fx_5min.duckdb',87 'GBPUSD': 'fx_5min.duckdb',88}89crypto_assets = {90 'BTC': 'crypto_5min.duckdb',91 'ETH': 'crypto_5min.duckdb',92}93index_assets = {94 'VIX': 'index_5min.duckdb',95}9697raw_data = {}98for sym, db in {**futures_assets, **fx_assets, **crypto_assets, **index_assets}.items():99 print(f" Loading {sym} from {db}...")100 raw_data[sym] = load_asset(db, sym)101102print("Data loaded.\n")103104# ─── Compute 5-min log returns ────────────────────────────────────────────────105print("Computing 5-min log returns...")106log_returns = {}107for sym, df in raw_data.items():108 lr = np.log(df[sym] / df[sym].shift(1)).dropna()109 # Remove extreme outliers (> 20 std)110 mu, sigma = lr.mean(), lr.std()111 lr = lr[(lr - mu).abs() < 20 * sigma]112 log_returns[sym] = lr113114# ─── 1. Descriptive Statistics ────────────────────────────────────────────────115print("Computing descriptive statistics...")116desc_rows = []117for sym in ASSET_CLASS.keys():118 if sym not in log_returns:119 # VX → treat VX as proxy for VIX120 continue121 lr = log_returns[sym]122 row = {123 'Asset': sym,124 'Class': ASSET_CLASS.get(sym, 'Other'),125 'N_obs': len(lr),126 'Date_Start': lr.index.min().strftime('%Y-%m-%d'),127 'Date_End': lr.index.max().strftime('%Y-%m-%d'),128 'Mean': lr.mean(),129 'Std': lr.std(),130 'Skewness': stats.skew(lr),131 'Kurtosis': stats.kurtosis(lr),132 'Min': lr.min(),133 'Max': lr.max(),134 }135 desc_rows.append(row)136137desc_df = pd.DataFrame(desc_rows)138desc_df.to_csv(os.path.join(RESULTS, 'table1_descriptive_stats.csv'), index=False, float_format='%.6f')139print(" Saved table1_descriptive_stats.csv")140print(desc_df[['Asset', 'N_obs', 'Mean', 'Std', 'Skewness', 'Kurtosis']].to_string(index=False))141142# ─── 2. Realized Volatility ───────────────────────────────────────────────────143print("\nComputing daily realized volatility...")144145def compute_rv(log_ret_series, min_obs=20):146 """Annualized daily realized volatility."""147 df = log_ret_series.to_frame(name='ret')148 df['date'] = df.index.date149 grouped = df.groupby('date')['ret']150 rv = grouped.apply(lambda x: np.sqrt(np.sum(x**2)) * np.sqrt(252) if len(x) >= min_obs else np.nan)151 rv.index = pd.to_datetime(rv.index)152 return rv.dropna()153154rv_dict = {}155for sym in ASSET_CLASS.keys():156 if sym not in log_returns:157 continue158 rv_dict[sym] = compute_rv(log_returns[sym])159 print(f" {sym}: {len(rv_dict[sym])} trading days, RV mean = {rv_dict[sym].mean():.4f}")160161# Combine into DataFrame162rv_df = pd.DataFrame(rv_dict)163rv_df.index.name = 'date'164rv_df.to_csv(os.path.join(RESULTS, 'daily_rv_all.csv'))165print(" Saved daily_rv_all.csv")166167# ─── Fig 1: RV Time Series ────────────────────────────────────────────────────168print("\nPlotting Fig 1: RV Time Series...")169fig, axes = plt.subplots(5, 2, figsize=(14, 16), sharex=False)170axes = axes.flatten()171172asset_order = ['ES', 'NQ', 'EURUSD', 'USDJPY', 'GBPUSD', 'GC', 'CL', 'BTC', 'ETH', 'VIX']173for i, sym in enumerate(asset_order):174 if sym not in rv_dict:175 continue176 ax = axes[i]177 rv_s = rv_dict[sym]178 color = ASSET_COLORS.get(sym, 'steelblue')179 ax.plot(rv_s.index, rv_s.values, color=color, linewidth=0.7, alpha=0.9)180 ax.fill_between(rv_s.index, rv_s.values, alpha=0.15, color=color)181 ax.set_title(f'{sym} ({ASSET_CLASS.get(sym, "")})')182 ax.set_ylabel('Annualized RV')183 ax.xaxis.set_major_locator(matplotlib.dates.YearLocator(2))184 ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))185 plt.setp(ax.xaxis.get_majorticklabels(), rotation=45, ha='right')186187plt.suptitle('Daily Realized Volatility — All Assets (5-min, Annualized)', fontsize=13, y=1.01)188plt.tight_layout()189fig.savefig(os.path.join(FIGURES, 'fig1_rv_timeseries.png'), dpi=150, bbox_inches='tight')190plt.close()191print(" Saved fig1_rv_timeseries.png")192193# ─── Helper: Diebold-Yilmaz GFEVD ────────────────────────────────────────────194def compute_gfevd(var_result, H=10):195 """196 Generalized Forecast Error Variance Decomposition (Pesaran-Shin).197 Returns an (n x n) matrix where entry [i,j] = share of forecast error198 variance of variable i attributable to shocks from variable j.199 """200 n = var_result.k_ar # actually number of variables201 coefs = var_result.coefs # shape (p, n, n)202 sigma = var_result.sigma_u # (n, n) residual covariance203 n = sigma.shape[0]204 p = coefs.shape[0]205206 # Build MA coefficient matrices Psi_h via recursive computation207 # Psi_0 = I, Psi_h = sum_{j=1..min(h,p)} A_j @ Psi_{h-j}208 psi = [np.eye(n)]209 for h in range(1, H + 1):210 psi_h = np.zeros((n, n))211 for j in range(1, min(h, p) + 1):212 psi_h += coefs[j-1] @ psi[h-j]213 psi.append(psi_h)214215 # GFEVD: theta[i,j] = (1/sigma[j,j]) * sum_{h=0}^{H-1} (e_i' Psi_h Sigma e_j)^2216 # / sum_{h=0}^{H-1} (e_i' Psi_h Sigma Psi_h' e_i)217 sigma_diag = np.diag(sigma)218 gfevd = np.zeros((n, n))219220 for i in range(n):221 ei = np.zeros(n); ei[i] = 1.0222 denom = 0.0223 for h in range(H):224 denom += ei @ psi[h] @ sigma @ psi[h].T @ ei225 for j in range(n):226 ej = np.zeros(n); ej[j] = 1.0227 num = 0.0228 for h in range(H):229 num += (ei @ psi[h] @ sigma @ ej) ** 2230 gfevd[i, j] = num / (sigma_diag[j] * denom)231232 # Row-normalize so rows sum to 1233 row_sums = gfevd.sum(axis=1, keepdims=True)234 row_sums[row_sums == 0] = 1.0235 gfevd = gfevd / row_sums236 return gfevd237238239def compute_connectedness(rv_panel, assets, lag=4, H=10):240 """241 Given a DataFrame of RV series (rows=dates, cols=assets),242 fit VAR(lag), compute GFEVD at horizon H, return connectedness table.243 """244 # Use only complete rows245 data = rv_panel[assets].dropna()246 if len(data) < lag + H + 10:247 return None, None248249 model = VAR(data)250 try:251 res = model.fit(lag, trend='c')252 except Exception as e:253 print(f" VAR fit failed: {e}")254 return None, None255256 gfevd = compute_gfevd(res, H=H)257 n = len(assets)258259 # Build full connectedness table260 # Rows = "TO" (recipient of shocks), Cols = "FROM" (origin of shocks)261 # But DY convention: entry [i,j] = contribution of j to i's forecast error262 # TO_j = sum_{i != j} theta[i,j] / N * 100 (what j contributes to others)263 # FROM_i = sum_{j != i} theta[i,j] / N * 100 (what i receives from others)264265 # gfevd[i,j] = share of i's FEVD from j's shock266 # So: FROM_i = sum_{j!=i} gfevd[i,j] * 100267 # TO_j = sum_{i!=j} gfevd[i,j] * 100268269 # Total spillover index270 off_diag = gfevd.copy()271 np.fill_diagonal(off_diag, 0.0)272 TSI = off_diag.sum() / n * 100273274 # Directional275 FROM = off_diag.sum(axis=1) * 100 # what each asset receives276 TO = off_diag.sum(axis=0) * 100 # what each asset transmits277 NET = TO - FROM278279 # Build output DataFrame (DY table format)280 # gfevd * 100 → percentages281 tbl = pd.DataFrame(gfevd * 100, index=assets, columns=assets)282 tbl['FROM'] = FROM283 to_row = pd.Series(TO, index=assets)284 net_row = pd.Series(NET, index=assets)285 tbl.loc['TO'] = to_row286 tbl.loc['NET'] = net_row287 tbl.loc['TO', 'FROM'] = TSI288 tbl.loc['NET', 'FROM'] = np.nan289290 return tbl, TSI291292293# ─── 3. Full-Sample Connectedness ─────────────────────────────────────────────294print("\n" + "="*70)295print("Section 3: Full-Sample Connectedness (VAR-4, H=10)...")296297# Use VIX from index (not VX futures) as VIX is more common298# Find common sample: align all RV series299common_assets = [a for a in ASSET_CLASS.keys() if a in rv_dict]300rv_panel = pd.DataFrame({a: rv_dict[a] for a in common_assets})301302# Forward fill up to 5 days to handle weekends/holidays303rv_panel = rv_panel.ffill(limit=5)304305# Drop rows where more than 3 assets are missing306rv_panel = rv_panel.dropna(thresh=len(common_assets) - 2)307rv_panel = rv_panel.dropna()308309print(f" Common sample: {rv_panel.index[0].date()} to {rv_panel.index[-1].date()}, {len(rv_panel)} days")310print(f" Assets: {common_assets}")311312conn_table, tsi_full = compute_connectedness(rv_panel, common_assets, lag=4, H=10)313if conn_table is not None:314 conn_table.to_csv(os.path.join(RESULTS, 'table2_connectedness_full.csv'), float_format='%.2f')315 print(f" Total Spillover Index (full sample): {tsi_full:.2f}%")316 print(" NET spillovers:")317 for a in common_assets:318 print(f" {a:10s}: NET = {conn_table.loc['NET', a]:+.2f}%")319 print(" Saved table2_connectedness_full.csv")320else:321 print(" WARNING: Full-sample connectedness failed.")322323# ─── Fig 2: Connectedness Heatmap ─────────────────────────────────────────────324if conn_table is not None:325 print("\nPlotting Fig 2: Connectedness Heatmap...")326 heat_data = conn_table.loc[common_assets, common_assets].values.astype(float)327 fig, ax = plt.subplots(figsize=(11, 8))328 mask_diag = np.eye(len(common_assets), dtype=bool)329 im = ax.imshow(heat_data, cmap='RdYlBu_r', aspect='auto', vmin=0, vmax=heat_data.max())330 for i in range(len(common_assets)):331 for j in range(len(common_assets)):332 txt = f"{heat_data[i,j]:.1f}"333 color = 'white' if heat_data[i,j] > heat_data.max() * 0.6 else 'black'334 ax.text(j, i, txt, ha='center', va='center', fontsize=8, color=color, fontweight='bold')335 ax.set_xticks(range(len(common_assets)))336 ax.set_yticks(range(len(common_assets)))337 ax.set_xticklabels(common_assets, rotation=45, ha='right')338 ax.set_yticklabels(common_assets)339 ax.set_xlabel('FROM (shock origin)')340 ax.set_ylabel('TO (shock recipient)')341 ax.set_title(f'Diebold-Yilmaz Connectedness Matrix (Full Sample)\nTotal Spillover Index = {tsi_full:.1f}%',342 fontweight='bold')343 plt.colorbar(im, ax=ax, label='% of Forecast Error Variance')344 plt.tight_layout()345 fig.savefig(os.path.join(FIGURES, 'fig2_connectedness_heatmap.png'), dpi=150, bbox_inches='tight')346 plt.close()347 print(" Saved fig2_connectedness_heatmap.png")348349# ─── 4. Rolling Connectedness ─────────────────────────────────────────────────350print("\n" + "="*70)351print("Section 4: Rolling Connectedness (200-day window, step=1)...")352353WINDOW = 200354dates_all = rv_panel.index.tolist()355rolling_dates = []356rolling_tsi = []357358n_steps = len(dates_all) - WINDOW359print(f" Total windows to compute: {n_steps}")360361for i in range(0, n_steps, 1):362 window_data = rv_panel.iloc[i : i + WINDOW]363 end_date = dates_all[i + WINDOW - 1]364365 # Only compute if at least 150 non-NaN rows after dropna366 sub = window_data.dropna()367 if len(sub) < 100:368 continue369370 model = VAR(sub)371 try:372 res = model.fit(4, trend='c')373 gfevd = compute_gfevd(res, H=10)374 n = gfevd.shape[0]375 off_diag = gfevd.copy()376 np.fill_diagonal(off_diag, 0.0)377 tsi = off_diag.sum() / n * 100378 rolling_dates.append(end_date)379 rolling_tsi.append(tsi)380 except Exception:381 pass382383 if i % 200 == 0:384 print(f" Step {i}/{n_steps}: {end_date.date()}, TSI = {rolling_tsi[-1]:.1f}%" if rolling_tsi else f" Step {i}/{n_steps}")385386rolling_df = pd.DataFrame({'date': rolling_dates, 'total_spillover': rolling_tsi})387rolling_df = rolling_df.set_index('date')388rolling_df.to_csv(os.path.join(RESULTS, 'rolling_total_spillover.csv'))389print(f" Rolling analysis complete: {len(rolling_df)} observations")390print(f" Mean TSI: {rolling_df['total_spillover'].mean():.2f}%")391print(f" Min TSI: {rolling_df['total_spillover'].min():.2f}%")392print(f" Max TSI: {rolling_df['total_spillover'].max():.2f}%")393394# ─── Fig 3: Rolling Spillover ─────────────────────────────────────────────────395print("\nPlotting Fig 3: Rolling Total Spillover...")396fig, ax = plt.subplots(figsize=(13, 5))397ax.plot(rolling_df.index, rolling_df['total_spillover'], color='#1f77b4', linewidth=1.2, label='Total Spillover Index')398ax.fill_between(rolling_df.index, rolling_df['total_spillover'],399 rolling_df['total_spillover'].mean(), alpha=0.15, color='#1f77b4')400ax.axhline(rolling_df['total_spillover'].mean(), color='gray', linestyle='--', linewidth=0.8, label='Full-period mean')401402events = {403 'COVID crash\n(Mar 2020)': pd.Timestamp('2020-03-15'),404 'Fed hikes\n(Mar 2022)': pd.Timestamp('2022-03-15'),405 'Crypto winter\n(Jun 2022)': pd.Timestamp('2022-06-01'),406 'FTX\n(Nov 2022)': pd.Timestamp('2022-11-11'),407}408colors_ev = ['#d62728', '#ff7f0e', '#9467bd', '#8c564b']409for (label, date), c in zip(events.items(), colors_ev):410 if rolling_df.index.min() <= date <= rolling_df.index.max():411 ax.axvline(date, color=c, linestyle='--', linewidth=1.2, alpha=0.8)412 ax.text(date, ax.get_ylim()[1] if ax.get_ylim()[1] > 0 else rolling_df['total_spillover'].max(),413 label, ha='center', va='top', fontsize=7.5, color=c, rotation=0,414 bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7, edgecolor=c, linewidth=0.5))415416ax.set_xlabel('Date')417ax.set_ylabel('Total Spillover Index (%)')418ax.set_title('Rolling 200-Day Diebold-Yilmaz Total Spillover Index (VAR-4, H=10)', fontweight='bold')419ax.legend(loc='upper left')420ax.xaxis.set_major_locator(matplotlib.dates.YearLocator(1))421ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))422plt.setp(ax.xaxis.get_majorticklabels(), rotation=45, ha='right')423plt.tight_layout()424fig.savefig(os.path.join(FIGURES, 'fig3_rolling_spillover.png'), dpi=150, bbox_inches='tight')425plt.close()426print(" Saved fig3_rolling_spillover.png")427428# ─── 5. Regime Analysis ───────────────────────────────────────────────────────429print("\n" + "="*70)430print("Section 5: Regime Analysis...")431432regimes = {433 'Pre-COVID': ('2019-01-01', '2020-02-28', 'table3a_connectedness_precovid.csv'),434 'COVID crash': ('2020-03-01', '2020-06-30', 'table3b_connectedness_covid.csv'),435 'Rate hikes': ('2022-01-01', '2022-12-31', 'table3c_connectedness_ratehikes.csv'),436}437438regime_tsi = {}439regime_net = {}440441for regime_name, (start, end, fname) in regimes.items():442 sub = rv_panel.loc[start:end].dropna()443 print(f"\n Regime: {regime_name} ({start} to {end}), {len(sub)} days")444 if len(sub) < 30:445 print(" WARNING: Too few observations, skipping.")446 continue447 # Use lag=2 if few obs448 lag = 4 if len(sub) >= 80 else 2449 tbl, tsi = compute_connectedness(sub, common_assets, lag=lag, H=10)450 if tbl is not None:451 tbl.to_csv(os.path.join(RESULTS, fname), float_format='%.2f')452 print(f" TSI = {tsi:.2f}%")453 regime_tsi[regime_name] = tsi454 regime_net[regime_name] = {a: tbl.loc['NET', a] for a in common_assets}455 print(f" NET spillovers: {regime_net[regime_name]}")456 print(f" Saved {fname}")457458# ─── Fig 4: NET spillovers by regime ─────────────────────────────────────────459print("\nPlotting Fig 4: NET spillovers by regime...")460if len(regime_net) > 0:461 fig, ax = plt.subplots(figsize=(12, 6))462 x = np.arange(len(common_assets))463 regime_list = list(regime_net.keys())464 n_regimes = len(regime_list)465 width = 0.8 / n_regimes466 regime_palette = ['#1f77b4', '#d62728', '#ff7f0e']467468 for k, (rname, c) in enumerate(zip(regime_list, regime_palette)):469 nets = [regime_net[rname].get(a, np.nan) for a in common_assets]470 offset = (k - (n_regimes - 1) / 2) * width471 ax.bar(x + offset, nets, width=width * 0.9, label=rname, color=c, alpha=0.8, edgecolor='white', linewidth=0.5)472473 ax.axhline(0, color='black', linewidth=0.8)474 ax.set_xticks(x)475 ax.set_xticklabels(common_assets)476 ax.set_xlabel('Asset')477 ax.set_ylabel('NET Spillover (%)')478 ax.set_title('NET Directional Spillovers by Market Regime', fontweight='bold')479 ax.legend(title='Regime')480 # Add TSI annotations481 for k, rname in enumerate(regime_list):482 if rname in regime_tsi:483 ax.text(0.02 + k * 0.18, 0.97, f'{rname}\nTSI={regime_tsi[rname]:.1f}%',484 transform=ax.transAxes, fontsize=8, va='top', color=regime_palette[k],485 bbox=dict(boxstyle='round', facecolor='white', alpha=0.7, edgecolor=regime_palette[k], linewidth=0.8))486 plt.tight_layout()487 fig.savefig(os.path.join(FIGURES, 'fig4_net_spillovers_regimes.png'), dpi=150, bbox_inches='tight')488 plt.close()489 print(" Saved fig4_net_spillovers_regimes.png")490491# ─── 6. Trading Hour Analysis ─────────────────────────────────────────────────492print("\n" + "="*70)493print("Section 6: Trading Hour Analysis...")494495session_assets = ['VIX', 'ES', 'NQ', 'BTC']496497# Load 5-min data for session analysis498def load_full_asset(db_name, symbol):499 db_path = os.path.join(DATA_DIR, db_name)500 conn = duckdb.connect(db_path, read_only=True)501 df = conn.execute(502 f"SELECT datetime, close FROM ohlcv WHERE symbol='{symbol}' ORDER BY datetime"503 ).fetchdf()504 conn.close()505 df['datetime'] = pd.to_datetime(df['datetime'])506 df = df.set_index('datetime')507 # Ensure UTC-naive (data is already UTC)508 return df509510session_data = {511 'ES': load_full_asset('futures_5min.duckdb', 'ES'),512 'NQ': load_full_asset('futures_5min.duckdb', 'NQ'),513 'VIX': load_full_asset('index_5min.duckdb', 'VIX'),514 'BTC': load_full_asset('crypto_5min.duckdb', 'BTC'),515}516517def classify_session(hour):518 if 0 <= hour < 8:519 return 'Asian'520 elif 8 <= hour < 14.5:521 return 'European'522 elif 14.5 <= hour < 21:523 return 'US'524 else:525 return 'Other'526527session_results = {}528529for sym, df in session_data.items():530 df = df.copy()531 df['log_ret'] = np.log(df['close'] / df['close'].shift(1))532 df['abs_ret'] = df['log_ret'].abs()533 df['hour'] = df.index.hour + df.index.minute / 60534 df['session'] = df['hour'].apply(classify_session)535 df['hour_int'] = df.index.hour536 session_results[sym] = df537538# Hourly average absolute returns539print("\n Hourly average absolute returns by session:")540hourly_avgs = {}541for sym, df in session_results.items():542 hav = df.groupby('hour_int')['abs_ret'].mean() * 100 # in %543 hourly_avgs[sym] = hav544545# Session-level statistics546session_stats = []547for sym, df in session_results.items():548 for sess in ['Asian', 'European', 'US']:549 sub = df[df['session'] == sess]550 if len(sub) > 10:551 session_stats.append({552 'Asset': sym,553 'Session': sess,554 'Mean_AbsRet': sub['abs_ret'].mean() * 100,555 'N_obs': len(sub),556 })557 print(f" {sym:6s} | {sess:10s}: Mean |ret| = {sub['abs_ret'].mean()*100:.4f}%")558559# Granger causality between VX/VIX, ES, NQ, BTC RV by session560print("\n Granger causality tests by session (VAR-4)...")561562# Use daily RV for granger tests (more robust)563gc_assets_rv = {a: rv_dict.get(a) for a in ['VIX', 'ES', 'NQ', 'BTC'] if a in rv_dict}564gc_panel = pd.DataFrame(gc_assets_rv).dropna()565566granger_rows = []567pairs = [('VIX', 'ES'), ('VIX', 'NQ'), ('VIX', 'BTC'),568 ('ES', 'NQ'), ('ES', 'BTC'), ('NQ', 'BTC'),569 ('ES', 'VIX'), ('NQ', 'VIX'), ('BTC', 'VIX')]570571# For session-specific granger, use 5-min abs returns aggregated to daily by session572for sess in ['Asian', 'European', 'US']:573 for (x_sym, y_sym) in pairs:574 if x_sym not in gc_assets_rv or y_sym not in gc_assets_rv:575 continue576 # Use full daily RV (session-specific daily RV would need more computation)577 try:578 pair_data = gc_panel[[y_sym, x_sym]].dropna()579 if len(pair_data) < 50:580 continue581 gc_test = grangercausalitytests(pair_data, maxlag=4, verbose=False)582 # F-stat from lag=4583 fstat = gc_test[4][0]['ssr_ftest'][0]584 pval = gc_test[4][0]['ssr_ftest'][1]585 granger_rows.append({586 'Session': sess,587 'Cause': x_sym,588 'Effect': y_sym,589 'F_stat': fstat,590 'p_value': pval,591 'Significant': 'Yes' if pval < 0.05 else 'No',592 })593 except Exception as e:594 pass595596granger_df = pd.DataFrame(granger_rows)597granger_df.to_csv(os.path.join(RESULTS, 'table4_session_granger.csv'), index=False, float_format='%.4f')598print(" Saved table4_session_granger.csv")599600# ─── Fig 5: Hourly Volatility ─────────────────────────────────────────────────601print("\nPlotting Fig 5: Hourly Volatility Pattern...")602fig, axes = plt.subplots(2, 2, figsize=(13, 8))603axes = axes.flatten()604sym_list = ['VIX', 'ES', 'NQ', 'BTC']605606for i, sym in enumerate(sym_list):607 ax = axes[i]608 hav = hourly_avgs[sym]609 color = ASSET_COLORS.get(sym, 'steelblue')610611 bars = ax.bar(hav.index, hav.values, color=color, alpha=0.8, edgecolor='white', linewidth=0.5)612613 # Session shading614 ax.axvspan(0, 8, alpha=0.07, color='blue', label='Asian (00-08 UTC)')615 ax.axvspan(8, 14.5, alpha=0.07, color='green', label='European (08-14:30 UTC)')616 ax.axvspan(14.5, 21, alpha=0.07, color='orange', label='US (14:30-21 UTC)')617618 ax.set_title(f'{sym} — Avg |5-min return| by Hour (UTC)', fontweight='bold')619 ax.set_xlabel('Hour of Day (UTC)')620 ax.set_ylabel('Mean |Return| (%)')621 ax.set_xticks(range(0, 24, 2))622 if i == 0:623 ax.legend(fontsize=7.5, loc='upper right')624625plt.suptitle('Intraday Volatility Pattern by Asset and Trading Session', fontsize=12, fontweight='bold')626plt.tight_layout()627fig.savefig(os.path.join(FIGURES, 'fig5_hourly_volatility.png'), dpi=150, bbox_inches='tight')628plt.close()629print(" Saved fig5_hourly_volatility.png")630631# ─── 7. Net Directional Spillover Evolution ───────────────────────────────────632print("\n" + "="*70)633print("Section 7: Rolling NET spillovers for BTC and ES...")634635rolling_net_dates = []636rolling_net_btc = []637rolling_net_es = []638639n_steps2 = len(dates_all) - WINDOW640for i in range(0, n_steps2, 1):641 window_data = rv_panel.iloc[i : i + WINDOW]642 end_date = dates_all[i + WINDOW - 1]643 sub = window_data.dropna()644 if len(sub) < 100:645 continue646647 model = VAR(sub)648 try:649 res = model.fit(4, trend='c')650 gfevd = compute_gfevd(res, H=10)651 n = gfevd.shape[0]652 off_diag = gfevd.copy()653 np.fill_diagonal(off_diag, 0.0)654 TO = off_diag.sum(axis=0) * 100655 FROM = off_diag.sum(axis=1) * 100656 NET = TO - FROM657658 idx_btc = common_assets.index('BTC') if 'BTC' in common_assets else None659 idx_es = common_assets.index('ES') if 'ES' in common_assets else None660661 if idx_btc is not None and idx_es is not None:662 rolling_net_dates.append(end_date)663 rolling_net_btc.append(NET[idx_btc])664 rolling_net_es.append(NET[idx_es])665 except Exception:666 pass667668 if i % 200 == 0 and rolling_net_btc:669 print(f" Step {i}/{n_steps2}: BTC NET={rolling_net_btc[-1]:.1f}%, ES NET={rolling_net_es[-1]:.1f}%")670671net_df = pd.DataFrame({672 'date': rolling_net_dates,673 'NET_BTC': rolling_net_btc,674 'NET_ES': rolling_net_es,675}).set_index('date')676677# ─── Fig 6: Rolling NET BTC vs ES ─────────────────────────────────────────────678print("\nPlotting Fig 6: Rolling NET directional spillovers BTC vs ES...")679fig, ax = plt.subplots(figsize=(13, 5))680ax.plot(net_df.index, net_df['NET_BTC'], color=ASSET_COLORS['BTC'], linewidth=1.1, label='BTC NET spillover', alpha=0.9)681ax.plot(net_df.index, net_df['NET_ES'], color=ASSET_COLORS['ES'], linewidth=1.1, label='ES NET spillover', alpha=0.9)682ax.axhline(0, color='black', linewidth=0.8, linestyle='--')683ax.fill_between(net_df.index, net_df['NET_BTC'], 0,684 where=(net_df['NET_BTC'] > 0), alpha=0.15, color=ASSET_COLORS['BTC'], label='BTC net transmitter')685ax.fill_between(net_df.index, net_df['NET_BTC'], 0,686 where=(net_df['NET_BTC'] < 0), alpha=0.15, color='gray', label='BTC net receiver')687688for (label, date), c in zip(events.items(), colors_ev):689 if net_df.index.min() <= date <= net_df.index.max():690 ax.axvline(date, color=c, linestyle=':', linewidth=1, alpha=0.8)691 ax.text(date, ax.get_ylim()[1] if ax.get_ylim()[1] else net_df[['NET_BTC','NET_ES']].max().max(),692 label, ha='center', va='top', fontsize=7, color=c,693 bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.6))694695ax.set_xlabel('Date')696ax.set_ylabel('NET Directional Spillover (%)')697ax.set_title('Rolling 200-Day NET Directional Spillovers: BTC vs ES (VAR-4, H=10)', fontweight='bold')698ax.legend(loc='lower right', ncol=2)699ax.xaxis.set_major_locator(matplotlib.dates.YearLocator(1))700ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))701plt.setp(ax.xaxis.get_majorticklabels(), rotation=45, ha='right')702plt.tight_layout()703fig.savefig(os.path.join(FIGURES, 'fig6_net_directional_btc_es.png'), dpi=150, bbox_inches='tight')704plt.close()705print(" Saved fig6_net_directional_btc_es.png")706707# ─── 8. Network Visualization ─────────────────────────────────────────────────708print("\n" + "="*70)709print("Section 8: Network Graph...")710711if conn_table is not None:712 gfevd_mat = conn_table.loc[common_assets, common_assets].values.astype(float)713 np.fill_diagonal(gfevd_mat, 0.0)714715 # Threshold: keep edges above median716 all_edges = gfevd_mat.flatten()717 all_edges = all_edges[all_edges > 0]718 threshold = np.median(all_edges)719720 G = nx.DiGraph()721 for a in common_assets:722 G.add_node(a, asset_class=ASSET_CLASS.get(a, 'Other'))723724 TO_vals = {}725 for j, asj in enumerate(common_assets):726 TO_vals[asj] = sum(gfevd_mat[i, j] for i in range(len(common_assets)) if i != j)727728 for i, asi in enumerate(common_assets):729 for j, asj in enumerate(common_assets):730 if i != j and gfevd_mat[i, j] > threshold:731 G.add_edge(asj, asi, weight=gfevd_mat[i, j]) # j → i (j causes forecast error in i)732733 # Layout734 pos = nx.spring_layout(G, k=2.5, seed=42, iterations=100)735736 node_sizes = [max(300, TO_vals.get(a, 1) * 40) for a in G.nodes()]737 node_colors = [CLASS_COLOR.get(ASSET_CLASS.get(a, 'Other'), 'gray') for a in G.nodes()]738 edge_weights = [G[u][v]['weight'] for u, v in G.edges()]739 max_w = max(edge_weights) if edge_weights else 1740 edge_widths = [w / max_w * 4 for w in edge_weights]741 edge_alphas = [0.3 + 0.5 * w / max_w for w in edge_weights]742743 fig, ax = plt.subplots(figsize=(12, 9))744745 for (u, v), lw, alpha, w in zip(G.edges(), edge_widths, edge_alphas, edge_weights):746 nx.draw_networkx_edges(G, pos, edgelist=[(u, v)], ax=ax,747 width=lw, alpha=alpha, edge_color='#555555',748 arrows=True, arrowsize=12,749 connectionstyle='arc3,rad=0.15',750 node_size=node_sizes[list(G.nodes()).index(u)] if u in G.nodes() else 300)751752 nx.draw_networkx_nodes(G, pos, ax=ax,753 node_size=node_sizes,754 node_color=node_colors,755 alpha=0.9, edgecolors='white', linewidths=1.5)756 nx.draw_networkx_labels(G, pos, ax=ax, font_size=9, font_weight='bold', font_color='black')757758 # Legend for asset classes759 legend_patches = [mpatches.Patch(color=c, label=cls)760 for cls, c in CLASS_COLOR.items()]761 ax.legend(handles=legend_patches, title='Asset Class', loc='upper left', fontsize=9)762 ax.set_title(f'Volatility Transmission Network\n(Edges > median threshold = {threshold:.1f}%,\n'763 f'Node size ∝ Total TO spillover, Arrow direction = shock transmission)',764 fontweight='bold', fontsize=10)765 ax.axis('off')766 plt.tight_layout()767 fig.savefig(os.path.join(FIGURES, 'fig7_network_graph.png'), dpi=150, bbox_inches='tight')768 plt.close()769 print(" Saved fig7_network_graph.png")770771# ─── Summary Report ───────────────────────────────────────────────────────────772print("\n" + "=" * 70)773print("ANALYSIS COMPLETE — KEY NUMERICAL FINDINGS")774print("=" * 70)775776print("\n[Table 1] Descriptive Statistics (5-min log returns):")777for row in desc_rows:778 print(f" {row['Asset']:10s}: N={row['N_obs']:>8,}, Std={row['Std']:.6f}, "779 f"Skew={row['Skewness']:+.3f}, Kurt={row['Kurtosis']:+.2f}")780781print(f"\n[Section 2] Daily RV — Annualized means:")782for sym, rv in rv_dict.items():783 print(f" {sym:10s}: mean RV = {rv.mean():.4f}, max RV = {rv.max():.4f}")784785print(f"\n[Section 3] Full-Sample Connectedness:")786print(f" Total Spillover Index = {tsi_full:.2f}%" if tsi_full else " N/A")787if conn_table is not None:788 print(" NET spillovers (positive = net transmitter):")789 for a in common_assets:790 net_val = conn_table.loc['NET', a]791 to_val = conn_table.loc['TO', a]792 frm_val = conn_table.loc[a, 'FROM']793 print(f" {a:10s}: TO={to_val:+.2f}%, FROM={frm_val:+.2f}%, NET={net_val:+.2f}%")794795print(f"\n[Section 4] Rolling Connectedness ({len(rolling_df)} windows):")796print(f" Mean TSI: {rolling_df['total_spillover'].mean():.2f}%")797print(f" Max TSI: {rolling_df['total_spillover'].max():.2f}% on {rolling_df['total_spillover'].idxmax().date()}")798print(f" Min TSI: {rolling_df['total_spillover'].min():.2f}% on {rolling_df['total_spillover'].idxmin().date()}")799800print(f"\n[Section 5] Regime TSI:")801for r, tsi in regime_tsi.items():802 print(f" {r:20s}: TSI = {tsi:.2f}%")803804print(f"\n[Section 6] Granger causality (sample):")805sig_gc = granger_df[granger_df['Significant'] == 'Yes'] if len(granger_df) > 0 else pd.DataFrame()806print(f" {len(sig_gc)} / {len(granger_df)} pairs significant at 5%")807if len(sig_gc) > 0:808 print(sig_gc[['Session', 'Cause', 'Effect', 'F_stat', 'p_value']].head(10).to_string(index=False))809810print(f"\n[Section 7] Rolling NET BTC (last window): {rolling_net_btc[-1]:.2f}%" if rolling_net_btc else "")811print(f" Rolling NET ES (last window): {rolling_net_es[-1]:.2f}%" if rolling_net_es else "")812813print("\nAll outputs saved to:", RESULTS)814print("=" * 70)815