#!/usr/bin/env python3
"""
=============================================================================
SUPPLEMENTARY MATERIALS — Replication Code
"Cryptocurrency Market Maturation: Evidence from Bitcoin and Ethereum"

Data source: Crypto_Data_2026.xlsx (provided alongside this script)
Sheets used: BTC_Raw_Data, ETH_Raw_Data

Requirements:
    pip install pandas numpy scipy statsmodels matplotlib openpyxl

Usage:
    python supplementary_analysis.py

All tables are printed to console and all figures are saved to ./figures/
=============================================================================
"""
import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
from scipy import stats
from scipy.stats import kruskal, mannwhitneyu, spearmanr, pearsonr, jarque_bera
from statsmodels.regression.linear_model import OLS
from statsmodels.tools import add_constant
from statsmodels.tsa.stattools import adfuller
import json, os

# ── Output directory ──────────────────────────────────────────────────────
OUT = './figures'
os.makedirs(OUT, exist_ok=True)

plt.rcParams.update({
    'font.family': 'serif',
    'font.size': 10,
    'axes.labelsize': 11,
    'axes.titlesize': 12,
    'figure.dpi': 300,
    'savefig.dpi': 300,
    'savefig.bbox': 'tight',
    'figure.figsize': (10, 6),
})

WINDOW = 365   # 365-day rolling window (crypto trades 365 days/year)

# =============================================================================
# 1. LOAD DATA FROM EXCEL
# =============================================================================
DATA_FILE = 'Crypto_Data_2026.xlsx'

print("=" * 60)
print("LOADING DATA FROM", DATA_FILE)
print("=" * 60)

btc_raw = pd.read_excel(DATA_FILE, sheet_name='BTC_Raw_Data', parse_dates=['Date'])
eth_raw = pd.read_excel(DATA_FILE, sheet_name='ETH_Raw_Data', parse_dates=['Date'])

btc_raw = btc_raw.sort_values('Date').set_index('Date')
eth_raw = eth_raw.sort_values('Date').set_index('Date')

# Log returns from Close prices (recomputed for full reproducibility)
btc_ret = (np.log(btc_raw['Close'] / btc_raw['Close'].shift(1)) * 100).dropna()
eth_ret = (np.log(eth_raw['Close'] / eth_raw['Close'].shift(1)) * 100).dropna()

print(f"BTC: {len(btc_ret)} daily log-returns, "
      f"{btc_ret.index[0].date()} to {btc_ret.index[-1].date()}")
print(f"ETH: {len(eth_ret)} daily log-returns, "
      f"{eth_ret.index[0].date()} to {eth_ret.index[-1].date()}")

# =============================================================================
# 2. MARKET STATE CLASSIFICATION (SMA-50)
# =============================================================================
def classify_states(prices):
    """Bull / Bear / Neutral via 50-day SMA."""
    sma50 = prices.rolling(50).mean()
    sma50_rising = sma50 > sma50.shift(1)
    bull = (prices > sma50) & sma50_rising
    bear = (prices < sma50) & (~sma50_rising)
    states = pd.Series('Neutral', index=prices.index)
    states[bull] = 'Bull'
    states[bear] = 'Bear'
    return states

btc_states = classify_states(btc_raw['Close'])
eth_states = classify_states(eth_raw['Close'])

# 30-day rolling volatility
btc_vol30 = btc_ret.rolling(30).std()
eth_vol30 = eth_ret.rolling(30).std()

# =============================================================================
# 3. ROLLING RISK METRICS (365-day window)
# =============================================================================
print(f"\nComputing {WINDOW}-day rolling VaR / CVaR / MDD ...")

def rolling_risk(returns, window=WINDOW):
    """
    Compute rolling VaR (1%, 5%), CVaR (1%, 5%), and Maximum Drawdown.
    Returns five pandas Series aligned to the returns index.
    """
    n = len(returns)
    idx = returns.index
    var1  = pd.Series(np.nan, index=idx)
    var5  = pd.Series(np.nan, index=idx)
    cvar1 = pd.Series(np.nan, index=idx)
    cvar5 = pd.Series(np.nan, index=idx)
    mdd   = pd.Series(np.nan, index=idx)

    vals = returns.values
    for i in range(window, n):
        w  = vals[i - window:i]
        q1 = np.percentile(w, 1)
        q5 = np.percentile(w, 5)
        var1.iloc[i]  = -q1
        var5.iloc[i]  = -q5
        cvar1.iloc[i] = -w[w <= q1].mean() if np.any(w <= q1) else -q1
        cvar5.iloc[i] = -w[w <= q5].mean() if np.any(w <= q5) else -q5
        cum         = np.cumprod(1 + w / 100)
        running_max = np.maximum.accumulate(cum)
        dd          = (running_max - cum) / running_max
        mdd.iloc[i] = dd.max() * 100

    return var1, var5, cvar1, cvar5, mdd

