UQO Working Paper No. 9 — A grand hedonic model of the Canadian housing market: decomposing structure and location value.
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 06 — Additional analyses added in the October 2026 revision.45All results are written to ``results/reproduced/`` and consumed by scripts 046and 05. Each block is a self-contained function.78 shapley.csv Shapley / Owen decomposition of R^2 across regressor blocks9 se_comparison.csv M5 key coefficients under alternative variance estimators,10 incl. a wild-cluster (Rademacher) bootstrap at the province level11 spline.csv piecewise-linear (spline) living-area elasticities vs. quadratic12 spatial_holdout.json FSA-level spatial holdout + repeated random holdouts13 moran_sensitivity.csv Moran's I for k in {5, 10, 20} and two subsamples14 gradient_metro.json urban gradient with metro-specific intercepts1516Usage: python scripts/06_revision_analyses.py17"""18import itertools19import json20import math21import sys22import warnings23from pathlib import Path2425import numpy as np26import pandas as pd27import statsmodels.api as sm28from sklearn.neighbors import NearestNeighbors2930sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))3132from wp9 import models, sample # noqa: E40233from wp9.config import METROS, REPRODUCED, SEED_MORAN, SEED_OOS, STRUCT, ensure_dirs # noqa: E4023435warnings.filterwarnings("ignore")36SEED_REV = 2026373839# --------------------------------------------------------------------------40def _r2(y, X):41 """OLS R^2 with an intercept-only baseline (X may already contain a constant)."""42 if "const" not in X.columns:43 X = sm.add_constant(X)44 return float(sm.OLS(y, X).fit().rsquared)454647def _r2_absorb(y, X, fsa):48 """R^2 of a model that absorbs FSA fixed effects (full-dummy equivalent)."""49 from linearmodels.iv.absorbing import AbsorbingLS50 X = X.copy()51 if "const" not in X.columns:52 X = sm.add_constant(X)53 res = AbsorbingLS(y, X, absorb=fsa.to_frame().astype("category"),54 drop_absorbed=True).fit(cov_type="unadjusted")55 return float(res.rsquared)565758def shapley_decomposition(s: pd.DataFrame, label: str) -> list[dict]:59 """Shapley values of R^2 for the blocks {structure, type/ownership, province, FSA}.6061 Province dummies are nested in the FSA effects, so any coalition that62 contains FSA has the same R^2 with or without the province block; the63 Shapley value of the province block then measures what provinces add when64 FSAs are *not* present, averaged over orderings. The Owen value groups65 province and FSA into a single 'location' block.66 """67 y = s["ln_price"]68 blocks = {69 "S": s[STRUCT].astype(float),70 "T": pd.concat([pd.get_dummies(s[c], prefix=p, drop_first=True).astype(float)71 for c, p in (("btype_c", "bt"), ("own_c", "ow"), ("cat", "cat"))], axis=1),72 "P": pd.get_dummies(s["prov"], prefix="pv", drop_first=True).astype(float),73 }74 names = ["S", "T", "P", "F"]75 cache = {}7677 def value(coal):78 key = tuple(sorted(coal))79 if key in cache:80 return cache[key]81 cols = [blocks[b] for b in key if b in blocks]82 X = pd.concat(cols, axis=1) if cols else pd.DataFrame(index=s.index)83 if "F" in key:84 X = X.drop(columns=[c for c in X.columns if c.startswith("pv_")], errors="ignore")85 v = _r2_absorb(y, X, s["fsa_c"]) if len(X.columns) else _r2_absorb(86 y, pd.DataFrame({"const": np.ones(len(s))}, index=s.index), s["fsa_c"])87 else:88 v = _r2(y, X) if len(X.columns) else 0.089 cache[key] = v90 return v9192 n = len(names)93 shap = {b: 0.0 for b in names}94 for b in names:95 others = [o for o in names if o != b]96 for k in range(len(others) + 1):97 for coal in itertools.combinations(others, k):98 w = math.factorial(k) * math.factorial(n - k - 1) / math.factorial(n)99 shap[b] += w * (value(coal + (b,)) - value(coal))100 total = value(tuple(names))101102 # Owen value: groups {S}, {T}, {L = P, F}; Shapley over the 3 groups.103 groups = {"S": ("S",), "T": ("T",), "L": ("P", "F")}104 gnames = list(groups)105 owen = {g: 0.0 for g in gnames}106 for g in gnames:107 others = [o for o in gnames if o != g]108 for k in range(len(others) + 1):109 for coal in itertools.combinations(others, k):110 w = math.factorial(k) * math.factorial(3 - k - 1) / math.factorial(3)111 with_g = sum((groups[x] for x in coal + (g,)), ())112 without = sum((groups[x] for x in coal), ())113 owen[g] += w * (value(with_g) - value(without))114115 seq = {"S": value(("S",)),116 "T": value(("S", "T")) - value(("S",)),117 "P": value(("S", "T", "P")) - value(("S", "T")),118 "F": value(("S", "T", "P", "F")) - value(("S", "T", "P"))}119 rows = []120 for b, lab in (("S", "Structure"), ("T", "Dwelling type and ownership"),121 ("P", "Province"), ("F", "Neighbourhood (FSA)")):122 rows.append({"sample": label, "block": lab, "sequential": seq[b], "shapley": shap[b],123 "owen": owen["L"] if b == "F" else (np.nan if b == "P" else owen[b]),124 "total_r2": total, "n": len(s)})125 return rows126127128# --------------------------------------------------------------------------129def _within(df, cols, by):130 """Demean columns by group (one-way fixed-effect transformation)."""131 return df[cols] - df.groupby(by)[cols].transform("mean")132133134def se_comparison(s: pd.DataFrame) -> pd.DataFrame:135 """M5 key coefficients under HC1, FSA-cluster, province-cluster and a136 province-level wild-cluster bootstrap (Rademacher, full enumeration)."""137 X = models.design(s).drop(columns="const")138 cols = X.columns.tolist()139 data = pd.concat([s[["ln_price", "fsa_c", "prov"]], X], axis=1)140 dd = _within(data, ["ln_price"] + cols, "fsa_c")141 y, Xd = dd["ln_price"].values, dd[cols].values142 keep = Xd.std(axis=0) > 1e-12 # drop columns absorbed by the FSA effects143 Xd, cols = Xd[:, keep], [c for c, k in zip(cols, keep) if k]144 n_fe = s["fsa_c"].nunique()145 ols = sm.OLS(y, Xd)146 base = ols.fit()147 dof_corr = (len(y) - 1) / (len(y) - Xd.shape[1] - n_fe) # absorbed-FE dof adjustment148 hc1 = ols.fit(cov_type="HC1")149 cl_fsa = ols.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(s["fsa_c"])[0]})150 cl_prov = ols.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(s["prov"])[0]})151 cl_2way = ols.fit(cov_type="cluster",152 cov_kwds={"groups": np.column_stack([pd.factorize(s["fsa_c"])[0],153 pd.factorize(s["prov"])[0]])})154 rows = []155 for key in ("ln_living", "bathrooms", "half_baths", "bedrooms", "ln_lot"):156 j = cols.index(key)157 rows.append({"variable": key, "coef": base.params[j],158 "se_hc1": hc1.bse[j] * math.sqrt(dof_corr),159 "se_fsa": cl_fsa.bse[j] * math.sqrt(dof_corr),160 "se_prov": cl_prov.bse[j],161 "se_twoway": cl_2way.bse[j] * math.sqrt(dof_corr)})162 # Wild-cluster restricted bootstrap-t at the province level (G = 9 -> 512 draws)163 groups = pd.factorize(s["prov"])[0]164 G = groups.max() + 1165 signs = np.array(list(itertools.product([-1.0, 1.0], repeat=G)))166 XtX_inv = np.linalg.inv(Xd.T @ Xd)167 for r in rows:168 j = cols.index(r["variable"])169 t_obs = r["coef"] / r["se_prov"]170 mask = np.ones(Xd.shape[1], bool); mask[j] = False171 Xr = Xd[:, mask]172 br = np.linalg.lstsq(Xr, y, rcond=None)[0]173 yhat_r, e_r = Xr @ br, y - Xr @ br174 t_boot = np.empty(len(signs))175 for b, w in enumerate(signs):176 yb = yhat_r + e_r * w[groups]177 bb = XtX_inv @ (Xd.T @ yb)178 ub = yb - Xd @ bb179 # cluster-robust variance of bb[j] by province180 meat = np.zeros((Xd.shape[1], Xd.shape[1]))181 for g in range(G):182 sg = Xd[groups == g].T @ ub[groups == g]183 meat += np.outer(sg, sg)184 V = XtX_inv @ meat @ XtX_inv * (G / (G - 1))185 t_boot[b] = bb[j] / math.sqrt(V[j, j])186 r["p_wild_prov"] = float(np.mean(np.abs(t_boot) >= abs(t_obs)))187 r["t_prov"] = float(t_obs)188 out = pd.DataFrame(rows)189 out.attrs["n_fsa_clusters"] = int(n_fe)190 out["n_fsa_clusters"] = n_fe191 out["n_prov_clusters"] = G192 return out193194195# --------------------------------------------------------------------------196def spline_elasticity(s: pd.DataFrame) -> pd.DataFrame:197 """Piecewise-linear spline in ln(area) with knots at the size quintiles."""198 knots = s["ln_living"].quantile([0.2, 0.4, 0.6, 0.8]).values199 extra = pd.DataFrame({f"hinge{k}": np.maximum(s["ln_living"] - kn, 0.0)200 for k, kn in enumerate(knots)}, index=s.index)201 res, cols = models.fit_absorbing(s, extra=extra)202 params = pd.Series(np.asarray(res.params).ravel(), index=res.params.index)203 cov = pd.DataFrame(np.asarray(res.cov), index=res.params.index, columns=res.params.index)204 names = ["ln_living"] + [f"hinge{k}" for k in range(len(knots))]205 edges = np.concatenate([[s["ln_living"].min()], knots, [s["ln_living"].max()]])206 rows = []207 for seg in range(len(knots) + 1):208 active = names[: seg + 1]209 g = np.array([1.0 if nm in active else 0.0 for nm in names])210 el = float(g @ params[names].values)211 se = float(np.sqrt(g @ cov.loc[names, names].values @ g))212 rows.append({"segment": seg + 1, "lo_m2": math.exp(edges[seg]), "hi_m2": math.exp(edges[seg + 1]),213 "elasticity": el, "se": se, "r2_spline": float(res.rsquared)})214 # quadratic and linear fits for comparison215 res_q, _ = models.fit_absorbing(s, extra=(s["ln_living"] ** 2).to_frame("ln_living2"))216 res_l, _ = models.fit_absorbing(s)217 out = pd.DataFrame(rows)218 out["r2_quadratic"] = float(res_q.rsquared)219 out["r2_linear"] = float(res_l.rsquared)220 return out221222223# --------------------------------------------------------------------------224def _score(test, pred_log, smear):225 err = test["ln_price"].values - pred_log226 ss_tot = np.sum((test["ln_price"].values - test["ln_price"].mean()) ** 2)227 ape = np.abs(np.exp(pred_log) * smear - test["price_cad"].values) / test["price_cad"].values228 return {"n_test": int(len(test)), "oos_r2": float(1 - np.sum(err ** 2) / ss_tot),229 "rmse_log": float(np.sqrt(np.mean(err ** 2))),230 "median_ape": float(np.median(ape) * 100), "mean_ape": float(np.mean(ape) * 100),231 # mean APE excluding the top 1% of errors: a handful of listings with232 # implausible parking/storey counts (e.g. 3,620 parking spaces) produce233 # astronomically large predictions in some splits234 "mean_ape_trim": float(np.mean(np.sort(ape)[: int(0.99 * len(ape))]) * 100),235 "within10": float(np.mean(ape <= 0.10) * 100), "within20": float(np.mean(ape <= 0.20) * 100)}236237238def spatial_holdout(s: pd.DataFrame) -> dict:239 """Hold out 20% of FSAs entirely; predict their listings with (a) the240 province-mean fixed effect and (b) an inverse-distance average of the five241 nearest training FSAs' effects. Also: five repeated random 80/20 splits."""242 out = {}243 rng = np.random.default_rng(SEED_REV)244 real = s[~s["fsa_c"].str.endswith("_other")]245 fsas = real["fsa_c"].unique()246 reps = []247 for rep in range(5):248 held = rng.choice(fsas, size=int(0.2 * len(fsas)), replace=False)249 train = real[~real["fsa_c"].isin(held)]250 test = real[real["fsa_c"].isin(held)]251 res, cols = models.fit_absorbing(train)252 fe = models.fsa_fixed_effects(train, res, cols)253 beta = pd.Series(np.asarray(res.params).ravel(), index=res.params.index)254 xb_test = models.design(test).reindex(columns=beta.index, fill_value=0.0).values @ beta.values255 pred_train = models.predict_with_fe(train, res, cols, fe)256 smear = models.duan_smearing(train["ln_price"].values - pred_train)257 # (a) province mean of the training FSA effects258 fe_prov = (fe.to_frame("fe").join(train.groupby("fsa_c")["prov"].first())259 .groupby("prov")["fe"].mean())260 pred_a = xb_test + test["prov"].map(fe_prov).fillna(fe.mean()).values261 # (b) inverse-distance weighted mean of the 5 nearest training FSAs262 cent = train.groupby("fsa_c")[["lat", "lon"]].mean()263 test_cent = test.groupby("fsa_c")[["lat", "lon"]].mean()264 nn = NearestNeighbors(n_neighbors=5).fit(np.radians(cent.values))265 dist, idx = nn.kneighbors(np.radians(test_cent.values))266 dist_km = 6371.0 * dist267 w = 1.0 / np.maximum(dist_km, 1.0)268 fe_train = fe.reindex(cent.index).values269 fe_nn = pd.Series((fe_train[idx] * w).sum(1) / w.sum(1), index=test_cent.index)270 pred_b = xb_test + test["fsa_c"].map(fe_nn).values271 reps.append({"rep": rep, "n_fsa_held": int(len(held)),272 "prov_mean": _score(test, pred_a, smear), "nn_idw": _score(test, pred_b, smear)})273 out["fsa_holdout"] = reps274 for key in ("prov_mean", "nn_idw"):275 out[f"fsa_holdout_mean_{key}"] = {m: float(np.mean([r[key][m] for r in reps]))276 for m in reps[0][key]}277 out[f"fsa_holdout_sd_{key}"] = {m: float(np.std([r[key][m] for r in reps], ddof=1))278 for m in reps[0][key]}279 # repeated random 80/20 splits on common support (5 seeds)280 rand = []281 for seed in (SEED_OOS, 1, 2, 3, 4):282 rs = np.random.RandomState(seed)283 idx = rs.permutation(len(s)); cut = int(0.8 * len(s))284 train, test = s.iloc[idx[:cut]], s.iloc[idx[cut:]]285 common = set(train["fsa_c"]).intersection(test["fsa_c"])286 train, test = train[train["fsa_c"].isin(common)], test[test["fsa_c"].isin(common)]287 res, cols = models.fit_absorbing(train)288 fe = models.fsa_fixed_effects(train, res, cols)289 smear = models.duan_smearing(train["ln_price"].values - models.predict_with_fe(train, res, cols, fe))290 rand.append(_score(test, models.predict_with_fe(test, res, cols, fe), smear))291 out["random_repeats"] = rand292 out["random_mean"] = {m: float(np.mean([r[m] for r in rand])) for m in rand[0]}293 out["random_sd"] = {m: float(np.std([r[m] for r in rand], ddof=1)) for m in rand[0]}294 return out295296297# --------------------------------------------------------------------------298def moran_sensitivity(s: pd.DataFrame) -> pd.DataFrame:299 """Moran's I of structure-only and grand-model residuals for several k and subsamples."""300 grand = pd.read_parquet(REPRODUCED / "grand_model.parquet")301 X = models.design(s)302 struct_resid = s["ln_price"] - sm.OLS(s["ln_price"], X).fit().predict(X)303 rows = []304 for seed in (SEED_MORAN, SEED_REV):305 sub = s.sample(15000, random_state=seed)306 coords = np.radians(sub[["lat", "lon"]].values)307 for k in (5, 10, 20):308 nn = NearestNeighbors(n_neighbors=k + 1, metric="haversine").fit(coords)309 _, idx = nn.kneighbors(coords); idx = idx[:, 1:]310 for label, vals in (("structure_only", struct_resid.loc[sub.index].values),311 ("grand_model", grand.loc[sub.index, "resid_grand"].values)):312 z = vals - vals.mean(); lag = z[idx].mean(axis=1)313 stat = np.sum(z * lag) / np.sum(z ** 2)314 rng = np.random.default_rng(seed + k)315 perms = np.array([np.sum((zp := rng.permutation(z)) * zp[idx].mean(axis=1)) / np.sum(zp ** 2)316 for _ in range(199)])317 rows.append({"seed": seed, "k": k, "model": label, "I": float(stat),318 "z": float((stat - perms.mean()) / perms.std()),319 "p_perm": float((np.sum(np.abs(perms) >= abs(stat)) + 1) / 200)})320 return pd.DataFrame(rows)321322323# --------------------------------------------------------------------------324def gradient_with_metro_fe(s: pd.DataFrame) -> dict:325 """Location premium on ln(distance) with nearest-metro intercepts, and a quadratic form."""326 X = models.design(s)327 prem = s["ln_price"] - sm.OLS(s["ln_price"], X).fit().predict(X)328 lat, lon = s["lat"].values, s["lon"].values329 rad = np.pi / 180330 dists = {}331 for name, (mlat, mlon) in METROS.items():332 a = (np.sin((mlat - lat) * rad / 2) ** 2333 + np.cos(lat * rad) * np.cos(mlat * rad) * np.sin((mlon - lon) * rad / 2) ** 2)334 dists[name] = 2 * 6371.0 * np.arcsin(np.sqrt(a))335 D = pd.DataFrame(dists, index=s.index)336 dist = D.min(axis=1); metro = D.idxmin(axis=1)337 keep = dist > 0.5338 ld = np.log(dist[keep])339 base = sm.OLS(prem[keep], sm.add_constant(ld)).fit(cov_type="cluster",340 cov_kwds={"groups": pd.factorize(s["fsa_c"][keep])[0]})341 Xm = pd.concat([ld.rename("ln_dist"), pd.get_dummies(metro[keep], drop_first=True).astype(float)], axis=1)342 mfe = sm.OLS(prem[keep], sm.add_constant(Xm)).fit(cov_type="cluster",343 cov_kwds={"groups": pd.factorize(s["fsa_c"][keep])[0]})344 Xq = pd.DataFrame({"ln_dist": ld, "ln_dist2": ld ** 2})345 quad = sm.OLS(prem[keep], sm.add_constant(Xq)).fit(cov_type="cluster",346 cov_kwds={"groups": pd.factorize(s["fsa_c"][keep])[0]})347 # metro-specific slopes348 slopes = {}349 for m, g in pd.DataFrame({"prem": prem[keep], "ld": ld, "metro": metro[keep]}).groupby("metro"):350 f = sm.OLS(g["prem"], sm.add_constant(g["ld"])).fit()351 slopes[m] = {"beta": float(f.params.iloc[1]), "n": int(len(g))}352 return {"pooled": {"beta": float(base.params.iloc[1]), "se_fsa": float(base.bse.iloc[1]), "r2": float(base.rsquared)},353 "metro_fe": {"beta": float(mfe.params["ln_dist"]), "se_fsa": float(mfe.bse["ln_dist"]), "r2": float(mfe.rsquared)},354 "quadratic": {"b1": float(quad.params["ln_dist"]), "b2": float(quad.params["ln_dist2"]),355 "se1": float(quad.bse["ln_dist"]), "se2": float(quad.bse["ln_dist2"]), "r2": float(quad.rsquared)},356 "metro_slopes": slopes, "share_within_40km": float((dist <= 40).mean()),357 "share_beyond_160km": float((dist > 160).mean())}358359360# --------------------------------------------------------------------------361def main() -> None:362 ensure_dirs()363 s = sample.load_sample()364365 rows = shapley_decomposition(s, "All residential") + \366 shapley_decomposition(s[s["cat"] == "house"], "Houses")367 pd.DataFrame(rows).to_csv(REPRODUCED / "shapley.csv", index=False)368 print(pd.DataFrame(rows).round(3).to_string())369370 sec = se_comparison(s)371 sec.to_csv(REPRODUCED / "se_comparison.csv", index=False)372 print(sec.round(4).to_string())373374 spl = spline_elasticity(s)375 spl.to_csv(REPRODUCED / "spline.csv", index=False)376 print(spl.round(3).to_string())377378 sh = spatial_holdout(s)379 json.dump(sh, open(REPRODUCED / "spatial_holdout.json", "w"), indent=2)380 print("spatial holdout:", {k: {m: round(v, 3) for m, v in d.items()} for k, d in sh.items() if k.startswith(("fsa_holdout_mean", "random_mean"))})381382 ms = moran_sensitivity(s)383 ms.to_csv(REPRODUCED / "moran_sensitivity.csv", index=False)384 print(ms.round(3).to_string())385386 gm = gradient_with_metro_fe(s)387 json.dump(gm, open(REPRODUCED / "gradient_metro.json", "w"), indent=2)388 print("gradient:", gm["pooled"], gm["metro_fe"], gm["quadratic"])389 print("ALL REVISION ANALYSES DONE")390391392if __name__ == "__main__":393 main()394