#!/usr/bin/env python3 # Author: Simon-Pierre Boucher — contact@spboucher.ai """Step 06 — Additional analyses added in the October 2026 revision. All results are written to ``results/reproduced/`` and consumed by scripts 04 and 05. Each block is a self-contained function. shapley.csv Shapley / Owen decomposition of R^2 across regressor blocks se_comparison.csv M5 key coefficients under alternative variance estimators, incl. a wild-cluster (Rademacher) bootstrap at the province level spline.csv piecewise-linear (spline) living-area elasticities vs. quadratic spatial_holdout.json FSA-level spatial holdout + repeated random holdouts moran_sensitivity.csv Moran's I for k in {5, 10, 20} and two subsamples gradient_metro.json urban gradient with metro-specific intercepts Usage: python scripts/06_revision_analyses.py """ import itertools import json import math import sys import warnings from pathlib import Path import numpy as np import pandas as pd import statsmodels.api as sm from sklearn.neighbors import NearestNeighbors sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) from wp9 import models, sample # noqa: E402 from wp9.config import METROS, REPRODUCED, SEED_MORAN, SEED_OOS, STRUCT, ensure_dirs # noqa: E402 warnings.filterwarnings("ignore") SEED_REV = 2026 # -------------------------------------------------------------------------- def _r2(y, X): """OLS R^2 with an intercept-only baseline (X may already contain a constant).""" if "const" not in X.columns: X = sm.add_constant(X) return float(sm.OLS(y, X).fit().rsquared) def _r2_absorb(y, X, fsa): """R^2 of a model that absorbs FSA fixed effects (full-dummy equivalent).""" from linearmodels.iv.absorbing import AbsorbingLS X = X.copy() if "const" not in X.columns: X = sm.add_constant(X) res = AbsorbingLS(y, X, absorb=fsa.to_frame().astype("category"), drop_absorbed=True).fit(cov_type="unadjusted") return float(res.rsquared) def shapley_decomposition(s: pd.DataFrame, label: str) -> list[dict]: """Shapley values of R^2 for the blocks {structure, type/ownership, province, FSA}. Province dummies are nested in the FSA effects, so any coalition that contains FSA has the same R^2 with or without the province block; the Shapley value of the province block then measures what provinces add when FSAs are *not* present, averaged over orderings. The Owen value groups province and FSA into a single 'location' block. """ y = s["ln_price"] blocks = { "S": s[STRUCT].astype(float), "T": pd.concat([pd.get_dummies(s[c], prefix=p, drop_first=True).astype(float) for c, p in (("btype_c", "bt"), ("own_c", "ow"), ("cat", "cat"))], axis=1), "P": pd.get_dummies(s["prov"], prefix="pv", drop_first=True).astype(float), } names = ["S", "T", "P", "F"] cache = {} def value(coal): key = tuple(sorted(coal)) if key in cache: return cache[key] cols = [blocks[b] for b in key if b in blocks] X = pd.concat(cols, axis=1) if cols else pd.DataFrame(index=s.index) if "F" in key: X = X.drop(columns=[c for c in X.columns if c.startswith("pv_")], errors="ignore") v = _r2_absorb(y, X, s["fsa_c"]) if len(X.columns) else _r2_absorb( y, pd.DataFrame({"const": np.ones(len(s))}, index=s.index), s["fsa_c"]) else: v = _r2(y, X) if len(X.columns) else 0.0 cache[key] = v return v n = len(names) shap = {b: 0.0 for b in names} for b in names: others = [o for o in names if o != b] for k in range(len(others) + 1): for coal in itertools.combinations(others, k): w = math.factorial(k) * math.factorial(n - k - 1) / math.factorial(n) shap[b] += w * (value(coal + (b,)) - value(coal)) total = value(tuple(names)) # Owen value: groups {S}, {T}, {L = P, F}; Shapley over the 3 groups. groups = {"S": ("S",), "T": ("T",), "L": ("P", "F")} gnames = list(groups) owen = {g: 0.0 for g in gnames} for g in gnames: others = [o for o in gnames if o != g] for k in range(len(others) + 1): for coal in itertools.combinations(others, k): w = math.factorial(k) * math.factorial(3 - k - 1) / math.factorial(3) with_g = sum((groups[x] for x in coal + (g,)), ()) without = sum((groups[x] for x in coal), ()) owen[g] += w * (value(with_g) - value(without)) seq = {"S": value(("S",)), "T": value(("S", "T")) - value(("S",)), "P": value(("S", "T", "P")) - value(("S", "T")), "F": value(("S", "T", "P", "F")) - value(("S", "T", "P"))} rows = [] for b, lab in (("S", "Structure"), ("T", "Dwelling type and ownership"), ("P", "Province"), ("F", "Neighbourhood (FSA)")): rows.append({"sample": label, "block": lab, "sequential": seq[b], "shapley": shap[b], "owen": owen["L"] if b == "F" else (np.nan if b == "P" else owen[b]), "total_r2": total, "n": len(s)}) return rows # -------------------------------------------------------------------------- def _within(df, cols, by): """Demean columns by group (one-way fixed-effect transformation).""" return df[cols] - df.groupby(by)[cols].transform("mean") def se_comparison(s: pd.DataFrame) -> pd.DataFrame: """M5 key coefficients under HC1, FSA-cluster, province-cluster and a province-level wild-cluster bootstrap (Rademacher, full enumeration).""" X = models.design(s).drop(columns="const") cols = X.columns.tolist() data = pd.concat([s[["ln_price", "fsa_c", "prov"]], X], axis=1) dd = _within(data, ["ln_price"] + cols, "fsa_c") y, Xd = dd["ln_price"].values, dd[cols].values keep = Xd.std(axis=0) > 1e-12 # drop columns absorbed by the FSA effects Xd, cols = Xd[:, keep], [c for c, k in zip(cols, keep) if k] n_fe = s["fsa_c"].nunique() ols = sm.OLS(y, Xd) base = ols.fit() dof_corr = (len(y) - 1) / (len(y) - Xd.shape[1] - n_fe) # absorbed-FE dof adjustment hc1 = ols.fit(cov_type="HC1") cl_fsa = ols.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(s["fsa_c"])[0]}) cl_prov = ols.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(s["prov"])[0]}) cl_2way = ols.fit(cov_type="cluster", cov_kwds={"groups": np.column_stack([pd.factorize(s["fsa_c"])[0], pd.factorize(s["prov"])[0]])}) rows = [] for key in ("ln_living", "bathrooms", "half_baths", "bedrooms", "ln_lot"): j = cols.index(key) rows.append({"variable": key, "coef": base.params[j], "se_hc1": hc1.bse[j] * math.sqrt(dof_corr), "se_fsa": cl_fsa.bse[j] * math.sqrt(dof_corr), "se_prov": cl_prov.bse[j], "se_twoway": cl_2way.bse[j] * math.sqrt(dof_corr)}) # Wild-cluster restricted bootstrap-t at the province level (G = 9 -> 512 draws) groups = pd.factorize(s["prov"])[0] G = groups.max() + 1 signs = np.array(list(itertools.product([-1.0, 1.0], repeat=G))) XtX_inv = np.linalg.inv(Xd.T @ Xd) for r in rows: j = cols.index(r["variable"]) t_obs = r["coef"] / r["se_prov"] mask = np.ones(Xd.shape[1], bool); mask[j] = False Xr = Xd[:, mask] br = np.linalg.lstsq(Xr, y, rcond=None)[0] yhat_r, e_r = Xr @ br, y - Xr @ br t_boot = np.empty(len(signs)) for b, w in enumerate(signs): yb = yhat_r + e_r * w[groups] bb = XtX_inv @ (Xd.T @ yb) ub = yb - Xd @ bb # cluster-robust variance of bb[j] by province meat = np.zeros((Xd.shape[1], Xd.shape[1])) for g in range(G): sg = Xd[groups == g].T @ ub[groups == g] meat += np.outer(sg, sg) V = XtX_inv @ meat @ XtX_inv * (G / (G - 1)) t_boot[b] = bb[j] / math.sqrt(V[j, j]) r["p_wild_prov"] = float(np.mean(np.abs(t_boot) >= abs(t_obs))) r["t_prov"] = float(t_obs) out = pd.DataFrame(rows) out.attrs["n_fsa_clusters"] = int(n_fe) out["n_fsa_clusters"] = n_fe out["n_prov_clusters"] = G return out # -------------------------------------------------------------------------- def spline_elasticity(s: pd.DataFrame) -> pd.DataFrame: """Piecewise-linear spline in ln(area) with knots at the size quintiles.""" knots = s["ln_living"].quantile([0.2, 0.4, 0.6, 0.8]).values extra = pd.DataFrame({f"hinge{k}": np.maximum(s["ln_living"] - kn, 0.0) for k, kn in enumerate(knots)}, index=s.index) res, cols = models.fit_absorbing(s, extra=extra) params = pd.Series(np.asarray(res.params).ravel(), index=res.params.index) cov = pd.DataFrame(np.asarray(res.cov), index=res.params.index, columns=res.params.index) names = ["ln_living"] + [f"hinge{k}" for k in range(len(knots))] edges = np.concatenate([[s["ln_living"].min()], knots, [s["ln_living"].max()]]) rows = [] for seg in range(len(knots) + 1): active = names[: seg + 1] g = np.array([1.0 if nm in active else 0.0 for nm in names]) el = float(g @ params[names].values) se = float(np.sqrt(g @ cov.loc[names, names].values @ g)) rows.append({"segment": seg + 1, "lo_m2": math.exp(edges[seg]), "hi_m2": math.exp(edges[seg + 1]), "elasticity": el, "se": se, "r2_spline": float(res.rsquared)}) # quadratic and linear fits for comparison res_q, _ = models.fit_absorbing(s, extra=(s["ln_living"] ** 2).to_frame("ln_living2")) res_l, _ = models.fit_absorbing(s) out = pd.DataFrame(rows) out["r2_quadratic"] = float(res_q.rsquared) out["r2_linear"] = float(res_l.rsquared) return out # -------------------------------------------------------------------------- def _score(test, pred_log, smear): err = test["ln_price"].values - pred_log ss_tot = np.sum((test["ln_price"].values - test["ln_price"].mean()) ** 2) ape = np.abs(np.exp(pred_log) * smear - test["price_cad"].values) / test["price_cad"].values return {"n_test": int(len(test)), "oos_r2": float(1 - np.sum(err ** 2) / ss_tot), "rmse_log": float(np.sqrt(np.mean(err ** 2))), "median_ape": float(np.median(ape) * 100), "mean_ape": float(np.mean(ape) * 100), # mean APE excluding the top 1% of errors: a handful of listings with # implausible parking/storey counts (e.g. 3,620 parking spaces) produce # astronomically large predictions in some splits "mean_ape_trim": float(np.mean(np.sort(ape)[: int(0.99 * len(ape))]) * 100), "within10": float(np.mean(ape <= 0.10) * 100), "within20": float(np.mean(ape <= 0.20) * 100)} def spatial_holdout(s: pd.DataFrame) -> dict: """Hold out 20% of FSAs entirely; predict their listings with (a) the province-mean fixed effect and (b) an inverse-distance average of the five nearest training FSAs' effects. Also: five repeated random 80/20 splits.""" out = {} rng = np.random.default_rng(SEED_REV) real = s[~s["fsa_c"].str.endswith("_other")] fsas = real["fsa_c"].unique() reps = [] for rep in range(5): held = rng.choice(fsas, size=int(0.2 * len(fsas)), replace=False) train = real[~real["fsa_c"].isin(held)] test = real[real["fsa_c"].isin(held)] res, cols = models.fit_absorbing(train) fe = models.fsa_fixed_effects(train, res, cols) beta = pd.Series(np.asarray(res.params).ravel(), index=res.params.index) xb_test = models.design(test).reindex(columns=beta.index, fill_value=0.0).values @ beta.values pred_train = models.predict_with_fe(train, res, cols, fe) smear = models.duan_smearing(train["ln_price"].values - pred_train) # (a) province mean of the training FSA effects fe_prov = (fe.to_frame("fe").join(train.groupby("fsa_c")["prov"].first()) .groupby("prov")["fe"].mean()) pred_a = xb_test + test["prov"].map(fe_prov).fillna(fe.mean()).values # (b) inverse-distance weighted mean of the 5 nearest training FSAs cent = train.groupby("fsa_c")[["lat", "lon"]].mean() test_cent = test.groupby("fsa_c")[["lat", "lon"]].mean() nn = NearestNeighbors(n_neighbors=5).fit(np.radians(cent.values)) dist, idx = nn.kneighbors(np.radians(test_cent.values)) dist_km = 6371.0 * dist w = 1.0 / np.maximum(dist_km, 1.0) fe_train = fe.reindex(cent.index).values fe_nn = pd.Series((fe_train[idx] * w).sum(1) / w.sum(1), index=test_cent.index) pred_b = xb_test + test["fsa_c"].map(fe_nn).values reps.append({"rep": rep, "n_fsa_held": int(len(held)), "prov_mean": _score(test, pred_a, smear), "nn_idw": _score(test, pred_b, smear)}) out["fsa_holdout"] = reps for key in ("prov_mean", "nn_idw"): out[f"fsa_holdout_mean_{key}"] = {m: float(np.mean([r[key][m] for r in reps])) for m in reps[0][key]} out[f"fsa_holdout_sd_{key}"] = {m: float(np.std([r[key][m] for r in reps], ddof=1)) for m in reps[0][key]} # repeated random 80/20 splits on common support (5 seeds) rand = [] for seed in (SEED_OOS, 1, 2, 3, 4): rs = np.random.RandomState(seed) idx = rs.permutation(len(s)); cut = int(0.8 * len(s)) train, test = s.iloc[idx[:cut]], s.iloc[idx[cut:]] common = set(train["fsa_c"]).intersection(test["fsa_c"]) train, test = train[train["fsa_c"].isin(common)], test[test["fsa_c"].isin(common)] res, cols = models.fit_absorbing(train) fe = models.fsa_fixed_effects(train, res, cols) smear = models.duan_smearing(train["ln_price"].values - models.predict_with_fe(train, res, cols, fe)) rand.append(_score(test, models.predict_with_fe(test, res, cols, fe), smear)) out["random_repeats"] = rand out["random_mean"] = {m: float(np.mean([r[m] for r in rand])) for m in rand[0]} out["random_sd"] = {m: float(np.std([r[m] for r in rand], ddof=1)) for m in rand[0]} return out # -------------------------------------------------------------------------- def moran_sensitivity(s: pd.DataFrame) -> pd.DataFrame: """Moran's I of structure-only and grand-model residuals for several k and subsamples.""" grand = pd.read_parquet(REPRODUCED / "grand_model.parquet") X = models.design(s) struct_resid = s["ln_price"] - sm.OLS(s["ln_price"], X).fit().predict(X) rows = [] for seed in (SEED_MORAN, SEED_REV): sub = s.sample(15000, random_state=seed) coords = np.radians(sub[["lat", "lon"]].values) for k in (5, 10, 20): nn = NearestNeighbors(n_neighbors=k + 1, metric="haversine").fit(coords) _, idx = nn.kneighbors(coords); idx = idx[:, 1:] for label, vals in (("structure_only", struct_resid.loc[sub.index].values), ("grand_model", grand.loc[sub.index, "resid_grand"].values)): z = vals - vals.mean(); lag = z[idx].mean(axis=1) stat = np.sum(z * lag) / np.sum(z ** 2) rng = np.random.default_rng(seed + k) perms = np.array([np.sum((zp := rng.permutation(z)) * zp[idx].mean(axis=1)) / np.sum(zp ** 2) for _ in range(199)]) rows.append({"seed": seed, "k": k, "model": label, "I": float(stat), "z": float((stat - perms.mean()) / perms.std()), "p_perm": float((np.sum(np.abs(perms) >= abs(stat)) + 1) / 200)}) return pd.DataFrame(rows) # -------------------------------------------------------------------------- def gradient_with_metro_fe(s: pd.DataFrame) -> dict: """Location premium on ln(distance) with nearest-metro intercepts, and a quadratic form.""" X = models.design(s) prem = s["ln_price"] - sm.OLS(s["ln_price"], X).fit().predict(X) lat, lon = s["lat"].values, s["lon"].values rad = np.pi / 180 dists = {} for name, (mlat, mlon) in METROS.items(): a = (np.sin((mlat - lat) * rad / 2) ** 2 + np.cos(lat * rad) * np.cos(mlat * rad) * np.sin((mlon - lon) * rad / 2) ** 2) dists[name] = 2 * 6371.0 * np.arcsin(np.sqrt(a)) D = pd.DataFrame(dists, index=s.index) dist = D.min(axis=1); metro = D.idxmin(axis=1) keep = dist > 0.5 ld = np.log(dist[keep]) base = sm.OLS(prem[keep], sm.add_constant(ld)).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(s["fsa_c"][keep])[0]}) Xm = pd.concat([ld.rename("ln_dist"), pd.get_dummies(metro[keep], drop_first=True).astype(float)], axis=1) mfe = sm.OLS(prem[keep], sm.add_constant(Xm)).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(s["fsa_c"][keep])[0]}) Xq = pd.DataFrame({"ln_dist": ld, "ln_dist2": ld ** 2}) quad = sm.OLS(prem[keep], sm.add_constant(Xq)).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(s["fsa_c"][keep])[0]}) # metro-specific slopes slopes = {} for m, g in pd.DataFrame({"prem": prem[keep], "ld": ld, "metro": metro[keep]}).groupby("metro"): f = sm.OLS(g["prem"], sm.add_constant(g["ld"])).fit() slopes[m] = {"beta": float(f.params.iloc[1]), "n": int(len(g))} return {"pooled": {"beta": float(base.params.iloc[1]), "se_fsa": float(base.bse.iloc[1]), "r2": float(base.rsquared)}, "metro_fe": {"beta": float(mfe.params["ln_dist"]), "se_fsa": float(mfe.bse["ln_dist"]), "r2": float(mfe.rsquared)}, "quadratic": {"b1": float(quad.params["ln_dist"]), "b2": float(quad.params["ln_dist2"]), "se1": float(quad.bse["ln_dist"]), "se2": float(quad.bse["ln_dist2"]), "r2": float(quad.rsquared)}, "metro_slopes": slopes, "share_within_40km": float((dist <= 40).mean()), "share_beyond_160km": float((dist > 160).mean())} # -------------------------------------------------------------------------- def main() -> None: ensure_dirs() s = sample.load_sample() rows = shapley_decomposition(s, "All residential") + \ shapley_decomposition(s[s["cat"] == "house"], "Houses") pd.DataFrame(rows).to_csv(REPRODUCED / "shapley.csv", index=False) print(pd.DataFrame(rows).round(3).to_string()) sec = se_comparison(s) sec.to_csv(REPRODUCED / "se_comparison.csv", index=False) print(sec.round(4).to_string()) spl = spline_elasticity(s) spl.to_csv(REPRODUCED / "spline.csv", index=False) print(spl.round(3).to_string()) sh = spatial_holdout(s) json.dump(sh, open(REPRODUCED / "spatial_holdout.json", "w"), indent=2) 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"))}) ms = moran_sensitivity(s) ms.to_csv(REPRODUCED / "moran_sensitivity.csv", index=False) print(ms.round(3).to_string()) gm = gradient_with_metro_fe(s) json.dump(gm, open(REPRODUCED / "gradient_metro.json", "w"), indent=2) print("gradient:", gm["pooled"], gm["metro_fe"], gm["quadratic"]) print("ALL REVISION ANALYSES DONE") if __name__ == "__main__": main()