TeX 55.1%
Python 44.2%
1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Hedonic engine: feature construction, robust (Huber-IRLS) fixed-effects regressions.34Headline specification (structural + location only; no market variables):5 log P = s(log floor area) + s(log land area) + age bins + storeys + genre + subtype [+ n_units for plex]6 + location FE (fine: neighbourhood unit / 1-km cell) + time effects7Missing characteristics: conditional-median imputation within (class × municipality) → (class × region) → class,8with time-stable missingness indicators (floor-area and land-area missingness are checked for drift in9results/profile/completeness_by_month.csv).10"""11from __future__ import annotations1213from dataclasses import dataclass1415import warnings1617import numpy as np18import polars as pl19import pyfixest as pf2021warnings.filterwarnings("ignore", message=".*singleton fixed effect.*")22warnings.filterwarnings("ignore", message=".*multicollinearity.*")23warnings.filterwarnings("ignore", category=RuntimeWarning, module="sharp_hpi")2425AGE_EDGES = [1, 5, 10, 20, 40, 60, 80, 120]26AGE_LABELS = ["a00_0-1", "a01_2-5", "a02_6-10", "a03_11-20", "a04_21-40", "a05_41-60", "a06_61-80", "a07_81-120", "a08_120+"]27HUBER_K = 1.34528N_IRLS = 3293031def _impute(df: pl.DataFrame, col: str, filled: str, cls: str = "property_class") -> pl.DataFrame:32 m1 = df.group_by([cls, "mun_code"]).agg(pl.col(col).median().alias("_m1"))33 m2 = df.group_by([cls, "region_code"]).agg(pl.col(col).median().alias("_m2"))34 m3 = df.group_by(cls).agg(pl.col(col).median().alias("_m3"))35 return (df.join(m1, on=[cls, "mun_code"], how="left").join(m2, on=[cls, "region_code"], how="left").join(m3, on=cls, how="left")36 .with_columns(pl.coalesce(pl.col(col), pl.col("_m1"), pl.col("_m2"), pl.col("_m3")).alias(filled)).drop("_m1", "_m2", "_m3"))373839def build_features(df: pl.DataFrame) -> pl.DataFrame:40 """Attach modelling features to the analytic layer (idempotent)."""41 out = df.with_columns(42 pl.when(pl.col("floor_area_m2") > 10).then(pl.col("floor_area_m2")).otherwise(None).alias("_fa"),43 pl.when(pl.col("land_area_m2") > 10).then(pl.col("land_area_m2")).otherwise(None).alias("_la"),44 pl.when(pl.col("age").is_null() | (pl.col("age") > 200)).then(None).otherwise(pl.col("age")).alias("_age"),45 ).with_columns(46 pl.col("_fa").is_null().alias("fa_missing"), pl.col("_la").is_null().alias("la_missing"), pl.col("_age").is_null().alias("age_missing"),47 )48 out = _impute(out, "_fa", "fa_f")49 out = _impute(out, "_la", "la_f")50 out = _impute(out, "_age", "age_f")51 age = pl.col("age_f")52 age_bin = pl.when(age <= 1).then(pl.lit(AGE_LABELS[0]))53 lo = 154 for edge, lab in zip(AGE_EDGES[1:], AGE_LABELS[1:-1]):55 age_bin = age_bin.when((age > lo) & (age <= edge)).then(pl.lit(lab))56 lo = edge57 age_bin = age_bin.otherwise(pl.lit(AGE_LABELS[-1]))58 out = out.with_columns(59 pl.col("fa_f").log().alias("log_fa"), pl.col("la_f").log().alias("log_la"), age_bin.alias("age_bin"),60 pl.col("storeys").fill_null(1.0).clip(1, 4).round(0).cast(pl.Int32).cast(pl.Utf8).alias("storeys_c"),61 pl.col("subtype").fill_null("na").alias("subtype_c"),62 pl.when(pl.col("property_class") == "plex").then(pl.col("n_units").clip(2, 5)).otherwise(1).cast(pl.Utf8).alias("units_c"),63 pl.col("genre").alias("genre_c"),64 pl.concat_str([pl.col("property_class"), pl.lit(":"), pl.col("month")]).alias("class_month"),65 ).drop("_fa", "_la", "_age")66 # log floor area: linear + squared + cubic around the class median (smooth nonlinearity); land area linear + squared67 med = out.group_by("property_class").agg(pl.col("log_fa").median().alias("_mfa"), pl.col("log_la").median().alias("_mla"))68 out = out.join(med, on="property_class", how="left").with_columns(69 (pl.col("log_fa") - pl.col("_mfa")).alias("lfa"), (pl.col("log_la") - pl.col("_mla")).alias("lla"),70 ).with_columns((pl.col("lfa") ** 2).alias("lfa2"), (pl.col("lfa") ** 3).alias("lfa3"), (pl.col("lla") ** 2).alias("lla2")).drop("_mfa", "_mla")71 return out727374def assign_loc_fe(df: pl.DataFrame, min_n: int = 15, levels: tuple[str, ...] = ("grid500", "grid1k", "grid2k", "grid5k", "mun_code", "mrc_code"),75 group: tuple[str, ...] = ("property_class",)) -> pl.DataFrame:76 """Finest nested spatial cell with at least `min_n` rows (within `group`), coarsening upward.7778 The partition depends only on sale locations (not prices) and is fixed for the estimation sample it is built on;79 real-time vintages rebuild it from the data available at the cutoff.80 """81 g = list(group)82 key = pl.lit(None, dtype=pl.Utf8)83 chosen = pl.lit(None, dtype=pl.Utf8)84 out = df85 for lvl in levels:86 n = pl.len().over(g + [lvl])87 chosen = pl.when(chosen.is_null() & (n >= min_n)).then(pl.concat_str([pl.lit(lvl + ":"), pl.col(lvl)])).otherwise(chosen)88 out = out.with_columns(chosen.alias("loc_fe"))89 out = out.with_columns(pl.when(pl.col("loc_fe").is_null()).then(pl.concat_str([pl.lit("mrc_code:"), pl.col("mrc_code")])).otherwise(pl.col("loc_fe")).alias("loc_fe"))90 # a cell that received rows only through fall-back might still be < min_n at a coarse level; that is acceptable (no singletons by91 # construction only if the group has >= min_n rows at the top level) — guard below.92 return out.with_columns(pl.col("loc_fe").str.slice(0, 8).alias("loc_fe_level"))939495# Right-hand sides -----------------------------------------------------------96RHS_STRUCT = "lfa + lfa2 + lfa3 + fa_missing + lla + lla2 + la_missing + C(age_bin) + age_missing + C(storeys_c) + C(genre_c) + C(subtype_c) + C(units_c)"97RHS_STRUCT_BY_CLASS = "lfa + lfa2 + lfa3 + fa_missing + lla + lla2 + la_missing + C(age_bin) + age_missing + C(storeys_c) + C(genre_c) + C(subtype_c) + C(units_c)"9899100@dataclass101class FitResult:102 coef: dict103 resid: np.ndarray104 fitted: np.ndarray105 weights: np.ndarray106 sigma2_robust: float107 n: int108 r2_within: float | None109110111def huber_weights(resid: np.ndarray, k: float = HUBER_K) -> tuple[np.ndarray, float]:112 s = 1.4826 * np.median(np.abs(resid - np.median(resid)))113 s = max(s, 1e-6)114 u = np.abs(resid) / s115 w = np.where(u <= k, 1.0, k / u)116 return w, s117118119def fit_fe(df: pl.DataFrame, rhs: str, fe: str, y: str = "log_price", robust: bool = True, n_irls: int = N_IRLS,120 base_weights: np.ndarray | None = None) -> tuple[FitResult, pf.estimation.Feols]:121 """Huber-IRLS fixed-effects regression (pyfixest). Returns residuals aligned with df rows."""122 n_all = df.height123 w_all = np.ones(n_all) if base_weights is None else np.asarray(base_weights, float)124 # iterative singleton removal on every FE factor so residuals stay aligned (dropped rows get NaN)125 keep = np.ones(n_all, bool)126 fe_cols = [c.strip() for c in fe.split("+")] if fe else []127 if fe_cols:128 sub = df.select(fe_cols).with_row_index("_i")129 while True:130 cur = sub.filter(pl.Series(keep)[sub["_i"]])131 drop = np.zeros(n_all, bool)132 for c in fe_cols:133 cnt = cur.group_by(c).agg(pl.len().alias("_n"))134 ones = cur.join(cnt.filter(pl.col("_n") == 1), on=c, how="semi")["_i"].to_numpy()135 drop[ones] = True136 if not drop.any():137 break138 keep &= ~drop139 pdf = df.filter(pl.Series(keep)).to_pandas()140 w = w_all[keep]141 fml = f"{y} ~ {rhs} | {fe}" if fe else f"{y} ~ {rhs}"142 pdf["_w"] = w143 m = pf.feols(fml, data=pdf, weights="_w", vcov="iid")144 resid = np.asarray(m.resid(), float)145 if len(resid) != len(pdf):146 raise ValueError(f"rows dropped by the FE regression ({len(pdf) - len(resid)}) after singleton removal")147 if robust:148 for _ in range(n_irls):149 hw, _s = huber_weights(resid)150 pdf["_w"] = w * hw151 m = pf.feols(fml, data=pdf, weights="_w", vcov="iid")152 resid = np.asarray(m.resid(), float)153 hw, s = huber_weights(resid)154 resid_all = np.full(n_all, np.nan); resid_all[keep] = resid155 fitted_all = np.full(n_all, np.nan); fitted_all[keep] = np.asarray(pdf[y].values, float) - resid156 w_out = np.zeros(n_all); w_out[keep] = pdf["_w"].values157 r2w = float(m._r2_within) if hasattr(m, "_r2_within") and m._r2_within is not None else None158 return FitResult(coef=dict(m.coef()), resid=resid_all, fitted=fitted_all, weights=w_out, sigma2_robust=float(s ** 2), n=int(keep.sum()), r2_within=r2w), m159