1"""Evaluation metrics, paired block-bootstrap intervals, Diebold-Mariano / Giacomini-White tests, Holm correction,2calibration, and two-way clustered panel regressions."""3from __future__ import annotations45import numpy as np6import pandas as pd7import polars as pl8from scipy import stats9from sklearn.metrics import average_precision_score, brier_score_loss, log_loss, roc_auc_score101112# --------------------------------------------------------------------- metrics13def binary_metrics(y: np.ndarray, p: np.ndarray) -> dict:14 p = np.clip(p, 1e-6, 1 - 1e-6)15 out = {"n": int(len(y)), "base_rate": float(y.mean()), "logloss": float(log_loss(y, p)), "brier": float(brier_score_loss(y, p))}16 if 0 < y.mean() < 1:17 out["auc"] = float(roc_auc_score(y, p))18 out["pr_auc"] = float(average_precision_score(y, p))19 # calibration intercept / slope (logistic recalibration of the logit)20 z = np.log(p / (1 - p))21 try:22 from sklearn.linear_model import LogisticRegression23 m = LogisticRegression(C=1e6, solver="lbfgs", max_iter=1000).fit(z.reshape(-1, 1), y)24 out["cal_slope"], out["cal_intercept"] = float(m.coef_[0][0]), float(m.intercept_[0])25 except Exception:26 pass27 return out282930def continuous_metrics(y: np.ndarray, yhat: np.ndarray, y_train_mean: float) -> dict:31 e = y - yhat32 ss = np.sum((y - y_train_mean) ** 2)33 return {"n": int(len(y)), "mae": float(np.mean(np.abs(e))), "rmse": float(np.sqrt(np.mean(e ** 2))),34 "oos_r2": float(1 - np.sum(e ** 2) / ss) if ss > 0 else float("nan")}353637def pointwise_logloss(y: np.ndarray, p: np.ndarray) -> np.ndarray:38 p = np.clip(p, 1e-6, 1 - 1e-6)39 return -(y * np.log(p) + (1 - y) * np.log(1 - p))404142def reliability_table(y: np.ndarray, p: np.ndarray, bins: int = 10) -> pd.DataFrame:43 q = np.quantile(p, np.linspace(0, 1, bins + 1))44 q[0], q[-1] = -np.inf, np.inf45 b = np.digitize(p, q[1:-1])46 rows = []47 for k in range(bins):48 m = b == k49 if m.sum():50 rows.append({"decile": k + 1, "n": int(m.sum()), "mean_pred": float(p[m].mean()), "obs_rate": float(y[m].mean())})51 return pd.DataFrame(rows)525354# ------------------------------------------------------------ paired inference55def quarter_blocks(dates: np.ndarray) -> np.ndarray:56 d = pd.to_datetime(dates)57 return (d.year * 4 + (d.month - 1) // 3).to_numpy()585960def paired_block_bootstrap(d: np.ndarray, blocks: np.ndarray, n_boot: int = 1000, seed: int = 20261102) -> dict:61 """Mean of a per-observation loss differential with a cluster (block) bootstrap over calendar quarters.62 Returns the mean, a percentile 95 % CI and a two-sided bootstrap p-value for H0: mean = 0."""63 rng = np.random.default_rng(seed)64 ub = np.unique(blocks)65 # per-block sums and counts -> resample blocks with replacement (equal-weight by observation)66 sums = np.array([d[blocks == b].sum() for b in ub]); cnts = np.array([(blocks == b).sum() for b in ub])67 mean = float(d.mean())68 draws = np.empty(n_boot)69 for i in range(n_boot):70 idx = rng.integers(0, len(ub), len(ub))71 draws[i] = sums[idx].sum() / cnts[idx].sum()72 lo, hi = np.percentile(draws, [2.5, 97.5])73 p = float(2 * min((draws - mean <= -abs(mean)).mean(), (draws - mean >= abs(mean)).mean()))74 return {"mean": mean, "ci_lo": float(lo), "ci_hi": float(hi), "p_boot": p, "n_blocks": int(len(ub)), "n": int(len(d))}757677def dm_gw_test(d: np.ndarray, blocks: np.ndarray) -> dict:78 """Unconditional Giacomini-White (= Diebold-Mariano with HAC) on the per-quarter mean loss differential.79 HAC: Newey-West with lag = floor(4 (T/100)^(2/9))."""80 ub = np.unique(blocks)81 dq = np.array([d[blocks == b].mean() for b in ub])82 T = len(dq)83 if T < 4:84 return {"t": float("nan"), "p": float("nan"), "T": T}85 m = dq.mean(); u = dq - m86 L = int(np.floor(4 * (T / 100) ** (2 / 9)))87 s2 = u @ u / T88 for l in range(1, L + 1):89 w = 1 - l / (L + 1)90 s2 += 2 * w * (u[l:] @ u[:-l]) / T91 se = np.sqrt(max(s2, 1e-18) / T)92 t = m / se93 return {"t": float(t), "p": float(2 * (1 - stats.t.cdf(abs(t), df=T - 1))), "T": int(T), "mean_quarterly": float(m)}949596def holm(pvals: dict[str, float]) -> dict[str, float]:97 items = sorted(pvals.items(), key=lambda kv: kv[1])98 m = len(items); out = {}; running = 0.099 for i, (k, p) in enumerate(items):100 adj = min(1.0, (m - i) * p); running = max(running, adj); out[k] = running101 return out102103104# ------------------------------------------------------------- panel regression105def demean_two_way(df: pd.DataFrame, cols: list[str], fe1: str, fe2: str, iters: int = 25) -> pd.DataFrame:106 """Alternating projections for two fixed effects (firm, sector x quarter)."""107 x = df[cols].astype(float).copy()108 for _ in range(iters):109 x = x - x.groupby(df[fe1]).transform("mean")110 x = x - x.groupby(df[fe2]).transform("mean")111 return x112113114def two_way_cluster_ols(df: pd.DataFrame, y: str, xs: list[str], fe1: str, fe2: str, cl1: str, cl2: str) -> pd.DataFrame:115 """OLS with two-way FE (within transformation) and Cameron-Gelbach-Miller two-way clustered SEs."""116 d = df.dropna(subset=[y] + xs + [fe1, fe2, cl1, cl2]).copy()117 z = demean_two_way(d, [y] + xs, fe1, fe2)118 X = z[xs].to_numpy(); Y = z[y].to_numpy()119 XtX_inv = np.linalg.pinv(X.T @ X)120 beta = XtX_inv @ X.T @ Y121 e = Y - X @ beta122 def meat(groups):123 S = np.zeros((X.shape[1], X.shape[1]))124 for _, idx in pd.Series(np.arange(len(groups))).groupby(groups.to_numpy()).groups.items():125 idx = np.asarray(list(idx)); s = X[idx].T @ e[idx]; S += np.outer(s, s)126 return S127 g12 = (d[cl1].astype(str) + "|" + d[cl2].astype(str))128 V = XtX_inv @ (meat(d[cl1]) + meat(d[cl2]) - meat(g12)) @ XtX_inv129 se = np.sqrt(np.clip(np.diag(V), 0, None))130 t = beta / se131 G = min(d[cl1].nunique(), d[cl2].nunique())132 p = 2 * (1 - stats.t.cdf(np.abs(t), df=max(G - 1, 1)))133 return pd.DataFrame({"var": xs, "coef": beta, "se": se, "t": t, "p": p, "n": len(d), "n_firms": d[cl1].nunique(),134 "n_quarters": d[cl2].nunique()})135