1"""Proper scoring rules and secondary metrics. Conventions: y = realised cumulative log return at horizon H (or the2per-minute return), forecasts are either parametric Student-t (closed form) or sample arrays [n, N]."""3from __future__ import annotations45import numpy as np6from scipy import special, stats789# ----------------------------------------------------------------------------- closed-form Student-t CRPS10def crps_student_t(y: np.ndarray, mu: np.ndarray, scale: np.ndarray, nu: np.ndarray) -> np.ndarray:11 """CRPS of a location-scale Student-t (Jordan, Krüger & Lerch 2019, scoringRules), valid for nu > 1."""12 y, mu, scale, nu = (np.asarray(a, dtype=np.float64) for a in (y, mu, scale, nu))13 z = (y - mu) / scale14 F = stats.t.cdf(z, nu); f = stats.t.pdf(z, nu)15 b1 = special.beta(0.5, nu - 0.5); b2 = special.beta(0.5, nu / 2) ** 216 term = z * (2 * F - 1) + 2 * f * (nu + z * z) / (nu - 1) - 2 * np.sqrt(nu) * b1 / ((nu - 1) * b2)17 return scale * term181920def crps_gaussian(y: np.ndarray, mu: np.ndarray, sd: np.ndarray) -> np.ndarray:21 z = (y - mu) / sd22 return sd * (z * (2 * stats.norm.cdf(z) - 1) + 2 * stats.norm.pdf(z) - 1 / np.sqrt(np.pi))232425# ----------------------------------------------------------------------------- sample based26def crps_samples(y: np.ndarray, s: np.ndarray) -> np.ndarray:27 """Energy-form estimator from samples s [n, N] for observations y [N] (unbiased pairwise term)."""28 n = s.shape[0]29 t1 = np.abs(s - y[None, :]).mean(0)30 ss = np.sort(s, axis=0)31 # E|X - X'| via order statistics: (2/(n(n-1))) sum_i (2i - n - 1) x_(i)32 i = np.arange(1, n + 1)[:, None]33 t2 = (2.0 / (n * (n - 1))) * ((2 * i - n - 1) * ss).sum(0)34 return t1 - 0.5 * t2353637def quantiles_from_samples(s: np.ndarray, levels: np.ndarray) -> np.ndarray:38 return np.quantile(s, levels, axis=0) # [Q, N]394041def quantile_loss(y: np.ndarray, q: np.ndarray, levels: np.ndarray) -> np.ndarray:42 """Mean pinball loss over levels; q [Q, N]."""43 d = y[None, :] - q44 return np.maximum(levels[:, None] * d, (levels[:, None] - 1) * d).mean(0)454647def interval_score(y: np.ndarray, lo: np.ndarray, hi: np.ndarray, alpha: float) -> np.ndarray:48 return (hi - lo) + 2 / alpha * (lo - y) * (y < lo) + 2 / alpha * (y - hi) * (y > hi)495051def weighted_interval_score(y: np.ndarray, s: np.ndarray, alphas=(0.02, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)) -> np.ndarray:52 """WIS (Bracher et al. 2021) with K central intervals + median."""53 med = np.median(s, axis=0)54 K = len(alphas)55 tot = 0.5 * np.abs(y - med)56 for a in alphas:57 lo, hi = np.quantile(s, [a / 2, 1 - a / 2], axis=0)58 tot = tot + a / 2 * interval_score(y, lo, hi, a)59 return tot / (K + 0.5)606162def coverage(y: np.ndarray, s: np.ndarray, level: float) -> tuple[float, float]:63 """Empirical coverage and mean width of the central `level` interval."""64 a = 1 - level65 lo, hi = np.quantile(s, [a / 2, 1 - a / 2], axis=0)66 return float(((y >= lo) & (y <= hi)).mean()), float((hi - lo).mean())676869def pit_samples(y: np.ndarray, s: np.ndarray) -> np.ndarray:70 """Randomised PIT from samples (ties broken uniformly)."""71 n = s.shape[0]72 below = (s < y[None, :]).sum(0); equal = (s == y[None, :]).sum(0)73 u = np.random.default_rng(0).random(len(y))74 return (below + u * (equal + 1)) / (n + 1)757677def log_score_student_t(y, mu, scale, nu) -> np.ndarray:78 return -stats.t.logpdf((y - mu) / scale, nu) + np.log(scale)798081def log_score_kde(y: np.ndarray, s: np.ndarray) -> np.ndarray:82 """Negative log density from a Gaussian KDE of the samples (comparable only across sample-based forecasts)."""83 n = s.shape[0]84 sd = s.std(0, ddof=1) + 1e-1285 h = 1.06 * sd * n ** (-0.2)86 dens = stats.norm.pdf((y[None, :] - s) / h[None, :]).mean(0) / h87 return -np.log(dens + 1e-300)888990# ----------------------------------------------------------------------------- secondary91def qlike(rv: np.ndarray, var_hat: np.ndarray) -> np.ndarray:92 v = np.maximum(var_hat, 1e-14)93 return np.log(v) + rv / v949596def brier(event: np.ndarray, p: np.ndarray) -> np.ndarray:97 return (p - event.astype(float)) ** 29899100def directional_accuracy(y: np.ndarray, mu: np.ndarray) -> float:101 nz = y != 0102 return float((np.sign(mu[nz]) == np.sign(y[nz])).mean()) if nz.any() else float("nan")103104105# ----------------------------------------------------------------------------- multivariate (paths)106def energy_score(Y: np.ndarray, S: np.ndarray) -> np.ndarray:107 """Y [N, H] observations, S [n, N, H] samples. ES = E||S - Y|| - 0.5 E||S - S'|| (pairs via a random permutation)."""108 n = S.shape[0]109 t1 = np.linalg.norm(S - Y[None], axis=-1).mean(0)110 perm = np.random.default_rng(1).permutation(n)111 t2 = np.linalg.norm(S - S[perm], axis=-1).mean(0)112 return t1 - 0.5 * t2113114115def variogram_score(Y: np.ndarray, S: np.ndarray, p: float = 0.5, max_pairs: int = 400) -> np.ndarray:116 """Variogram score of order p (Scheuerer & Hamill 2015) on a random subset of index pairs, equal weights."""117 N, H = Y.shape118 rng = np.random.default_rng(2)119 i = rng.integers(0, H, max_pairs); j = rng.integers(0, H, max_pairs)120 keep = i != j; i, j = i[keep], j[keep]121 obs = np.abs(Y[:, i] - Y[:, j]) ** p122 sim = (np.abs(S[:, :, i] - S[:, :, j]) ** p).mean(0)123 return ((obs - sim) ** 2).mean(-1)124125126def path_stats(S: np.ndarray) -> dict[str, np.ndarray]:127 """From sample paths of per-minute returns S [n, N, H]: cumulative path, max drawdown, realised vol."""128 cum = np.cumsum(S, axis=-1)129 peak = np.maximum.accumulate(np.concatenate([np.zeros_like(cum[..., :1]), cum], -1), axis=-1)[..., 1:]130 mdd = (peak - cum).max(-1)131 rv = (S ** 2).sum(-1)132 return {"cum": cum, "mdd": mdd, "rv": rv}133