1"""Dependence-aware inference: Diebold–Mariano with HAC (overlap-aware lag), cluster bootstrap by trading day2(pooling all assets of a day into one cluster), Holm and Benjamini–Hochberg controls."""3from __future__ import annotations45import numpy as np6from scipy import stats789def hac_variance(d: np.ndarray, lag: int) -> float:10 d = d - d.mean()11 n = len(d)12 v = (d * d).sum() / n13 for k in range(1, lag + 1):14 w = 1 - k / (lag + 1)15 v += 2 * w * (d[k:] * d[:-k]).sum() / n16 return float(max(v, 1e-30))171819def diebold_mariano(loss_a: np.ndarray, loss_b: np.ndarray, order: np.ndarray, lag: int) -> dict:20 """DM on the loss differential ordered chronologically (`order` sorts the observations), HAC lag = overlap."""21 d = (loss_a - loss_b)[np.argsort(order, kind="stable")]22 n = len(d)23 v = hac_variance(d, lag)24 t = d.mean() / np.sqrt(v / n)25 return {"mean_diff": float(d.mean()), "t": float(t), "p": float(2 * (1 - stats.norm.cdf(abs(t)))), "n": int(n),26 "lag": int(lag)}272829def cluster_bootstrap_mean(x: np.ndarray, clusters: np.ndarray, n_boot: int = 2000, seed: int = 0,30 alpha: float = 0.05) -> dict:31 """Percentile CI of mean(x) resampling whole clusters (trading days) with replacement."""32 rng = np.random.default_rng(seed)33 uc, inv = np.unique(clusters, return_inverse=True)34 sums = np.bincount(inv, weights=x, minlength=len(uc)); cnt = np.bincount(inv, minlength=len(uc)).astype(float)35 C = len(uc)36 boots = np.empty(n_boot)37 for b in range(n_boot):38 idx = rng.integers(0, C, C)39 boots[b] = sums[idx].sum() / cnt[idx].sum()40 return {"mean": float(x.mean()), "lo": float(np.quantile(boots, alpha / 2)), "hi": float(np.quantile(boots, 1 - alpha / 2)),41 "se": float(boots.std(ddof=1)), "n_clusters": int(C), "n": int(len(x))}424344def holm(pvals: list[float]) -> list[float]:45 p = np.asarray(pvals, dtype=float); m = len(p)46 order = np.argsort(p); adj = np.empty(m)47 running = 0.048 for rank, i in enumerate(order):49 running = max(running, (m - rank) * p[i])50 adj[i] = min(1.0, running)51 return adj.tolist()525354def benjamini_hochberg(pvals: list[float]) -> list[float]:55 p = np.asarray(pvals, dtype=float); m = len(p)56 order = np.argsort(p)[::-1]; adj = np.empty(m); running = 1.057 for rank, i in zip(range(m, 0, -1), order):58 running = min(running, m / rank * p[i]); adj[i] = running59 return adj.tolist()60