"""
WP 6 -- Revision of the connectedness analysis (editorial-board review, 2026-10-17).

Starts from the daily realized-volatility panel produced by 01_analysis.py
(results/daily_rv_all.csv, annualized realized volatility per instrument and calendar day)
and re-estimates every VAR/GFEVD object on a correctly aligned panel:

  * trading-day panel: a day enters only if all ten instruments have an own realized
    volatility (no forward filling; the original script forward-filled weekends and
    holidays, so that 1,137 of its 3,651 rows were stale copies);
  * log realized volatility (Andersen et al., 2003) instead of levels;
  * VAR(4) -- HQ criterion -- with H = 10, as in the original paper;
  * residual-bootstrap confidence intervals for TSI, TO, FROM, NET (Kilian, 1998);
  * robustness grid (transformation, p, H, sample, weekend-inclusive crypto RV, subsets);
  * rolling windows of 200 / 250 / 500 trading days; calendar-year sub-samples;
  * regime sub-samples (pre-COVID, COVID crash, 2022 hiking cycle, spot-ETF era);
  * Granger-causality tests with conventional F and HAC-robust Wald statistics;
  * replication of the original (flawed) construction for the audit trail.

Outputs: results/rev_*.csv, paper/tables/*.tex (generated, not typed), figures/fig*.png.
Run time: about 3 minutes on a laptop (bootstrap included).
"""
import os
import time
import warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import networkx as nx
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import adfuller, kpss, acf
import statsmodels.api as sm

warnings.filterwarnings('ignore')
T0 = time.time()

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
RESULTS = os.path.join(ROOT, 'results')
FIGURES = os.path.join(ROOT, 'figures')
TABLES = os.path.join(ROOT, 'paper', 'tables')
PAPER_FIG = os.path.join(ROOT, 'paper', 'figures')
for d in (RESULTS, FIGURES, TABLES, PAPER_FIG):
    os.makedirs(d, exist_ok=True)

ASSETS = ['ES', 'NQ', 'EURUSD', 'USDJPY', 'GBPUSD', 'GC', 'CL', 'BTC', 'ETH', 'VIX']
ASSET_CLASS = {'ES': 'Equity', 'NQ': 'Equity', 'EURUSD': 'FX', 'USDJPY': 'FX', 'GBPUSD': 'FX',
               'GC': 'Commodity', 'CL': 'Commodity', 'BTC': 'Crypto', 'ETH': 'Crypto', 'VIX': 'Volatility'}
CLASS_COLOR = {'Equity': '#1f77b4', 'FX': '#2ca02c', 'Commodity': '#ff7f0e', 'Crypto': '#9467bd',
               'Volatility': '#8c564b'}
P_BASE, H_BASE, W_BASE = 4, 10, 200
B_BOOT = 500
RNG = np.random.default_rng(20261017)
SAMPLE_START, SAMPLE_END = '2016-03-11', '2026-03-11'

plt.rcParams.update({'font.size': 10, 'axes.titlesize': 11, 'axes.labelsize': 10, 'figure.dpi': 150,
                     'savefig.dpi': 150, 'font.family': 'DejaVu Sans',
                     'axes.spines.top': False, 'axes.spines.right': False})


def log(msg):
    print(f"[{time.time() - T0:6.1f}s] {msg}", flush=True)


def savefig(fig, name):
    for d in (FIGURES, PAPER_FIG):
        fig.savefig(os.path.join(d, name), dpi=150, bbox_inches='tight')
    plt.close(fig)


# ─── 1. Panels ────────────────────────────────────────────────────────────────
rv_all = pd.read_csv(os.path.join(RESULTS, 'daily_rv_all.csv'), parse_dates=['date']).set_index('date')
rv_all = rv_all[ASSETS]

# (a) trading-day panel: all ten instruments observed (VIX trades only on CBOE days)
panel = rv_all.loc[SAMPLE_START:SAMPLE_END].dropna()
log(f"trading-day panel: {panel.index[0].date()} -> {panel.index[-1].date()}, T = {len(panel)}; "
    f"weekday counts = {panel.index.dayofweek.value_counts().sort_index().to_dict()}")

# (b) replication of the original construction (calendar days, forward fill up to 5 days, levels)
orig_panel = rv_all.ffill(limit=5).dropna(thresh=len(ASSETS) - 2).dropna()
n_weekend = int((orig_panel.index.dayofweek >= 5).sum())
n_ffilled_rows = int((orig_panel.index.isin(rv_all.dropna().index) == False).sum())
log(f"original construction: T = {len(orig_panel)}, weekend rows = {n_weekend}, "
    f"rows with at least one forward-filled value = {n_ffilled_rows}")

# (c) weekend-inclusive crypto RV: average daily variance over the calendar days since the
#     previous trading day (so Monday carries Saturday, Sunday and Monday variance)
def weekend_inclusive(series, trading_days):
    var = series.dropna() ** 2
    out = {}
    prev = None
    for d in trading_days:
        if prev is None:
            out[d] = var.get(d, np.nan)
        else:
            win = var.loc[prev + pd.Timedelta(days=1):d]
            out[d] = win.mean() if len(win) else np.nan
        prev = d
    return pd.Series(out) ** 0.5


panel_wk = panel.copy()
for c in ('BTC', 'ETH'):
    panel_wk[c] = weekend_inclusive(rv_all[c], panel.index)
panel_wk = panel_wk.dropna()


# ─── 2. Descriptive statistics of daily RV ────────────────────────────────────
rows = []
for a in ASSETS:
    s = panel[a]
    ls = np.log(s)
    adf_l = adfuller(ls, regression='c', autolag='AIC')
    kp_l = kpss(ls, regression='c', nlags='auto')
    rows.append({'Asset': a, 'Class': ASSET_CLASS[a], 'Mean_RV_pct': 100 * s.mean(), 'Median_RV_pct': 100 * s.median(),
                 'Max_RV_pct': 100 * s.max(), 'Skew_level': s.skew(), 'Kurt_level': s.kurt(),
                 'Skew_log': ls.skew(), 'Kurt_log': ls.kurt(), 'AC1_log': acf(ls, nlags=22)[1],
                 'AC22_log': acf(ls, nlags=22)[22], 'ADF_log': adf_l[0], 'ADF_p': adf_l[1],
                 'KPSS_log': kp_l[0], 'KPSS_p': kp_l[1]})
