spb/wp11_uqo

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 06 — Inference and robustness (added in the 2026-10-17 revision).45Reads the test-set predictions stored by step 02 and runs the additional6analyses requested by the editorial review:781. cluster-bootstrap confidence intervals (municipality clusters) for the9   MdAPE of every model, and Diebold–Mariano / Giacomini–White-type tests10   of equal accuracy against the spline hedonic (A6) and gradient11   boosting (D4), under both holdouts            -> inference_<scheme>.csv122. same-period comparison: random-split models evaluated on the 2025–2613   sales of the random test set versus forward-in-time models on the14   full 2025–26 test set                            -> same_period.csv153. forward-in-time accuracy by test quarter        -> by_quarter.csv164. leakage decomposition of the random holdout: accuracy on test sales17   whose parcel / coordinate pair also appears in training versus the18   others                                           -> leakage.csv195. leakage-free holdouts: the full race under parcel-grouped and20   coordinate-grouped 80/20 splits   -> horserace_group_parcel.csv,21                                        horserace_group_coord.csv226. nested spatial fixed effects (municipality + grid)  -> nested_fe.csv237. Box–Cox on a finer lambda grid                       -> boxcox_fine.csv248. gradient-boosting hyper-parameter sensitivity       -> hgb_sensitivity.csv2510. targeted pairwise comparisons along each axis (bootstrap CI of the26    MdAPE difference and clustered DM tests)              -> pairs.csv279. level-bias decomposition: median log residual (aggregate level error)28   and MdAPE after re-levelling predictions by the test-set median ratio29   (relative pricing skill net of the level error)       -> relevel.csv3031Usage:  python scripts/06_robustness.py          (about 15 min)32"""33import sys34import time35from pathlib import Path3637import numpy as np38import pandas as pd3940sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))4142from wp11 import config, inference, models, sample  # noqa: E4024344OUT = config.REPRODUCED45MODEL_NAMES = [m["name"] for m in models.model_registry()]46A6, D4, A5 = ("A6 Semi-log + splines", "D4 Gradient boosting",47              "A5 Semi-log + quadratics")484950def _mdape(y, p):51    return float(np.median(np.abs(p - y) / y) * 100)525354def step_inference(P: dict):55    for scheme, Ps in P.items():56        t0 = time.time()57        boot = inference.cluster_bootstrap_mdape(Ps, MODEL_NAMES, ref=A6)58        boot_d4 = inference.cluster_bootstrap_mdape(Ps, MODEL_NAMES, ref=D4)59        boot = boot.merge(boot_d4[["model", "diff_vs_ref", "diff_lo",60                                   "diff_hi"]]61                          .rename(columns={"diff_vs_ref": "diff_vs_d4",62                                           "diff_lo": "diff_d4_lo",63                                           "diff_hi": "diff_d4_hi"}),64                          on="model", how="left")65        rows = []66        for m in MODEL_NAMES:67            r = {"model": m}68            for ref, tag in ((A6, "a6"), (D4, "d4")):69                if m == ref:70                    continue71                for loss, ltag in ((inference.ape, "ape"),72                                   (inference.sq_log, "sqlog")):73                    t = inference.dm_test(Ps, m, ref, loss=loss)74                    r[f"dm_{ltag}_vs_{tag}_t"] = t["t"]75                    r[f"dm_{ltag}_vs_{tag}_p"] = t["p"]76                    r[f"dm_{ltag}_vs_{tag}_diff"] = t["mean_diff"]77            rows.append(r)78        res = boot.merge(pd.DataFrame(rows), on="model")79        res["n_clusters"] = Ps["muni"].nunique()80        res.to_csv(OUT / f"inference_{scheme}.csv", index=False)81        print(f"  inference [{scheme}] done in {time.time() - t0:.0f}s")828384def step_same_period(P: dict):85    cut = pd.Timestamp(config.TEMPORAL_CUTOFF)86    Pr = P["random"][P["random"]["sale_date"] >= cut]87    Pt = P["temporal"]88    rows = []89    for m in MODEL_NAMES:90        rows.append({"model": m,91                     "mdape_random_all": _mdape(P["random"]["price"].to_numpy(),92                                                P["random"][m].to_numpy()),93                     "mdape_random_2025_26": _mdape(Pr["price"].to_numpy(),94                                                    Pr[m].to_numpy()),95                     "mdape_forward_2025_26": _mdape(Pt["price"].to_numpy(),96                                                     Pt[m].to_numpy()),97                     "n_random_2025_26": len(Pr), "n_forward": len(Pt)})98    pd.DataFrame(rows).to_csv(OUT / "same_period.csv", index=False)99    print("  same-period comparison done")100101102def step_by_quarter(P: dict):103    Pt = P["temporal"].copy()104    Pt["q"] = Pt["sale_date"].dt.to_period("Q").astype(str)105    rows = []106    for q, g in Pt.groupby("q"):107        for m in (A5, A6, "D2 Ridge (cross-validated)", "D3 Random forest",108                  D4):109            rows.append({"quarter": q, "model": m, "n": len(g),110                         "mdape": _mdape(g["price"].to_numpy(),111                                         g[m].to_numpy())})112    pd.DataFrame(rows).to_csv(OUT / "by_quarter.csv", index=False)113    print("  by-quarter accuracy done")114115116def step_leakage(P: dict):117    Pr = P["random"]118    rows = []119    for key in ("in_train_parcel", "in_train_coord"):120        for flag in (True, False):121            g = Pr[Pr[key] == flag]122            for m in (A5, A6, "D2 Ridge (cross-validated)",123                      "D3 Random forest", D4,124                      "D5 Gradient boosting, no coordinates",125                      "D6 Spatial k-NN comparables"):126                rows.append({"key": key, "seen_in_training": flag,127                             "segment": "all", "n": len(g), "model": m,128                             "mdape": _mdape(g["price"].to_numpy(),129                                             g[m].to_numpy())})130            gc = g[g["prop_class"] == "condo"]131            for m in (A5, D4):132                rows.append({"key": key, "seen_in_training": flag,133                             "segment": "condo", "n": len(gc), "model": m,134                             "mdape": _mdape(gc["price"].to_numpy(),135                                             gc[m].to_numpy())})136    pd.DataFrame(rows).to_csv(OUT / "leakage.csv", index=False)137    print("  leakage decomposition done")138139140def step_relevel(P: dict):141    """Decompose forward-in-time error into an aggregate level error and142    relative pricing skill. For each model the median log residual143    b = median(ln y - ln p) measures the common level error (staleness);144    the re-levelled MdAPE uses p * exp(b), i.e. an oracle time adjustment145    of the kind practitioners obtain from a published index."""146    rows = []147    for scheme, Ps in P.items():148        y = Ps["price"].to_numpy(float)149        for m in MODEL_NAMES:150            p = Ps[m].to_numpy(float)151            b = float(np.median(np.log(y) - np.log(p)))152            rows.append({"scheme": scheme, "model": m,153                         "median_log_resid": b,154                         "mdape": _mdape(y, p),155                         "mdape_relevelled": _mdape(y, p * np.exp(b)),156                         "mape_relevelled": float(np.mean(157                             np.abs(p * np.exp(b) - y) / y) * 100),158                         "rmse_ln_centered": float(np.sqrt(np.mean(159                             (np.log(y) - np.log(p) - b) ** 2)))})160    pd.DataFrame(rows).to_csv(OUT / "relevel.csv", index=False)161    print("  level-bias decomposition done")162163164PAIRS = [("A2 Semi-log", "A1 Linear (levels)"),165         ("A6 Semi-log + splines", "A2 Semi-log"),166         ("A4 Box-Cox", "A2 Semi-log"),167         ("A3 Log-log", "A2 Semi-log"),168         ("B2 Year FE", "B1 No time effects"),169         ("B4 Month FE", "B3 Quarter FE"),170         ("C3 Grid-cell FE (~5.5 km)", "C2 Municipality FE"),171         ("C4 Grid-cell FE (~1.1 km)", "C2 Municipality FE"),172         ("D2 Ridge (cross-validated)", "D1 OLS (muni + month FE)"),173         ("D4 Gradient boosting", "D3 Random forest"),174         ("D4 Gradient boosting", "D5 Gradient boosting, no coordinates"),175         ("D4 Gradient boosting", "A1 Linear (levels)"),176         ("A6 Semi-log + splines", "A1 Linear (levels)")]177178179def step_pairs(P: dict):180    """Targeted comparisons: for each (A, B) the MdAPE difference A - B181    with its cluster-bootstrap interval and the clustered DM statistics182    on the APE and squared-log losses (positive t = B more accurate)."""183    rows = []184    for scheme, Ps in P.items():185        for a, b in PAIRS:186            boot = inference.cluster_bootstrap_mdape(Ps, [a, b], ref=b)187            x = boot.set_index("model").loc[a]188            t1 = inference.dm_test(Ps, a, b, loss=inference.ape)189            t2 = inference.dm_test(Ps, a, b, loss=inference.sq_log)190            rows.append({"scheme": scheme, "model_a": a, "model_b": b,191                         "mdape_a": x["mdape"],192                         "mdape_b": boot.set_index("model").loc[b, "mdape"],193                         "diff": x["diff_vs_ref"], "diff_lo": x["diff_lo"],194                         "diff_hi": x["diff_hi"],195                         "dm_ape_t": t1["t"], "dm_ape_p": t1["p"],196                         "dm_sqlog_t": t2["t"], "dm_sqlog_p": t2["p"]})197    pd.DataFrame(rows).to_csv(OUT / "pairs.csv", index=False)198    print("  targeted pairwise comparisons done")199200201def step_group_splits(df):202    for scheme in ("group_parcel", "group_coord"):203        res = models.run_horserace(df, scheme)204        res.to_csv(OUT / f"horserace_{scheme}.csv", index=False)205206207def step_nested_fe(df):208    rows = []209    for scheme in ("random", "temporal"):210        res = models.run_horserace(df, scheme,211                                   registry=models.robustness_registry())212        res["scheme"] = scheme213        rows.append(res)214    pd.concat(rows).to_csv(OUT / "nested_fe.csv", index=False)215216217def step_boxcox_fine(df):218    d = models.add_derived(df)219    d["sale_year_c"] = d["sale_year"].astype(str)220    train, test = models.split(d, "random")221    grid = [round(x, 2) for x in np.arange(-0.25, 0.751, 0.05)]222    lam, prof = models.fit_boxcox(train, models.NUM_LEVELS,223                                  models.CATS + ["muni", "quarter"],224                                  grid=grid)225    prof["chosen"] = prof["lambda"] == lam226    # out-of-sample accuracy of the Box–Cox model at the refined lambda227    spec = dict(name="A4 fine", group="x", kind="boxcox",228                numeric=models.NUM_LEVELS, fe=["muni", "quarter"],229                boxcox_grid=[lam])230    pred, _ = models.fit_predict(spec, train, test)231    prof.attrs = {}232    prof["mdape_at_chosen"] = np.where(prof["chosen"],233                                       _mdape(test["price"].to_numpy(float),234                                              np.clip(pred, config.PRED_FLOOR,235                                                      None)), np.nan)236    prof.to_csv(OUT / "boxcox_fine.csv", index=False)237    print(f"  Box–Cox fine grid: lambda* = {lam}")238239240def step_hgb_sensitivity(df):241    rows = []242    for scheme in ("random", "temporal"):243        reg = [dict(name="D4 Gradient boosting (baseline)", group="S",244                    kind="hgb")]245        reg += [dict(name=f"D4 {label}", group="S", kind="hgb",246                     hgb_params=params)247                for label, params in models.HGB_ALTERNATIVES.items()]248        res = models.run_horserace(df, scheme, registry=reg)249        res["scheme"] = scheme250        rows.append(res)251    pd.concat(rows).to_csv(OUT / "hgb_sensitivity.csv", index=False)252253254def main() -> None:255    config.ensure_dirs()256    df = sample.load()257    P = {s: pd.read_parquet(config.PROCESSED / f"predictions_{s}.parquet")258         for s in ("random", "temporal")}259    t0 = time.time()260    step_inference(P)261    step_same_period(P)262    step_by_quarter(P)263    step_leakage(P)264    step_relevel(P)265    step_pairs(P)266    step_group_splits(df)267    step_nested_fe(df)268    step_boxcox_fine(df)269    step_hgb_sensitivity(df)270    print(f"Robustness step finished in {(time.time() - t0) / 60:.1f} min")271272273if __name__ == "__main__":274    main()275