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"""Extended geographic-validation analyses added in the October 2026 revision.56All analyses start from the stored analytical sample (data/processed/7analytical_sample.pkl) and use the XGBoost configuration of src/wp3/models.py.8Each XGBoost fit takes 6-8 s on an 18-core Apple Silicon machine, so the whole9script runs in roughly 10-15 minutes.1011Produces results/revision/ with:12 geo_trees_curve.csv held-out-state R2 and early-stopping-set R2 as a13 function of the number of boosting rounds14 size_matched_split.csv random split with the geographic training-set size15 variants_two_designs.csv full model / no-school / no-school-no-region /16 missing-indicator variants under both designs,17 plus within-state R2 for the 10-state holdout18 loso_by_state.csv leave-one-state-out R2, RMSE, within-state R2 and19 MdAPE for each of the 51 jurisdictions20 grouped_cv_by_state.csv 5-fold grouped cross-validation (folds = sets of21 states), pooled and per-fold metrics22 ols_geo_holdout.csv OLS under the 10-state holdout, raw and after23 removing implausible attribute values2425Usage: python scripts/05_revision_geo_validation.py [--quick]26"""2728import argparse29import sys30import time31from pathlib import Path3233import numpy as np34import pandas as pd35from sklearn.linear_model import LinearRegression36from sklearn.metrics import mean_squared_error, r2_score37from sklearn.model_selection import GroupKFold, train_test_split38from xgboost import XGBRegressor3940sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))41from wp3 import config, data # noqa: E40242from wp3.models import XGB_PARAMS # noqa: E4024344OUT = config.RESULTS_DIR / "revision"45SCHOOL_FEATS = ["avg_school_rating", "school_count", "nearest_school_distance"]464748def fit_xgb(X_tr, y_tr, n_estimators=None):49 params = dict(XGB_PARAMS)50 if n_estimators is not None:51 params["n_estimators"] = n_estimators52 X_fit, X_eval, y_fit, y_eval = train_test_split(53 X_tr, y_tr, test_size=0.1, random_state=config.RANDOM_STATE)54 model = XGBRegressor(**params)55 model.fit(X_fit, y_fit, eval_set=[(X_eval, y_eval)], verbose=0)56 return model, (X_eval, y_eval)575859def metrics(y, yhat, groups=None):60 out = {61 "r2": r2_score(y, yhat),62 "rmse": float(np.sqrt(mean_squared_error(y, yhat))),63 "mdape": float(np.median(np.abs(np.exp(yhat) / np.exp(y) - 1)) * 100),64 }65 if groups is not None:66 # within-group R2: demean actual and predicted log prices by group67 g = pd.Series(groups)68 yd = y - pd.Series(y).groupby(g.values).transform("mean").values69 pd_ = yhat - pd.Series(yhat).groupby(g.values).transform("mean").values70 out["r2_within"] = r2_score(yd, pd_)71 # share of the test-set variance that is between groups72 out["between_share"] = 1 - yd.var() / y.var()73 return out747576def main():77 parser = argparse.ArgumentParser(description=__doc__)78 parser.add_argument("--quick", action="store_true",79 help="skip the leave-one-state-out loop")80 args = parser.parse_args()81 OUT.mkdir(parents=True, exist_ok=True)82 t0 = time.time()8384 df = data.load_analytical_sample()85 X_full, parts = data.build_feature_matrix(df)86 cols = list(X_full.columns)87 X = X_full.values.astype(np.float32)88 y = df["ln_price"].values89 states = df["address_state"].values90 geo_mask = df["address_state"].isin(config.HOLDOUT_STATES).values91 idx_train, idx_test = train_test_split(92 np.arange(len(df)), test_size=0.2, random_state=config.RANDOM_STATE)9394 # ---------------------------------------------------------------- 1. trees curve95 print("1. Held-out-state R2 as a function of boosting rounds")96 model, (X_eval, y_eval) = fit_xgb(X[~geo_mask], y[~geo_mask])97 rows = []98 for k in [25, 50, 75, 100, 150, 200, 300, 400, 500, 600, 700, 800, 900, 1000]:99 rows.append({100 "n_trees": k,101 "r2_geo": r2_score(y[geo_mask], model.predict(X[geo_mask], iteration_range=(0, k))),102 "r2_eval_within_training_states": r2_score(103 y_eval, model.predict(X_eval, iteration_range=(0, k))),104 })105 curve = pd.DataFrame(rows)106 # random-split model, same grid, for comparison107 model_r, (Xe_r, ye_r) = fit_xgb(X[idx_train], y[idx_train])108 curve["r2_random_test"] = [109 r2_score(y[idx_test], model_r.predict(X[idx_test], iteration_range=(0, k)))110 for k in curve["n_trees"]]111 curve.to_csv(OUT / "geo_trees_curve.csv", index=False)112 print(curve.to_string(index=False))113114 # ---------------------------------------------------------------- 2. size-matched split115 print("\n2. Size-matched random split")116 rng = np.random.RandomState(config.RANDOM_STATE)117 perm = rng.permutation(len(df))118 n_geo_train = int((~geo_mask).sum())119 i_tr, i_te = perm[:n_geo_train], perm[n_geo_train:]120 m_sm, _ = fit_xgb(X[i_tr], y[i_tr])121 sm = pd.DataFrame([122 {"design": "random 80/20 (631,073 train)", "n_train": len(idx_train),123 "n_test": len(idx_test), **metrics(y[idx_test], model_r.predict(X[idx_test]))},124 {"design": "random, geographic training size", "n_train": len(i_tr),125 "n_test": len(i_te), **metrics(y[i_te], m_sm.predict(X[i_te]))},126 {"design": "10-state geographic holdout", "n_train": n_geo_train,127 "n_test": int(geo_mask.sum()),128 **metrics(y[geo_mask], model.predict(X[geo_mask]), groups=states[geo_mask])},129 ])130 sm.to_csv(OUT / "size_matched_split.csv", index=False)131 print(sm.to_string(index=False))132133 # ---------------------------------------------------------------- 3. feature variants134 print("\n3. Feature variants under both designs")135 region_cols = parts["region_dummies"]136 yb_missing = (df["age"].round(6).isin(137 df.groupby("address_state")["age"].transform("median").round(6)) &138 (df["age"] > 0)).values # proxy not used; replaced below by exact flags139 # exact missingness flags come from the raw database when available140 flags = None141 try:142 import duckdb143 con = duckdb.connect(str(config.DUCKDB_PATH), read_only=True)144 flags = con.execute("""145 SELECT zpid,146 CASE WHEN year_built IS NULL THEN 1 ELSE 0 END AS miss_year_built,147 CASE WHEN lot_size_sqft IS NULL THEN 1 ELSE 0 END AS miss_lot,148 CASE WHEN avg_school_rating IS NULL THEN 1 ELSE 0 END AS miss_school,149 CASE WHEN walk_score IS NULL THEN 1 ELSE 0 END AS miss_walk,150 CASE WHEN has_garage IS NULL THEN 1 ELSE 0 END AS miss_garage,151 CASE WHEN transit_score IS NULL THEN 1 ELSE 0 END AS miss_transit152 FROM properties""").fetchdf()153 con.close()154 flags = df[["zpid"]].merge(flags, on="zpid", how="left").fillna(0)155 except Exception as exc: # pragma: no cover156 print(" (raw database unavailable, missing-indicator variant skipped):", exc)157158 variants = {159 "Full model (62 features)": cols,160 "Without school variables (59)": [c for c in cols if c not in SCHOOL_FEATS],161 "Without region dummies (59)": [c for c in cols if c not in region_cols],162 "Without school variables and region dummies (56)":163 [c for c in cols if c not in SCHOOL_FEATS and c not in region_cols],164 }165 frames = {name: X_full[c].values.astype(np.float32) for name, c in variants.items()}166 if flags is not None:167 ind = flags.drop(columns="zpid").values.astype(np.float32)168 frames["Full model + 6 missingness indicators (68)"] = np.hstack([X, ind])169 rows = []170 for name, Xv in frames.items():171 m_rand, _ = fit_xgb(Xv[idx_train], y[idx_train])172 m_geo, _ = fit_xgb(Xv[~geo_mask], y[~geo_mask])173 r = metrics(y[idx_test], m_rand.predict(Xv[idx_test]))174 g = metrics(y[geo_mask], m_geo.predict(Xv[geo_mask]), groups=states[geo_mask])175 rows.append({"variant": name, "n_features": Xv.shape[1],176 "r2_random": r["r2"], "mdape_random": r["mdape"],177 "r2_geo": g["r2"], "r2_geo_within_state": g["r2_within"],178 "mdape_geo": g["mdape"],179 "r2_geo_100_trees": r2_score(180 y[geo_mask], m_geo.predict(Xv[geo_mask], iteration_range=(0, 100)))})181 print(f" {name:<55} random={r['r2']:.4f} geo={g['r2']:.4f} "182 f"within={g['r2_within']:.4f} geo@100={rows[-1]['r2_geo_100_trees']:.4f}")183 pd.DataFrame(rows).to_csv(OUT / "variants_two_designs.csv", index=False)184185 # ---------------------------------------------------------------- 4. OLS geo holdout186 print("\n4. OLS under the 10-state holdout")187 def ols_geo(mask_keep, label):188 tr = (~geo_mask) & mask_keep189 te = geo_mask & mask_keep190 lr = LinearRegression(n_jobs=-1).fit(X_full.values[tr], y[tr])191 pred = lr.predict(X_full.values[te])192 m = metrics(y[te], pred, groups=states[te])193 return {"specification": label, "n_train": int(tr.sum()), "n_test": int(te.sum()), **m}194 implausible = ((df["parking_spaces"] > 20) | (df["stories"] > 10) |195 (df["lot_size_sqft"] > 43560 * 640) | (df["age"] >= 200)).values196 nogeo_cols = [c for c in cols if c not in region_cols]197 ols_rows = [ols_geo(np.ones(len(df), bool), "OLS, 62 regressors, all observations"),198 ols_geo(~implausible, "OLS, 62 regressors, implausible values removed")]199 tr = ~geo_mask & ~implausible200 te = geo_mask & ~implausible201 lr = LinearRegression(n_jobs=-1).fit(X_full[nogeo_cols].values[tr], y[tr])202 ols_rows.append({"specification": "OLS, no region dummies, implausible values removed",203 "n_train": int(tr.sum()), "n_test": int(te.sum()),204 **metrics(y[te], lr.predict(X_full[nogeo_cols].values[te]), groups=states[te])})205 ols_df = pd.DataFrame(ols_rows)206 ols_df["n_implausible_removed"] = int(implausible.sum())207 ols_df.to_csv(OUT / "ols_geo_holdout.csv", index=False)208 print(ols_df.to_string(index=False))209210 # ---------------------------------------------------------------- 5. grouped 5-fold CV211 print("\n5. Grouped 5-fold cross-validation by state (full model and no-school model)")212 gkf = GroupKFold(n_splits=5)213 cv_rows = []214 for name in ["Full model (62 features)", "Without school variables (59)"]:215 Xv = frames[name]216 pooled_pred = np.empty(len(df))217 for k, (tr, te) in enumerate(gkf.split(Xv, y, groups=states)):218 m, _ = fit_xgb(Xv[tr], y[tr])219 pred = m.predict(Xv[te])220 pooled_pred[te] = pred221 mm = metrics(y[te], pred, groups=states[te])222 cv_rows.append({"variant": name, "fold": k + 1, "n_test": len(te),223 "n_states": len(np.unique(states[te])), **mm})224 print(f" {name[:30]:<30} fold {k + 1}: r2={mm['r2']:.4f} within={mm['r2_within']:.4f}")225 mm = metrics(y, pooled_pred, groups=states)226 cv_rows.append({"variant": name, "fold": "pooled", "n_test": len(df), "n_states": 51, **mm})227 pd.DataFrame(cv_rows).to_csv(OUT / "grouped_cv_by_state.csv", index=False)228229 # ---------------------------------------------------------------- 6. leave-one-state-out230 if not args.quick:231 print("\n6. Leave-one-state-out (51 fits per variant)")232 loso_rows = []233 for name in ["Full model (62 features)", "Without school variables (59)"]:234 Xv = frames[name]235 for st in sorted(np.unique(states)):236 te = states == st237 m, _ = fit_xgb(Xv[~te], y[~te])238 pred = m.predict(Xv[te])239 mm = metrics(y[te], pred)240 mm["mean_error"] = float(np.mean(pred - y[te])) # level bias in log points241 mm["r2_100_trees"] = r2_score(y[te], m.predict(Xv[te], iteration_range=(0, 100)))242 loso_rows.append({"variant": name, "state": st, "n_test": int(te.sum()), **mm})243 print(f" {name[:12]:<12} {st}: n={te.sum():>7,} r2={mm['r2']:.3f} "244 f"bias={mm['mean_error']:+.3f} mdape={mm['mdape']:.1f}")245 pd.DataFrame(loso_rows).to_csv(OUT / "loso_by_state.csv", index=False)246247 print(f"\nDone in {time.time() - t0:.0f}s; outputs in {OUT}")248249250if __name__ == "__main__":251 main()252