UQO Working Paper No. 3 — Hedonic housing price models for the US: parametric, quantile, and machine-learning approaches.
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3#4"""Inference on the corrected design (November 2026 review, batch L2; items5WP03-DATA-06, EDIT-04, ADVE-05, ADVE-08, ECON-02, THEO-06, ECON-04, THEO-05,6EDIT-07, ADVE-14, ADVE-12, ECON-10, THEO-08, ECON-07, EDIT-09, ECON-08, ADVE-23,7EDIT-02). Refactor of 06_revision_inference.py on the v2 design (trimmed sample).89Steps (--steps to select; default all)10 assert zpid / split-index assertions (DATA-06)11 cluster OLS with HC3 / state-CR1 / ZIP3-CR1 SEs; counts over the 62 slopes with12 |t| > 1.96, t(G-1), Holm, BH; CSS effective number of clusters; '000' cluster13 reported and reassigned to the nearest ZIP3 centroid (SEs compared);14 clustered Kennedy effects for the dummies15 wild wild cluster bootstrap by state (wildboottest, Rademacher and Webb,16 B = 9,999, seed 42) for every slope17 states state-by-state OLS with ZIP3-cluster SEs, reversal counts (|t| > 2; state CI18 excludes the pooled value), cluster-robust homogeneity test, EB shrinkage19 lot lot-size codings on the v2 sample (legacy vs v2, trimmed and untrimmed)20 moran Moran's I on the full random test set (cKDTree k = 8, 999 permutations) and21 by distance band (0-1, 1-5, 5-25 km) for OLS / LightGBM / XGBoost residuals22 (stored v1.2 predictions and, when present, v2 predictions from batch N2)23 conley Conley HAC (Bartlett, 25 / 50 / 100 km) on a 100,000-listing subsample24 paired ZIP3-cluster SE of paired squared-error differences on the random split2526Outputs: results/inference_v2/*.csv2728Usage: python scripts/06_inference_v2.py [--steps assert cluster wild states lot moran conley paired]29 [--B 9999] [--workers 8]30"""3132import argparse33import json34import sys35import time36import warnings37from pathlib import Path3839import numpy as np40import pandas as pd41import statsmodels.api as sm42from scipy import sparse, stats43from scipy.spatial import cKDTree44from sklearn.metrics import r2_score45from sklearn.model_selection import train_test_split4647sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))48from wp3 import config, data, design # noqa: E4024950OUT = config.RESULTS_DIR / "inference_v2"51KEY24 = ["ln_sqft", "bedrooms", "bathrooms", "age", "age_sq", "ln_lot", "lot_missing", "has_pool",52 "has_garage", "on_waterfront", "has_central_air", "has_hardwood", "luxury_score",53 "walk_score", "bike_score", "transit_score", "avg_school_rating", "property_tax_rate",54 "is_condo", "has_hoa", "tag_new_construction", "tag_foreclosure",55 "region_Northeast", "region_South", "region_West"]56REPORT13 = ["ln_sqft", "bedrooms", "bathrooms", "age", "ln_lot", "has_pool", "has_garage",57 "has_central_air", "has_hardwood", "luxury_score", "tag_foreclosure", "is_condo", "has_hoa"]58DUMMIES = ["on_waterfront", "has_pool", "has_garage", "has_central_air", "has_hardwood", "is_condo",59 "tag_foreclosure", "tag_new_construction", "has_hoa", "has_basement", "has_spa", "has_fireplace", "lot_missing"]606162def holm(p):63 p = np.asarray(p, float); m = len(p); order = np.argsort(p); adj = np.empty(m); prev = 0.064 for i, j in enumerate(order):65 v = max(min(1.0, (m - i) * p[j]), prev); adj[j] = v; prev = v66 return adj676869def bh(p):70 p = np.asarray(p, float); m = len(p); order = np.argsort(p); adj = np.empty(m); prev = 1.071 for i in range(m - 1, -1, -1):72 j = order[i]; v = min(min(1.0, m / (i + 1) * p[j]), prev); adj[j] = v; prev = v73 return adj747576def to_xy_km(lat, lon):77 """Equirectangular projection (km), adequate for neighbour searches below ~100 km."""78 lat0 = np.deg2rad(np.mean(lat))79 return np.column_stack([np.deg2rad(lon) * 6371.0 * np.cos(lat0), np.deg2rad(lat) * 6371.0])808182# ----------------------------------------------------------------------------- assert83def step_assert():84 sample = pd.read_parquet(config.DATA_PROCESSED / "analytical_sample_v2.parquet", columns=["zpid"])85 d = pd.read_parquet(design.DESIGN_V2, columns=["zpid"])86 assert (sample["zpid"].values == d["zpid"].values).all(), "design_v2 / analytical_sample_v2 zpid order differs"87 rows = {"design_vs_sample_zpid_order": True}88 if config.ANALYTICAL_SAMPLE.exists() and config.MODEL_DATA.exists():89 pk = pd.read_pickle(config.ANALYTICAL_SAMPLE)90 md = data.load_model_data()91 assert (md["df_clean"]["zpid"].values == pk["zpid"].values).all(), "model_data.pkl / analytical_sample.pkl zpid order differs"92 assert (pk["zpid"].values == d["zpid"].values).all()93 rows["model_data_vs_pickle_zpid_order"] = True94 try:95 ml = data.load_result("ml_results.pkl")96 idx_tr, idx_te = train_test_split(np.arange(len(pk)), test_size=0.2, random_state=config.RANDOM_STATE)97 assert (ml["X_test"].index.values == idx_te).all(), "ml_results X_test index != train_test_split(seed 42)"98 assert ml["X_test"].index.isin(pk.index).all()99 rows["ml_results_test_index_equals_seed42_split"] = True100 except FileNotFoundError:101 rows["ml_results_test_index_equals_seed42_split"] = "ml_results.pkl absent"102 pd.Series(rows).to_csv(OUT / "assertions.csv", header=["passed"])103 print("assertions:", rows)104105106# ----------------------------------------------------------------------------- cluster107def css_effective_clusters(X, groups):108 """Carter-Schnepel-Steigerwald (2017) feasible effective number of clusters per coefficient109 (homoskedastic within-cluster errors, Omega_g = I): gamma_g = e_k' A X_g'X_g A e_k."""110 A = np.linalg.inv(X.T @ X)111 codes, uniq = pd.factorize(groups)112 G = len(uniq); K = X.shape[1]113 gam = np.zeros((G, K))114 for g in range(G):115 Xg = X[codes == g]116 M = A @ (Xg.T @ Xg) @ A117 gam[g] = np.diag(M)118 gbar = gam.mean(axis=0)119 Gamma = ((gam - gbar) ** 2 / gbar ** 2).mean(axis=0)120 return G / (1 + Gamma)121122123def step_cluster():124 X, y, info, ex = design.load_design("v2", trimmed=True, standardize=True)125 Xc = sm.add_constant(X, has_constant="add")126 names = list(Xc.columns)127 st = pd.factorize(info["address_state"].values)[0]128 z3 = pd.factorize(info["zip3"].values)[0]129 base = sm.OLS(y, Xc.values)130 r_hc3 = base.fit(cov_type="HC3")131 r_st = base.fit(cov_type="cluster", cov_kwds={"groups": st})132 r_z3 = base.fit(cov_type="cluster", cov_kwds={"groups": z3})133 out = pd.DataFrame({"coef": r_hc3.params, "se_hc3": r_hc3.bse, "se_cluster_state": r_st.bse,134 "se_cluster_zip3": r_z3.bse}, index=names)135 out["ratio_zip3_hc3"] = out.se_cluster_zip3 / out.se_hc3; out["ratio_state_hc3"] = out.se_cluster_state / out.se_hc3136 for k in ["hc3", "cluster_state", "cluster_zip3"]:137 out[f"t_{k}"] = out.coef / out[f"se_{k}"]138 out["p_state_t50"] = 2 * (1 - stats.t.cdf(out.t_cluster_state.abs(), 50))139 out["p_zip3_t885"] = 2 * (1 - stats.t.cdf(out.t_cluster_zip3.abs(), len(np.unique(z3)) - 1))140 out["p_hc3_normal"] = 2 * (1 - stats.norm.cdf(out.t_hc3.abs()))141 slopes = out.index != "const"142 for k, pcol in [("state", "p_state_t50"), ("zip3", "p_zip3_t885"), ("hc3", "p_hc3_normal")]:143 out.loc[slopes, f"p_holm_{k}"] = holm(out.loc[slopes, pcol]); out.loc[slopes, f"p_bh_{k}"] = bh(out.loc[slopes, pcol])144 # CSS effective clusters (state)145 out["G_eff_css_state"] = css_effective_clusters(Xc.values, info["address_state"].values)146 out["G_eff_css_zip3"] = css_effective_clusters(Xc.values, info["zip3"].values)147 out["n_clusters_state"] = len(np.unique(st)); out["n_clusters_zip3"] = len(np.unique(z3))148 out.index.name = "variable"149 out.to_csv(OUT / "cluster_se_v2.csv")150 s = out[slopes]151 counts = []152 for k, G in [("hc3", None), ("cluster_zip3", len(np.unique(z3))), ("cluster_state", 51)]:153 t = s[f"t_{k}"].abs()154 crit = stats.t.ppf(0.975, G - 1) if G else 1.96155 pc = {"hc3": "p_hc3_normal", "cluster_zip3": "p_zip3_t885", "cluster_state": "p_state_t50"}[k]156 kk = k.replace("cluster_", "")157 counts.append({"se": k, "n_slopes": int(slopes.sum()), "n_abs_t_gt_196": int((t > 1.96).sum()),158 "critical_t": crit, "n_abs_t_gt_crit": int((t > crit).sum()),159 "n_holm_5pct": int((s[f"p_holm_{kk}"] < 0.05).sum()), "n_bh_5pct": int((s[f"p_bh_{kk}"] < 0.05).sum()),160 "median_G_eff_css": float(s["G_eff_css_state" if kk == "state" else "G_eff_css_zip3"].median()) if G else np.nan,161 "min_G_eff_css": float(s["G_eff_css_state" if kk == "state" else "G_eff_css_zip3"].min()) if G else np.nan})162 pd.DataFrame(counts).to_csv(OUT / "significance_counts_v2.csv", index=False)163 print(pd.DataFrame(counts).to_string(index=False))164 # '000' cluster: size and nearest-ZIP3 reassignment165 miss = info["zip3"].values == "000"166 cent = info.loc[~miss].groupby("zip3")[["latitude", "longitude"]].mean()167 tree = cKDTree(to_xy_km(cent.latitude.values, cent.longitude.values))168 _, j = tree.query(to_xy_km(info.loc[miss, "latitude"].values, info.loc[miss, "longitude"].values))169 z3b = info["zip3"].values.copy(); z3b[miss] = cent.index.values[j]170 r_z3b = base.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(z3b)[0]})171 cmp = pd.DataFrame({"se_cluster_zip3_with_000": r_z3.bse, "se_cluster_zip3_reassigned": r_z3b.bse}, index=names)172 cmp["rel_change"] = cmp.se_cluster_zip3_reassigned / cmp.se_cluster_zip3_with_000 - 1173 cmp.index.name = "variable"; cmp.to_csv(OUT / "zip3_missing_reassignment.csv")174 pd.DataFrame([{"n_listings_missing_zip": int(miss.sum()), "n_zip3_prefixes": int(len(cent)),175 "n_clusters_with_000": int(len(np.unique(z3))), "max_abs_rel_se_change_after_reassignment": float(cmp.rel_change.abs().max()),176 "reassigned_to": " ".join(sorted(set(cent.index.values[j])))}]).to_csv(OUT / "zip3_missing_cluster.csv", index=False)177 # Kennedy with clustered variance (natural units)178 Xu, yu, infou, _ = design.load_design("v2", trimmed=True, standardize=False)179 Xcu = sm.add_constant(Xu, has_constant="add")180 bu = sm.OLS(yu, Xcu.values)181 ru = {"hc3": bu.fit(cov_type="HC3"), "cluster_zip3": bu.fit(cov_type="cluster", cov_kwds={"groups": z3}),182 "cluster_state": bu.fit(cov_type="cluster", cov_kwds={"groups": st})}183 kr = []184 for v in DUMMIES:185 jv = list(Xcu.columns).index(v)186 for k, r in ru.items():187 b, se_ = r.params[jv], r.bse[jv]188 kr.append({"variable": v, "se": k, "coef": b, "se_value": se_, "naive_pct": 100 * (np.exp(b) - 1),189 "kennedy_pct": 100 * (np.exp(b - 0.5 * se_ ** 2) - 1), "t": b / se_})190 pd.DataFrame(kr).to_csv(OUT / "kennedy_clustered_v2.csv", index=False)191 unstd = pd.DataFrame({"coef": ru["hc3"].params, "se_hc3": ru["hc3"].bse, "se_cluster_zip3": ru["cluster_zip3"].bse,192 "se_cluster_state": ru["cluster_state"].bse}, index=list(Xcu.columns))193 unstd.index.name = "variable"; unstd.to_csv(OUT / "ols_v2_unstandardized_cluster_se.csv")194195196# ----------------------------------------------------------------------------- wild197def step_wild(B, workers):198 from wildboottest.wildboottest import wildboottest199 X, y, info, ex = design.load_design("v2", trimmed=True, standardize=True)200 Xc = sm.add_constant(X, has_constant="add")201 model = sm.OLS(y, Xc)202 cl = pd.Series(pd.factorize(info["address_state"].values)[0])203 rows = []204 t0 = time.time()205 for wt in ["rademacher", "webb"]:206 for v in [c for c in Xc.columns if c != "const"]:207 res = wildboottest(model, cluster=cl, B=B, param=v, weights_type=wt, seed=config.RANDOM_STATE, show=False)208 r = res.iloc[0].to_dict()209 rows.append({"variable": v, "weights": wt, **{k: r[k] for k in r if k in ("statistic", "p-value", "param")}, "B": B})210 print(f" wild bootstrap ({wt}) done in {time.time() - t0:.0f}s", flush=True)211 out = pd.DataFrame(rows).rename(columns={"p-value": "p_wild", "statistic": "t_cr1"})212 wide = out.pivot(index="variable", columns="weights", values="p_wild").rename(columns=lambda c: f"p_wild_{c}")213 tt = out[out.weights == "webb"].set_index("variable")["t_cr1"]214 wide["t_cr1"] = tt215 wide["p_holm_webb"] = holm(wide["p_wild_webb"]); wide["p_bh_webb"] = bh(wide["p_wild_webb"])216 wide["p_holm_rademacher"] = holm(wide["p_wild_rademacher"]); wide["p_bh_rademacher"] = bh(wide["p_wild_rademacher"])217 wide.to_csv(OUT / "wild_bootstrap_state.csv")218 cnt = {"n_slopes": len(wide), "n_sig_5pct_webb": int((wide.p_wild_webb < 0.05).sum()),219 "n_sig_5pct_rademacher": int((wide.p_wild_rademacher < 0.05).sum()),220 "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()),221 "n_sig_5pct_rademacher_holm": int((wide.p_holm_rademacher < 0.05).sum()), "B": B}222 pd.DataFrame([cnt]).to_csv(OUT / "wild_bootstrap_counts.csv", index=False)223 print(cnt)224225226# ----------------------------------------------------------------------------- states227def step_states():228 Xu, y, info, _ = design.load_design("v2", trimmed=True, standardize=False)229 cols = [c for c in Xu.columns if not c.startswith("region_")]230 states = info["address_state"].values; z3 = info["zip3"].values231 idx_tr, idx_te = train_test_split(np.arange(788_842), test_size=0.2, random_state=config.RANDOM_STATE)232 # map the full-sample split onto the trimmed sample via zpid233 full = pd.read_parquet(design.DESIGN_V2, columns=["zpid"])234 test_z = set(full["zpid"].values[idx_te])235 is_test = info["zpid"].isin(test_z).values236 rows, se_rows, dropped = [], [], []237 pred_inter = np.full(len(y), np.nan)238 for s_ in sorted(np.unique(states)):239 m = states == s_240 Xs = Xu.loc[m, cols]241 keep = [c for c in cols if Xs[c].nunique() > 1]242 dropped.append({"state": s_, "n": int(m.sum()), "dropped_constant_columns": " ".join(c for c in cols if c not in keep)})243 Xc = sm.add_constant(Xs[keep].values, has_constant="add")244 g = pd.factorize(z3[m])[0]245 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")246 cond = np.linalg.cond(Xc)247 row = {"state": s_, "n": int(m.sum()), "n_zip3": int(len(np.unique(g))), "r2": r.rsquared, "cond_number": cond,248 "median_price": float(np.exp(np.median(y[m])))}249 for v in REPORT13:250 if v in keep:251 j = keep.index(v) + 1252 row[v] = r.params[j]; row[f"se_{v}"] = r.bse[j]253 else:254 row[v] = np.nan; row[f"se_{v}"] = np.nan255 rows.append(row)256 tr = m & ~is_test; te = m & is_test257 if tr.sum() > len(keep) + 5 and te.sum() > 0:258 from sklearn.linear_model import LinearRegression259 pred_inter[te] = LinearRegression().fit(Xu.loc[tr, keep].values, y[tr]).predict(Xu.loc[te, keep].values)260 het = pd.DataFrame(rows); het.to_csv(OUT / "state_heterogeneity_v2.csv", index=False)261 pd.DataFrame(dropped).to_csv(OUT / "state_regressions_dropped_columns.csv", index=False)262 # pooled state-FE reference with ZIP3-cluster SEs263 sfe = pd.get_dummies(pd.Series(states), drop_first=True).astype(float).values264 Xp = np.column_stack([np.ones(len(y)), Xu[cols].values, sfe])265 pooled = sm.OLS(y, Xp).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(z3)[0]})266 summ = []267 for v in REPORT13:268 jp = cols.index(v) + 1269 b0, s0 = pooled.params[jp], pooled.bse[jp]270 b = het[v].values; se_ = het[f"se_{v}"].values; ok = ~np.isnan(b) & ~np.isnan(se_) & (se_ > 0)271 b, se_ = b[ok], se_[ok]272 w = 1 / se_ ** 2273 b_pw = np.sum(w * b) / np.sum(w)274 # cluster-robust homogeneity (states independent; ZIP3 clusters nested in states)275 Q = float(np.sum(w * (b - b_pw) ** 2)); dfq = len(b) - 1276 tau2 = max(0.0, (Q - dfq) / (np.sum(w) - np.sum(w ** 2) / np.sum(w))) # DerSimonian-Laird277 summ.append({"variable": v, "pooled_state_fe": b0, "pooled_se_zip3": s0, "n_states": int(len(b)),278 "precision_weighted_mean": b_pw, "listing_weighted_mean": np.average(het[v].values[ok], weights=het["n"].values[ok]),279 "p10": np.percentile(b, 10), "median": np.median(b), "p90": np.percentile(b, 90), "min": b.min(), "max": b.max(),280 "n_sign_differs_from_pooled": int((np.sign(b) != np.sign(b0)).sum()),281 "n_sign_differs_and_abs_t_gt_2": int(((np.sign(b) != np.sign(b0)) & (np.abs(b / se_) > 2)).sum()),282 "n_state_ci_excludes_pooled": int(((b - 1.96 * se_ > b0) | (b + 1.96 * se_ < b0)).sum()),283 "n_abs_t_gt_2": int((np.abs(b / se_) > 2).sum()),284 "homogeneity_Q": Q, "homogeneity_df": dfq, "homogeneity_p": float(1 - stats.chi2.cdf(Q, dfq)),285 "tau2_dl": tau2, "tau_dl": np.sqrt(tau2), "I2_pct": max(0.0, 100 * (Q - dfq) / Q) if Q > 0 else 0.0})286 summ = pd.DataFrame(summ); summ.to_csv(OUT / "state_heterogeneity_summary_v2.csv", index=False)287 print(summ[["variable", "pooled_state_fe", "p10", "median", "p90", "n_sign_differs_from_pooled", "n_sign_differs_and_abs_t_gt_2",288 "n_state_ci_excludes_pooled", "homogeneity_p", "tau_dl"]].round(3).to_string(index=False))289 # empirical-Bayes shrinkage for Figure 3 (normal-normal, tau2 DL)290 eb = []291 for v in REPORT13:292 t2 = float(summ.set_index("variable").loc[v, "tau2_dl"]); mu = float(summ.set_index("variable").loc[v, "precision_weighted_mean"])293 for _, r in het.iterrows():294 b, s_ = r[v], r[f"se_{v}"]295 if np.isnan(b) or np.isnan(s_) or s_ <= 0:296 continue297 shrink = t2 / (t2 + s_ ** 2)298 post = shrink * b + (1 - shrink) * mu; pv = shrink * s_ ** 2299 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_,300 "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),301 "shrinkage_factor": shrink})302 pd.DataFrame(eb).to_csv(OUT / "state_heterogeneity_eb.csv", index=False)303 # pooled vs state FE vs fully interacted on the random test split (trimmed)304 from sklearn.linear_model import LinearRegression305 te = is_test; tr = ~is_test306 p_pooled = LinearRegression().fit(Xu.loc[tr, cols].values, y[tr]).predict(Xu.loc[te, cols].values)307 gtr = pd.Series(states[tr]); Xd = Xu.loc[tr, cols].values - pd.DataFrame(Xu.loc[tr, cols].values).groupby(gtr.values).transform("mean").values308 yd = y[tr] - pd.Series(y[tr]).groupby(gtr.values).transform("mean").values309 beta, *_ = np.linalg.lstsq(Xd, yd, rcond=None)310 alpha = pd.Series(y[tr] - Xu.loc[tr, cols].values @ beta).groupby(gtr.values).mean()311 p_fe = Xu.loc[te, cols].values @ beta + pd.Series(states[te]).map(alpha).fillna(alpha.mean()).values312 ok = ~np.isnan(pred_inter[te])313 comp = pd.DataFrame([{"model": "Pooled OLS, no geographic controls", "r2_test": r2_score(y[te], p_pooled)},314 {"model": "Pooled OLS + 50 state fixed effects", "r2_test": r2_score(y[te], p_fe)},315 {"model": "State-by-state OLS (fully interacted)", "r2_test": r2_score(y[te][ok], pred_inter[te][ok])}])316 comp["n_train"] = int(tr.sum()); comp["n_test"] = int(te.sum())317 comp.to_csv(OUT / "pooled_vs_interacted_v2.csv", index=False); print(comp.to_string(index=False))318319320# ----------------------------------------------------------------------------- lot321def step_lot():322 rows = []323 for trimmed in [False, True]:324 Xl, y, info, _ = design.load_design("legacy", trimmed=trimmed, standardize=False, lot_indicator=False)325 Xv, _, _, _ = design.load_design("v2", trimmed=trimmed, standardize=False, lot_indicator=False)326 Xvi, _, _, _ = design.load_design("v2", trimmed=trimmed, standardize=False, lot_indicator=True)327 lm = info["lot_missing"].values == 1328 def fit(X, mask, label):329 Xc = sm.add_constant(X.values[mask], has_constant="add")330 r = sm.OLS(y[mask], Xc).fit(cov_type="HC3")331 names = list(X.columns)332 return {"sample": "trimmed (783,872)" if trimmed else "all (788,842)", "specification": label, "N": int(mask.sum()),333 "R2": r.rsquared, "beta_ln_lot": r.params[names.index("ln_lot") + 1], "se_ln_lot": r.bse[names.index("ln_lot") + 1],334 "beta_ln_sqft": r.params[names.index("ln_sqft") + 1],335 "beta_lot_missing": r.params[names.index("lot_missing") + 1] if "lot_missing" in names else np.nan}336 allm = np.ones(len(y), bool)337 rows.append(fit(Xl, allm, "v1.1 coding: ln(lot) = 0 when missing, zero, 1 sq ft or negative"))338 rows.append(fit(Xv, allm, "v2: ln(state-median lot) when missing/zero/<= 1"))339 rows.append(fit(Xvi, allm, "v2 + lot-missing indicator (the design used throughout)"))340 rows.append(fit(Xl, ~lm, "Complete cases (raw lot > 1 sq ft)"))341 p995 = info["lot_size_sqft"].quantile(0.995)342 Xw = Xvi.copy(); lot_v2 = np.exp(Xvi["ln_lot"].values); Xw["ln_lot"] = np.log(np.minimum(lot_v2, p995))343 rows.append(fit(Xw, allm, "v2 + indicator, lot winsorised at p99.5"))344 out = pd.DataFrame(rows); out.to_csv(OUT / "lot_recoding_v2.csv", index=False); print(out.round(4).to_string(index=False))345346347# ----------------------------------------------------------------------------- moran348def moran_knn(e, xy, k=8, permutations=999, seed=42):349 from esda.moran import Moran350 from libpysal.weights import WSP351 tree = cKDTree(xy)352 d, j = tree.query(xy, k=k + 1)353 n = len(xy)354 rows_ = np.repeat(np.arange(n), k); cols_ = j[:, 1:].ravel()355 W = sparse.csr_matrix((np.ones(n * k) / k, (rows_, cols_)), shape=(n, n))356 w = WSP(W).to_W()357 np.random.seed(seed)358 mi = Moran(e, w, transformation="r", permutations=permutations)359 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,360 "median_kth_neighbour_km": float(np.median(d[:, k]))}361362363def moran_band(e, pairs_i, pairs_j, n):364 """Moran's I for a symmetric binary band weight matrix given as pair lists; normal-approx z (randomisation)."""365 z = e - e.mean()366 S0 = 2 * len(pairs_i)367 num = 2 * np.sum(z[pairs_i] * z[pairs_j])368 I = (n / S0) * num / np.sum(z ** 2)369 EI = -1 / (n - 1)370 deg = np.bincount(np.concatenate([pairs_i, pairs_j]), minlength=n).astype(float)371 S1 = 0.5 * np.sum((2.0) ** 2) * len(pairs_i) # each symmetric pair contributes (w_ij + w_ji)^2 = 4, times 1/2372 S2 = np.sum((2 * deg) ** 2)373 b2 = n * np.sum(z ** 4) / np.sum(z ** 2) ** 2374 VI = (n * ((n ** 2 - 3 * n + 3) * S1 - n * S2 + 3 * S0 ** 2) - b2 * ((n ** 2 - n) * S1 - 2 * n * S2 + 6 * S0 ** 2)) / \375 ((n - 1) * (n - 2) * (n - 3) * S0 ** 2) - EI ** 2376 zI = (I - EI) / np.sqrt(VI)377 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)),378 "share_units_with_neighbour": float((deg > 0).mean())}379380381def step_moran(permutations):382 rows = []383 sets = []384 # v1.2 stored predictions (random test set, legacy design)385 try:386 ml = data.load_result("ml_results.pkl"); pk = pd.read_pickle(config.ANALYTICAL_SAMPLE)387 idx = ml["X_test"].index.values388 xy = to_xy_km(pk.loc[idx, "latitude"].values, pk.loc[idx, "longitude"].values)389 y = ml["y_test"].values390 for name, p in [("OLS", ml["y_pred_ols"]), ("LightGBM", ml["y_pred_lgb"]), ("XGBoost", ml["y_pred_xgb"])]:391 sets.append(("v1.2 stored predictions, random test set", name, y - p, xy))392 except FileNotFoundError:393 pass394 # v2 predictions from batch N2 when available395 prim = config.RESULTS_DIR / "primary"396 d = pd.read_parquet(design.DESIGN_V2, columns=["zpid", "latitude", "longitude"]).set_index("zpid")397 for name, pat in [("OLS", "predictions_random_OLS_62_regressors_+_lot_indicator.parquet"),398 ("LightGBM", "predictions_random_LightGBM_1000_rounds_early_stopping_seed42.parquet"),399 ("XGBoost", "predictions_random_XGBoost_1000_rounds_early_stopping_seed42.parquet")]:400 f = prim / pat401 if f.exists():402 p = pd.read_parquet(f)403 xy = to_xy_km(d.loc[p.zpid, "latitude"].values, d.loc[p.zpid, "longitude"].values)404 sets.append(("v2 predictions (N2), random test set, trimmed", name, p["y"].values - p["pred"].values, xy))405 bands = [(0, 1), (1, 5), (5, 25)]406 for label, name, e, xy in sets:407 t0 = time.time()408 r = moran_knn(e, xy, k=8, permutations=permutations)409 rows.append({"predictions": label, "model": name, "weights": "kNN k=8 row-standardised", "n": len(e), **r})410 print(f" {label} {name} kNN I={r['I']:.4f} z_sim={r['z_sim']:.1f} ({time.time() - t0:.0f}s)", flush=True)411 tree = cKDTree(xy)412 prev = None413 for lo, hi in bands:414 pairs = tree.query_pairs(hi, output_type="ndarray")415 if lo > 0:416 dd = np.linalg.norm(xy[pairs[:, 0]] - xy[pairs[:, 1]], axis=1)417 pairs = pairs[dd > lo]418 rb = moran_band(e, pairs[:, 0], pairs[:, 1], len(e))419 rows.append({"predictions": label, "model": name, "weights": f"distance band {lo}-{hi} km (binary)", "n": len(e), **rb})420 print(f" band {lo}-{hi} km: I={rb['I']:.4f} z={rb['z_rand']:.1f} pairs={rb['n_pairs']:,}", flush=True)421 out = pd.DataFrame(rows); out.to_csv(OUT / "moran_full_test_v2.csv", index=False)422423424# ----------------------------------------------------------------------------- conley425def step_conley(n_sub=100_000, cutoffs=(25, 50, 100)):426 X, y, info, ex = design.load_design("v2", trimmed=True, standardize=True)427 rng = np.random.RandomState(config.RANDOM_STATE)428 idx = np.sort(rng.choice(len(y), n_sub, replace=False))429 Xc = sm.add_constant(X.iloc[idx], has_constant="add").values; yy = y[idx]430 names = ["const"] + list(X.columns)431 r = sm.OLS(yy, Xc).fit(cov_type="HC0")432 e = r.resid433 A = np.linalg.inv(Xc.T @ Xc)434 xy = to_xy_km(info["latitude"].values[idx], info["longitude"].values[idx])435 tree = cKDTree(xy)436 out = pd.DataFrame({"coef": r.params, "se_hc0": r.bse}, index=names)437 st = pd.factorize(info["address_state"].values[idx])[0]438 out["se_cluster_state"] = sm.OLS(yy, Xc).fit(cov_type="cluster", cov_kwds={"groups": st}).bse439 Xe = Xc * e[:, None]440 for c in cutoffs:441 t0 = time.time()442 pairs = tree.query_pairs(c, output_type="ndarray")443 dd = np.linalg.norm(xy[pairs[:, 0]] - xy[pairs[:, 1]], axis=1)444 k = 1 - dd / c # Bartlett kernel445 n = len(yy)446 K = sparse.coo_matrix((np.concatenate([k, k]), (np.concatenate([pairs[:, 0], pairs[:, 1]]), np.concatenate([pairs[:, 1], pairs[:, 0]]))), shape=(n, n)).tocsr()447 M = Xe.T @ Xe + Xe.T @ (K @ Xe)448 V = A @ M @ A449 out[f"se_conley_{c}km"] = np.sqrt(np.diag(V))450 print(f" Conley {c} km: {len(pairs):,} pairs ({time.time() - t0:.0f}s)", flush=True)451 for c in cutoffs:452 out[f"ratio_conley_{c}km_hc0"] = out[f"se_conley_{c}km"] / out["se_hc0"]453 out["n_subsample"] = n_sub; out.index.name = "variable"454 out.to_csv(OUT / "conley_hac_v2.csv")455 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())456457458# ----------------------------------------------------------------------------- paired459def step_paired():460 prim = config.RESULTS_DIR / "primary"461 files = {p.name: p for p in prim.glob("predictions_random_*.parquet")}462 if not files:463 print(" no v2 random-split predictions yet (batch N2) — using stored v1.2 predictions")464 ml = data.load_result("ml_results.pkl"); pk = pd.read_pickle(config.ANALYTICAL_SAMPLE)465 z3 = pk.loc[ml["X_test"].index, "address_zipcode"].fillna("000").astype(str).str[:3].values466 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"467 else:468 d = pd.read_parquet(design.DESIGN_V2, columns=["zpid", "zip3"]).set_index("zpid")469 preds = {}470 for name, pat in [("OLS", "predictions_random_OLS_62_regressors_+_lot_indicator.parquet"),471 ("LightGBM", "predictions_random_LightGBM_1000_rounds_early_stopping_seed42.parquet"),472 ("XGBoost", "predictions_random_XGBoost_1000_rounds_early_stopping_seed42.parquet")]:473 p = pd.read_parquet(prim / pat); preds[name] = p["pred"].values; y = p["y"].values; zp = p["zpid"].values474 z3 = d.loc[zp, "zip3"].values; label = "v2 (N2)"475 g = pd.factorize(z3)[0]476 rows = []477 for a, c in [("XGBoost", "OLS"), ("XGBoost", "LightGBM"), ("LightGBM", "OLS")]:478 dl = (y - preds[c]) ** 2 - (y - preds[a]) ** 2479 r = sm.OLS(dl, np.ones(len(dl))).fit(cov_type="cluster", cov_kwds={"groups": g})480 r0 = sm.OLS(dl, np.ones(len(dl))).fit(cov_type="HC0")481 rows.append({"predictions": label, "quantity": f"{c} minus {a} squared-error difference", "mean": float(dl.mean()),482 "se_iid": float(dl.std(ddof=1) / np.sqrt(len(dl))), "se_hc0": float(r0.bse[0]), "se_cluster_zip3": float(r.bse[0]),483 "t_cluster_zip3": float(r.params[0] / r.bse[0]), "n_clusters": int(len(np.unique(g)))})484 pd.DataFrame(rows).to_csv(OUT / "paired_loss_cluster_se.csv", index=False); print(pd.DataFrame(rows).round(5).to_string(index=False))485486487def main():488 ap = argparse.ArgumentParser(description=__doc__)489 ap.add_argument("--steps", nargs="+", default=["assert", "cluster", "wild", "states", "lot", "moran", "conley", "paired"])490 ap.add_argument("--B", type=int, default=9999)491 ap.add_argument("--workers", type=int, default=8)492 ap.add_argument("--permutations", type=int, default=999)493 args = ap.parse_args()494 OUT.mkdir(parents=True, exist_ok=True)495 warnings.filterwarnings("ignore")496 import os497 os.environ.setdefault("OMP_NUM_THREADS", str(args.workers))498 t0 = time.time()499 for s in args.steps:500 print(f"== {s} ({time.time() - t0:.0f}s)", flush=True)501 {"assert": step_assert, "cluster": step_cluster, "wild": lambda: step_wild(args.B, args.workers), "states": step_states,502 "lot": step_lot, "moran": lambda: step_moran(args.permutations), "conley": step_conley, "paired": step_paired}[s]()503 (OUT / "run_meta.json").write_text(json.dumps({"steps": args.steps, "B": args.B, "elapsed_s": round(time.time() - t0)}, indent=1))504 print(f"Done in {time.time() - t0:.0f}s")505506507if __name__ == "__main__":508 main()509