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/09_monthly.py7# Purpose : v2.1 — monthly hybrid indexes: RTD province anchor, DIRECT robust8#           local time-dummy for liquid cells, hierarchical shrinkage for9#           thin cells only.10# =============================================================================11"""Monthly index build (v2.1 headline).1213Architecture (arbitrated against repeat sales, local time dummies and14stratified matched-cell medians — see methodology_v2_research.md):1516- Province × type ....... rolling-time-dummy, Huber IRLS, mean splice.17- Liquid cell (median ≥ LIQUID_MIN tx/month) ... DIRECT local robust18  time-dummy (own β, own FSA FE), then a light state-space smoother19  (obs var σ²/nₘ) for the real-time variant. Deviation-from-pooled-surface20  methods compress strong local divergences (Québec City condo: +50–55%21  by every composition-free arbiter vs +42% by deviations), so liquid22  cells are estimated on their own data.23- Thin cell ............. deviation from its parent's published path via the24  heteroskedastic local-level Kalman model (unchanged machinery).25"""2627from __future__ import annotations2829import sys30import time31import unicodedata32from pathlib import Path3334sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))3536import numpy as np37import polars as pl3839from qwhpi.clean import CLEAN_PARQUET, research_sample40from qwhpi.config import CITY_TARGETS, INTERIM_DIR, PROPERTY_TYPES, TABLES_DIR, ensure_dirs41from qwhpi.features import build_features42from qwhpi.hierarchy import CellDeviation, base_constant, cell_deviation, compose_cell_table43from qwhpi.rtd import fullpool_residualize, local_time_dummy, run_rtd44from qwhpi.state_space import fit_local_level4546MONTHLY_PARQUET = INTERIM_DIR / "monthly_indexes.parquet"47FP_RESID = INTERIM_DIR / "rtd_residuals.parquet"48RTD_PATHS = INTERIM_DIR / "rtd_paths.parquet"4950LIQUID_MIN = 40  # median tx/month for direct local estimation515253def slug(name: str) -> str:54    s = unicodedata.normalize("NFKD", name).encode("ascii", "ignore").decode()55    return s.lower().replace(" ", "-").replace("'", "")565758def city_slug(name: str) -> str:59    return "quebec-city" if name == "Québec" else slug(name)606162def compose_direct(grid, raw_path, n, sigma2, gid, gname, glevel, ptype):63    """Published table for a directly-estimated liquid cell."""64    fit = fit_local_level(raw_path, np.maximum(n, 1.0), sigma2)65    base = base_constant(grid, fit.filtered)66    log_raw = raw_path - base67    log_filt = fit.filtered - base68    log_smooth = fit.smoothed - base69    se = np.sqrt(np.maximum(fit.filtered_var, 0.0))70    return pl.DataFrame({71        "week": grid,72        "geography_id": gid, "geography_name": gname, "geography_level": glevel,73        "property_type": ptype,74        "index": np.exp(log_raw) * 100.0,75        "index_smoothed": np.exp(log_filt) * 100.0,76        "index_research": np.exp(log_smooth) * 100.0,77        "log_index": log_filt,78        "se_log": se,79        "lower_95": np.exp(log_filt - 1.96 * se) * 100.0,80        "upper_95": np.exp(log_filt + 1.96 * se) * 100.0,81        "transactions": n,82        "shrinkage_weight": 1.0 - fit.gain,  # temporal smoothing weight83        "research_var": np.maximum(fit.smoothed_var, 0.0),84    }), fit858687def combine_all(tables, weights, gid, gname, glevel):88    """Fixed-share composition of the three type indexes into 'all'."""89    wsum = sum(weights.values())90    w = {k: v / wsum for k, v in weights.items()}91    base = tables[next(iter(tables))].sort("week")92    grid = base["week"].to_list()9394    def wavg(col, log=True):95        acc = np.zeros(len(grid))96        for k, tbl in tables.items():97            vals = tbl.sort("week")[col].to_numpy()98            acc += w[k] * (np.log(vals / 100.0) if log else vals)99        return acc100101    log_raw, log_filt, log_res = wavg("index"), wavg("index_smoothed"), wavg("index_research")102    var = sum((w[k] ** 2) * (tbl.sort("week")["se_log"].to_numpy() ** 2)103              for k, tbl in tables.items())104    se = np.sqrt(var)105    n = sum(np.nan_to_num(tbl.sort("week")["transactions"].to_numpy())106            for tbl in tables.values())107    return pl.DataFrame({108        "week": grid,109        "geography_id": gid, "geography_name": gname, "geography_level": glevel,110        "property_type": "all",111        "index": np.exp(log_raw) * 100.0,112        "index_smoothed": np.exp(log_filt) * 100.0,113        "index_research": np.exp(log_res) * 100.0,114        "log_index": log_filt,115        "se_log": se,116        "lower_95": np.exp(log_filt - 1.96 * se) * 100.0,117        "upper_95": np.exp(log_filt + 1.96 * se) * 100.0,118        "transactions": n,119        "shrinkage_weight": wavg("shrinkage_weight", log=False),120        "research_var": var,121    })122123124def main() -> None:125    ensure_dirs()126    feat = build_features(research_sample(pl.read_parquet(CLEAN_PARQUET)))127128    print("Stage 1: province RTD (13-month windows, Huber IRLS) ...")129    t0 = time.time()130    rtd = run_rtd(feat, verbose=False)131    print(f"  done in {time.time() - t0:.0f}s")132    rtd.drift.write_csv(TABLES_DIR / "rtd_coefficient_drift.csv")133134    print("Stage 1b: full-pool robust residualization (thin-cell hierarchy) ...")135    fp_resid = fullpool_residualize(feat)136    fp_resid.write_parquet(FP_RESID)137    sigma2 = {t: float(fp_resid.filter(pl.col("propertyType") == t)["resid"].var())138              for t in PROPERTY_TYPES}139140    grid = rtd.months141    T = len(grid)142    pl.concat([143        pl.DataFrame({"month": grid, "property_type": t,144                      "mean_splice": rtd.paths[t],145                      "movement_splice": rtd.paths_movement[t],146                      "window_splice": rtd.paths_window[t]})147        for t in PROPERTY_TYPES148    ]).write_parquet(RTD_PATHS)149150    resid = fp_resid.rename({"month": "week_str"})151    region_names = dict(feat.select("region_code", "region").unique().iter_rows())152    frames: list[pl.DataFrame] = []153    estimation_log: list[dict] = []154155    # Parent references for thin cells:156    region_resid_dev: dict[tuple[str, str], CellDeviation] = {}157    region_anchor: dict[tuple[str, str], np.ndarray] = {}  # published filtered log path158159    for ptype in PROPERTY_TYPES:160        path = rtd.paths[ptype]161        r_t = resid.filter(pl.col("propertyType") == ptype)162        f_t = feat.filter(pl.col("propertyType") == ptype)163        counts = {r["week_str"]: r["n"] for r in164                  r_t.group_by("week_str").agg(pl.len().alias("n")).iter_rows(named=True)}165        n_prov = np.array([counts.get(m, 0) for m in grid], dtype=float)166        path_var = np.where(n_prov > 0, sigma2[ptype] / np.maximum(n_prov, 1), np.nan)167        path_var = np.nan_to_num(path_var, nan=float(np.nanmax(path_var)))168169        prov = compose_cell_table(170            grid=grid, province_path=path, deviations=[],171            geography_id="quebec", geography_name="Québec (province)",172            geography_level="province", property_type=ptype, path_var=path_var,173        ).with_columns(pl.Series("transactions", n_prov),174                       pl.Series("shrinkage_weight", np.zeros(T)))175        frames.append(prov)176        estimation_log.append({"cell": f"quebec/{ptype}", "method": "rtd_mean_splice"})177178        def monthly_counts(sub: pl.DataFrame) -> np.ndarray:179            c = dict(sub.group_by("month").agg(pl.len().alias("n")).iter_rows())180            return np.array([c.get(m, 0) for m in grid], dtype=float)181182        # ---------------- regions ---------------- #183        for rcode in sorted(region_names):184            cell_r = r_t.filter(pl.col("region_code") == rcode)185            cell_f = f_t.filter(pl.col("region_code") == rcode)186            n_m = monthly_counts(cell_f)187            med = float(np.median(n_m))188            rdev = cell_deviation(cell_r.rename({"resid": "dev"})189                                  .select("week_str", "dev"), grid, sigma2[ptype])190            region_resid_dev[(rcode, ptype)] = rdev191            gid = f"region-{rcode}"192193            direct = local_time_dummy(cell_f, grid) if med >= LIQUID_MIN else None194            if direct is not None:195                raw_path, n_cell, s2 = direct196                tbl, fit = compose_direct(grid, raw_path, n_cell, s2,197                                          gid, region_names[rcode], "region", ptype)198                frames.append(tbl)199                region_anchor[(rcode, ptype)] = raw_path * 0 + fit.filtered  # filtered log path200                estimation_log.append({"cell": f"{gid}/{ptype}", "method": "direct_local_td",201                                       "median_tx": med})202            else:203                tbl = compose_cell_table(204                    grid=grid, province_path=path, deviations=[rdev],205                    geography_id=gid, geography_name=region_names[rcode],206                    geography_level="region", property_type=ptype, path_var=path_var)207                frames.append(tbl)208                region_anchor[(rcode, ptype)] = path + rdev.fit.filtered209                estimation_log.append({"cell": f"{gid}/{ptype}", "method": "hierarchical",210                                       "median_tx": med})211212        # ---------------- cities ---------------- #213        for city in CITY_TARGETS:214            cell_f = f_t.filter(pl.col("municipality") == city)215            if cell_f.height == 0:216                continue217            rcode = cell_f["region_code"][0]218            gid = city_slug(city)219            n_m = monthly_counts(cell_f)220            med = float(np.median(n_m))221222            direct = local_time_dummy(cell_f, grid) if med >= LIQUID_MIN else None223            if direct is not None:224                raw_path, n_cell, s2 = direct225                tbl, _ = compose_direct(grid, raw_path, n_cell, s2,226                                        gid, city, "municipality", ptype)227                frames.append(tbl)228                estimation_log.append({"cell": f"{gid}/{ptype}", "method": "direct_local_td",229                                       "median_tx": med})230            else:231                # deviation of the city vs its REGION's residual path, anchored232                # on the region's published (filtered) log path233                rdev = region_resid_dev[(rcode, ptype)]234                rpath = {m: v for m, v in zip(grid, rdev.fit.smoothed)}235                cell_r = r_t.filter(pl.col("municipality") == city).with_columns(236                    (pl.col("resid") - pl.col("week_str").replace_strict(rpath, default=0.0))237                    .alias("dev"))238                mdev = cell_deviation(cell_r.select("week_str", "dev"), grid, sigma2[ptype])239                tbl = compose_cell_table(240                    grid=grid, province_path=region_anchor[(rcode, ptype)],241                    deviations=[mdev],242                    geography_id=gid, geography_name=city,243                    geography_level="municipality", property_type=ptype,244                    path_var=path_var)245                frames.append(tbl)246                estimation_log.append({"cell": f"{gid}/{ptype}", "method": "hierarchical",247                                       "median_tx": med})248249    out = pl.concat(frames)250251    all_frames = []252    for gid, gname, glevel in out.select(253            "geography_id", "geography_name", "geography_level").unique().iter_rows():254        tables = {t: out.filter((pl.col("geography_id") == gid)255                                & (pl.col("property_type") == t))256                  for t in PROPERTY_TYPES}257        tables = {t: tbl for t, tbl in tables.items() if tbl.height}258        weights = {t: max(float(np.nansum(tbl["transactions"].to_numpy())), 1.0)259                   for t, tbl in tables.items()}260        all_frames.append(combine_all(tables, weights, gid, gname, glevel))261262    out = pl.concat([out, *all_frames]).rename({"week": "period"}).sort(263        ["geography_level", "geography_id", "property_type", "period"])264    out.write_parquet(MONTHLY_PARQUET)265266    log_df = pl.DataFrame(estimation_log)267    log_df.write_csv(TABLES_DIR / "estimation_methods.csv")268    n_direct = log_df.filter(pl.col("method") == "direct_local_td").height269    n_hier = log_df.filter(pl.col("method") == "hierarchical").height270    print(f"\nCells: {n_direct} direct, {n_hier} hierarchical, 3 RTD province")271    print(f"Saved {out.height:,} observations to {MONTHLY_PARQUET}")272273274if __name__ == "__main__":275    main()276