TeX 55.1%
Python 44.2%
1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Benchmark index estimators (protocol §5, A–J) on one geography × class sample.34Every estimator returns an `IndexSeries` of **log** index levels on a continuous month grid5(NaN where not estimable). Normalisation (2021 average = 100 or January 2021 = 100) is applied6downstream in `panel.py`, so it never changes growth rates.7"""8from __future__ import annotations910from dataclasses import dataclass, field1112import numpy as np13import polars as pl14import pyfixest as pf15from scipy import sparse16from scipy.sparse.linalg import lsqr1718from sharp_hpi.hedonic import RHS_STRUCT, fit_fe, huber_weights192021@dataclass22class IndexSeries:23 months: list[str]24 log_index: np.ndarray25 se: np.ndarray | None = None26 n: np.ndarray | None = None27 method: str = ""28 meta: dict = field(default_factory=dict)2930 def to_frame(self) -> pl.DataFrame:31 d = {"month": self.months, "log_index": self.log_index, "method": self.method}32 if self.se is not None:33 d["se"] = self.se34 if self.n is not None:35 d["n"] = self.n36 return pl.DataFrame(d)373839def month_grid(first: str, last: str) -> list[str]:40 y, m = int(first[:4]), int(first[5:7])41 out = []42 while True:43 cur = f"{y:04d}-{m:02d}"44 out.append(cur)45 if cur == last:46 return out47 m += 148 if m == 13:49 y, m = y + 1, 1505152def _counts(df: pl.DataFrame, grid: list[str]) -> np.ndarray:53 c = df.group_by("month").agg(pl.len().alias("n"))54 d = dict(zip(c["month"].to_list(), c["n"].to_list()))55 return np.array([d.get(m, 0) for m in grid], float)565758def _rhs_for(df: pl.DataFrame) -> str:59 rhs = RHS_STRUCT60 if df["property_class"].n_unique() == 1 and df["property_class"][0] != "plex":61 rhs = rhs.replace(" + C(units_c)", "")62 if df["property_class"].n_unique() == 1 and df["property_class"][0] != "house":63 rhs = rhs.replace(" + C(subtype_c)", "")64 return rhs656667# A. raw median / mean -------------------------------------------------------68def median_index(df: pl.DataFrame, grid: list[str], stat: str = "median") -> IndexSeries:69 agg = pl.col("price").median() if stat == "median" else pl.col("price").mean()70 g = df.group_by("month").agg(agg.alias("v")).sort("month")71 d = dict(zip(g["month"].to_list(), g["v"].to_list()))72 li = np.array([np.log(d[m]) if m in d else np.nan for m in grid])73 return IndexSeries(grid, li, n=_counts(df, grid), method=f"A_{stat}")747576# B. stratified median (fixed-weight Laspeyres over strata) ------------------77def stratified_median(df: pl.DataFrame, grid: list[str], n_strata_fa: int = 3, n_strata_age: int = 3) -> IndexSeries:78 qfa = [df["fa_f"].quantile(q) for q in np.linspace(0, 1, n_strata_fa + 1)[1:-1]]79 qage = [df["age_f"].quantile(q) for q in np.linspace(0, 1, n_strata_age + 1)[1:-1]]80 d = df.with_columns(81 pl.col("fa_f").cut(qfa, labels=[f"fa{i}" for i in range(n_strata_fa)]).cast(pl.Utf8).alias("_sfa"),82 pl.col("age_f").cut(qage, labels=[f"ag{i}" for i in range(n_strata_age)]).cast(pl.Utf8).alias("_sag"),83 ).with_columns(pl.concat_str([pl.col("_sfa"), pl.col("_sag"), pl.col("subtype_c")]).alias("stratum"))84 w = d.group_by("stratum").agg(pl.len().alias("w")).with_columns(pl.col("w") / pl.col("w").sum())85 med = d.group_by(["stratum", "month"]).agg(pl.col("price").median().log().alias("lm"))86 strata = w["stratum"].to_list()87 wv = dict(zip(strata, w["w"].to_list()))88 lookup = {(s, m): v for s, m, v in zip(med["stratum"].to_list(), med["month"].to_list(), med["lm"].to_list())}89 out = np.full(len(grid), np.nan)90 last = {s: np.nan for s in strata}91 for i, m in enumerate(grid):92 acc, wsum = 0.0, 0.093 for s in strata:94 v = lookup.get((s, m), np.nan)95 if np.isnan(v):96 v = last[s] # carry forward an empty stratum's last median97 else:98 last[s] = v99 if not np.isnan(v):100 acc += wv[s] * v101 wsum += wv[s]102 out[i] = acc / wsum if wsum > 0.5 else np.nan103 return IndexSeries(grid, out, n=_counts(df, grid), method="B_stratified_median", meta={"n_strata": len(strata)})104105106# C. pooled hedonic time dummy (Huber-IRLS, location FE) ----------------------107def pooled_td(df: pl.DataFrame, grid: list[str], robust: bool = True, fe: str = "loc_fe") -> IndexSeries:108 rhs = _rhs_for(df) + " + C(month)"109 res, m = fit_fe(df, rhs, fe, robust=robust)110 coef = res.coef111 se = dict(m.se())112 base = grid[0]113 li = np.array([0.0 if g == base else coef.get(f"C(month)[T.{g}]", np.nan) for g in grid])114 s = np.array([0.0 if g == base else se.get(f"C(month)[T.{g}]", np.nan) for g in grid])115 return IndexSeries(grid, li, se=s, n=_counts(df, grid), method="C_pooled_td" + ("" if robust else "_ols"),116 meta={"sigma_robust": res.sigma2_robust ** 0.5, "n": res.n})117118119# D. rolling time dummy with window / movement / mean splice ------------------120def rolling_td(df: pl.DataFrame, grid: list[str], window: int = 13, splice: str = "mean", robust: bool = True, fe: str = "loc_fe",121 min_obs: int = 50) -> IndexSeries:122 """Hill et al.-type rolling time dummy. Each window of `window` months is estimated; the published path is123 extended one month at a time. splice ∈ {window, movement, mean, half}."""124 rhs = _rhs_for(df) + " + C(month)"125 T = len(grid)126 path = np.full(T, np.nan)127 paths_w: list[np.ndarray] = []128 first_end = min(window, T)129 for end in range(first_end, T + 1):130 w_months = grid[end - window: end] if end >= window else grid[:end]131 sub = df.filter(pl.col("month").is_in(w_months))132 if sub.height < min_obs or sub["month"].n_unique() < 2:133 paths_w.append(np.full(len(w_months), np.nan)); continue134 res, m = fit_fe(sub, rhs, fe, robust=robust)135 coef = res.coef136 b = w_months[0]137 pw = np.array([0.0 if g == b else coef.get(f"C(month)[T.{g}]", np.nan) for g in w_months])138 paths_w.append(pw)139 if end == first_end:140 path[:end] = pw141 continue142 new_pos = end - 1143 if splice == "window":144 path[new_pos] = path[end - window] + (pw[-1] - pw[0])145 elif splice == "movement":146 path[new_pos] = path[new_pos - 1] + (pw[-1] - pw[-2])147 elif splice == "half":148 k = max(1, (window - 1) // 2)149 path[new_pos] = path[new_pos - k] + (pw[-1] - pw[-1 - k])150 else: # mean splice (Diewert–Fox / Hill et al.): average over all overlap anchors151 cands = []152 for j in range(window - 1):153 anchor = new_pos - (window - 1) + j154 if anchor >= 0 and not np.isnan(path[anchor]) and not np.isnan(pw[j]):155 cands.append(path[anchor] + (pw[-1] - pw[j]))156 path[new_pos] = np.mean(cands) if cands else np.nan157 return IndexSeries(grid, path, n=_counts(df, grid), method=f"D_rtd_w{window}_{splice}")158159160# E/F. hedonic imputation (Laspeyres / Paasche / Fisher) and GEKS ------------161def _period_fits(df: pl.DataFrame, grid: list[str], loc_col: str, min_obs: int = 40, robust: bool = False):162 """Per-month hedonic fits (OLS or Huber) with a coarse location FE; returns {month: (model, pandas data)}."""163 rhs = _rhs_for(df)164 fits = {}165 for mth in grid:166 sub = df.filter(pl.col("month") == mth)167 if sub.height < min_obs:168 continue169 try:170 res, m = fit_fe(sub, rhs, loc_col, robust=robust, n_irls=2)171 except Exception:172 continue173 fits[mth] = (m, sub)174 return fits175176177def _impute(m, data_pl: pl.DataFrame) -> np.ndarray:178 pdf = data_pl.to_pandas()179 try:180 pred = np.asarray(m.predict(newdata=pdf), float)181 except Exception:182 return np.full(len(pdf), np.nan)183 return pred184185186def hedonic_imputation(df: pl.DataFrame, grid: list[str], loc_col: str = "loc_coarse", kind: str = "fisher", chained: bool = True,187 min_obs: int = 40, robust: bool = False) -> IndexSeries:188 fits = _period_fits(df, grid, loc_col, min_obs, robust)189 T = len(grid)190 li = np.full(T, np.nan)191 if not fits:192 return IndexSeries(grid, li, n=_counts(df, grid), method=f"E_hi_{kind}")193 months = [g for g in grid if g in fits]194 li[grid.index(months[0])] = 0.0195196 def bilateral(s: str, t: str) -> float:197 ms, ds = fits[s]; mt, dt = fits[t]198 # Laspeyres: basket s, prices t vs s (double imputation)199 ps_s = _impute(ms, ds); pt_s = _impute(mt, ds)200 ok = ~np.isnan(ps_s) & ~np.isnan(pt_s)201 lasp = np.mean(pt_s[ok] - ps_s[ok]) if ok.sum() >= 10 else np.nan202 ps_t = _impute(ms, dt); pt_t = _impute(mt, dt)203 ok2 = ~np.isnan(ps_t) & ~np.isnan(pt_t)204 paas = np.mean(pt_t[ok2] - ps_t[ok2]) if ok2.sum() >= 10 else np.nan205 return {"laspeyres": lasp, "paasche": paas, "fisher": 0.5 * (lasp + paas)}[kind]206207 if chained:208 for a, b in zip(months[:-1], months[1:]):209 li[grid.index(b)] = li[grid.index(a)] + bilateral(a, b)210 else:211 for b in months[1:]:212 li[grid.index(b)] = bilateral(months[0], b)213 return IndexSeries(grid, li, n=_counts(df, grid), method=f"E_hi_{kind}{'_chained' if chained else '_direct'}", meta={"months_estimated": len(months)})214215216def geks(df: pl.DataFrame, grid: list[str], loc_col: str = "loc_coarse", window: int = 25, min_obs: int = 40, robust: bool = False) -> IndexSeries:217 """Rolling-window GEKS on bilateral double-imputation Fisher links, extended by movement splice (RYGEKS)."""218 fits = _period_fits(df, grid, loc_col, min_obs, robust)219 months = [g for g in grid if g in fits]220 T = len(months)221 li = np.full(len(grid), np.nan)222 if T < 2:223 return IndexSeries(grid, li, n=_counts(df, grid), method=f"F_geks_w{window}")224 # cache imputed prices: pred[t][s] = f_t applied to basket s225 pred: dict[tuple[str, str], np.ndarray] = {}226227 def P(t: str, s: str) -> np.ndarray:228 if (t, s) not in pred:229 pred[(t, s)] = _impute(fits[t][0], fits[s][1])230 return pred[(t, s)]231232 def fisher(s: str, t: str) -> float:233 ok = ~np.isnan(P(t, s)) & ~np.isnan(P(s, s)); lasp = np.mean((P(t, s) - P(s, s))[ok]) if ok.sum() >= 10 else np.nan234 ok2 = ~np.isnan(P(t, t)) & ~np.isnan(P(s, t)); paas = np.mean((P(t, t) - P(s, t))[ok2]) if ok2.sum() >= 10 else np.nan235 return 0.5 * (lasp + paas)236237 def geks_window(ws: list[str]) -> np.ndarray:238 k = len(ws)239 F = np.zeros((k, k))240 for i in range(k):241 for j in range(i + 1, k):242 F[i, j] = fisher(ws[i], ws[j]); F[j, i] = -F[i, j]243 F = np.where(np.isnan(F), 0.0, F)244 return F.mean(axis=0) # mean_j F[j,t] = log P_t − mean_j log P_j (GEKS, up to a constant)245 out = np.full(T, np.nan)246 first = min(window, T)247 g0 = geks_window(months[:first]); out[:first] = g0 - g0[0]248 for end in range(first + 1, T + 1):249 ws = months[end - window: end]250 g = geks_window(ws)251 out[end - 1] = out[end - 2] + (g[-1] - g[-2]) # movement splice252 for i, m in enumerate(months):253 li[grid.index(m)] = out[i]254 return IndexSeries(grid, li, n=_counts(df, grid), method=f"F_geks_w{window}_movement")255256257# H. repeat sales -------------------------------------------------------------258def repeat_sales_pairs(df: pl.DataFrame, min_hold_months: int = 6) -> pl.DataFrame:259 d = df.filter(pl.col("parcel_id").is_not_null()).sort(["parcel_id", "date"]).with_columns(260 pl.col("month").shift(1).over("parcel_id").alias("m1"), pl.col("log_price").shift(1).over("parcel_id").alias("lp1"),261 pl.col("date").shift(1).over("parcel_id").alias("d1"),262 ).filter(pl.col("m1").is_not_null())263 d = d.with_columns(((pl.col("date") - pl.col("d1")).dt.total_days() / 30.4375).alias("hold_months"))264 return d.filter(pl.col("hold_months") >= min_hold_months).select(["parcel_id", "m1", pl.col("month").alias("m2"), "lp1", pl.col("log_price").alias("lp2"), "hold_months"])265266267def repeat_sales(df: pl.DataFrame, grid: list[str], weighted: bool = True, min_hold_months: int = 6) -> IndexSeries:268 pairs = repeat_sales_pairs(df, min_hold_months)269 T = len(grid)270 li = np.full(T, np.nan)271 if pairs.height < 30:272 return IndexSeries(grid, li, method="H_repeat_sales")273 idx = {m: i for i, m in enumerate(grid)}274 i1 = np.array([idx[m] for m in pairs["m1"].to_list()]); i2 = np.array([idx[m] for m in pairs["m2"].to_list()])275 y = (pairs["lp2"] - pairs["lp1"]).to_numpy()276 n = len(y)277 rows = np.concatenate([np.arange(n), np.arange(n)]); cols = np.concatenate([i1, i2]); vals = np.concatenate([-np.ones(n), np.ones(n)])278 X = sparse.csr_matrix((vals, (rows, cols)), shape=(n, T))[:, 1:] # base = first month279 beta = lsqr(X, y, atol=1e-10, btol=1e-10)[0]280 if weighted: # Case–Shiller: var(resid) = a + b*hold281 r = y - X @ beta282 h = pairs["hold_months"].to_numpy()283 A = np.column_stack([np.ones(n), h]); ab = np.linalg.lstsq(A, r ** 2, rcond=None)[0]284 v = np.maximum(A @ ab, 1e-4); w = 1 / np.sqrt(v)285 beta = lsqr(sparse.diags(w) @ X, w * y, atol=1e-10, btol=1e-10)[0]286 li[0] = 0.0; li[1:] = beta287 cnt = np.bincount(np.concatenate([i1, i2]), minlength=T).astype(float)288 # months with no pair endpoints are not identified289 li[cnt == 0] = np.nan290 return IndexSeries(grid, li, n=cnt, method="H_repeat_sales" + ("_cs" if weighted else "_bmn"), meta={"pairs": n})291292293# I. static empirical-Bayes partial pooling of cell deviations (mixed-effects benchmark) ---------294def static_eb_deviation(resid_mean: np.ndarray, n: np.ndarray, sigma2: float) -> tuple[np.ndarray, np.ndarray, float]:295 """Random-intercept shrinkage per period: deviation_t ~ N(0, tau2); ybar_t | dev ~ N(dev, sigma2/n_t).296 tau2 by marginal likelihood (method of moments fallback). Returns (posterior mean, weight, tau2)."""297 ok = n > 0298 if ok.sum() < 2:299 return np.zeros_like(resid_mean), np.zeros_like(resid_mean), 0.0300 y = resid_mean[ok]; v = sigma2 / n[ok]301 tau2 = max(np.var(y) - np.mean(v), 1e-8)302 for _ in range(50): # DerSimonian–Laird style iteration303 w = 1 / (v + tau2); mu = np.sum(w * y) / np.sum(w)304 tau2_new = max(np.sum(w ** 2 * ((y - mu) ** 2 - v)) / np.sum(w ** 2), 1e-8)305 if abs(tau2_new - tau2) < 1e-10:306 break307 tau2 = tau2_new308 wgt = np.zeros_like(resid_mean); post = np.zeros_like(resid_mean)309 wgt[ok] = tau2 / (tau2 + v); post[ok] = wgt[ok] * y310 return post, wgt, tau2311312313# J. out-of-time gradient boosting residual index (diagnostic only) ------------------------------314def gbm_residual_index(df: pl.DataFrame, grid: list[str], retrain_every: int = 3, min_train: int = 2000, seed: int = 20261021) -> IndexSeries:315 """Train LightGBM on log price ~ characteristics + coordinates using only months < t (expanding window,316 retrained every `retrain_every` months), then index_t = mean out-of-sample residual (chain of monthly317 quality-adjusted price levels relative to the training-period average). Predictive accuracy is NOT index accuracy."""318 import lightgbm as lgb319 feats = ["lfa", "lla", "age_f", "fa_missing", "la_missing", "age_missing", "x_m", "y_m"]320 cat = ["storeys_c", "genre_c", "subtype_c", "units_c"]321 d = df.with_columns([pl.col(c).cast(pl.Categorical).to_physical().alias(c + "_i") for c in cat])322 X_cols = feats + [c + "_i" for c in cat]323 T = len(grid)324 li = np.full(T, np.nan)325 model = None326 train_mean_month = None327 for i, mth in enumerate(grid):328 past = d.filter(pl.col("month") < mth)329 if past.height < min_train:330 continue331 if model is None or i % retrain_every == 0:332 # detrend training target with month means so the model learns quality, not time333 mm = past.group_by("month").agg(pl.col("log_price").mean().alias("_mm"))334 tr = past.join(mm, on="month")335 yt = (tr["log_price"] - tr["_mm"]).to_numpy()336 model = lgb.LGBMRegressor(n_estimators=400, learning_rate=0.05, num_leaves=63, min_child_samples=20, subsample=0.8, subsample_freq=1,337 colsample_bytree=0.8, objective="l1", random_state=seed, verbose=-1)338 model.fit(tr.select(X_cols).to_numpy(), yt)339 train_mean_month = mm340 cur = d.filter(pl.col("month") == mth)341 if cur.height == 0:342 continue343 pred = model.predict(cur.select(X_cols).to_numpy())344 li[i] = float(np.mean(cur["log_price"].to_numpy() - pred))345 return IndexSeries(grid, li, n=_counts(df, grid), method="J_gbm_oot_residual")346347348# G. HMTS — hedonic multilateral time-series (Ishaak, Ouwehand & Remøy 2024, JOS 40(3)) ---------------------349def hmts(df: pl.DataFrame, grid: list[str], loc_col: str = "loc_coarse", window: int = 13, min_obs: int = 40, start_periods: int = 20) -> IndexSeries:350 """For every base period b, the series p̂_{t(b)} = exp(α̂_t + β̂_t' x̄_b), t = b … b+w−1 (period-t hedonic coefficients on the351 base-period mean characteristics; here evaluated as the mean imputed log price on the base basket, eq. 4), is smoothed by a352 smooth-trend state-space model (local linear trend with zero level disturbance, ML hyperparameters, Kalman smoother; eqs. 9–12)353 fitted on log p̂; each smoothed series is indexed on its own first period (eq. 17) and the published index is the geometric mean354 over the w most recent base-period series (eq. 18). Periods beyond the first `start_periods` are 'limited-revisable' (w−1 preliminary355 periods). See docs/literature/hmts_method.md. Monthly adaptation: w = `window` (paper: w = 4 quarters)."""356 from scipy.optimize import minimize357358 from sharp_hpi.statespace import _filter_llt359 fits = _period_fits(df, grid, loc_col, min_obs, robust=False)360 months = [g for g in grid if g in fits]361 T = len(months)362 li = np.full(len(grid), np.nan)363 if T < 3:364 return IndexSeries(grid, li, n=_counts(df, grid), method=f"G_hmts_w{window}")365 # imputed log prices of basket b at period-t prices (row vectors, cached); links use the common support of366 # consecutive periods so that rows unpredictable in one period (unseen location levels) do not shift the mean367 cache: dict[tuple[str, str], np.ndarray] = {}368369 def P(t: str, b: str) -> np.ndarray:370 if (t, b) not in cache:371 cache[(t, b)] = _impute(fits[t][0], fits[b][1])372 return cache[(t, b)]373374 def link(b: str, t0: str, t1: str) -> float:375 a0, a1 = P(t0, b), P(t1, b); ok = np.isfinite(a0) & np.isfinite(a1)376 return float(np.mean(a1[ok] - a0[ok])) if ok.sum() >= 10 else np.nan377378 def smooth_trend(y: np.ndarray) -> np.ndarray:379 ok = ~np.isnan(y)380 if ok.sum() < 3:381 return y382 yy = np.where(ok, y, 0.0)383384 def nll(th):385 th = np.clip(th, -18.0, 0.0)386 se2, sz2 = np.exp(th)387 r = np.where(ok, se2, np.inf)388 v = -_filter_llt(yy, r, 0.0, sz2, 1.0)[5]389 return v if np.isfinite(v) else 1e12390 best = minimize(nll, x0=[np.log(1e-3), np.log(1e-4)], method="Nelder-Mead", options={"maxiter": 300, "xatol": 1e-3, "fatol": 1e-4})391 se2, sz2 = np.exp(np.clip(best.x, -18.0, 0.0))392 s = _filter_llt(yy, np.where(ok, se2, np.inf), 0.0, sz2, 1.0)[2]393 if not np.all(np.isfinite(s[ok])) or np.nanmax(np.abs(s - np.nanmean(s))) > 1.0:394 return y # degenerate fit on a very short / sparse series: keep the raw imputed series395 return s396397 # base-period series: each extended to at most `window` periods, then frozen (Table 5)398 series: dict[int, np.ndarray] = {}399 for bi in range(T):400 end = min(bi + window, T)401 lk = [0.0] + [link(months[bi], months[t - 1], months[t]) for t in range(bi + 1, end)]402 y = np.cumsum(np.nan_to_num(lk)); y[np.isnan(lk)] = np.nan403 s = smooth_trend(y)404 series[bi] = s - s[0] # indexed on its first period (eq. 17)405 # final index: for period t, geometric mean over the base-period series that cover t; chain across the moving set (eq. 18–19)406 out = np.full(T, np.nan); out[0] = 0.0407 for t in range(1, T):408 links = []409 for bi, s in series.items():410 if bi <= t - 1 and t - bi < len(s):411 links.append(s[t - bi] - s[t - 1 - bi])412 out[t] = out[t - 1] + (np.nanmean(links) if links else np.nan)413 for i, m in enumerate(months):414 li[grid.index(m)] = out[i]415 return IndexSeries(grid, li, n=_counts(df, grid), method=f"G_hmts_w{window}", meta={"months_estimated": T})416