desc = pd.DataFrame(rows).set_index('Asset')
desc.to_csv(os.path.join(RESULTS, 'rev_table_rv_descriptives.csv'), float_format='%.4f')
log("descriptives done")


# ─── 3. GFEVD machinery ───────────────────────────────────────────────────────
def gfevd_from_var(coefs, sigma, H):
    """Normalized generalized FEVD (Koop et al., 1996; Pesaran and Shin, 1998)."""
    n = sigma.shape[0]
    p = coefs.shape[0]
    psi = [np.eye(n)]
    for h in range(1, H):
        acc = np.zeros((n, n))
        for j in range(1, min(h, p) + 1):
            acc += coefs[j - 1] @ psi[h - j]
        psi.append(acc)
    num = np.zeros((n, n))
    den = np.zeros(n)
    for h in range(H):
        ps = psi[h] @ sigma
        num += ps ** 2
        den += np.einsum('ij,ij->i', ps, psi[h])
    theta = num / (np.diag(sigma)[None, :] * den[:, None])
    return theta / theta.sum(axis=1, keepdims=True)


def connectedness(y, p=P_BASE, H=H_BASE):
    """Fit VAR(p) with constant, return dict of connectedness objects (in percent)."""
    res = VAR(y.values).fit(p, trend='c')
    th = gfevd_from_var(res.coefs, res.sigma_u, H) * 100
    off = th.copy()
    np.fill_diagonal(off, 0)
    out = {'theta': th, 'FROM': off.sum(axis=1), 'TO': off.sum(axis=0), 'TSI': off.sum() / th.shape[0],
           'res': res}
    out['NET'] = out['TO'] - out['FROM']
    out['NETPAIR'] = th.T - th  # [i,j] = theta_ji - theta_ij : net from i to j
    return out


def bootstrap_connectedness(y, p=P_BASE, H=H_BASE, B=B_BOOT, rng=RNG):
    """Residual (i.i.d.) bootstrap of the VAR; returns arrays of TSI, TO, FROM, NET."""
    base = connectedness(y, p, H)
    res = base['res']
    T, n = y.shape
    u = res.resid - res.resid.mean(axis=0)
    A = res.coefs
    c = res.intercept
    Y = y.values
    tsi_b, to_b, from_b, net_b = [], [], [], []
    for b in range(B):
        idx = rng.integers(0, len(u), len(u))
        ub = u[idx]
        yb = np.empty_like(Y)
        yb[:p] = Y[:p]
        for t in range(p, T):
            acc = c.copy()
            for l in range(p):
                acc += A[l] @ yb[t - l - 1]
            yb[t] = acc + ub[t - p]
        rb = VAR(yb).fit(p, trend='c')
        th = gfevd_from_var(rb.coefs, rb.sigma_u, H) * 100
        off = th.copy()
        np.fill_diagonal(off, 0)
        fr, to = off.sum(axis=1), off.sum(axis=0)
        tsi_b.append(off.sum() / n)
        to_b.append(to)
        from_b.append(fr)
        net_b.append(to - fr)
    return base, np.array(tsi_b), np.array(to_b), np.array(from_b), np.array(net_b)


def ci(arr, q=(5, 95)):
    return np.percentile(arr, q, axis=0)


# ─── 4. Full-sample connectedness (log RV, trading days) with bootstrap ───────
ylog = np.log(panel)
log("full-sample bootstrap ...")
base, tsi_b, to_b, from_b, net_b = bootstrap_connectedness(ylog)
tsi_ci, to_ci, from_ci, net_ci = ci(tsi_b), ci(to_b), ci(from_b), ci(net_b)
log(f"full-sample TSI = {base['TSI']:.2f} [{tsi_ci[0]:.2f}, {tsi_ci[1]:.2f}]")

full_tbl = pd.DataFrame(base['theta'], index=ASSETS, columns=ASSETS)
full_tbl['FROM'] = base['FROM']
full_tbl.loc['TO'] = list(base['TO']) + [base['TSI']]
full_tbl.loc['NET'] = list(base['NET']) + [np.nan]
full_tbl.to_csv(os.path.join(RESULTS, 'rev_table_connectedness_full.csv'), float_format='%.2f')
ci_tbl = pd.DataFrame({'TO': base['TO'], 'TO_lo': to_ci[0], 'TO_hi': to_ci[1], 'FROM': base['FROM'],
                       'FROM_lo': from_ci[0], 'FROM_hi': from_ci[1], 'NET': base['NET'], 'NET_lo': net_ci[0],
                       'NET_hi': net_ci[1]}, index=ASSETS)
ci_tbl.loc['TSI', ['NET', 'NET_lo', 'NET_hi']] = [base['TSI'], tsi_ci[0], tsi_ci[1]]
ci_tbl.to_csv(os.path.join(RESULTS, 'rev_table_connectedness_full_ci.csv'), float_format='%.2f')

