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/06_validation.py7# Purpose : Step 6 — repeat-sales comparison, downsampling experiment,8#           composition-shock simulation, weekly-vs-monthly comparison.9# =============================================================================10"""Validation suite (Execution Order step 6). Writes tables + report."""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 REPORTS_DIR, TABLES_DIR, ensure_dirs24from qwhpi.features import build_features25from qwhpi.hierarchy import cell_deviation, week_grid26from qwhpi.repeat_sales import bmn_index, build_pairs2728RNG = np.random.default_rng(42)29LINES: list[str] = []303132def note(line: str = "") -> None:33    LINES.append(line)34    print(line)353637# ------------------------------------------------------------------------- #38# 1. Repeat-sales comparison39# ------------------------------------------------------------------------- #4041def repeat_sales_comparison(cl: pl.DataFrame, hier: pl.DataFrame) -> None:42    note("## 1. Repeat-sales (BMN) vs hedonic index")43    note("")44    cells = [45        ("quebec", None, "unifamilial"), ("quebec", None, "condo"),46        ("quebec", None, "plex"),47        ("montreal", "Montréal", "condo"), ("montreal", "Montréal", "unifamilial"),48        ("quebec-city", "Québec", "unifamilial"),49    ]50    rows = []51    for gid, mun, ptype in cells:52        sub = cl.filter(pl.col("propertyType") == ptype)53        if mun:54            sub = sub.filter(pl.col("municipality") == mun)55        rs = bmn_index(build_pairs(sub))56        if rs is None:57            continue58        h = (59            hier.filter((pl.col("geography_id") == gid)60                        & (pl.col("property_type") == ptype))61            .with_columns(pl.col("week").str.slice(0, 7).alias("month"))62            .group_by("month").agg(pl.col("index_research").mean().alias("hedonic"))63        )64        j = rs.join(h, on="month", how="inner").sort("month")65        lg_rs = np.log(j["rs_index"].to_numpy())66        lg_he = np.log(j["hedonic"].to_numpy())67        corr = float(np.corrcoef(np.diff(lg_rs), np.diff(lg_he))[0, 1])68        rows.append({69            "cell": f"{gid}/{ptype}", "n_pairs": build_pairs(sub).height,70            "corr_monthly_changes": round(corr, 3),71            "end_gap_pct": round((lg_he[-1] - lg_rs[-1]) * 100, 2),72            "rs_end": round(float(j["rs_index"][-1]), 1),73            "hedonic_end": round(float(j["hedonic"][-1]), 1),74        })75        note(f"- {gid}/{ptype}: Δcorr={corr:.3f}, end levels RS={j['rs_index'][-1]:.1f} "76             f"vs hedonic={j['hedonic'][-1]:.1f}")77    pl.DataFrame(rows).write_csv(TABLES_DIR / "repeat_sales_comparison.csv")787980# ------------------------------------------------------------------------- #81# 2. Downsampling experiment (Montreal condo lab)82# ------------------------------------------------------------------------- #8384def downsampling(resid: pl.DataFrame, n_reps: int = 20) -> None:85    note("")86    note("## 2. Downsampling experiment — Montréal condo lab")87    note("")88    cell = resid.filter((pl.col("municipality") == "Montréal")89                        & (pl.col("propertyType") == "condo"))90    grid = week_grid(resid)91    sigma2 = float(cell["resid"].var())9293    ref = cell_deviation(cell.rename({"resid": "dev"}).select("week_str", "dev"),94                         grid, sigma2)95    ref_path = ref.fit.smoothed96    ref_turn = np.sign(np.diff(_ma(ref_path, 4)))9798    rows = []99    for target in [100, 50, 25, 15, 10, 5]:100        rmse_l, bias_l, vol_l, turn_l, cov_l = [], [], [], [], []101        for _ in range(n_reps):102            thin = (103                cell.with_columns(pl.lit(RNG.random(cell.height)).alias("_u"))104                .with_columns(pl.col("_u").rank().over("week_str").alias("_rk"))105                .filter(pl.col("_rk") <= target)106            )107            fit = cell_deviation(thin.rename({"resid": "dev"}).select("week_str", "dev"),108                                 grid, sigma2)109            path = fit.fit.smoothed110            err = path - ref_path111            rmse_l.append(np.sqrt(np.mean(err ** 2)))112            bias_l.append(np.mean(err))113            vol_l.append(np.std(np.diff(path)) / max(np.std(np.diff(ref_path)), 1e-12))114            turn = np.sign(np.diff(_ma(path, 4)))115            turn_l.append(float(np.mean(turn == ref_turn)))116            band = 1.96 * np.sqrt(fit.fit.smoothed_var)117            cov_l.append(float(np.mean(np.abs(err) <= band)))118        rows.append({119            "target_tx_per_week": target,120            "rmse_log_pct": round(float(np.mean(rmse_l)) * 100, 3),121            "bias_log_pct": round(float(np.mean(bias_l)) * 100, 3),122            "volatility_ratio": round(float(np.mean(vol_l)), 3),123            "turning_point_agreement": round(float(np.mean(turn_l)), 3),124            "ci95_coverage": round(float(np.mean(cov_l)), 3),125        })126        note(f"- {target:>3} tx/wk: RMSE={rows[-1]['rmse_log_pct']}% "127             f"bias={rows[-1]['bias_log_pct']}% vol_ratio={rows[-1]['volatility_ratio']} "128             f"turn_agree={rows[-1]['turning_point_agreement']} "129             f"CI95_cov={rows[-1]['ci95_coverage']}")130    pl.DataFrame(rows).write_csv(TABLES_DIR / "downsampling_results.csv")131    note("")132    note("Reading: RMSE vs the full-sample (~145 tx/wk) smoothed path; CI "133         "coverage of the deviation posterior; turning points on 4-week MA.")134135136def _ma(x: np.ndarray, w: int) -> np.ndarray:137    return np.convolve(x, np.ones(w) / w, mode="same")138139140# ------------------------------------------------------------------------- #141# 3. Composition-shock simulation142# ------------------------------------------------------------------------- #143144def composition_shock(feat: pl.DataFrame) -> None:145    note("")146    note("## 3. Composition-shock simulation (province unifamilial)")147    note("")148    import pyfixest as pf149150    uni = feat.filter(pl.col("propertyType") == "unifamilial")151    med_fa = float(uni["floor_area_filled"].median())152    shock_weeks = [w for w in uni["week_str"].unique().to_list()153                   if w.startswith("2024") and w < "2024-07"]154155    # Drop 70% of below-median-floorArea sales in shocked weeks.156    u = RNG.random(uni.height)157    uni = uni.with_columns(pl.lit(u).alias("_u"))158    shocked = uni.filter(159        ~(160            pl.col("week_str").is_in(shock_weeks)161            & (pl.col("floor_area_filled") < med_fa)162            & (pl.col("_u") < 0.7)163        )164    )165    note(f"- Shock: drop 70% of below-median floorArea sales in 2024-H1 "166         f"({uni.height - shocked.height:,} of {uni.height:,} rows removed)")167168    def weekly_median(df: pl.DataFrame) -> dict[str, float]:169        return {r["week_str"]: r["m"] for r in170                df.group_by("week_str").agg(pl.col("amount").median().alias("m"))171                .iter_rows(named=True)}172173    def hedonic_path(df: pl.DataFrame) -> dict[str, float]:174        pdf = df.select("log_amount", "week_str", "log_fa", "fa_missing",175                        "age_bin", "building_type", "loc_fine").to_pandas()176        pdf["fa_missing"] = pdf["fa_missing"].astype(int)177        fit = pf.feols("log_amount ~ log_fa + fa_missing + C(age_bin) "178                       "| loc_fine + week_str", data=pdf)179        fx = fit.fixef()180        key = next(k for k in fx if "week_str" in k)181        return {k: float(v) for k, v in fx[key].items()}182183    med_b, med_s = weekly_median(uni), weekly_median(shocked)184    hed_b, hed_s = hedonic_path(uni), hedonic_path(shocked)185186    common = [w for w in shock_weeks if w in med_b and w in med_s187              and w in hed_b and w in hed_s]188    dmed = np.mean([np.log(med_s[w] / med_b[w]) for w in common]) * 100189    dhed = np.mean([hed_s[w] - hed_b[w] for w in common]) * 100190    note(f"- Raw weekly median moves by {dmed:+.2f}% in shocked weeks (must move)")191    note(f"- Hedonic time-dummy moves by {dhed:+.2f}% (must stay ~0)")192    pl.DataFrame({193        "metric": ["raw_median_shift_pct", "hedonic_shift_pct"],194        "value": [round(float(dmed), 3), round(float(dhed), 3)],195    }).write_csv(TABLES_DIR / "composition_shock.csv")196197198# ------------------------------------------------------------------------- #199# 4. Weekly vs monthly comparison200# ------------------------------------------------------------------------- #201202def weekly_vs_monthly(hier: pl.DataFrame) -> None:203    note("")204    note("## 4. Weekly vs monthly (same methodology)")205    note("")206    rows = []207    for gid, ptype in [("quebec", "unifamilial"), ("quebec", "condo"),208                       ("montreal", "condo")]:209        h = hier.filter((pl.col("geography_id") == gid)210                        & (pl.col("property_type") == ptype)).sort("week")211        wk_raw = np.log(h["index"].to_numpy())212        wk_smooth = np.log(h["index_smoothed"].to_numpy())213        wk_res = np.log(h["index_research"].to_numpy())214        monthly = (215            h.with_columns(pl.col("week").str.slice(0, 7).alias("m"))216            .group_by("m").agg(pl.col("index_research").mean().alias("ix"))217            .sort("m")218        )219        mo = np.log(monthly["ix"].to_numpy())220        # noise = raw minus final smoothed; signal = final smoothed changes221        noise_sd = float(np.nanstd(wk_raw - wk_res))222        signal_sd = float(np.std(np.diff(wk_res)))223        revision_sd = float(np.std(wk_smooth - wk_res))  # real-time vs final224        rows.append({225            "cell": f"{gid}/{ptype}",226            "weekly_signal_sd_pct": round(signal_sd * 100, 3),227            "weekly_noise_sd_pct": round(noise_sd * 100, 3),228            "monthly_change_sd_pct": round(float(np.std(np.diff(mo))) * 100, 3),229            "realtime_revision_sd_pct": round(revision_sd * 100, 3),230            "weeks_to_monthly_ratio": round(signal_sd * np.sqrt(52 / 12)231                                            / max(float(np.std(np.diff(mo))), 1e-9), 2),232        })233        note(f"- {gid}/{ptype}: weekly signal SD={rows[-1]['weekly_signal_sd_pct']}%, "234             f"noise SD={rows[-1]['weekly_noise_sd_pct']}%, monthly ΔSD="235             f"{rows[-1]['monthly_change_sd_pct']}%, real-time revision SD="236             f"{rows[-1]['realtime_revision_sd_pct']}%")237    pl.DataFrame(rows).write_csv(TABLES_DIR / "weekly_vs_monthly.csv")238239240def main() -> None:241    ensure_dirs()242    cl = research_sample(pl.read_parquet(CLEAN_PARQUET))243    feat = build_features(cl)244    hier = pl.read_parquet("data/interim/hierarchical_indexes.parquet")245    resid = pl.read_parquet("data/interim/stage1_residuals.parquet")246247    note("# QWHPI validation report")248    note("")249    repeat_sales_comparison(cl, hier)250    downsampling(resid)251    composition_shock(feat)252    weekly_vs_monthly(hier)253254    header = (255        "<!--\n"256        "=============================================================================\n"257        "QWHPI — Quebec Weekly Housing Price Index\n"258        "Author  : Simon-Pierre Boucher\n"259        "Contact : contact@spboucher.ai\n"260        "File    : outputs/reports/validation_report.md\n"261        "Purpose : Step 6 validation report (generated by 06_validation.py)\n"262        "=============================================================================\n"263        "-->\n\n"264    )265    (REPORTS_DIR / "validation_report.md").write_text(266        header + "\n".join(LINES) + "\n", encoding="utf-8")267    print("\nReport written to outputs/reports/validation_report.md")268269270if __name__ == "__main__":271    main()272