btc_var1, btc_var5, btc_cvar1, btc_cvar5, btc_mdd = rolling_risk(btc_ret)
eth_var1, eth_var5, eth_cvar1, eth_cvar5, eth_mdd = rolling_risk(eth_ret)

print(f"BTC rolling metrics: {btc_var1.dropna().shape[0]} observations")
print(f"ETH rolling metrics: {eth_var1.dropna().shape[0]} observations")

# =============================================================================
# 4. TABLE 1: DESCRIPTIVE STATISTICS
# =============================================================================
def desc_stats(r):
    jb, jbp   = jarque_bera(r.dropna())
    adf_res   = adfuller(r.dropna(), maxlag=20)
    return {
        'Obs':              len(r),
        'Mean (%)':         r.mean(),
        'Median (%)':       r.median(),
        'Std Dev (%)':      r.std(),
        'Skewness':         r.skew(),
        'Excess Kurtosis':  r.kurtosis(),
        'Min (%)':          r.min(),
        'Max (%)':          r.max(),
        'JB Statistic':     jb,
        'JB p-value':       jbp,
        'ADF Statistic':    adf_res[0],
        'ADF p-value':      adf_res[1],
    }

bs = desc_stats(btc_ret)
es = desc_stats(eth_ret)

print("\n" + "=" * 60)
print("TABLE 1: DESCRIPTIVE STATISTICS OF DAILY LOG-RETURNS")
print("=" * 60)
print(f"{'Statistic':<25} {'BTC':>15} {'ETH':>15}")
print("-" * 55)
for k in bs:
    fmt = f"{bs[k]:>15.4f}" if isinstance(bs[k], float) else f"{bs[k]:>15}"
    fmt2 = f"{es[k]:>15.4f}" if isinstance(es[k], float) else f"{es[k]:>15}"
    print(f"{k:<25} {fmt} {fmt2}")

# =============================================================================
# 5. TABLE 2: ANNUAL STATE DISTRIBUTIONS
# =============================================================================
def annual_table(prices, returns, states, start_year, end_year=2026):
    rows = []
    for yr in range(start_year, end_year):
        mask = returns.index.year == yr
        if mask.sum() == 0:
            continue
        r  = returns[mask]
        s  = states.reindex(r.index).dropna()
        p  = prices.reindex(r.index).dropna()
        v  = r.rolling(30).std().dropna()
        rows.append({
            'Year':       yr,
            'N':          len(r),
            'P(Bull)':    (s == 'Bull').mean(),
            'P(Bear)':    (s == 'Bear').mean(),
            'P(Neutral)': (s == 'Neutral').mean(),
            'Mean Ret%':  r.mean(),
            'Avg Vol%':   v.mean() if len(v) > 0 else np.nan,
            'Min Price':  p.min(),
            'Max Price':  p.max(),
        })
    return pd.DataFrame(rows)

tbl2_btc = annual_table(btc_raw['Close'], btc_ret, btc_states, 2015)
tbl2_eth = annual_table(eth_raw['Close'], eth_ret, eth_states, 2018)

print("\n" + "=" * 60)
print("TABLE 2a: BTC ANNUAL STATE DISTRIBUTIONS")
print("=" * 60)
print(tbl2_btc.to_string(index=False, float_format=lambda x: f"{x:.4f}"))

