#!/usr/bin/env python3 # ============================================================================= # QWHPI — Quebec Weekly Housing Price Index # Author : Simon-Pierre Boucher # Contact : contact@spboucher.ai # File : engine/scripts/05_hierarchical.py # Purpose : Step 5 — hierarchical weekly indexes: province, 17 regions and # target cities × type, with Kalman shrinkage for thin cells. # ============================================================================= """Hierarchical index build (Execution Order step 5). Persists stage-1 artifacts (paths, residuals, variances) for reuse by the validation and scaling steps, and writes ``data/interim/hierarchical_indexes.parquet``. """ from __future__ import annotations import sys import time import unicodedata 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 CITY_TARGETS, INTERIM_DIR, PROPERTY_TYPES, TABLES_DIR, ensure_dirs from qwhpi.features import build_features from qwhpi.hierarchy import ( CellDeviation, cell_deviation, compose_cell_table, rebase_log, stage1_residualize, week_grid, ) HIER_PARQUET = INTERIM_DIR / "hierarchical_indexes.parquet" STAGE1_RESID = INTERIM_DIR / "stage1_residuals.parquet" STAGE1_PATHS = INTERIM_DIR / "stage1_paths.parquet" def slug(name: str) -> str: s = unicodedata.normalize("NFKD", name).encode("ascii", "ignore").decode() return s.lower().replace(" ", "-").replace("'", "") def city_slug(name: str) -> str: # Municipality 'Québec' must not collide with the province id. return "quebec-city" if name == "Québec" else slug(name) def aligned_path(paths: pl.DataFrame, ptype: str, grid: list[str]) -> np.ndarray: sub = {r["week"]: r["delta"] for r in paths.filter(pl.col("property_type") == ptype).iter_rows(named=True)} arr = np.array([sub.get(w, np.nan) for w in grid]) # Interpolate the rare missing week levels (singleton-dropped in stage 1). if np.isnan(arr).any(): idx = np.arange(len(arr)) ok = np.isfinite(arr) arr = np.interp(idx, idx[ok], arr[ok]) return arr def combine_all(tables: dict[str, pl.DataFrame], weights: dict[str, float], geography_id: str, geography_name: str, geography_level: str) -> pl.DataFrame: """Fixed-share composition of the three type indexes into 'all'.""" wsum = sum(weights.values()) w = {k: v / wsum for k, v in weights.items()} base = tables[next(iter(tables))].sort("week") grid = base["week"].to_list() def wavg(col: str, log: bool = True) -> np.ndarray: acc = np.zeros(len(grid)) for k, tbl in tables.items(): vals = tbl.sort("week")[col].to_numpy() acc += w[k] * (np.log(vals / 100.0) if log else vals) return acc log_raw = wavg("index") log_filt = wavg("index_smoothed") log_res = wavg("index_research") var = np.zeros(len(grid)) for k, tbl in tables.items(): var += (w[k] ** 2) * (tbl.sort("week")["se_log"].to_numpy() ** 2) se = np.sqrt(var) n = sum(np.nan_to_num(tbl.sort("week")["transactions"].to_numpy()) for tbl in tables.values()) shrink = wavg("shrinkage_weight", log=False) return pl.DataFrame({ "week": grid, "geography_id": geography_id, "geography_name": geography_name, "geography_level": geography_level, "property_type": "all", "index": np.exp(log_raw) * 100.0, "index_smoothed": np.exp(log_filt) * 100.0, "index_research": np.exp(log_res) * 100.0, "log_index": log_filt, "se_log": se, "lower_95": np.exp(log_filt - 1.96 * se) * 100.0, "upper_95": np.exp(log_filt + 1.96 * se) * 100.0, "transactions": n, "shrinkage_weight": shrink, "research_var": var, }) def main() -> None: ensure_dirs() feat = build_features(research_sample(pl.read_parquet(CLEAN_PARQUET))) print("Stage 1: pooled residualization ...") t0 = time.time() s1 = stage1_residualize(feat) print(f" n={s1.n_obs:,} r2={s1.r2:.4f} ({time.time() - t0:.0f}s)") s1.residuals.write_parquet(STAGE1_RESID) s1.paths.write_parquet(STAGE1_PATHS) resid = s1.residuals grid = week_grid(resid) T = len(grid) # Type-level residual variances (fallback for tiny cells) and weekly # counts for the province path variance. sigma2_type = {t: float(resid.filter(pl.col("propertyType") == t)["resid"].var()) for t in PROPERTY_TYPES} frames: list[pl.DataFrame] = [] region_dev_store: dict[tuple[str, str], CellDeviation] = {} region_names = dict( feat.select("region_code", "region").unique().iter_rows() ) for ptype in PROPERTY_TYPES: path = aligned_path(s1.paths, ptype, grid) r_t = resid.filter(pl.col("propertyType") == ptype) counts = {r["week_str"]: r["n"] for r in r_t.group_by("week_str").agg(pl.len().alias("n")).iter_rows(named=True)} n_prov = np.array([counts.get(w, 0) for w in grid], dtype=float) path_var = np.where(n_prov > 0, sigma2_type[ptype] / np.maximum(n_prov, 1), np.nan) path_var = np.nan_to_num(path_var, nan=np.nanmax(path_var)) # ---- province cell ---- # prov = compose_cell_table( grid=grid, province_path=path, deviations=[], geography_id="quebec", geography_name="Québec (province)", geography_level="province", property_type=ptype, path_var=path_var, ).with_columns( pl.Series("transactions", n_prov), pl.Series("shrinkage_weight", np.zeros(T)), ) frames.append(prov) # ---- 17 regions ---- # for rcode in sorted(region_names): cell = r_t.filter(pl.col("region_code") == rcode).rename({"resid": "dev"}) dev = cell_deviation(cell.select("week_str", "dev"), grid, sigma2_type[ptype]) region_dev_store[(rcode, ptype)] = dev frames.append(compose_cell_table( grid=grid, province_path=path, deviations=[dev], geography_id=f"region-{rcode}", geography_name=region_names[rcode], geography_level="region", property_type=ptype, path_var=path_var, )) # ---- target cities ---- # for city in CITY_TARGETS: cell = r_t.filter(pl.col("municipality") == city) if cell.height == 0: continue rcode = cell["region_code"][0] rdev = region_dev_store[(rcode, ptype)] rpath = {w: v for w, v in zip(grid, rdev.fit.smoothed)} cell = cell.with_columns( (pl.col("resid") - pl.col("week_str").replace_strict(rpath, default=0.0)) .alias("dev") ) mdev = cell_deviation(cell.select("week_str", "dev"), grid, sigma2_type[ptype]) frames.append(compose_cell_table( grid=grid, province_path=path, deviations=[rdev, mdev], geography_id=city_slug(city), geography_name=city, geography_level="municipality", property_type=ptype, path_var=path_var, )) out = pl.concat(frames) # ---- 'all' composition per geography (fixed full-sample type shares) ---- # all_frames: list[pl.DataFrame] = [] geos = out.select("geography_id", "geography_name", "geography_level").unique() for gid, gname, glevel in geos.iter_rows(): tables = { t: out.filter((pl.col("geography_id") == gid) & (pl.col("property_type") == t)) for t in PROPERTY_TYPES } tables = {t: tbl for t, tbl in tables.items() if tbl.height} weights = {t: float(np.nansum(tbl["transactions"].to_numpy())) for t, tbl in tables.items()} weights = {t: max(v, 1.0) for t, v in weights.items()} all_frames.append(combine_all(tables, weights, gid, gname, glevel)) out = pl.concat([out, *all_frames]).sort( ["geography_level", "geography_id", "property_type", "week"]) out.write_parquet(HIER_PARQUET) # Summary: volatility comparison raw vs smoothed for thin cells. rows = [] for (gid, ptype), grp in out.to_pandas().groupby(["geography_id", "property_type"]): grp = grp.sort_values("week") li_raw = np.log(grp["index"].to_numpy()) li_smo = np.log(grp["index_smoothed"].to_numpy()) ok = np.isfinite(li_raw) rows.append({ "geography": gid, "property_type": ptype, "median_weekly_tx": float(np.nanmedian(grp["transactions"])), "sd_raw_pct": round(float(np.nanstd(np.diff(li_raw[ok]))) * 100, 2), "sd_smoothed_pct": round(float(np.std(np.diff(li_smo))) * 100, 2), "mean_shrinkage": round(float(np.nanmean(grp["shrinkage_weight"])), 3), "last_index": round(float(grp["index_smoothed"].iloc[-1]), 2), }) summary = pl.DataFrame(rows).sort(["geography", "property_type"]) summary.write_csv(TABLES_DIR / "hierarchical_summary.csv") print(summary) print(f"\nSaved {out.height:,} observations to {HIER_PARQUET}") if __name__ == "__main__": main()