#!/usr/bin/env python3 # Author: Simon-Pierre Boucher — contact@spboucher.ai # """Inference on the corrected design (November 2026 review, batch L2; items WP03-DATA-06, EDIT-04, ADVE-05, ADVE-08, ECON-02, THEO-06, ECON-04, THEO-05, EDIT-07, ADVE-14, ADVE-12, ECON-10, THEO-08, ECON-07, EDIT-09, ECON-08, ADVE-23, EDIT-02). Refactor of 06_revision_inference.py on the v2 design (trimmed sample). Steps (--steps to select; default all) assert zpid / split-index assertions (DATA-06) cluster OLS with HC3 / state-CR1 / ZIP3-CR1 SEs; counts over the 62 slopes with |t| > 1.96, t(G-1), Holm, BH; CSS effective number of clusters; '000' cluster reported and reassigned to the nearest ZIP3 centroid (SEs compared); clustered Kennedy effects for the dummies wild wild cluster bootstrap by state (wildboottest, Rademacher and Webb, B = 9,999, seed 42) for every slope states state-by-state OLS with ZIP3-cluster SEs, reversal counts (|t| > 2; state CI excludes the pooled value), cluster-robust homogeneity test, EB shrinkage lot lot-size codings on the v2 sample (legacy vs v2, trimmed and untrimmed) moran Moran's I on the full random test set (cKDTree k = 8, 999 permutations) and by distance band (0-1, 1-5, 5-25 km) for OLS / LightGBM / XGBoost residuals (stored v1.2 predictions and, when present, v2 predictions from batch N2) conley Conley HAC (Bartlett, 25 / 50 / 100 km) on a 100,000-listing subsample paired ZIP3-cluster SE of paired squared-error differences on the random split Outputs: results/inference_v2/*.csv Usage: python scripts/06_inference_v2.py [--steps assert cluster wild states lot moran conley paired] [--B 9999] [--workers 8] """ import argparse import json import sys import time import warnings from pathlib import Path import numpy as np import pandas as pd import statsmodels.api as sm from scipy import sparse, stats from scipy.spatial import cKDTree from sklearn.metrics import r2_score from sklearn.model_selection import train_test_split sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) from wp3 import config, data, design # noqa: E402 OUT = config.RESULTS_DIR / "inference_v2" KEY24 = ["ln_sqft", "bedrooms", "bathrooms", "age", "age_sq", "ln_lot", "lot_missing", "has_pool", "has_garage", "on_waterfront", "has_central_air", "has_hardwood", "luxury_score", "walk_score", "bike_score", "transit_score", "avg_school_rating", "property_tax_rate", "is_condo", "has_hoa", "tag_new_construction", "tag_foreclosure", "region_Northeast", "region_South", "region_West"] REPORT13 = ["ln_sqft", "bedrooms", "bathrooms", "age", "ln_lot", "has_pool", "has_garage", "has_central_air", "has_hardwood", "luxury_score", "tag_foreclosure", "is_condo", "has_hoa"] DUMMIES = ["on_waterfront", "has_pool", "has_garage", "has_central_air", "has_hardwood", "is_condo", "tag_foreclosure", "tag_new_construction", "has_hoa", "has_basement", "has_spa", "has_fireplace", "lot_missing"] def holm(p): p = np.asarray(p, float); m = len(p); order = np.argsort(p); adj = np.empty(m); prev = 0.0 for i, j in enumerate(order): v = max(min(1.0, (m - i) * p[j]), prev); adj[j] = v; prev = v return adj def bh(p): p = np.asarray(p, float); m = len(p); order = np.argsort(p); adj = np.empty(m); prev = 1.0 for i in range(m - 1, -1, -1): j = order[i]; v = min(min(1.0, m / (i + 1) * p[j]), prev); adj[j] = v; prev = v return adj def to_xy_km(lat, lon): """Equirectangular projection (km), adequate for neighbour searches below ~100 km.""" lat0 = np.deg2rad(np.mean(lat)) return np.column_stack([np.deg2rad(lon) * 6371.0 * np.cos(lat0), np.deg2rad(lat) * 6371.0]) # ----------------------------------------------------------------------------- assert def step_assert(): sample = pd.read_parquet(config.DATA_PROCESSED / "analytical_sample_v2.parquet", columns=["zpid"]) d = pd.read_parquet(design.DESIGN_V2, columns=["zpid"]) assert (sample["zpid"].values == d["zpid"].values).all(), "design_v2 / analytical_sample_v2 zpid order differs" rows = {"design_vs_sample_zpid_order": True} if config.ANALYTICAL_SAMPLE.exists() and config.MODEL_DATA.exists(): pk = pd.read_pickle(config.ANALYTICAL_SAMPLE) md = data.load_model_data() assert (md["df_clean"]["zpid"].values == pk["zpid"].values).all(), "model_data.pkl / analytical_sample.pkl zpid order differs" assert (pk["zpid"].values == d["zpid"].values).all() rows["model_data_vs_pickle_zpid_order"] = True try: ml = data.load_result("ml_results.pkl") idx_tr, idx_te = train_test_split(np.arange(len(pk)), test_size=0.2, random_state=config.RANDOM_STATE) assert (ml["X_test"].index.values == idx_te).all(), "ml_results X_test index != train_test_split(seed 42)" assert ml["X_test"].index.isin(pk.index).all() rows["ml_results_test_index_equals_seed42_split"] = True except FileNotFoundError: rows["ml_results_test_index_equals_seed42_split"] = "ml_results.pkl absent" pd.Series(rows).to_csv(OUT / "assertions.csv", header=["passed"]) print("assertions:", rows) # ----------------------------------------------------------------------------- cluster def css_effective_clusters(X, groups): """Carter-Schnepel-Steigerwald (2017) feasible effective number of clusters per coefficient (homoskedastic within-cluster errors, Omega_g = I): gamma_g = e_k' A X_g'X_g A e_k.""" A = np.linalg.inv(X.T @ X) codes, uniq = pd.factorize(groups) G = len(uniq); K = X.shape[1] gam = np.zeros((G, K)) for g in range(G): Xg = X[codes == g] M = A @ (Xg.T @ Xg) @ A gam[g] = np.diag(M) gbar = gam.mean(axis=0) Gamma = ((gam - gbar) ** 2 / gbar ** 2).mean(axis=0) return G / (1 + Gamma) def step_cluster(): X, y, info, ex = design.load_design("v2", trimmed=True, standardize=True) Xc = sm.add_constant(X, has_constant="add") names = list(Xc.columns) st = pd.factorize(info["address_state"].values)[0] z3 = pd.factorize(info["zip3"].values)[0] base = sm.OLS(y, Xc.values) r_hc3 = base.fit(cov_type="HC3") r_st = base.fit(cov_type="cluster", cov_kwds={"groups": st}) r_z3 = base.fit(cov_type="cluster", cov_kwds={"groups": z3}) out = pd.DataFrame({"coef": r_hc3.params, "se_hc3": r_hc3.bse, "se_cluster_state": r_st.bse, "se_cluster_zip3": r_z3.bse}, index=names) out["ratio_zip3_hc3"] = out.se_cluster_zip3 / out.se_hc3; out["ratio_state_hc3"] = out.se_cluster_state / out.se_hc3 for k in ["hc3", "cluster_state", "cluster_zip3"]: out[f"t_{k}"] = out.coef / out[f"se_{k}"] out["p_state_t50"] = 2 * (1 - stats.t.cdf(out.t_cluster_state.abs(), 50)) out["p_zip3_t885"] = 2 * (1 - stats.t.cdf(out.t_cluster_zip3.abs(), len(np.unique(z3)) - 1)) out["p_hc3_normal"] = 2 * (1 - stats.norm.cdf(out.t_hc3.abs())) slopes = out.index != "const" for k, pcol in [("state", "p_state_t50"), ("zip3", "p_zip3_t885"), ("hc3", "p_hc3_normal")]: out.loc[slopes, f"p_holm_{k}"] = holm(out.loc[slopes, pcol]); out.loc[slopes, f"p_bh_{k}"] = bh(out.loc[slopes, pcol]) # CSS effective clusters (state) out["G_eff_css_state"] = css_effective_clusters(Xc.values, info["address_state"].values) out["G_eff_css_zip3"] = css_effective_clusters(Xc.values, info["zip3"].values) out["n_clusters_state"] = len(np.unique(st)); out["n_clusters_zip3"] = len(np.unique(z3)) out.index.name = "variable" out.to_csv(OUT / "cluster_se_v2.csv") s = out[slopes] counts = [] for k, G in [("hc3", None), ("cluster_zip3", len(np.unique(z3))), ("cluster_state", 51)]: t = s[f"t_{k}"].abs() crit = stats.t.ppf(0.975, G - 1) if G else 1.96 pc = {"hc3": "p_hc3_normal", "cluster_zip3": "p_zip3_t885", "cluster_state": "p_state_t50"}[k] kk = k.replace("cluster_", "") counts.append({"se": k, "n_slopes": int(slopes.sum()), "n_abs_t_gt_196": int((t > 1.96).sum()), "critical_t": crit, "n_abs_t_gt_crit": int((t > crit).sum()), "n_holm_5pct": int((s[f"p_holm_{kk}"] < 0.05).sum()), "n_bh_5pct": int((s[f"p_bh_{kk}"] < 0.05).sum()), "median_G_eff_css": float(s["G_eff_css_state" if kk == "state" else "G_eff_css_zip3"].median()) if G else np.nan, "min_G_eff_css": float(s["G_eff_css_state" if kk == "state" else "G_eff_css_zip3"].min()) if G else np.nan}) pd.DataFrame(counts).to_csv(OUT / "significance_counts_v2.csv", index=False) print(pd.DataFrame(counts).to_string(index=False)) # '000' cluster: size and nearest-ZIP3 reassignment miss = info["zip3"].values == "000" cent = info.loc[~miss].groupby("zip3")[["latitude", "longitude"]].mean() tree = cKDTree(to_xy_km(cent.latitude.values, cent.longitude.values)) _, j = tree.query(to_xy_km(info.loc[miss, "latitude"].values, info.loc[miss, "longitude"].values)) z3b = info["zip3"].values.copy(); z3b[miss] = cent.index.values[j] r_z3b = base.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(z3b)[0]}) cmp = pd.DataFrame({"se_cluster_zip3_with_000": r_z3.bse, "se_cluster_zip3_reassigned": r_z3b.bse}, index=names) cmp["rel_change"] = cmp.se_cluster_zip3_reassigned / cmp.se_cluster_zip3_with_000 - 1 cmp.index.name = "variable"; cmp.to_csv(OUT / "zip3_missing_reassignment.csv") pd.DataFrame([{"n_listings_missing_zip": int(miss.sum()), "n_zip3_prefixes": int(len(cent)), "n_clusters_with_000": int(len(np.unique(z3))), "max_abs_rel_se_change_after_reassignment": float(cmp.rel_change.abs().max()), "reassigned_to": " ".join(sorted(set(cent.index.values[j])))}]).to_csv(OUT / "zip3_missing_cluster.csv", index=False) # Kennedy with clustered variance (natural units) Xu, yu, infou, _ = design.load_design("v2", trimmed=True, standardize=False) Xcu = sm.add_constant(Xu, has_constant="add") bu = sm.OLS(yu, Xcu.values) ru = {"hc3": bu.fit(cov_type="HC3"), "cluster_zip3": bu.fit(cov_type="cluster", cov_kwds={"groups": z3}), "cluster_state": bu.fit(cov_type="cluster", cov_kwds={"groups": st})} kr = [] for v in DUMMIES: jv = list(Xcu.columns).index(v) for k, r in ru.items(): b, se_ = r.params[jv], r.bse[jv] kr.append({"variable": v, "se": k, "coef": b, "se_value": se_, "naive_pct": 100 * (np.exp(b) - 1), "kennedy_pct": 100 * (np.exp(b - 0.5 * se_ ** 2) - 1), "t": b / se_}) pd.DataFrame(kr).to_csv(OUT / "kennedy_clustered_v2.csv", index=False) unstd = pd.DataFrame({"coef": ru["hc3"].params, "se_hc3": ru["hc3"].bse, "se_cluster_zip3": ru["cluster_zip3"].bse, "se_cluster_state": ru["cluster_state"].bse}, index=list(Xcu.columns)) unstd.index.name = "variable"; unstd.to_csv(OUT / "ols_v2_unstandardized_cluster_se.csv") # ----------------------------------------------------------------------------- wild def step_wild(B, workers): from wildboottest.wildboottest import wildboottest X, y, info, ex = design.load_design("v2", trimmed=True, standardize=True) Xc = sm.add_constant(X, has_constant="add") model = sm.OLS(y, Xc) cl = pd.Series(pd.factorize(info["address_state"].values)[0]) rows = [] t0 = time.time() for wt in ["rademacher", "webb"]: for v in [c for c in Xc.columns if c != "const"]: res = wildboottest(model, cluster=cl, B=B, param=v, weights_type=wt, seed=config.RANDOM_STATE, show=False) r = res.iloc[0].to_dict() rows.append({"variable": v, "weights": wt, **{k: r[k] for k in r if k in ("statistic", "p-value", "param")}, "B": B}) print(f" wild bootstrap ({wt}) done in {time.time() - t0:.0f}s", flush=True) out = pd.DataFrame(rows).rename(columns={"p-value": "p_wild", "statistic": "t_cr1"}) wide = out.pivot(index="variable", columns="weights", values="p_wild").rename(columns=lambda c: f"p_wild_{c}") tt = out[out.weights == "webb"].set_index("variable")["t_cr1"] wide["t_cr1"] = tt wide["p_holm_webb"] = holm(wide["p_wild_webb"]); wide["p_bh_webb"] = bh(wide["p_wild_webb"]) wide["p_holm_rademacher"] = holm(wide["p_wild_rademacher"]); wide["p_bh_rademacher"] = bh(wide["p_wild_rademacher"]) wide.to_csv(OUT / "wild_bootstrap_state.csv") cnt = {"n_slopes": len(wide), "n_sig_5pct_webb": int((wide.p_wild_webb < 0.05).sum()), "n_sig_5pct_rademacher": int((wide.p_wild_rademacher < 0.05).sum()), "n_sig_5pct_webb_holm": int((wide.p_holm_webb < 0.05).sum()), "n_sig_5pct_webb_bh": int((wide.p_bh_webb < 0.05).sum()), "n_sig_5pct_rademacher_holm": int((wide.p_holm_rademacher < 0.05).sum()), "B": B} pd.DataFrame([cnt]).to_csv(OUT / "wild_bootstrap_counts.csv", index=False) print(cnt) # ----------------------------------------------------------------------------- states def step_states(): Xu, y, info, _ = design.load_design("v2", trimmed=True, standardize=False) cols = [c for c in Xu.columns if not c.startswith("region_")] states = info["address_state"].values; z3 = info["zip3"].values idx_tr, idx_te = train_test_split(np.arange(788_842), test_size=0.2, random_state=config.RANDOM_STATE) # map the full-sample split onto the trimmed sample via zpid full = pd.read_parquet(design.DESIGN_V2, columns=["zpid"]) test_z = set(full["zpid"].values[idx_te]) is_test = info["zpid"].isin(test_z).values rows, se_rows, dropped = [], [], [] pred_inter = np.full(len(y), np.nan) for s_ in sorted(np.unique(states)): m = states == s_ Xs = Xu.loc[m, cols] keep = [c for c in cols if Xs[c].nunique() > 1] dropped.append({"state": s_, "n": int(m.sum()), "dropped_constant_columns": " ".join(c for c in cols if c not in keep)}) Xc = sm.add_constant(Xs[keep].values, has_constant="add") g = pd.factorize(z3[m])[0] r = sm.OLS(y[m], Xc).fit(cov_type="cluster", cov_kwds={"groups": g}) if len(np.unique(g)) > 1 else sm.OLS(y[m], Xc).fit(cov_type="HC3") cond = np.linalg.cond(Xc) row = {"state": s_, "n": int(m.sum()), "n_zip3": int(len(np.unique(g))), "r2": r.rsquared, "cond_number": cond, "median_price": float(np.exp(np.median(y[m])))} for v in REPORT13: if v in keep: j = keep.index(v) + 1 row[v] = r.params[j]; row[f"se_{v}"] = r.bse[j] else: row[v] = np.nan; row[f"se_{v}"] = np.nan rows.append(row) tr = m & ~is_test; te = m & is_test if tr.sum() > len(keep) + 5 and te.sum() > 0: from sklearn.linear_model import LinearRegression pred_inter[te] = LinearRegression().fit(Xu.loc[tr, keep].values, y[tr]).predict(Xu.loc[te, keep].values) het = pd.DataFrame(rows); het.to_csv(OUT / "state_heterogeneity_v2.csv", index=False) pd.DataFrame(dropped).to_csv(OUT / "state_regressions_dropped_columns.csv", index=False) # pooled state-FE reference with ZIP3-cluster SEs sfe = pd.get_dummies(pd.Series(states), drop_first=True).astype(float).values Xp = np.column_stack([np.ones(len(y)), Xu[cols].values, sfe]) pooled = sm.OLS(y, Xp).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(z3)[0]}) summ = [] for v in REPORT13: jp = cols.index(v) + 1 b0, s0 = pooled.params[jp], pooled.bse[jp] b = het[v].values; se_ = het[f"se_{v}"].values; ok = ~np.isnan(b) & ~np.isnan(se_) & (se_ > 0) b, se_ = b[ok], se_[ok] w = 1 / se_ ** 2 b_pw = np.sum(w * b) / np.sum(w) # cluster-robust homogeneity (states independent; ZIP3 clusters nested in states) Q = float(np.sum(w * (b - b_pw) ** 2)); dfq = len(b) - 1 tau2 = max(0.0, (Q - dfq) / (np.sum(w) - np.sum(w ** 2) / np.sum(w))) # DerSimonian-Laird summ.append({"variable": v, "pooled_state_fe": b0, "pooled_se_zip3": s0, "n_states": int(len(b)), "precision_weighted_mean": b_pw, "listing_weighted_mean": np.average(het[v].values[ok], weights=het["n"].values[ok]), "p10": np.percentile(b, 10), "median": np.median(b), "p90": np.percentile(b, 90), "min": b.min(), "max": b.max(), "n_sign_differs_from_pooled": int((np.sign(b) != np.sign(b0)).sum()), "n_sign_differs_and_abs_t_gt_2": int(((np.sign(b) != np.sign(b0)) & (np.abs(b / se_) > 2)).sum()), "n_state_ci_excludes_pooled": int(((b - 1.96 * se_ > b0) | (b + 1.96 * se_ < b0)).sum()), "n_abs_t_gt_2": int((np.abs(b / se_) > 2).sum()), "homogeneity_Q": Q, "homogeneity_df": dfq, "homogeneity_p": float(1 - stats.chi2.cdf(Q, dfq)), "tau2_dl": tau2, "tau_dl": np.sqrt(tau2), "I2_pct": max(0.0, 100 * (Q - dfq) / Q) if Q > 0 else 0.0}) summ = pd.DataFrame(summ); summ.to_csv(OUT / "state_heterogeneity_summary_v2.csv", index=False) print(summ[["variable", "pooled_state_fe", "p10", "median", "p90", "n_sign_differs_from_pooled", "n_sign_differs_and_abs_t_gt_2", "n_state_ci_excludes_pooled", "homogeneity_p", "tau_dl"]].round(3).to_string(index=False)) # empirical-Bayes shrinkage for Figure 3 (normal-normal, tau2 DL) eb = [] for v in REPORT13: t2 = float(summ.set_index("variable").loc[v, "tau2_dl"]); mu = float(summ.set_index("variable").loc[v, "precision_weighted_mean"]) for _, r in het.iterrows(): b, s_ = r[v], r[f"se_{v}"] if np.isnan(b) or np.isnan(s_) or s_ <= 0: continue shrink = t2 / (t2 + s_ ** 2) post = shrink * b + (1 - shrink) * mu; pv = shrink * s_ ** 2 eb.append({"variable": v, "state": r.state, "n": r.n, "estimate": b, "se_zip3": s_, "ci_low": b - 1.96 * s_, "ci_high": b + 1.96 * s_, "eb_posterior_mean": post, "eb_posterior_sd": np.sqrt(pv), "eb_low": post - 1.96 * np.sqrt(pv), "eb_high": post + 1.96 * np.sqrt(pv), "shrinkage_factor": shrink}) pd.DataFrame(eb).to_csv(OUT / "state_heterogeneity_eb.csv", index=False) # pooled vs state FE vs fully interacted on the random test split (trimmed) from sklearn.linear_model import LinearRegression te = is_test; tr = ~is_test p_pooled = LinearRegression().fit(Xu.loc[tr, cols].values, y[tr]).predict(Xu.loc[te, cols].values) gtr = pd.Series(states[tr]); Xd = Xu.loc[tr, cols].values - pd.DataFrame(Xu.loc[tr, cols].values).groupby(gtr.values).transform("mean").values yd = y[tr] - pd.Series(y[tr]).groupby(gtr.values).transform("mean").values beta, *_ = np.linalg.lstsq(Xd, yd, rcond=None) alpha = pd.Series(y[tr] - Xu.loc[tr, cols].values @ beta).groupby(gtr.values).mean() p_fe = Xu.loc[te, cols].values @ beta + pd.Series(states[te]).map(alpha).fillna(alpha.mean()).values ok = ~np.isnan(pred_inter[te]) comp = pd.DataFrame([{"model": "Pooled OLS, no geographic controls", "r2_test": r2_score(y[te], p_pooled)}, {"model": "Pooled OLS + 50 state fixed effects", "r2_test": r2_score(y[te], p_fe)}, {"model": "State-by-state OLS (fully interacted)", "r2_test": r2_score(y[te][ok], pred_inter[te][ok])}]) comp["n_train"] = int(tr.sum()); comp["n_test"] = int(te.sum()) comp.to_csv(OUT / "pooled_vs_interacted_v2.csv", index=False); print(comp.to_string(index=False)) # ----------------------------------------------------------------------------- lot def step_lot(): rows = [] for trimmed in [False, True]: Xl, y, info, _ = design.load_design("legacy", trimmed=trimmed, standardize=False, lot_indicator=False) Xv, _, _, _ = design.load_design("v2", trimmed=trimmed, standardize=False, lot_indicator=False) Xvi, _, _, _ = design.load_design("v2", trimmed=trimmed, standardize=False, lot_indicator=True) lm = info["lot_missing"].values == 1 def fit(X, mask, label): Xc = sm.add_constant(X.values[mask], has_constant="add") r = sm.OLS(y[mask], Xc).fit(cov_type="HC3") names = list(X.columns) return {"sample": "trimmed (783,872)" if trimmed else "all (788,842)", "specification": label, "N": int(mask.sum()), "R2": r.rsquared, "beta_ln_lot": r.params[names.index("ln_lot") + 1], "se_ln_lot": r.bse[names.index("ln_lot") + 1], "beta_ln_sqft": r.params[names.index("ln_sqft") + 1], "beta_lot_missing": r.params[names.index("lot_missing") + 1] if "lot_missing" in names else np.nan} allm = np.ones(len(y), bool) rows.append(fit(Xl, allm, "v1.1 coding: ln(lot) = 0 when missing, zero, 1 sq ft or negative")) rows.append(fit(Xv, allm, "v2: ln(state-median lot) when missing/zero/<= 1")) rows.append(fit(Xvi, allm, "v2 + lot-missing indicator (the design used throughout)")) rows.append(fit(Xl, ~lm, "Complete cases (raw lot > 1 sq ft)")) p995 = info["lot_size_sqft"].quantile(0.995) Xw = Xvi.copy(); lot_v2 = np.exp(Xvi["ln_lot"].values); Xw["ln_lot"] = np.log(np.minimum(lot_v2, p995)) rows.append(fit(Xw, allm, "v2 + indicator, lot winsorised at p99.5")) out = pd.DataFrame(rows); out.to_csv(OUT / "lot_recoding_v2.csv", index=False); print(out.round(4).to_string(index=False)) # ----------------------------------------------------------------------------- moran def moran_knn(e, xy, k=8, permutations=999, seed=42): from esda.moran import Moran from libpysal.weights import WSP tree = cKDTree(xy) d, j = tree.query(xy, k=k + 1) n = len(xy) rows_ = np.repeat(np.arange(n), k); cols_ = j[:, 1:].ravel() W = sparse.csr_matrix((np.ones(n * k) / k, (rows_, cols_)), shape=(n, n)) w = WSP(W).to_W() np.random.seed(seed) mi = Moran(e, w, transformation="r", permutations=permutations) return {"I": mi.I, "E_I": mi.EI, "z_norm": mi.z_norm, "p_norm": mi.p_norm, "z_sim": mi.z_sim, "p_sim": mi.p_sim, "median_kth_neighbour_km": float(np.median(d[:, k]))} def moran_band(e, pairs_i, pairs_j, n): """Moran's I for a symmetric binary band weight matrix given as pair lists; normal-approx z (randomisation).""" z = e - e.mean() S0 = 2 * len(pairs_i) num = 2 * np.sum(z[pairs_i] * z[pairs_j]) I = (n / S0) * num / np.sum(z ** 2) EI = -1 / (n - 1) deg = np.bincount(np.concatenate([pairs_i, pairs_j]), minlength=n).astype(float) S1 = 0.5 * np.sum((2.0) ** 2) * len(pairs_i) # each symmetric pair contributes (w_ij + w_ji)^2 = 4, times 1/2 S2 = np.sum((2 * deg) ** 2) b2 = n * np.sum(z ** 4) / np.sum(z ** 2) ** 2 VI = (n * ((n ** 2 - 3 * n + 3) * S1 - n * S2 + 3 * S0 ** 2) - b2 * ((n ** 2 - n) * S1 - 2 * n * S2 + 6 * S0 ** 2)) / \ ((n - 1) * (n - 2) * (n - 3) * S0 ** 2) - EI ** 2 zI = (I - EI) / np.sqrt(VI) return {"I": I, "E_I": EI, "z_rand": zI, "p_rand": float(2 * (1 - stats.norm.cdf(abs(zI)))), "n_pairs": int(len(pairs_i)), "share_units_with_neighbour": float((deg > 0).mean())} def step_moran(permutations): rows = [] sets = [] # v1.2 stored predictions (random test set, legacy design) try: ml = data.load_result("ml_results.pkl"); pk = pd.read_pickle(config.ANALYTICAL_SAMPLE) idx = ml["X_test"].index.values xy = to_xy_km(pk.loc[idx, "latitude"].values, pk.loc[idx, "longitude"].values) y = ml["y_test"].values for name, p in [("OLS", ml["y_pred_ols"]), ("LightGBM", ml["y_pred_lgb"]), ("XGBoost", ml["y_pred_xgb"])]: sets.append(("v1.2 stored predictions, random test set", name, y - p, xy)) except FileNotFoundError: pass # v2 predictions from batch N2 when available prim = config.RESULTS_DIR / "primary" d = pd.read_parquet(design.DESIGN_V2, columns=["zpid", "latitude", "longitude"]).set_index("zpid") for name, pat in [("OLS", "predictions_random_OLS_62_regressors_+_lot_indicator.parquet"), ("LightGBM", "predictions_random_LightGBM_1000_rounds_early_stopping_seed42.parquet"), ("XGBoost", "predictions_random_XGBoost_1000_rounds_early_stopping_seed42.parquet")]: f = prim / pat if f.exists(): p = pd.read_parquet(f) xy = to_xy_km(d.loc[p.zpid, "latitude"].values, d.loc[p.zpid, "longitude"].values) sets.append(("v2 predictions (N2), random test set, trimmed", name, p["y"].values - p["pred"].values, xy)) bands = [(0, 1), (1, 5), (5, 25)] for label, name, e, xy in sets: t0 = time.time() r = moran_knn(e, xy, k=8, permutations=permutations) rows.append({"predictions": label, "model": name, "weights": "kNN k=8 row-standardised", "n": len(e), **r}) print(f" {label} {name} kNN I={r['I']:.4f} z_sim={r['z_sim']:.1f} ({time.time() - t0:.0f}s)", flush=True) tree = cKDTree(xy) prev = None for lo, hi in bands: pairs = tree.query_pairs(hi, output_type="ndarray") if lo > 0: dd = np.linalg.norm(xy[pairs[:, 0]] - xy[pairs[:, 1]], axis=1) pairs = pairs[dd > lo] rb = moran_band(e, pairs[:, 0], pairs[:, 1], len(e)) rows.append({"predictions": label, "model": name, "weights": f"distance band {lo}-{hi} km (binary)", "n": len(e), **rb}) print(f" band {lo}-{hi} km: I={rb['I']:.4f} z={rb['z_rand']:.1f} pairs={rb['n_pairs']:,}", flush=True) out = pd.DataFrame(rows); out.to_csv(OUT / "moran_full_test_v2.csv", index=False) # ----------------------------------------------------------------------------- conley def step_conley(n_sub=100_000, cutoffs=(25, 50, 100)): X, y, info, ex = design.load_design("v2", trimmed=True, standardize=True) rng = np.random.RandomState(config.RANDOM_STATE) idx = np.sort(rng.choice(len(y), n_sub, replace=False)) Xc = sm.add_constant(X.iloc[idx], has_constant="add").values; yy = y[idx] names = ["const"] + list(X.columns) r = sm.OLS(yy, Xc).fit(cov_type="HC0") e = r.resid A = np.linalg.inv(Xc.T @ Xc) xy = to_xy_km(info["latitude"].values[idx], info["longitude"].values[idx]) tree = cKDTree(xy) out = pd.DataFrame({"coef": r.params, "se_hc0": r.bse}, index=names) st = pd.factorize(info["address_state"].values[idx])[0] out["se_cluster_state"] = sm.OLS(yy, Xc).fit(cov_type="cluster", cov_kwds={"groups": st}).bse Xe = Xc * e[:, None] for c in cutoffs: t0 = time.time() pairs = tree.query_pairs(c, output_type="ndarray") dd = np.linalg.norm(xy[pairs[:, 0]] - xy[pairs[:, 1]], axis=1) k = 1 - dd / c # Bartlett kernel n = len(yy) K = sparse.coo_matrix((np.concatenate([k, k]), (np.concatenate([pairs[:, 0], pairs[:, 1]]), np.concatenate([pairs[:, 1], pairs[:, 0]]))), shape=(n, n)).tocsr() M = Xe.T @ Xe + Xe.T @ (K @ Xe) V = A @ M @ A out[f"se_conley_{c}km"] = np.sqrt(np.diag(V)) print(f" Conley {c} km: {len(pairs):,} pairs ({time.time() - t0:.0f}s)", flush=True) for c in cutoffs: out[f"ratio_conley_{c}km_hc0"] = out[f"se_conley_{c}km"] / out["se_hc0"] out["n_subsample"] = n_sub; out.index.name = "variable" out.to_csv(OUT / "conley_hac_v2.csv") print(out.loc[[v for v in KEY24 if v in out.index], ["coef", "se_hc0", "se_cluster_state"] + [f"se_conley_{c}km" for c in cutoffs]].round(4).to_string()) # ----------------------------------------------------------------------------- paired def step_paired(): prim = config.RESULTS_DIR / "primary" files = {p.name: p for p in prim.glob("predictions_random_*.parquet")} if not files: print(" no v2 random-split predictions yet (batch N2) — using stored v1.2 predictions") ml = data.load_result("ml_results.pkl"); pk = pd.read_pickle(config.ANALYTICAL_SAMPLE) z3 = pk.loc[ml["X_test"].index, "address_zipcode"].fillna("000").astype(str).str[:3].values y = ml["y_test"].values; preds = {"OLS": ml["y_pred_ols"], "LightGBM": ml["y_pred_lgb"], "XGBoost": ml["y_pred_xgb"]}; label = "v1.2 stored" else: d = pd.read_parquet(design.DESIGN_V2, columns=["zpid", "zip3"]).set_index("zpid") preds = {} for name, pat in [("OLS", "predictions_random_OLS_62_regressors_+_lot_indicator.parquet"), ("LightGBM", "predictions_random_LightGBM_1000_rounds_early_stopping_seed42.parquet"), ("XGBoost", "predictions_random_XGBoost_1000_rounds_early_stopping_seed42.parquet")]: p = pd.read_parquet(prim / pat); preds[name] = p["pred"].values; y = p["y"].values; zp = p["zpid"].values z3 = d.loc[zp, "zip3"].values; label = "v2 (N2)" g = pd.factorize(z3)[0] rows = [] for a, c in [("XGBoost", "OLS"), ("XGBoost", "LightGBM"), ("LightGBM", "OLS")]: dl = (y - preds[c]) ** 2 - (y - preds[a]) ** 2 r = sm.OLS(dl, np.ones(len(dl))).fit(cov_type="cluster", cov_kwds={"groups": g}) r0 = sm.OLS(dl, np.ones(len(dl))).fit(cov_type="HC0") rows.append({"predictions": label, "quantity": f"{c} minus {a} squared-error difference", "mean": float(dl.mean()), "se_iid": float(dl.std(ddof=1) / np.sqrt(len(dl))), "se_hc0": float(r0.bse[0]), "se_cluster_zip3": float(r.bse[0]), "t_cluster_zip3": float(r.params[0] / r.bse[0]), "n_clusters": int(len(np.unique(g)))}) pd.DataFrame(rows).to_csv(OUT / "paired_loss_cluster_se.csv", index=False); print(pd.DataFrame(rows).round(5).to_string(index=False)) def main(): ap = argparse.ArgumentParser(description=__doc__) ap.add_argument("--steps", nargs="+", default=["assert", "cluster", "wild", "states", "lot", "moran", "conley", "paired"]) ap.add_argument("--B", type=int, default=9999) ap.add_argument("--workers", type=int, default=8) ap.add_argument("--permutations", type=int, default=999) args = ap.parse_args() OUT.mkdir(parents=True, exist_ok=True) warnings.filterwarnings("ignore") import os os.environ.setdefault("OMP_NUM_THREADS", str(args.workers)) t0 = time.time() for s in args.steps: print(f"== {s} ({time.time() - t0:.0f}s)", flush=True) {"assert": step_assert, "cluster": step_cluster, "wild": lambda: step_wild(args.B, args.workers), "states": step_states, "lot": step_lot, "moran": lambda: step_moran(args.permutations), "conley": step_conley, "paired": step_paired}[s]() (OUT / "run_meta.json").write_text(json.dumps({"steps": args.steps, "B": args.B, "elapsed_s": round(time.time() - t0)}, indent=1)) print(f"Done in {time.time() - t0:.0f}s") if __name__ == "__main__": main()