print("\n" + "=" * 60)
print("TABLE 2b: ETH ANNUAL STATE DISTRIBUTIONS")
print("=" * 60)
print(tbl2_eth.to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# =============================================================================
# 6. TABLE 3: TREND TESTS (Newey-West HAC)
# =============================================================================
def trend_test(series, asset, metric):
    s    = series.dropna()
    t    = np.arange(len(s))
    X    = add_constant(t)
    nlag = int(len(s) ** (1 / 3))
    res  = OLS(s.values, X).fit(cov_type='HAC', cov_kwds={'maxlags': nlag})
    b1, tstat, pval, r2 = res.params[1], res.tvalues[1], res.pvalues[1], res.rsquared
    sig = '***' if pval < 0.01 else ('**' if pval < 0.05 else ('*' if pval < 0.10 else ''))
    return {
        'Asset':       asset,
        'Metric':      metric,
        'β₁ (×10⁻⁴)': b1 * 1e4,
        'NW t-stat':   tstat,
        'p-value':     pval,
        'R²':          r2,
        'N':           len(s),
        'Dir':         '↓' if b1 < 0 else '↑',
        'Sig':         sig,
    }

trend_rows = []
for asset, v1, v5, c1, c5, md in [
    ('BTC', btc_var1, btc_var5, btc_cvar1, btc_cvar5, btc_mdd),
    ('ETH', eth_var1, eth_var5, eth_cvar1, eth_cvar5, eth_mdd),
]:
    for name, s in [('VaR1', v1), ('VaR5', v5), ('CVaR1', c1), ('CVaR5', c5), ('MDD', md)]:
        trend_rows.append(trend_test(s, asset, name))

tbl3 = pd.DataFrame(trend_rows)
print("\n" + "=" * 60)
print("TABLE 3: LINEAR TREND TESTS ON ROLLING RISK METRICS (HAC s.e.)")
print("=" * 60)
print(tbl3.to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# =============================================================================
# 7. TABLE 4: RISK BY MARKET STATE & VOLATILITY REGIME
# =============================================================================
def risk_by_state(v1, v5, c1, c5, md, states, asset):
    rows = []
    for g in ['Bull', 'Bear', 'Neutral']:
        mask  = (states == g).reindex(v1.index, fill_value=False)
        valid = mask & v1.notna()
        if valid.sum() == 0:
            continue
        rows.append({
            'Asset': asset, 'State': g, 'N': int(valid.sum()),
            'VaR1 Mean':  v1[valid].mean(),  'VaR1 Median':  v1[valid].median(),  'VaR1 Std':  v1[valid].std(),
            'VaR5 Mean':  v5[valid].mean(),  'VaR5 Median':  v5[valid].median(),  'VaR5 Std':  v5[valid].std(),
            'CVaR1 Mean': c1[valid].mean(),  'CVaR1 Median': c1[valid].median(),  'CVaR1 Std': c1[valid].std(),
            'CVaR5 Mean': c5[valid].mean(),  'CVaR5 Median': c5[valid].median(),  'CVaR5 Std': c5[valid].std(),
            'MDD Mean':   md[valid].mean(),  'MDD Median':   md[valid].median(),  'MDD Std':   md[valid].std(),
        })
    return pd.DataFrame(rows)

tbl4a_btc = risk_by_state(btc_var1, btc_var5, btc_cvar1, btc_cvar5, btc_mdd, btc_states, 'BTC')
tbl4a_eth = risk_by_state(eth_var1, eth_var5, eth_cvar1, eth_cvar5, eth_mdd, eth_states, 'ETH')

print("\n" + "=" * 60)
print("TABLE 4a: RISK METRICS BY MARKET STATE")
print("=" * 60)
print(pd.concat([tbl4a_btc, tbl4a_eth]).to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# Kruskal-Wallis tests
print("\n  Kruskal-Wallis tests (Bull vs Bear vs Neutral):")
for asset, v1, v5, c1, c5, md, st in [
    ('BTC', btc_var1, btc_var5, btc_cvar1, btc_cvar5, btc_mdd, btc_states),
    ('ETH', eth_var1, eth_var5, eth_cvar1, eth_cvar5, eth_mdd, eth_states),
]:
    for name, series in [('VaR1', v1), ('VaR5', v5), ('CVaR1', c1), ('CVaR5', c5), ('MDD', md)]:
        s = series.dropna()
        groups = [s[(st.reindex(s.index) == g)].values
                  for g in ['Bull', 'Bear', 'Neutral']
                  if (st.reindex(s.index) == g).sum() > 0]
        if len(groups) >= 2:
            H, p = kruskal(*groups)
            sig = '***' if p < 0.01 else ('**' if p < 0.05 else ('*' if p < 0.10 else ''))
            print(f"    {asset} {name}: H={H:.2f}, p={p:.4e} {sig}")

# Volatility regime (High vs Low)
def risk_by_vol(v1, v5, c1, c5, md, vol30, asset):
    vol = vol30.reindex(v1.index)
    med = vol.median()
    rows = []
    for label, mask in [('High σ', vol > med), ('Low σ', vol <= med)]:
        valid = mask & v1.notna()
        rows.append({
            'Asset': asset, 'Regime': label, 'N': int(valid.sum()),
            'VaR1 Mean': v1[valid].mean(), 'VaR5 Mean': v5[valid].mean(),
            'CVaR1 Mean': c1[valid].mean(), 'CVaR5 Mean': c5[valid].mean(),
            'MDD Mean': md[valid].mean(),
        })
    return pd.DataFrame(rows), med

tbl4b_btc, btc_vol_med = risk_by_vol(btc_var1, btc_var5, btc_cvar1, btc_cvar5, btc_mdd, btc_vol30, 'BTC')
tbl4b_eth, eth_vol_med = risk_by_vol(eth_var1, eth_var5, eth_cvar1, eth_cvar5, eth_mdd, eth_vol30, 'ETH')

print("\n" + "=" * 60)
print("TABLE 4b: RISK METRICS BY VOLATILITY REGIME")
print(f"  (BTC median 30-day vol = {btc_vol_med:.4f}%, ETH = {eth_vol_med:.4f}%)")
print("=" * 60)
print(pd.concat([tbl4b_btc, tbl4b_eth]).to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# Mann-Whitney U tests
print("\n  Mann-Whitney U tests (High σ vs Low σ):")
for asset, v1, v5, c1, c5, md, vol30 in [
    ('BTC', btc_var1, btc_var5, btc_cvar1, btc_cvar5, btc_mdd, btc_vol30),
    ('ETH', eth_var1, eth_var5, eth_cvar1, eth_cvar5, eth_mdd, eth_vol30),
]:
    vol = vol30.reindex(v1.index)
    med = vol.median()
    for name, series in [('VaR1', v1), ('VaR5', v5), ('CVaR1', c1), ('CVaR5', c5), ('MDD', md)]:
        s  = series.dropna()
        hi = s[(vol.reindex(s.index) > med).fillna(False)]
        lo = s[(vol.reindex(s.index) <= med).fillna(False)]
        if len(hi) > 10 and len(lo) > 10:
            U, p = mannwhitneyu(hi, lo, alternative='two-sided')
            sig = '***' if p < 0.01 else ('**' if p < 0.05 else ('*' if p < 0.10 else ''))
            print(f"    {asset} {name}: U={U:.0f}, p={p:.4e} {sig}")

# =============================================================================
# 8. TABLE 5: SUB-PERIOD COMPARISON (Early vs Late)
# =============================================================================
def sub_period_test(series, asset, metric, early_end, late_start):
    s     = series.dropna()
    early = s[s.index.year <= early_end]
    late  = s[s.index.year >= late_start]
    if len(early) < 10 or len(late) < 10:
        return None
    U, p = mannwhitneyu(early, late, alternative='two-sided')
    sig  = '***' if p < 0.01 else ('**' if p < 0.05 else ('*' if p < 0.10 else ''))
    return {
        'Asset': asset, 'Metric': metric,
        'Early mean': early.mean(), 'Late mean': late.mean(),
        'Change (%)': (late.mean() - early.mean()) / early.mean() * 100,
        'MW U-stat': U, 'p-value': p, 'Sig': sig,
    }

# BTC: Early 2015–2019, Late 2020–2026
# ETH: Early 2018–2020, Late 2021–2026
tbl5_btc = pd.DataFrame([r for r in [
    sub_period_test(btc_var1,  'BTC', 'VaR1',  2019, 2020),
    sub_period_test(btc_var5,  'BTC', 'VaR5',  2019, 2020),
    sub_period_test(btc_cvar1, 'BTC', 'CVaR1', 2019, 2020),
    sub_period_test(btc_cvar5, 'BTC', 'CVaR5', 2019, 2020),
    sub_period_test(btc_mdd,   'BTC', 'MDD',   2019, 2020),
] if r])

tbl5_eth = pd.DataFrame([r for r in [
    sub_period_test(eth_var1,  'ETH', 'VaR1',  2020, 2021),
    sub_period_test(eth_var5,  'ETH', 'VaR5',  2020, 2021),
    sub_period_test(eth_cvar1, 'ETH', 'CVaR1', 2020, 2021),
    sub_period_test(eth_cvar5, 'ETH', 'CVaR5', 2020, 2021),
    sub_period_test(eth_mdd,   'ETH', 'MDD',   2020, 2021),
] if r])

print("\n" + "=" * 60)
print("TABLE 5: SUB-PERIOD COMPARISON (Early vs Late)")
print("  BTC: Early = 2015–2019, Late = 2020–2026")
print("  ETH: Early = 2018–2020, Late = 2021–2026")
print("=" * 60)
print(pd.concat([tbl5_btc, tbl5_eth]).to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# =============================================================================
# 9. TABLE 6: VOLATILITY REGIME × SUB-PERIOD
# =============================================================================
def regime_sub_test(series, vol30, asset, metric, early_end, late_start):
    s    = series.dropna()
    vol  = vol30.reindex(s.index)
    med  = vol.median()
    rows = []
    for regime, mask in [('High σ', vol > med), ('Low σ', vol <= med)]:
        s_r   = s[mask.reindex(s.index, fill_value=False)]
        early = s_r[s_r.index.year <= early_end]
        late  = s_r[s_r.index.year >= late_start]
        if len(early) < 5 or len(late) < 5:
            continue
        U, p = mannwhitneyu(early, late, alternative='two-sided')
        sig  = '***' if p < 0.01 else ('**' if p < 0.05 else ('*' if p < 0.10 else ''))
        rows.append({
            'Asset': asset, 'Regime': regime, 'Metric': metric,
            'Early mean': early.mean(), 'Late mean': late.mean(),
            'Change (%)': (late.mean() - early.mean()) / early.mean() * 100,
            'p-value': p, 'Sig': sig,
        })
    return rows

tbl6_rows = []
for m, s in [('VaR1', btc_var1), ('CVaR1', btc_cvar1), ('MDD', btc_mdd)]:
    tbl6_rows.extend(regime_sub_test(s, btc_vol30, 'BTC', m, 2019, 2020))
for m, s in [('VaR1', eth_var1), ('CVaR1', eth_cvar1), ('MDD', eth_mdd)]:
    tbl6_rows.extend(regime_sub_test(s, eth_vol30, 'ETH', m, 2020, 2021))

print("\n" + "=" * 60)
print("TABLE 6: VOLATILITY REGIME × SUB-PERIOD INTERACTION")
print("=" * 60)
print(pd.DataFrame(tbl6_rows).to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# =============================================================================
# 10. TABLE 7: EXCEEDANCE CORRELATIONS (BTC–ETH)
# =============================================================================
common = btc_ret.index.intersection(eth_ret.index)
br = btc_ret.reindex(common).dropna()
er = eth_ret.reindex(common).dropna()
common2 = br.index.intersection(er.index)
br, er = br.reindex(common2), er.reindex(common2)

full_pearson, _ = pearsonr(br, er)
print(f"\nFull-sample BTC–ETH Pearson ρ = {full_pearson:.4f}")

def exceedance_corr(r1, r2, quantiles=(75, 90, 95)):
    pooled = pd.concat([r1.abs(), r2.abs()])
    rows   = []
    for q in quantiles:
        theta  = np.percentile(pooled, q)
        lo     = (r1 < -theta) & (r2 < -theta)
        hi     = (r1 > theta)  & (r2 > theta)
        rho_lo = np.corrcoef(r1[lo], r2[lo])[0, 1] if lo.sum() > 5 else np.nan
        rho_hi = np.corrcoef(r1[hi], r2[hi])[0, 1] if hi.sum() > 5 else np.nan
        rows.append({
            'Quantile': f'{q}th',
            'ρ⁻ (lower)': rho_lo,
            'ρ⁺ (upper)': rho_hi,
            'ρ⁻ − ρ⁺': rho_lo - rho_hi if not (np.isnan(rho_lo) or np.isnan(rho_hi)) else np.nan,
            'N_lower': int(lo.sum()),
            'N_upper': int(hi.sum()),
        })
    return pd.DataFrame(rows)

tbl7 = exceedance_corr(br, er)
print("\n" + "=" * 60)
print("TABLE 7: EXCEEDANCE CORRELATIONS (full sample)")
print("=" * 60)
print(tbl7.to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# Sub-period exceedance
early_mask = common2.year <= 2020
late_mask  = common2.year >= 2021

print("\n" + "=" * 60)
print("TABLE 8: EXCEEDANCE CORRELATIONS BY SUB-PERIOD")
print("  Early = 2018–2020, Late = 2021–2026")
print("=" * 60)
for label, mask in [('Early (2018–2020)', early_mask), ('Late (2021–2026)', late_mask)]:
    print(f"\n  {label}:")
    df = exceedance_corr(br[mask], er[mask])
    print(df.to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# Correlation by market state
tbl9_rows = []
bs_c = btc_states.reindex(common2)
for state in ['Bull', 'Bear', 'Neutral']:
    mask = bs_c == state
    if mask.sum() < 30:
        continue
    r1, r2 = br[mask], er[mask]
    pr, pp = pearsonr(r1, r2)
    sr, sp = spearmanr(r1, r2)
    tbl9_rows.append({
        'State': state, 'N': int(mask.sum()),
        'Pearson ρ': pr, 'Pearson p': pp,
        'Spearman ρ': sr, 'Spearman p': sp,
    })

print("\n" + "=" * 60)
print("TABLE 9: BTC–ETH CORRELATIONS BY MARKET STATE")
print("=" * 60)
print(pd.DataFrame(tbl9_rows).to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# =============================================================================
# 11. FIGURES
# =============================================================================
print("\nGenerating figures ...")

def pct_fmt(x, pos): return f'{x:.0f}%'

# ── Figure 1: Rolling VaR ─────────────────────────────────────────────────
fig, axes = plt.subplots(2, 1, figsize=(10, 7))
for ax, asset, v1, v5, lbl in [
    (axes[0], 'Bitcoin',  btc_var1, btc_var5, 'A'),
    (axes[1], 'Ethereum', eth_var1, eth_var5, 'B'),
]:
    s1 = v1.dropna(); s5 = v5.reindex(s1.index)
    ax.fill_between(s1.index, s5, s1, alpha=0.2, color='steelblue')
    ax.plot(s1.index, s1, 'b-', lw=1,   label='VaR 1%')
    ax.plot(s5.index, s5, 'c-', lw=0.8, label='VaR 5%')
    z = np.polyfit(np.arange(len(s1)), s1.values, 1)
    ax.plot(s1.index, z[0] * np.arange(len(s1)) + z[1], 'k--', lw=1.2, label='Linear trend')
    ax.set_title(f'Panel {lbl}: {asset}', fontsize=11, fontweight='bold')
    ax.set_ylabel('VaR (%)')
    ax.yaxis.set_major_formatter(FuncFormatter(pct_fmt))
    ax.legend(fontsize=8, loc='upper right')
    ax.grid(alpha=0.3)
plt.suptitle(f'Figure 1. {WINDOW}-Day Rolling VaR (1% and 5%)', y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig(f'{OUT}/fig1_rolling_var.png')
plt.close()

# ── Figure 2: Rolling CVaR ────────────────────────────────────────────────
fig, axes = plt.subplots(2, 1, figsize=(10, 7))
for ax, asset, c1, c5, lbl in [
    (axes[0], 'Bitcoin',  btc_cvar1, btc_cvar5, 'A'),
    (axes[1], 'Ethereum', eth_cvar1, eth_cvar5, 'B'),
]:
    s1 = c1.dropna(); s5 = c5.reindex(s1.index)
    ax.fill_between(s1.index, s5, s1, alpha=0.2, color='coral')
    ax.plot(s1.index, s1, 'r-', lw=1, label='CVaR 1%')
    ax.plot(s5.index, s5, color='salmon', lw=0.8, label='CVaR 5%')
    z = np.polyfit(np.arange(len(s1)), s1.values, 1)
    ax.plot(s1.index, z[0] * np.arange(len(s1)) + z[1], 'k--', lw=1.2, label='Linear trend')
    ax.set_title(f'Panel {lbl}: {asset}', fontsize=11, fontweight='bold')
    ax.set_ylabel('CVaR (%)')
    ax.yaxis.set_major_formatter(FuncFormatter(pct_fmt))
    ax.legend(fontsize=8, loc='upper right')
    ax.grid(alpha=0.3)
plt.suptitle(f'Figure 2. {WINDOW}-Day Rolling CVaR (1% and 5%)', y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig(f'{OUT}/fig2_rolling_cvar.png')
plt.close()

# ── Figure 3: Rolling MDD ─────────────────────────────────────────────────
fig, axes = plt.subplots(2, 1, figsize=(10, 7))
for ax, asset, md, lbl in [
    (axes[0], 'Bitcoin',  btc_mdd, 'A'),
    (axes[1], 'Ethereum', eth_mdd, 'B'),
]:
    s = md.dropna()
    ax.fill_between(s.index, 0, s, alpha=0.3, color='purple')
    ax.plot(s.index, s, color='purple', lw=0.8)
    z = np.polyfit(np.arange(len(s)), s.values, 1)
    ax.plot(s.index, z[0] * np.arange(len(s)) + z[1], 'k--', lw=1.2, label='Linear trend')
    ax.set_title(f'Panel {lbl}: {asset}', fontsize=11, fontweight='bold')
    ax.set_ylabel('MDD (%)')
    ax.yaxis.set_major_formatter(FuncFormatter(pct_fmt))
    ax.legend(fontsize=8)
    ax.grid(alpha=0.3)
plt.suptitle(f'Figure 3. {WINDOW}-Day Rolling Maximum Drawdown', y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig(f'{OUT}/fig3_rolling_mdd.png')
plt.close()

# ── Figure 4: Box plots by market state ──────────────────────────────────
fig, axes = plt.subplots(2, 3, figsize=(12, 7))
clr_map = {'Bull': '#2ecc71', 'Bear': '#e74c3c', 'Neutral': '#95a5a6'}
for row, (asset, v1, c1, md, st) in enumerate([
    ('BTC', btc_var1, btc_cvar1, btc_mdd, btc_states),
    ('ETH', eth_var1, eth_cvar1, eth_mdd, eth_states),
]):
    for col, (metric, series) in enumerate([('VaR 1%', v1), ('CVaR 1%', c1), ('MDD', md)]):
        ax  = axes[row, col]
        dat, lbls, clrs = [], [], []
        for state in ['Bull', 'Neutral', 'Bear']:
            s    = series.dropna()
            vals = s[st.reindex(s.index) == state]
            if len(vals) > 0:
                dat.append(vals.values); lbls.append(state); clrs.append(clr_map[state])
        bp = ax.boxplot(dat, labels=lbls, patch_artist=True,
                        flierprops={'marker': 'D', 'markersize': 2, 'alpha': 0.3})
        for patch, c in zip(bp['boxes'], clrs):
            patch.set_facecolor(c); patch.set_alpha(0.6)
        ax.set_title(f'{asset} — {metric}', fontsize=10, fontweight='bold')
        ax.set_ylabel('%')
        ax.grid(alpha=0.3, axis='y')
plt.suptitle('Figure 4. Risk Metrics by Market State', y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig(f'{OUT}/fig4_boxplots_state.png')
plt.close()

# ── Figure 5: Volatility regime × sub-period ─────────────────────────────
fig, axes = plt.subplots(2, 3, figsize=(12, 7))
for row, (asset, v1, c1, md, vol30, ee, ls) in enumerate([
    ('BTC', btc_var1, btc_cvar1, btc_mdd, btc_vol30, 2019, 2020),
    ('ETH', eth_var1, eth_cvar1, eth_mdd, eth_vol30, 2020, 2021),
]):
    med = vol30.reindex(v1.dropna().index).median()
    for col, (metric, series) in enumerate([('VaR 1%', v1), ('CVaR 1%', c1), ('MDD', md)]):
        ax   = axes[row, col]
        s    = series.dropna()
        va   = vol30.reindex(s.index)
        hi   = va > med; lo = va <= med
        early = s.index.year <= ee; late = s.index.year >= ls
        vals = [s[hi & early].mean(), s[hi & late].mean(),
                s[lo & early].mean(), s[lo & late].mean()]
        bars = ax.bar([0, 1, 3, 4], vals,
                      color=['#e74c3c', '#ff9999', '#3498db', '#99ccff'],
                      width=0.8, edgecolor='black', linewidth=0.5)
        ax.set_xticks([0.5, 3.5])
        ax.set_xticklabels(['High σ', 'Low σ'], fontsize=9)
        ax.set_title(f'{asset} — {metric}', fontsize=10, fontweight='bold')
        ax.set_ylabel('%')
        ax.grid(alpha=0.3, axis='y')
        if col == 0 and row == 0:
            ax.legend([bars[0], bars[1]], ['Early period', 'Late period'],
                      fontsize=8, loc='upper right')
plt.suptitle('Figure 5. Risk Metrics by Volatility Regime and Sub-Period', y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig(f'{OUT}/fig5_regime_subperiod.png')
plt.close()

# ── Figure 6: Rolling correlations ───────────────────────────────────────
fig, axes = plt.subplots(2, 1, figsize=(10, 7))

roll_p = br.rolling(90).corr(er).dropna()
roll_s_vals = [spearmanr(br.iloc[i-90:i], er.iloc[i-90:i])[0] for i in range(90, len(br))]
roll_s = pd.Series(roll_s_vals, index=br.index[90:])

ax = axes[0]
ax.plot(roll_p.index, roll_p, 'b-', lw=0.8, label='Pearson (90-day)')
ax.plot(roll_s.index, roll_s, color='orange', lw=0.6, alpha=0.7, label='Spearman (90-day)')
ax.axhline(full_pearson, color='navy', ls='--', lw=1,
           label=f'Full-sample ρ = {full_pearson:.3f}')
bs_c2 = btc_states.reindex(roll_p.index)
for state, color in [('Bull', 'green'), ('Bear', 'red')]:
    ax.fill_between(roll_p.index, -0.2, 1.1,
                    where=(bs_c2 == state), color=color, alpha=0.08)
ax.set_ylim(-0.2, 1.05)
ax.set_title('Panel A: 90-Day Rolling Correlations', fontweight='bold')
ax.set_ylabel('Correlation')
ax.legend(fontsize=8, loc='lower left')
ax.grid(alpha=0.3)

ax = axes[1]
rlo, rhi, dates_exc = [], [], []
for i in range(WINDOW, len(br)):
    b, e = br.iloc[i-WINDOW:i], er.iloc[i-WINDOW:i]
    theta = np.percentile(pd.concat([b.abs(), e.abs()]), 90)
    lo = (b < -theta) & (e < -theta)
    hi = (b > theta)  & (e > theta)
    rlo.append(np.corrcoef(b[lo], e[lo])[0, 1] if lo.sum() > 5 else np.nan)
    rhi.append(np.corrcoef(b[hi], e[hi])[0, 1] if hi.sum() > 5 else np.nan)
    dates_exc.append(br.index[i])
ax.plot(dates_exc, rlo, 'r-', lw=1, label='ρ⁻ (lower tail)')
ax.plot(dates_exc, rhi, 'g-', lw=1, label='ρ⁺ (upper tail)')
ax.axhline(0, color='black', lw=0.5)
ax.set_title(f'Panel B: {WINDOW}-Day Rolling Exceedance Correlations (90th percentile)',
             fontweight='bold')
ax.set_ylabel('Exceedance Correlation')
ax.legend(fontsize=8)
ax.grid(alpha=0.3)

plt.suptitle('Figure 6. Rolling BTC–ETH Correlations', y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig(f'{OUT}/fig6_rolling_correlations.png')
plt.close()

# ── Figure 7: Scatter with exceedance zones ───────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
for ax, label, mask in [
    (axes[0], 'Early (2018–2020)', early_mask),
    (axes[1], 'Late (2021–2026)',  late_mask),
]:
    b, e   = br[mask], er[mask]
    theta  = np.percentile(pd.concat([b.abs(), e.abs()]), 90)
    lo     = (b < -theta) & (e < -theta)
    hi     = (b > theta)  & (e > theta)
    other  = ~lo & ~hi
    ax.scatter(b[other], e[other], s=3,  alpha=0.15, color='grey')
    ax.scatter(b[lo],    e[lo],    s=15, alpha=0.6,  color='red',
               label=f'Lower tail (n={lo.sum()})', zorder=5)
    ax.scatter(b[hi],    e[hi],    s=15, alpha=0.6,  color='green',
               label=f'Upper tail (n={hi.sum()})', zorder=5)
    ax.axhline(0, color='black', lw=0.3)
    ax.axvline(0, color='black', lw=0.3)
    ax.set_xlabel('BTC return (%)')
    ax.set_ylabel('ETH return (%)')
    ax.set_title(label, fontweight='bold')
    ax.legend(fontsize=8, loc='upper left')
    ax.grid(alpha=0.2)
    lim = max(b.abs().max(), e.abs().max()) * 1.05
    ax.set_xlim(-lim, lim); ax.set_ylim(-lim, lim)
plt.suptitle('Figure 7. BTC–ETH Return Scatter with Exceedance Zones (90th percentile)',
             y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig(f'{OUT}/fig7_scatter_exceedance.png')
plt.close()

# ── Appendix Figure A1: Price series with states ──────────────────────────
fig, axes = plt.subplots(2, 1, figsize=(10, 7))
for ax, asset, prices, states, lbl in [
    (axes[0], 'Bitcoin',  btc_raw['Close'], btc_states, 'A'),
    (axes[1], 'Ethereum', eth_raw['Close'], eth_states, 'B'),
]:
    ax.semilogy(prices.index, prices, 'k-', lw=0.8)
    for state, color in [('Bull', 'green'), ('Bear', 'red')]:
        m = states.reindex(prices.index) == state
        ax.fill_between(prices.index,
                        prices.min() * 0.8, prices.max() * 1.2,
                        where=m, color=color, alpha=0.1)
    ax.set_title(f'Panel {lbl}: {asset}', fontweight='bold')
    ax.set_ylabel('Price (USD, log scale)')
    ax.grid(alpha=0.3)
plt.suptitle('Appendix Figure A1. Price Series with Market State Regimes',
             y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig(f'{OUT}/figA1_price_states.png')
plt.close()

# ── Appendix Figure A2: 30-day rolling volatility ─────────────────────────
fig, axes = plt.subplots(2, 1, figsize=(10, 6))
for ax, asset, vol, med, lbl in [
    (axes[0], 'Bitcoin',  btc_vol30, btc_vol_med, 'A'),
    (axes[1], 'Ethereum', eth_vol30, eth_vol_med, 'B'),
]:
    v  = vol.dropna()
    hi = v > med
    ax.fill_between(v.index, 0, v, where=hi,  color='red',  alpha=0.3, label='High uncertainty')
    ax.fill_between(v.index, 0, v, where=~hi, color='blue', alpha=0.3, label='Low uncertainty')
    ax.axhline(med, color='black', ls='--', lw=1, label=f'Median = {med:.2f}%')
    ax.set_title(f'Panel {lbl}: {asset}', fontweight='bold')
    ax.set_ylabel('30-day σ (%)')
    ax.legend(fontsize=8)
    ax.grid(alpha=0.3)
plt.suptitle('Appendix Figure A2. 30-Day Rolling Volatility with Uncertainty Regimes',
             y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig(f'{OUT}/figA2_rolling_vol.png')
plt.close()

# =============================================================================
# 12. SAVE SUMMARY TO JSON
# =============================================================================
summary = {
    'data_source':    DATA_FILE,
    'rolling_window': WINDOW,
    'btc_period':     [str(btc_ret.index[0].date()), str(btc_ret.index[-1].date())],
    'eth_period':     [str(eth_ret.index[0].date()), str(eth_ret.index[-1].date())],
    'btc_obs':        int(bs['Obs']),
    'eth_obs':        int(es['Obs']),
    'btc_vol_median': float(btc_vol_med),
    'eth_vol_median': float(eth_vol_med),
    'btc_eth_pearson_full': float(full_pearson),
}
with open('supplementary_results.json', 'w') as f:
    json.dump(summary, f, indent=2)

print(f"\nAll figures saved to ./{OUT}/")
print("Summary saved to supplementary_results.json")
print("=" * 60)
print("REPLICATION COMPLETE")
print("=" * 60)
