#!/usr/bin/env python3 # Author: Simon-Pierre Boucher — contact@spboucher.ai """Step 03 — Extended analyses. Out-of-sample validation, robustness across samples, province heterogeneity, quantile regressions, nonlinearity in floor space, the urban price gradient, Moran's I spatial diagnostics, leave-one-province-out cross-validation. Writes to ``results/reproduced/``: oos.json, robustness.csv, heterogeneity.csv, quantile.csv, lopo.csv, ext2.json, nonlinear_band.csv, gradient_bins.csv Usage: python scripts/03_estimate_extended.py """ import json 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 from statsmodels.regression.quantile_regression import QuantReg 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_EXT2, SEED_MORAN, # noqa: E402 SEED_OOS, STRUCT, ensure_dirs) warnings.filterwarnings("ignore") # -------------------------------------------------------------------------- def out_of_sample(s: pd.DataFrame) -> dict: """80/20 split validation of the grand model, on common-FSA support.""" np.random.seed(SEED_OOS) # legacy global seed, kept identical to the original idx = np.arange(len(s)) np.random.shuffle(idx) cut = int(0.8 * len(s)) train, test = s.iloc[idx[:cut]].copy(), s.iloc[idx[cut:]].copy() common = set(train["fsa_c"]).intersection(test["fsa_c"]) train = train[train["fsa_c"].isin(common)].copy() test = test[test["fsa_c"].isin(common)].copy() res, cols = models.fit_absorbing(train) fe = models.fsa_fixed_effects(train, res, cols) pred_test = models.predict_with_fe(test, res, cols, fe) err = test["ln_price"].values - pred_test ss_res = np.sum(err ** 2) ss_tot = np.sum((test["ln_price"].values - test["ln_price"].mean()) ** 2) pred_train = models.predict_with_fe(train, res, cols, fe) smear = models.duan_smearing(train["ln_price"].values - pred_train) ape = np.abs(np.exp(pred_test) * smear - test["price_cad"].values) / test["price_cad"].values return { "n_train": int(len(train)), "n_test": int(len(test)), "oos_r2": float(1 - ss_res / ss_tot), "rmse_log": float(np.sqrt(np.mean(err ** 2))), "mae_log": float(np.mean(np.abs(err))), "median_ape": float(np.median(ape) * 100), "mean_ape": float(np.mean(ape) * 100), "within10": float(np.mean(ape <= 0.10) * 100), "within20": float(np.mean(ape <= 0.20) * 100), } # -------------------------------------------------------------------------- def robustness(s: pd.DataFrame) -> pd.DataFrame: """Key implicit prices of the grand model across alternative samples.""" def key_coefs(data, label): res, cols = models.fit_absorbing(data) tab = models.coef_table(res, cols) return {"label": label, "n": int(res.nobs), "r2": float(res.rsquared), "ln_living": tab.loc["ln_living", "coef"], "se_living": tab.loc["ln_living", "se"], "bathrooms": tab.loc["bathrooms", "coef"], "se_bath": tab.loc["bathrooms", "se"], "ln_lot": tab.loc["ln_lot", "coef"], "se_lot": tab.loc["ln_lot", "se"]} rows = [key_coefs(s, "Baseline (all residential)"), key_coefs(s[s["cat"] == "house"], "Houses only"), key_coefs(s[s["cat"] == "condo"], "Condominiums only")] lo, hi = s["price_cad"].quantile([0.005, 0.995]) rows.append(key_coefs(s[s["price_cad"].between(lo, hi)], "Price trimmed 0.5/99.5\\%")) counts = s["fsa_c"].value_counts() rows.append(key_coefs(s[s["fsa_c"].isin(counts[counts >= 50].index)], "FSAs with $\\geq$50 listings")) rows.append(key_coefs(s[s["prov"].isin(["ON", "QC", "BC"])], "ON/QC/BC only")) return pd.DataFrame(rows) # -------------------------------------------------------------------------- def heterogeneity(s: pd.DataFrame) -> pd.DataFrame: """Within-FSA living-area elasticity, estimated province by province.""" rows = [] for prov, g in s.groupby("prov"): if len(g) < 400: continue res, cols = models.fit_absorbing(g) tab = models.coef_table(res, cols) rows.append({"prov": prov, "n": len(g), "elast": tab.loc["ln_living", "coef"], "se_el": tab.loc["ln_living", "se"], "bath": tab.loc["bathrooms", "coef"], "se_bath": tab.loc["bathrooms", "se"]}) return pd.DataFrame(rows).sort_values("elast") # -------------------------------------------------------------------------- def quantile_regressions(s: pd.DataFrame) -> pd.DataFrame: """Quantile hedonic estimates (house subsample, province FE).""" houses = s[s["cat"] == "house"] prov_d = pd.get_dummies(houses["prov"], prefix="p", drop_first=True).astype(float) X = sm.add_constant(pd.concat([houses[STRUCT].astype(float), prov_d], axis=1)) y = houses["ln_price"] rows = [] for tau in (0.1, 0.25, 0.5, 0.75, 0.9): fit = QuantReg(y, X).fit(q=tau, max_iter=2000) rows.append({"tau": tau, "ln_living": fit.params["ln_living"], "se_living": fit.bse["ln_living"], "bathrooms": fit.params["bathrooms"], "se_bath": fit.bse["bathrooms"], "ln_lot": fit.params["ln_lot"], "se_lot": fit.bse["ln_lot"]}) ols = sm.OLS(y, X).fit(cov_type="HC1") rows.append({"tau": np.nan, "ln_living": ols.params["ln_living"], "se_living": ols.bse["ln_living"], "bathrooms": ols.params["bathrooms"], "se_bath": ols.bse["bathrooms"], "ln_lot": ols.params["ln_lot"], "se_lot": ols.bse["ln_lot"]}) return pd.DataFrame(rows) # -------------------------------------------------------------------------- def nonlinearity(s: pd.DataFrame) -> tuple[dict, pd.DataFrame]: """Quadratic-in-log-area grand model and the implied marginal elasticity.""" extra = (s["ln_living"] ** 2).to_frame("ln_living2") res, cols = models.fit_absorbing(s, extra=extra) params = pd.Series(np.asarray(res.params).ravel(), index=cols) cov = pd.DataFrame(np.asarray(res.cov), index=cols, columns=cols) b1, b2 = params["ln_living"], params["ln_living2"] grid = np.linspace(s["ln_living"].quantile(0.02), s["ln_living"].quantile(0.98), 60) marginal = b1 + 2 * b2 * grid v11 = cov.loc["ln_living", "ln_living"] v22 = cov.loc["ln_living2", "ln_living2"] v12 = cov.loc["ln_living", "ln_living2"] se = np.sqrt(v11 + 4 * grid ** 2 * v22 + 4 * grid * v12) band = pd.DataFrame({"ln_area": grid, "elasticity": marginal, "se": se}) return {"b1": float(b1), "b2": float(b2), "r2": float(res.rsquared)}, band # -------------------------------------------------------------------------- def haversine_km(lat1, lon1, lat2, lon2): """Great-circle distance in kilometres.""" rad = np.pi / 180 a = (np.sin((lat2 - lat1) * rad / 2) ** 2 + np.cos(lat1 * rad) * np.cos(lat2 * rad) * np.sin((lon2 - lon1) * rad / 2) ** 2) return 2 * 6371.0 * np.arcsin(np.sqrt(a)) def distance_to_metro(s: pd.DataFrame) -> np.ndarray: """Distance from each listing to the nearest of the nine major metros.""" dist = np.full(len(s), np.inf) for lat, lon in METROS.values(): dist = np.minimum(dist, haversine_km(s["lat"].values, s["lon"].values, lat, lon)) return dist def urban_gradient(s: pd.DataFrame) -> tuple[dict, pd.DataFrame, pd.Series]: """Location premium (residual from a structure-only model) vs. metro distance.""" X = models.design(s) structural = sm.OLS(s["ln_price"], X).fit() loc_premium = s["ln_price"] - structural.predict(X) dist = distance_to_metro(s) bins = [0, 5, 10, 20, 40, 80, 160, 320, 2000] banded = (pd.DataFrame({"dist": dist, "prem": loc_premium}) .assign(band=lambda x: pd.cut(x["dist"], bins)) .groupby("band", observed=True) .agg(x=("dist", "median"), prem=("prem", "mean"), n=("prem", "size")) .reset_index(drop=True)) keep = dist > 0.5 fit = sm.OLS(loc_premium[keep], sm.add_constant(np.log(dist[keep]))).fit(cov_type="HC1") stats = {"beta_logdist": float(fit.params.iloc[1]), "se": float(fit.bse.iloc[1]), "r2": float(fit.rsquared)} return stats, banded, loc_premium # -------------------------------------------------------------------------- def morans_i(coords: np.ndarray, values: np.ndarray, k: int = 10, permutations: int = 199) -> tuple[float, float, float]: """Moran's I with row-standardised kNN weights and a permutation test.""" nn = NearestNeighbors(n_neighbors=k + 1).fit(coords) _, idx = nn.kneighbors(coords) idx = idx[:, 1:] z = values - values.mean() lag = z[idx].mean(axis=1) stat = np.sum(z * lag) / np.sum(z ** 2) perms = np.empty(permutations) for b in range(permutations): zp = np.random.permutation(z) perms[b] = np.sum(zp * zp[idx].mean(axis=1)) / np.sum(zp ** 2) zscore = (stat - perms.mean()) / perms.std() pvalue = (np.sum(np.abs(perms) >= abs(stat)) + 1) / (permutations + 1) return float(stat), float(zscore), float(pvalue) def spatial_diagnostics(s: pd.DataFrame, loc_premium: pd.Series) -> dict: """Moran's I on structure-only vs. grand-model residuals (15k sample).""" grand = pd.read_parquet(REPRODUCED / "grand_model.parquet") subsample = s.sample(15000, random_state=SEED_MORAN) coords = subsample[["lat", "lon"]].values resid_struct = loc_premium.loc[subsample.index].values resid_grand = grand.loc[subsample.index, "resid_grand"].values np.random.seed(SEED_EXT2) i_s, z_s, p_s = morans_i(coords, resid_struct) i_g, z_g, p_g = morans_i(coords, resid_grand) np.save(REPRODUCED / "moran_coords.npy", coords) np.save(REPRODUCED / "moran_resid_struct.npy", resid_struct) np.save(REPRODUCED / "moran_resid_grand.npy", resid_grand) return {"struct_I": i_s, "struct_z": z_s, "struct_p": p_s, "grand_I": i_g, "grand_z": z_g, "grand_p": p_g} # -------------------------------------------------------------------------- def leave_one_province_out(s: pd.DataFrame) -> pd.DataFrame: """Estimate the structural model out-of-province; within-province R^2.""" X = models.design(s) rows = [] for prov in s["prov"].value_counts().index: held = s["prov"] == prov if held.sum() < 800: continue fit = sm.OLS(s["ln_price"].values[~held], X.values[~held]).fit() err = s["ln_price"].values[held] - X.values[held] @ fit.params err = err - err.mean() # province-specific intercept allowed y_held = s["ln_price"].values[held] r2 = 1 - np.sum(err ** 2) / np.sum((y_held - y_held.mean()) ** 2) rows.append({"prov": prov, "n": int(held.sum()), "r2_within": float(r2)}) return pd.DataFrame(rows).sort_values("r2_within", ascending=False) # -------------------------------------------------------------------------- def main() -> None: ensure_dirs() s = sample.load_sample() results = {} oos = out_of_sample(s) json.dump(oos, open(REPRODUCED / "oos.json", "w"), indent=2) print("OOS:", {k: round(v, 3) for k, v in oos.items()}) robustness(s).to_csv(REPRODUCED / "robustness.csv", index=False) print("robustness done") heterogeneity(s).to_csv(REPRODUCED / "heterogeneity.csv", index=False) print("heterogeneity done") quantile_regressions(s).to_csv(REPRODUCED / "quantile.csv", index=False) print("quantile done") results["nonlin"], band = nonlinearity(s) band.to_csv(REPRODUCED / "nonlinear_band.csv", index=False) print("nonlinearity:", results["nonlin"]) results["gradient"], banded, loc_premium = urban_gradient(s) banded.to_csv(REPRODUCED / "gradient_bins.csv", index=False) print("gradient:", results["gradient"]) results["moran"] = spatial_diagnostics(s, loc_premium) print("moran:", {k: round(v, 3) for k, v in results["moran"].items()}) lopo = leave_one_province_out(s) lopo.to_csv(REPRODUCED / "lopo.csv", index=False) results["lopo_mean"] = float(lopo["r2_within"].mean()) print("LOPO mean within-province R2:", round(results["lopo_mean"], 3)) json.dump(results, open(REPRODUCED / "ext2.json", "w"), indent=2) print("ALL EXTENDED ANALYSES DONE") if __name__ == "__main__": main()