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, data-quality and heterogeneity analyses added in the October 20265revision. Everything starts from the stored artifacts (model_data.pkl,6analytical_sample.pkl, ml_results.pkl) and the raw DuckDB extract (for the7original missingness flags). Runtime: about 5 minutes.89Produces results/revision/ with:10 missingness_full.csv original missingness of every regressor input and11 the imputation rule actually applied12 implausible_values.csv counts of implausible attribute values13 cluster_se.csv OLS coefficients with HC3, state-clustered and14 ZIP3-clustered standard errors15 lot_recoding.csv OLS under alternative codings of lot size16 state_heterogeneity.csv state-by-state OLS coefficients (natural units)17 state_heterogeneity_summary.csv18 pooled_vs_interacted.csv random-split R2 of pooled, state-FE and19 state-by-state OLS20 r2_bootstrap.csv bootstrap CIs for test-set R2 and paired21 differences (observation and ZIP3-cluster bootstrap)22 moran_xgb_vs_ols.csv Moran's I on OLS and XGBoost test residuals23 list_to_zestimate.csv distribution of list price / Zestimate2425Usage: python scripts/06_revision_inference.py26"""2728import sys29import time30from pathlib import Path3132import duckdb33import numpy as np34import pandas as pd35import statsmodels.api as sm36from sklearn.metrics import r2_score37from sklearn.model_selection import train_test_split3839sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))40from wp3 import config, data # noqa: E4024142OUT = config.RESULTS_DIR / "revision"43KEY = ["ln_sqft", "bedrooms", "bathrooms", "age", "age_sq", "ln_lot", "has_pool",44 "has_garage", "on_waterfront", "has_central_air", "has_hardwood", "luxury_score",45 "walk_score", "bike_score", "transit_score", "avg_school_rating",46 "property_tax_rate", "is_condo", "has_hoa", "tag_new_construction",47 "tag_foreclosure", "region_Northeast", "region_South", "region_West"]4849RAW_INPUTS = { # regressor input -> raw column, imputation rule actually applied50 "year_built (age, age^2, interactions)": ("year_built", "state median year built"),51 "lot_size_sqft (ln lot)": ("lot_size_sqft", "state median in lot_size_sqft; ln(lot) coded 0"),52 "avg_school_rating": ("avg_school_rating", "state median"),53 "school_count": ("school_count", "state median (= 3 in every state)"),54 "nearest_school_distance": ("nearest_school_distance", "state median"),55 "walk_score": ("walk_score", "state median"),56 "bike_score": ("bike_score", "state median"),57 "transit_score": ("transit_score", "set to 0"),58 "has_garage": ("has_garage", "set to False"),59 "parking_spaces": ("parking_spaces", "set to 0"),60 "stories": ("stories", "set to 1"),61 "property_tax_rate": ("property_tax_rate", "state median"),62 "hoa_fee_value (ln HOA)": ("hoa_fee_value", "set to 0"),63 "has_pool / has_spa / has_basement / has_fireplace / on_waterfront":64 ("has_pool", "set to False"),65 "roof_type / construction_material / foundation_type":66 ("roof_type", "'Unknown' category"),67}686970def missingness(df):71 con = duckdb.connect(str(config.DUCKDB_PATH), read_only=True)72 cols = ["zpid", "year_built", "lot_size_sqft", "avg_school_rating", "school_count",73 "nearest_school_distance", "walk_score", "bike_score", "transit_score",74 "has_garage", "parking_spaces", "stories", "property_tax_rate", "hoa_fee_value",75 "has_pool", "has_spa", "has_basement", "has_fireplace", "on_waterfront",76 "roof_type", "construction_material", "foundation_type", "zestimate",77 "cooling_type", "heating_type", "flooring"]78 raw = con.execute(f"SELECT {', '.join(cols)} FROM properties").fetchdf()79 con.close()80 m = df[["zpid", "address_state", "price", "ln_lot", "lot_size_sqft", "age",81 "parking_spaces", "stories", "has_pool", "has_spa", "has_basement",82 "has_fireplace", "on_waterfront", "has_central_air", "has_forced_air",83 "has_hardwood", "roof_cat", "construction_cat", "foundation_cat",84 "avg_school_rating", "region"]].merge(raw, on="zpid", how="left",85 suffixes=("", "_raw"))86 rows = []87 for label, (col, rule) in RAW_INPUTS.items():88 src = col + "_raw" if col + "_raw" in m.columns else col89 rows.append({"input": label, "raw_column": col,90 "missing_share": float(m[src].isna().mean()), "imputation": rule})91 for col in ["cooling_type", "heating_type", "flooring"]:92 rows.append({"input": f"{col} (central air / forced air / hardwood dummies)",93 "raw_column": col, "missing_share": float(m[col].isna().mean()),94 "imputation": "dummy set to False"})95 miss = pd.DataFrame(rows).sort_values("missing_share", ascending=False)96 miss.to_csv(OUT / "missingness_full.csv", index=False)9798 imp = {99 "parking_spaces > 20": int((df["parking_spaces"] > 20).sum()),100 "parking_spaces > 1,000": int((df["parking_spaces"] > 1000).sum()),101 "stories > 10": int((df["stories"] > 10).sum()),102 "lot_size > 1 sq. mile (640 acres)": int((df["lot_size_sqft"] > 43560 * 640).sum()),103 "lot_size > 100 acres": int((df["lot_size_sqft"] > 43560 * 100).sum()),104 "lot_size recorded as 0 (raw)": int((m["lot_size_sqft_raw"] == 0).sum()),105 "ln(lot) coded 0 (missing or zero lot)": int((df["ln_lot"] == 0).sum()),106 "year_built > 2026 (age set to 0)": int((m["year_built"] > 2026).sum()),107 "year_built < 1826 (age capped at 200)": int((m["year_built"] < 1826).sum()),108 "age == 0": int((df["age"] == 0).sum()),109 "bike_score > 100": int((m["bike_score"] > 100).sum()),110 "school rating imputed (state median)": int(m["avg_school_rating_raw"].isna().sum()),111 }112 pd.Series(imp, name="n").rename_axis("check").to_csv(OUT / "implausible_values.csv")113114 z = m["price"] / m["zestimate"]115 z = z[m["zestimate"].notna() & (z > 0.2) & (z < 5)]116 q = z.quantile([0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95])117 out = pd.DataFrame({"all": q})118 for reg in ["Northeast", "Midwest", "South", "West"]:119 zz = z[m.loc[z.index, "region"] == reg]120 out[reg] = zz.quantile(q.index)121 out.loc["n"] = [len(z)] + [int((m.loc[z.index, "region"] == r).sum())122 for r in ["Northeast", "Midwest", "South", "West"]]123 out.loc["share_above_1"] = [float((z > 1).mean())] + [124 float((z[m.loc[z.index, "region"] == r] > 1).mean())125 for r in ["Northeast", "Midwest", "South", "West"]]126 out.to_csv(OUT / "list_to_zestimate.csv")127 return m128129130def cluster_se(md, df):131 X, y = md["X_const"], md["y_clean"]132 zip3 = df["address_zipcode"].fillna("000").astype(str).str[:3].values133 state = df["address_state"].fillna("NA").values134 base = sm.OLS(y.values, X.values)135 res_hc3 = base.fit(cov_type="HC3")136 res_st = base.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(state)[0]})137 res_z3 = base.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(zip3)[0]})138 cols = list(X.columns)139 out = pd.DataFrame({140 "coef": res_hc3.params, "se_hc3": res_hc3.bse, "se_cluster_state": res_st.bse,141 "se_cluster_zip3": res_z3.bse}, index=cols)142 out["ratio_zip3_hc3"] = out["se_cluster_zip3"] / out["se_hc3"]143 out["ratio_state_hc3"] = out["se_cluster_state"] / out["se_hc3"]144 out["t_hc3"] = out["coef"] / out["se_hc3"]145 out["t_cluster_state"] = out["coef"] / out["se_cluster_state"]146 out["t_cluster_zip3"] = out["coef"] / out["se_cluster_zip3"]147 out["n_clusters_state"] = len(np.unique(state))148 out["n_clusters_zip3"] = len(np.unique(zip3))149 out.index.name = "variable"150 out.to_csv(OUT / "cluster_se.csv")151 print(out.loc[[k for k in KEY if k in out.index],152 ["coef", "se_hc3", "se_cluster_zip3", "se_cluster_state",153 "ratio_zip3_hc3", "ratio_state_hc3"]].round(4).to_string())154 return out155156157def unstandardized_matrix(df_clean, X_const):158 """Rebuild the 62-regressor design matrix in natural units (as script 02)."""159 df_m = df_clean.copy()160 xcols = [c for c in X_const.columns if c != "const"]161 for cat_col in config.CATEGORICAL_COLS:162 dummies = pd.get_dummies(df_m[cat_col], prefix=cat_col)163 for dc in dummies.columns:164 if dc in xcols:165 df_m[dc] = dummies[dc].astype(float)166 return df_m, xcols167168169def lot_recoding(df_m, xcols, raw_m):170 y = df_m["ln_price"].values171 lot_missing = raw_m["lot_size_sqft_raw"].isna().values172 lot_zero = (raw_m["lot_size_sqft_raw"] == 0).values173 state_median = df_m.groupby("address_state")["lot_size_sqft"].transform("median")174 lot_imp = df_m["lot_size_sqft"].where(~(lot_missing | lot_zero), state_median)175 ln_lot_imp = np.log(lot_imp.clip(lower=1)).values176 ind = (lot_missing | lot_zero).astype(float)177178 def fit(Xmat, names, mask, label):179 Xc = sm.add_constant(Xmat[mask])180 r = sm.OLS(y[mask], Xc).fit(cov_type="HC3")181 j = names.index("ln_lot") + 1182 k = names.index("ln_sqft") + 1183 return {"specification": label, "N": int(mask.sum()), "R2": r.rsquared,184 "beta_ln_lot": r.params[j], "se_ln_lot": r.bse[j],185 "beta_ln_sqft": r.params[k],186 "beta_lot_missing_indicator": (r.params[names.index("lot_missing") + 1]187 if "lot_missing" in names else np.nan)}188189 X0 = df_m[xcols].values.astype(float)190 allm = np.ones(len(df_m), bool)191 X1 = X0.copy(); X1[:, xcols.index("ln_lot")] = ln_lot_imp192 X2 = np.column_stack([X1, ind]); names2 = xcols + ["lot_missing"]193 rows = [194 fit(X0, xcols, allm, "Baseline coding: ln(lot) = 0 when lot missing or zero (as published)"),195 fit(X1, xcols, allm, "ln(lot) = ln(state-median lot) when missing or zero"),196 fit(X2, names2, allm, "ln(lot) = ln(state-median lot) + missing-lot indicator"),197 fit(X0, xcols, ~(lot_missing | lot_zero), "Complete cases (lot reported and positive)"),198 ]199 # winsorized version of the recoded variable, for comparison with the published check200 p995 = df_m["lot_size_sqft"].quantile(0.995)201 X3 = X1.copy()202 X3[:, xcols.index("ln_lot")] = np.log(lot_imp.clip(lower=1, upper=p995)).values203 rows.append(fit(X3, xcols, allm, "Recoded + lot winsorized at p99.5"))204 out = pd.DataFrame(rows)205 out.to_csv(OUT / "lot_recoding.csv", index=False)206 print(out.round(4).to_string(index=False))207 return ln_lot_imp208209210def state_heterogeneity(df_m, xcols, ln_lot_imp):211 """State-by-state OLS in natural units (region dummies dropped: constant within state)."""212 y = df_m["ln_price"].values213 cols = [c for c in xcols if not c.startswith("region_")]214 Xdf = df_m[cols].astype(float).copy()215 Xdf["ln_lot"] = ln_lot_imp216 # trim implausible values so that state regressions are not driven by a few records217 ok = ((df_m["parking_spaces"] <= 20) & (df_m["stories"] <= 10) &218 (df_m["lot_size_sqft"] <= 43560 * 640) & (df_m["age"] < 200)).values219 states = df_m["address_state"].values220 report = ["ln_sqft", "bedrooms", "bathrooms", "age", "ln_lot", "has_pool",221 "has_garage", "has_central_air", "has_hardwood", "luxury_score",222 "tag_foreclosure", "is_condo", "has_hoa"]223 rows = []224 idx_train, idx_test = train_test_split(np.arange(len(df_m)), test_size=0.2,225 random_state=config.RANDOM_STATE)226 is_test = np.zeros(len(df_m), bool); is_test[idx_test] = True227 pred_inter = np.full(len(df_m), np.nan)228 for st in sorted(np.unique(states)):229 m_st = (states == st) & ok230 tr = m_st & ~is_test231 te = m_st & is_test232 Xc = sm.add_constant(Xdf.values[tr], has_constant="add")233 beta, *_ = np.linalg.lstsq(Xc, y[tr], rcond=None)234 pred_inter[te] = sm.add_constant(Xdf.values[te], has_constant="add") @ beta235 full = sm.OLS(y[m_st], sm.add_constant(Xdf.values[m_st], has_constant="add")).fit()236 row = {"state": st, "n": int(m_st.sum()), "r2": full.rsquared,237 "median_price": float(np.exp(np.median(y[m_st])))}238 for v in report:239 row[v] = full.params[cols.index(v) + 1]240 rows.append(row)241 het = pd.DataFrame(rows)242 het.to_csv(OUT / "state_heterogeneity.csv", index=False)243244 # pooled with state FE (same trimmed sample, natural units) for reference245 sfe = pd.get_dummies(df_m["address_state"], drop_first=True).astype(float)246 Xp = np.column_stack([Xdf.values, sfe.values])247 tr = ok & ~is_test; te = ok & is_test248 pooled_fe = sm.OLS(y[tr], sm.add_constant(Xp[tr], has_constant="add")).fit()249 pred_fe = sm.add_constant(Xp[te], has_constant="add") @ pooled_fe.params250 pooled = sm.OLS(y[tr], sm.add_constant(Xdf.values[tr], has_constant="add")).fit()251 pred_pooled = sm.add_constant(Xdf.values[te], has_constant="add") @ pooled.params252 comp = pd.DataFrame([253 {"model": "Pooled OLS, no geographic controls", "r2_test": r2_score(y[te], pred_pooled)},254 {"model": "Pooled OLS + 50 state fixed effects", "r2_test": r2_score(y[te], pred_fe)},255 {"model": "State-by-state OLS (fully interacted)", "r2_test": r2_score(y[te], pred_inter[te])},256 ])257 comp["n_train"] = int(tr.sum()); comp["n_test"] = int(te.sum())258 comp.to_csv(OUT / "pooled_vs_interacted.csv", index=False)259 print(comp.to_string(index=False))260261 summ = []262 w = het["n"].values263 for v in report:264 vals = het[v].values265 wm = np.average(vals, weights=w)266 summ.append({"variable": v, "pooled_state_fe": pooled_fe.params[cols.index(v) + 1],267 "weighted_mean": wm, "min": vals.min(), "p10": np.percentile(vals, 10),268 "median": np.median(vals), "p90": np.percentile(vals, 90), "max": vals.max(),269 "share_positive": float((vals > 0).mean()),270 "n_states_sign_differs_from_pooled":271 int((np.sign(vals) != np.sign(pooled_fe.params[cols.index(v) + 1])).sum())})272 summ = pd.DataFrame(summ)273 summ.to_csv(OUT / "state_heterogeneity_summary.csv", index=False)274 print(summ.round(4).to_string(index=False))275276277def xgb_geo_trimmed(df):278 """XGBoost (full model) under the 10-state holdout, scored on the same trimmed279 test set as the OLS comparison in script 05 (implausible records removed)."""280 from sklearn.model_selection import train_test_split as tts281 from xgboost import XGBRegressor282 from wp3.models import XGB_PARAMS283 X_full, _ = data.build_feature_matrix(df)284 X = X_full.values.astype(np.float32); y = df["ln_price"].values285 geo = df["address_state"].isin(config.HOLDOUT_STATES).values286 implausible = ((df["parking_spaces"] > 20) | (df["stories"] > 10) |287 (df["lot_size_sqft"] > 43560 * 640) | (df["age"] >= 200)).values288 X_fit, X_eval, y_fit, y_eval = tts(X[~geo], y[~geo], test_size=0.1,289 random_state=config.RANDOM_STATE)290 m = XGBRegressor(**XGB_PARAMS).fit(X_fit, y_fit, eval_set=[(X_eval, y_eval)], verbose=0)291 te = geo & ~implausible292 rows = []293 for k in [100, 1000]:294 pred = m.predict(X[te], iteration_range=(0, k))295 rows.append({"model": f"XGBoost, {k} trees, trimmed test set", "n_test": int(te.sum()),296 "r2": r2_score(y[te], pred),297 "rmse": float(np.sqrt(np.mean((y[te] - pred) ** 2))),298 "mdape": float(np.median(np.abs(np.exp(pred - y[te]) - 1)) * 100)})299 out = pd.DataFrame(rows)300 out.to_csv(OUT / "xgb_geo_trimmed.csv", index=False)301 print(out.round(4).to_string(index=False))302303304def r2_bootstrap(ml, df, n_boot=1000):305 y = ml["y_test"].values306 preds = {"OLS": ml["y_pred_ols"], "LightGBM": ml["y_pred_lgb"], "XGBoost": ml["y_pred_xgb"]}307 zip3 = df.loc[ml["X_test"].index, "address_zipcode"].fillna("000").astype(str).str[:3].values308 codes, uniq = pd.factorize(zip3)309 rng = np.random.RandomState(config.RANDOM_STATE)310 n = len(y)311312 def r2s(idx):313 yy = y[idx]314 return {k: r2_score(yy, v[idx]) for k, v in preds.items()}315316 obs, clu = [], []317 members = [np.where(codes == g)[0] for g in range(len(uniq))]318 for _ in range(n_boot):319 obs.append(r2s(rng.randint(0, n, n)))320 draw = rng.randint(0, len(uniq), len(uniq))321 clu.append(r2s(np.concatenate([members[g] for g in draw])))322 rows = []323 for label, boots in [("observation bootstrap", obs), ("ZIP3-cluster bootstrap", clu)]:324 b = pd.DataFrame(boots)325 for k in preds:326 rows.append({"scheme": label, "quantity": f"R2 {k}", "point": r2_score(y, preds[k]),327 "ci_low": b[k].quantile(0.025), "ci_high": b[k].quantile(0.975),328 "boot_se": b[k].std()})329 for a, c in [("XGBoost", "OLS"), ("XGBoost", "LightGBM"), ("LightGBM", "OLS")]:330 d = b[a] - b[c]331 rows.append({"scheme": label, "quantity": f"R2 {a} - R2 {c}",332 "point": r2_score(y, preds[a]) - r2_score(y, preds[c]),333 "ci_low": d.quantile(0.025), "ci_high": d.quantile(0.975),334 "boot_se": d.std()})335 # paired Diebold-Mariano-type test on squared errors (random split => iid losses)336 for a, c in [("XGBoost", "OLS"), ("XGBoost", "LightGBM")]:337 dl = (y - preds[c]) ** 2 - (y - preds[a]) ** 2338 rows.append({"scheme": "paired t on squared-error differences", "quantity": f"{c} minus {a} loss",339 "point": dl.mean(), "ci_low": np.nan, "ci_high": np.nan,340 "boot_se": dl.std() / np.sqrt(n)})341 out = pd.DataFrame(rows)342 out["n_boot"] = n_boot343 out.to_csv(OUT / "r2_bootstrap.csv", index=False)344 print(out.round(4).to_string(index=False))345346347def moran_xgb_vs_ols(ml, df, n_sub=3, size=5000, k=8):348 from esda.moran import Moran349 from libpysal.weights import KNN350 idx = ml["X_test"].index.values351 coords = df.loc[idx, ["latitude", "longitude"]].values352 y = ml["y_test"].values353 res = {"OLS": y - ml["y_pred_ols"], "LightGBM": y - ml["y_pred_lgb"],354 "XGBoost": y - ml["y_pred_xgb"]}355 rows = []356 for i in range(n_sub):357 rng = np.random.RandomState(config.RANDOM_STATE + i)358 sub = rng.choice(len(idx), size=size, replace=False)359 w = KNN.from_array(coords[sub], k=k); w.transform = "r"360 for name, e in res.items():361 mi = Moran(e[sub], w, permutations=999)362 rows.append({"subsample": i + 1, "model": name, "I": mi.I, "z": mi.z_sim,363 "p_sim": mi.p_sim})364 out = pd.DataFrame(rows)365 out.to_csv(OUT / "moran_xgb_vs_ols.csv", index=False)366 print(out.groupby("model")["I"].agg(["mean", "std"]).round(4))367368369def main():370 OUT.mkdir(parents=True, exist_ok=True)371 t0 = time.time()372 df = data.load_analytical_sample()373 md = data.load_model_data()374 ml = data.load_result("ml_results.pkl")375376 print("1. Missingness, implausible values, list/Zestimate")377 raw_m = missingness(df)378 print("\n2. Clustered standard errors")379 cluster_se(md, df)380 print("\n3. Lot-size coding")381 df_m, xcols = unstandardized_matrix(md["df_clean"], md["X_const"])382 ln_lot_imp = lot_recoding(df_m, xcols, raw_m)383 print("\n4. State-by-state heterogeneity")384 state_heterogeneity(df_m, xcols, ln_lot_imp)385 print("\n4b. XGBoost on the trimmed geographic test set")386 xgb_geo_trimmed(df)387 print("\n5. Bootstrap CIs for test-set R2")388 r2_bootstrap(ml, df)389 print("\n6. Moran's I on test residuals, OLS vs tree ensembles")390 moran_xgb_vs_ols(ml, df)391 print(f"\nDone in {time.time() - t0:.0f}s")392393394if __name__ == "__main__":395 main()396