# original construction replicated (audit) and intermediate variants
orig = connectedness(orig_panel, P_BASE, H_BASE)
lvl_td = connectedness(panel, P_BASE, H_BASE)
audit = pd.DataFrame([
    {'Construction': 'Original: calendar days, forward-filled, RV levels', 'T': len(orig_panel), 'TSI': orig['TSI'],
     'NET_ES': orig['NET'][0], 'NET_NQ': orig['NET'][1], 'NET_VIX': orig['NET'][9], 'NET_BTC': orig['NET'][7],
     'NET_USDJPY': orig['NET'][3]},
    {'Construction': 'Trading days, RV levels', 'T': len(panel), 'TSI': lvl_td['TSI'], 'NET_ES': lvl_td['NET'][0],
     'NET_NQ': lvl_td['NET'][1], 'NET_VIX': lvl_td['NET'][9], 'NET_BTC': lvl_td['NET'][7],
     'NET_USDJPY': lvl_td['NET'][3]},
    {'Construction': 'Trading days, log RV (revised baseline)', 'T': len(panel), 'TSI': base['TSI'],
     'NET_ES': base['NET'][0], 'NET_NQ': base['NET'][1], 'NET_VIX': base['NET'][9], 'NET_BTC': base['NET'][7],
     'NET_USDJPY': base['NET'][3]},
])
audit.to_csv(os.path.join(RESULTS, 'rev_audit_construction.csv'), index=False, float_format='%.2f')
log("audit table done")

# ─── 5. Robustness grid ───────────────────────────────────────────────────────
rob_rows = []


def rob(label, y, p=P_BASE, H=H_BASE):
    r = connectedness(y, p, H)
    d = {'Specification': label, 'T': len(y), 'p': p, 'H': H, 'TSI': r['TSI']}
    for a in ['ES', 'NQ', 'USDJPY', 'GC', 'CL', 'BTC', 'ETH', 'VIX']:
        if a in y.columns:
            d[f'NET_{a}'] = r['NET'][list(y.columns).index(a)]
        else:
            d[f'NET_{a}'] = np.nan
    if 'BTC' in y.columns:
        d['Own_BTC'] = r['theta'][list(y.columns).index('BTC'), list(y.columns).index('BTC')]
    rob_rows.append(d)
    return r


rob('Baseline: log RV, trading days, VAR(4), H = 10', ylog)
rob('RV in levels (trading days)', panel)
for p in (1, 2, 3, 6):
    rob(f'VAR({p})', ylog, p=p)
for H in (5, 20):
    rob(f'H = {H}', ylog, H=H)
rob('Weekend-inclusive crypto RV', np.log(panel_wk))
rob('Sample 2018-01 to 2026-03', ylog.loc['2018-01-01':])
rob('Sample 2016-03 to 2019-12', ylog.loc[:'2019-12-31'])
rob('Sample 2020-01 to 2026-03', ylog.loc['2020-01-01':])
rob('Excluding VIX', ylog.drop(columns=['VIX']))
rob('Excluding BTC and ETH', ylog.drop(columns=['BTC', 'ETH']))
rob('Excluding ETH (BTC only)', ylog.drop(columns=['ETH']))
rob_df = pd.DataFrame(rob_rows)
rob_df.to_csv(os.path.join(RESULTS, 'rev_table_robustness.csv'), index=False, float_format='%.2f')
log("robustness grid done")

# ─── 6. Rolling windows ───────────────────────────────────────────────────────
def rolling(y, w, p=P_BASE, H=H_BASE):
    dates, tsi, net = [], [], []
    Y = y.values
    for i in range(len(y) - w + 1):
        r = VAR(Y[i:i + w]).fit(p, trend='c')
        th = gfevd_from_var(r.coefs, r.sigma_u, H) * 100
        off = th.copy()
        np.fill_diagonal(off, 0)
        dates.append(y.index[i + w - 1])
        tsi.append(off.sum() / th.shape[0])
        net.append(off.sum(axis=0) - off.sum(axis=1))
    out = pd.DataFrame(net, index=dates, columns=[f'NET_{a}' for a in y.columns])
    out.insert(0, 'TSI', tsi)
    out.index.name = 'date'
    return out


roll = {}
for w in (200, 250, 500):
    roll[w] = rolling(ylog, w)
    roll[w].to_csv(os.path.join(RESULTS, f'rev_rolling_w{w}.csv'), float_format='%.3f')
    log(f"rolling w={w}: {len(roll[w])} windows, TSI mean {roll[w]['TSI'].mean():.2f}, "
        f"min {roll[w]['TSI'].min():.2f} ({roll[w]['TSI'].idxmin().date()}), "
        f"max {roll[w]['TSI'].max():.2f} ({roll[w]['TSI'].idxmax().date()})")
roll_orig = pd.read_csv(os.path.join(RESULTS, 'rolling_total_spillover.csv'), parse_dates=['date']).set_index('date')
roll_summary = pd.DataFrame({
    'window': ['original (200 calendar days, levels, forward-filled)', '200 trading days', '250 trading days',
               '500 trading days'],
    'n_windows': [len(roll_orig), len(roll[200]), len(roll[250]), len(roll[500])],
    'mean': [roll_orig['total_spillover'].mean()] + [roll[w]['TSI'].mean() for w in (200, 250, 500)],
    'min': [roll_orig['total_spillover'].min()] + [roll[w]['TSI'].min() for w in (200, 250, 500)],
    'min_date': [roll_orig['total_spillover'].idxmin().date()] + [roll[w]['TSI'].idxmin().date() for w in (200, 250, 500)],
    'max': [roll_orig['total_spillover'].max()] + [roll[w]['TSI'].max() for w in (200, 250, 500)],
    'max_date': [roll_orig['total_spillover'].idxmax().date()] + [roll[w]['TSI'].idxmax().date() for w in (200, 250, 500)],
})
roll_summary.to_csv(os.path.join(RESULTS, 'rev_rolling_summary.csv'), index=False, float_format='%.2f')

