spb/wp9_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 02 — Estimate the M1–M5 hedonic specification ladder.45Writes to ``results/reproduced/``:6  fit.json          R^2 and N for each specification7  coef_M{1,2,3,5}.csv  coefficient tables (coef, se, p)8  grand_model.parquet  per-listing grand-model (M5) fitted values/residuals9  fsa_premia.csv    FSA location premia (grand-model fixed effects, >=50 listings)1011Usage:  python scripts/02_estimate_core.py12"""13import json14import sys15from pathlib import Path1617import numpy as np18import pandas as pd1920sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2122from wp9 import models, sample  # noqa: E40223from wp9.config import REPRODUCED, ensure_dirs  # noqa: E402242526def main() -> None:27    ensure_dirs()28    s = sample.load_sample()29    ladder = models.specification_ladder(s)3031    fit = {}32    for name, (res, cols, data) in ladder.items():33        r2 = float(res.rsquared)34        fit[name] = {"r2": r2, "n": int(res.nobs)}35        print(f"{name}: N={int(res.nobs):,}  R2={r2:.4f}")36        if name != "M4":37            models.coef_table(res, cols).to_csv(REPRODUCED / f"coef_{name}.csv")38    json.dump(fit, open(REPRODUCED / "fit.json", "w"), indent=2)3940    # Grand model (M5): per-listing predictions, residuals and FSA premia41    res5, cols5, _ = ladder["M5"]42    fe = models.fsa_fixed_effects(s, res5, cols5)43    pred = models.predict_with_fe(s, res5, cols5, fe)44    frame = s[["fsa_c", "prov", "lat", "lon", "ln_price"]].copy()45    frame["pred_grand"] = pred46    frame["resid_grand"] = frame["ln_price"] - frame["pred_grand"]47    frame.to_parquet(REPRODUCED / "grand_model.parquet", index=False)4849    premia = (frame.assign(fe=s["fsa_c"].map(fe))50              .groupby("fsa_c")51              .agg(fe=("fe", "first"), n=("fe", "size"), prov=("prov", "first"),52                   lat=("lat", "mean"), lon=("lon", "mean")))53    premia = premia[premia["n"] >= 50].copy()54    reference_level = premia["fe"].median()55    premia["premium_pct"] = (np.exp(premia["fe"] - reference_level) - 1) * 10056    premia.index.name = "fsa_c"57    premia.sort_values("premium_pct", ascending=False).to_csv(REPRODUCED / "fsa_premia.csv")58    print(f"FSA premia written for {len(premia):,} neighbourhoods (>=50 listings)")596061if __name__ == "__main__":62    main()63