spb/wp9_uqo

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 03 — Extended analyses.45Out-of-sample validation, robustness across samples, province heterogeneity,6quantile regressions, nonlinearity in floor space, the urban price gradient,7Moran's I spatial diagnostics, leave-one-province-out cross-validation.89Writes to ``results/reproduced/``:10  oos.json, robustness.csv, heterogeneity.csv, quantile.csv, lopo.csv,11  ext2.json, nonlinear_band.csv, gradient_bins.csv1213Usage:  python scripts/03_estimate_extended.py14"""15import json16import sys17import warnings18from pathlib import Path1920import numpy as np21import pandas as pd22import statsmodels.api as sm23from sklearn.neighbors import NearestNeighbors24from statsmodels.regression.quantile_regression import QuantReg2526sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2728from wp9 import models, sample  # noqa: E40229from wp9.config import (METROS, REPRODUCED, SEED_EXT2, SEED_MORAN,  # noqa: E40230                        SEED_OOS, STRUCT, ensure_dirs)3132warnings.filterwarnings("ignore")333435# --------------------------------------------------------------------------36def out_of_sample(s: pd.DataFrame) -> dict:37    """80/20 split validation of the grand model, on common-FSA support."""38    np.random.seed(SEED_OOS)  # legacy global seed, kept identical to the original39    idx = np.arange(len(s))40    np.random.shuffle(idx)41    cut = int(0.8 * len(s))42    train, test = s.iloc[idx[:cut]].copy(), s.iloc[idx[cut:]].copy()43    common = set(train["fsa_c"]).intersection(test["fsa_c"])44    train = train[train["fsa_c"].isin(common)].copy()45    test = test[test["fsa_c"].isin(common)].copy()4647    res, cols = models.fit_absorbing(train)48    fe = models.fsa_fixed_effects(train, res, cols)49    pred_test = models.predict_with_fe(test, res, cols, fe)50    err = test["ln_price"].values - pred_test51    ss_res = np.sum(err ** 2)52    ss_tot = np.sum((test["ln_price"].values - test["ln_price"].mean()) ** 2)5354    pred_train = models.predict_with_fe(train, res, cols, fe)55    smear = models.duan_smearing(train["ln_price"].values - pred_train)56    ape = np.abs(np.exp(pred_test) * smear - test["price_cad"].values) / test["price_cad"].values57    return {58        "n_train": int(len(train)), "n_test": int(len(test)),59        "oos_r2": float(1 - ss_res / ss_tot),60        "rmse_log": float(np.sqrt(np.mean(err ** 2))),61        "mae_log": float(np.mean(np.abs(err))),62        "median_ape": float(np.median(ape) * 100),63        "mean_ape": float(np.mean(ape) * 100),64        "within10": float(np.mean(ape <= 0.10) * 100),65        "within20": float(np.mean(ape <= 0.20) * 100),66    }676869# --------------------------------------------------------------------------70def robustness(s: pd.DataFrame) -> pd.DataFrame:71    """Key implicit prices of the grand model across alternative samples."""72    def key_coefs(data, label):73        res, cols = models.fit_absorbing(data)74        tab = models.coef_table(res, cols)75        return {"label": label, "n": int(res.nobs), "r2": float(res.rsquared),76                "ln_living": tab.loc["ln_living", "coef"], "se_living": tab.loc["ln_living", "se"],77                "bathrooms": tab.loc["bathrooms", "coef"], "se_bath": tab.loc["bathrooms", "se"],78                "ln_lot": tab.loc["ln_lot", "coef"], "se_lot": tab.loc["ln_lot", "se"]}7980    rows = [key_coefs(s, "Baseline (all residential)"),81            key_coefs(s[s["cat"] == "house"], "Houses only"),82            key_coefs(s[s["cat"] == "condo"], "Condominiums only")]83    lo, hi = s["price_cad"].quantile([0.005, 0.995])84    rows.append(key_coefs(s[s["price_cad"].between(lo, hi)], "Price trimmed 0.5/99.5\\%"))85    counts = s["fsa_c"].value_counts()86    rows.append(key_coefs(s[s["fsa_c"].isin(counts[counts >= 50].index)],87                          "FSAs with $\\geq$50 listings"))88    rows.append(key_coefs(s[s["prov"].isin(["ON", "QC", "BC"])], "ON/QC/BC only"))89    # Revision 2026-10: two further cuts. (i) Storeys and parking counts are90    # noisy (max 593 storeys, 3,620 parking spaces); cap both at their 99th91    # percentile. (ii) Drop the province-level residual categories that pool92    # FSAs with fewer than 25 listings, so every absorbed effect is a true FSA.93    capped = s.copy()94    for col in ("stories_n", "parking_n"):95        capped[col] = capped[col].clip(upper=capped[col].quantile(0.99))96    rows.append(key_coefs(capped, "Storeys/parking capped at P99"))97    rows.append(key_coefs(s[~s["fsa_c"].str.endswith("_other")], "Excl. pooled residual FSAs"))98    return pd.DataFrame(rows)99100101# --------------------------------------------------------------------------102def heterogeneity(s: pd.DataFrame) -> pd.DataFrame:103    """Within-FSA living-area elasticity, estimated province by province."""104    rows = []105    for prov, g in s.groupby("prov"):106        if len(g) < 400:107            continue108        res, cols = models.fit_absorbing(g)109        tab = models.coef_table(res, cols)110        rows.append({"prov": prov, "n": len(g),111                     "elast": tab.loc["ln_living", "coef"], "se_el": tab.loc["ln_living", "se"],112                     "bath": tab.loc["bathrooms", "coef"], "se_bath": tab.loc["bathrooms", "se"]})113    return pd.DataFrame(rows).sort_values("elast")114115116# --------------------------------------------------------------------------117def quantile_regressions(s: pd.DataFrame) -> pd.DataFrame:118    """Quantile hedonic estimates (house subsample, province FE)."""119    houses = s[s["cat"] == "house"]120    prov_d = pd.get_dummies(houses["prov"], prefix="p", drop_first=True).astype(float)121    X = sm.add_constant(pd.concat([houses[STRUCT].astype(float), prov_d], axis=1))122    y = houses["ln_price"]123    rows = []124    for tau in (0.1, 0.25, 0.5, 0.75, 0.9):125        # statsmodels default: Koenker--Bassett kernel (Epanechnikov) sandwich126        # variance with Hall--Sheather bandwidth, robust to heteroskedasticity.127        fit = QuantReg(y, X).fit(q=tau, vcov="robust", kernel="epa", bandwidth="hsheather",128                                 max_iter=2000)129        rows.append({"tau": tau,130                     "ln_living": fit.params["ln_living"], "se_living": fit.bse["ln_living"],131                     "bathrooms": fit.params["bathrooms"], "se_bath": fit.bse["bathrooms"],132                     "ln_lot": fit.params["ln_lot"], "se_lot": fit.bse["ln_lot"]})133    ols = sm.OLS(y, X).fit(cov_type="HC1")134    rows.append({"tau": np.nan,135                 "ln_living": ols.params["ln_living"], "se_living": ols.bse["ln_living"],136                 "bathrooms": ols.params["bathrooms"], "se_bath": ols.bse["bathrooms"],137                 "ln_lot": ols.params["ln_lot"], "se_lot": ols.bse["ln_lot"]})138    return pd.DataFrame(rows)139140141# --------------------------------------------------------------------------142def nonlinearity(s: pd.DataFrame) -> tuple[dict, pd.DataFrame]:143    """Quadratic-in-log-area grand model and the implied marginal elasticity."""144    extra = (s["ln_living"] ** 2).to_frame("ln_living2")145    res, cols = models.fit_absorbing(s, extra=extra)146    params = pd.Series(np.asarray(res.params).ravel(), index=cols)147    cov = pd.DataFrame(np.asarray(res.cov), index=cols, columns=cols)148    b1, b2 = params["ln_living"], params["ln_living2"]149    grid = np.linspace(s["ln_living"].quantile(0.02), s["ln_living"].quantile(0.98), 60)150    marginal = b1 + 2 * b2 * grid151    v11 = cov.loc["ln_living", "ln_living"]152    v22 = cov.loc["ln_living2", "ln_living2"]153    v12 = cov.loc["ln_living", "ln_living2"]154    se = np.sqrt(v11 + 4 * grid ** 2 * v22 + 4 * grid * v12)155    band = pd.DataFrame({"ln_area": grid, "elasticity": marginal, "se": se})156    return {"b1": float(b1), "b2": float(b2), "r2": float(res.rsquared)}, band157158159# --------------------------------------------------------------------------160def haversine_km(lat1, lon1, lat2, lon2):161    """Great-circle distance in kilometres."""162    rad = np.pi / 180163    a = (np.sin((lat2 - lat1) * rad / 2) ** 2164         + np.cos(lat1 * rad) * np.cos(lat2 * rad) * np.sin((lon2 - lon1) * rad / 2) ** 2)165    return 2 * 6371.0 * np.arcsin(np.sqrt(a))166167168def distance_to_metro(s: pd.DataFrame) -> np.ndarray:169    """Distance from each listing to the nearest of the nine major metros."""170    dist = np.full(len(s), np.inf)171    for lat, lon in METROS.values():172        dist = np.minimum(dist, haversine_km(s["lat"].values, s["lon"].values, lat, lon))173    return dist174175176def urban_gradient(s: pd.DataFrame) -> tuple[dict, pd.DataFrame, pd.Series]:177    """Location premium (residual from a structure-only model) vs. metro distance."""178    X = models.design(s)179    structural = sm.OLS(s["ln_price"], X).fit()180    loc_premium = s["ln_price"] - structural.predict(X)181    dist = distance_to_metro(s)182183    bins = [0, 5, 10, 20, 40, 80, 160, 320, 2000]184    banded = (pd.DataFrame({"dist": dist, "prem": loc_premium})185              .assign(band=lambda x: pd.cut(x["dist"], bins))186              .groupby("band", observed=True)187              .agg(x=("dist", "median"), prem=("prem", "mean"), n=("prem", "size"))188              .reset_index(drop=True))189190    keep = dist > 0.5191    fit = sm.OLS(loc_premium[keep],192                 sm.add_constant(np.log(dist[keep]))).fit(cov_type="HC1")193    stats = {"beta_logdist": float(fit.params.iloc[1]), "se": float(fit.bse.iloc[1]),194             "r2": float(fit.rsquared)}195    return stats, banded, loc_premium196197198# --------------------------------------------------------------------------199def morans_i(coords: np.ndarray, values: np.ndarray, k: int = 10,200             permutations: int = 199) -> tuple[float, float, float]:201    """Moran's I with row-standardised kNN weights and a permutation test."""202    nn = NearestNeighbors(n_neighbors=k + 1).fit(coords)203    _, idx = nn.kneighbors(coords)204    idx = idx[:, 1:]205    z = values - values.mean()206    lag = z[idx].mean(axis=1)207    stat = np.sum(z * lag) / np.sum(z ** 2)208    perms = np.empty(permutations)209    for b in range(permutations):210        zp = np.random.permutation(z)211        perms[b] = np.sum(zp * zp[idx].mean(axis=1)) / np.sum(zp ** 2)212    zscore = (stat - perms.mean()) / perms.std()213    pvalue = (np.sum(np.abs(perms) >= abs(stat)) + 1) / (permutations + 1)214    return float(stat), float(zscore), float(pvalue)215216217def spatial_diagnostics(s: pd.DataFrame, loc_premium: pd.Series) -> dict:218    """Moran's I on structure-only vs. grand-model residuals (15k sample)."""219    grand = pd.read_parquet(REPRODUCED / "grand_model.parquet")220    subsample = s.sample(15000, random_state=SEED_MORAN)221    coords = subsample[["lat", "lon"]].values222    resid_struct = loc_premium.loc[subsample.index].values223    resid_grand = grand.loc[subsample.index, "resid_grand"].values224225    np.random.seed(SEED_EXT2)226    i_s, z_s, p_s = morans_i(coords, resid_struct)227    i_g, z_g, p_g = morans_i(coords, resid_grand)228    np.save(REPRODUCED / "moran_coords.npy", coords)229    np.save(REPRODUCED / "moran_resid_struct.npy", resid_struct)230    np.save(REPRODUCED / "moran_resid_grand.npy", resid_grand)231    return {"struct_I": i_s, "struct_z": z_s, "struct_p": p_s,232            "grand_I": i_g, "grand_z": z_g, "grand_p": p_g}233234235# --------------------------------------------------------------------------236def leave_one_province_out(s: pd.DataFrame) -> pd.DataFrame:237    """Estimate the structural model out-of-province; within-province R^2."""238    X = models.design(s)239    rows = []240    for prov in s["prov"].value_counts().index:241        held = s["prov"] == prov242        if held.sum() < 800:243            continue244        fit = sm.OLS(s["ln_price"].values[~held], X.values[~held]).fit()245        err = s["ln_price"].values[held] - X.values[held] @ fit.params246        err = err - err.mean()  # province-specific intercept allowed247        y_held = s["ln_price"].values[held]248        r2 = 1 - np.sum(err ** 2) / np.sum((y_held - y_held.mean()) ** 2)249        rows.append({"prov": prov, "n": int(held.sum()), "r2_within": float(r2)})250    return pd.DataFrame(rows).sort_values("r2_within", ascending=False)251252253# --------------------------------------------------------------------------254def main() -> None:255    ensure_dirs()256    s = sample.load_sample()257    results = {}258259    oos = out_of_sample(s)260    json.dump(oos, open(REPRODUCED / "oos.json", "w"), indent=2)261    print("OOS:", {k: round(v, 3) for k, v in oos.items()})262263    robustness(s).to_csv(REPRODUCED / "robustness.csv", index=False)264    print("robustness done")265266    heterogeneity(s).to_csv(REPRODUCED / "heterogeneity.csv", index=False)267    print("heterogeneity done")268269    quantile_regressions(s).to_csv(REPRODUCED / "quantile.csv", index=False)270    print("quantile done")271272    results["nonlin"], band = nonlinearity(s)273    band.to_csv(REPRODUCED / "nonlinear_band.csv", index=False)274    print("nonlinearity:", results["nonlin"])275276    results["gradient"], banded, loc_premium = urban_gradient(s)277    banded.to_csv(REPRODUCED / "gradient_bins.csv", index=False)278    print("gradient:", results["gradient"])279280    results["moran"] = spatial_diagnostics(s, loc_premium)281    print("moran:", {k: round(v, 3) for k, v in results["moran"].items()})282283    lopo = leave_one_province_out(s)284    lopo.to_csv(REPRODUCED / "lopo.csv", index=False)285    results["lopo_mean"] = float(lopo["r2_within"].mean())286    print("LOPO mean within-province R2:", round(results["lopo_mean"], 3))287288    json.dump(results, open(REPRODUCED / "ext2.json", "w"), indent=2)289    print("ALL EXTENDED ANALYSES DONE")290291292if __name__ == "__main__":293    main()294