spb/wp11_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 03 — Extensions beyond the scoreboard.451. Constant-quality price indices implied by the month fixed effects of6   four linear functional forms, plus a gradient-boosting index obtained7   by repricing a fixed reference sample at each month.82. Implicit-price comparison: the age and floor-area profiles implied by9   the quadratic OLS versus the partial-dependence profile of the10   gradient-boosting model.113. Learning curves: accuracy versus training-set size for OLS and GB.124. Segment analysis: MdAPE by property class and municipality size.1314Writes: index.csv, profiles.csv, learning.csv, segments.csv1516Usage:  python scripts/03_extensions.py17"""18import sys19from pathlib import Path2021import numpy as np22import pandas as pd2324sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2526from wp11 import config, models, sample  # noqa: E40227from sklearn.ensemble import HistGradientBoostingRegressor  # noqa: E402282930def month_index_from_linear(spec, train, months):31    """Fit a month-FE linear model and read the index off the dummies.3233    Numeric regressors are z-scored first: month-dummy coefficients are34    invariant to that reparametrization, but the conjugate-gradient35    solver needs the improved conditioning for the individual36    coefficients (not just the fit) to be accurate when levels like lot37    area (up to 10^5 m2) share the design with dummies.38    """39    d2 = train.copy()40    for c in spec.get("numeric", []):41        sd = d2[c].std()42        if sd > 0:43            d2[c] = (d2[c] - d2[c].mean()) / sd44    for c in spec.get("spline", []):45        d2[c] = train[c]  # splines keep their raw support46    train = d247    pred, pipe = models.fit_predict(dict(spec), train, train.head(50))48    ct = pipe.named_steps["ct"]49    names = ct.get_feature_names_out()50    coefs = pipe.named_steps["reg"].coef_51    idx = {}52    for nm, c in zip(names, coefs):53        if "month_" in nm:54            idx[nm.split("month_")[-1]] = float(c)55    base = sorted(months)[0]56    ref = idx.get(base, 0.0)57    return {m: 100 * float(np.exp(idx.get(m, 0.0) - ref)) for m in months}585960def main() -> None:61    config.ensure_dirs()62    df = sample.load()63    out = config.REPRODUCED64    d = models.add_derived(df)65    d["sale_year_c"] = d["sale_year"].astype(str)66    months = sorted(d["month"].unique())6768    # ------------------------------------------------------------ 1. indices69    forms = {70        "Semi-log": dict(kind="linear_log", numeric=models.NUM_LEVELS,71                         fe=["muni", "month"]),72        "Log-log": dict(kind="linear_log", numeric=models.NUM_LOGS,73                        fe=["muni", "month"]),74        "Semi-log + quadratics": dict(kind="linear_log",75                                      numeric=models.NUM_QUAD,76                                      fe=["muni", "month"]),77        "Semi-log + splines": dict(kind="linear_log",78                                   numeric=models.NUM_LEVELS,79                                   fe=["muni", "month"],80                                   spline=["area", "age", "lot"]),81    }82    rows = []83    for label, spec in forms.items():84        spec["name"], spec["group"] = label, "index"85        idx = month_index_from_linear(spec, d, months)86        rows += [{"month": m, "form": label, "index": v}87                 for m, v in idx.items()]88        print(f"  index from {label}: done")8990    # gradient boosting: reprice a fixed 20k reference sample each month91    rng = np.random.default_rng(config.SEED)92    hgb = HistGradientBoostingRegressor(max_iter=600, learning_rate=0.08,93                                        max_leaf_nodes=63, min_samples_leaf=20,94                                        l2_regularization=1e-2,95                                        random_state=config.SEED)96    hgb.fit(models._ml_frame(d), d["ln_price"].to_numpy())97    ref = d.sample(20_000, random_state=config.SEED).copy()98    t_by_month = d.groupby("month")["t"].mean()99    base_val = None100    for m in months:101        ref["t"] = t_by_month[m]102        v = float(np.exp(hgb.predict(models._ml_frame(ref))).mean())103        base_val = base_val or v104        rows.append({"month": m, "form": "Gradient boosting",105                     "index": 100 * v / base_val})106    pd.DataFrame(rows).to_csv(out / "index.csv", index=False)107    print("  index from Gradient boosting: done")108109    # ------------------------------------------------------------ 2. profiles110    train, _ = models.split(d, "random")111    spec = dict(name="A5", group="x", kind="linear_log",112                numeric=models.NUM_QUAD, fe=["muni", "quarter"])113    _, pipe = models.fit_predict(spec, train, train.head(50))114    names = pipe.named_steps["ct"].get_feature_names_out()115    coefs = dict(zip(names, pipe.named_steps["reg"].coef_))116    prof_rows = []117    ages = np.arange(0, 121, 5)118    b_age = coefs.get("num__age", 0.0)119    b_age2 = coefs.get("num__age2", 0.0)120    for a in ages:121        y = b_age * a + b_age2 * (a / 10.0) ** 2122        prof_rows.append({"var": "age", "x": a, "model": "OLS quadratic",123                          "y": y - (b_age * ages[0])})124    areas = np.arange(50, 401, 25)125    b_ar = coefs.get("num__area", 0.0)126    b_ar2 = coefs.get("num__area2", 0.0)127    y0 = b_ar * areas[0] + b_ar2 * (areas[0] / 100.0) ** 2128    for a in areas:129        y = b_ar * a + b_ar2 * (a / 100.0) ** 2130        prof_rows.append({"var": "area", "x": a, "model": "OLS quadratic",131                          "y": y - y0})132133    # partial dependence of the GB model (average prediction on a grid)134    sub = models._ml_frame(train.sample(30_000, random_state=config.SEED))135    for var, grid in (("age", ages), ("area", areas)):136        vals = []137        for g in grid:138            s2 = sub.copy()139            s2[var] = g140            vals.append(float(hgb.predict(s2).mean()))141        v0 = vals[0]142        prof_rows += [{"var": var, "x": g, "model": "Gradient boosting",143                       "y": v - v0} for g, v in zip(grid, vals)]144    pd.DataFrame(prof_rows).to_csv(out / "profiles.csv", index=False)145    print("  implicit-price profiles: done")146147    # ------------------------------------------------------------ 3. learning curves148    train_full, test = models.split(d, "random")149    y_test = test["price"].to_numpy(float)150    rows = []151    for n in config.LEARNING_SIZES:152        if n > len(train_full):153            continue154        tr = train_full.sample(n, random_state=config.SEED)155        for label, spec in (156                ("OLS quadratic (muni+quarter FE)",157                 dict(name="A5", group="x", kind="linear_log",158                      numeric=models.NUM_QUAD, fe=["muni", "quarter"])),159                ("Gradient boosting",160                 dict(name="D4", group="x", kind="hgb"))):161            pred, _ = models.fit_predict(spec, tr, test)162            met = models._metrics(y_test, pred)163            rows.append({"n_train": n, "model": label, **met})164            print(f"  learning n={n:>7,} {label:<34} "165                  f"MdAPE={met['mdape']:.1f}%")166    pd.DataFrame(rows).to_csv(out / "learning.csv", index=False)167168    # ------------------------------------------------------------ 4. segments169    rows = []170    muni_sales = d.groupby("muni")["muni"].transform("size")171    segs = {f"Class: {c}": d["prop_class"] == c172            for c in ("single_family", "condo", "plex", "cottage")}173    segs.update({"Muni < 1k sales": muni_sales < 1_000,174                 "Muni 1k–10k sales": muni_sales.between(1_000, 10_000),175                 "Muni > 10k sales": muni_sales > 10_000})176    train, test = models.split(d, "random")177    preds = {}178    for label, spec in (179            ("OLS quadratic", dict(name="A5", group="x", kind="linear_log",180                                   numeric=models.NUM_QUAD,181                                   fe=["muni", "quarter"])),182            ("Gradient boosting", dict(name="D4", group="x", kind="hgb"))):183        preds[label], _ = models.fit_predict(spec, train, test)184    y_test = test["price"].to_numpy(float)185    for seg, mask in segs.items():186        m = mask.loc[test.index].to_numpy()187        if m.sum() < 500:188            continue189        for label, p in preds.items():190            met = models._metrics(y_test[m], p[m])191            rows.append({"segment": seg, "model": label, "n": int(m.sum()),192                         **met})193    pd.DataFrame(rows).to_csv(out / "segments.csv", index=False)194    print("  segment analysis: done")195196197if __name__ == "__main__":198    main()199