# ─── 7. Regimes (with bootstrap CIs on NET) ───────────────────────────────────
REGIMES = {
    'Pre-COVID (2019-01 to 2020-02)': ('2019-01-01', '2020-02-28'),
    'COVID crash (2020-03 to 2020-06)': ('2020-03-01', '2020-06-30'),
    'Rate hikes (2022-01 to 2022-12)': ('2022-01-01', '2022-12-31'),
    'Spot-ETF era (2024-01 to 2026-03)': ('2024-01-11', '2026-03-11'),
}
reg_rows, reg_full = [], {}
for name, (s, e) in REGIMES.items():
    y = ylog.loc[s:e]
    p = P_BASE if len(y) >= 150 else 1
    b, tb, tob, frb, nb = bootstrap_connectedness(y, p=p, B=B_BOOT)
    tci, nci = ci(tb), ci(nb)
    row = {'Regime': name, 'T': len(y), 'p': p, 'TSI': b['TSI'], 'TSI_lo': tci[0], 'TSI_hi': tci[1]}
    for k, a in enumerate(ASSETS):
        row[f'NET_{a}'] = b['NET'][k]
        row[f'NET_{a}_lo'] = nci[0][k]
        row[f'NET_{a}_hi'] = nci[1][k]
    reg_rows.append(row)
    tbl = pd.DataFrame(b['theta'], index=ASSETS, columns=ASSETS)
    tbl['FROM'] = b['FROM']
    tbl.loc['TO'] = list(b['TO']) + [b['TSI']]
    tbl.loc['NET'] = list(b['NET']) + [np.nan]
    reg_full[name] = tbl
    slug = name.split(' (')[0].lower().replace(' ', '_').replace('-', '')
    tbl.to_csv(os.path.join(RESULTS, f'rev_table_regime_{slug}.csv'), float_format='%.2f')
    log(f"regime {name}: T={len(y)}, p={p}, TSI={b['TSI']:.1f} [{tci[0]:.1f},{tci[1]:.1f}], "
        f"NET BTC={b['NET'][7]:+.1f} [{nci[0][7]:+.1f},{nci[1][7]:+.1f}], NET USDJPY={b['NET'][3]:+.1f}")
    # same regime with p = 4 for comparison when p = 1 was used
    if p != P_BASE:
        alt = connectedness(y, P_BASE, H_BASE)
        log(f"   (VAR(4) on the same window: TSI={alt['TSI']:.1f}, NET USDJPY={alt['NET'][3]:+.1f}, "
            f"NET BTC={alt['NET'][7]:+.1f})")
        reg_rows[-1]['TSI_var4'] = alt['TSI']
        reg_rows[-1]['NET_USDJPY_var4'] = alt['NET'][3]
        reg_rows[-1]['NET_BTC_var4'] = alt['NET'][7]
reg_df = pd.DataFrame(reg_rows)
reg_df.to_csv(os.path.join(RESULTS, 'rev_table_regimes.csv'), index=False, float_format='%.2f')

# ─── 8. Calendar-year sub-samples (crypto integration) ────────────────────────
yr_rows = []
for yr in range(2016, 2027):
    y = ylog.loc[str(yr)]
    if len(y) < 60:
        continue
    b, tb, tob, frb, nb = bootstrap_connectedness(y, p=1 if len(y) < 150 else 2, B=300)
    tci, nci = ci(tb), ci(nb)
    i_btc, i_es, i_eth = ASSETS.index('BTC'), ASSETS.index('ES'), ASSETS.index('ETH')
    th = b['theta']
    yr_rows.append({'Year': yr, 'T': len(y), 'TSI': b['TSI'], 'TSI_lo': tci[0], 'TSI_hi': tci[1],
                    'NET_BTC': b['NET'][i_btc], 'NET_BTC_lo': nci[0][i_btc], 'NET_BTC_hi': nci[1][i_btc],
                    'NET_ETH': b['NET'][i_eth], 'NET_ETH_lo': nci[0][i_eth], 'NET_ETH_hi': nci[1][i_eth],
                    'FROM_BTC': b['FROM'][i_btc], 'FROM_BTC_lo': ci(frb)[0][i_btc], 'FROM_BTC_hi': ci(frb)[1][i_btc],
                    'FROM_BTC_excl_ETH': b['FROM'][i_btc] - th[i_btc, i_eth],
                    'BTC_from_ES': th[i_btc, i_es], 'ES_from_BTC': th[i_es, i_btc],
                    'NET_ES': b['NET'][i_es], 'NET_ES_lo': nci[0][i_es], 'NET_ES_hi': nci[1][i_es]})
yr_df = pd.DataFrame(yr_rows)
yr_df.to_csv(os.path.join(RESULTS, 'rev_table_yearly.csv'), index=False, float_format='%.2f')
log("yearly sub-samples done")

# ─── 9. Granger causality (log RV, p lags, F and HAC Wald) ────────────────────
PAIRS = [('ES', 'VIX'), ('NQ', 'VIX'), ('VIX', 'ES'), ('VIX', 'NQ'), ('ES', 'NQ'), ('NQ', 'ES'),
         ('ES', 'BTC'), ('NQ', 'BTC'), ('VIX', 'BTC'), ('BTC', 'ES'), ('BTC', 'NQ'), ('BTC', 'VIX'),
         ('USDJPY', 'ES'), ('ES', 'USDJPY'), ('GC', 'BTC'), ('BTC', 'GC')]


def granger(y, cause, effect, p=P_BASE, hac_lags=10):
    df = pd.DataFrame({'y': y[effect], 'x': y[cause]})
    X = pd.concat([df['y'].shift(l).rename(f'y_l{l}') for l in range(1, p + 1)] +
                  [df['x'].shift(l).rename(f'x_l{l}') for l in range(1, p + 1)], axis=1)
    X = sm.add_constant(X)
    data = pd.concat([df['y'], X], axis=1).dropna()
    yy, XX = data['y'], data.drop(columns='y')
    m_u = sm.OLS(yy, XX).fit()
    m_r = sm.OLS(yy, XX[[c for c in XX.columns if not c.startswith('x_')]]).fit()
    T = len(yy)
    F = ((m_r.ssr - m_u.ssr) / p) / (m_u.ssr / (T - 2 * p - 1))
    from scipy import stats
    pF = 1 - stats.f.cdf(F, p, T - 2 * p - 1)
    m_hac = sm.OLS(yy, XX).fit(cov_type='HAC', cov_kwds={'maxlags': hac_lags})
    R = np.zeros((p, XX.shape[1]))
    for k, c in enumerate([c for c in XX.columns if c.startswith('x_')]):
        R[k, list(XX.columns).index(c)] = 1
    w = m_hac.wald_test(R, use_f=False, scalar=True)
    return {'Cause': cause, 'Effect': effect, 'T': T, 'F': F, 'p_F': pF, 'Wald_HAC': float(w.statistic),
            'p_HAC': float(w.pvalue), 'sum_coef': float(m_u.params[[c for c in XX.columns if c.startswith('x_')]].sum())}


