#!/usr/bin/env python3 # Author: Simon-Pierre Boucher — contact@spboucher.ai # """Extended geographic-validation analyses added in the October 2026 revision. All analyses start from the stored analytical sample (data/processed/ analytical_sample.pkl) and use the XGBoost configuration of src/wp3/models.py. Each XGBoost fit takes 6-8 s on an 18-core Apple Silicon machine, so the whole script runs in roughly 10-15 minutes. Produces results/revision/ with: geo_trees_curve.csv held-out-state R2 and early-stopping-set R2 as a function of the number of boosting rounds size_matched_split.csv random split with the geographic training-set size variants_two_designs.csv full model / no-school / no-school-no-region / missing-indicator variants under both designs, plus within-state R2 for the 10-state holdout loso_by_state.csv leave-one-state-out R2, RMSE, within-state R2 and MdAPE for each of the 51 jurisdictions grouped_cv_by_state.csv 5-fold grouped cross-validation (folds = sets of states), pooled and per-fold metrics ols_geo_holdout.csv OLS under the 10-state holdout, raw and after removing implausible attribute values Usage: python scripts/05_revision_geo_validation.py [--quick] """ import argparse import sys import time from pathlib import Path import numpy as np import pandas as pd from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_squared_error, r2_score from sklearn.model_selection import GroupKFold, train_test_split from xgboost import XGBRegressor sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) from wp3 import config, data # noqa: E402 from wp3.models import XGB_PARAMS # noqa: E402 OUT = config.RESULTS_DIR / "revision" SCHOOL_FEATS = ["avg_school_rating", "school_count", "nearest_school_distance"] def fit_xgb(X_tr, y_tr, n_estimators=None): params = dict(XGB_PARAMS) if n_estimators is not None: params["n_estimators"] = n_estimators X_fit, X_eval, y_fit, y_eval = train_test_split( X_tr, y_tr, test_size=0.1, random_state=config.RANDOM_STATE) model = XGBRegressor(**params) model.fit(X_fit, y_fit, eval_set=[(X_eval, y_eval)], verbose=0) return model, (X_eval, y_eval) def metrics(y, yhat, groups=None): out = { "r2": r2_score(y, yhat), "rmse": float(np.sqrt(mean_squared_error(y, yhat))), "mdape": float(np.median(np.abs(np.exp(yhat) / np.exp(y) - 1)) * 100), } if groups is not None: # within-group R2: demean actual and predicted log prices by group g = pd.Series(groups) yd = y - pd.Series(y).groupby(g.values).transform("mean").values pd_ = yhat - pd.Series(yhat).groupby(g.values).transform("mean").values out["r2_within"] = r2_score(yd, pd_) # share of the test-set variance that is between groups out["between_share"] = 1 - yd.var() / y.var() return out def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--quick", action="store_true", help="skip the leave-one-state-out loop") args = parser.parse_args() OUT.mkdir(parents=True, exist_ok=True) t0 = time.time() df = data.load_analytical_sample() X_full, parts = data.build_feature_matrix(df) cols = list(X_full.columns) X = X_full.values.astype(np.float32) y = df["ln_price"].values states = df["address_state"].values geo_mask = df["address_state"].isin(config.HOLDOUT_STATES).values idx_train, idx_test = train_test_split( np.arange(len(df)), test_size=0.2, random_state=config.RANDOM_STATE) # ---------------------------------------------------------------- 1. trees curve print("1. Held-out-state R2 as a function of boosting rounds") model, (X_eval, y_eval) = fit_xgb(X[~geo_mask], y[~geo_mask]) rows = [] for k in [25, 50, 75, 100, 150, 200, 300, 400, 500, 600, 700, 800, 900, 1000]: rows.append({ "n_trees": k, "r2_geo": r2_score(y[geo_mask], model.predict(X[geo_mask], iteration_range=(0, k))), "r2_eval_within_training_states": r2_score( y_eval, model.predict(X_eval, iteration_range=(0, k))), }) curve = pd.DataFrame(rows) # random-split model, same grid, for comparison model_r, (Xe_r, ye_r) = fit_xgb(X[idx_train], y[idx_train]) curve["r2_random_test"] = [ r2_score(y[idx_test], model_r.predict(X[idx_test], iteration_range=(0, k))) for k in curve["n_trees"]] curve.to_csv(OUT / "geo_trees_curve.csv", index=False) print(curve.to_string(index=False)) # ---------------------------------------------------------------- 2. size-matched split print("\n2. Size-matched random split") rng = np.random.RandomState(config.RANDOM_STATE) perm = rng.permutation(len(df)) n_geo_train = int((~geo_mask).sum()) i_tr, i_te = perm[:n_geo_train], perm[n_geo_train:] m_sm, _ = fit_xgb(X[i_tr], y[i_tr]) sm = pd.DataFrame([ {"design": "random 80/20 (631,073 train)", "n_train": len(idx_train), "n_test": len(idx_test), **metrics(y[idx_test], model_r.predict(X[idx_test]))}, {"design": "random, geographic training size", "n_train": len(i_tr), "n_test": len(i_te), **metrics(y[i_te], m_sm.predict(X[i_te]))}, {"design": "10-state geographic holdout", "n_train": n_geo_train, "n_test": int(geo_mask.sum()), **metrics(y[geo_mask], model.predict(X[geo_mask]), groups=states[geo_mask])}, ]) sm.to_csv(OUT / "size_matched_split.csv", index=False) print(sm.to_string(index=False)) # ---------------------------------------------------------------- 3. feature variants print("\n3. Feature variants under both designs") region_cols = parts["region_dummies"] yb_missing = (df["age"].round(6).isin( df.groupby("address_state")["age"].transform("median").round(6)) & (df["age"] > 0)).values # proxy not used; replaced below by exact flags # exact missingness flags come from the raw database when available flags = None try: import duckdb con = duckdb.connect(str(config.DUCKDB_PATH), read_only=True) flags = con.execute(""" SELECT zpid, CASE WHEN year_built IS NULL THEN 1 ELSE 0 END AS miss_year_built, CASE WHEN lot_size_sqft IS NULL THEN 1 ELSE 0 END AS miss_lot, CASE WHEN avg_school_rating IS NULL THEN 1 ELSE 0 END AS miss_school, CASE WHEN walk_score IS NULL THEN 1 ELSE 0 END AS miss_walk, CASE WHEN has_garage IS NULL THEN 1 ELSE 0 END AS miss_garage, CASE WHEN transit_score IS NULL THEN 1 ELSE 0 END AS miss_transit FROM properties""").fetchdf() con.close() flags = df[["zpid"]].merge(flags, on="zpid", how="left").fillna(0) except Exception as exc: # pragma: no cover print(" (raw database unavailable, missing-indicator variant skipped):", exc) variants = { "Full model (62 features)": cols, "Without school variables (59)": [c for c in cols if c not in SCHOOL_FEATS], "Without region dummies (59)": [c for c in cols if c not in region_cols], "Without school variables and region dummies (56)": [c for c in cols if c not in SCHOOL_FEATS and c not in region_cols], } frames = {name: X_full[c].values.astype(np.float32) for name, c in variants.items()} if flags is not None: ind = flags.drop(columns="zpid").values.astype(np.float32) frames["Full model + 6 missingness indicators (68)"] = np.hstack([X, ind]) rows = [] for name, Xv in frames.items(): m_rand, _ = fit_xgb(Xv[idx_train], y[idx_train]) m_geo, _ = fit_xgb(Xv[~geo_mask], y[~geo_mask]) r = metrics(y[idx_test], m_rand.predict(Xv[idx_test])) g = metrics(y[geo_mask], m_geo.predict(Xv[geo_mask]), groups=states[geo_mask]) rows.append({"variant": name, "n_features": Xv.shape[1], "r2_random": r["r2"], "mdape_random": r["mdape"], "r2_geo": g["r2"], "r2_geo_within_state": g["r2_within"], "mdape_geo": g["mdape"], "r2_geo_100_trees": r2_score( y[geo_mask], m_geo.predict(Xv[geo_mask], iteration_range=(0, 100)))}) print(f" {name:<55} random={r['r2']:.4f} geo={g['r2']:.4f} " f"within={g['r2_within']:.4f} geo@100={rows[-1]['r2_geo_100_trees']:.4f}") pd.DataFrame(rows).to_csv(OUT / "variants_two_designs.csv", index=False) # ---------------------------------------------------------------- 4. OLS geo holdout print("\n4. OLS under the 10-state holdout") def ols_geo(mask_keep, label): tr = (~geo_mask) & mask_keep te = geo_mask & mask_keep lr = LinearRegression(n_jobs=-1).fit(X_full.values[tr], y[tr]) pred = lr.predict(X_full.values[te]) m = metrics(y[te], pred, groups=states[te]) return {"specification": label, "n_train": int(tr.sum()), "n_test": int(te.sum()), **m} implausible = ((df["parking_spaces"] > 20) | (df["stories"] > 10) | (df["lot_size_sqft"] > 43560 * 640) | (df["age"] >= 200)).values nogeo_cols = [c for c in cols if c not in region_cols] ols_rows = [ols_geo(np.ones(len(df), bool), "OLS, 62 regressors, all observations"), ols_geo(~implausible, "OLS, 62 regressors, implausible values removed")] tr = ~geo_mask & ~implausible te = geo_mask & ~implausible lr = LinearRegression(n_jobs=-1).fit(X_full[nogeo_cols].values[tr], y[tr]) ols_rows.append({"specification": "OLS, no region dummies, implausible values removed", "n_train": int(tr.sum()), "n_test": int(te.sum()), **metrics(y[te], lr.predict(X_full[nogeo_cols].values[te]), groups=states[te])}) ols_df = pd.DataFrame(ols_rows) ols_df["n_implausible_removed"] = int(implausible.sum()) ols_df.to_csv(OUT / "ols_geo_holdout.csv", index=False) print(ols_df.to_string(index=False)) # ---------------------------------------------------------------- 5. grouped 5-fold CV print("\n5. Grouped 5-fold cross-validation by state (full model and no-school model)") gkf = GroupKFold(n_splits=5) cv_rows = [] for name in ["Full model (62 features)", "Without school variables (59)"]: Xv = frames[name] pooled_pred = np.empty(len(df)) for k, (tr, te) in enumerate(gkf.split(Xv, y, groups=states)): m, _ = fit_xgb(Xv[tr], y[tr]) pred = m.predict(Xv[te]) pooled_pred[te] = pred mm = metrics(y[te], pred, groups=states[te]) cv_rows.append({"variant": name, "fold": k + 1, "n_test": len(te), "n_states": len(np.unique(states[te])), **mm}) print(f" {name[:30]:<30} fold {k + 1}: r2={mm['r2']:.4f} within={mm['r2_within']:.4f}") mm = metrics(y, pooled_pred, groups=states) cv_rows.append({"variant": name, "fold": "pooled", "n_test": len(df), "n_states": 51, **mm}) pd.DataFrame(cv_rows).to_csv(OUT / "grouped_cv_by_state.csv", index=False) # ---------------------------------------------------------------- 6. leave-one-state-out if not args.quick: print("\n6. Leave-one-state-out (51 fits per variant)") loso_rows = [] for name in ["Full model (62 features)", "Without school variables (59)"]: Xv = frames[name] for st in sorted(np.unique(states)): te = states == st m, _ = fit_xgb(Xv[~te], y[~te]) pred = m.predict(Xv[te]) mm = metrics(y[te], pred) mm["mean_error"] = float(np.mean(pred - y[te])) # level bias in log points mm["r2_100_trees"] = r2_score(y[te], m.predict(Xv[te], iteration_range=(0, 100))) loso_rows.append({"variant": name, "state": st, "n_test": int(te.sum()), **mm}) print(f" {name[:12]:<12} {st}: n={te.sum():>7,} r2={mm['r2']:.3f} " f"bias={mm['mean_error']:+.3f} mdape={mm['mdape']:.1f}") pd.DataFrame(loso_rows).to_csv(OUT / "loso_by_state.csv", index=False) print(f"\nDone in {time.time() - t0:.0f}s; outputs in {OUT}") if __name__ == "__main__": main()