spb/wp11_uqo Public
UQO Working Paper No. 11 — Half a million prices, twenty models: a systematic assessment of hedonic specifications.
TeX 54.7%
Python 45.2%
1# Author: Simon-Pierre Boucher — contact@spboucher.ai2"""The WP11 horse-race harness.34Twenty hedonic models organised along four design axes:56 A. functional form — linear, semi-log, log-log, Box–Cox, quadratic,7 splines (linear estimator, muni + quarter FE)8 B. time effects — none / year / quarter / month FE (quadratic backbone)9 C. spatial controls — none / municipality / ~1.1 km grid / neighbourhood10 unit FE (quadratic backbone, quarter FE)11 D. estimation method — OLS, ridge, random forest, gradient boosting12 (with and without coordinates), spatial k-NN13 comparables1415Every linear model is estimated as a sparse ridge with a vanishing penalty16(numerically OLS, but it (i) survives collinear dummies and (ii) predicts17unseen categories at the reference level via ``handle_unknown='ignore'``).18Every model is evaluated on identical train/test splits, and all metrics19are computed on *price levels* — log models are retransformed with Duan's20smearing factor, Box–Cox by analytic inversion — so that functional forms21compete on the same scoreboard.22"""23from __future__ import annotations2425import time2627import numpy as np28import pandas as pd29from sklearn.compose import ColumnTransformer30from sklearn.ensemble import (HistGradientBoostingRegressor,31 RandomForestRegressor)32from sklearn.linear_model import Ridge, RidgeCV33from sklearn.neighbors import KNeighborsRegressor34from sklearn.pipeline import Pipeline35from sklearn.preprocessing import OneHotEncoder, SplineTransformer3637from . import config3839NUM_LEVELS = ["area", "lot", "has_lot", "age", "floors", "units"]40NUM_LOGS = ["ln_area", "ln_lot", "has_lot", "age", "floors", "units"]41NUM_QUAD = NUM_LEVELS + ["area2", "age2"]42CATS = ["prop_class", "link"]43ML_NUM = ["area", "lot", "has_lot", "age", "floors", "units",44 "lat", "lng", "t"]454647def add_derived(df: pd.DataFrame) -> pd.DataFrame:48 """Quadratic terms, scaled to keep the design well conditioned."""49 out = df.copy()50 out["area2"] = (out["area"] / 100.0) ** 251 out["age2"] = (out["age"] / 10.0) ** 252 return out535455# ---------------------------------------------------------------- splits56def split(df: pd.DataFrame, scheme: str):57 """Return (train, test) under the random or forward-in-time scheme."""58 if scheme == "random":59 rng = np.random.default_rng(config.SEED)60 mask = rng.random(len(df)) < config.TEST_SHARE61 return df[~mask], df[mask]62 if scheme == "temporal":63 cut = pd.Timestamp(config.TEMPORAL_CUTOFF)64 return df[df["sale_date"] < cut], df[df["sale_date"] >= cut]65 raise ValueError(scheme)666768# ---------------------------------------------------------------- pipelines69def _linear_pipeline(numeric: list[str], cats: list[str],70 spline_cols: list[str] | None = None,71 alpha: float = 1e-6) -> Pipeline:72 transformers = [("cat",73 OneHotEncoder(handle_unknown="ignore", drop="first"),74 cats)]75 if spline_cols:76 transformers.append(77 ("spl", SplineTransformer(degree=3, n_knots=6), spline_cols))78 passthrough = [c for c in numeric if c not in spline_cols]79 else:80 passthrough = numeric81 if passthrough:82 transformers.append(("num", "passthrough", passthrough))83 ct = ColumnTransformer(transformers, sparse_threshold=1.0)84 return Pipeline([("ct", ct),85 ("reg", Ridge(alpha=alpha, solver="sparse_cg"))])868788def _ml_frame(df: pd.DataFrame, geo: bool = True) -> pd.DataFrame:89 cols = ML_NUM if geo else [c for c in ML_NUM if c not in ("lat", "lng")]90 X = df[cols].copy()91 for c in CATS:92 X[c] = df[c].astype("category").cat.codes93 return X949596# ---------------------------------------------------------------- metrics97def _metrics(y_level: np.ndarray, pred_level: np.ndarray) -> dict:98 pred_level = np.clip(pred_level, config.PRED_FLOOR, None)99 ape = np.abs(pred_level - y_level) / y_level100 ln_y, ln_p = np.log(y_level), np.log(pred_level)101 err = ln_y - ln_p102 return {"mdape": float(np.median(ape) * 100),103 "mape": float(np.mean(ape) * 100),104 "rmse_ln": float(np.sqrt(np.mean(err ** 2))),105 "r2_ln": float(1 - np.sum(err ** 2)106 / np.sum((ln_y - ln_y.mean()) ** 2))}107108109# ---------------------------------------------------------------- Box–Cox110def _boxcox(y: np.ndarray, lam: float) -> np.ndarray:111 return np.log(y) if lam == 0 else (y ** lam - 1.0) / lam112113114def _inv_boxcox(z: np.ndarray, lam: float) -> np.ndarray:115 if lam == 0:116 return np.exp(z)117 base = np.clip(lam * z + 1.0, 1e-6, None)118 return base ** (1.0 / lam)119120121def fit_boxcox(train: pd.DataFrame, numeric: list[str], cats: list[str]):122 """Profile-likelihood choice of λ on the training sample."""123 y = train["price"].to_numpy(float)124 n = len(y)125 sum_ln_y = np.log(y).sum()126 rows = []127 for lam in config.BOXCOX_GRID:128 pipe = _linear_pipeline(numeric, cats)129 z = _boxcox(y, lam)130 pipe.fit(train, z)131 sse = float(np.sum((z - pipe.predict(train)) ** 2))132 ll = -0.5 * n * np.log(sse / n) + (lam - 1.0) * sum_ln_y133 rows.append({"lambda": lam, "loglik": ll})134 prof = pd.DataFrame(rows)135 lam_star = float(prof.loc[prof["loglik"].idxmax(), "lambda"])136 return lam_star, prof137138139# ---------------------------------------------------------------- registry140def model_registry() -> list[dict]:141 """The twenty models. ``fe`` lists categorical FE columns appended to142 the base dummies; ``numeric``/``spline`` define the design."""143 A = [144 dict(name="A1 Linear (levels)", group="A. Functional form",145 kind="linear_level", numeric=NUM_LEVELS, fe=["muni", "quarter"]),146 dict(name="A2 Semi-log", group="A. Functional form",147 kind="linear_log", numeric=NUM_LEVELS, fe=["muni", "quarter"]),148 dict(name="A3 Log-log", group="A. Functional form",149 kind="linear_log", numeric=NUM_LOGS, fe=["muni", "quarter"]),150 dict(name="A4 Box-Cox", group="A. Functional form",151 kind="boxcox", numeric=NUM_LEVELS, fe=["muni", "quarter"]),152 dict(name="A5 Semi-log + quadratics", group="A. Functional form",153 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "quarter"]),154 dict(name="A6 Semi-log + splines", group="A. Functional form",155 kind="linear_log", numeric=NUM_LEVELS, fe=["muni", "quarter"],156 spline=["area", "age", "lot"]),157 ]158 B = [159 dict(name="B1 No time effects", group="B. Time effects",160 kind="linear_log", numeric=NUM_QUAD, fe=["muni"]),161 dict(name="B2 Year FE", group="B. Time effects",162 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "sale_year_c"]),163 dict(name="B3 Quarter FE", group="B. Time effects",164 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "quarter"]),165 dict(name="B4 Month FE", group="B. Time effects",166 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "month"]),167 ]168 C = [169 dict(name="C1 No spatial controls", group="C. Spatial controls",170 kind="linear_log", numeric=NUM_QUAD, fe=["quarter"]),171 dict(name="C2 Municipality FE", group="C. Spatial controls",172 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "quarter"]),173 dict(name="C3 Grid-cell FE (~5.5 km)", group="C. Spatial controls",174 kind="linear_log", numeric=NUM_QUAD, fe=["grid5", "quarter"]),175 dict(name="C4 Grid-cell FE (~1.1 km)", group="C. Spatial controls",176 kind="linear_log", numeric=NUM_QUAD, fe=["grid", "quarter"]),177 ]178 D = [179 dict(name="D1 OLS (muni + month FE)", group="D. Estimation method",180 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "month"]),181 dict(name="D2 Ridge (cross-validated)", group="D. Estimation method",182 kind="ridge", numeric=NUM_QUAD, fe=["muni", "month"]),183 dict(name="D3 Random forest", group="D. Estimation method",184 kind="rf"),185 dict(name="D4 Gradient boosting", group="D. Estimation method",186 kind="hgb"),187 dict(name="D5 Gradient boosting, no coordinates",188 group="D. Estimation method", kind="hgb_nogeo"),189 dict(name="D6 Spatial k-NN comparables", group="D. Estimation method",190 kind="knn"),191 ]192 return A + B + C + D193194195# ---------------------------------------------------------------- fit/eval196def fit_predict(spec: dict, train: pd.DataFrame, test: pd.DataFrame):197 """Fit one registry entry, return level predictions on the test set."""198 kind = spec["kind"]199 y_train = train["price"].to_numpy(float)200201 if kind in ("linear_level", "linear_log", "boxcox", "ridge"):202 cats = CATS + spec.get("fe", [])203 pipe = _linear_pipeline(spec["numeric"], cats,204 spline_cols=spec.get("spline"))205 if kind == "linear_level":206 pipe.fit(train, y_train)207 return pipe.predict(test), pipe208 if kind in ("linear_log", "ridge"):209 if kind == "ridge":210 from sklearn.preprocessing import StandardScaler211 pipe.steps[-1:] = [212 ("sc", StandardScaler(with_mean=False)),213 ("reg", RidgeCV(alphas=np.logspace(-6, 2, 17),214 cv=3))]215 z = np.log(y_train)216 pipe.fit(train, z)217 resid = z - pipe.predict(train)218 smear = float(np.mean(np.exp(resid))) # Duan (1983)219 return np.exp(pipe.predict(test)) * smear, pipe220 lam, _ = fit_boxcox(train, spec["numeric"], cats)221 pipe = _linear_pipeline(spec["numeric"], cats)222 pipe.fit(train, _boxcox(y_train, lam))223 spec["lambda"] = lam224 return _inv_boxcox(pipe.predict(test), lam), pipe225226 if kind == "rf":227 m = RandomForestRegressor(n_estimators=120, min_samples_leaf=5,228 max_features=0.5, n_jobs=-1,229 random_state=config.SEED)230 m.fit(_ml_frame(train), np.log(y_train))231 return np.exp(m.predict(_ml_frame(test))), m232233 if kind in ("hgb", "hgb_nogeo"):234 geo = kind == "hgb"235 m = HistGradientBoostingRegressor(236 max_iter=600, learning_rate=0.08, max_leaf_nodes=63,237 min_samples_leaf=20, l2_regularization=1e-2,238 random_state=config.SEED)239 m.fit(_ml_frame(train, geo), np.log(y_train))240 return np.exp(m.predict(_ml_frame(test, geo))), m241242 if kind == "knn":243 # price-per-m2 of the 10 nearest sold neighbours (lat/lng, km-scaled)244 Xtr = train[["lat", "lng"]].to_numpy() * np.array([111.0, 78.0])245 Xte = test[["lat", "lng"]].to_numpy() * np.array([111.0, 78.0])246 m = KNeighborsRegressor(n_neighbors=10, weights="distance", n_jobs=-1)247 m.fit(Xtr, np.log(y_train) - train["ln_area"].to_numpy())248 return np.exp(m.predict(Xte) + test["ln_area"].to_numpy()), m249250 raise ValueError(kind)251252253def run_horserace(df: pd.DataFrame, scheme: str) -> pd.DataFrame:254 """Evaluate the full registry under one split scheme.255256 Under the forward-in-time split, test-period time categories (months,257 quarters, years) were never seen in training. A one-hot encoder would258 silently price them at the *baseline* period; the honest forecasting259 rule for a static model is to freeze the price level at the last260 period observed in training, so unseen time labels are remapped to the261 latest training label before prediction.262 """263 d = add_derived(df)264 d["sale_year_c"] = d["sale_year"].astype(str)265 train, test = split(d, scheme)266 if scheme == "temporal":267 test = test.copy()268 for col in ("quarter", "month", "sale_year_c"):269 last = train[col].max()270 test.loc[~test[col].isin(train[col].unique()), col] = last271 y_test = test["price"].to_numpy(float)272 rows = []273 for spec in model_registry():274 t0 = time.time()275 pred, _ = fit_predict(dict(spec), train, test)276 met = _metrics(y_test, pred)277 rows.append({"name": spec["name"], "group": spec["group"],278 **met, "seconds": round(time.time() - t0, 1),279 "n_train": len(train), "n_test": len(test)})280 print(f" [{scheme}] {spec['name']:<38} "281 f"MdAPE={met['mdape']:5.1f}% R2_ln={met['r2_ln']:.3f} "282 f"({rows[-1]['seconds']}s)")283 return pd.DataFrame(rows)284