gr = pd.DataFrame([granger(ylog, c, e) for c, e in PAIRS])
gr.to_csv(os.path.join(RESULTS, 'rev_table_granger.csv'), index=False, float_format='%.4f')
log("Granger tests done")

# ─── 10. Lag-order criteria ───────────────────────────────────────────────────
sel = VAR(ylog.values).select_order(maxlags=10)
lag_tbl = pd.DataFrame({'p': range(0, 11), 'AIC': sel.ics['aic'], 'BIC': sel.ics['bic'], 'HQ': sel.ics['hqic']})
lag_tbl.to_csv(os.path.join(RESULTS, 'rev_table_lag_selection.csv'), index=False, float_format='%.3f')
log(f"lag selection: {sel.selected_orders}")


# ─── 11. LaTeX tables (generated) ─────────────────────────────────────────────
def fmt(x, d=1, sign=False):
    if pd.isna(x):
        return ''
    s = f"{x:+.{d}f}" if sign else f"{x:.{d}f}"
    return s.replace('-', '$-$').replace('+', '$+$') if sign else s


def write(name, text):
    with open(os.path.join(TABLES, name), 'w') as f:
        f.write(text)


# 11a. RV descriptives
lines = [r"\begin{tabular}{@{}llrrrrrrrrr@{}}", r"\toprule",
         r"Asset & Class & Mean & Median & Max & Skew & Skew & Kurt & $\rho_1$ & ADF & KPSS \\",
         r" & & \multicolumn{3}{c}{annualized RV (\%)} & level & log & log & log & log & log \\",
         r"\cmidrule(lr){3-5}\cmidrule(lr){6-6}\cmidrule(lr){7-9}\cmidrule(lr){10-11}"]
lines[-1] = r"\midrule"
for a in ASSETS:
    r = desc.loc[a]
    kp = fmt(r['KPSS_log'], 2) + ('' if r['KPSS_p'] > 0.05 else r'$^{*}$')
    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 ''
    lines.append(f"{a} & {r['Class']} & {fmt(r['Mean_RV_pct'])} & {fmt(r['Median_RV_pct'])} & {fmt(r['Max_RV_pct'], 0)} & "
                 f"{fmt(r['Skew_level'], 2)} & {fmt(r['Skew_log'], 2)} & {fmt(r['Kurt_log'], 2)} & {fmt(r['AC1_log'], 2)} & "
                 f"{fmt(r['ADF_log'], 1)}{adf_st} & {kp} \\\\")
lines += [r"\bottomrule", r"\end{tabular}"]
write('tab_rv_descriptives.tex', "\n".join(lines))

# 11b. full connectedness table (with CI rows)
def conn_table_tex(tbl, to_ci_=None, from_ci_=None, net_ci_=None, tsi_ci_=None):
    cols = ASSETS
    L = [r"\begin{tabular}{@{}l*{10}{r}r@{}}", r"\toprule",
         r"& \multicolumn{10}{c}{Shock origin $j$} & \\", r"\cmidrule(lr){2-11}",
         "Recipient $i$ & " + " & ".join(cols) + r" & FROM \\", r"\midrule"]
    for i, a in enumerate(cols):
        vals = [fmt(tbl.loc[a, c], 1) for c in cols]
        vals = [f"\\textbf{{{v}}}" if c == a else v for v, c in zip(vals, cols)]
        L.append(f"{a} & " + " & ".join(vals) + f" & {fmt(tbl.loc[a, 'FROM'], 1)} \\\\")
        if from_ci_ is not None:
            L.append(" & " + " & ".join([''] * 10) + f" & {{\\scriptsize[{fmt(from_ci_[0][i], 1)}, {fmt(from_ci_[1][i], 1)}]}} \\\\")
    L.append(r"\midrule")
    L.append("TO & " + " & ".join(fmt(tbl.loc['TO', c], 1) for c in cols) + f" & \\textbf{{{fmt(tbl.loc['TO', 'FROM'], 1)}}} \\\\")
    if to_ci_ is not None:
        L.append(" & " + " & ".join(f"{{\\scriptsize[{fmt(to_ci_[0][k], 1)}, {fmt(to_ci_[1][k], 1)}]}}" for k in range(10)) +
                 f" & {{\\scriptsize[{fmt(tsi_ci_[0], 1)}, {fmt(tsi_ci_[1], 1)}]}} \\\\")
    L.append("NET & " + " & ".join(fmt(tbl.loc['NET', c], 1, sign=True) for c in cols) + r" & \\")
    if net_ci_ is not None:
        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" & \\")
    L += [r"\bottomrule", r"\end{tabular}"]
    return "\n".join(L)


write('tab_connectedness_full.tex', conn_table_tex(full_tbl, to_ci, from_ci, net_ci, tsi_ci))
for name, tbl in reg_full.items():
    slug = name.split(' (')[0].lower().replace(' ', '_').replace('-', '')
    write(f'tab_regime_{slug}.tex', conn_table_tex(tbl))

# 11c. audit / construction table
L = [r"\begin{tabular}{@{}lrrrrrrr@{}}", r"\toprule",
     r"Construction of the panel & $T$ & TSI & ES & NQ & VIX & BTC & USDJPY \\", r"\midrule"]
for _, r in audit.iterrows():
    L.append(f"{r['Construction']} & {int(r['T'])} & {fmt(r['TSI'])} & {fmt(r['NET_ES'], 1, True)} & {fmt(r['NET_NQ'], 1, True)} & "
             f"{fmt(r['NET_VIX'], 1, True)} & {fmt(r['NET_BTC'], 1, True)} & {fmt(r['NET_USDJPY'], 1, True)} \\\\")
