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 geographic designs from the out-of-fold predictions of505_geo_designs.py (November 2026 review, batch L2; items WP03-ADVE-17, EDIT-01,6THEO-03, THEO-01, THEO-02, ECON-05, EDIT-14, THEO-12, ADVE-22, THEO-14, EDIT-13,7THEO-15, ADVE-03, ECON-06, EDIT-06).89From results/geo/predictions/<design>_<model>_<variant>_seed<k>.parquet:10 gap ZIP3-cluster bootstrap (1,000 draws, seed 42) of R2(OLS) - R2(XGB nested),11 - R2(XGB 1,000 rounds), - R2(XGB 100 rounds), - R2(LightGBM nested) for every12 design (pooled out-of-fold R2); seed dispersion (mean, SD, min, max over seeds)13 loso per-state paired differences (51 states), sign test and Wilcoxon; LOSO14 summaries (median, listing-weighted mean, p10-p90, states > 0.5)15 decomposition SSE within / between (level) shares and oracle within-group R2 per design16 feasible feasible level correction: level estimated from k in {5, 20, 50, 100, 200}17 random listings of the held-out state, 100 repetitions (LOSO and ten-state)18 two_by_two region in/out x rounds {nested, 100, 1,000} under ten-state and LOSO19 duplicates random-split R2 without the test rows that have an exact duplicate in train20 range the OLS-vs-XGB gap as a range across designs and seeds (cluster A headline)2122Outputs: results/geo_inference/*.csv2324Usage: python scripts/09_geo_inference.py [--draws 1000]25"""2627import argparse28import sys29import time30from pathlib import Path3132import numpy as np33import pandas as pd34from scipy import stats35from sklearn.metrics import r2_score3637sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))38from wp3 import config, design # noqa: E4023940GEO = config.RESULTS_DIR / "geo"41PRED = GEO / "predictions"42OUT = config.RESULTS_DIR / "geo_inference"43DESIGNS = ["random", "size_matched", "ten_state", "gkf_state", "gkf_zip3", "gkf_zip5", "loso"]44COMPARATORS = [("xgb", "full", "pred_nested"), ("xgb", "full", "pred_fixed1000"), ("xgb", "full", "pred_fixed100"),45 ("xgb", "plus_missing", "pred_nested"), ("lgb", "full", "pred_nested")]464748def load_pred(dsg, model, variant, seed):49 f = PRED / f"{dsg}_{model}_{variant}_seed{seed}.parquet"50 return pd.read_parquet(f) if f.exists() else None515253def seeds_available(dsg, model, variant):54 return sorted(int(p.stem.split("seed")[-1]) for p in PRED.glob(f"{dsg}_{model}_{variant}_seed*.parquet"))555657def within_r2(y, p, g):58 gs = pd.Series(g)59 yd = y - pd.Series(y).groupby(gs.values).transform("mean").values60 pdm = p - pd.Series(p).groupby(gs.values).transform("mean").values61 return r2_score(yd, pdm)626364def step_gap(info, draws):65 rng = np.random.RandomState(config.RANDOM_STATE)66 rows = []67 for dsg in DESIGNS:68 ols = load_pred(dsg, "ols", "full", 42)69 if ols is None:70 continue71 ols = ols.set_index("zpid")72 y = ols["y"].values73 z3 = info.loc[ols.index, "zip3"].values74 codes, uniq = pd.factorize(z3)75 members = [np.where(codes == g)[0] for g in range(len(uniq))]76 boot_idx = [np.concatenate([members[g] for g in rng.randint(0, len(uniq), len(uniq))]) for _ in range(draws)]77 base = {"design": dsg, "n_test": len(y), "n_zip3": len(uniq), "r2_ols": r2_score(y, ols["pred"].values)}78 for model, variant, col in COMPARATORS:79 seeds = seeds_available(dsg, model, variant)80 if not seeds:81 continue82 r2s, gaps = [], []83 for s in seeds:84 p = load_pred(dsg, model, variant, s).set_index("zpid").loc[ols.index]85 r2s.append(r2_score(y, p[col].values)); gaps.append(base["r2_ols"] - r2s[-1])86 # ZIP3-cluster bootstrap of the gap for the first seed (paired: same resampled listings)87 p = load_pred(dsg, model, variant, seeds[0]).set_index("zpid").loc[ols.index]88 po, pm = ols["pred"].values, p[col].values89 bg = np.array([r2_score(y[b], po[b]) - r2_score(y[b], pm[b]) for b in boot_idx])90 rows.append({**base, "comparator": f"{model}_{variant}_{col.replace('pred_', '')}", "n_seeds": len(seeds),91 "r2_comparator_mean": np.mean(r2s), "r2_comparator_sd": np.std(r2s, ddof=1) if len(r2s) > 1 else np.nan,92 "r2_comparator_min": np.min(r2s), "r2_comparator_max": np.max(r2s),93 "gap_ols_minus_comparator_mean": np.mean(gaps), "gap_sd_over_seeds": np.std(gaps, ddof=1) if len(gaps) > 1 else np.nan,94 "gap_min_over_seeds": np.min(gaps), "gap_max_over_seeds": np.max(gaps),95 "gap_seed42": gaps[0], "gap_boot_se_zip3": bg.std(ddof=1), "gap_ci_low_zip3": np.percentile(bg, 2.5),96 "gap_ci_high_zip3": np.percentile(bg, 97.5), "n_draws": draws})97 print(f" {dsg:12s} {rows[-1]['comparator']:28s} OLS {base['r2_ols']:.4f} comp {np.mean(r2s):.4f} "98 f"gap {np.mean(gaps):+.4f} [{rows[-1]['gap_ci_low_zip3']:+.4f}, {rows[-1]['gap_ci_high_zip3']:+.4f}]", flush=True)99 out = pd.DataFrame(rows); out.to_csv(OUT / "gap_bootstrap.csv", index=False)100 return out101102103def step_loso(info):104 ols = load_pred("loso", "ols", "full", 42)105 if ols is None:106 return107 ols = ols.set_index("zpid"); y = ols["y"].values; st = info.loc[ols.index, "address_state"].values108 per = []109 for s_ in sorted(np.unique(st)):110 m = st == s_111 row = {"state": s_, "n_test": int(m.sum()), "r2_ols": r2_score(y[m], ols["pred"].values[m]),112 "r2_within_ols": r2_score(y[m] - y[m].mean(), ols["pred"].values[m] - ols["pred"].values[m].mean()),113 "mean_error_ols": float(np.mean(ols["pred"].values[m] - y[m]))}114 for model, variant, col in COMPARATORS:115 seeds = seeds_available("loso", model, variant)116 if not seeds:117 continue118 vals, within, me = [], [], []119 for s in seeds:120 p = load_pred("loso", model, variant, s).set_index("zpid").loc[ols.index][col].values121 vals.append(r2_score(y[m], p[m])); within.append(r2_score(y[m] - y[m].mean(), p[m] - p[m].mean())); me.append(np.mean(p[m] - y[m]))122 k = f"{model}_{variant}_{col.replace('pred_', '')}"123 row[f"r2_{k}"] = np.mean(vals); row[f"r2_{k}_sd"] = np.std(vals, ddof=1) if len(vals) > 1 else np.nan124 row[f"r2_within_{k}"] = np.mean(within); row[f"mean_error_{k}"] = np.mean(me)125 row[f"diff_ols_minus_{k}"] = row["r2_ols"] - np.mean(vals)126 per.append(row)127 per = pd.DataFrame(per); per.to_csv(OUT / "loso_per_state.csv", index=False)128 tests, summ = [], []129 cols = [c for c in per.columns if c.startswith("r2_") and not c.endswith("_sd") and "within" not in c]130 for c in cols:131 v = per[c].values; w = per["n_test"].values132 summ.append({"model": c.replace("r2_", ""), "median": np.median(v), "listing_weighted_mean": np.average(v, weights=w),133 "unweighted_mean": v.mean(), "p10": np.percentile(v, 10), "p90": np.percentile(v, 90), "min": v.min(), "max": v.max(),134 "n_states_gt_05": int((v > 0.5).sum()), "n_negative": int((v < 0).sum()),135 "median_abs_mean_error": float(np.median(np.abs(per[c.replace("r2_", "mean_error_")].values))) if c.replace("r2_", "mean_error_") in per else np.nan})136 if c != "r2_ols":137 d = per["r2_ols"].values - v138 nz = d[d != 0]139 tests.append({"comparison": f"OLS minus {c.replace('r2_', '')}", "n_states": len(d), "mean_diff": d.mean(), "median_diff": np.median(d),140 "n_states_ols_better": int((d > 0).sum()), "sign_test_p": float(stats.binomtest(int((nz > 0).sum()), len(nz)).pvalue),141 "wilcoxon_p": float(stats.wilcoxon(d).pvalue), "paired_t_p": float(stats.ttest_1samp(d, 0).pvalue),142 "listing_weighted_mean_diff": float(np.average(d, weights=w))})143 pd.DataFrame(summ).to_csv(OUT / "loso_summaries.csv", index=False)144 pd.DataFrame(tests).to_csv(OUT / "loso_paired_tests.csv", index=False)145 print(pd.DataFrame(summ).round(3).to_string(index=False)); print(pd.DataFrame(tests).round(4).to_string(index=False))146147148def step_decomposition(info):149 rows = []150 for dsg in DESIGNS:151 for model, variant, col in [("ols", "full", "pred"), ("ols", "no_region", "pred"), ("ols_state_fe", "no_region", "pred")] + COMPARATORS:152 seeds = seeds_available(dsg, model, variant)153 if not seeds:154 continue155 p = load_pred(dsg, model, variant, seeds[0])156 y = p["y"].values; pr = p[col].values157 st = info.loc[p["zpid"].values, "address_state"].values158 groups = {"state": st}159 if dsg == "gkf_zip3":160 groups["zip3"] = info.loc[p["zpid"].values, "zip3"].values161 if dsg == "gkf_zip5":162 groups["zip5"] = info.loc[p["zpid"].values, "zip5"].values163 row = {"design": dsg, "model": f"{model}_{variant}_{col.replace('pred_', '').replace('pred', 'ols')}", "seed": seeds[0], "n": len(y), "r2": r2_score(y, pr)}164 e = pr - y; sst = ((y - y.mean()) ** 2).sum()165 for gname, g in groups.items():166 eb = pd.Series(e).groupby(pd.Series(g).values).transform("mean").values167 yd = y - pd.Series(y).groupby(pd.Series(g).values).transform("mean").values168 row[f"r2_within_{gname}_oracle"] = within_r2(y, pr, g)169 row[f"sse_within_share_{gname}"] = ((e - eb) ** 2).sum() / sst170 row[f"sse_between_share_{gname}"] = (eb ** 2).sum() / sst171 row[f"between_share_of_variance_{gname}"] = 1 - yd.var() / y.var()172 rows.append(row)173 out = pd.DataFrame(rows); out.to_csv(OUT / "sse_decomposition.csv", index=False)174 print(out[["design", "model", "r2", "r2_within_state_oracle", "sse_within_share_state", "sse_between_share_state"]].round(4).to_string(index=False))175176177def step_feasible(info, ks=(5, 20, 50, 100, 200), reps=100):178 rng = np.random.RandomState(config.RANDOM_STATE)179 rows = []180 for dsg in ["loso", "ten_state"]:181 for model, variant, col in [("ols", "full", "pred"), ("xgb", "full", "pred_nested"), ("xgb", "full", "pred_fixed1000")]:182 seeds = seeds_available(dsg, model, variant)183 if not seeds:184 continue185 p = load_pred(dsg, model, variant, seeds[0]); y = p["y"].values; pr = p[col].values186 st = info.loc[p["zpid"].values, "address_state"].values187 states = sorted(np.unique(st))188 for k in ks:189 r2_pooled, r2_states = [], {s_: [] for s_ in states}190 for rep in range(reps):191 adj = pr.copy()192 for s_ in states:193 m = np.where(st == s_)[0]194 if len(m) <= k:195 continue196 cal = rng.choice(m, k, replace=False)197 level = np.mean(y[cal] - pr[cal])198 adj[m] = pr[m] + level199 rest = np.setdiff1d(m, cal)200 r2_states[s_].append(r2_score(y[rest], adj[rest]))201 r2_pooled.append(r2_score(y, adj))202 # oracle = all listings used for the level203 adj_o = pr + (pd.Series(y - pr).groupby(st).transform("mean").values)204 rows.append({"design": dsg, "model": f"{model}_{variant}_{col.replace('pred_', '').replace('pred', 'ols')}", "k_calibration": k, "reps": reps,205 "r2_raw": r2_score(y, pr), "r2_oracle_level": r2_score(y, adj_o),206 "r2_feasible_median": np.median(r2_pooled), "r2_feasible_p25": np.percentile(r2_pooled, 25), "r2_feasible_p75": np.percentile(r2_pooled, 75),207 "r2_feasible_mean": np.mean(r2_pooled),208 "median_state_r2_feasible": np.median([np.median(v) for v in r2_states.values() if v])})209 print(f" {dsg} {rows[-1]['model']} k={k}: raw {rows[-1]['r2_raw']:.4f} oracle {rows[-1]['r2_oracle_level']:.4f} feasible {rows[-1]['r2_feasible_median']:.4f}", flush=True)210 pd.DataFrame(rows).to_csv(OUT / "feasible_level_correction.csv", index=False)211212213def step_two_by_two(info):214 met = pd.read_csv(GEO / "metrics.csv")215 pooled = pd.read_csv(GEO / "metrics_pooled.csv") if (GEO / "metrics_pooled.csv").exists() else None216 rows = []217 for dsg in ["ten_state", "loso", "gkf_state"]:218 for variant in ["full", "no_region", "no_school", "no_school_no_region", "plus_missing"]:219 for rule in ["nested", "fixed100", "fixed1000"]:220 if dsg == "ten_state":221 sub = met[(met.design == dsg) & (met.model == "xgb") & (met.variant == variant) & (met.rounds_rule == rule)]222 if sub.empty:223 continue224 rows.append({"design": dsg, "variant": variant, "rounds_rule": rule, "region_dummies": "no_region" not in variant,225 "school_vars": "no_school" not in variant, "r2_mean": sub.r2.mean(), "r2_sd": sub.r2.std(ddof=1) if len(sub) > 1 else np.nan,226 "n_seeds": len(sub), "r2_within_state_oracle": sub.r2_within_state.mean(), "mdape": sub.mdape.mean(),227 "selected_rounds": sub.selected_rounds.iloc[0] if "selected_rounds" in sub else np.nan})228 else:229 sub = pooled[(pooled.design == dsg) & (pooled.model == "xgb") & (pooled.variant == variant) & (pooled.rounds_rule == rule)] if pooled is not None else pd.DataFrame()230 if sub.empty:231 continue232 per = met[(met.design == dsg) & (met.model == "xgb") & (met.variant == variant) & (met.rounds_rule == rule)]233 rows.append({"design": dsg, "variant": variant, "rounds_rule": rule, "region_dummies": "no_region" not in variant,234 "school_vars": "no_school" not in variant, "r2_mean": sub.r2.mean(), "r2_sd": sub.r2.std(ddof=1) if len(sub) > 1 else np.nan,235 "n_seeds": len(sub), "r2_within_state_oracle": sub.r2_within_state.mean(), "mdape": sub.mdape.mean(),236 "median_fold_r2": per.groupby("seed").r2.median().mean()})237 out = pd.DataFrame(rows); out.to_csv(OUT / "two_by_two_region_rounds.csv", index=False)238 print(out.round(4).to_string(index=False))239240241def step_duplicates(info):242 dup = info["dup_group_id"].values243 rows = []244 for model, variant, col in [("ols", "full", "pred")] + COMPARATORS[:2]:245 p = load_pred("random", model, variant, 42)246 if p is None:247 continue248 y = p["y"].values; pr = p[col].values249 d = info.loc[p["zpid"].values, "dup_group_id"].values250 # a test row has a duplicate in train iff its group has a member not in the test set251 test_z = set(p["zpid"].values)252 grp_all = info[info["dup_group_id"] >= 0]253 in_train = grp_all[~grp_all.index.isin(test_z)].groupby("dup_group_id").size()254 has_train_dup = (d >= 0) & pd.Series(d).map(in_train).fillna(0).values.astype(bool)255 rows.append({"model": f"{model}_{variant}_{col.replace('pred_', '').replace('pred', 'ols')}", "n_test": len(y),256 "n_test_with_train_duplicate": int(has_train_dup.sum()), "r2_all": r2_score(y, pr),257 "r2_without_duplicates": r2_score(y[~has_train_dup], pr[~has_train_dup]),258 "r2_duplicates_only": r2_score(y[has_train_dup], pr[has_train_dup]) if has_train_dup.sum() > 10 else np.nan})259 pd.DataFrame(rows).to_csv(OUT / "duplicates_random_split_v2.csv", index=False); print(pd.DataFrame(rows).round(4).to_string(index=False))260261262def step_range(gap):263 if gap is None or gap.empty:264 return265 sub = gap[gap.comparator.isin(["xgb_full_nested", "xgb_full_fixed1000", "xgb_full_fixed100", "lgb_full_nested"])]266 rows = []267 for comp, g in sub.groupby("comparator"):268 rows.append({"comparator": comp, "designs": " ".join(g.design), "gap_min_design": g.gap_ols_minus_comparator_mean.min(),269 "gap_max_design": g.gap_ols_minus_comparator_mean.max(),270 "design_of_min": g.loc[g.gap_ols_minus_comparator_mean.idxmin(), "design"],271 "design_of_max": g.loc[g.gap_ols_minus_comparator_mean.idxmax(), "design"],272 "min_ci_low": g.gap_ci_low_zip3.min(), "max_ci_high": g.gap_ci_high_zip3.max(),273 "n_designs_ci_excludes_zero": int(((g.gap_ci_low_zip3 > 0) | (g.gap_ci_high_zip3 < 0)).sum()), "n_designs": len(g)})274 pd.DataFrame(rows).to_csv(OUT / "gap_range.csv", index=False); print(pd.DataFrame(rows).round(4).to_string(index=False))275276277def main():278 ap = argparse.ArgumentParser(description=__doc__)279 ap.add_argument("--draws", type=int, default=1000)280 ap.add_argument("--steps", nargs="+", default=["gap", "loso", "decomposition", "feasible", "two_by_two", "duplicates", "range"])281 args = ap.parse_args()282 OUT.mkdir(parents=True, exist_ok=True)283 t0 = time.time()284 info = pd.read_parquet(design.DESIGN_V2, columns=["zpid", "address_state", "zip3", "zip5", "dup_group_id"]).set_index("zpid")285 gap = None286 for s in args.steps:287 print(f"== {s} ({time.time() - t0:.0f}s)", flush=True)288 if s == "gap":289 gap = step_gap(info, args.draws)290 elif s == "range":291 step_range(gap if gap is not None else (pd.read_csv(OUT / "gap_bootstrap.csv") if (OUT / "gap_bootstrap.csv").exists() else None))292 else:293 {"loso": step_loso, "decomposition": step_decomposition, "feasible": step_feasible, "two_by_two": step_two_by_two,294 "duplicates": step_duplicates}[s](info)295 print(f"Done in {time.time() - t0:.0f}s")296297298if __name__ == "__main__":299 main()300