UQO Working Paper No. 11 — Half a million prices, twenty models: a systematic assessment of hedonic specifications.
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 / ~5.5 km grid / ~1.1 km10 grid 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* so that functional forms compete on the20same scoreboard. Two retransformation rules are available for every model21whose dependent variable is a transformation of the price:2223 * ``"median"`` (main scoreboard) — the fitted transformed value is24 inverted as is: exp(.) for log models, the analytic inverse for25 Box–Cox. Under symmetric transformed-scale errors this is the26 conditional *median* price, the target that MdAPE/MAPE reward and the27 quantity the tree ensembles (fit on ln P) deliver natively. Log-scale28 metrics (RMSE_ln, R2_ln) are always computed on this prediction.29 * ``"mean"`` (robustness) — Duan's (1983) smearing estimator of the30 conditional *mean*: exp(z_hat) * mean(exp(e_i)) for log models and31 mean_i g^{-1}(z_hat + e_i) over training residuals e_i for Box–Cox;32 tree ensembles receive the same correction from their (in-sample)33 training residuals.3435Version 1.0 of the paper smeared the linear log models but not the tree36ensembles, which penalised the linear models on every metric; the37``"median"`` rule applied uniformly is the scoreboard of version 1.1.38"""39from __future__ import annotations4041import time4243import numpy as np44import pandas as pd45from sklearn.compose import ColumnTransformer46from sklearn.ensemble import (HistGradientBoostingRegressor,47 RandomForestRegressor)48from sklearn.linear_model import Ridge, RidgeCV49from sklearn.neighbors import KNeighborsRegressor50from sklearn.pipeline import Pipeline51from sklearn.preprocessing import OneHotEncoder, SplineTransformer5253from . import config5455NUM_LEVELS = ["area", "lot", "has_lot", "age", "floors", "units"]56NUM_LOGS = ["ln_area", "ln_lot", "has_lot", "age", "floors", "units"]57NUM_QUAD = NUM_LEVELS + ["area2", "age2"]58CATS = ["prop_class", "link"]59ML_NUM = ["area", "lot", "has_lot", "age", "floors", "units",60 "lat", "lng", "t"]616263def add_derived(df: pd.DataFrame) -> pd.DataFrame:64 """Quadratic terms, scaled to keep the design well conditioned."""65 out = df.copy()66 out["area2"] = (out["area"] / 100.0) ** 267 out["age2"] = (out["age"] / 10.0) ** 268 return out697071# ---------------------------------------------------------------- splits72def split(df: pd.DataFrame, scheme: str):73 """Return (train, test) under the random or forward-in-time scheme."""74 if scheme == "random":75 rng = np.random.default_rng(config.SEED)76 mask = rng.random(len(df)) < config.TEST_SHARE77 return df[~mask], df[mask]78 if scheme == "temporal":79 cut = pd.Timestamp(config.TEMPORAL_CUTOFF)80 return df[df["sale_date"] < cut], df[df["sale_date"] >= cut]81 if scheme in ("group_parcel", "group_coord"):82 # Leakage-free random holdouts: 20% of *groups* (parcels, or exact83 # coordinate pairs = buildings/complexes) go to the test set, so84 # that no property (resp. building) is seen on both sides.85 key = "parcel" if scheme == "group_parcel" else "coord"86 rng = np.random.default_rng(config.SEED)87 groups = df[key].unique()88 test_groups = set(groups[rng.random(len(groups)) < config.TEST_SHARE])89 mask = df[key].isin(test_groups).to_numpy()90 return df[~mask], df[mask]91 raise ValueError(scheme)929394# ---------------------------------------------------------------- pipelines95def _linear_pipeline(numeric: list[str], cats: list[str],96 spline_cols: list[str] | None = None,97 alpha: float = 1e-6) -> Pipeline:98 transformers = [("cat",99 OneHotEncoder(handle_unknown="ignore", drop="first"),100 cats)]101 if spline_cols:102 transformers.append(103 ("spl", SplineTransformer(degree=3, n_knots=6), spline_cols))104 passthrough = [c for c in numeric if c not in spline_cols]105 else:106 passthrough = numeric107 if passthrough:108 transformers.append(("num", "passthrough", passthrough))109 ct = ColumnTransformer(transformers, sparse_threshold=1.0)110 return Pipeline([("ct", ct),111 ("reg", Ridge(alpha=alpha, solver="sparse_cg"))])112113114def _ml_frame(df: pd.DataFrame, geo: bool = True) -> pd.DataFrame:115 cols = ML_NUM if geo else [c for c in ML_NUM if c not in ("lat", "lng")]116 X = df[cols].copy()117 for c in CATS:118 X[c] = df[c].astype("category").cat.codes119 return X120121122# ---------------------------------------------------------------- metrics123def _metrics(y_level: np.ndarray, pred_level: np.ndarray) -> dict:124 pred_level = np.clip(pred_level, config.PRED_FLOOR, None)125 ape = np.abs(pred_level - y_level) / y_level126 ln_y, ln_p = np.log(y_level), np.log(pred_level)127 err = ln_y - ln_p128 return {"mdape": float(np.median(ape) * 100),129 "mape": float(np.mean(ape) * 100),130 "rmse_ln": float(np.sqrt(np.mean(err ** 2))),131 "r2_ln": float(1 - np.sum(err ** 2)132 / np.sum((ln_y - ln_y.mean()) ** 2))}133134135# ---------------------------------------------------------------- Box–Cox136def _boxcox(y: np.ndarray, lam: float) -> np.ndarray:137 return np.log(y) if lam == 0 else (y ** lam - 1.0) / lam138139140def _inv_boxcox(z: np.ndarray, lam: float) -> np.ndarray:141 if lam == 0:142 return np.exp(z)143 base = np.clip(lam * z + 1.0, 1e-6, None)144 return base ** (1.0 / lam)145146147def _smear_boxcox(z_hat: np.ndarray, resid: np.ndarray, lam: float,148 n_resid: int = 1000) -> np.ndarray:149 """Duan's smearing estimator of E[P | X] for a Box–Cox model:150 the average of g^{-1}(z_hat + e_i) over (a subsample of) the training151 residuals e_i. Reduces to exp(z_hat) * mean(exp(e)) when lam = 0."""152 rng = np.random.default_rng(config.SEED)153 e = resid if len(resid) <= n_resid else rng.choice(resid, n_resid,154 replace=False)155 out = np.zeros_like(z_hat)156 for chunk in np.array_split(np.arange(len(z_hat)),157 max(1, len(z_hat) // 5000)):158 out[chunk] = _inv_boxcox(z_hat[chunk, None] + e[None, :],159 lam).mean(axis=1)160 return out161162163def fit_boxcox(train: pd.DataFrame, numeric: list[str], cats: list[str],164 grid: list[float] | None = None):165 """Profile-likelihood choice of λ on the training sample."""166 y = train["price"].to_numpy(float)167 n = len(y)168 sum_ln_y = np.log(y).sum()169 rows = []170 for lam in (grid if grid is not None else config.BOXCOX_GRID):171 pipe = _linear_pipeline(numeric, cats)172 z = _boxcox(y, lam)173 pipe.fit(train, z)174 sse = float(np.sum((z - pipe.predict(train)) ** 2))175 ll = -0.5 * n * np.log(sse / n) + (lam - 1.0) * sum_ln_y176 rows.append({"lambda": lam, "loglik": ll})177 prof = pd.DataFrame(rows)178 lam_star = float(prof.loc[prof["loglik"].idxmax(), "lambda"])179 return lam_star, prof180181182HGB_PARAMS = dict(max_iter=600, learning_rate=0.08, max_leaf_nodes=63,183 min_samples_leaf=20, l2_regularization=1e-2)184185# Alternative gradient-boosting configurations for the hyper-parameter186# sensitivity analysis (Section "Robustness"). ``early_stopping`` holds out187# the last 10% of the training *rows in time order* as a validation set,188# i.e. a time-honest tuning rule.189HGB_ALTERNATIVES = {190 "shallow (15 leaves, lr 0.05, 1500 iter)":191 dict(max_iter=1500, learning_rate=0.05, max_leaf_nodes=15,192 min_samples_leaf=20, l2_regularization=1e-2),193 "deep (255 leaves, lr 0.05, 600 iter)":194 dict(max_iter=600, learning_rate=0.05, max_leaf_nodes=255,195 min_samples_leaf=20, l2_regularization=1e-2),196 "regularised (min leaf 200, L2 = 1)":197 dict(max_iter=600, learning_rate=0.08, max_leaf_nodes=63,198 min_samples_leaf=200, l2_regularization=1.0),199 "early stopping on last 10% of training months":200 dict(max_iter=2000, learning_rate=0.05, max_leaf_nodes=63,201 min_samples_leaf=20, l2_regularization=1e-2,202 early_stopping=True),203}204205206def robustness_registry() -> list[dict]:207 """Additional specifications: nested spatial FE (unseen grid cells fall208 back to their municipality instead of an arbitrary reference cell)."""209 return [210 dict(name="C3' Municipality + grid-cell FE (~5.5 km)",211 group="R. Nested spatial FE", kind="linear_log",212 numeric=NUM_QUAD, fe=["muni", "grid5", "quarter"]),213 dict(name="C4' Municipality + grid-cell FE (~1.1 km)",214 group="R. Nested spatial FE", kind="linear_log",215 numeric=NUM_QUAD, fe=["muni", "grid", "quarter"]),216 ]217218219# ---------------------------------------------------------------- registry220def model_registry() -> list[dict]:221 """The twenty models. ``fe`` lists categorical FE columns appended to222 the base dummies; ``numeric``/``spline`` define the design."""223 A = [224 dict(name="A1 Linear (levels)", group="A. Functional form",225 kind="linear_level", numeric=NUM_LEVELS, fe=["muni", "quarter"]),226 dict(name="A2 Semi-log", group="A. Functional form",227 kind="linear_log", numeric=NUM_LEVELS, fe=["muni", "quarter"]),228 dict(name="A3 Log-log", group="A. Functional form",229 kind="linear_log", numeric=NUM_LOGS, fe=["muni", "quarter"]),230 dict(name="A4 Box-Cox", group="A. Functional form",231 kind="boxcox", numeric=NUM_LEVELS, fe=["muni", "quarter"]),232 dict(name="A5 Semi-log + quadratics", group="A. Functional form",233 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "quarter"]),234 dict(name="A6 Semi-log + splines", group="A. Functional form",235 kind="linear_log", numeric=NUM_LEVELS, fe=["muni", "quarter"],236 spline=["area", "age", "lot"]),237 ]238 B = [239 dict(name="B1 No time effects", group="B. Time effects",240 kind="linear_log", numeric=NUM_QUAD, fe=["muni"]),241 dict(name="B2 Year FE", group="B. Time effects",242 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "sale_year_c"]),243 dict(name="B3 Quarter FE", group="B. Time effects",244 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "quarter"]),245 dict(name="B4 Month FE", group="B. Time effects",246 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "month"]),247 ]248 C = [249 dict(name="C1 No spatial controls", group="C. Spatial controls",250 kind="linear_log", numeric=NUM_QUAD, fe=["quarter"]),251 dict(name="C2 Municipality FE", group="C. Spatial controls",252 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "quarter"]),253 dict(name="C3 Grid-cell FE (~5.5 km)", group="C. Spatial controls",254 kind="linear_log", numeric=NUM_QUAD, fe=["grid5", "quarter"]),255 dict(name="C4 Grid-cell FE (~1.1 km)", group="C. Spatial controls",256 kind="linear_log", numeric=NUM_QUAD, fe=["grid", "quarter"]),257 ]258 D = [259 dict(name="D1 OLS (muni + month FE)", group="D. Estimation method",260 kind="linear_log", numeric=NUM_QUAD, fe=["muni", "month"]),261 dict(name="D2 Ridge (cross-validated)", group="D. Estimation method",262 kind="ridge", numeric=NUM_QUAD, fe=["muni", "month"]),263 dict(name="D3 Random forest", group="D. Estimation method",264 kind="rf"),265 dict(name="D4 Gradient boosting", group="D. Estimation method",266 kind="hgb"),267 dict(name="D5 Gradient boosting, no coordinates",268 group="D. Estimation method", kind="hgb_nogeo"),269 dict(name="D6 Spatial k-NN comparables", group="D. Estimation method",270 kind="knn"),271 ]272 return A + B + C + D273274275# ---------------------------------------------------------------- fit/eval276def fit_predict(spec: dict, train: pd.DataFrame, test: pd.DataFrame,277 retransform: str = "median"):278 """Fit one registry entry, return level predictions on the test set.279280 ``retransform`` is ``"median"`` (plain inversion of the transformed281 prediction) or ``"mean"`` (Duan smearing). ``spec`` is updated in282 place with ``smear`` (log models), ``lambda`` (Box–Cox) and, when283 ``retransform="mean"`` is requested, both predictions are returned284 through ``spec["pred_median"]`` / ``spec["pred_mean"]``.285 """286 kind = spec["kind"]287 y_train = train["price"].to_numpy(float)288 both = retransform == "mean"289290 def _log_out(z_hat_test, z_hat_train, model):291 resid = np.log(y_train) - z_hat_train292 smear = float(np.mean(np.exp(resid))) # Duan (1983)293 spec["smear"] = smear294 p_med = np.exp(z_hat_test)295 spec["pred_median"], spec["pred_mean"] = p_med, p_med * smear296 return (p_med * smear if both else p_med), model297298 if kind in ("linear_level", "linear_log", "boxcox", "ridge"):299 cats = CATS + spec.get("fe", [])300 pipe = _linear_pipeline(spec["numeric"], cats,301 spline_cols=spec.get("spline"))302 if kind == "linear_level":303 pipe.fit(train, y_train)304 p = pipe.predict(test)305 spec["pred_median"] = spec["pred_mean"] = p306 return p, pipe307 if kind in ("linear_log", "ridge"):308 if kind == "ridge":309 from sklearn.preprocessing import StandardScaler310 pipe.steps[-1:] = [311 ("sc", StandardScaler(with_mean=False)),312 ("reg", RidgeCV(alphas=np.logspace(-6, 2, 17),313 cv=3))]314 pipe.fit(train, np.log(y_train))315 return _log_out(pipe.predict(test), pipe.predict(train), pipe)316 lam, _ = fit_boxcox(train, spec["numeric"], cats,317 grid=spec.get("boxcox_grid"))318 pipe = _linear_pipeline(spec["numeric"], cats)319 z = _boxcox(y_train, lam)320 pipe.fit(train, z)321 spec["lambda"] = lam322 z_hat = pipe.predict(test)323 p_med = _inv_boxcox(z_hat, lam)324 spec["pred_median"] = p_med325 if both:326 resid = z - pipe.predict(train)327 spec["pred_mean"] = _smear_boxcox(z_hat, resid, lam)328 return spec["pred_mean"], pipe329 return p_med, pipe330331 if kind == "rf":332 m = RandomForestRegressor(n_estimators=120, min_samples_leaf=5,333 max_features=0.5, n_jobs=-1,334 random_state=config.SEED)335 m.fit(_ml_frame(train), np.log(y_train))336 return _log_out(m.predict(_ml_frame(test)),337 m.predict(_ml_frame(train)) if both338 else np.log(y_train), m)339340 if kind in ("hgb", "hgb_nogeo"):341 geo = kind == "hgb"342 params = dict(spec.get("hgb_params", HGB_PARAMS))343 tr = train344 if params.pop("early_stopping", False):345 # time-honest validation: the last 10% of training rows in346 # sale-date order serve as the early-stopping set347 tr = train.sort_values("sale_date")348 n_val = int(0.1 * len(tr))349 Xv, yv = _ml_frame(tr.tail(n_val), geo), np.log(350 tr.tail(n_val)["price"].to_numpy(float))351 tr = tr.head(len(tr) - n_val)352 m = HistGradientBoostingRegressor(random_state=config.SEED,353 **params)354 # scikit-learn's early stopping draws a random validation355 # split; emulate a time-ordered one by staged fitting356 m.set_params(early_stopping=False)357 best, best_iter = np.inf, params["max_iter"]358 m.set_params(warm_start=True, max_iter=100)359 for it in range(100, params["max_iter"] + 1, 100):360 m.set_params(max_iter=it)361 m.fit(_ml_frame(tr, geo),362 np.log(tr["price"].to_numpy(float)))363 mse = float(np.mean((yv - m.predict(Xv)) ** 2))364 if mse < best - 1e-5:365 best, best_iter = mse, it366 elif it - best_iter >= 300:367 break368 spec["best_iter"] = best_iter369 m = HistGradientBoostingRegressor(random_state=config.SEED,370 **{**params,371 "max_iter": best_iter})372 tr = train373 else:374 m = HistGradientBoostingRegressor(random_state=config.SEED,375 **params)376 m.fit(_ml_frame(tr, geo), np.log(tr["price"].to_numpy(float)))377 return _log_out(m.predict(_ml_frame(test, geo)),378 m.predict(_ml_frame(train, geo)) if both379 else np.log(y_train), m)380381 if kind == "knn":382 # price-per-m2 of the 10 nearest sold neighbours (lat/lng, km-scaled)383 Xtr = train[["lat", "lng"]].to_numpy() * np.array([111.0, 78.0])384 Xte = test[["lat", "lng"]].to_numpy() * np.array([111.0, 78.0])385 m = KNeighborsRegressor(n_neighbors=10, weights="distance", n_jobs=-1)386 m.fit(Xtr, np.log(y_train) - train["ln_area"].to_numpy())387 z_te = m.predict(Xte) + test["ln_area"].to_numpy()388 z_tr = (m.predict(Xtr) + train["ln_area"].to_numpy() if both389 else np.log(y_train))390 return _log_out(z_te, z_tr, m)391392 raise ValueError(kind)393394395def prepare(df: pd.DataFrame, scheme: str):396 """Derived columns, split, and the carry-forward rule for the397 forward-in-time scheme. Returns (train, test)."""398 d = add_derived(df)399 d["sale_year_c"] = d["sale_year"].astype(str)400 train, test = split(d, scheme)401 if scheme == "temporal":402 test = test.copy()403 for col in ("quarter", "month", "sale_year_c"):404 last = train[col].max()405 test.loc[~test[col].isin(train[col].unique()), col] = last406 return train, test407408409def run_horserace(df: pd.DataFrame, scheme: str,410 registry: list[dict] | None = None,411 return_predictions: bool = False):412 """Evaluate a registry (default: the twenty models) under one split.413414 Under the forward-in-time split, test-period time categories (months,415 quarters, years) were never seen in training. A one-hot encoder would416 silently price them at the *baseline* period; the honest forecasting417 rule for a static model is to freeze the price level at the last418 period observed in training, so unseen time labels are remapped to the419 latest training label before prediction. Tree ensembles are piecewise420 constant in the month counter, so a test month beyond the training421 range is priced by the terminal leaf — the tree analogue of the same422 carry-forward rule.423424 Metrics are reported for the median-consistent prediction (main425 scoreboard, suffix none) and for the Duan-smeared mean prediction426 (suffix ``_mean``). With ``return_predictions`` the test-set427 predictions of every model are returned as a second DataFrame.428 """429 train, test = prepare(df, scheme)430 y_test = test["price"].to_numpy(float)431 rows, preds = [], {}432 for spec in (registry if registry is not None else model_registry()):433 spec = dict(spec)434 t0 = time.time()435 fit_predict(spec, train, test, retransform="mean")436 p_med, p_mean = spec["pred_median"], spec["pred_mean"]437 met = _metrics(y_test, p_med)438 met_mean = {f"{k}_mean": v for k, v in439 _metrics(y_test, p_mean).items()}440 rows.append({"name": spec["name"], "group": spec["group"],441 **met, **met_mean,442 "smear": spec.get("smear", np.nan),443 "lambda": spec.get("lambda", np.nan),444 "best_iter": spec.get("best_iter", np.nan),445 "seconds": round(time.time() - t0, 1),446 "n_train": len(train), "n_test": len(test)})447 preds[spec["name"]] = np.clip(p_med, config.PRED_FLOOR, None)448 print(f" [{scheme}] {spec['name']:<44} "449 f"MdAPE={met['mdape']:5.1f}% R2_ln={met['r2_ln']:.3f} "450 f"({rows[-1]['seconds']}s)")451 res = pd.DataFrame(rows)452 if not return_predictions:453 return res454 P = test[["id", "sale_date", "price", "muni", "prop_class", "grid",455 "parcel", "coord"]].reset_index(drop=True)456 P["in_train_parcel"] = P["parcel"].isin(set(train["parcel"]))457 P["in_train_coord"] = P["coord"].isin(set(train["coord"]))458 for k, v in preds.items():459 P[k] = v460 return res, P461