L += [r"\bottomrule", r"\end{tabular}"]
write('tab_audit_construction.tex', "\n".join(L))

# 11d. robustness grid
L = [r"\begin{tabular}{@{}lrrrrrrrrr@{}}", r"\toprule",
     r"Specification & $T$ & TSI & ES & NQ & USDJPY & CL & BTC & VIX & BTC own \\", r"\midrule"]
for _, r in rob_df.iterrows():
    L.append(f"{r['Specification']} & {int(r['T'])} & {fmt(r['TSI'])} & {fmt(r['NET_ES'], 1, True)} & {fmt(r['NET_NQ'], 1, True)} & "
             f"{fmt(r['NET_USDJPY'], 1, True)} & {fmt(r['NET_CL'], 1, True)} & {fmt(r['NET_BTC'], 1, True)} & "
             f"{fmt(r['NET_VIX'], 1, True)} & {fmt(r.get('Own_BTC', np.nan))} \\\\")
L += [r"\bottomrule", r"\end{tabular}"]
write('tab_robustness.tex', "\n".join(L))

# 11e. regimes summary with CIs
L = [r"\begin{tabular}{@{}l" + "r" * len(REGIMES) + "@{}}", r"\toprule",
     " & " + " & ".join(n.split(' (')[0] for n in REGIMES) + r" \\",
     " & " + " & ".join("{\\scriptsize " + n.split(' (')[1].rstrip(')') + "}" for n in REGIMES) + r" \\", r"\midrule",
     "$T$ (trading days) & " + " & ".join(str(int(r['T'])) for _, r in reg_df.iterrows()) + r" \\",
     "VAR lag $p$ & " + " & ".join(str(int(r['p'])) for _, r in reg_df.iterrows()) + r" \\",
     "TSI (\\%) & " + " & ".join(fmt(r['TSI']) for _, r in reg_df.iterrows()) + r" \\",
     " & " + " & ".join(f"{{\\scriptsize[{fmt(r['TSI_lo'])}, {fmt(r['TSI_hi'])}]}}" for _, r in reg_df.iterrows()) + r" \\",
     r"\midrule", r"\multicolumn{" + str(1 + len(REGIMES)) + r"}{@{}l}{\textit{Net directional spillovers (pp)}} \\"]
for a in ASSETS:
    L.append(f"{a} & " + " & ".join(fmt(r[f'NET_{a}'], 1, True) for _, r in reg_df.iterrows()) + r" \\")
    L.append(" & " + " & ".join(f"{{\\scriptsize[{fmt(r[f'NET_{a}_lo'], 1, True)}, {fmt(r[f'NET_{a}_hi'], 1, True)}]}}"
                                for _, r in reg_df.iterrows()) + r" \\")
L += [r"\bottomrule", r"\end{tabular}"]
write('tab_regimes.tex', "\n".join(L))

# 11f. yearly table
L = [r"\begin{tabular}{@{}lrrrrrrr@{}}", r"\toprule",
     r"Year & $T$ & TSI & NET BTC & NET ETH & FROM BTC & BTC $\leftarrow$ ES & ES $\leftarrow$ BTC \\",
     r" & & (\%) & (pp) & (pp) & excl.\ ETH (\%) & (\%) & (\%) \\", r"\midrule"]
for _, r in yr_df.iterrows():
    L.append(f"{int(r['Year'])} & {int(r['T'])} & {fmt(r['TSI'])} {{\\scriptsize[{fmt(r['TSI_lo'])}, {fmt(r['TSI_hi'])}]}} & "
             f"{fmt(r['NET_BTC'], 1, True)} {{\\scriptsize[{fmt(r['NET_BTC_lo'], 1, True)}, {fmt(r['NET_BTC_hi'], 1, True)}]}} & "
             f"{fmt(r['NET_ETH'], 1, True)} & {fmt(r['FROM_BTC_excl_ETH'])} & {fmt(r['BTC_from_ES'])} & {fmt(r['ES_from_BTC'])} \\\\")
L += [r"\bottomrule", r"\end{tabular}"]
write('tab_yearly.tex', "\n".join(L))

# 11g. Granger
def stars(p):
    return r'$^{***}$' if p < 0.01 else r'$^{**}$' if p < 0.05 else r'$^{*}$' if p < 0.10 else ''


L = [r"\begin{tabular}{@{}llrrrrr@{}}", r"\toprule",
     r"Cause & Effect & $T$ & $F$ & $p$ & Wald (HAC) & $p$ \\", r"\midrule"]
for _, r in gr.iterrows():
    L.append(f"{r['Cause']} & {r['Effect']} & {int(r['T'])} & {fmt(r['F'], 2)}{stars(r['p_F'])} & {fmt(r['p_F'], 3)} & "
             f"{fmt(r['Wald_HAC'], 2)}{stars(r['p_HAC'])} & {fmt(r['p_HAC'], 3)} \\\\")
L += [r"\bottomrule", r"\end{tabular}"]
write('tab_granger.tex', "\n".join(L))

# 11h. lag selection
L = [r"\begin{tabular}{@{}lrrrrrrrrrr@{}}", r"\toprule", "Criterion & " + " & ".join(f"$p={p}$" for p in range(1, 11)) + r" \\", r"\midrule"]
for crit in ('AIC', 'BIC', 'HQ'):
    vals = lag_tbl[crit].values[1:]
    best = np.argmin(vals)
    L.append(f"{crit} & " + " & ".join((f"\\textbf{{{v:.2f}}}" if k == best else f"{v:.2f}") for k, v in enumerate(vals)) + r" \\")
L += [r"\bottomrule", r"\end{tabular}"]
write('tab_lag_selection.tex', "\n".join(L))

