TeX 55.1%
Python 44.2%
1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""SHARP-HPI: Spatially Hierarchical, Adaptive and Robust housing price index.34Decomposition (per class k):5 log P_i = f(X_i; θ_t) + α_loc(i) + μ_{k,t} + δ^R_{r(i),t} + δ^M_{m(i),t} + δ^U_{u(i),t} + ε_i67Stage 1 — robust pooled hedonic (Huber-IRLS) with fine-location fixed effects α_loc (nested metric grid),8 province × month effects μ_t, and smoothly time-varying quality prices (interactions of the9 structural terms with a centred linear time index). Residuals e_i are quality- and10 location-adjusted log prices net of the provincial path.11Stage 2 — nested deviations. For each level ℓ ∈ {region, MRC, municipality} the monthly mean residual of a12 cell (net of its parent's smoothed/filtered deviation) is modelled by a heteroskedastic local-level13 (or local-linear-trend) state-space model with observation variance σ²_cell / n_t. Innovation14 variances τ² are calibrated by empirical Bayes across the cells of a level (grid EM on a normal prior for15 log τ²), so a thin cell's responsiveness is borrowed from its peers instead of being estimated from16 its own handful of sales (the adaptive pooling). Innovations are Huber down-weighted (robust).17Real-time mode uses the filtered path at every level; retrospective mode the smoothed path.1819Reliability bundle per cell-month: actual count, effective sample size (σ²_cell / Var(state)), local data weight20(Kalman gain), posterior sd, covariate overlap, residual autocorrelation, one-step-ahead standardized error.21"""22from __future__ import annotations2324from dataclasses import dataclass, field2526import numpy as np27import polars as pl2829from sharp_hpi.hedonic import RHS_STRUCT, assign_loc_fe, fit_fe30from sharp_hpi.statespace import LOG_TAU2_GRID, SSFit, eb_prior, fit_cell, profile_loglik3132LEVELS = [("region", "region_code"), ("mrc", "mrc_code"), ("municipality", "mun_code")]33TV_TERMS = ["lfa", "lfa2", "lla"] # quality prices allowed to drift linearly in time (ablation: none)343536@dataclass37class Stage1:38 months: list[str]39 mu: np.ndarray # provincial log path (month effects, first month = 0)40 mu_se: np.ndarray41 resid: pl.DataFrame # tx_id, month, region_code, mrc_code, mun_code, e42 sigma2: float43 coef: dict44 n: int45 r2_within: float | None464748@dataclass49class CellFit:50 level: str51 geo_id: str52 parent_id: str | None53 months: list[str]54 n: np.ndarray55 ybar: np.ndarray # observed monthly mean deviation (NaN if none)56 fit: SSFit57 sigma2_cell: float58 overlap: np.ndarray # share of the cell's sales inside the class common support59 rho1: float # lag-1 autocorrelation of observed deviations (residual dependence)60 oos_z: np.ndarray # one-step-ahead standardized innovation (NaN if none)616263@dataclass64class SharpResult:65 cls: str66 stage1: Stage167 cells: dict[tuple[str, str], CellFit] = field(default_factory=dict)68 priors: dict[str, tuple[float, float]] = field(default_factory=dict)69 config: dict = field(default_factory=dict)707172def _rhs(df: pl.DataFrame, time_varying: bool) -> str:73 rhs = RHS_STRUCT74 cls = df["property_class"][0]75 if cls != "plex":76 rhs = rhs.replace(" + C(units_c)", "")77 if cls != "house":78 rhs = rhs.replace(" + C(subtype_c)", "")79 if time_varying:80 rhs += " + " + " + ".join(f"{t}:tt" for t in TV_TERMS)81 return rhs + " + C(month)"828384def stage1(df: pl.DataFrame, grid: list[str], time_varying: bool = True, min_loc_n: int = 15, robust: bool = True, n_irls: int = 2) -> Stage1:85 T = len(grid)86 tt = {m: (i - (T - 1) / 2) / max(T - 1, 1) for i, m in enumerate(grid)}87 d = assign_loc_fe(df, min_n=min_loc_n).with_columns(pl.col("month").replace_strict(tt, default=0.0).alias("tt"))88 res, m = fit_fe(d, _rhs(d, time_varying), "loc_fe", robust=robust, n_irls=n_irls)89 se = dict(m.se())90 base = grid[0]91 mu = np.array([0.0 if g == base else res.coef.get(f"C(month)[T.{g}]", np.nan) for g in grid])92 mu_se = np.array([0.0 if g == base else se.get(f"C(month)[T.{g}]", np.nan) for g in grid])93 resid = d.select(["tx_id", "month", "region_code", "mrc_code", "mun_code", "lfa", "age_f"]).with_columns(pl.Series("e", res.resid))94 return Stage1(grid, mu, mu_se, resid, res.sigma2_robust, res.coef, res.n, res.r2_within)959697def _cell_series(df: pl.DataFrame, key: str, grid: list[str], ycol: str = "e") -> dict[str, tuple[np.ndarray, np.ndarray, float, np.ndarray]]:98 """Per cell id: (ybar by month, n by month, robust sigma2, overlap by month)."""99 agg = df.group_by([key, "month"]).agg(pl.col(ycol).mean().alias("ybar"), pl.len().alias("n"), pl.col("_in_support").mean().alias("ov"))100 sig = df.group_by(key).agg(((pl.col(ycol) - pl.col(ycol).median()).abs().median() * 1.4826).pow(2).alias("s2"), pl.len().alias("N"))101 idx = {m: i for i, m in enumerate(grid)}102 out = {}103 s2 = dict(zip(sig[key].to_list(), zip(sig["s2"].to_list(), sig["N"].to_list())))104 for gid, sub in agg.group_by(key):105 gid = gid[0]106 y = np.full(len(grid), np.nan); n = np.zeros(len(grid)); ov = np.full(len(grid), np.nan)107 for mth, yb, nn, o in zip(sub["month"].to_list(), sub["ybar"].to_list(), sub["n"].to_list(), sub["ov"].to_list()):108 if mth in idx:109 y[idx[mth]] = yb; n[idx[mth]] = nn; ov[idx[mth]] = o110 out[gid] = (y, n, s2[gid], ov)111 return out112113114def _acf1(y: np.ndarray) -> float:115 v = y[~np.isnan(y)]116 if len(v) < 6:117 return np.nan118 v = v - v.mean()119 return float(np.sum(v[1:] * v[:-1]) / np.sum(v * v)) if np.sum(v * v) > 0 else np.nan120121122def stage2(s1: Stage1, df_geo: pl.DataFrame, grid: list[str], model: str = "ll", eb: bool = True, robust: bool = True,123 realtime_parent: bool = False, sigma2_shrink_n: int = 100, levels=LEVELS) -> tuple[dict, dict]:124 """Nested state-space deviations. Returns (cells, priors)."""125 e = s1.resid126 # common support: class-level 1–99 % band of (centred) log floor area and age127 q = e.select(pl.col("lfa").quantile(0.01).alias("l1"), pl.col("lfa").quantile(0.99).alias("l99"), pl.col("age_f").quantile(0.99).alias("a99"))128 e = e.with_columns((pl.col("lfa").is_between(q["l1"][0], q["l99"][0]) & (pl.col("age_f") <= q["a99"][0])).alias("_in_support"))129 e = e.with_columns(pl.lit(0.0).alias("parent_dev"))130 cells: dict[tuple[str, str], CellFit] = {}131 priors: dict[str, tuple[float, float]] = {}132 parent_key = None133 for lname, key in levels:134 e = e.with_columns((pl.col("e") - pl.col("parent_dev")).alias("e_net"))135 series = _cell_series(e, key, grid, "e_net")136 # cell sigma2 shrunk toward the class robust variance for small cells137 for gid, (y, n, (s2, N), ov) in series.items():138 w = N / (N + sigma2_shrink_n)139 series[gid] = (y, n, (w * s2 + (1 - w) * s1.sigma2, N), ov)140 prior = None141 if eb:142 profs = [profile_loglik(y, n, s2[0], model) for (y, n, s2, ov) in series.values() if (n > 0).sum() >= 6]143 prior = eb_prior(profs) if len(profs) >= 3 else None144 priors[lname] = prior if prior else (float(np.log(1e-5)), 1.0)145 dev_map = {}146 for gid, (y, n, (s2, N), ov) in series.items():147 fit = fit_cell(y, n, s2, model=model, prior=prior if eb else None, robust=robust)148 parent = e.filter(pl.col(key) == gid)[parent_key][0] if parent_key else None149 z = np.where(n > 0, (np.nan_to_num(y) - np.concatenate([[0.0], fit.filtered[:-1]])) / np.sqrt(np.where(n > 0, s2 / np.maximum(n, 1), 1) + fit.tau2 + np.concatenate([[0.25], fit.filtered_var[:-1]])), np.nan)150 cells[(lname, gid)] = CellFit(lname, gid, parent, grid, n, y, fit, s2, ov, _acf1(y), z)151 dev_map[gid] = fit.filtered if realtime_parent else fit.smoothed152 # pass this level's deviation to the children as their parent path153 idx = {m: i for i, m in enumerate(grid)}154 dev_tbl = pl.DataFrame({key: list(dev_map.keys()), "_dev": [list(v) for v in dev_map.values()]}).explode("_dev").with_columns(155 pl.Series("month", grid * len(dev_map)))156 e = e.drop("parent_dev").join(dev_tbl, on=[key, "month"], how="left").with_columns(pl.col("_dev").fill_null(0.0).alias("_dev"))157 e = e.with_columns((pl.col("parent_dev_prev") if "parent_dev_prev" in e.columns else pl.lit(0.0)).alias("_pp"))158 e = e.with_columns((pl.col("_dev") + pl.col("_pp")).alias("parent_dev")).drop("_dev", "_pp")159 e = e.with_columns(pl.col("parent_dev").alias("parent_dev_prev"))160 parent_key = key161 return cells, priors162163164def run_sharp(df: pl.DataFrame, grid: list[str], time_varying: bool = True, model: str = "ll", eb: bool = True, robust: bool = True,165 realtime_parent: bool = False, min_loc_n: int = 15, levels=LEVELS) -> SharpResult:166 cls = df["property_class"][0]167 s1 = stage1(df, grid, time_varying=time_varying, min_loc_n=min_loc_n, robust=robust)168 cells, priors = stage2(s1, df, grid, model=model, eb=eb, robust=robust, realtime_parent=realtime_parent, levels=levels)169 return SharpResult(cls, s1, cells, priors, {"time_varying": time_varying, "model": model, "eb": eb, "robust": robust, "realtime_parent": realtime_parent, "min_loc_n": min_loc_n})170171172def cell_log_index(res: SharpResult, level: str, geo_id: str, mode: str = "smoothed") -> tuple[np.ndarray, np.ndarray]:173 """Log index of a cell = provincial path + sum of its ancestors' and own deviations; variance = sum of state variances174 (independence across levels is an approximation; stage-1 μ uncertainty is added)."""175 path = res.stage1.mu.copy(); var = res.stage1.mu_se ** 2176 chain = []177 lvl_names = [l for l, _ in LEVELS]178 cur = (level, geo_id)179 while cur is not None and cur in res.cells:180 c = res.cells[cur]; chain.append(c)181 li = lvl_names.index(c.level)182 cur = (lvl_names[li - 1], c.parent_id) if li > 0 and c.parent_id is not None else None183 for c in chain:184 f = c.fit185 path = path + (f.smoothed if mode == "smoothed" else f.filtered)186 var = var + (f.smoothed_var if mode == "smoothed" else f.filtered_var)187 return path, var188