Repository Public Research preview
1# Author: Simon-Pierre Boucher — contact@spboucher.ai2"""Inference for out-of-sample accuracy comparisons.34Two tools, both operating on the stored test-set predictions:56* a *cluster bootstrap* (resampling municipalities with replacement) for7 the MdAPE of each model and for the MdAPE difference between two8 models evaluated on the same test sales;9* a Diebold–Mariano / Giacomini–White-type test of equal predictive10 accuracy: the mean of the per-sale loss differential d_i = L(e_i^A) -11 L(e_i^B) is tested against zero with standard errors clustered by12 municipality (the cross-sectional analogue of the HAC-robust DM13 statistic). Losses: absolute percentage error (the MAPE loss) and the14 squared log error (the RMSE_ln loss).15"""16from __future__ import annotations1718import numpy as np19import pandas as pd2021from . import config222324def ape(y: np.ndarray, p: np.ndarray) -> np.ndarray:25 return np.abs(p - y) / y262728def sq_log(y: np.ndarray, p: np.ndarray) -> np.ndarray:29 return (np.log(y) - np.log(p)) ** 2303132def cluster_bootstrap_mdape(P: pd.DataFrame, models: list[str],33 cluster: str = "muni", B: int = 500,34 ref: str | None = None,35 seed: int = config.SEED) -> pd.DataFrame:36 """Percentile 95% intervals for MdAPE (and MdAPE minus the reference37 model's MdAPE) from B cluster-bootstrap replicates.3839 Clusters are municipalities; each replicate draws clusters with40 replacement (so a drawn municipality contributes all its test sales,41 possibly several times) and recomputes the median APE.42 """43 y = P["price"].to_numpy(float)44 A = {m: ape(y, P[m].to_numpy(float)) for m in models}45 codes, uniq = pd.factorize(P[cluster])46 G = len(uniq)47 rng = np.random.default_rng(seed)48 # pre-compute row indices per cluster49 order = np.argsort(codes, kind="stable")50 bounds = np.searchsorted(codes[order], np.arange(G + 1))51 rows_by_cluster = [order[bounds[g]:bounds[g + 1]] for g in range(G)]52 reps = {m: np.empty(B) for m in models}53 for b in range(B):54 draw = rng.integers(0, G, G)55 idx = np.concatenate([rows_by_cluster[g] for g in draw])56 for m in models:57 reps[m][b] = np.median(A[m][idx])58 out = []59 for m in models:60 point = float(np.median(A[m]) * 100)61 lo, hi = np.percentile(reps[m], [2.5, 97.5]) * 10062 row = {"model": m, "mdape": point, "ci_lo": lo, "ci_hi": hi,63 "boot_se": float(reps[m].std(ddof=1) * 100)}64 if ref is not None and m != ref:65 d = (reps[m] - reps[ref]) * 10066 row.update({"diff_vs_ref": point - float(np.median(A[ref]) * 100),67 "diff_lo": float(np.percentile(d, 2.5)),68 "diff_hi": float(np.percentile(d, 97.5))})69 out.append(row)70 return pd.DataFrame(out)717273def dm_test(P: pd.DataFrame, model_a: str, model_b: str,74 loss=ape, cluster: str = "muni") -> dict:75 """Clustered test of E[d_i] = 0 for d_i = L_A(i) - L_B(i).7677 Returns the mean differential (positive = A has larger loss, i.e.78 B is more accurate), its cluster-robust standard error, the t79 statistic and the two-sided normal p-value.80 """81 y = P["price"].to_numpy(float)82 d = loss(y, P[model_a].to_numpy(float)) - loss(y, P[model_b].to_numpy(float))83 codes, _ = pd.factorize(P[cluster])84 n, G = len(d), codes.max() + 185 dbar = d.mean()86 # cluster-robust variance of the mean: sum over clusters of squared87 # within-cluster sums of demeaned d, divided by n^288 sums = np.bincount(codes, weights=d - dbar, minlength=G)89 var = (sums ** 2).sum() / n ** 2 * G / (G - 1)90 se = float(np.sqrt(var))91 t = float(dbar / se) if se > 0 else np.nan92 from math import erf, sqrt93 p = 2 * (1 - 0.5 * (1 + erf(abs(t) / sqrt(2)))) if np.isfinite(t) else np.nan94 return {"model_a": model_a, "model_b": model_b, "mean_diff": float(dbar),95 "se": se, "t": t, "p": float(p), "n_clusters": int(G)}96