1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Arm A: hedonic ridge regression on log price with municipality x quarter and 1-km grid fixed effects.3Arm B: hierarchical hedonic (empirical-Bayes shrinkage of region > municipality > grid residual means)4with a repeat-sales adjustment when the parcel has a prior sale."""5from __future__ import annotations67import numpy as np8import polars as pl9import scipy.sparse as sp10from sklearn.linear_model import Ridge11from sklearn.preprocessing import OneHotEncoder, StandardScaler1213from .common import Prediction, QUANTILES, numeric_features, residual_quantile_intervals141516def _quarter(df: pl.DataFrame) -> pl.Series:17 return (df["date"].dt.year().cast(pl.Utf8) + "Q" + ((df["date"].dt.month() - 1) // 3 + 1).cast(pl.Utf8))181920class HedonicRidge:21 def __init__(self, alpha: float = 1.0, fe: tuple[str, ...] = ("mun_code", "grid1k", "property_class", "subtype", "mun_class")):22 # Out-of-time prediction: calendar-quarter fixed effects (and polynomial trends) cannot be extrapolated23 # safely to 2025-26 (a cubic trend fitted on 2021-24 produced a +0.5 log bias on the 2026 holdout). Time24 # therefore enters only through month-of-year effects and the as-of market-level features of the matrix25 # (12-month municipality and 3-km neighbourhood medians), which are legitimately observable at the sale date.26 self.alpha = alpha27 self.fe = fe28 self.num_cols: list[str] = []29 self.enc = OneHotEncoder(handle_unknown="ignore", min_frequency=5)30 self.scaler = StandardScaler()31 self.model = Ridge(alpha=alpha, solver="sparse_cg", random_state=0)32 self.resid_valid: np.ndarray | None = None3334 def _design(self, df: pl.DataFrame, fit: bool) -> sp.csr_matrix:35 d = df.with_columns(36 (pl.col("mun_code").cast(pl.Utf8) + "_" + pl.col("property_class").cast(pl.Utf8)).alias("mun_class"),37 ((pl.col("date").dt.year() - 2021) * 12 + pl.col("date").dt.month() - 1).cast(pl.Float64).alias("_t"),38 pl.col("date").dt.month().cast(pl.Utf8).alias("_moy"),39 )40 cat = d.select([c for c in self.fe if c in d.columns] + ["_moy"]).fill_null("NA").to_pandas().astype(str)41 if fit:42 self.num_cols = [c for c in numeric_features(df) if c not in ("lat", "lng", "x_m", "y_m", "month_index")]43 X = d.select(self.num_cols).fill_null(strategy="mean").fill_nan(0.0).to_numpy().astype(float)44 X = np.nan_to_num(X)45 self._mun_trend = sp.csr_matrix((len(X), 0))46 if fit:47 Xc = self.enc.fit_transform(cat)48 Xn = self.scaler.fit_transform(X)49 else:50 Xc = self.enc.transform(cat)51 Xn = self.scaler.transform(X)52 return sp.hstack([sp.csr_matrix(Xn), Xc, self._mun_trend]).tocsr()5354 def fit(self, train: pl.DataFrame, valid: pl.DataFrame | None = None) -> "HedonicRidge":55 X = self._design(train, fit=True)56 self.model.fit(X, train["log_price"].to_numpy())57 if valid is not None:58 self.resid_valid = valid["log_price"].to_numpy() - self.model.predict(self._design(valid, fit=False))59 return self6061 def predict_log(self, df: pl.DataFrame) -> np.ndarray:62 return self.model.predict(self._design(df, fit=False))6364 def predict(self, df: pl.DataFrame) -> Prediction:65 p = self.predict_log(df)66 q = residual_quantile_intervals(p, self.resid_valid) if self.resid_valid is not None else None67 return Prediction(tx_id=df["tx_id"].to_numpy(), point=np.exp(p), quantiles=q)686970class HierarchicalHedonic:71 """Ridge hedonic without geographic fixed effects + empirical-Bayes shrinkage of residual means at region,72 municipality and 1-km grid levels (estimated on train), + repeat-sales blend when a prior sale exists.7374 The blend weight w on the repeat-sales forecast (prior log price + municipality index change) is chosen on75 validation by least squares, never on test."""7677 def __init__(self, alpha: float = 1.0, k_shrink: float = 20.0):78 self.base = HedonicRidge(alpha=alpha, fe=("property_class", "subtype"))79 self.k = k_shrink80 self.levels = ("region_name", "mun_code", "grid1k")81 self.effects: dict[str, dict[str, float]] = {}82 self.w_repeat = 0.083 self.resid_valid: np.ndarray | None = None8485 def _shrunk_effects(self, df: pl.DataFrame, resid: np.ndarray) -> np.ndarray:86 r = resid.copy()87 total = np.zeros_like(r)88 for lvl in self.levels:89 key = df[lvl].cast(pl.Utf8).fill_null("NA").to_numpy()90 tab = pl.DataFrame({"k": key, "r": r}).group_by("k").agg(pl.col("r").mean().alias("m"), pl.len().alias("n"))91 eff = {k: m * n / (n + self.k) for k, m, n in zip(tab["k"], tab["m"], tab["n"])}92 self.effects[lvl] = eff93 e = np.array([eff.get(k, 0.0) for k in key])94 r = r - e95 total += e96 return total9798 def _apply_effects(self, df: pl.DataFrame) -> np.ndarray:99 total = np.zeros(df.height)100 for lvl in self.levels:101 key = df[lvl].cast(pl.Utf8).fill_null("NA").to_numpy()102 total += np.array([self.effects[lvl].get(k, 0.0) for k in key])103 return total104105 def _repeat_forecast(self, df: pl.DataFrame) -> np.ndarray:106 if "prior_sale_log_price" in df.columns and "prior_sale_index_change" in df.columns:107 return df.select((pl.col("prior_sale_log_price") + pl.col("prior_sale_index_change").fill_null(0.0))).to_numpy().ravel()108 if "prior_sale_log_price" in df.columns:109 return df["prior_sale_log_price"].to_numpy()110 return np.full(df.height, np.nan)111112 def fit(self, train: pl.DataFrame, valid: pl.DataFrame | None = None) -> "HierarchicalHedonic":113 self.base.fit(train)114 resid = train["log_price"].to_numpy() - self.base.predict_log(train)115 self._shrunk_effects(train, resid)116 if valid is not None:117 h = self.base.predict_log(valid) + self._apply_effects(valid)118 rs = self._repeat_forecast(valid)119 y = valid["log_price"].to_numpy()120 m = ~np.isnan(rs)121 if m.sum() > 50:122 d = rs[m] - h[m]123 self.w_repeat = float(np.clip(np.sum(d * (y[m] - h[m])) / max(np.sum(d * d), 1e-9), 0.0, 1.0))124 self.resid_valid = y - self._combine(h, rs)125 return self126127 def _combine(self, h: np.ndarray, rs: np.ndarray) -> np.ndarray:128 out = h.copy()129 m = ~np.isnan(rs)130 out[m] = (1 - self.w_repeat) * h[m] + self.w_repeat * rs[m]131 return out132133 def predict_log(self, df: pl.DataFrame) -> np.ndarray:134 return self._combine(self.base.predict_log(df) + self._apply_effects(df), self._repeat_forecast(df))135136 def predict(self, df: pl.DataFrame) -> Prediction:137 p = self.predict_log(df)138 q = residual_quantile_intervals(p, self.resid_valid) if self.resid_valid is not None else None139 return Prediction(tx_id=df["tx_id"].to_numpy(), point=np.exp(p), quantiles=q)140