# 11i. rolling summary
L = [r"\begin{tabular}{@{}lrrrrrr@{}}", r"\toprule", r"Window & Windows & Mean & Min & (date) & Max & (date) \\", r"\midrule"]
for _, r in roll_summary.iterrows():
    L.append(f"{r['window']} & {int(r['n_windows'])} & {fmt(r['mean'])} & {fmt(r['min'])} & {r['min_date']} & {fmt(r['max'])} & {r['max_date']} \\\\")
L += [r"\bottomrule", r"\end{tabular}"]
write('tab_rolling_summary.tex', "\n".join(L))
log("LaTeX tables written")

# ─── 12. Figures ──────────────────────────────────────────────────────────────
# Fig 2: heatmap
fig, ax = plt.subplots(figsize=(10, 7.5))
mat = base['theta']
im = ax.imshow(mat, cmap='Blues', vmin=0, vmax=mat.max())
for i in range(10):
    for j in range(10):
        ax.text(j, i, f"{mat[i, j]:.1f}", ha='center', va='center', fontsize=8,
                color='white' if mat[i, j] > 0.55 * mat.max() else 'black')
ax.set_xticks(range(10)); ax.set_yticks(range(10))
ax.set_xticklabels(ASSETS, rotation=45, ha='right'); ax.set_yticklabels(ASSETS)
ax.set_xlabel('Shock origin $j$'); ax.set_ylabel('Recipient $i$')
ax.set_title(f'Generalized FEVD shares (%), log RV, VAR(4), H = 10 — TSI = {base["TSI"]:.1f}%')
plt.colorbar(im, ax=ax, label='% of 10-day forecast-error variance')
plt.tight_layout(); savefig(fig, 'fig2_connectedness_heatmap.png')

# Fig 3: rolling TSI (three windows) + original for reference
fig, ax = plt.subplots(figsize=(13, 5))
ax.plot(roll[200].index, roll[200]['TSI'], color='#1f77b4', lw=1.0, label='200 trading days')
ax.plot(roll[250].index, roll[250]['TSI'], color='#ff7f0e', lw=0.9, alpha=0.9, label='250 trading days')
ax.plot(roll[500].index, roll[500]['TSI'], color='#2ca02c', lw=1.3, label='500 trading days')
ax.axhline(base['TSI'], color='gray', ls='--', lw=0.8, label=f'full sample ({base["TSI"]:.1f}%)')
events = {'COVID crash\n(Mar 2020)': '2020-03-16', 'First Fed hike\n(Mar 2022)': '2022-03-16',
          'FTX\n(Nov 2022)': '2022-11-11', 'Spot BTC ETFs\n(Jan 2024)': '2024-01-11', 'Volatility spike\n(Aug 2024)': '2024-08-05'}
for k, (lab, d) in enumerate(events.items()):
    d = pd.Timestamp(d)
    ax.axvline(d, color='k', ls=':', lw=0.8, alpha=0.6)
    ax.text(d, ax.get_ylim()[1] - (k % 2) * 3, lab, ha='center', va='top', fontsize=7)
ax.set_ylabel('Total spillover index (%)'); ax.set_xlabel('End of window')
ax.set_title('Rolling total spillover index — log RV, VAR(4), H = 10, trading-day panel')
ax.legend(loc='lower left', ncol=4, fontsize=8)
plt.tight_layout(); savefig(fig, 'fig3_rolling_spillover.png')

# Fig 4: NET by regime with CI bars
fig, ax = plt.subplots(figsize=(12, 5.5))
x = np.arange(10); nreg = len(REGIMES); width = 0.8 / nreg
pal = ['#1f77b4', '#d62728', '#ff7f0e', '#9467bd']
for k, (_, r) in enumerate(reg_df.iterrows()):
    vals = np.array([r[f'NET_{a}'] for a in ASSETS])
    lo = np.array([r[f'NET_{a}_lo'] for a in ASSETS]); hi = np.array([r[f'NET_{a}_hi'] for a in ASSETS])
    off = (k - (nreg - 1) / 2) * width
    ax.bar(x + off, vals, width * 0.9, color=pal[k], alpha=0.85, label=f"{r['Regime']} (TSI {r['TSI']:.1f}%)")
    ax.errorbar(x + off, vals, yerr=[vals - lo, hi - vals], fmt='none', ecolor='k', elinewidth=0.6, capsize=1.5)
ax.axhline(0, color='k', lw=0.8); ax.set_xticks(x); ax.set_xticklabels(ASSETS)
ax.set_ylabel('Net directional spillover (pp)'); ax.legend(fontsize=8)
ax.set_title('Net directional spillovers by regime, with 90% bootstrap intervals')
plt.tight_layout(); savefig(fig, 'fig4_net_spillovers_regimes.png')

# Fig 6: rolling NET BTC vs ES (w = 200 and 500)
fig, ax = plt.subplots(figsize=(13, 5))
ax.plot(roll[200].index, roll[200]['NET_BTC'], color='#9467bd', lw=0.8, alpha=0.6, label='BTC, 200 days')
ax.plot(roll[500].index, roll[500]['NET_BTC'], color='#9467bd', lw=1.6, label='BTC, 500 days')
ax.plot(roll[200].index, roll[200]['NET_ES'], color='#1f77b4', lw=0.8, alpha=0.6, label='ES, 200 days')
ax.plot(roll[500].index, roll[500]['NET_ES'], color='#1f77b4', lw=1.6, label='ES, 500 days')
ax.axhline(0, color='k', lw=0.8, ls='--')
for lab, d in events.items():
    ax.axvline(pd.Timestamp(d), color='k', ls=':', lw=0.8, alpha=0.6)
ax.set_ylabel('Net directional spillover (pp)'); ax.set_xlabel('End of window')
ax.set_title('Rolling net directional spillovers: Bitcoin and E-mini S&P 500')
ax.legend(ncol=4, fontsize=8, loc='lower left')
plt.tight_layout(); savefig(fig, 'fig6_net_directional_btc_es.png')

# Fig 7: network of net pairwise directional spillovers
netpair = base['NETPAIR']
thr = np.percentile(np.abs(netpair[np.triu_indices(10, 1)]), 60)
G = nx.DiGraph()
for a in ASSETS:
    G.add_node(a)
