spb/wp3_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3#4"""Regenerate all 11 paper figures from the stored pickles.56Every figure is rebuilt from the data actually used in the paper7(analytical_sample.pkl, model_data.pkl, ml_results.pkl, qr_results.pkl,8shap_data.pkl) in a unified publication style. The script prints9verification statistics (Jarque-Bera, skewness, kurtosis, model R2, OLS10reference coefficients) so the regenerated content can be checked against11the numbers reported in the paper. The six original PNGs remain untouched12in the source archive (immo-wp3-spb-20260519/figures/).1314Usage:15    python scripts/03_make_figures.py               # all 11 figures16    python scripts/03_make_figures.py --figs fig2 fig417"""1819import argparse20import sys21from pathlib import Path2223import matplotlib2425matplotlib.use("Agg")26import matplotlib.pyplot as plt27import numpy as np28import pandas as pd29import scipy.stats as st3031sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))32from wp3 import config, data33from wp3.plotting import (ACCENT, CMAP_DIV, CMAP_SEQ, NEUTRAL, PRIMARY,34                          PRIMARY_LIGHT, apply_paper_style, label, panel_title)3536RNG = np.random.RandomState(config.RANDOM_STATE)3738FILE_NAMES = {39    "fig1": "fig1_price_distribution.png",40    "fig2": "fig2_ols_diagnostics.png",41    "fig3": "fig3_quantile_coefficients.png",42    "fig4": "fig4_model_comparison.png",43    "fig5": "fig5_shap_summary.png",44    "fig6": "fig6_shap_importance.png",45    "fig7": "fig7_geographic_prices.png",46    "fig8": "fig8_regional_prices.png",47    "fig9": "fig9_marginal_effects.png",48    "fig10": "fig10_shap_dependence.png",49    "fig11": "fig11_shap_interactions.png",50}515253def fig1_price_distribution(path):54    """Histograms of price (thousands) and log-price with median markers."""55    df = data.load_analytical_sample()56    price_k = df["price"].values / 1000.05758    fig, axes = plt.subplots(1, 2, figsize=(16, 6))59    ax = axes[0]60    ax.hist(price_k[price_k <= 3000], bins=100, color=PRIMARY, alpha=0.85,61            edgecolor="white", linewidth=0.2)62    med = np.median(price_k)63    ax.axvline(med, color=ACCENT, linestyle="--", linewidth=1.8,64               label=f"Median: \\${med:,.0f}K")65    ax.set_xlabel("Listing Price (\\$ thousands)")66    ax.set_ylabel("Frequency")67    ax.yaxis.set_major_formatter(plt.matplotlib.ticker.StrMethodFormatter("{x:,.0f}"))68    panel_title(ax, "(a) Price Distribution")69    ax.legend()7071    ax = axes[1]72    ax.hist(df["ln_price"], bins=100, color=PRIMARY, alpha=0.85,73            edgecolor="white", linewidth=0.2)74    med_ln = df["ln_price"].median()75    ax.axvline(med_ln, color=ACCENT, linestyle="--", linewidth=1.8,76               label=f"Median: {med_ln:.2f}")77    ax.set_xlabel("Log(Price)")78    ax.set_ylabel("Frequency")79    ax.yaxis.set_major_formatter(plt.matplotlib.ticker.StrMethodFormatter("{x:,.0f}"))80    panel_title(ax, "(b) Log-Price Distribution")81    ax.legend()8283    fig.tight_layout()84    fig.savefig(path)85    plt.close(fig)868788def fig2_ols_diagnostics(path):89    """OLS diagnostics: residuals vs fitted, Q-Q, distribution, scale-location."""90    md = data.load_model_data()91    resid = np.asarray(md["residuals"], dtype=float)92    fitted = np.asarray(md["y_clean"], dtype=float) - resid9394    jb = st.jarque_bera(resid)[0]95    print(f"  [verify] skew={st.skew(resid):.2f} (paper: 0.18)  "96          f"kurtosis={st.kurtosis(resid, fisher=False):.2f} (paper: 5.55)  "97          f"JB={jb:,.0f} (paper: 217,078)")9899    idx = RNG.choice(len(resid), size=100_000, replace=False)100101    fig, axes = plt.subplots(2, 2, figsize=(14, 10))102103    ax = axes[0, 0]104    ax.hexbin(fitted[idx], resid[idx], gridsize=60, cmap=CMAP_SEQ, mincnt=1)105    ax.axhline(0, color=ACCENT, linestyle="--", linewidth=1.5)106    ax.set_xlabel("Fitted Values")107    ax.set_ylabel("Residuals")108    panel_title(ax, "(a) Residuals vs. Fitted")109110    ax = axes[0, 1]111    (osm, osr), (slope, intercept, _) = st.probplot(resid[idx], dist="norm")112    ax.plot(osm, osr, ".", color=PRIMARY, markersize=2.5, rasterized=True)113    ax.plot(osm, slope * osm + intercept, color=ACCENT, linestyle="--", linewidth=1.5)114    ax.set_xlabel("Theoretical Quantiles")115    ax.set_ylabel("Sample Quantiles")116    panel_title(ax, "(b) Normal Q-Q Plot")117118    ax = axes[1, 0]119    ax.hist(resid, bins=100, color=PRIMARY, alpha=0.85, density=True,120            edgecolor="white", linewidth=0.2)121    x = np.linspace(resid.min(), resid.max(), 400)122    ax.plot(x, st.norm.pdf(x, resid.mean(), resid.std()), color=ACCENT,123            linestyle="--", linewidth=1.6, label="Normal density")124    ax.set_xlabel("Residuals")125    ax.set_ylabel("Density")126    panel_title(ax, "(c) Residual Distribution")127    ax.legend()128129    ax = axes[1, 1]130    ax.hexbin(fitted[idx], np.sqrt(np.abs(resid[idx] / resid.std())),131              gridsize=60, cmap=CMAP_SEQ, mincnt=1)132    ax.set_xlabel("Fitted Values")133    ax.set_ylabel(r"$\sqrt{|\mathrm{Standardized\ Residuals}|}$")134    panel_title(ax, "(d) Scale-Location")135136    fig.tight_layout()137    fig.savefig(path)138    plt.close(fig)139140141def fig3_quantile_coefficients(path):142    """Quantile-regression coefficient paths with the OLS benchmark."""143    qr = data.load_result("qr_results.pkl")144    md = data.load_model_data()145    ols_params = md["ols_results"]["params"]146    taus = sorted(qr.keys())147148    plot_vars = ["ln_sqft", "bathrooms", "age", "ln_lot",149                 "has_pool", "has_garage", "luxury_score", "tag_foreclosure"]150    print(f"  [verify] OLS ref ln_sqft={ols_params['ln_sqft']:.3f} "151          f"(original figure: 0.312)")152153    fig, axes = plt.subplots(2, 4, figsize=(19, 9))154    for k, (ax, v) in enumerate(zip(axes.ravel(), plot_vars)):155        coefs = [qr[t]["params"][v] for t in taus]156        ax.plot(taus, coefs, color=PRIMARY, marker="o", markersize=6,157                linewidth=2, zorder=3)158        ax.axhline(ols_params[v], color=ACCENT, linestyle="--", linewidth=1.5,159                   label="OLS")160        ax.set_xticks(taus)161        ax.set_xlabel(r"Quantile $\tau$")162        if k % 4 == 0:163            ax.set_ylabel("Coefficient")164        panel_title(ax, f"({chr(97 + k)}) {label(v)}")165        if k == 0:166            ax.legend(loc="best")167168    fig.tight_layout()169    fig.savefig(path)170    plt.close(fig)171172173def fig4_model_comparison(path):174    """Predicted vs actual log-prices on the random test set: OLS, XGB, LGBM."""175    from sklearn.metrics import r2_score176177    ml = data.load_result("ml_results.pkl")178    y_test = np.asarray(ml["y_test"], dtype=float)179    panels = [180        ("(a) OLS", ml["y_pred_ols"]),181        ("(b) XGBoost", ml["y_pred_xgb"]),182        ("(c) LightGBM", ml["y_pred_lgb"]),183    ]184185    idx = RNG.choice(len(y_test), size=50_000, replace=False)186    lims = (y_test.min() - 0.1, y_test.max() + 0.1)187188    fig, axes = plt.subplots(1, 3, figsize=(18, 6), sharex=True, sharey=True)189    for ax, (title, y_pred) in zip(axes, panels):190        y_pred = np.asarray(y_pred, dtype=float)191        r2 = r2_score(y_test, y_pred)192        print(f"  [verify] {title.split(') ')[1]}: R2={r2:.4f} "193              f"(paper: OLS 0.630 / XGBoost 0.833 / LightGBM 0.809)")194        ax.hexbin(y_test[idx], y_pred[idx], gridsize=70, cmap=CMAP_SEQ, mincnt=1)195        ax.plot(lims, lims, color=ACCENT, linestyle="--", linewidth=1.5)196        ax.set_xlim(lims)197        ax.set_ylim(lims)198        ax.set_xlabel("Actual Log(Price)")199        panel_title(ax, title)200        ax.text(0.05, 0.92, f"$R^2 = {r2:.3f}$", transform=ax.transAxes,201                fontsize=13)202    axes[0].set_ylabel("Predicted Log(Price)")203204    fig.tight_layout()205    fig.savefig(path)206    plt.close(fig)207208209def _shap_frame():210    """SHAP values and features with human-readable column labels."""211    sd = data.load_shap_data()212    X = sd["X_shap"].copy()213    return sd, X214215216def fig5_shap_summary(path):217    """SHAP beeswarm summary plot for the XGBoost model (top 20 features)."""218    import shap219220    sd, X = _shap_frame()221    X_lab = X.rename(columns={c: label(c) for c in X.columns})222    plt.figure(figsize=(12, 10))223    shap.summary_plot(sd["shap_values_xgb"], X_lab, max_display=20, show=False,224                      plot_size=None)225    ax = plt.gca()226    ax.set_xlabel("SHAP Value (impact on predicted log-price)")227    ax.grid(axis="y", visible=False)228    plt.tight_layout()229    plt.savefig(path)230    plt.close()231232233def fig6_shap_importance(path):234    """Mean |SHAP| feature-importance bar chart (top 20)."""235    sd, _ = _shap_frame()236    mean_shap = pd.Series(sd["mean_shap"]).sort_values(ascending=True).tail(20)237238    fig, ax = plt.subplots(figsize=(11, 9))239    ax.barh([label(v) for v in mean_shap.index], mean_shap.values,240            color=PRIMARY, alpha=0.9, height=0.72)241    ax.set_xlabel(r"Mean $|$SHAP Value$|$")242    ax.grid(axis="x", linestyle=":", linewidth=0.7, alpha=0.45, color="#999999")243    ax.grid(axis="y", visible=False)244    fig.tight_layout()245    fig.savefig(path)246    plt.close(fig)247248249def fig7_geographic_prices(path):250    """Map of log listing prices (50,000-listing random subsample)."""251    df = data.load_analytical_sample()252    idx = RNG.choice(len(df), size=50_000, replace=False)253    sub = df.iloc[idx]254255    fig, ax = plt.subplots(figsize=(15, 9))256    sc = ax.scatter(sub["longitude"], sub["latitude"], c=sub["ln_price"],257                    s=3.5, cmap="viridis", alpha=0.75, linewidths=0,258                    rasterized=True)259    cbar = fig.colorbar(sc, ax=ax, shrink=0.75, pad=0.02)260    cbar.set_label("Log(Price)")261    cbar.outline.set_visible(False)262    ax.set_xlabel("Longitude")263    ax.set_ylabel("Latitude")264    ax.grid(False)265    ax.set_axisbelow(True)266267    fig.tight_layout()268    fig.savefig(path)269    plt.close(fig)270271272def fig8_regional_prices(path):273    """Distribution of log listing prices by Census region."""274    df = data.load_analytical_sample()275    order = ["Northeast", "Midwest", "South", "West"]276    groups = [df.loc[df["region"] == r, "ln_price"].values for r in order]277    labels_n = [f"{r}\n(N = {len(g):,})" for r, g in zip(order, groups)]278279    fig, ax = plt.subplots(figsize=(10, 6.5))280    ax.boxplot(groups, tick_labels=labels_n, showfliers=False, patch_artist=True,281               widths=0.55,282               medianprops=dict(color=ACCENT, linewidth=1.8),283               boxprops=dict(facecolor=PRIMARY, alpha=0.75, edgecolor="#333333"),284               whiskerprops=dict(color="#333333", linewidth=1.0),285               capprops=dict(color="#333333", linewidth=1.0))286    ax.set_ylabel("Log(Price)")287    fig.tight_layout()288    fig.savefig(path)289    plt.close(fig)290291292def fig9_marginal_effects(path):293    """Unconditional bivariate relationships: scatter plus binned means."""294    df = data.load_analytical_sample()295    specs = [296        ("living_area_sqft", "Living Area (sqft)", (0, 8000)),297        ("age", "Property Age (years)", (0, 150)),298        ("lot_size_sqft", "Lot Size (sqft)", (0, 45000)),299        ("bedrooms", "Bedrooms", None),300        ("bathrooms", "Bathrooms", (0, 8)),301        ("avg_school_rating", "Avg. School Rating", None),302    ]303    idx = RNG.choice(len(df), size=30_000, replace=False)304    sub = df.iloc[idx]305306    fig, axes = plt.subplots(2, 3, figsize=(18, 10))307    for k, (ax, (col, xlabel, xlim)) in enumerate(zip(axes.ravel(), specs)):308        x, y = sub[col].values, sub["ln_price"].values309        if xlim is not None:310            keep = (x >= xlim[0]) & (x <= xlim[1])311            x, y = x[keep], y[keep]312        ax.scatter(x, y, s=2.5, color=PRIMARY_LIGHT, alpha=0.2, linewidths=0,313                   rasterized=True)314315        # Binned means computed on the full sample316        xf, yf = df[col].values, df["ln_price"].values317        if xlim is not None:318            keepf = (xf >= xlim[0]) & (xf <= xlim[1])319            xf, yf = xf[keepf], yf[keepf]320        bins = np.linspace(np.nanmin(xf), np.nanmax(xf), 25)321        which = np.digitize(xf, bins)322        centers = [xf[which == b].mean() for b in range(1, len(bins))323                   if (which == b).sum() > 50]324        means = [yf[which == b].mean() for b in range(1, len(bins))325                 if (which == b).sum() > 50]326        ax.plot(centers, means, color=ACCENT, linewidth=2.2, marker="o",327                markersize=4.5, label="Binned mean", zorder=3)328329        ax.set_xlabel(xlabel)330        if k % 3 == 0:331            ax.set_ylabel("Log(Price)")332        if col in ("living_area_sqft", "lot_size_sqft"):333            ax.xaxis.set_major_formatter(334                plt.matplotlib.ticker.StrMethodFormatter("{x:,.0f}"))335        panel_title(ax, f"({chr(97 + k)}) {xlabel}")336        if k == 0:337            ax.legend(loc="lower right")338339    fig.tight_layout()340    fig.savefig(path)341    plt.close(fig)342343344def fig10_shap_dependence(path):345    """SHAP dependence plots for the eight features shown in the paper."""346    sd, X = _shap_frame()347    shap_vals = sd["shap_values_xgb"]348    feats = ["ln_sqft", "bathrooms", "ln_lot", "avg_school_rating",349             "property_tax_rate", "age", "walk_score", "luxury_score"]350351    fig, axes = plt.subplots(2, 4, figsize=(20, 9.5))352    for k, (ax, feat) in enumerate(zip(axes.ravel(), feats)):353        j = list(X.columns).index(feat)354        xv = X[feat].values355        ax.scatter(xv, shap_vals[:, j], c=xv, cmap=CMAP_DIV, s=4, alpha=0.6,356                   linewidths=0, rasterized=True)357        ax.axhline(0, color=NEUTRAL, linestyle="--", linewidth=1.1)358        ax.set_xlabel(label(feat))359        if k % 4 == 0:360            ax.set_ylabel("SHAP Value")361        panel_title(ax, f"({chr(97 + k)}) {label(feat)}")362363    fig.tight_layout()364    fig.savefig(path)365    plt.close(fig)366367368def fig11_shap_interactions(path):369    """SHAP values of the two interaction features highlighted in the paper."""370    sd, X = _shap_frame()371    shap_vals = sd["shap_values_xgb"]372    panels = [373        ("waterfront_x_sqft", "(a) Waterfront $\\times$ Living Area Interaction"),374        ("pool_x_south", "(b) Pool $\\times$ South Region Interaction"),375    ]376377    fig, axes = plt.subplots(1, 2, figsize=(16, 6.5))378    for ax, (feat, title) in zip(axes, panels):379        j = list(X.columns).index(feat)380        ax.scatter(X[feat].values, shap_vals[:, j], s=6, color=PRIMARY,381                   alpha=0.35, linewidths=0, rasterized=True)382        ax.axhline(0, color=NEUTRAL, linestyle="--", linewidth=1.1)383        ax.set_xlabel(label(feat))384        ax.set_ylabel("SHAP Value")385        panel_title(ax, title)386387    fig.tight_layout()388    fig.savefig(path)389    plt.close(fig)390391392BUILDERS = {393    "fig1": fig1_price_distribution,394    "fig2": fig2_ols_diagnostics,395    "fig3": fig3_quantile_coefficients,396    "fig4": fig4_model_comparison,397    "fig5": fig5_shap_summary,398    "fig6": fig6_shap_importance,399    "fig7": fig7_geographic_prices,400    "fig8": fig8_regional_prices,401    "fig9": fig9_marginal_effects,402    "fig10": fig10_shap_dependence,403    "fig11": fig11_shap_interactions,404}405406407def main():408    parser = argparse.ArgumentParser(description=__doc__)409    parser.add_argument("--figs", nargs="+", choices=sorted(BUILDERS),410                        help="build only these figures")411    args = parser.parse_args()412413    apply_paper_style()414    config.FIGURES_DIR.mkdir(parents=True, exist_ok=True)415416    names = args.figs if args.figs else list(BUILDERS)417    for name in names:418        path = config.FIGURES_DIR / FILE_NAMES[name]419        print(f"Building {name} -> {path.relative_to(config.PROJECT_ROOT)}")420        BUILDERS[name](path)421422    print("Done.")423424425if __name__ == "__main__":426    main()427