spb/wp9_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 04 — Generate every figure used in the paper (17 PNG files).45Number-bearing summary figures (R^2 ladder, variance decomposition, OOS6accuracy, FSA premia, heterogeneity, quantile) are drawn from the results7tier selected with ``--results`` so the paper's published numbers are used by8default; distribution/scatter/map figures are drawn from the micro sample.910Usage:  python scripts/04_make_figures.py [--results reference|reproduced]11"""12import argparse13import json14import sys15from pathlib import Path1617import numpy as np18import pandas as pd1920sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2122from wp9 import models, sample  # noqa: E40223from wp9.config import FIGURES, REFERENCE, REPRODUCED, RESULTS, ensure_dirs  # noqa: E40224from wp9.plotstyle import (ACCENT, ACCENT2, ACCENT3, GREEN, GREY, LIGHT,  # noqa: E40225                           NEUTRAL, ORANGE, apply_style)2627apply_style()28import matplotlib.pyplot as plt  # noqa: E40229import statsmodels.api as sm  # noqa: E4023031PROVINCE_NAMES = {"ON": "Ontario", "QC": "Quebec", "BC": "British Columbia",32                  "AB": "Alberta", "SK": "Saskatchewan", "MB": "Manitoba",33                  "NS": "Nova Scotia", "NL": "Nfld. & Labrador", "NB": "New Brunswick"}343536def save(fig, name):37    fig.tight_layout()38    fig.savefig(FIGURES / name, bbox_inches="tight")39    plt.close(fig)40    print(" fig", name)414243# ---------------------------------------------------------------- micro-data44def fig_price_dist(s):45    fig, ax = plt.subplots(1, 2, figsize=(11, 4))46    shown = s[s["price_cad"] <= 3e6]["price_cad"]47    ax[0].hist(shown / 1e3, bins=80, color=ACCENT, alpha=0.9)48    ax[0].axvline(s["price_cad"].median() / 1e3, color=ACCENT2, ls="--", lw=1.2,49                  label=f"median \\${s['price_cad'].median():,.0f}")50    ax[0].set_xlabel("List price (thousand CAD, truncated at \\$3M)")51    ax[0].set_ylabel("Listings")52    ax[0].set_title("(a) Raw list price")53    ax[0].legend()54    ax[1].hist(s["ln_price"], bins=80, color=ACCENT, alpha=0.9)55    ax[1].set_xlabel("ln(price)")56    ax[1].set_title("(b) Log price (dependent variable)")57    save(fig, "fig_price_dist.png")585960def fig_province_ppm2(s):61    g = (s.groupby("prov")["ppm2"].agg(["median", "mean"])62         .sort_values("median"))63    fig, ax = plt.subplots(figsize=(7.5, 4.2))64    y = np.arange(len(g))65    ax.barh(y, g["median"], color=ACCENT, alpha=0.9, label="median")66    ax.barh(y, (g["mean"] - g["median"]).clip(lower=0), left=g["median"],67            color=GREY, alpha=0.55, label="mean$-$median gap")68    ax.set_yticks(y)69    ax.set_yticklabels([PROVINCE_NAMES.get(p, p) for p in g.index])70    ax.set_xlabel("Price per m$^2$ of living area (CAD)")71    ax.set_title("Price per square metre by province")72    ax.legend()73    save(fig, "fig_province_ppm2.png")747576def fig_size_gradient(s):77    bins = np.exp(np.linspace(np.log(45), np.log(470), 14))78    banded = (s.assign(band=pd.cut(s["living_m2"], bins))79              .groupby(["band", "cat"], observed=True)80              .agg(x=("living_m2", "median"), price=("price_cad", "median"),81                   n=("price_cad", "size"))82              .reset_index())83    fig, ax = plt.subplots(figsize=(7, 4.4))84    for cat, colour, label in (("house", ACCENT, "Houses"), ("condo", ACCENT2, "Condominiums")):85        g = banded[(banded["cat"] == cat) & (banded["n"] >= 30)]86        ax.plot(g["x"], g["price"] / 1e3, "o-", color=colour, lw=2, ms=5, label=label)87    ax.set_xscale("log")88    ax.set_yscale("log")89    ax.set_xlabel("Living area (m$^2$, log scale)")90    ax.set_ylabel("Median list price (thousand CAD, log scale)")91    ax.set_title("Median price by living-area bin and dwelling type")92    ax.legend()93    save(fig, "fig_size_gradient.png")949596def _canada_axes(ax, s):97    ax.set_xlim(-140, -50)98    ax.set_ylim(41, 61)99    ax.set_xlabel("Longitude")100    ax.set_ylabel("Latitude")101    for prov, g in s.groupby("prov"):102        if len(g) < 400:103            continue104        ax.annotate(prov, (g["lon"].median(), g["lat"].quantile(0.9) + 1.2),105                    fontsize=8, fontweight="bold", color="#333333", ha="center")106107108def fig_maps(s):109    colour = np.log(s["ppm2"])110    fig, ax = plt.subplots(figsize=(7, 5))111    sc = ax.scatter(s["lon"], s["lat"], c=colour, s=2, alpha=0.35, cmap="viridis")112    _canada_axes(ax, s)113    ax.set_title("Listing locations, coloured by log price per m$^2$")114    fig.colorbar(sc, ax=ax, shrink=0.75, label="ln(price per m$^2$)")115    save(fig, "fig_map.png")116117    g = (s.groupby("fsa")118         .agg(lat=("lat", "mean"), lon=("lon", "mean"),119              ppm2=("ppm2", "median"), n=("ppm2", "size")))120    g = g[g["n"] >= 25]121    fig, ax = plt.subplots(figsize=(7, 5))122    sc = ax.scatter(g["lon"], g["lat"], c=np.log(g["ppm2"]),123                    s=np.sqrt(g["n"]) * 1.8, alpha=0.75, cmap="viridis",124                    edgecolors="white", linewidths=0.2)125    _canada_axes(ax, s)126    ax.set_title("FSA neighbourhood medians (bubble area $\\propto\\sqrt{\\mathrm{listings}}$)")127    fig.colorbar(sc, ax=ax, shrink=0.75, label="ln(median price per m$^2$)")128    save(fig, "fig_fsa_map.png")129130131def fig_fit_resid():132    grand = pd.read_parquet(REPRODUCED / "grand_model.parquet")133    fig, ax = plt.subplots(figsize=(5.4, 5.2))134    ax.hexbin(grand["pred_grand"], grand["ln_price"], gridsize=90, cmap="Blues",135              mincnt=1, bins="log")136    lims = [grand["ln_price"].min(), grand["ln_price"].max()]137    ax.plot(lims, lims, color=ACCENT2, lw=1.4, ls="--")138    ax.set_xlabel("Predicted ln(price)")139    ax.set_ylabel("Actual ln(price)")140    ax.set_title("Grand model: predicted vs. actual")141    save(fig, "fig_fit.png")142143    resid = grand["resid_grand"]144    fig, ax = plt.subplots(1, 2, figsize=(9.6, 4.1))145    ax[0].hist(resid, bins=100, color=ACCENT, alpha=0.9, density=True)146    grid = np.linspace(resid.quantile(0.001), resid.quantile(0.999), 200)147    ax[0].plot(grid, np.exp(-0.5 * ((grid - resid.mean()) / resid.std()) ** 2)148               / (resid.std() * np.sqrt(2 * np.pi)), color=ACCENT2, lw=1.5,149               label="Normal density")150    ax[0].set_xlabel("Residual")151    ax[0].set_title("(a) Residual distribution")152    ax[0].legend()153    sm.qqplot(resid, line="45", fit=True, ax=ax[1], markerfacecolor=ACCENT,154              markeredgecolor=ACCENT, markersize=2, alpha=0.4)155    ax[1].set_title("(b) Normal Q--Q plot")156    save(fig, "fig_resid.png")157158159def fig_moran(res_dir):160    ext2 = json.load(open(res_dir / "ext2.json"))161    coords = np.load(REPRODUCED / "moran_coords.npy")162    panels = [("Structural-only residuals", np.load(REPRODUCED / "moran_resid_struct.npy"),163               ext2["moran"]["struct_I"]),164              ("Grand model (FSA FE) residuals", np.load(REPRODUCED / "moran_resid_grand.npy"),165               ext2["moran"]["grand_I"])]166    from sklearn.neighbors import NearestNeighbors167    fig, ax = plt.subplots(1, 2, figsize=(11, 4.4))168    for j, (title, resid, moran) in enumerate(panels):169        nn = NearestNeighbors(n_neighbors=11).fit(coords)170        _, idx = nn.kneighbors(coords)171        z = resid - resid.mean()172        lag = z[idx[:, 1:]].mean(axis=1)173        ax[j].scatter(z, lag, s=3, alpha=0.15, color=ACCENT)174        slope, intercept = np.polyfit(z, lag, 1)175        xs = np.linspace(z.min(), z.max(), 10)176        ax[j].plot(xs, slope * xs + intercept, color=ACCENT2, lw=1.8)177        ax[j].axhline(0, color="k", lw=0.5)178        ax[j].axvline(0, color="k", lw=0.5)179        ax[j].set_xlabel("Residual ($z$)")180        ax[j].set_ylabel("Spatial lag of residual")181        ax[j].set_title(f"{title}\nMoran's I = {moran:.3f}")182    save(fig, "fig_moran.png")183184185def fig_gradient(res_dir):186    banded = pd.read_csv(REPRODUCED / "gradient_bins.csv")187    fig, ax = plt.subplots(figsize=(7, 4.2))188    ax.plot(banded["x"], (np.exp(banded["prem"]) - 1) * 100, "o-", color=ACCENT, lw=2)189    ax.axhline(0, color="k", lw=0.7, ls=":")190    ax.set_xscale("symlog")191    ax.set_xlabel("Distance to nearest major metro (km, symlog)")192    ax.set_ylabel("Location premium vs structure-only (%)")193    ax.set_title("The urban price gradient: value falls with distance to metro")194    save(fig, "fig_gradient.png")195196197def fig_nonlinear(res_dir):198    band = pd.read_csv(REPRODUCED / "nonlinear_band.csv")199    baseline = pd.read_csv(res_dir / "robustness.csv").iloc[0]["ln_living"]200    fig, ax = plt.subplots(figsize=(7, 4.2))201    area = np.exp(band["ln_area"])202    ax.plot(area, band["elasticity"], color=ACCENT, lw=2)203    ax.fill_between(area, band["elasticity"] - 1.96 * band["se"],204                    band["elasticity"] + 1.96 * band["se"], color=ACCENT, alpha=0.18)205    ax.axhline(baseline, color=ACCENT2, ls="--", lw=1,206               label=f"linear-model elasticity ({baseline:.2f})")207    ax.set_xscale("log")208    ax.set_xlabel("Living area (m$^2$)")209    ax.set_ylabel("Marginal elasticity of price w.r.t. area")210    ax.set_title("Diminishing returns to floor space (quadratic spec., grand model)")211    ax.legend()212    save(fig, "fig_nonlinear.png")213214215# ----------------------------------------------------------- results-driven216def fig_r2(res_dir):217    fit = json.load(open(res_dir / "fit.json"))218    names = ["M1", "M2", "M3", "M4", "M5"]219    labels = ["M1\nStructural", "M2\n+Type/Own.", "M3\n+Province",220              "M4\nHouses+FSA", "M5\nGrand+FSA"]221    values = [fit[m]["r2"] for m in names]222    fig, ax = plt.subplots(figsize=(7.5, 4))223    bars = ax.bar(labels, values, color=[ACCENT, ACCENT, ACCENT, ACCENT3, ACCENT2], alpha=0.92)224    for bar, val in zip(bars, values):225        ax.text(bar.get_x() + bar.get_width() / 2, val + 0.012, f"{val:.3f}",226                ha="center", fontsize=9)227    ax.set_ylim(0, 0.9)228    ax.set_ylabel("$R^2$ (share of log-price variance explained)")229    ax.set_title("Explanatory power across the specification ladder")230    save(fig, "fig_r2.png")231232233def fig_decomp(res_dir):234    fit = json.load(open(res_dir / "fit.json"))235    blocks = [("Structure", fit["M1"]["r2"]),236              ("Dwelling type & ownership", fit["M2"]["r2"] - fit["M1"]["r2"]),237              ("Province", fit["M3"]["r2"] - fit["M2"]["r2"]),238              ("Neighbourhood (FSA)", fit["M5"]["r2"] - fit["M3"]["r2"]),239              ("Unexplained", 1 - fit["M5"]["r2"])]240    colours = [ACCENT, ACCENT3, GREEN, ACCENT2, NEUTRAL]241    fig, ax = plt.subplots(figsize=(8, 2.4))242    left = 0.0243    for (label, share), colour in zip(blocks, colours):244        ax.barh(0, share, left=left, color=colour, edgecolor="white")245        if share > 0.03:246            ax.text(left + share / 2, 0, f"{label}\n{share * 100:.1f}%",247                    ha="center", va="center", fontsize=8,248                    color="black" if colour == NEUTRAL else "white")249        left += share250    ax.set_xlim(0, 1)251    ax.set_ylim(-0.5, 0.5)252    ax.set_yticks([])253    ax.set_xlabel("Share of log-price variance")254    ax.set_title("Variance decomposition of Canadian house prices")255    save(fig, "fig_decomp.png")256257258def fig_forest(res_dir):259    coefs = pd.read_csv(res_dir / "coef_M3.csv", index_col=0)260    order = [("ln_living", "$\\ln$ living area"), ("bathrooms", "Full bathrooms"),261             ("half_baths", "Half bathrooms"), ("bedrooms", "Bedrooms"),262             ("parking_n", "Parking spaces"), ("stories_n", "Storeys"),263             ("ln_lot", "$\\ln(1+$lot m$^2)$")]264    rows = [(label, coefs.loc[key, "coef"], coefs.loc[key, "se"])265            for key, label in order if key in coefs.index]266    fig, ax = plt.subplots(figsize=(7, 4))267    y = np.arange(len(rows))[::-1]268    ax.errorbar([r[1] for r in rows], y, xerr=[1.96 * r[2] for r in rows],269                fmt="o", color=ACCENT, ecolor=GREY, capsize=3, ms=6)270    ax.axvline(0, color=ACCENT2, ls="--", lw=1)271    ax.set_yticks(y)272    ax.set_yticklabels([r[0] for r in rows])273    ax.set_xlabel("Coefficient on ln(price), 95% cluster-robust CI")274    ax.set_title("Structural implicit prices (houses, province-FE model M3)")275    save(fig, "fig_forest.png")276277278def fig_heterogeneity(res_dir):279    het = pd.read_csv(res_dir / "heterogeneity.csv").sort_values("elast")280    fig, ax = plt.subplots(figsize=(7, 4))281    y = np.arange(len(het))282    ax.errorbar(het["elast"], y, xerr=1.96 * het["se_el"], fmt="o",283                color=ACCENT, ecolor=GREY, capsize=3, ms=6)284    ax.axvline(het["elast"].mean(), color=ACCENT2, ls="--", lw=1,285               label="cross-province mean")286    ax.set_yticks(y)287    ax.set_yticklabels([f"{p} (n={n:,})" for p, n in zip(het["prov"], het["n"])])288    ax.set_xlabel("Living-area elasticity (within-FSA)")289    ax.set_title("Heterogeneity in the size elasticity of price across provinces")290    ax.legend()291    save(fig, "fig_heterogeneity.png")292293294def fig_premia(res_dir):295    premia = pd.read_csv(res_dir / "fsa_premia.csv")296    premia = premia.rename(columns={premia.columns[0]: "fsa"})297    shown = pd.concat([premia.nsmallest(12, "premium_pct"),298                       premia.nlargest(12, "premium_pct")]).sort_values("premium_pct")299    colours = [ACCENT2 if v < 0 else ACCENT for v in shown["premium_pct"]]300    fig, ax = plt.subplots(figsize=(7.5, 5))301    ax.barh(np.arange(len(shown)), shown["premium_pct"], color=colours, alpha=0.9)302    ax.set_yticks(np.arange(len(shown)))303    ax.set_yticklabels([f"{f} ({p})" for f, p in zip(shown["fsa"], shown["prov"])],304                       fontsize=8)305    ax.axvline(0, color="k", lw=0.8)306    ax.set_xlabel("Neighbourhood (FSA) price premium vs. national median (%), net of structure")307    ax.set_title("Highest- and lowest-valued neighbourhoods in Canada")308    save(fig, "fig_premia.png")309310311def fig_quantile(res_dir):312    q = pd.read_csv(res_dir / "quantile.csv")313    taus = q[q["tau"].notna()]314    ols = q[q["tau"].isna()].iloc[0]315    fig, ax = plt.subplots(1, 2, figsize=(11, 4.2))316    for j, (var, se_var, title) in enumerate(317            [("ln_living", "se_living", "$\\ln$ living-area elasticity"),318             ("bathrooms", "se_bath", "Full-bathroom premium")]):319        ax[j].errorbar(taus["tau"], taus[var], yerr=1.96 * taus[se_var],320                       fmt="o-", color=ACCENT, capsize=3, label="Quantile")321        ax[j].axhline(ols[var], color=ACCENT2, ls="--", label="OLS")322        ax[j].set_xlabel("Quantile of price ($\\tau$)")323        ax[j].set_title(title)324        ax[j].legend()325    save(fig, "fig_quantile.png")326327328def fig_oos(res_dir):329    oos = json.load(open(res_dir / "oos.json"))330    within10 = oos["within10"]331    mid = oos["within20"] - oos["within10"]332    beyond = 100 - oos["within20"]333    fig, ax = plt.subplots(figsize=(6, 4))334    ax.bar([5, 15, 30], [within10, mid, beyond], width=8,335           color=[ACCENT, ACCENT3, LIGHT], alpha=0.9)336    for x, v, label in ((5, within10, "≤10%"), (15, mid, "10–20%"), (30, beyond, ">20%")):337        ax.text(x, v + 1, f"{label}\n{v:.0f}%", ha="center", fontsize=9)338    ax.axvline(oos["median_ape"], color=ACCENT2, ls="--",339               label=f"median APE {oos['median_ape']:.1f}%")340    ax.set_xticks([5, 15, 30])341    ax.set_xticklabels(["≤10%", "10–20%", ">20%"])342    ax.set_xlabel("Absolute percentage error (out-of-sample)")343    ax.set_ylabel("Share of test listings (%)")344    ax.set_ylim(0, 52)345    ax.set_title(f"Out-of-sample valuation accuracy (OOS $R^2$={oos['oos_r2']:.3f})", pad=12)346    ax.legend(loc="upper right", fontsize=8)347    save(fig, "fig_oos.png")348349350def main() -> None:351    parser = argparse.ArgumentParser()352    parser.add_argument("--results", choices=["reference", "reproduced"],353                        default="reference",354                        help="results tier for number-bearing figures (default: reference)")355    args = parser.parse_args()356    res_dir = REFERENCE if args.results == "reference" else REPRODUCED357    ensure_dirs()358    s = sample.load_sample()359360    fig_price_dist(s)361    fig_province_ppm2(s)362    fig_size_gradient(s)363    fig_maps(s)364    fig_fit_resid()365    fig_r2(res_dir)366    fig_decomp(res_dir)367    fig_forest(res_dir)368    fig_heterogeneity(res_dir)369    fig_premia(res_dir)370    fig_quantile(res_dir)371    fig_oos(res_dir)372    fig_nonlinear(res_dir)373    fig_gradient(res_dir)374    fig_moran(res_dir)375    print("all figures written to", FIGURES)376377378if __name__ == "__main__":379    main()380