spb/wp22_uqo

1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Benchmark K — re-implementation of the QWHPI v2.1 monthly architecture (PhD/recherche/QWHPI, 2026-08):3province × class = 13-month rolling time dummy, Huber IRLS, mean splice; liquid cells (median ≥ 40 tx/month) = direct4local time dummy lightly smoothed by a local-level filter; thin cells = local-level Kalman deviation from the parent's5published path with per-cell maximum-likelihood τ² (no pooling across cells)."""6from __future__ import annotations78import numpy as np9import polars as pl1011from sharp_hpi.estimators import IndexSeries, pooled_td, rolling_td12from sharp_hpi.hedonic import RHS_STRUCT, assign_loc_fe, fit_fe13from sharp_hpi.statespace import fit_cell1415LIQUID_MIN = 40161718def qhpi_province(df: pl.DataFrame, grid: list[str]) -> IndexSeries:19    d = assign_loc_fe(df, min_n=15)20    return rolling_td(d, grid, window=13, splice="mean", robust=True)212223def qhpi_fullpool_residuals(df: pl.DataFrame, grid: list[str]) -> tuple[pl.DataFrame, float]:24    d = assign_loc_fe(df, min_n=15)25    rhs = RHS_STRUCT.replace(" + C(units_c)", "") if d["property_class"][0] != "plex" else RHS_STRUCT26    if d["property_class"][0] != "house":27        rhs = rhs.replace(" + C(subtype_c)", "")28    res, m = fit_fe(d, rhs + " + C(month)", "loc_fe", robust=True, n_irls=2)29    return d.select(["tx_id", "month", "region_code", "mrc_code", "mun_code"]).with_columns(pl.Series("e", res.resid)), res.sigma2_robust303132def qhpi_cell(df_cell: pl.DataFrame, resid_cell: pl.DataFrame, parent_path: np.ndarray, grid: list[str], sigma2: float) -> tuple[np.ndarray, np.ndarray, str, np.ndarray]:33    """Returns (log path, variance, method, gain) for a sub-provincial cell."""34    cnt = resid_cell.group_by("month").agg(pl.len().alias("n"))35    nmap = dict(zip(cnt["month"].to_list(), cnt["n"].to_list()))36    n = np.array([nmap.get(m, 0) for m in grid], float)37    T = len(grid)38    if np.median(n[n > 0]) >= LIQUID_MIN if (n > 0).any() else False:39        d = assign_loc_fe(df_cell, min_n=15)40        s = pooled_td(d, grid, robust=True)41        y = s.log_index42        ok = ~np.isnan(y)43        # light smoothing: local level with obs var sigma2/n on the direct TD path44        f = fit_cell(np.where(ok, y, np.nan), np.where(ok, n, 0), s.meta.get("sigma_robust", np.sqrt(sigma2)) ** 2, model="ll", robust=False)45        return f.smoothed, f.smoothed_var, "direct_local_td", f.gain46    agg = resid_cell.group_by("month").agg(pl.col("e").mean().alias("yb"))47    ymap = dict(zip(agg["month"].to_list(), agg["yb"].to_list()))48    y = np.array([ymap.get(m, np.nan) for m in grid])49    f = fit_cell(y, n, sigma2, model="ll", robust=False)50    return parent_path + f.smoothed, f.smoothed_var, "hierarchical", f.gain51