spb/wp9_uqo Public

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    return pd.DataFrame(rows)909192# --------------------------------------------------------------------------93def heterogeneity(s: pd.DataFrame) -> pd.DataFrame:94    """Within-FSA living-area elasticity, estimated province by province."""95    rows = []96    for prov, g in s.groupby("prov"):97        if len(g) < 400:98            continue99        res, cols = models.fit_absorbing(g)100        tab = models.coef_table(res, cols)101        rows.append({"prov": prov, "n": len(g),102                     "elast": tab.loc["ln_living", "coef"], "se_el": tab.loc["ln_living", "se"],103                     "bath": tab.loc["bathrooms", "coef"], "se_bath": tab.loc["bathrooms", "se"]})104    return pd.DataFrame(rows).sort_values("elast")105106107# --------------------------------------------------------------------------108def quantile_regressions(s: pd.DataFrame) -> pd.DataFrame:109    """Quantile hedonic estimates (house subsample, province FE)."""110    houses = s[s["cat"] == "house"]111    prov_d = pd.get_dummies(houses["prov"], prefix="p", drop_first=True).astype(float)112    X = sm.add_constant(pd.concat([houses[STRUCT].astype(float), prov_d], axis=1))113    y = houses["ln_price"]114    rows = []115    for tau in (0.1, 0.25, 0.5, 0.75, 0.9):116        fit = QuantReg(y, X).fit(q=tau, max_iter=2000)117        rows.append({"tau": tau,118                     "ln_living": fit.params["ln_living"], "se_living": fit.bse["ln_living"],119                     "bathrooms": fit.params["bathrooms"], "se_bath": fit.bse["bathrooms"],120                     "ln_lot": fit.params["ln_lot"], "se_lot": fit.bse["ln_lot"]})121    ols = sm.OLS(y, X).fit(cov_type="HC1")122    rows.append({"tau": np.nan,123                 "ln_living": ols.params["ln_living"], "se_living": ols.bse["ln_living"],124                 "bathrooms": ols.params["bathrooms"], "se_bath": ols.bse["bathrooms"],125                 "ln_lot": ols.params["ln_lot"], "se_lot": ols.bse["ln_lot"]})126    return pd.DataFrame(rows)127128129# --------------------------------------------------------------------------130def nonlinearity(s: pd.DataFrame) -> tuple[dict, pd.DataFrame]:131    """Quadratic-in-log-area grand model and the implied marginal elasticity."""132    extra = (s["ln_living"] ** 2).to_frame("ln_living2")133    res, cols = models.fit_absorbing(s, extra=extra)134    params = pd.Series(np.asarray(res.params).ravel(), index=cols)135    cov = pd.DataFrame(np.asarray(res.cov), index=cols, columns=cols)136    b1, b2 = params["ln_living"], params["ln_living2"]137    grid = np.linspace(s["ln_living"].quantile(0.02), s["ln_living"].quantile(0.98), 60)138    marginal = b1 + 2 * b2 * grid139    v11 = cov.loc["ln_living", "ln_living"]140    v22 = cov.loc["ln_living2", "ln_living2"]141    v12 = cov.loc["ln_living", "ln_living2"]142    se = np.sqrt(v11 + 4 * grid ** 2 * v22 + 4 * grid * v12)143    band = pd.DataFrame({"ln_area": grid, "elasticity": marginal, "se": se})144    return {"b1": float(b1), "b2": float(b2), "r2": float(res.rsquared)}, band145146147# --------------------------------------------------------------------------148def haversine_km(lat1, lon1, lat2, lon2):149    """Great-circle distance in kilometres."""150    rad = np.pi / 180151    a = (np.sin((lat2 - lat1) * rad / 2) ** 2152         + np.cos(lat1 * rad) * np.cos(lat2 * rad) * np.sin((lon2 - lon1) * rad / 2) ** 2)153    return 2 * 6371.0 * np.arcsin(np.sqrt(a))154155156def distance_to_metro(s: pd.DataFrame) -> np.ndarray:157    """Distance from each listing to the nearest of the nine major metros."""158    dist = np.full(len(s), np.inf)159    for lat, lon in METROS.values():160        dist = np.minimum(dist, haversine_km(s["lat"].values, s["lon"].values, lat, lon))161    return dist162163164def urban_gradient(s: pd.DataFrame) -> tuple[dict, pd.DataFrame, pd.Series]:165    """Location premium (residual from a structure-only model) vs. metro distance."""166    X = models.design(s)167    structural = sm.OLS(s["ln_price"], X).fit()168    loc_premium = s["ln_price"] - structural.predict(X)169    dist = distance_to_metro(s)170171    bins = [0, 5, 10, 20, 40, 80, 160, 320, 2000]172    banded = (pd.DataFrame({"dist": dist, "prem": loc_premium})173              .assign(band=lambda x: pd.cut(x["dist"], bins))174              .groupby("band", observed=True)175              .agg(x=("dist", "median"), prem=("prem", "mean"), n=("prem", "size"))176              .reset_index(drop=True))177178    keep = dist > 0.5179    fit = sm.OLS(loc_premium[keep],180                 sm.add_constant(np.log(dist[keep]))).fit(cov_type="HC1")181    stats = {"beta_logdist": float(fit.params.iloc[1]), "se": float(fit.bse.iloc[1]),182             "r2": float(fit.rsquared)}183    return stats, banded, loc_premium184185186# --------------------------------------------------------------------------187def morans_i(coords: np.ndarray, values: np.ndarray, k: int = 10,188             permutations: int = 199) -> tuple[float, float, float]:189    """Moran's I with row-standardised kNN weights and a permutation test."""190    nn = NearestNeighbors(n_neighbors=k + 1).fit(coords)191    _, idx = nn.kneighbors(coords)192    idx = idx[:, 1:]193    z = values - values.mean()194    lag = z[idx].mean(axis=1)195    stat = np.sum(z * lag) / np.sum(z ** 2)196    perms = np.empty(permutations)197    for b in range(permutations):198        zp = np.random.permutation(z)199        perms[b] = np.sum(zp * zp[idx].mean(axis=1)) / np.sum(zp ** 2)200    zscore = (stat - perms.mean()) / perms.std()201    pvalue = (np.sum(np.abs(perms) >= abs(stat)) + 1) / (permutations + 1)202    return float(stat), float(zscore), float(pvalue)203204205def spatial_diagnostics(s: pd.DataFrame, loc_premium: pd.Series) -> dict:206    """Moran's I on structure-only vs. grand-model residuals (15k sample)."""207    grand = pd.read_parquet(REPRODUCED / "grand_model.parquet")208    subsample = s.sample(15000, random_state=SEED_MORAN)209    coords = subsample[["lat", "lon"]].values210    resid_struct = loc_premium.loc[subsample.index].values211    resid_grand = grand.loc[subsample.index, "resid_grand"].values212213    np.random.seed(SEED_EXT2)214    i_s, z_s, p_s = morans_i(coords, resid_struct)215    i_g, z_g, p_g = morans_i(coords, resid_grand)216    np.save(REPRODUCED / "moran_coords.npy", coords)217    np.save(REPRODUCED / "moran_resid_struct.npy", resid_struct)218    np.save(REPRODUCED / "moran_resid_grand.npy", resid_grand)219    return {"struct_I": i_s, "struct_z": z_s, "struct_p": p_s,220            "grand_I": i_g, "grand_z": z_g, "grand_p": p_g}221222223# --------------------------------------------------------------------------224def leave_one_province_out(s: pd.DataFrame) -> pd.DataFrame:225    """Estimate the structural model out-of-province; within-province R^2."""226    X = models.design(s)227    rows = []228    for prov in s["prov"].value_counts().index:229        held = s["prov"] == prov230        if held.sum() < 800:231            continue232        fit = sm.OLS(s["ln_price"].values[~held], X.values[~held]).fit()233        err = s["ln_price"].values[held] - X.values[held] @ fit.params234        err = err - err.mean()  # province-specific intercept allowed235        y_held = s["ln_price"].values[held]236        r2 = 1 - np.sum(err ** 2) / np.sum((y_held - y_held.mean()) ** 2)237        rows.append({"prov": prov, "n": int(held.sum()), "r2_within": float(r2)})238    return pd.DataFrame(rows).sort_values("r2_within", ascending=False)239240241# --------------------------------------------------------------------------242def main() -> None:243    ensure_dirs()244    s = sample.load_sample()245    results = {}246247    oos = out_of_sample(s)248    json.dump(oos, open(REPRODUCED / "oos.json", "w"), indent=2)249    print("OOS:", {k: round(v, 3) for k, v in oos.items()})250251    robustness(s).to_csv(REPRODUCED / "robustness.csv", index=False)252    print("robustness done")253254    heterogeneity(s).to_csv(REPRODUCED / "heterogeneity.csv", index=False)255    print("heterogeneity done")256257    quantile_regressions(s).to_csv(REPRODUCED / "quantile.csv", index=False)258    print("quantile done")259260    results["nonlin"], band = nonlinearity(s)261    band.to_csv(REPRODUCED / "nonlinear_band.csv", index=False)262    print("nonlinearity:", results["nonlin"])263264    results["gradient"], banded, loc_premium = urban_gradient(s)265    banded.to_csv(REPRODUCED / "gradient_bins.csv", index=False)266    print("gradient:", results["gradient"])267268    results["moran"] = spatial_diagnostics(s, loc_premium)269    print("moran:", {k: round(v, 3) for k, v in results["moran"].items()})270271    lopo = leave_one_province_out(s)272    lopo.to_csv(REPRODUCED / "lopo.csv", index=False)273    results["lopo_mean"] = float(lopo["r2_within"].mean())274    print("LOPO mean within-province R2:", round(results["lopo_mean"], 3))275276    json.dump(results, open(REPRODUCED / "ext2.json", "w"), indent=2)277    print("ALL EXTENDED ANALYSES DONE")278279280if __name__ == "__main__":281    main()282