spb/qwhpi Public

1#!/usr/bin/env python32# =============================================================================3# QWHPI — Quebec Weekly Housing Price Index4# Author  : Simon-Pierre Boucher5# Contact : contact@spboucher.ai6# File    : engine/scripts/10_monthly_validation.py7# Purpose : v2.1 validation — repeat sales, downsampling (direct vs8#           hierarchical threshold), composition shock, drift & seasonality.9# =============================================================================10"""Monthly validation suite (v2.1)."""1112from __future__ import annotations1314import sys15from pathlib import Path1617sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))1819import numpy as np20import polars as pl2122from qwhpi.clean import CLEAN_PARQUET, research_sample23from qwhpi.config import INTERIM_DIR, REPORTS_DIR, TABLES_DIR, ensure_dirs24from qwhpi.features import build_features25from qwhpi.hierarchy import cell_deviation26from qwhpi.repeat_sales import bmn_index, build_pairs27from qwhpi.rtd import local_time_dummy, month_grid28from qwhpi.seasonal import seasonality_test29from qwhpi.state_space import fit_local_level3031RNG = np.random.default_rng(42)32LINES: list[str] = []333435def note(line: str = "") -> None:36    LINES.append(line)37    print(line)383940def main() -> None:41    ensure_dirs()42    cl = research_sample(pl.read_parquet(CLEAN_PARQUET))43    feat = build_features(cl)44    monthly = pl.read_parquet(INTERIM_DIR / "monthly_indexes.parquet")45    fp = pl.read_parquet(INTERIM_DIR / "rtd_residuals.parquet")46    grid = month_grid(sorted(fp["month"].unique().to_list()))4748    note("# QHPI v2.1 validation report (monthly, robust)")49    note("")5051    # ------------------------------------------------------------------ #52    # 1. Repeat-sales comparison53    # ------------------------------------------------------------------ #54    note("## 1. Repeat sales (BMN, unit-signature pairs) vs published index")55    note("")56    rows = []57    for gid, mun, ptype in [58        ("quebec", None, "unifamilial"), ("quebec", None, "condo"),59        ("quebec", None, "plex"), ("montreal", "Montréal", "condo"),60        ("montreal", "Montréal", "unifamilial"),61        ("quebec-city", "Québec", "condo"),62        ("quebec-city", "Québec", "unifamilial"),63    ]:64        sub = cl.filter(pl.col("propertyType") == ptype)65        if mun:66            sub = sub.filter(pl.col("municipality") == mun)67        rs = bmn_index(build_pairs(sub))68        if rs is None:69            continue70        h = monthly.filter((pl.col("geography_id") == gid)71                           & (pl.col("property_type") == ptype)) \72            .select(pl.col("period").alias("month"), "index_research")73        j = rs.join(h, on="month", how="inner").sort("month")74        lg_rs = np.log(j["rs_index"].to_numpy())75        lg_he = np.log(j["index_research"].to_numpy())76        corr = float(np.corrcoef(np.diff(lg_rs), np.diff(lg_he))[0, 1])77        gap = float((lg_he[-2] - lg_rs[-2]) * 100)78        rows.append({"cell": f"{gid}/{ptype}", "n_pairs": build_pairs(sub).height,79                     "corr_monthly_changes": round(corr, 3),80                     "end_gap_pct": round(gap, 2)})81        note(f"- {gid}/{ptype}: Δcorr={corr:.3f}, écart final hedonic−RS = {gap:+.1f}%")82    pl.DataFrame(rows).write_csv(TABLES_DIR / "repeat_sales_comparison_monthly.csv")8384    # ------------------------------------------------------------------ #85    # 2. Downsampling — validates LIQUID_MIN and the grade scale86    # ------------------------------------------------------------------ #87    note("")88    note("## 2. Downsampling (Montréal condo, ~620 tx/month)")89    note("")90    mc_feat = feat.filter((pl.col("municipality") == "Montréal")91                          & (pl.col("propertyType") == "condo"))92    ref = local_time_dummy(mc_feat, grid)93    assert ref is not None94    ref_path = ref[0] - ref[0].mean()9596    mc_resid = fp.filter((pl.col("municipality") == "Montréal")97                         & (pl.col("propertyType") == "condo")) \98        .rename({"month": "week_str"})99    sig2 = float(mc_resid["resid"].var())100    prov_path_tbl = pl.read_parquet(INTERIM_DIR / "rtd_paths.parquet")101    prov = prov_path_tbl.filter(pl.col("property_type") == "condo").sort("month")102    prov_path = prov["mean_splice"].to_numpy()103104    rows = []105    for target, n_reps in [(300, 8), (150, 8), (75, 8), (40, 10), (20, 10), (10, 10)]:106        rmse_d, rmse_h = [], []107        for _ in range(n_reps):108            thin_f = (mc_feat.with_columns(pl.lit(RNG.random(mc_feat.height)).alias("_u"))109                      .with_columns(pl.col("_u").rank().over("month").alias("_rk"))110                      .filter(pl.col("_rk") <= target))111            # direct estimator on thinned cell112            d = local_time_dummy(thin_f, grid)113            if d is not None:114                sm = fit_local_level(d[0], np.maximum(d[1], 1.0), d[2]).smoothed115                err = (sm - sm.mean()) - ref_path116                rmse_d.append(np.sqrt(np.mean(err ** 2)))117            # hierarchical estimator on thinned residuals118            thin_r = (mc_resid.with_columns(pl.lit(RNG.random(mc_resid.height)).alias("_u"))119                      .with_columns(pl.col("_u").rank().over("week_str").alias("_rk"))120                      .filter(pl.col("_rk") <= target))121            dev = cell_deviation(thin_r.rename({"resid": "dev"})122                                 .select("week_str", "dev"), grid, sig2)123            hier = prov_path + dev.fit.smoothed124            err = (hier - hier.mean()) - ref_path125            rmse_h.append(np.sqrt(np.mean(err ** 2)))126        rows.append({"target_tx_per_month": target,127                     "rmse_direct_pct": round(float(np.mean(rmse_d)) * 100, 2) if rmse_d else None,128                     "rmse_hierarchical_pct": round(float(np.mean(rmse_h)) * 100, 2)})129        note(f"- {target:>3} tx/mois: RMSE direct={rows[-1]['rmse_direct_pct']}% "130             f"| hiérarchique={rows[-1]['rmse_hierarchical_pct']}%")131    pl.DataFrame(rows).write_csv(TABLES_DIR / "downsampling_monthly.csv")132    note("")133    note("Lecture: en RMSE pur (référence = direct plein échantillon, cellule à "134         "faible divergence) le croisement se situe vers 75–150 tx/mois. Mais le "135         "RMSE du hiérarchique contient un biais de compression qui explose sur "136         "les cellules fortement divergentes (Québec-ville condo: ~8–10 pts, cf. "137         "arbitrage RS/stratifié) alors que l'erreur du direct est du bruit pur "138         "(~1.9% à 40 tx/mois), réduit par le lisseur d'état. LIQUID_MIN=40 "139         "échange donc un bruit borné contre un biais non borné.")140141    # ------------------------------------------------------------------ #142    # 3. Composition shock (monthly)143    # ------------------------------------------------------------------ #144    note("")145    note("## 3. Choc de composition (province unifamilial, mensuel)")146    import pyfixest as pf147    uni = feat.filter(pl.col("propertyType") == "unifamilial")148    med_fa = float(uni["floor_area_filled"].median())149    shock_months = [m for m in grid if m.startswith("2024")][:6]150    u = RNG.random(uni.height)151    uni = uni.with_columns(pl.lit(u).alias("_u"))152    shocked = uni.filter(~(pl.col("month").is_in(shock_months)153                           & (pl.col("floor_area_filled") < med_fa)154                           & (pl.col("_u") < 0.7)))155156    def med_path(df):157        return {r["month"]: r["m"] for r in158                df.group_by("month").agg(pl.col("amount").median().alias("m"))159                .iter_rows(named=True)}160161    def hed_path(df):162        pdf = df.select("log_amount", "month", "log_fa", "fa_missing",163                        "age_bin", "loc_fine").to_pandas()164        pdf["fa_missing"] = pdf["fa_missing"].astype(int)165        fit = pf.feols("log_amount ~ log_fa + fa_missing + C(age_bin) "166                       "| loc_fine + month", data=pdf)167        fx = fit.fixef()168        key = next(k for k in fx if "month" in k)169        return {k: float(v) for k, v in fx[key].items()}170171    mb, ms = med_path(uni), med_path(shocked)172    hb, hs = hed_path(uni), hed_path(shocked)173    common = [m for m in shock_months if m in mb and m in ms and m in hb and m in hs]174    dmed = float(np.mean([np.log(ms[m] / mb[m]) for m in common])) * 100175    dhed = float(np.mean([hs[m] - hb[m] for m in common])) * 100176    note(f"- Médiane brute: {dmed:+.2f}% dans les mois choqués (doit bouger)")177    note(f"- Hédonique: {dhed:+.2f}% (doit rester ~0)")178    pl.DataFrame({"metric": ["raw_median_shift_pct", "hedonic_shift_pct"],179                  "value": [round(dmed, 3), round(dhed, 3)]}) \180        .write_csv(TABLES_DIR / "composition_shock_monthly.csv")181182    # ------------------------------------------------------------------ #183    # 4. Drift, splice, seasonality184    # ------------------------------------------------------------------ #185    note("")186    note("## 4. Diagnostics de robustesse")187    drift = pl.read_csv(TABLES_DIR / "rtd_coefficient_drift.csv")188    b = drift["beta_log_fa"]189    note(f"- Dérive β(log surface) sur 55 fenêtres: {b.min():.3f} → {b.max():.3f} "190         f"(amplitude {(b.max() - b.min()) / b.mean() * 100:.1f}% de la moyenne) — "191         "justifie le RTD contre le pool fixe.")192    seas_rows = []193    for gid, t in [("quebec", "all"), ("quebec", "unifamilial"), ("montreal", "condo")]:194        s = monthly.filter((pl.col("geography_id") == gid)195                           & (pl.col("property_type") == t)).sort("period")196        res = seasonality_test(np.log(s["index_research"].to_numpy()), n_harmonics=2)197        seas_rows.append({"cell": f"{gid}/{t}", **res})198        note(f"- Saisonnalité {gid}/{t}: F={res['f_stat']} p={res['p_value']} → "199             + ("significative" if res["significant_5pct"] else "non significative (NSA)"))200    pl.DataFrame(seas_rows).write_csv(TABLES_DIR / "seasonality_monthly.csv")201202    header = (203        "<!--\n"204        "=============================================================================\n"205        "QWHPI — Quebec Weekly Housing Price Index\n"206        "Author  : Simon-Pierre Boucher\n"207        "Contact : contact@spboucher.ai\n"208        "File    : outputs/reports/validation_v2_report.md\n"209        "Purpose : v2.1 monthly validation report (10_monthly_validation.py)\n"210        "=============================================================================\n"211        "-->\n\n"212    )213    (REPORTS_DIR / "validation_v2_report.md").write_text(214        header + "\n".join(LINES) + "\n", encoding="utf-8")215    print("\nReport: outputs/reports/validation_v2_report.md")216217218if __name__ == "__main__":219    main()220