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/05_hierarchical.py7# Purpose : Step 5 — hierarchical weekly indexes: province, 17 regions and8#           target cities × type, with Kalman shrinkage for thin cells.9# =============================================================================10"""Hierarchical index build (Execution Order step 5).1112Persists stage-1 artifacts (paths, residuals, variances) for reuse by the13validation and scaling steps, and writes14``data/interim/hierarchical_indexes.parquet``.15"""1617from __future__ import annotations1819import sys20import time21import unicodedata22from pathlib import Path2324sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2526import numpy as np27import polars as pl2829from qwhpi.clean import CLEAN_PARQUET, research_sample30from qwhpi.config import CITY_TARGETS, INTERIM_DIR, PROPERTY_TYPES, TABLES_DIR, ensure_dirs31from qwhpi.features import build_features32from qwhpi.hierarchy import (33    CellDeviation, cell_deviation, compose_cell_table, rebase_log,34    stage1_residualize, week_grid,35)3637HIER_PARQUET = INTERIM_DIR / "hierarchical_indexes.parquet"38STAGE1_RESID = INTERIM_DIR / "stage1_residuals.parquet"39STAGE1_PATHS = INTERIM_DIR / "stage1_paths.parquet"404142def slug(name: str) -> str:43    s = unicodedata.normalize("NFKD", name).encode("ascii", "ignore").decode()44    return s.lower().replace(" ", "-").replace("'", "")454647def city_slug(name: str) -> str:48    # Municipality 'Québec' must not collide with the province id.49    return "quebec-city" if name == "Québec" else slug(name)505152def aligned_path(paths: pl.DataFrame, ptype: str, grid: list[str]) -> np.ndarray:53    sub = {r["week"]: r["delta"] for r in54           paths.filter(pl.col("property_type") == ptype).iter_rows(named=True)}55    arr = np.array([sub.get(w, np.nan) for w in grid])56    # Interpolate the rare missing week levels (singleton-dropped in stage 1).57    if np.isnan(arr).any():58        idx = np.arange(len(arr))59        ok = np.isfinite(arr)60        arr = np.interp(idx, idx[ok], arr[ok])61    return arr626364def combine_all(tables: dict[str, pl.DataFrame], weights: dict[str, float],65                geography_id: str, geography_name: str,66                geography_level: str) -> pl.DataFrame:67    """Fixed-share composition of the three type indexes into 'all'."""68    wsum = sum(weights.values())69    w = {k: v / wsum for k, v in weights.items()}70    base = tables[next(iter(tables))].sort("week")71    grid = base["week"].to_list()7273    def wavg(col: str, log: bool = True) -> np.ndarray:74        acc = np.zeros(len(grid))75        for k, tbl in tables.items():76            vals = tbl.sort("week")[col].to_numpy()77            acc += w[k] * (np.log(vals / 100.0) if log else vals)78        return acc7980    log_raw = wavg("index")81    log_filt = wavg("index_smoothed")82    log_res = wavg("index_research")83    var = np.zeros(len(grid))84    for k, tbl in tables.items():85        var += (w[k] ** 2) * (tbl.sort("week")["se_log"].to_numpy() ** 2)86    se = np.sqrt(var)87    n = sum(np.nan_to_num(tbl.sort("week")["transactions"].to_numpy())88            for tbl in tables.values())89    shrink = wavg("shrinkage_weight", log=False)9091    return pl.DataFrame({92        "week": grid,93        "geography_id": geography_id,94        "geography_name": geography_name,95        "geography_level": geography_level,96        "property_type": "all",97        "index": np.exp(log_raw) * 100.0,98        "index_smoothed": np.exp(log_filt) * 100.0,99        "index_research": np.exp(log_res) * 100.0,100        "log_index": log_filt,101        "se_log": se,102        "lower_95": np.exp(log_filt - 1.96 * se) * 100.0,103        "upper_95": np.exp(log_filt + 1.96 * se) * 100.0,104        "transactions": n,105        "shrinkage_weight": shrink,106        "research_var": var,107    })108109110def main() -> None:111    ensure_dirs()112    feat = build_features(research_sample(pl.read_parquet(CLEAN_PARQUET)))113114    print("Stage 1: pooled residualization ...")115    t0 = time.time()116    s1 = stage1_residualize(feat)117    print(f"  n={s1.n_obs:,} r2={s1.r2:.4f} ({time.time() - t0:.0f}s)")118    s1.residuals.write_parquet(STAGE1_RESID)119    s1.paths.write_parquet(STAGE1_PATHS)120121    resid = s1.residuals122    grid = week_grid(resid)123    T = len(grid)124125    # Type-level residual variances (fallback for tiny cells) and weekly126    # counts for the province path variance.127    sigma2_type = {t: float(resid.filter(pl.col("propertyType") == t)["resid"].var())128                   for t in PROPERTY_TYPES}129130    frames: list[pl.DataFrame] = []131    region_dev_store: dict[tuple[str, str], CellDeviation] = {}132133    region_names = dict(134        feat.select("region_code", "region").unique().iter_rows()135    )136137    for ptype in PROPERTY_TYPES:138        path = aligned_path(s1.paths, ptype, grid)139        r_t = resid.filter(pl.col("propertyType") == ptype)140        counts = {r["week_str"]: r["n"] for r in141                  r_t.group_by("week_str").agg(pl.len().alias("n")).iter_rows(named=True)}142        n_prov = np.array([counts.get(w, 0) for w in grid], dtype=float)143        path_var = np.where(n_prov > 0, sigma2_type[ptype] / np.maximum(n_prov, 1), np.nan)144        path_var = np.nan_to_num(path_var, nan=np.nanmax(path_var))145146        # ---- province cell ---- #147        prov = compose_cell_table(148            grid=grid, province_path=path, deviations=[],149            geography_id="quebec", geography_name="Québec (province)",150            geography_level="province", property_type=ptype, path_var=path_var,151        ).with_columns(152            pl.Series("transactions", n_prov),153            pl.Series("shrinkage_weight", np.zeros(T)),154        )155        frames.append(prov)156157        # ---- 17 regions ---- #158        for rcode in sorted(region_names):159            cell = r_t.filter(pl.col("region_code") == rcode).rename({"resid": "dev"})160            dev = cell_deviation(cell.select("week_str", "dev"), grid, sigma2_type[ptype])161            region_dev_store[(rcode, ptype)] = dev162            frames.append(compose_cell_table(163                grid=grid, province_path=path, deviations=[dev],164                geography_id=f"region-{rcode}",165                geography_name=region_names[rcode],166                geography_level="region", property_type=ptype, path_var=path_var,167            ))168169        # ---- target cities ---- #170        for city in CITY_TARGETS:171            cell = r_t.filter(pl.col("municipality") == city)172            if cell.height == 0:173                continue174            rcode = cell["region_code"][0]175            rdev = region_dev_store[(rcode, ptype)]176            rpath = {w: v for w, v in zip(grid, rdev.fit.smoothed)}177            cell = cell.with_columns(178                (pl.col("resid") - pl.col("week_str").replace_strict(rpath, default=0.0))179                .alias("dev")180            )181            mdev = cell_deviation(cell.select("week_str", "dev"), grid, sigma2_type[ptype])182            frames.append(compose_cell_table(183                grid=grid, province_path=path, deviations=[rdev, mdev],184                geography_id=city_slug(city), geography_name=city,185                geography_level="municipality", property_type=ptype, path_var=path_var,186            ))187188    out = pl.concat(frames)189190    # ---- 'all' composition per geography (fixed full-sample type shares) ---- #191    all_frames: list[pl.DataFrame] = []192    geos = out.select("geography_id", "geography_name", "geography_level").unique()193    for gid, gname, glevel in geos.iter_rows():194        tables = {195            t: out.filter((pl.col("geography_id") == gid) & (pl.col("property_type") == t))196            for t in PROPERTY_TYPES197        }198        tables = {t: tbl for t, tbl in tables.items() if tbl.height}199        weights = {t: float(np.nansum(tbl["transactions"].to_numpy()))200                   for t, tbl in tables.items()}201        weights = {t: max(v, 1.0) for t, v in weights.items()}202        all_frames.append(combine_all(tables, weights, gid, gname, glevel))203204    out = pl.concat([out, *all_frames]).sort(205        ["geography_level", "geography_id", "property_type", "week"])206    out.write_parquet(HIER_PARQUET)207208    # Summary: volatility comparison raw vs smoothed for thin cells.209    rows = []210    for (gid, ptype), grp in out.to_pandas().groupby(["geography_id", "property_type"]):211        grp = grp.sort_values("week")212        li_raw = np.log(grp["index"].to_numpy())213        li_smo = np.log(grp["index_smoothed"].to_numpy())214        ok = np.isfinite(li_raw)215        rows.append({216            "geography": gid, "property_type": ptype,217            "median_weekly_tx": float(np.nanmedian(grp["transactions"])),218            "sd_raw_pct": round(float(np.nanstd(np.diff(li_raw[ok]))) * 100, 2),219            "sd_smoothed_pct": round(float(np.std(np.diff(li_smo))) * 100, 2),220            "mean_shrinkage": round(float(np.nanmean(grp["shrinkage_weight"])), 3),221            "last_index": round(float(grp["index_smoothed"].iloc[-1]), 2),222        })223    summary = pl.DataFrame(rows).sort(["geography", "property_type"])224    summary.write_csv(TABLES_DIR / "hierarchical_summary.csv")225    print(summary)226    print(f"\nSaved {out.height:,} observations to {HIER_PARQUET}")227228229if __name__ == "__main__":230    main()231