1"""r3k_inference.py — Stage D (FIRC-0002, package L4, clusters E / M / N) inference helpers, imported by the r3k_15+ scripts.2Everything here is self-contained (numpy / scipy / statsmodels): no remembered critical values — the fixed-b and sup-Wald3reference distributions are simulated with a fixed seed and the simulation parameters are returned with the statistic.4 cgm_cov(res, g1, g2) two-way clustered covariance (Cameron-Gelbach-Miller 2011) for a fitted statsmodels OLS5 ols_cluster(y, X, groups) OLS with one-way cluster-robust SE (small-sample G/(G-1) x (N-1)/(N-k) adjustment)6 wild_cluster_boot(y, X, g, j, ...) restricted wild-cluster bootstrap p-value and bootstrap t draws (Rademacher / Webb)7 nw_t(x, lags) Newey-West (Bartlett, `lags`) t of a mean; plain_t(x) its i.i.d. counterpart; ci958 fixed_b_cv(T, bw, ...) Kiefer-Vogelsang (2005) fixed-b critical value for the Bartlett HAC t, simulated9 block_boot_ci(x, B, seed) percentile CI of a mean by resampling the series elements (cohort block bootstrap)10 mde(se, power=0.8) minimum detectable effect at 5 % two-sided11 holm(p), bh(p), romano_wolf(t_obs, t_boot) multiple-testing adjustments (RW: step-down max-t on shared bootstrap draws)12 deflated_sharpe(sr, T, skew, kurt, sr_trials) Bailey-Lopez de Prado (2014) DSR probability13 permutation_three_cohorts(spreads, post_years) exact permutation p over all C(n,3) subsets and over contiguous blocks14 sup_wald(x, trim, nsim, seed) Andrews (1993) sup-Wald mean-shift test with a simulated finite-sample p-value15"""16import itertools, numpy as np, pandas as pd17from scipy import stats18from statsmodels.stats.sandwich_covariance import cov_hac19import statsmodels.api as sm20RNG_SEED = 202611052122def cgm_cov(res, g1, g2):23 v1 = res.get_robustcov_results(cov_type="cluster", groups=np.asarray(g1)).cov_params()24 v2 = res.get_robustcov_results(cov_type="cluster", groups=np.asarray(g2)).cov_params()25 g12 = pd.factorize(pd.Series(list(zip(np.asarray(g1), np.asarray(g2)))))[0]26 v12 = res.get_robustcov_results(cov_type="cluster", groups=g12).cov_params()27 return np.asarray(v1) + np.asarray(v2) - np.asarray(v12)2829def se_from_cov(V):30 """Standard errors from a covariance matrix; the CGM two-way estimator V1 + V2 - V12 is not guaranteed positive semi-definite31 with few clusters (3 post cohorts), so a negative diagonal element yields NaN (printed as '--'), never a complex number."""32 d = np.diag(np.asarray(V, float)); return np.where(d > 0, np.sqrt(np.where(d > 0, d, 1.0)), np.nan)3334def ols_cluster(y, X, groups):35 y = np.asarray(y, float); X = np.asarray(X, float); g = pd.factorize(np.asarray(groups))[0]; n, k = X.shape; G = g.max() + 136 XtX_inv = np.linalg.pinv(X.T @ X); b = XtX_inv @ X.T @ y; e = y - X @ b37 meat = np.zeros((k, k))38 for gi in range(G):39 m = g == gi; s = X[m].T @ e[m]; meat += np.outer(s, s)40 adj = (G / (G - 1)) * ((n - 1) / (n - k)); V = adj * XtX_inv @ meat @ XtX_inv41 return b, np.sqrt(np.diag(V)), e, XtX_inv, g4243def wild_cluster_boot(y, X, groups, j, B=9999, weights="rademacher", seed=RNG_SEED):44 """Restricted (null-imposed) wild-cluster bootstrap for H0: beta_j = 0 (Cameron-Gelbach-Miller 2008; Webb 2023 six-point45 weights for few clusters). Returns dict(p, t_obs, t_boot (B,), G)."""46 y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape; g = pd.factorize(np.asarray(groups))[0]; G = g.max() + 147 b, se, e, XtX_inv, _ = ols_cluster(y, X, g); t_obs = b[j] / se[j]48 Xr = np.delete(X, j, axis=1); br = np.linalg.pinv(Xr.T @ Xr) @ Xr.T @ y; er = y - Xr @ br; yr = Xr @ br49 rng = np.random.default_rng(seed); adj = (G / (G - 1)) * ((n - 1) / (n - k))50 blocks = [np.where(g == gi)[0] for gi in range(G)]; Xg = [X[ix] for ix in blocks]51 if weights == "webb":52 pts = np.array([-np.sqrt(1.5), -1, -np.sqrt(0.5), np.sqrt(0.5), 1, np.sqrt(1.5)]); W = rng.choice(pts, size=(B, G))53 else:54 W = rng.choice([-1.0, 1.0], size=(B, G))55 XtX_inv_Xt = XtX_inv @ X.T; t_boot = np.empty(B)56 for bi in range(B):57 ys = yr + er * W[bi][g]; bs = XtX_inv_Xt @ ys; es = ys - X @ bs; meat = np.zeros((k, k))58 for gi in range(G):59 s = Xg[gi].T @ es[blocks[gi]]; meat += np.outer(s, s)60 V = adj * XtX_inv @ meat @ XtX_inv; t_boot[bi] = bs[j] / np.sqrt(V[j, j])61 p = float((np.abs(t_boot) >= abs(t_obs)).mean()); return {"p": round(p, 4), "t_obs": round(float(t_obs), 2), "t_boot": t_boot, "G": int(G), "B": B, "weights": weights}6263def nw_t(x, lags=3):64 x = pd.Series(x).dropna().values65 if len(x) < 3: return (float(x.mean()) if len(x) else np.nan), np.nan66 r = sm.OLS(x, np.ones((len(x), 1))).fit(); return float(x.mean()), float(x.mean() / np.sqrt(cov_hac(r, nlags=lags)[0, 0]))67def plain_t(x):68 x = pd.Series(x).dropna().values69 if len(x) < 2: return np.nan70 return float(x.mean() / (x.std(ddof=1) / np.sqrt(len(x))))71def ci95(x):72 x = pd.Series(x).dropna().values73 if len(x) < 2: return [None, None]74 se = x.std(ddof=1) / np.sqrt(len(x)); q = stats.t.ppf(0.975, len(x) - 1); return [float(x.mean() - q * se), float(x.mean() + q * se)]75def nw_ci95(x, lags=3):76 x = pd.Series(x).dropna().values77 if len(x) < 3: return [None, None]78 r = sm.OLS(x, np.ones((len(x), 1))).fit(); se = np.sqrt(cov_hac(r, nlags=lags)[0, 0]); return [float(x.mean() - 1.96 * se), float(x.mean() + 1.96 * se)]7980_FB_CACHE = {}81def fixed_b_cv(T, bw, alpha=0.05, nsim=20000, steps=1000, seed=RNG_SEED):82 """Kiefer-Vogelsang (2005) fixed-b limit of the Bartlett-kernel HAC t-statistic of a mean, b = (bw+1)/T as in statsmodels83 cov_hac(nlags=bw) (Bartlett weights 1 - l/(bw+1)). Limit: W(1)/sqrt(Q_b), Q_b = (2/b) int_0^1 B(r)^2 dr - (2/b) int_0^{1-b} B(r)B(r+b) dr,84 B = Brownian bridge. Two-sided critical value at `alpha`, simulated once per (T, bw) with a fixed seed."""85 b = (bw + 1) / T; key = (round(b, 4), alpha)86 if key in _FB_CACHE: return _FB_CACHE[key]87 rng = np.random.default_rng(seed); dt = 1.0 / steps; cvs = []88 for chunk in range(nsim // 2000):89 dW = rng.normal(0, np.sqrt(dt), size=(2000, steps)); W = np.cumsum(dW, axis=1); r = np.arange(1, steps + 1) / steps90 Bb = W - r * W[:, -1:]; lag = int(round(b * steps)); Q = (2 / b) * (Bb ** 2).sum(1) * dt - (2 / b) * (Bb[:, :steps - lag] * Bb[:, lag:]).sum(1) * dt91 cvs.append(np.abs(W[:, -1]) / np.sqrt(np.maximum(Q, 1e-12)))92 cv = float(np.quantile(np.concatenate(cvs), 1 - alpha)); _FB_CACHE[key] = {"cv": round(cv, 3), "b": round(b, 4), "nsim": nsim, "steps": steps}; return _FB_CACHE[key]9394def block_boot_ci(x, B=9999, seed=RNG_SEED):95 x = pd.Series(x).dropna().values; rng = np.random.default_rng(seed)96 if len(x) < 2: return [None, None]97 m = rng.choice(x, size=(B, len(x)), replace=True).mean(1); return [float(np.quantile(m, 0.025)), float(np.quantile(m, 0.975))]98def mde(se, power=0.8, alpha=0.05): return float((stats.norm.ppf(1 - alpha / 2) + stats.norm.ppf(power)) * se)99100def holm(p):101 p = np.asarray(p, float); m = len(p); o = np.argsort(p); adj = np.empty(m); run = 0.0102 for i, idx in enumerate(o): run = max(run, (m - i) * p[idx]); adj[idx] = min(1.0, run)103 return adj104def bh(p):105 p = np.asarray(p, float); m = len(p); o = np.argsort(p)[::-1]; adj = np.empty(m); run = 1.0106 for i, idx in enumerate(o): rank = m - i; run = min(run, m / rank * p[idx]); adj[idx] = min(1.0, run)107 return adj108def romano_wolf(t_obs, t_boot):109 """Step-down max-|t| adjustment (Romano-Wolf 2005) from a (B, m) matrix of bootstrap t's computed under the null with the SAME110 draws for all m hypotheses (so cross-sectional dependence is preserved). Returns adjusted p-values."""111 t_obs = np.abs(np.asarray(t_obs, float)); tb = np.abs(np.asarray(t_boot, float)); m = len(t_obs); o = np.argsort(-t_obs); adj = np.empty(m); prev = 0.0112 for i, idx in enumerate(o):113 rest = o[i:]; maxt = tb[:, rest].max(1); p = float((maxt >= t_obs[idx]).mean()); prev = max(prev, p); adj[idx] = prev114 return adj115116def deflated_sharpe(sr_ann, T, skew, kurt, sr_trials_ann):117 """Bailey & Lopez de Prado (2014). Inputs: annualised SR of the candidate, number of monthly observations T, skewness and118 (non-excess) kurtosis of its monthly returns, the annualised SR of every strategy tried (list). Monthly SR = SR_ann / sqrt(12)."""119 sr = sr_ann / np.sqrt(12); trials = np.asarray(sr_trials_ann, float) / np.sqrt(12); N = len(trials); gamma = 0.5772156649120 v = float(trials.var(ddof=1)) if N > 1 else 0.0121 sr0 = np.sqrt(v) * ((1 - gamma) * stats.norm.ppf(1 - 1 / N) + gamma * stats.norm.ppf(1 - 1 / (N * np.e))) if N > 1 else 0.0122 denom = np.sqrt(max(1 - skew * sr + (kurt - 1) / 4 * sr ** 2, 1e-9)); z = (sr - sr0) * np.sqrt(T - 1) / denom123 return {"dsr_prob": round(float(stats.norm.cdf(z)), 4), "sr0_ann": round(float(sr0 * np.sqrt(12)), 3), "n_trials": int(N), "psr_z": round(float(z), 3)}124125def permutation_three_cohorts(spreads, post_years):126 """spreads: dict year -> spread; post_years: the observed post block. Exact permutation over all C(n, k) subsets of k cohorts127 (statistic: mean(post) - mean(rest)) and over the n-k+1 contiguous blocks of length k."""128 yrs = sorted(spreads); v = np.array([spreads[y] for y in yrs]); k = len(post_years); obs = v[[yrs.index(y) for y in post_years]].mean() - v[[i for i, y in enumerate(yrs) if y not in post_years]].mean()129 diffs = []130 for comb in itertools.combinations(range(len(yrs)), k):131 m = np.zeros(len(yrs), bool); m[list(comb)] = True; diffs.append(v[m].mean() - v[~m].mean())132 diffs = np.array(diffs); contig = []133 for s in range(len(yrs) - k + 1):134 m = np.zeros(len(yrs), bool); m[s:s + k] = True; contig.append((yrs[s], float(v[m].mean() - v[~m].mean())))135 cvals = np.array([c[1] for c in contig])136 return {"observed_diff": float(obs), "n_subsets": int(len(diffs)), "p_two_sided": float((np.abs(diffs) >= abs(obs) - 1e-12).mean()), "p_one_sided": float((diffs >= obs - 1e-12).mean()), "rank_among_subsets": int((diffs > obs + 1e-12).sum() + 1),137 "contiguous_blocks": {str(y): round(d, 4) for y, d in contig}, "p_contiguous_one_sided": float((cvals >= obs - 1e-12).mean()), "rank_among_contiguous": int((cvals > obs + 1e-12).sum() + 1)}138139def sup_wald(x, trim=0.15, nsim=20000, seed=RNG_SEED):140 """Andrews (1993) sup-Wald test for a single mean shift at an unknown date in a short series, trimming `trim` at each end;141 the Wald statistic uses the full-sample residual variance (homoskedastic, no HAC: T is 15). The finite-sample null distribution142 is simulated (i.i.d. normal series of the same length, same trimming), so the p-value needs no tabulated critical value."""143 x = np.asarray(pd.Series(x).dropna(), float); T = len(x); lo, hi = max(2, int(np.floor(trim * T))), T - max(2, int(np.floor(trim * T)))144 def sw(v):145 best, bk = -1.0, None146 for k in range(lo, hi + 1):147 a, b = v[:k], v[k:]; s2 = (((a - a.mean()) ** 2).sum() + ((b - b.mean()) ** 2).sum()) / (T - 2); w = (a.mean() - b.mean()) ** 2 / (s2 * (1 / k + 1 / (T - k)))148 if w > best: best, bk = w, k149 return best, bk150 obs, k = sw(x); rng = np.random.default_rng(seed); sims = np.array([sw(rng.normal(size=T))[0] for _ in range(nsim)])151 return {"sup_wald": round(float(obs), 3), "break_after_index": int(k), "p_simulated": round(float((sims >= obs).mean()), 4), "T": T, "trim": trim, "nsim": nsim, "cv05_simulated": round(float(np.quantile(sims, 0.95)), 3), "cv10_simulated": round(float(np.quantile(sims, 0.90)), 3)}152