spb/wp9_uqo

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        if name in ("M4", "M5"):36            fit[name]["n_fe"] = int(data["fsa_c"].nunique())37        print(f"{name}: N={int(res.nobs):,}  R2={r2:.4f}")38        if name != "M4":39            models.coef_table(res, cols).to_csv(REPRODUCED / f"coef_{name}.csv")40    json.dump(fit, open(REPRODUCED / "fit.json", "w"), indent=2)4142    # Grand model (M5): per-listing predictions, residuals and FSA premia43    res5, cols5, _ = ladder["M5"]44    fe = models.fsa_fixed_effects(s, res5, cols5)45    pred = models.predict_with_fe(s, res5, cols5, fe)46    frame = s[["fsa_c", "prov", "lat", "lon", "ln_price"]].copy()47    frame["pred_grand"] = pred48    frame["resid_grand"] = frame["ln_price"] - frame["pred_grand"]49    frame.to_parquet(REPRODUCED / "grand_model.parquet", index=False)5051    premia = (frame.assign(fe=s["fsa_c"].map(fe))52              .groupby("fsa_c")53              .agg(fe=("fe", "first"), n=("fe", "size"), prov=("prov", "first"),54                   lat=("lat", "mean"), lon=("lon", "mean")))55    premia = premia[premia["n"] >= 50].copy()56    reference_level = premia["fe"].median()57    premia["premium_pct"] = (np.exp(premia["fe"] - reference_level) - 1) * 10058    premia.index.name = "fsa_c"59    premia.sort_values("premium_pct", ascending=False).to_csv(REPRODUCED / "fsa_premia.csv")60    print(f"FSA premia written for {len(premia):,} neighbourhoods (>=50 listings)")616263if __name__ == "__main__":64    main()65