spb/wp22_uqo

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Phase 2 — benchmark indices A–K and SHARP at province × class and region × class (retrospective, full sample).45Output: results/benchmarks/series.parquet (long: class, level, geo_id, method, month, log_index, se, n) + timing.csv.6GBM (J) is run at province level and for the three largest regions only (cost), see docs/decision_log.md D2.1.7"""8from __future__ import annotations910import sys11import time12from pathlib import Path1314sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))15import numpy as np16import polars as pl1718from sharp_hpi import config as C19from sharp_hpi import estimators as E20from sharp_hpi.hedonic import assign_loc_fe, build_features21from sharp_hpi.qhpi import qhpi_cell, qhpi_fullpool_residuals, qhpi_province22from sharp_hpi.sharp import cell_log_index, run_sharp2324OUT = C.RESULTS / "benchmarks"25OUT.mkdir(parents=True, exist_ok=True)26GRID = E.month_grid("2021-01", "2026-10")27GBM_REGIONS = {"06", "03", "16"}282930def rec(rows, timing, cls, level, gid, s: E.IndexSeries, t0):31    f = s.to_frame().with_columns(pl.lit(cls).alias("class"), pl.lit(level).alias("level"), pl.lit(gid).alias("geo_id"))32    if "se" not in f.columns:33        f = f.with_columns(pl.lit(None, dtype=pl.Float64).alias("se"))34    if "n" not in f.columns:35        f = f.with_columns(pl.lit(None, dtype=pl.Float64).alias("n"))36    rows.append(f.select(["class", "level", "geo_id", "method", "month", "log_index", "se", "n"]))37    timing.append({"class": cls, "level": level, "geo_id": gid, "method": s.method, "seconds": round(time.time() - t0, 2)})38    print(cls, level, gid, s.method, f"{time.time() - t0:.1f}s", flush=True)394041def main() -> None:42    a = pl.read_parquet(C.ANALYTIC_TX)43    f = build_features(a)44    rows, timing = [], []45    for cls in ("house", "condo", "plex"):46        d = f.filter(pl.col("property_class") == cls)47        prov = assign_loc_fe(d, min_n=15).with_columns(pl.col("mrc_code").alias("loc_coarse"))48        # --- province ---49        for fn, kw in [(E.median_index, {"stat": "median"}), (E.median_index, {"stat": "mean"}), (E.stratified_median, {}),50                       (E.pooled_td, {"robust": True}), (E.pooled_td, {"robust": False}),51                       (E.rolling_td, {"window": 13, "splice": "mean"}), (E.rolling_td, {"window": 13, "splice": "movement"}), (E.rolling_td, {"window": 13, "splice": "window"}),52                       (E.rolling_td, {"window": 25, "splice": "mean"}), (E.rolling_td, {"window": 7, "splice": "mean"}),53                       (E.hedonic_imputation, {"kind": "laspeyres"}), (E.hedonic_imputation, {"kind": "paasche"}), (E.hedonic_imputation, {"kind": "fisher"}),54                       (E.hedonic_imputation, {"kind": "fisher", "chained": False}), (E.geks, {"window": 25}), (E.geks, {"window": 13}),55                       (E.repeat_sales, {"weighted": True}), (E.repeat_sales, {"weighted": False}), (E.gbm_residual_index, {})]:56            t0 = time.time(); rec(rows, timing, cls, "province", "QC", fn(prov, GRID, **kw), t0)57        t0 = time.time(); sK = qhpi_province(prov, GRID); sK.method = "K_qhpi"; rec(rows, timing, cls, "province", "QC", sK, t0)58        t0 = time.time(); r = run_sharp(d, GRID); s1 = r.stage159        rec(rows, timing, cls, "province", "QC", E.IndexSeries(GRID, s1.mu, se=s1.mu_se, n=E._counts(d, GRID), method="S_sharp"), t0)60        t0 = time.time(); r0 = run_sharp(d, GRID, time_varying=False)61        rec(rows, timing, cls, "province", "QC", E.IndexSeries(GRID, r0.stage1.mu, se=r0.stage1.mu_se, n=E._counts(d, GRID), method="S_sharp_notv"), t0)62        # QHPI hierarchy inputs63        qres, qsig = qhpi_fullpool_residuals(prov, GRID)64        # --- regions ---65        for reg in sorted(d["region_code"].unique().to_list()):66            dr = d.filter(pl.col("region_code") == reg)67            if dr.height < 500:68                continue69            drr = assign_loc_fe(dr, min_n=15).with_columns(pl.col("mun_code").alias("loc_coarse"))70            for fn, kw in [(E.median_index, {"stat": "median"}), (E.stratified_median, {}), (E.pooled_td, {"robust": True}),71                           (E.rolling_td, {"window": 13, "splice": "mean"}), (E.hedonic_imputation, {"kind": "fisher"}), (E.geks, {"window": 25}),72                           (E.repeat_sales, {"weighted": True})]:73                t0 = time.time()74                try:75                    rec(rows, timing, cls, "region", reg, fn(drr, GRID, **kw), t0)76                except Exception as ex:77                    print("FAILED", cls, reg, fn.__name__, ex, flush=True)78            if reg in GBM_REGIONS and cls == "house":79                t0 = time.time(); rec(rows, timing, cls, "region", reg, E.gbm_residual_index(drr, GRID), t0)80            t0 = time.time()81            path, var, meth, gain = qhpi_cell(drr, qres.filter(pl.col("region_code") == reg), sK.log_index, GRID, qsig)82            rec(rows, timing, cls, "region", reg, E.IndexSeries(GRID, path, se=np.sqrt(var), n=E._counts(dr, GRID), method=f"K_qhpi_{meth}"), t0)83            if ("region", reg) in r.cells:84                p, v = cell_log_index(r, "region", reg, "smoothed")85                rec(rows, timing, cls, "region", reg, E.IndexSeries(GRID, p, se=np.sqrt(v), n=r.cells[("region", reg)].n, method="S_sharp"), time.time())86                p, v = cell_log_index(r, "region", reg, "filtered")87                rec(rows, timing, cls, "region", reg, E.IndexSeries(GRID, p, se=np.sqrt(v), n=r.cells[("region", reg)].n, method="S_sharp_filtered"), time.time())88        pl.concat(rows).write_parquet(OUT / "series.parquet")89        pl.DataFrame(timing).write_csv(OUT / "timing.csv")90    print("done")919293if __name__ == "__main__":94    main()95