| """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 |
|
|
| DATA = ROOT / "outputs" / "sentiment_replication" |
| EPS = 1e-12 |
| DEVICE = torch.device("cuda" if torch.cuda.is_available() else "cpu") |
|
|
| |
| HORIZON = {"BTCUSDT": 12, "ETHUSDT": 12, "ADAUSDT": 1} |
|
|
|
|
| |
| 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): |
| 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): |
| 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, |
| } |
|
|
|
|
| |
| 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 |
|
|
|
|
| |
| 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 |
| 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) |
| |
| 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 |
| if name == "majority": |
| probs = np.full(len(test_idx), float(y[train_idx].mean())) |
| elif name == "persistence": |
| |
| |
| 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": |
| |
| |
| |
| |
| try: |
| from sklearn.ensemble import HistGradientBoostingClassifier |
| except Exception: |
| return None |
| 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() |
|
|
| |
| |
| |
| 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}") |
|
|
| |
| 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 |
| |
| 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() |
|
|
