#!/usr/bin/env python3 # Author: Simon-Pierre Boucher — contact@spboucher.ai """Step 03 — Extensions beyond the scoreboard. 1. Constant-quality price indices implied by the month fixed effects of four linear functional forms, plus a gradient-boosting index obtained by repricing a fixed reference sample at each month. 2. Implicit-price comparison: the age and floor-area profiles implied by the quadratic OLS versus the partial-dependence profile of the gradient-boosting model. 3. Learning curves: accuracy versus training-set size for OLS and GB. 4. Segment analysis: MdAPE by property class and municipality size. Writes: index.csv, profiles.csv, learning.csv, segments.csv Usage: python scripts/03_extensions.py """ import sys from pathlib import Path import numpy as np import pandas as pd sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) from wp11 import config, models, sample # noqa: E402 from sklearn.ensemble import HistGradientBoostingRegressor # noqa: E402 def month_index_from_linear(spec, train, months): """Fit a month-FE linear model and read the index off the dummies. Numeric regressors are z-scored first: month-dummy coefficients are invariant to that reparametrization, but the conjugate-gradient solver needs the improved conditioning for the individual coefficients (not just the fit) to be accurate when levels like lot area (up to 10^5 m2) share the design with dummies. """ d2 = train.copy() for c in spec.get("numeric", []): sd = d2[c].std() if sd > 0: d2[c] = (d2[c] - d2[c].mean()) / sd for c in spec.get("spline", []): d2[c] = train[c] # splines keep their raw support train = d2 pred, pipe = models.fit_predict(dict(spec), train, train.head(50)) ct = pipe.named_steps["ct"] names = ct.get_feature_names_out() coefs = pipe.named_steps["reg"].coef_ idx = {} for nm, c in zip(names, coefs): if "month_" in nm: idx[nm.split("month_")[-1]] = float(c) base = sorted(months)[0] ref = idx.get(base, 0.0) return {m: 100 * float(np.exp(idx.get(m, 0.0) - ref)) for m in months} def main() -> None: config.ensure_dirs() df = sample.load() out = config.REPRODUCED d = models.add_derived(df) d["sale_year_c"] = d["sale_year"].astype(str) months = sorted(d["month"].unique()) # ------------------------------------------------------------ 1. indices forms = { "Semi-log": dict(kind="linear_log", numeric=models.NUM_LEVELS, fe=["muni", "month"]), "Log-log": dict(kind="linear_log", numeric=models.NUM_LOGS, fe=["muni", "month"]), "Semi-log + quadratics": dict(kind="linear_log", numeric=models.NUM_QUAD, fe=["muni", "month"]), "Semi-log + splines": dict(kind="linear_log", numeric=models.NUM_LEVELS, fe=["muni", "month"], spline=["area", "age", "lot"]), } rows = [] for label, spec in forms.items(): spec["name"], spec["group"] = label, "index" idx = month_index_from_linear(spec, d, months) rows += [{"month": m, "form": label, "index": v} for m, v in idx.items()] print(f" index from {label}: done") # gradient boosting: reprice a fixed 20k reference sample each month rng = np.random.default_rng(config.SEED) hgb = HistGradientBoostingRegressor(max_iter=600, learning_rate=0.08, max_leaf_nodes=63, min_samples_leaf=20, l2_regularization=1e-2, random_state=config.SEED) hgb.fit(models._ml_frame(d), d["ln_price"].to_numpy()) ref = d.sample(20_000, random_state=config.SEED).copy() t_by_month = d.groupby("month")["t"].mean() base_val = None for m in months: ref["t"] = t_by_month[m] v = float(np.exp(hgb.predict(models._ml_frame(ref))).mean()) base_val = base_val or v rows.append({"month": m, "form": "Gradient boosting", "index": 100 * v / base_val}) pd.DataFrame(rows).to_csv(out / "index.csv", index=False) print(" index from Gradient boosting: done") # ------------------------------------------------------------ 2. profiles train, _ = models.split(d, "random") spec = dict(name="A5", group="x", kind="linear_log", numeric=models.NUM_QUAD, fe=["muni", "quarter"]) _, pipe = models.fit_predict(spec, train, train.head(50)) names = pipe.named_steps["ct"].get_feature_names_out() coefs = dict(zip(names, pipe.named_steps["reg"].coef_)) prof_rows = [] ages = np.arange(0, 121, 5) b_age = coefs.get("num__age", 0.0) b_age2 = coefs.get("num__age2", 0.0) for a in ages: y = b_age * a + b_age2 * (a / 10.0) ** 2 prof_rows.append({"var": "age", "x": a, "model": "OLS quadratic", "y": y - (b_age * ages[0])}) areas = np.arange(50, 401, 25) b_ar = coefs.get("num__area", 0.0) b_ar2 = coefs.get("num__area2", 0.0) y0 = b_ar * areas[0] + b_ar2 * (areas[0] / 100.0) ** 2 for a in areas: y = b_ar * a + b_ar2 * (a / 100.0) ** 2 prof_rows.append({"var": "area", "x": a, "model": "OLS quadratic", "y": y - y0}) # partial dependence of the GB model (average prediction on a grid) sub = models._ml_frame(train.sample(30_000, random_state=config.SEED)) for var, grid in (("age", ages), ("area", areas)): vals = [] for g in grid: s2 = sub.copy() s2[var] = g vals.append(float(hgb.predict(s2).mean())) v0 = vals[0] prof_rows += [{"var": var, "x": g, "model": "Gradient boosting", "y": v - v0} for g, v in zip(grid, vals)] pd.DataFrame(prof_rows).to_csv(out / "profiles.csv", index=False) print(" implicit-price profiles: done") # ------------------------------------------------------------ 3. learning curves train_full, test = models.split(d, "random") y_test = test["price"].to_numpy(float) rows = [] for n in config.LEARNING_SIZES: if n > len(train_full): continue tr = train_full.sample(n, random_state=config.SEED) for label, spec in ( ("OLS quadratic (muni+quarter FE)", dict(name="A5", group="x", kind="linear_log", numeric=models.NUM_QUAD, fe=["muni", "quarter"])), ("Gradient boosting", dict(name="D4", group="x", kind="hgb"))): pred, _ = models.fit_predict(spec, tr, test) met = models._metrics(y_test, pred) rows.append({"n_train": n, "model": label, **met}) print(f" learning n={n:>7,} {label:<34} " f"MdAPE={met['mdape']:.1f}%") pd.DataFrame(rows).to_csv(out / "learning.csv", index=False) # ------------------------------------------------------------ 4. segments rows = [] muni_sales = d.groupby("muni")["muni"].transform("size") segs = {f"Class: {c}": d["prop_class"] == c for c in ("single_family", "condo", "plex", "cottage")} segs.update({"Muni < 1k sales": muni_sales < 1_000, "Muni 1k–10k sales": muni_sales.between(1_000, 10_000), "Muni > 10k sales": muni_sales > 10_000}) train, test = models.split(d, "random") preds = {} for label, spec in ( ("OLS quadratic", dict(name="A5", group="x", kind="linear_log", numeric=models.NUM_QUAD, fe=["muni", "quarter"])), ("Gradient boosting", dict(name="D4", group="x", kind="hgb"))): preds[label], _ = models.fit_predict(spec, train, test) y_test = test["price"].to_numpy(float) for seg, mask in segs.items(): m = mask.loc[test.index].to_numpy() if m.sum() < 500: continue for label, p in preds.items(): met = models._metrics(y_test[m], p[m]) rows.append({"segment": seg, "model": label, "n": int(m.sum()), **met}) pd.DataFrame(rows).to_csv(out / "segments.csv", index=False) print(" segment analysis: done") if __name__ == "__main__": main()