#!/usr/bin/env python3 # ============================================================================= # QWHPI — Quebec Weekly Housing Price Index # Author : Simon-Pierre Boucher # Contact : contact@spboucher.ai # File : engine/scripts/10_monthly_validation.py # Purpose : v2.1 validation — repeat sales, downsampling (direct vs # hierarchical threshold), composition shock, drift & seasonality. # ============================================================================= """Monthly validation suite (v2.1).""" from __future__ import annotations import sys from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) import numpy as np import polars as pl from qwhpi.clean import CLEAN_PARQUET, research_sample from qwhpi.config import INTERIM_DIR, REPORTS_DIR, TABLES_DIR, ensure_dirs from qwhpi.features import build_features from qwhpi.hierarchy import cell_deviation from qwhpi.repeat_sales import bmn_index, build_pairs from qwhpi.rtd import local_time_dummy, month_grid from qwhpi.seasonal import seasonality_test from qwhpi.state_space import fit_local_level RNG = np.random.default_rng(42) LINES: list[str] = [] def note(line: str = "") -> None: LINES.append(line) print(line) def main() -> None: ensure_dirs() cl = research_sample(pl.read_parquet(CLEAN_PARQUET)) feat = build_features(cl) monthly = pl.read_parquet(INTERIM_DIR / "monthly_indexes.parquet") fp = pl.read_parquet(INTERIM_DIR / "rtd_residuals.parquet") grid = month_grid(sorted(fp["month"].unique().to_list())) note("# QHPI v2.1 validation report (monthly, robust)") note("") # ------------------------------------------------------------------ # # 1. Repeat-sales comparison # ------------------------------------------------------------------ # note("## 1. Repeat sales (BMN, unit-signature pairs) vs published index") note("") rows = [] for gid, mun, ptype in [ ("quebec", None, "unifamilial"), ("quebec", None, "condo"), ("quebec", None, "plex"), ("montreal", "Montréal", "condo"), ("montreal", "Montréal", "unifamilial"), ("quebec-city", "Québec", "condo"), ("quebec-city", "Québec", "unifamilial"), ]: sub = cl.filter(pl.col("propertyType") == ptype) if mun: sub = sub.filter(pl.col("municipality") == mun) rs = bmn_index(build_pairs(sub)) if rs is None: continue h = monthly.filter((pl.col("geography_id") == gid) & (pl.col("property_type") == ptype)) \ .select(pl.col("period").alias("month"), "index_research") j = rs.join(h, on="month", how="inner").sort("month") lg_rs = np.log(j["rs_index"].to_numpy()) lg_he = np.log(j["index_research"].to_numpy()) corr = float(np.corrcoef(np.diff(lg_rs), np.diff(lg_he))[0, 1]) gap = float((lg_he[-2] - lg_rs[-2]) * 100) rows.append({"cell": f"{gid}/{ptype}", "n_pairs": build_pairs(sub).height, "corr_monthly_changes": round(corr, 3), "end_gap_pct": round(gap, 2)}) note(f"- {gid}/{ptype}: Δcorr={corr:.3f}, écart final hedonic−RS = {gap:+.1f}%") pl.DataFrame(rows).write_csv(TABLES_DIR / "repeat_sales_comparison_monthly.csv") # ------------------------------------------------------------------ # # 2. Downsampling — validates LIQUID_MIN and the grade scale # ------------------------------------------------------------------ # note("") note("## 2. Downsampling (Montréal condo, ~620 tx/month)") note("") mc_feat = feat.filter((pl.col("municipality") == "Montréal") & (pl.col("propertyType") == "condo")) ref = local_time_dummy(mc_feat, grid) assert ref is not None ref_path = ref[0] - ref[0].mean() mc_resid = fp.filter((pl.col("municipality") == "Montréal") & (pl.col("propertyType") == "condo")) \ .rename({"month": "week_str"}) sig2 = float(mc_resid["resid"].var()) prov_path_tbl = pl.read_parquet(INTERIM_DIR / "rtd_paths.parquet") prov = prov_path_tbl.filter(pl.col("property_type") == "condo").sort("month") prov_path = prov["mean_splice"].to_numpy() rows = [] for target, n_reps in [(300, 8), (150, 8), (75, 8), (40, 10), (20, 10), (10, 10)]: rmse_d, rmse_h = [], [] for _ in range(n_reps): thin_f = (mc_feat.with_columns(pl.lit(RNG.random(mc_feat.height)).alias("_u")) .with_columns(pl.col("_u").rank().over("month").alias("_rk")) .filter(pl.col("_rk") <= target)) # direct estimator on thinned cell d = local_time_dummy(thin_f, grid) if d is not None: sm = fit_local_level(d[0], np.maximum(d[1], 1.0), d[2]).smoothed err = (sm - sm.mean()) - ref_path rmse_d.append(np.sqrt(np.mean(err ** 2))) # hierarchical estimator on thinned residuals thin_r = (mc_resid.with_columns(pl.lit(RNG.random(mc_resid.height)).alias("_u")) .with_columns(pl.col("_u").rank().over("week_str").alias("_rk")) .filter(pl.col("_rk") <= target)) dev = cell_deviation(thin_r.rename({"resid": "dev"}) .select("week_str", "dev"), grid, sig2) hier = prov_path + dev.fit.smoothed err = (hier - hier.mean()) - ref_path rmse_h.append(np.sqrt(np.mean(err ** 2))) rows.append({"target_tx_per_month": target, "rmse_direct_pct": round(float(np.mean(rmse_d)) * 100, 2) if rmse_d else None, "rmse_hierarchical_pct": round(float(np.mean(rmse_h)) * 100, 2)}) note(f"- {target:>3} tx/mois: RMSE direct={rows[-1]['rmse_direct_pct']}% " f"| hiérarchique={rows[-1]['rmse_hierarchical_pct']}%") pl.DataFrame(rows).write_csv(TABLES_DIR / "downsampling_monthly.csv") note("") note("Lecture: en RMSE pur (référence = direct plein échantillon, cellule à " "faible divergence) le croisement se situe vers 75–150 tx/mois. Mais le " "RMSE du hiérarchique contient un biais de compression qui explose sur " "les cellules fortement divergentes (Québec-ville condo: ~8–10 pts, cf. " "arbitrage RS/stratifié) alors que l'erreur du direct est du bruit pur " "(~1.9% à 40 tx/mois), réduit par le lisseur d'état. LIQUID_MIN=40 " "échange donc un bruit borné contre un biais non borné.") # ------------------------------------------------------------------ # # 3. Composition shock (monthly) # ------------------------------------------------------------------ # note("") note("## 3. Choc de composition (province unifamilial, mensuel)") import pyfixest as pf uni = feat.filter(pl.col("propertyType") == "unifamilial") med_fa = float(uni["floor_area_filled"].median()) shock_months = [m for m in grid if m.startswith("2024")][:6] u = RNG.random(uni.height) uni = uni.with_columns(pl.lit(u).alias("_u")) shocked = uni.filter(~(pl.col("month").is_in(shock_months) & (pl.col("floor_area_filled") < med_fa) & (pl.col("_u") < 0.7))) def med_path(df): return {r["month"]: r["m"] for r in df.group_by("month").agg(pl.col("amount").median().alias("m")) .iter_rows(named=True)} def hed_path(df): pdf = df.select("log_amount", "month", "log_fa", "fa_missing", "age_bin", "loc_fine").to_pandas() pdf["fa_missing"] = pdf["fa_missing"].astype(int) fit = pf.feols("log_amount ~ log_fa + fa_missing + C(age_bin) " "| loc_fine + month", data=pdf) fx = fit.fixef() key = next(k for k in fx if "month" in k) return {k: float(v) for k, v in fx[key].items()} mb, ms = med_path(uni), med_path(shocked) hb, hs = hed_path(uni), hed_path(shocked) common = [m for m in shock_months if m in mb and m in ms and m in hb and m in hs] dmed = float(np.mean([np.log(ms[m] / mb[m]) for m in common])) * 100 dhed = float(np.mean([hs[m] - hb[m] for m in common])) * 100 note(f"- Médiane brute: {dmed:+.2f}% dans les mois choqués (doit bouger)") note(f"- Hédonique: {dhed:+.2f}% (doit rester ~0)") pl.DataFrame({"metric": ["raw_median_shift_pct", "hedonic_shift_pct"], "value": [round(dmed, 3), round(dhed, 3)]}) \ .write_csv(TABLES_DIR / "composition_shock_monthly.csv") # ------------------------------------------------------------------ # # 4. Drift, splice, seasonality # ------------------------------------------------------------------ # note("") note("## 4. Diagnostics de robustesse") drift = pl.read_csv(TABLES_DIR / "rtd_coefficient_drift.csv") b = drift["beta_log_fa"] note(f"- Dérive β(log surface) sur 55 fenêtres: {b.min():.3f} → {b.max():.3f} " f"(amplitude {(b.max() - b.min()) / b.mean() * 100:.1f}% de la moyenne) — " "justifie le RTD contre le pool fixe.") seas_rows = [] for gid, t in [("quebec", "all"), ("quebec", "unifamilial"), ("montreal", "condo")]: s = monthly.filter((pl.col("geography_id") == gid) & (pl.col("property_type") == t)).sort("period") res = seasonality_test(np.log(s["index_research"].to_numpy()), n_harmonics=2) seas_rows.append({"cell": f"{gid}/{t}", **res}) note(f"- Saisonnalité {gid}/{t}: F={res['f_stat']} p={res['p_value']} → " + ("significative" if res["significant_5pct"] else "non significative (NSA)")) pl.DataFrame(seas_rows).write_csv(TABLES_DIR / "seasonality_monthly.csv") header = ( "\n\n" ) (REPORTS_DIR / "validation_v2_report.md").write_text( header + "\n".join(LINES) + "\n", encoding="utf-8") print("\nReport: outputs/reports/validation_v2_report.md") if __name__ == "__main__": main()