#!/usr/bin/env python3 # ============================================================================= # QWHPI — Quebec Weekly Housing Price Index # Author : Simon-Pierre Boucher # Contact : contact@spboucher.ai # File : engine/scripts/06_validation.py # Purpose : Step 6 — repeat-sales comparison, downsampling experiment, # composition-shock simulation, weekly-vs-monthly comparison. # ============================================================================= """Validation suite (Execution Order step 6). Writes tables + report.""" 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 REPORTS_DIR, TABLES_DIR, ensure_dirs from qwhpi.features import build_features from qwhpi.hierarchy import cell_deviation, week_grid from qwhpi.repeat_sales import bmn_index, build_pairs RNG = np.random.default_rng(42) LINES: list[str] = [] def note(line: str = "") -> None: LINES.append(line) print(line) # ------------------------------------------------------------------------- # # 1. Repeat-sales comparison # ------------------------------------------------------------------------- # def repeat_sales_comparison(cl: pl.DataFrame, hier: pl.DataFrame) -> None: note("## 1. Repeat-sales (BMN) vs hedonic index") note("") cells = [ ("quebec", None, "unifamilial"), ("quebec", None, "condo"), ("quebec", None, "plex"), ("montreal", "Montréal", "condo"), ("montreal", "Montréal", "unifamilial"), ("quebec-city", "Québec", "unifamilial"), ] rows = [] for gid, mun, ptype in cells: 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 = ( hier.filter((pl.col("geography_id") == gid) & (pl.col("property_type") == ptype)) .with_columns(pl.col("week").str.slice(0, 7).alias("month")) .group_by("month").agg(pl.col("index_research").mean().alias("hedonic")) ) j = rs.join(h, on="month", how="inner").sort("month") lg_rs = np.log(j["rs_index"].to_numpy()) lg_he = np.log(j["hedonic"].to_numpy()) corr = float(np.corrcoef(np.diff(lg_rs), np.diff(lg_he))[0, 1]) rows.append({ "cell": f"{gid}/{ptype}", "n_pairs": build_pairs(sub).height, "corr_monthly_changes": round(corr, 3), "end_gap_pct": round((lg_he[-1] - lg_rs[-1]) * 100, 2), "rs_end": round(float(j["rs_index"][-1]), 1), "hedonic_end": round(float(j["hedonic"][-1]), 1), }) note(f"- {gid}/{ptype}: Δcorr={corr:.3f}, end levels RS={j['rs_index'][-1]:.1f} " f"vs hedonic={j['hedonic'][-1]:.1f}") pl.DataFrame(rows).write_csv(TABLES_DIR / "repeat_sales_comparison.csv") # ------------------------------------------------------------------------- # # 2. Downsampling experiment (Montreal condo lab) # ------------------------------------------------------------------------- # def downsampling(resid: pl.DataFrame, n_reps: int = 20) -> None: note("") note("## 2. Downsampling experiment — Montréal condo lab") note("") cell = resid.filter((pl.col("municipality") == "Montréal") & (pl.col("propertyType") == "condo")) grid = week_grid(resid) sigma2 = float(cell["resid"].var()) ref = cell_deviation(cell.rename({"resid": "dev"}).select("week_str", "dev"), grid, sigma2) ref_path = ref.fit.smoothed ref_turn = np.sign(np.diff(_ma(ref_path, 4))) rows = [] for target in [100, 50, 25, 15, 10, 5]: rmse_l, bias_l, vol_l, turn_l, cov_l = [], [], [], [], [] for _ in range(n_reps): thin = ( cell.with_columns(pl.lit(RNG.random(cell.height)).alias("_u")) .with_columns(pl.col("_u").rank().over("week_str").alias("_rk")) .filter(pl.col("_rk") <= target) ) fit = cell_deviation(thin.rename({"resid": "dev"}).select("week_str", "dev"), grid, sigma2) path = fit.fit.smoothed err = path - ref_path rmse_l.append(np.sqrt(np.mean(err ** 2))) bias_l.append(np.mean(err)) vol_l.append(np.std(np.diff(path)) / max(np.std(np.diff(ref_path)), 1e-12)) turn = np.sign(np.diff(_ma(path, 4))) turn_l.append(float(np.mean(turn == ref_turn))) band = 1.96 * np.sqrt(fit.fit.smoothed_var) cov_l.append(float(np.mean(np.abs(err) <= band))) rows.append({ "target_tx_per_week": target, "rmse_log_pct": round(float(np.mean(rmse_l)) * 100, 3), "bias_log_pct": round(float(np.mean(bias_l)) * 100, 3), "volatility_ratio": round(float(np.mean(vol_l)), 3), "turning_point_agreement": round(float(np.mean(turn_l)), 3), "ci95_coverage": round(float(np.mean(cov_l)), 3), }) note(f"- {target:>3} tx/wk: RMSE={rows[-1]['rmse_log_pct']}% " f"bias={rows[-1]['bias_log_pct']}% vol_ratio={rows[-1]['volatility_ratio']} " f"turn_agree={rows[-1]['turning_point_agreement']} " f"CI95_cov={rows[-1]['ci95_coverage']}") pl.DataFrame(rows).write_csv(TABLES_DIR / "downsampling_results.csv") note("") note("Reading: RMSE vs the full-sample (~145 tx/wk) smoothed path; CI " "coverage of the deviation posterior; turning points on 4-week MA.") def _ma(x: np.ndarray, w: int) -> np.ndarray: return np.convolve(x, np.ones(w) / w, mode="same") # ------------------------------------------------------------------------- # # 3. Composition-shock simulation # ------------------------------------------------------------------------- # def composition_shock(feat: pl.DataFrame) -> None: note("") note("## 3. Composition-shock simulation (province unifamilial)") note("") import pyfixest as pf uni = feat.filter(pl.col("propertyType") == "unifamilial") med_fa = float(uni["floor_area_filled"].median()) shock_weeks = [w for w in uni["week_str"].unique().to_list() if w.startswith("2024") and w < "2024-07"] # Drop 70% of below-median-floorArea sales in shocked weeks. u = RNG.random(uni.height) uni = uni.with_columns(pl.lit(u).alias("_u")) shocked = uni.filter( ~( pl.col("week_str").is_in(shock_weeks) & (pl.col("floor_area_filled") < med_fa) & (pl.col("_u") < 0.7) ) ) note(f"- Shock: drop 70% of below-median floorArea sales in 2024-H1 " f"({uni.height - shocked.height:,} of {uni.height:,} rows removed)") def weekly_median(df: pl.DataFrame) -> dict[str, float]: return {r["week_str"]: r["m"] for r in df.group_by("week_str").agg(pl.col("amount").median().alias("m")) .iter_rows(named=True)} def hedonic_path(df: pl.DataFrame) -> dict[str, float]: pdf = df.select("log_amount", "week_str", "log_fa", "fa_missing", "age_bin", "building_type", "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 + week_str", data=pdf) fx = fit.fixef() key = next(k for k in fx if "week_str" in k) return {k: float(v) for k, v in fx[key].items()} med_b, med_s = weekly_median(uni), weekly_median(shocked) hed_b, hed_s = hedonic_path(uni), hedonic_path(shocked) common = [w for w in shock_weeks if w in med_b and w in med_s and w in hed_b and w in hed_s] dmed = np.mean([np.log(med_s[w] / med_b[w]) for w in common]) * 100 dhed = np.mean([hed_s[w] - hed_b[w] for w in common]) * 100 note(f"- Raw weekly median moves by {dmed:+.2f}% in shocked weeks (must move)") note(f"- Hedonic time-dummy moves by {dhed:+.2f}% (must stay ~0)") pl.DataFrame({ "metric": ["raw_median_shift_pct", "hedonic_shift_pct"], "value": [round(float(dmed), 3), round(float(dhed), 3)], }).write_csv(TABLES_DIR / "composition_shock.csv") # ------------------------------------------------------------------------- # # 4. Weekly vs monthly comparison # ------------------------------------------------------------------------- # def weekly_vs_monthly(hier: pl.DataFrame) -> None: note("") note("## 4. Weekly vs monthly (same methodology)") note("") rows = [] for gid, ptype in [("quebec", "unifamilial"), ("quebec", "condo"), ("montreal", "condo")]: h = hier.filter((pl.col("geography_id") == gid) & (pl.col("property_type") == ptype)).sort("week") wk_raw = np.log(h["index"].to_numpy()) wk_smooth = np.log(h["index_smoothed"].to_numpy()) wk_res = np.log(h["index_research"].to_numpy()) monthly = ( h.with_columns(pl.col("week").str.slice(0, 7).alias("m")) .group_by("m").agg(pl.col("index_research").mean().alias("ix")) .sort("m") ) mo = np.log(monthly["ix"].to_numpy()) # noise = raw minus final smoothed; signal = final smoothed changes noise_sd = float(np.nanstd(wk_raw - wk_res)) signal_sd = float(np.std(np.diff(wk_res))) revision_sd = float(np.std(wk_smooth - wk_res)) # real-time vs final rows.append({ "cell": f"{gid}/{ptype}", "weekly_signal_sd_pct": round(signal_sd * 100, 3), "weekly_noise_sd_pct": round(noise_sd * 100, 3), "monthly_change_sd_pct": round(float(np.std(np.diff(mo))) * 100, 3), "realtime_revision_sd_pct": round(revision_sd * 100, 3), "weeks_to_monthly_ratio": round(signal_sd * np.sqrt(52 / 12) / max(float(np.std(np.diff(mo))), 1e-9), 2), }) note(f"- {gid}/{ptype}: weekly signal SD={rows[-1]['weekly_signal_sd_pct']}%, " f"noise SD={rows[-1]['weekly_noise_sd_pct']}%, monthly ΔSD=" f"{rows[-1]['monthly_change_sd_pct']}%, real-time revision SD=" f"{rows[-1]['realtime_revision_sd_pct']}%") pl.DataFrame(rows).write_csv(TABLES_DIR / "weekly_vs_monthly.csv") def main() -> None: ensure_dirs() cl = research_sample(pl.read_parquet(CLEAN_PARQUET)) feat = build_features(cl) hier = pl.read_parquet("data/interim/hierarchical_indexes.parquet") resid = pl.read_parquet("data/interim/stage1_residuals.parquet") note("# QWHPI validation report") note("") repeat_sales_comparison(cl, hier) downsampling(resid) composition_shock(feat) weekly_vs_monthly(hier) header = ( "\n\n" ) (REPORTS_DIR / "validation_report.md").write_text( header + "\n".join(LINES) + "\n", encoding="utf-8") print("\nReport written to outputs/reports/validation_report.md") if __name__ == "__main__": main()