for i, a in enumerate(ASSETS):
    for j, b in enumerate(ASSETS):
        if i < j:
            v = netpair[i, j]
            if abs(v) > thr:
                if v > 0:
                    G.add_edge(a, b, weight=v)
                else:
                    G.add_edge(b, a, weight=-v)
order = ['ES', 'NQ', 'VIX', 'GC', 'CL', 'EURUSD', 'USDJPY', 'GBPUSD', 'BTC', 'ETH']
pos = {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)}
fig, ax = plt.subplots(figsize=(10.5, 8))
sizes = [300 + 25 * base['TO'][ASSETS.index(a)] for a in G.nodes()]
colors = [CLASS_COLOR[ASSET_CLASS[a]] for a in G.nodes()]
wts = np.array([G[u][v]['weight'] for u, v in G.edges()])
nx.draw_networkx_edges(G, pos, ax=ax, width=0.6 + 3.5 * wts / wts.max(), alpha=0.7, edge_color='#555555',
                       arrows=True, arrowsize=14, connectionstyle='arc3,rad=0.12', node_size=sizes)
nx.draw_networkx_nodes(G, pos, ax=ax, node_size=sizes, node_color=colors, alpha=0.92, edgecolors='white', linewidths=1.5)
nx.draw_networkx_labels(G, pos, ax=ax, font_size=9, font_weight='bold')
ax.legend(handles=[mpatches.Patch(color=c, label=k) for k, c in CLASS_COLOR.items()], loc='upper left', fontsize=8)
ax.set_title(f'Net pairwise directional spillovers above the 60th percentile (|net| > {thr:.1f} pp)\n'
             'arrow from net transmitter to net receiver; node size proportional to TO spillover')
ax.axis('off'); plt.tight_layout(); savefig(fig, 'fig7_network_graph.png')

# Fig 8: yearly BTC net and FROM with CIs
fig, axes = plt.subplots(1, 2, figsize=(12, 4.2))
ax = axes[0]
ax.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']],
            fmt='o-', color='#9467bd', capsize=3, label='BTC')
ax.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']],
            fmt='s-', color='#1f77b4', capsize=3, label='ES')
ax.axhline(0, color='k', lw=0.8, ls='--'); ax.set_ylabel('Net directional spillover (pp)'); ax.set_title('(a) Net spillovers by calendar year')
ax.legend()
ax = axes[1]
ax.plot(yr_df['Year'], yr_df['FROM_BTC_excl_ETH'], 'o-', color='#9467bd', label='BTC: FROM traditional assets (excl. ETH)')
ax.plot(yr_df['Year'], yr_df['BTC_from_ES'], 's--', color='#1f77b4', label='BTC $\\leftarrow$ ES')
ax.plot(yr_df['Year'], yr_df['ES_from_BTC'], '^--', color='#ff7f0e', label='ES $\\leftarrow$ BTC')
ax.set_ylabel('% of forecast-error variance'); ax.set_title('(b) Bitcoin–traditional-asset shares by year'); ax.legend(fontsize=8)
plt.tight_layout(); savefig(fig, 'fig8_yearly_crypto_integration.png')
log("figures written")

# ─── 13. Key numbers for the paper ────────────────────────────────────────────
key = {
    'T_panel': len(panel), 'T_original_panel': len(orig_panel), 'weekend_rows_original': n_weekend,
    'ffilled_rows_original': n_ffilled_rows, 'TSI': base['TSI'], 'TSI_lo': tsi_ci[0], 'TSI_hi': tsi_ci[1],
    'TSI_original_replicated': orig['TSI'], 'TSI_levels_trading_days': lvl_td['TSI'],
}
for k, a in enumerate(ASSETS):
    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]
    key[f'TO_{a}'] = base['TO'][k]; key[f'FROM_{a}'] = base['FROM'][k]; key[f'OWN_{a}'] = base['theta'][k, k]
key['BTC_ETH'] = base['theta'][7, 8]; key['ETH_BTC'] = base['theta'][8, 7]
key['VIX_from_ES'] = base['theta'][9, 0]; key['VIX_from_NQ'] = base['theta'][9, 1]
key['BTC_from_traditional'] = base['FROM'][7] - base['theta'][7, 8]
key['ETH_from_traditional'] = base['FROM'][8] - base['theta'][8, 7]
for w in (200, 250, 500):
    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()
    key[f'roll{w}_min_date'] = str(roll[w]['TSI'].idxmin().date()); key[f'roll{w}_max_date'] = str(roll[w]['TSI'].idxmax().date())
    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]
    key[f'roll{w}_NET_BTC_2025_mean'] = roll[w].loc['2025-01-01':'2025-12-31', 'NET_BTC'].mean()
    r17 = roll[w].loc['2017-01-01':'2017-12-31', 'NET_BTC']
    key[f'roll{w}_NET_BTC_2017_mean'] = r17.mean() if len(r17) else np.nan
    key[f'roll{w}_share_BTC_positive'] = (roll[w]['NET_BTC'] > 0).mean()
# rolling TSI around Feb 2018 and Mar 2020
key['roll200_2018-02-02'] = roll[200]['TSI'].get(pd.Timestamp('2018-02-02'), np.nan)
key['roll200_2018-02-06'] = roll[200]['TSI'].get(pd.Timestamp('2018-02-06'), np.nan)
key['roll200_2020-02-20'] = roll[200]['TSI'].get(pd.Timestamp('2020-02-20'), np.nan)
key['roll200_2020-03-31'] = roll[200]['TSI'].get(pd.Timestamp('2020-03-31'), np.nan)
key['roll200_2020-06-30'] = roll[200]['TSI'].get(pd.Timestamp('2020-06-30'), np.nan)
pd.Series(key).to_csv(os.path.join(RESULTS, 'rev_key_numbers.csv'), header=['value'])
log("key numbers saved; done")
print(pd.Series(key).to_string())