scripts/run_ablation_study.py

"""Architecture and feature ablation on the best target, with naive baselines.

Answers the two questions the thesis methodology asks (3.1.5 baselines, 3.1.9
ablation) on a target that demonstrably carries signal, so that a component
contributing nothing is a finding rather than a failure to learn anything at
all.

**The ladder.** Every rung sees identical folds, identical features and
identical scaling, so a difference between rungs is attributable to the rung:

    majority          always answer the more common class
    persistence       carry the previous period's answer forward
    logistic          one linear layer on the latest row
    MLP               non-linear, still one row: no temporal structure
    GRU               recurrence over the window
    Transformer       attention over the window, no latent variable
    VAE+Transformer   the deployed architecture

The two baselines come first on purpose. On this project's own screening a
target reached 94.71% accuracy while losing to its majority class, and traded
volume had the highest ROC-AUC of eight candidates while losing to persistence.
A ladder without naive rungs can rank seven architectures and still not
establish that any of them works.

**The feature arm** repeats the deployed architecture on market data alone,
market plus Fear & Greed, and market plus on-chain, which is the historic
versus sentiment comparison. Both external sources are daily and are lagged a
full day, so no feature describes the hour it scores.

**Implementation note.** The models are defined here rather than imported from
src/pipeline, which pulls in scikit-learn estimators. That keeps this script
runnable where those estimators are unavailable, and keeps the ablation
self-contained: a reader can see every architecture being compared in one file.

Usage (Colab for the full grid, roughly 2 hours on an A100):
    python scripts/run_ablation_study.py
    python scripts/run_ablation_study.py --symbols BTCUSDT --folds 3
    python scripts/run_ablation_study.py --max-rows 6000 --epochs 2   # smoke test
    python scripts/run_ablation_study.py --target regime --symbols BTCUSDT --seq-len 168 --skip-features
        # the deployed regime question, with HAR and GARCH(1,1) rungs on the same folds
"""
from __future__ import annotations

import argparse
import ast
import json
import sys
import warnings
from pathlib import Path

import numpy as np
import pandas as pd
import torch
from torch import nn

ROOT = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(ROOT))
warnings.filterwarnings("ignore")

from sklearn.metrics import accuracy_score, matthews_corrcoef, roc_auc_score  # noqa: E402

DATA = ROOT / "outputs" / "sentiment_replication"
EPS = 1e-12
DEVICE = torch.device("cuda" if torch.cuda.is_available() else "cpu")

# horizons from scripts/screen_horizon_and_spillover.py
HORIZON = {"BTCUSDT": 12, "ETHUSDT": 12, "ADAUSDT": 1}


# ----------------------------------------------------------------- the models
class Logistic(nn.Module):
    """One linear layer on the most recent row. The simplest learned rung."""

    def __init__(self, n_features: int, seq_len: int):
        super().__init__()
        self.net = nn.Linear(n_features, 1)

    def forward(self, x):                      # x: (batch, seq, features)
        return self.net(x[:, -1, :]).squeeze(-1)


