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"""Generate the LaTeX table fragments (paper/tables/*.tex) and the three new5figures (figures/fig12-fig14) of the October 2026 revision from the CSV outputs6of scripts 05-07 in results/revision/. Nothing is typed by hand: every number7in these tables comes from a stored CSV.89Usage: python scripts/08_revision_tables_figures.py10"""1112import sys13from pathlib import Path1415import matplotlib16matplotlib.use("Agg")17import matplotlib.pyplot as plt # noqa: E40218import numpy as np # noqa: E40219import pandas as pd # noqa: E4022021sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))22from wp3 import config # noqa: E40223from wp3.plotting import ACCENT, NEUTRAL, PRIMARY, PRIMARY_LIGHT, apply_paper_style, label, panel_title # noqa: E4022425REV = config.RESULTS_DIR / "revision"26TAB = config.PROJECT_ROOT / "paper" / "tables"27FIG = config.FIGURES_DIR28TAB.mkdir(exist_ok=True, parents=True)2930HOLDOUT = set(config.HOLDOUT_STATES)313233def f(x, d=3, plus=False):34 if pd.isna(x):35 return "---"36 s = f"{x:+.{d}f}" if plus else f"{x:.{d}f}"37 return s.replace("-", "$-$")383940def n(x):41 return f"{int(x):,}"424344def write(name, body):45 r"""Write a table body. The final row is written without its trailing46 line break: the \\ is supplied by the LaTeX source after \input, because a47 row terminator at the very end of an \input file leaves tokens between48 \cr and the following \bottomrule (``Misplaced \noalign'')."""49 body = body.rstrip()50 if body.endswith("\\\\"):51 body = body[:-2].rstrip()52 (TAB / name).write_text(body + "\n")53 print("wrote", TAB / name)545556def tex_label(v):57 return label(v).replace("Log(", "$\\ln$(")585960# --------------------------------------------------------------------------- tables61def table_missingness():62 m = pd.read_csv(REV / "missingness_full.csv")63 imp = pd.read_csv(REV / "implausible_values.csv", index_col=0)["n"]64 nice = {"avg_school_rating": "Average school rating", "school_count": "School count",65 "nearest_school_distance": "Nearest-school distance", "transit_score": "Transit Score",66 "hoa_fee_value (ln HOA)": "HOA fee ($\\ln$ HOA)", "has_garage": "Garage indicator",67 "roof_type / construction_material / foundation_type": "Roof, construction, foundation type",68 "flooring (central air / forced air / hardwood dummies)": "Flooring (hardwood dummy)",69 "stories": "Stories", "parking_spaces": "Parking spaces",70 "year_built (age, age^2, interactions)": "Year built (age, age$^2$, interactions)",71 "lot_size_sqft (ln lot)": "Lot size ($\\ln$ lot)",72 "cooling_type (central air / forced air / hardwood dummies)": "Cooling type (central-air dummy)",73 "heating_type (central air / forced air / hardwood dummies)": "Heating type (forced-air dummy)",74 "bike_score": "Bike Score", "walk_score": "Walk Score"}75 rows = []76 for _, r in m.iterrows():77 if r["missing_share"] < 0.0005:78 continue79 rows.append(f"{nice.get(r['input'], r['input']).replace('_', '\\_')} & "80 f"{r['missing_share'] * 100:.1f} & {r['imputation'].replace('_', '\\_')} \\\\")81 body = "\n".join(rows)82 write("missingness.tex", body)83 rows = [f"{k.replace('_', '\\_').replace('>', '$>$')} & {n(v)} & {v / 788842 * 100:.2f} \\\\"84 for k, v in imp.items()]85 write("implausible.tex", "\n".join(rows))868788def table_cluster_se():89 c = pd.read_csv(REV / "cluster_se.csv", index_col=0)90 key = ["ln_sqft", "bedrooms", "bathrooms", "age", "age_sq", "ln_lot", "has_pool",91 "has_garage", "on_waterfront", "has_central_air", "has_hardwood", "luxury_score",92 "walk_score", "bike_score", "transit_score", "avg_school_rating",93 "property_tax_rate", "is_condo", "has_hoa", "tag_new_construction",94 "tag_foreclosure", "region_Northeast", "region_South", "region_West"]95 rows = []96 for v in key:97 r = c.loc[v]98 rows.append(f"{tex_label(v)} & {f(r.coef, 4)} & {f(r.se_hc3, 4)} & {f(r.se_cluster_zip3, 4)} & "99 f"{f(r.se_cluster_state, 4)} & {f(r.ratio_zip3_hc3, 1)} & {f(r.ratio_state_hc3, 1)} & "100 f"{f(r.t_cluster_state, 2)} \\\\")101 sig = {k: int((c[f"t_{k}"].abs() > 1.96).sum()) for k in ["hc3", "cluster_zip3", "cluster_state"]}102 rows.append("\\midrule")103 rows.append(f"Coefficients significant at 5\\% (of {len(c)}) & & {sig['hc3']} & {sig['cluster_zip3']} & "104 f"{sig['cluster_state']} & & & \\\\")105 rows.append(f"Median SE ratio across all coefficients & & & & & {c.ratio_zip3_hc3.median():.1f} & "106 f"{c.ratio_state_hc3.median():.1f} & \\\\")107 write("cluster_se.tex", "\n".join(rows))108 return c, sig109110111def table_lot_recoding():112 d = pd.read_csv(REV / "lot_recoding.csv")113 labels = ["Published: $\\ln(\\text{lot}) = 0$ if missing or zero",114 "Recoded: log of state-median lot if missing or zero",115 "Recoded + missing-lot indicator",116 "Complete cases (lot reported, positive)",117 "Recoded, lot winsorised at p99.5"]118 rows = []119 for lab, (_, r) in zip(labels, d.iterrows()):120 ind = f(r.beta_lot_missing_indicator, 3) if not pd.isna(r.beta_lot_missing_indicator) else "---"121 rows.append(f"{lab} & {n(r.N)} & {f(r.R2, 3)} & {f(r.beta_ln_lot, 3)} & ({f(r.se_ln_lot, 4)}) & "122 f"{f(r.beta_ln_sqft, 3)} & {ind} \\\\")123 write("lot_recoding.tex", "\n".join(rows))124 return d125126127def table_state_heterogeneity():128 s = pd.read_csv(REV / "state_heterogeneity_summary.csv")129 p = pd.read_csv(REV / "pooled_vs_interacted.csv")130 rows = []131 for _, r in s.iterrows():132 rows.append(f"{tex_label(r.variable)} & {f(r.pooled_state_fe, 3)} & {f(r.p10, 3)} & {f(r["median"], 3)} & "133 f"{f(r.p90, 3)} & {f(r['min'], 3)} & {f(r['max'], 3)} & {int(r.n_states_sign_differs_from_pooled)} \\\\")134 write("state_heterogeneity.tex", "\n".join(rows))135 rows = [f"{r.model} & {f(r.r2_test, 3)} \\\\" for _, r in p.iterrows()]136 write("pooled_vs_interacted.tex", "\n".join(rows))137 return s, p138139140def table_r2_bootstrap():141 b = pd.read_csv(REV / "r2_bootstrap.csv")142 for scheme, fname in [("observation bootstrap", "r2_bootstrap_obs.tex"),143 ("ZIP3-cluster bootstrap", "r2_bootstrap_zip3.tex")]:144 sub = b[b.scheme == scheme]145 rows = []146 for _, r in sub.iterrows():147 q = r.quantity.replace("R2", "$R^2$")148 rows.append(f"\\quad {q} & {f(r.point, 4)} & [{f(r.ci_low, 4)}, {f(r.ci_high, 4)}] & {f(r.boot_se, 4)} \\\\")149 write(fname, "\n".join(rows))150 return b151152153def table_moran_models():154 m = pd.read_csv(REV / "moran_xgb_vs_ols.csv")155 rows = []156 for model in ["OLS", "LightGBM", "XGBoost"]:157 sub = m[m.model == model].sort_values("subsample")158 vals = " & ".join(f(v, 4) for v in sub["I"])159 rows.append(f"{model} & {vals} & {f(sub['I'].mean(), 4)} & {f(sub['I'].std(ddof=0), 4)} & "160 f"{f(sub['z'].min(), 1)}--{f(sub['z'].max(), 1)} \\\\")161 write("moran_models.tex", "\n".join(rows))162 return m163164165def table_designs():166 sm = pd.read_csv(REV / "size_matched_split.csv")167 rows = []168 short = {"random 80/20 (631,073 train)": "Random 80/20 split",169 "random, geographic training size": "Random split, geographic training size",170 "10-state geographic holdout": "Ten-state geographic holdout"}171 for _, r in sm.iterrows():172 r.design = short.get(r.design, r.design)173 within = f(r.r2_within, 3) if not pd.isna(r.r2_within) else "---"174 rows.append(f"{r.design} & {n(r.n_train)} & {n(r.n_test)} & {f(r.r2, 3)} & {within} & "175 f"{f(r.rmse, 3)} & {r.mdape:.1f} \\\\")176 write("designs.tex", "\n".join(rows))177 return sm178179180def table_variants():181 v = pd.read_csv(REV / "variants_two_designs.csv")182 short = {"Full model (62 features)": "Full model",183 "Without school variables (59)": "No school variables",184 "Without region dummies (59)": "No region dummies",185 "Without school variables and region dummies (56)": "No school variables, no region dummies",186 "Full model + 6 missingness indicators (68)": "Full model + 6 missingness indicators"}187 rows = []188 for _, r in v.iterrows():189 rows.append(f"{short.get(r.variant, r.variant)} & {int(r.n_features)} & {f(r.r2_random, 3)} & {f(r.r2_geo, 3)} & "190 f"{f(r.r2_geo_100_trees, 3)} & {f(r.r2_geo_within_state, 3)} & {r.mdape_geo:.1f} \\\\")191 write("variants.tex", "\n".join(rows))192 return v193194195def table_ols_geo():196 o = pd.read_csv(REV / "ols_geo_holdout.csv")197 x = pd.read_csv(REV / "xgb_geo_trimmed.csv")198 rows = []199 short = {"OLS, 62 regressors, all observations": "OLS, all observations",200 "OLS, 62 regressors, implausible values removed": "OLS, implausible records removed",201 "OLS, no region dummies, implausible values removed": "OLS, no region dummies, implausible removed",202 "XGBoost, 100 trees, trimmed test set": "XGBoost, 100 rounds, same trimmed test set",203 "XGBoost, 1000 trees, trimmed test set": "XGBoost, 1,000 rounds, same trimmed test set"}204 for _, r in o.iterrows():205 r.specification = short.get(r.specification, r.specification)206 r2 = f(r.r2, 3) if r.r2 > -10 else f"$-{abs(r.r2) / 1e6:.1f}\\times 10^{{6}}$"207 rw = f(r.r2_within, 3) if r.r2_within > -10 else "---"208 rows.append(f"{r.specification} & {n(r.n_test)} & {r2} & {rw} & {f(r.rmse, 3) if r.rmse < 10 else f(r.rmse, 0)} & {r.mdape:.1f} \\\\")209 for _, r in x.iterrows():210 rows.append(f"{short.get(r.model, r.model)} & {n(r.n_test)} & {f(r.r2, 3)} & --- & {f(r.rmse, 3)} & {r.mdape:.1f} \\\\")211 write("ols_geo.tex", "\n".join(rows))212 return o, x213214215def table_grouped_cv():216 g = pd.read_csv(REV / "grouped_cv_by_state.csv")217 for variant, fname in [("Full model (62 features)", "grouped_cv_full.tex"),218 ("Without school variables (59)", "grouped_cv_noschool.tex")]:219 sub = g[g.variant == variant]220 rows = []221 for _, r in sub.iterrows():222 lab = f"Fold {r.fold}" if r.fold != "pooled" else "Pooled over the five folds"223 rows.append(f"\\quad {lab} & {int(r.n_states)} & {n(r.n_test)} & {f(r.r2, 3)} & "224 f"{f(r.r2_within, 3)} & {f(r.rmse, 3)} & {r.mdape:.1f} \\\\")225 write(fname, "\n".join(rows))226 return g227228229def table_loso():230 p = REV / "loso_by_state.csv"231 if not p.exists():232 print("loso_by_state.csv not available yet"); return None233 l = pd.read_csv(p)234 rows = []235 summ = {}236 for variant in l.variant.unique():237 sub = l[l.variant == variant]238 w = sub.n_test.values239 summ[variant] = {240 "median": sub.r2.median(), "p10": sub.r2.quantile(0.1), "p90": sub.r2.quantile(0.9),241 "min": sub.r2.min(), "max": sub.r2.max(), "wmean": np.average(sub.r2, weights=w),242 "neg": int((sub.r2 < 0).sum()), "bias_med": sub.mean_error.abs().median(),243 "mdape_med": sub.mdape.median(), "within_pooled": np.nan}244 s = summ[variant]245 short = {"Full model (62 features)": "Full model", "Without school variables (59)": "No school variables"}246 rows.append(f"{short.get(variant, variant)} & {f(s['median'], 3)} & {f(s['p10'], 3)} & {f(s['p90'], 3)} & "247 f"{f(s['min'], 3)} & {f(s['max'], 3)} & {f(s['wmean'], 3)} & {s['neg']} & "248 f"{f(s['bias_med'], 3)} & {s['mdape_med']:.1f} \\\\")249 write("loso_summary.tex", "\n".join(rows))250 # per-state table (full model) for the appendix251 full = l[l.variant.str.startswith("Full")].sort_values("r2", ascending=False)252 rows = []253 for _, r in full.iterrows():254 rows.append(f"{r.state} & {n(r.n_test)} & {f(r.r2, 3)} & {f(r.r2_100_trees, 3)} & "255 f"{f(r.mean_error, 3, plus=True)} & {r.mdape:.1f} \\\\")256 write("loso_states.tex", "\n".join(rows))257 return l, summ258259260def table_iqr_bootstrap():261 p = REV / "iqr_wald_bootstrap.csv"262 if not p.exists():263 print("iqr_wald_bootstrap.csv not available yet"); return None264 d = pd.read_csv(p)265 d["absz"] = d.z_boot_zip3.abs()266 d = d.sort_values("absz", ascending=False)267 rows = []268 for _, r in d.iterrows():269 star = "***" if abs(r.z_boot_zip3) > 3.29 else "**" if abs(r.z_boot_zip3) > 2.576 else "*" if abs(r.z_boot_zip3) > 1.96 else ""270 rows.append(f"{tex_label(r.variable)} & {f(r.difference, 3, plus=True)} & {f(r.z_independent, 2)} & "271 f"{f(r.z_boot_obs, 2)} & {f(r.z_boot_zip3, 2)} & [{f(r.ci_low_zip3, 3)}, {f(r.ci_high_zip3, 3)}] & {star} \\\\")272 write("iqr_bootstrap.tex", "\n".join(rows))273 return d274275276def table_zestimate():277 z = pd.read_csv(REV / "list_to_zestimate.csv", index_col=0)278 rows = []279 for q in ["0.1", "0.25", "0.5", "0.75", "0.9"]:280 r = z.loc[float(q)] if float(q) in z.index else z.loc[q]281 rows.append(f"$p_{{{int(float(q) * 100)}}}$ & " + " & ".join(f(r[c], 3) for c in z.columns) + " \\\\")282 rows.append("Share above 1 & " + " & ".join(f"{z.loc['share_above_1', c] * 100:.1f}\\%" for c in z.columns) + " \\\\")283 rows.append("$N$ & " + " & ".join(n(z.loc['n', c]) for c in z.columns) + " \\\\")284 write("zestimate.tex", "\n".join(rows))285 return z286287288def table_trees_curve():289 c = pd.read_csv(REV / "geo_trees_curve.csv")290 rows = [f"{int(r.n_trees)} & {f(r.r2_random_test, 3)} & {f(r.r2_eval_within_training_states, 3)} & {f(r.r2_geo, 3)} \\\\"291 for _, r in c.iterrows() if int(r.n_trees) in (25, 50, 75, 100, 200, 300, 500, 700, 1000)]292 write("trees_curve.tex", "\n".join(rows))293 return c294295296# --------------------------------------------------------------------------- figures297def fig12_trees_curve(c):298 apply_paper_style()299 fig, ax = plt.subplots(figsize=(8.5, 5))300 ax.plot(c.n_trees, c.r2_random_test, color=PRIMARY, lw=2, marker="o", ms=5,301 label="Random 20% test set")302 ax.plot(c.n_trees, c.r2_eval_within_training_states, color=NEUTRAL, lw=2, ls="--", marker="s", ms=4,303 label="Early-stopping set (within training states)")304 ax.plot(c.n_trees, c.r2_geo, color=ACCENT, lw=2, marker="D", ms=5,305 label="Ten held-out states")306 k = c.loc[c.r2_geo.idxmax()]307 ax.annotate(f"peak {k.r2_geo:.3f} at {int(k.n_trees)} rounds", xy=(k.n_trees, k.r2_geo),308 xytext=(k.n_trees + 120, k.r2_geo + 0.06), fontsize=11, color=ACCENT,309 arrowprops=dict(arrowstyle="-", color=ACCENT, lw=0.8))310 last = c.iloc[-1]311 ax.annotate(f"{last.r2_geo:.3f}", xy=(last.n_trees, last.r2_geo), xytext=(8, -4),312 textcoords="offset points", fontsize=11, color=ACCENT)313 ax.set_xlabel("Boosting rounds")314 ax.set_ylabel("Out-of-sample $R^2$")315 ax.set_ylim(0.3, 0.9)316 ax.legend(loc="lower right")317 fig.savefig(FIG / "fig12_geo_trees_curve.png")318 plt.close(fig)319320321def fig13_loso(l):322 apply_paper_style()323 full = l[l.variant.str.startswith("Full")].sort_values("r2")324 fig, ax = plt.subplots(figsize=(8.5, 10))325 colors = [ACCENT if s in HOLDOUT else PRIMARY for s in full.state]326 ax.barh(full.state, full.r2, color=colors, height=0.72)327 ax.axvline(0, color=NEUTRAL, lw=0.9)328 ax.axvline(full.r2.median(), color=NEUTRAL, lw=1, ls=":")329 ax.text(full.r2.median(), len(full) - 0.2, f"median {full.r2.median():.2f}", fontsize=11,330 color=NEUTRAL, ha="center", va="bottom")331 ax.set_xlabel("Leave-one-state-out $R^2$ (XGBoost, 62 features, 1,000 rounds)")332 ax.grid(axis="x", linestyle=":", alpha=0.45); ax.grid(axis="y", visible=False)333 ax.tick_params(axis="y", labelsize=9)334 from matplotlib.patches import Patch335 ax.legend(handles=[Patch(color=ACCENT, label="State in the original ten-state holdout"),336 Patch(color=PRIMARY, label="Other states and DC")], loc="lower right")337 fig.savefig(FIG / "fig13_loso_by_state.png")338 plt.close(fig)339340341def fig14_state_heterogeneity(summ_row_pooled):342 apply_paper_style()343 h = pd.read_csv(REV / "state_heterogeneity.csv").sort_values("ln_sqft")344 fig, axes = plt.subplots(1, 2, figsize=(12, 9), sharey=False)345 for ax, var, pooled, title in [346 (axes[0], "ln_sqft", summ_row_pooled["ln_sqft"], "(a) Living-area elasticity"),347 (axes[1], "has_garage", summ_row_pooled["has_garage"], "(b) Garage coefficient")]:348 hh = h.sort_values(var)349 sizes = 10 + 60 * np.sqrt(hh.n / hh.n.max())350 ax.scatter(hh[var], hh.state, s=sizes, color=PRIMARY, alpha=0.85, edgecolor="white", lw=0.6)351 ax.axvline(pooled, color=ACCENT, lw=1.5, ls="--")352 ax.text(pooled, -1.2, f"pooled, state FE: {pooled:.3f}", color=ACCENT, fontsize=10.5,353 ha="center", va="top")354 ax.axvline(0, color=NEUTRAL, lw=0.8)355 ax.tick_params(axis="y", labelsize=8.5)356 ax.grid(axis="x", linestyle=":", alpha=0.45); ax.grid(axis="y", visible=False)357 panel_title(ax, title)358 ax.set_xlabel("State-by-state OLS coefficient (natural units)")359 fig.tight_layout()360 fig.savefig(FIG / "fig14_state_heterogeneity.png")361 plt.close(fig)362363364def main():365 table_missingness()366 table_cluster_se()367 table_lot_recoding()368 s, p = table_state_heterogeneity()369 table_r2_bootstrap()370 table_moran_models()371 table_designs()372 table_variants()373 table_ols_geo()374 table_grouped_cv()375 table_zestimate()376 c = table_trees_curve()377 fig12_trees_curve(c)378 loso = table_loso()379 if loso is not None:380 fig13_loso(loso[0])381 pooled = dict(zip(s.variable, s.pooled_state_fe))382 fig14_state_heterogeneity(pooled)383 table_iqr_bootstrap()384385386if __name__ == "__main__":387 main()388