TeX 55.1%
Python 44.2%
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Phase 5 — complete monthly SHARP-HPI panel (retrospective research series, full sample) at province, region, MRC and4municipality level for house / condo / plex, plus coverage, reliability and representative-price tables and the5QHPI-architecture benchmark panel. Writes data/processed/index_panel.parquet|csv and results/panel/*.csv."""6from __future__ import annotations78import json9import sys10import time11from pathlib import Path1213sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))14import numpy as np15import polars as pl1617from sharp_hpi import config as C18from sharp_hpi.estimators import month_grid19from sharp_hpi.hedonic import build_features20from sharp_hpi.panel import build_rows21from sharp_hpi.sharp import LEVELS, cell_log_index, run_sharp2223OUT = C.RESULTS / "panel"24OUT.mkdir(parents=True, exist_ok=True)25GRID = month_grid("2021-01", "2026-10")26VINTAGE = C.SNAPSHOT_TAKEN[:10]272829def names(geo: pl.DataFrame) -> dict:30 d = {}31 for lvl, key, nm in [("region", "region_code", "region_name"), ("mrc", "mrc_code", "mrc_name"), ("municipality", "mun_code", "mun_name")]:32 for k, v in geo.select([key, nm]).unique().iter_rows():33 d[(lvl, k)] = v34 return d353637def main() -> None:38 a = pl.read_parquet(C.ANALYTIC_TX)39 f = build_features(a)40 nm = names(pl.read_parquet(C.GEO_TX, columns=["region_code", "region_name", "mrc_code", "mrc_name", "mun_code", "mun_name"]))41 rows, rep, cov, meta = [], [], [], {}42 for cls in C.PROPERTY_TYPES:43 d = f.filter(pl.col("property_class") == cls)44 t0 = time.time()45 r = run_sharp(d, GRID)46 meta[cls] = {"n": r.stage1.n, "sigma_robust": r.stage1.sigma2 ** 0.5, "r2_within": r.stage1.r2_within, "priors": {k: [float(v[0]), float(v[1])] for k, v in r.priors.items()},47 "seconds": round(time.time() - t0, 1), "coef": {k: float(v) for k, v in r.stage1.coef.items() if "month" not in k}}48 n_prov = np.array([d.filter(pl.col("month") == m).height for m in GRID], float)49 rows.append(build_rows("QC", "Québec (province)", "province", cls, GRID, r.stage1.mu, r.stage1.mu_se ** 2, n_prov, n_prov, np.ones(len(GRID)),50 "SHARP-HPI v%s (stage-1 robust TD)" % C.MODEL_VERSION, C.SAMPLE_END, VINTAGE, C.MODEL_VERSION))51 for (lvl, gid), c in r.cells.items():52 p, v = cell_log_index(r, lvl, gid, "smoothed")53 ess = np.where(c.fit.smoothed_var > 0, c.sigma2_cell / np.maximum(c.fit.smoothed_var, 1e-12), np.inf)54 ess = np.minimum(ess, 1e6)55 rows.append(build_rows(gid, nm.get((lvl, gid), gid), lvl, cls, GRID, p, v, c.n, ess, c.fit.gain, "SHARP-HPI v%s (hierarchical, smoothed)" % C.MODEL_VERSION,56 C.SAMPLE_END, VINTAGE, C.MODEL_VERSION))57 # reliability bundle58 cov.append({"geography_level": lvl, "geography_id": gid, "property_type": cls, "n_total": int(c.n.sum()), "months_with_sales": int((c.n > 0).sum()),59 "median_monthly_n": float(np.median(c.n[c.n > 0])) if (c.n > 0).any() else 0.0, "tau": float(np.sqrt(c.fit.tau2)), "sigma_cell": float(np.sqrt(c.sigma2_cell)),60 "mean_gain": float(np.nanmean(c.fit.gain[c.n > 0])) if (c.n > 0).any() else 0.0, "mean_ess": float(np.mean(ess)), "mean_overlap": float(np.nanmean(c.overlap)),61 "acf1_dev": c.rho1, "oos_z_rmse": float(np.sqrt(np.nanmean(c.oos_z ** 2))) if np.isfinite(c.oos_z).any() else np.nan,62 "innov_weight_min": float(np.min(c.fit.innov_weights))})63 # representative price at 2021 prices: exp(mean(log P - log index)) over the cell's sample64 sub = d.filter(pl.col(dict(LEVELS)[lvl]) == gid)65 li = dict(zip(GRID, p)); adj = np.array([li[m] for m in sub["month"].to_list()])66 base = [i for i, m in enumerate(GRID) if m.startswith("2021")]67 b2021 = float(np.mean(p[base]))68 rep.append({"geography_level": lvl, "geography_id": gid, "property_type": cls, "representative_price_2021": float(np.exp(np.mean(sub["log_price"].to_numpy() - adj + b2021))),69 "n": sub.height})70 panel = pl.concat(rows)71 panel.write_parquet(C.PROCESSED / "index_panel.parquet"); panel.write_csv(C.PROCESSED / "index_panel.csv")72 pl.DataFrame(cov).write_csv(OUT / "reliability_cells.csv")73 pl.DataFrame(rep).write_csv(OUT / "representative_prices.csv")74 (OUT / "stage1_meta.json").write_text(json.dumps(meta, indent=2))75 # coverage table: every geography × type × month76 panel.select(["geography_level", "geography_id", "geography_name", "property_type", "month", "actual_transaction_count", "effective_sample_size",77 "reliability_grade", "direct_or_model_assisted", "suppression_flag"]).write_csv(OUT / "coverage_cells_months.csv")78 panel.group_by(["geography_level", "property_type", "reliability_grade"]).agg(pl.len().alias("cell_months")).sort(["geography_level", "property_type", "reliability_grade"]).write_csv(OUT / "grade_distribution.csv")79 print(panel.shape, meta)808182if __name__ == "__main__":83 main()84