class MLP(nn.Module):
    """Non-linear but still one row: isolates non-linearity from memory."""

    def __init__(self, n_features: int, seq_len: int, hidden: int = 96):
        super().__init__()
        self.net = nn.Sequential(nn.Linear(n_features, hidden), nn.ReLU(),
                                 nn.Dropout(0.15), nn.Linear(hidden, hidden // 2),
                                 nn.ReLU(), nn.Linear(hidden // 2, 1))

    def forward(self, x):
        return self.net(x[:, -1, :]).squeeze(-1)


class GRUModel(nn.Module):
    """Recurrence over the window. The first rung that can use history."""

    def __init__(self, n_features: int, seq_len: int, hidden: int = 64):
        super().__init__()
        self.rnn = nn.GRU(n_features, hidden, batch_first=True)
        self.head = nn.Sequential(nn.Dropout(0.15), nn.Linear(hidden, 1))

    def forward(self, x):
        out, _ = self.rnn(x)
        return self.head(out[:, -1, :]).squeeze(-1)


class TransformerModel(nn.Module):
    """Attention over the window, with no latent variable.

    The rung that isolates what the VAE contributes: everything else is shared
    with the deployed architecture.
    """

    def __init__(self, n_features: int, seq_len: int, d_model: int = 64,
                 nhead: int = 4, layers: int = 2, dropout: float = 0.15):
        super().__init__()
        self.project = nn.Linear(n_features, d_model)
        self.position = nn.Parameter(torch.zeros(1, seq_len, d_model))
        encoder = nn.TransformerEncoderLayer(d_model, nhead, d_model * 2, dropout,
                                             batch_first=True)
        self.encoder = nn.TransformerEncoder(encoder, layers)
        self.head = nn.Sequential(nn.Dropout(dropout), nn.Linear(d_model, 1))

    def forward(self, x):
        h = self.encoder(self.project(x) + self.position)
        return self.head(h[:, -1, :]).squeeze(-1)


class VAETransformer(nn.Module):
    """The deployed architecture: a variational latent, then attention.

    Returns the logits and the reconstruction terms, because the training
    objective needs all three; only the logits are a prediction.
    """

    def __init__(self, n_features: int, seq_len: int, d_model: int = 64,
                 nhead: int = 4, layers: int = 2, latent: int = 24,
                 dropout: float = 0.15):
        super().__init__()
        self.encoder_in = nn.Sequential(nn.Linear(n_features, d_model), nn.ReLU())
        self.mu = nn.Linear(d_model, latent)
        self.logvar = nn.Linear(d_model, latent)
        self.decoder = nn.Sequential(nn.Linear(latent, d_model), nn.ReLU(),
                                     nn.Linear(d_model, n_features))
        self.project = nn.Linear(latent, d_model)
        self.position = nn.Parameter(torch.zeros(1, seq_len, d_model))
        encoder = nn.TransformerEncoderLayer(d_model, nhead, d_model * 2, dropout,
                                             batch_first=True)
        self.attn = nn.TransformerEncoder(encoder, layers)
        self.head = nn.Sequential(nn.Dropout(dropout), nn.Linear(d_model, 1))

    def forward(self, x):
        h = self.encoder_in(x)
        mu, logvar = self.mu(h), self.logvar(h).clamp(-8, 8)
        z = mu + torch.randn_like(mu) * torch.exp(0.5 * logvar) if self.training else mu
        recon = self.decoder(z)
        out = self.attn(self.project(z) + self.position)
        return self.head(out[:, -1, :]).squeeze(-1), recon, mu, logvar


class LSTMModel(nn.Module):
    """Recurrence with gating. Baseline 4 in the methodology alongside the GRU."""

    def __init__(self, n_features: int, seq_len: int, hidden: int = 64):
        super().__init__()
        self.rnn = nn.LSTM(n_features, hidden, batch_first=True)
        self.head = nn.Sequential(nn.Dropout(0.15), nn.Linear(hidden, 1))

    def forward(self, x):
        out, _ = self.rnn(x)
        return self.head(out[:, -1, :]).squeeze(-1)


class CNNEncoder(nn.Module):
    """1D convolution over the time axis: local shape, no long-range view.

    Shared by every CNN rung so the three of them differ only in what follows
    the convolution, which is what an ablation of the CNN is supposed to isolate.
    """

    def __init__(self, n_features: int, channels: int = 64):
        super().__init__()
        self.net = nn.Sequential(
            nn.Conv1d(n_features, channels, kernel_size=5, padding=2), nn.ReLU(),
            nn.Conv1d(channels, channels, kernel_size=3, padding=1), nn.ReLU())
        self.channels = channels

    def forward(self, x):                      # (batch, seq, features)
        return self.net(x.transpose(1, 2)).transpose(1, 2)


class CNNModel(nn.Module):
    """CNN alone, pooled over time. Baseline 4 in the methodology."""

    def __init__(self, n_features: int, seq_len: int, channels: int = 64):
        super().__init__()
        self.cnn = CNNEncoder(n_features, channels)
        self.head = nn.Sequential(nn.Dropout(0.15), nn.Linear(channels, 1))

    def forward(self, x):
        return self.head(self.cnn(x)[:, -1, :]).squeeze(-1)


class CNNTransformer(nn.Module):
    """Rung (b): convolutional front-end, then attention. No latent variable."""

    def __init__(self, n_features: int, seq_len: int, d_model: int = 64,
                 nhead: int = 4, layers: int = 2, dropout: float = 0.15):
        super().__init__()
        self.cnn = CNNEncoder(n_features, d_model)
        self.position = nn.Parameter(torch.zeros(1, seq_len, d_model))
        encoder = nn.TransformerEncoderLayer(d_model, nhead, d_model * 2, dropout,
                                             batch_first=True)
        self.encoder = nn.TransformerEncoder(encoder, layers)
        self.head = nn.Sequential(nn.Dropout(dropout), nn.Linear(d_model, 1))

    def forward(self, x):
        h = self.encoder(self.cnn(x) + self.position)
        return self.head(h[:, -1, :]).squeeze(-1)


class CNNVAETransformer(nn.Module):
    """Rung (d) and baseline 6: the architecture the thesis actually proposes.

    CNN for local feature extraction, a variational latent, then attention.
    Until now no rung in this study contained a CNN at all, so the proposed
    model had never been measured.
    """

    def __init__(self, n_features: int, seq_len: int, d_model: int = 64,
                 nhead: int = 4, layers: int = 2, latent: int = 24,
                 dropout: float = 0.15):
        super().__init__()
        self.cnn = CNNEncoder(n_features, d_model)
        self.mu = nn.Linear(d_model, latent)
        self.logvar = nn.Linear(d_model, latent)
        self.decoder = nn.Sequential(nn.Linear(latent, d_model), nn.ReLU(),
                                     nn.Linear(d_model, n_features))
        self.project = nn.Linear(latent, d_model)
        self.position = nn.Parameter(torch.zeros(1, seq_len, d_model))
        encoder = nn.TransformerEncoderLayer(d_model, nhead, d_model * 2, dropout,
                                             batch_first=True)
        self.attn = nn.TransformerEncoder(encoder, layers)
        self.head = nn.Sequential(nn.Dropout(dropout), nn.Linear(d_model, 1))

    def forward(self, x):
        h = self.cnn(x)
        mu, logvar = self.mu(h), self.logvar(h).clamp(-8, 8)
        z = mu + torch.randn_like(mu) * torch.exp(0.5 * logvar) if self.training else mu
        recon = self.decoder(z)
        out = self.attn(self.project(z) + self.position)
        return self.head(out[:, -1, :]).squeeze(-1), recon, mu, logvar


ARCHITECTURES = {
    "logistic": Logistic,
    "MLP": MLP,
    "CNN": CNNModel,
    "LSTM": LSTMModel,
    "GRU": GRUModel,
    "Transformer": TransformerModel,
    "CNN+Transformer": CNNTransformer,
    "VAE+Transformer": VAETransformer,
    "CNN+VAE+Transformer": CNNVAETransformer,
}


# ------------------------------------------------------------------- the data
def logret(frame) -> np.ndarray:
    close = frame["close"].to_numpy(dtype=float)
    return np.diff(np.log(np.clip(close, EPS, None)), prepend=np.nan)


def volatility_rises(frame, horizon: int) -> np.ndarray:
    """Next `horizon` hours more volatile than the last `horizon`."""
    r = logret(frame)
    n = len(r)
    forward = np.full((n, horizon), np.nan)
    for k in range(1, horizon + 1):
        forward[: n - k, k - 1] = r[k:]
    ahead = np.sqrt(np.nanmean(forward ** 2, axis=1))
    behind = (pd.Series(r ** 2).rolling(horizon, min_periods=horizon)
              .mean().pow(0.5).to_numpy())
    out = np.full(n, np.nan)
    ok = np.isfinite(ahead) & np.isfinite(behind)
    out[ok] = (ahead[ok] > behind[ok]).astype(float)
    return out


def daily_block(path: Path) -> pd.DataFrame:
    """A daily CSV lagged one full day, so nothing describes the hour it scores."""
    frame = pd.read_csv(path, parse_dates=["date"])
    frame["date"] = pd.to_datetime(frame["date"], utc=True).dt.normalize()
    frame = frame.set_index("date").sort_index()
    frame.index = frame.index + pd.Timedelta(days=1)
    return frame


def join_hourly(data: pd.DataFrame, block: pd.DataFrame) -> np.ndarray:
    days = pd.to_datetime(data["timestamp"], utc=True).dt.normalize()
    return (block.reindex(block.index.union(days.unique())).sort_index()
            .ffill().reindex(days).to_numpy(dtype=float))


def _fetch_klines(symbol: str):
    """Binance klines, without importing the module that owns the fetcher.

    scripts/run_volatility_experiment.py imports a scikit-learn estimator at
    module level, and this ablation needs none of them -- only the download. On
    hosts where those estimators are unavailable, importing the module to reach
    one function would take the whole script down, so the function is loaded
    from the file without executing the module body.
    """
    import importlib.util

    spec = importlib.util.spec_from_file_location(
        "_kl", ROOT / "scripts" / "run_volatility_experiment.py")
    module = importlib.util.module_from_spec(spec)
    source = (ROOT / "scripts" / "run_volatility_experiment.py").read_text(encoding="utf-8")
    tree = ast.parse(source)
    keep = [n for n in tree.body
            if not (isinstance(n, (ast.Import, ast.ImportFrom))
                    and "sklearn" in ast.unparse(n))]
    exec(compile(ast.Module(body=keep, type_ignores=[]), "<kl>", "exec"), module.__dict__)
    cache = ROOT / "data" / "cache" / f"{symbol}_1h.parquet"
    cache.parent.mkdir(parents=True, exist_ok=True)
    return module.fetch_klines("2019-01-01", cache, symbol=symbol)


def longest_contiguous(frame: pd.DataFrame) -> pd.DataFrame:
    """The longest run of consecutive hourly bars, so no window spans a gap."""
    stamps = pd.to_datetime(frame["timestamp"], utc=True)
    breaks = ([0] + list(np.flatnonzero(stamps.diff().to_numpy() != np.timedelta64(1, "h"))[1:])
              + [len(frame)])
    spans = [(breaks[i], breaks[i + 1]) for i in range(len(breaks) - 1)]
    lo, hi = max(spans, key=lambda s: s[1] - s[0])
    return frame.iloc[lo:hi].reset_index(drop=True)


REGIME_HORIZON = 24
HAR_WINDOWS = (24, 168, 720)


def load(symbol: str, max_rows: int | None, target: str = "rises"):
    """Features and label. `target="regime"` is the deployed volatility-regime
    question (next 24h above the trailing 90-day median), with the extra
    columns the HAR, GARCH and regime-persistence rungs need."""
    from src.volatility import add_volatility_target, build_volatility_features

    raw = _fetch_klines(symbol)
    featured, features = build_volatility_features(raw)
    if target == "regime":
        featured = add_volatility_target(featured, REGIME_HORIZON)
        r2 = featured["log_return"] ** 2
        for w in HAR_WINDOWS:
            featured[f"rv{w}"] = r2.rolling(w, min_periods=int(w * 0.9)).mean()
        featured = featured.assign(y=featured["target_volatility_regime"].astype(float))
        extra = ["forward_vol", "threshold_causal"] + [f"rv{w}" for w in HAR_WINDOWS]
    else:
        featured = featured.assign(y=volatility_rises(featured, HORIZON[symbol]))
        extra = []
    data = featured.dropna(subset=features + ["y"] + extra).reset_index(drop=True)
    data = longest_contiguous(data)
    if max_rows:
        data = data.iloc[-max_rows:].reset_index(drop=True)
    return data, features


# ------------------------------------------------------------------ evaluation
def windows(X: np.ndarray, index: np.ndarray, seq_len: int) -> torch.Tensor:
    """Stack the `seq_len` rows ending at each index, without copying the lot.

    Materialising every window for 30,000 rows at 168 hours would be hundreds of
    megabytes; this builds only the batch being used.
    """
    idx = index[:, None] - np.arange(seq_len - 1, -1, -1)[None, :]
    return torch.tensor(X[idx], dtype=torch.float32)


def train_one(name, X, y, train_idx, valid_idx, test_idx, seq_len, epochs, seed):
    torch.manual_seed(seed)
    model = ARCHITECTURES[name](X.shape[1], seq_len).to(DEVICE)
    opt = torch.optim.AdamW(model.parameters(), lr=1e-3, weight_decay=1e-4)
    loss_fn = nn.BCEWithLogitsLoss()
    best_state, best_score, patience = None, -np.inf, 0

    for _ in range(epochs):
        model.train()
        order = np.random.default_rng(seed).permutation(len(train_idx))
        for start in range(0, len(order), 256):
            batch = train_idx[order[start:start + 256]]
            xb = windows(X, batch, seq_len).to(DEVICE)
            yb = torch.tensor(y[batch], dtype=torch.float32, device=DEVICE)
            opt.zero_grad()
            out = model(xb)
            if isinstance(out, tuple):
                logits, recon, mu, logvar = out
                kl = -0.5 * torch.mean(1 + logvar - mu.pow(2) - logvar.exp())
                loss = loss_fn(logits, yb) + 0.1 * nn.functional.mse_loss(recon, xb) \
                    + 1e-4 * kl
            else:
                loss = loss_fn(out, yb)
            loss.backward()
            nn.utils.clip_grad_norm_(model.parameters(), 1.0)
            opt.step()

        score = roc_auc_score(y[valid_idx], predict(model, X, valid_idx, seq_len))
        if score > best_score:
            best_score, best_state, patience = score, \
                {k: v.detach().clone() for k, v in model.state_dict().items()}, 0
        else:
            patience += 1
            if patience >= 3:
                break

    if best_state:
        model.load_state_dict(best_state)
    return predict(model, X, test_idx, seq_len), predict(model, X, valid_idx, seq_len)


@torch.no_grad()
def predict(model, X, index, seq_len) -> np.ndarray:
    model.eval()
    out = []
    for start in range(0, len(index), 512):
        xb = windows(X, index[start:start + 512], seq_len).to(DEVICE)
        logits = model(xb)
        if isinstance(logits, tuple):
            logits = logits[0]
        out.append(torch.sigmoid(logits).cpu().numpy())
    return np.concatenate(out)


THRESHOLDS = np.round(np.arange(0.05, 0.951, 0.01), 2)


def tune_threshold(valid_truth, valid_probs) -> float:
    """The cut-off with the best MCC on the VALIDATION block (the last fifth of
    training), never on the test fold. Every rung gets the same rule, so a
    model whose probabilities are well ordered but shifted is not penalised
    by an arbitrary 0.5."""
    if valid_probs is None or len(np.unique(valid_truth)) < 2:
        return 0.5
    return float(max(THRESHOLDS, key=lambda t: matthews_corrcoef(
        valid_truth, (valid_probs >= t).astype(int))))


def score(truth, probs, threshold: float = 0.5) -> dict:
    predictions = (probs >= threshold).astype(int)
    at_half = (probs >= 0.5).astype(int)
    majority = float(max(truth.mean(), 1 - truth.mean()))
    accuracy = float(accuracy_score(truth, predictions))
    return {"roc_auc": float(roc_auc_score(truth, probs)), "accuracy": accuracy,
            "mcc": float(matthews_corrcoef(truth, predictions)),
            "accuracy_at_0.5": float(accuracy_score(truth, at_half)),
            "mcc_at_0.5": float(matthews_corrcoef(truth, at_half)),
            "threshold": threshold, "majority": majority, "margin": accuracy - majority}


def auc_fast(truth: np.ndarray, probs: np.ndarray) -> float:
    order = np.argsort(probs, kind="mergesort")
    ranks = np.empty(len(probs))
    ranks[order] = np.arange(1, len(probs) + 1)
    pos = truth.sum()
    return float((ranks[truth == 1].sum() - pos * (pos + 1) / 2) / (pos * (len(truth) - pos)))


def auc_difference(truth, a, b, rng, block: int = 168, n: int = 2000) -> dict:
    """AUC(a) - AUC(b) on pooled test hours, with a moving-block bootstrap
    (one-week blocks, since neighbouring hours share their 24h outcome)."""
    starts = np.arange(len(truth) - block)
    diffs = []
    for _ in range(n):
        idx = np.concatenate([np.arange(s, s + block)
                              for s in rng.choice(starts, len(truth) // block + 1)])[:len(truth)]
        t = truth[idx]
        if 0 < t.sum() < len(t):
            diffs.append(auc_fast(t, a[idx]) - auc_fast(t, b[idx]))
    diffs = np.array(diffs)
    point = auc_fast(truth, a) - auc_fast(truth, b)
    p = float(min(1.0, 2 * min((diffs <= 0).mean(), (diffs >= 0).mean())))
    return {"auc_difference": point, "ci95": [float(np.quantile(diffs, 0.025)),
                                              float(np.quantile(diffs, 0.975))], "p": p}


def _sigmoid(x: np.ndarray) -> np.ndarray:
    return 1.0 / (1.0 + np.exp(-np.clip(x, -30, 30)))


def _calibrated(train_score, train_y, test_score) -> np.ndarray:
    """Logistic calibration of a one-number volatility score on the training
    block (Newton steps on intercept and slope). A variance forecast is not a
    probability: GARCH forecasts the MEAN variance while the regime threshold
    is a MEDIAN, so a raw "forecast above threshold" cut is biased."""
    a, b = 0.0, 1.0
    for _ in range(50):
        p = _sigmoid(a + b * train_score)
        w = np.clip(p * (1 - p), 1e-6, None)
        g = np.array([np.sum(train_y - p), np.sum((train_y - p) * train_score)])
        H = np.array([[np.sum(w), np.sum(w * train_score)],
                      [np.sum(w * train_score), np.sum(w * train_score ** 2)]])
        step = np.linalg.solve(H + 1e-9 * np.eye(2), g)
        a, b = a + step[0], b + step[1]
        if np.abs(step).max() < 1e-8:
            break
    return _sigmoid(a + b * test_score)


def har_regime(ctx: dict, train_idx, test_idx) -> np.ndarray:
    """HAR (Corsi 2009) as a regime classifier: log next-24h variance on the
    log of the last day, week and month, fitted on the training block, then
    predicted over threshold variance calibrated into a probability."""
    X = np.column_stack([np.ones(len(ctx["har_x"])), ctx["har_x"]])
    ok = train_idx[np.isfinite(ctx["fwd_logvar"][train_idx])]
    beta, *_ = np.linalg.lstsq(X[ok], ctx["fwd_logvar"][ok], rcond=None)
    score = X @ beta - ctx["log_thr"]
    return _calibrated(score[train_idx], ctx["y"][train_idx], score[test_idx])


def garch_regime(ctx: dict, train_idx, test_idx, horizon: int) -> np.ndarray:
    """GARCH(1,1) by maximum likelihood on the training block's returns, then
    filtered causally through every hour; the mean variance forecast over the
    next `horizon` hours is compared with the threshold, as for HAR."""
    from scipy.optimize import minimize

    r = ctx["ret"] * 100.0                                   # percent, for conditioning
    fit = r[train_idx]
    var0 = float(np.var(fit))

    def filter_var(params, x):
        omega, alpha, beta = params
        s2 = np.empty(len(x))
        s2[0] = var0
        for t in range(1, len(x)):
            s2[t] = omega + alpha * x[t - 1] ** 2 + beta * s2[t - 1]
        return s2

    def nll(params):
        omega, alpha, beta = params
        if omega <= 0 or alpha < 0 or beta < 0 or alpha + beta >= 0.999:
            return 1e12
        s2 = filter_var(params, fit)
        return 0.5 * float(np.sum(np.log(s2) + fit ** 2 / s2))

    res = minimize(nll, x0=[var0 * 0.05, 0.08, 0.9], method="Nelder-Mead",
                   options={"maxiter": 2000, "xatol": 1e-6, "fatol": 1e-6})
    omega, alpha, beta = res.x
    s2 = filter_var(res.x, r)
    # s2[t] is the variance of hour t given t-1; one step on gives hour t+1.
    nxt = omega + alpha * r ** 2 + beta * s2
    persist = alpha + beta
    long_run = omega / max(1e-9, 1 - persist)
    k = np.arange(horizon)
    path_mean = long_run + (nxt[:, None] - long_run) * (persist ** k)[None, :]
    score = np.log(path_mean.mean(axis=1) / 1e4) - ctx["log_thr"]
    return _calibrated(score[train_idx], ctx["y"][train_idx], score[test_idx])


def evaluate_arm(name, X, y, seq_len, horizon, folds, epochs, ctx=None) -> dict | None:
    """Purged walk-forward. The last fifth of each training block is validation."""
    usable = np.arange(seq_len - 1, len(y))
    edges = np.linspace(0, len(usable), folds + 2, dtype=int)
    results, pooled = [], []
    for f in range(folds):
        train_end, test_end = edges[f + 1], edges[f + 2]
        train_all = usable[: max(0, train_end - horizon)]
        test_idx = usable[train_end:test_end]
        if len(train_all) < 500 or len(test_idx) < 100:
            continue
        cut = int(len(train_all) * 0.8)
        train_idx, valid_idx = train_all[:cut], train_all[cut:]
        truth = y[test_idx].astype(int)
        if len(np.unique(truth)) < 2:
            continue

        vprobs = None                      # validation-block probabilities, for the cut-off
        if name == "majority":
            probs = np.full(len(test_idx), float(y[train_idx].mean()))
        elif name == "persistence":
            # Regime: "the current 24h is already above the median". Rises:
            # the label of the last fully observed window. A 0/1 rule: no cut-off.
            probs = (ctx["regime_now"][test_idx] if ctx else y[test_idx - horizon]).astype(float)
        elif name == "HAR":
            probs = har_regime(ctx, train_all, test_idx)
            vprobs = har_regime(ctx, train_idx, valid_idx)
        elif name == "GARCH(1,1)":
            probs = garch_regime(ctx, train_all, test_idx, horizon)
            vprobs = garch_regime(ctx, train_idx, valid_idx, horizon)
        elif name == "gradient boosting":
            # Baseline 3 in the methodology, brought into the same ladder on
            # identical folds rather than quoted from a separate screen. It
            # reads the latest row only, like logistic and the MLP: the point
            # of its rung is the tree ensemble, not a window.
            try:
                from sklearn.ensemble import HistGradientBoostingClassifier
            except Exception:                                      # noqa: BLE001
                return None          # unavailable on this host; the rung is skipped
            model = HistGradientBoostingClassifier(
                max_iter=250, learning_rate=0.06, max_depth=4,
                l2_regularization=1.0, random_state=42)
            model.fit(X[train_idx], y[train_idx].astype(int))
            probs = model.predict_proba(X[test_idx])[:, 1]
            vprobs = model.predict_proba(X[valid_idx])[:, 1]
        else:
            mean = X[train_idx].mean(axis=0)
            std = X[train_idx].std(axis=0)
            Xs = (X - mean) / np.clip(std, 1e-8, None)
            probs, vprobs = train_one(name, Xs, y, train_idx, valid_idx, test_idx,
                                      seq_len, epochs, 42)
        threshold = tune_threshold(y[valid_idx].astype(int), vprobs)
        results.append(score(truth, probs, threshold) | {"fold": f + 1})
        pooled.append((test_idx, truth, probs))

    if not results:
        return None
    summary = {k: float(np.mean([r[k] for r in results])) for k in results[0] if k != "fold"}
    return summary | {"folds": len(results), "per_fold": results,
                      "_pooled": tuple(np.concatenate(p) for p in zip(*pooled))}


def main() -> None:
    ap = argparse.ArgumentParser()
    ap.add_argument("--symbols", nargs="+", default=list(HORIZON))
    ap.add_argument("--folds", type=int, default=5)
    ap.add_argument("--epochs", type=int, default=15)
    ap.add_argument("--seq-len", type=int, default=48,
                    help="held fixed across architectures so the comparison is fair")
    ap.add_argument("--max-rows", type=int, default=None)
    ap.add_argument("--target", choices=["rises", "regime"], default="rises",
                    help="regime = the deployed 24h volatility-regime question, with HAR and "
                         "GARCH(1,1) rungs; rises = the original screen's target")
    ap.add_argument("--skip-features", action="store_true",
                    help="architecture ladder only, no Fear & Greed / on-chain arm")
    args = ap.parse_args()

    print(f"device: {DEVICE}   seq_len {args.seq_len} (same for every architecture)")
    fear = daily_block(DATA / "fear_greed.csv")
    chain = np.log(daily_block(DATA / "onchain.csv").clip(lower=EPS)).diff().dropna()

    # Ordered from naive to proposed, which is the order the thesis tables use.
    # Gradient boosting sits with the other single-row models; it is skipped
    # with a note where scikit-learn is unavailable rather than failing the run.
    regime = args.target == "regime"
    ladder = (["majority", "persistence"] + (["HAR", "GARCH(1,1)"] if regime else [])
              + ["logistic", "gradient boosting", "MLP"]
              + [k for k in ARCHITECTURES if k not in ("logistic", "MLP")])
    results: dict = {}

    for symbol in args.symbols:
        data, features = load(symbol, args.max_rows, args.target)
        horizon = REGIME_HORIZON if regime else HORIZON[symbol]
        y = data["y"].to_numpy(dtype=float)
        market = data[features].to_numpy(dtype=float)
        ctx = None
        if regime:
            thr = data["threshold_causal"].to_numpy(dtype=float)
            ctx = {"y": y, "regime_now": (data["vol_24"].to_numpy() > thr).astype(float),
                   "har_x": np.log(data[[f"rv{w}" for w in HAR_WINDOWS]].to_numpy(dtype=float)),
                   "fwd_logvar": np.log(data["forward_vol"].to_numpy(dtype=float) ** 2),
                   "log_thr": np.log(thr ** 2),
                   "ret": np.nan_to_num(data["log_return"].to_numpy(dtype=float))}
        label = "volatility regime (next 24h above 90-day median)" if regime else "volatility rises"
        print(f"\n{'=' * 96}\n{symbol}   {len(data):,} rows   horizon {horizon}h"
              f"   target: {label}\n{'=' * 96}")
        print(f"{'model':<20}{'ROC-AUC':>9}{'acc tuned':>10}{'MCC tuned':>10}{'cut-off':>8}"
              f"{'acc @0.5':>9}{'MCC @0.5':>9}{'majority':>9}{'folds':>6}")
        results[symbol] = {"architecture": {}, "features": {}, "significance": {}}
        pooled = {}

        for name in ladder:
            got = evaluate_arm(name, market, y, args.seq_len, horizon,
                               args.folds, args.epochs, ctx)
            if got is None:
                print(f"{name:<20}  skipped (unavailable on this host or no usable fold)")
                continue
            pooled[name] = got.pop("_pooled")
            results[symbol]["architecture"][name] = got
            print(f"{name:<20}{got['roc_auc']:>9.4f}{got['accuracy']*100:>9.2f}%"
                  f"{got['mcc']:>+10.3f}{got['threshold']:>8.2f}"
                  f"{got['accuracy_at_0.5']*100:>8.2f}%{got['mcc_at_0.5']:>+9.3f}"
                  f"{got['majority']*100:>8.2f}%{got['folds']:>6}")

        # ---- is the proposed model's AUC significantly different from each rung's?
        proposed = "VAE+Transformer"
        if proposed in pooled:
            rng = np.random.default_rng(0)
            idx_p, truth_p, prob_p = pooled[proposed]
            print(f"\n{proposed} minus each rung, ROC-AUC on pooled test hours "
                  f"(moving-block bootstrap, one-week blocks):")
            for other, (idx_o, _, prob_o) in pooled.items():
                if other == proposed or other == "majority":
                    continue
                common, ia, ib = np.intersect1d(idx_p, idx_o, return_indices=True)
                got = auc_difference(truth_p[ia].astype(int), prob_p[ia], prob_o[ib], rng)
                results[symbol]["significance"][other] = got
                print(f"  vs {other:<20} {got['auc_difference']:>+8.4f}  "
                      f"95% CI [{got['ci95'][0]:+.4f}, {got['ci95'][1]:+.4f}]  p = {got['p']:.3f}")

        if args.skip_features:
            continue
        # ---- historic versus sentiment, on the deployed architecture
        arms = {"market only (33)": market,
                "+ fear & greed": np.column_stack([market, join_hourly(data, fear)])}
        if symbol == "BTCUSDT":
            arms["+ on-chain"] = np.column_stack([market, join_hourly(data, chain)])
            arms["+ both"] = np.column_stack(
                [market, join_hourly(data, fear), join_hourly(data, chain)])
        print(f"\n{'feature set':<20}{'ROC-AUC':>9}{'accuracy':>10}{'margin':>9}"
              f"{'MCC':>8}{'vs market':>11}")
        reference = None
        for label, matrix in arms.items():
            finite = np.isfinite(matrix).all(axis=1)
            if finite.sum() < 3000:
                print(f"{label:<20}  too few usable rows")
                continue
            got = evaluate_arm("VAE+Transformer", np.nan_to_num(matrix), y,
                               args.seq_len, horizon, args.folds, args.epochs)
            if got is None:
                continue
            got.pop("_pooled", None)
            if reference is None:
                reference = got["roc_auc"]
            got["delta_auc"] = got["roc_auc"] - reference
            results[symbol]["features"][label] = got
            print(f"{label:<20}{got['roc_auc']:>9.4f}{got['accuracy']*100:>9.2f}%"
                  f"{got['margin']*100:>+9.2f}{got['mcc']:>+8.3f}"
                  f"{got['delta_auc']:>+11.4f}")

    out = ROOT / "outputs" / ("ablation_study_regime.json" if regime else "ablation_study.json")
    out.write_text(json.dumps(results, indent=2), encoding="utf-8")
    print(f"\nwritten: {out.relative_to(ROOT)}")


if __name__ == "__main__":
    main()