TeX 55.1%
Python 44.2%
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Phase 6 — validation experiments on the real data (protocol §13–15). Sub-commands (each writes results/validation/*):4 repeat municipal/regional repeat-sales agreement (SHARP vs Case–Shiller repeat sales where ≥ 300 pairs)5 composition synthetic quality-composition shock (drop the smallest 40 % of houses from 2024-01): median vs TD vs SHARP6 placebo months permuted within municipality × class: estimated annual changes and CI false-positive rate7 moran residual spatial autocorrelation (Moran's I on k-nearest neighbours) of stage-1 residuals by month8 bootstrap cluster bootstrap (by fine-location cell) of SHARP condo indices vs model-based intervals9 aggregation stock- / transaction- / equal-weighted bottom-up provincial aggregates vs the direct provincial path10 ablation real-data ablations: variants vs full model (path differences, out-of-time innovation error), border sensitivity11 trim sensitivity of the provincial path to the assessment-ratio screen, 1/99 % price trims, type conflicts, new construction12"""13from __future__ import annotations1415import sys16import time17from pathlib import Path1819sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))20import numpy as np21import polars as pl2223from sharp_hpi import config as C24from sharp_hpi import estimators as E25from sharp_hpi.hedonic import assign_loc_fe, build_features26from sharp_hpi.sharp import LEVELS, cell_log_index, run_sharp, stage12728OUT = C.RESULTS / "validation"29OUT.mkdir(parents=True, exist_ok=True)30GRID = E.month_grid("2021-01", "2026-10")313233def load():34 return build_features(pl.read_parquet(C.ANALYTIC_TX))353637def panel_path(level, gid, cls):38 p = pl.read_parquet(C.PROCESSED / "index_panel.parquet").filter((pl.col("geography_level") == level) & (pl.col("geography_id") == gid) & (pl.col("property_type") == cls)).sort("month")39 return np.log(p["index_value"].to_numpy() / 100.0)404142def cmd_repeat(f):43 rows = []44 for cls in C.PROPERTY_TYPES:45 d = f.filter(pl.col("property_class") == cls)46 for level, key in [("province", None), ("region", "region_code"), ("municipality", "mun_code")]:47 ids = ["QC"] if key is None else d.group_by(key).agg(pl.len().alias("n")).filter(pl.col("n") >= 2000)[key].to_list()48 for gid in ids:49 dc = d if key is None else d.filter(pl.col(key) == gid)50 rs = E.repeat_sales(dc, GRID, weighted=True)51 if rs.meta.get("pairs", 0) < 300:52 continue53 s = panel_path(level, gid, cls)54 ok = ~np.isnan(rs.log_index) & ~np.isnan(s)55 a = rs.log_index[ok] - np.mean(rs.log_index[ok]); b = s[ok] - np.mean(s[ok])56 da, db = np.diff(rs.log_index), np.diff(s); okd = ~np.isnan(da) & ~np.isnan(db)57 rows.append({"class": cls, "level": level, "geo_id": gid, "pairs": rs.meta["pairs"], "months": int(ok.sum()), "rmse_level": float(np.sqrt(np.mean((a - b) ** 2))),58 "corr_level": float(np.corrcoef(a, b)[0, 1]), "total_change_rs": float(rs.log_index[ok][-1] - rs.log_index[ok][0]), "total_change_sharp": float(s[ok][-1] - s[ok][0]),59 "corr_12m_changes": float(np.corrcoef((rs.log_index[12:] - rs.log_index[:-12])[~np.isnan(rs.log_index[12:] - rs.log_index[:-12]) & ~np.isnan(s[12:] - s[:-12])], (s[12:] - s[:-12])[~np.isnan(rs.log_index[12:] - rs.log_index[:-12]) & ~np.isnan(s[12:] - s[:-12])])[0, 1])})60 pl.DataFrame(rows).write_csv(OUT / "repeat_sales_agreement.csv"); print(pl.DataFrame(rows))616263def cmd_composition(f):64 rows = []65 for reg in ["06", "03", "07", "16"]:66 d = assign_loc_fe(f.filter((pl.col("property_class") == "house") & (pl.col("region_code") == reg)), min_n=15)67 q40 = d["fa_f"].quantile(0.4)68 shocked = d.filter(~((pl.col("month") >= "2024-01") & (pl.col("fa_f") < q40)))69 for nm, dd in [("original", d), ("shock_drop_small_from_2024", shocked)]:70 med = E.median_index(dd, GRID).log_index; td = E.pooled_td(dd, GRID).log_index71 r = run_sharp(dd, GRID); sh = r.stage1.mu72 for meth, li in [("A_median", med), ("C_pooled_td", td), ("S_sharp", sh)]:73 pre = np.nanmean(li[[GRID.index(m) for m in GRID if "2023-01" <= m <= "2023-12"]]); post = np.nanmean(li[[GRID.index(m) for m in GRID if "2024-01" <= m <= "2024-12"]])74 rows.append({"region": reg, "sample": nm, "method": meth, "n": dd.height, "change_2023_to_2024": float(post - pre)})75 out = pl.DataFrame(rows)76 out = out.pivot(on="sample", index=["region", "method"], values="change_2023_to_2024").with_columns((pl.col("shock_drop_small_from_2024") - pl.col("original")).alias("shock_effect"))77 out.write_csv(OUT / "composition_shock.csv"); print(out)787980def cmd_placebo(f):81 rng = np.random.default_rng(C.SEED)82 rows = []83 for cls in C.PROPERTY_TYPES:84 d = f.filter(pl.col("property_class") == cls)85 for rep in range(3):86 # permute months within municipality (keeps composition, location and overall monthly volumes within cell)87 perm = d.with_columns(pl.Series("_r", rng.random(d.height))).sort(["mun_code", "_r"])88 months_sorted = d.sort(["mun_code", "date"])["month"].to_numpy()89 perm = perm.with_columns(pl.Series("month", months_sorted)).drop("_r")90 # but keep the province-wide path: re-attach the true month effect? No — the placebo tests the *deviation* machinery:91 r = run_sharp(perm, GRID)92 for (lvl, gid), c in r.cells.items():93 dev = c.fit.smoothed; v = c.fit.smoothed_var94 d12 = dev[12:] - dev[:-12]; s12 = np.sqrt(v[12:] + v[:-12])95 rows.append({"class": cls, "rep": rep, "level": lvl, "geo_id": gid, "n_total": int(c.n.sum()), "sd_dev": float(np.std(dev)), "max_abs_12m_change": float(np.max(np.abs(d12))),96 "share_12m_outside_95": float(np.mean(np.abs(d12) > 1.96 * s12)), "tau": float(np.sqrt(c.fit.tau2))})97 out = pl.DataFrame(rows); out.write_csv(OUT / "placebo_cells.csv")98 summ = out.group_by(["class", "level"]).agg(pl.len().alias("cells"), pl.col("sd_dev").median(), pl.col("max_abs_12m_change").median(), pl.col("share_12m_outside_95").mean(), pl.col("tau").median()).sort(["class", "level"])99 summ.write_csv(OUT / "placebo_summary.csv"); print(summ)100101102def cmd_moran(f):103 from scipy.spatial import cKDTree104 rows = []105 for cls in C.PROPERTY_TYPES:106 d = f.filter(pl.col("property_class") == cls)107 s1 = stage1(d, GRID)108 e = s1.resid.join(d.select(["tx_id", "x_m", "y_m"]), on="tx_id").filter(pl.col("e").is_not_null() & pl.col("e").is_not_nan())109 for scope, sub in [("province", e), ("Montréal", e.filter(pl.col("mun_code") == "66023")), ("Québec", e.filter(pl.col("mun_code") == "23027")), ("Gatineau", e.filter(pl.col("mun_code") == "81017"))]:110 for mth in ["2021-06", "2022-06", "2023-06", "2024-06", "2025-06", "2026-06"]:111 m = sub.filter(pl.col("month") == mth)112 if m.height < 200:113 continue114 xy = m.select(["x_m", "y_m"]).to_numpy(); z = m["e"].to_numpy(); z = z - z.mean()115 k = 10; tree = cKDTree(xy); dist, idx = tree.query(xy, k=k + 1); idx = idx[:, 1:]116 lag = z[idx].mean(axis=1); I = float(np.sum(z * lag) / np.sum(z * z))117 # permutation reference118 perm = [float(np.sum(z * rng_z[idx].mean(axis=1)) / np.sum(z * z)) for rng_z in [np.random.default_rng(s).permutation(z) for s in range(99)]]119 rows.append({"class": cls, "scope": scope, "month": mth, "n": m.height, "moran_I_knn10": I, "perm_mean": float(np.mean(perm)), "perm_sd": float(np.std(perm)), "z": float((I - np.mean(perm)) / np.std(perm))})120 out = pl.DataFrame(rows); out.write_csv(OUT / "moran_residuals.csv"); print(out)121122123def cmd_bootstrap(f, B: int = 40):124 rng = np.random.default_rng(C.SEED)125 d = assign_loc_fe(f.filter(pl.col("property_class") == "condo"), min_n=15)126 r0 = run_sharp(d, GRID)127 targets = [("province", "QC")] + [k for k in r0.cells if k[0] == "municipality"][:12] + [k for k in r0.cells if k[0] == "region"][:5]128 base = {k: (r0.stage1.mu if k[0] == "province" else cell_log_index(r0, *k)[0]) for k in targets}129 sev = {k: (r0.stage1.mu_se if k[0] == "province" else np.sqrt(cell_log_index(r0, *k)[1])) for k in targets}130 clusters = d["loc_fe"].unique().to_list(); boots = {k: [] for k in targets}131 for b in range(B):132 pick = rng.choice(clusters, size=len(clusters), replace=True)133 cnt = pl.DataFrame({"loc_fe": pick}).group_by("loc_fe").agg(pl.len().alias("_w"))134 db = d.join(cnt, on="loc_fe", how="inner")135 db = db.select(pl.exclude("_w").repeat_by(pl.col("_w")).explode()) if False else db.with_columns(pl.col("_w")).select(pl.all().repeat_by("_w")).explode(pl.all()).drop("_w")136 rb = run_sharp(db, GRID)137 for k in targets:138 if k[0] == "province":139 boots[k].append(rb.stage1.mu)140 elif k in rb.cells:141 boots[k].append(cell_log_index(rb, *k)[0])142 print("boot", b, flush=True)143 rows = []144 for k in targets:145 arr = np.array(boots[k]); arr = arr - arr[:, :12].mean(axis=1, keepdims=True); bse = arr.std(axis=0)146 ref = base[k] - base[k][:12].mean(); se = sev[k]147 d12b = (arr[:, 12:] - arr[:, :-12]).std(axis=0); d12m = np.sqrt(sev[k][12:] ** 2 + sev[k][:-12] ** 2)148 rows.append({"level": k[0], "geo_id": k[1], "n": int(d.filter(pl.col(dict(LEVELS).get(k[0], "mun_code")) == k[1]).height) if k[0] != "province" else d.height,149 "model_se_level_mean": float(np.mean(se[12:])), "boot_se_level_mean": float(np.mean(bse[12:])), "ratio_level": float(np.mean(bse[12:]) / np.mean(se[12:])),150 "model_se_12m_mean": float(np.mean(d12m)), "boot_se_12m_mean": float(np.mean(d12b)), "ratio_12m": float(np.mean(d12b) / np.mean(d12m))})151 out = pl.DataFrame(rows); out.write_csv(OUT / "bootstrap_vs_model_se.csv"); print(out)152153154def cmd_aggregation(f):155 panel = pl.read_parquet(C.PROCESSED / "index_panel.parquet")156 stock = pl.read_parquet(C.RAW_UNITS_STOCK).filter(pl.col("cubf") == 1000).with_columns(157 pl.when(pl.col("nb_logements") >= 6).then(pl.lit("multi")).when(pl.col("nb_logements").is_between(2, 5)).then(pl.lit("plex"))158 .when((pl.col("nb_logements") <= 1) & (pl.col("lien_physique") == "Intégré")).then(pl.lit("condo")).when(pl.col("nb_logements") <= 1).then(pl.lit("house")).otherwise(None).alias("cls"))159 sw = stock.group_by(["code_mun", "cls"]).agg(pl.len().alias("stock_units"), pl.col("valeur_role").sum().alias("stock_value")).rename({"code_mun": "geography_id"})160 sw.write_csv(OUT / "stock_weights_by_municipality_class.csv")161 rows = []162 for cls in C.PROPERTY_TYPES:163 mun = panel.filter((pl.col("geography_level") == "municipality") & (pl.col("property_type") == cls)).with_columns(np.log(pl.col("index_value") / 100).alias("li"))164 prov = panel.filter((pl.col("geography_level") == "province") & (pl.col("property_type") == cls)).sort("month").with_columns(np.log(pl.col("index_value") / 100).alias("li"))165 tw = mun.group_by("geography_id").agg(pl.col("actual_transaction_count").sum().alias("tx"))166 mun = mun.join(sw.filter(pl.col("cls") == cls).select(["geography_id", "stock_units", "stock_value"]), on="geography_id", how="left").join(tw, on="geography_id").with_columns(pl.lit(1.0).alias("eq"))167 for wname in ["stock_units", "stock_value", "tx", "eq"]:168 agg = mun.filter(pl.col(wname).is_not_null()).group_by("month").agg(((pl.col("li") * pl.col(wname)).sum() / pl.col(wname).sum()).alias("li"), pl.col(wname).sum().alias("w")).sort("month")169 a = agg.join(prov.select(["month", pl.col("li").alias("li_prov")]), on="month")170 li = a["li"].to_numpy().copy(); lp = a["li_prov"].to_numpy().copy(); li -= li[:12].mean(); lp -= lp[:12].mean()171 rows.append({"class": cls, "weights": wname, "cells": mun.filter(pl.col(wname).is_not_null())["geography_id"].n_unique(), "rmse_vs_direct": float(np.sqrt(np.mean((li - lp) ** 2))),172 "total_change_bottomup": float(li[-1] - li[0]), "total_change_direct": float(lp[-1] - lp[0]), "max_abs_gap": float(np.max(np.abs(li - lp)))})173 out = pl.DataFrame(rows); out.write_csv(OUT / "aggregation_sensitivity.csv"); print(out)174175176def cmd_ablation(f):177 rows = []178 variants = {"full": {}, "noEB": {"eb": False}, "noTV": {"time_varying": False}, "LLT": {"model": "llt"}, "norobust": {"robust": False},179 "coarseloc": {"min_loc_n": 10**9}, "skipMRC": {"levels": [("region", "region_code"), ("municipality", "mun_code")]}}180 for cls in C.PROPERTY_TYPES:181 d = f.filter(pl.col("property_class") == cls)182 runs = {}183 for v, kw in variants.items():184 t0 = time.time(); runs[v] = run_sharp(d, GRID, **kw); print(cls, v, f"{time.time() - t0:.0f}s", flush=True)185 full = runs["full"]186 for v, r in runs.items():187 for lvl in ["region", "mrc", "municipality"]:188 keys = [k for k in r.cells if k[0] == lvl and k in full.cells]189 diffs, oos, cov, gains = [], [], [], []190 for k in keys:191 p, var = cell_log_index(r, *k); p0, _ = cell_log_index(full, *k)192 diffs.append(np.sqrt(np.mean((np.diff(p) - np.diff(p0)) ** 2)))193 c = r.cells[k]; z = c.oos_z[np.isfinite(c.oos_z)]194 if len(z) > 5:195 oos.append(np.sqrt(np.mean(z ** 2))); cov.append(np.mean(np.abs(z) < 1.96))196 gains.append(np.nanmean(c.fit.gain[c.n > 0]) if (c.n > 0).any() else np.nan)197 rows.append({"class": cls, "variant": v, "level": lvl, "cells": len(keys), "rmse_monthly_change_vs_full": float(np.mean(diffs)) if diffs else np.nan,198 "oos_innovation_rmse_z": float(np.mean(oos)) if oos else np.nan, "oos_innovation_share_within_95": float(np.mean(cov)) if cov else np.nan,199 "mean_gain": float(np.nanmean(gains))})200 out = pl.DataFrame(rows); out.write_csv(OUT / "ablation_real.csv"); print(out)201202203def cmd_trim(f):204 clean = pl.read_parquet(C.CLEAN_TX)205 from sharp_hpi.cleaning import analytic206 rows = []207 for cls in ["house", "condo"]:208 base = f.filter(pl.col("property_class") == cls)209 variants = {"headline": base}210 # include assessment-ratio outliers back211 inc = clean.with_columns(pl.when(pl.col("exclude_reason_headline") == "assessment_ratio_outlier").then(None).otherwise(pl.col("exclude_reason_headline")).alias("exclude_reason_headline")).with_columns(pl.col("exclude_reason_headline").is_null().alias("in_headline_sample"))212 variants["with_ratio_outliers"] = build_features(analytic(inc)).filter(pl.col("property_class") == cls)213 variants["drop_type_conflicts"] = base.filter(~pl.col("type_conflict"))214 variants["drop_new_construction"] = base.filter(~pl.col("new_construction"))215 q = base.group_by("month").agg(pl.col("price").quantile(0.01).alias("q1"), pl.col("price").quantile(0.99).alias("q99"))216 variants["trim_1_99_monthly"] = base.join(q, on="month").filter((pl.col("price") > pl.col("q1")) & (pl.col("price") < pl.col("q99")))217 variants["provider_only_classes"] = base.filter(pl.col("class_source") == "roll")218 ref = None219 for nm, dd in variants.items():220 s = E.pooled_td(assign_loc_fe(dd, min_n=15), GRID, robust=True).log_index; s_ols = E.pooled_td(assign_loc_fe(dd, min_n=15), GRID, robust=False).log_index221 if ref is None:222 ref = s223 rows.append({"class": cls, "variant": nm, "n": dd.height, "total_change_huber": float(s[-1] - s[0]), "total_change_ols": float(s_ols[-1] - s_ols[0]),224 "rmse_monthly_vs_headline": float(np.sqrt(np.nanmean((np.diff(s) - np.diff(ref)) ** 2))), "max_abs_level_gap": float(np.nanmax(np.abs(s - ref)))})225 out = pl.DataFrame(rows); out.write_csv(OUT / "trimming_sensitivity.csv"); print(out)226227228if __name__ == "__main__":229 cmd = sys.argv[1]230 f = load()231 {"repeat": cmd_repeat, "composition": cmd_composition, "placebo": cmd_placebo, "moran": cmd_moran, "bootstrap": cmd_bootstrap,232 "aggregation": cmd_aggregation, "ablation": cmd_ablation, "trim": cmd_trim}[cmd](f)233