spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Figures du chapitre 18 — Modèles hédoniques (cas 1, Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).3Produit dans assets/figures/ch18/ :4  ch18-fig-formes.pdf            prix prédit selon la superficie pour quatre formes fonctionnelles (autres attributs = sujet)5  ch18-fig-diagnostics.pdf       résidus vs valeurs ajustées et diagramme quantile-quantile du modèle de référence6  ch18-fig-predit-observe.pdf    prix prédit (Duan) vs prix observé, 600 ventes, cinq comparables du ch. 16 mis en évidence7  ch18-fig-surface.pdf           surface du prix prédit selon la superficie habitable et l'âge (autres attributs = sujet)8  ch18-fig-prix-implicites.pdf   prix implicites au sujet en dollars, avec erreurs-types HC1 (méthode delta)9  ch18-fig-indice.pdf            indice trimestriel par indicatrices vs tendance linéaire du modèle de référence10Usage (racine) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch18_figures.py"""11from __future__ import annotations12import sys, pathlib, math13ROOT = pathlib.Path(__file__).resolve().parents[3]14sys.path.insert(0, str(ROOT / "analysis/python")); sys.path.insert(0, str(ROOT / "analysis/python/cases"))15import numpy as np, pandas as pd16import matplotlib.pyplot as plt17from matplotlib.ticker import FuncFormatter18from scipy import stats19from traite.style import setup, palette, savefig, LARGEUR, fmt_milliers20from traite import cas121import ch18_hedonique as c182223def _k(x, pos=None): return f"{x/1000:,.0f}".replace(",", " ")2425def fig_formes(res):26    v = res["ventes"]; s = res["sujet"]; vari = res["variantes"]; ref = res["modele"]27    S = np.linspace(90, 260, 120)28    base = cas1.sujet_df(s); grille = pd.concat([base] * len(S), ignore_index=True); grille["superficie_habitable_m2"] = S29    duan = lambda r: float(np.mean(np.exp(r["residus"])))30    # semi-log (référence)31    y_ref = cas1.predire(ref, grille)32    # linéaire33    rL = vari["linéaire"]["mco"]; XL = c18.matrice_lineaire(grille)[vari["linéaire"]["X"].columns].values; y_lin = XL @ rL["beta"]34    # log-log35    rG = vari["log-log"]["mco"]; XG = c18.matrice_loglog(grille)[vari["log-log"]["X"].columns].values; y_log = np.exp(XG @ rG["beta"]) * duan(rG)36    # Box–Cox37    rB = vari["Box–Cox"]["mco"]; lam = vari["Box–Cox"]["lambda"]; XB = cas1.matrice(grille)[ref["X"].columns].values38    yl = XB @ rB["beta"]39    inv = lambda z: np.exp(z) if abs(lam) < 1e-8 else np.power(np.maximum(lam * z + 1, 1e-9), 1 / lam)40    y_bc = np.array([np.mean(inv(z + rB["residus"])) for z in yl])41    fig, ax = plt.subplots(figsize=(LARGEUR, 3.4))42    hull = v[(v["secteur"] == "Hull") & (v["type"] == "à étages")]43    ax.scatter(hull["superficie_habitable_m2"], hull["prix_vente"], s=9, color=palette.ink2, alpha=0.3, lw=0, label="Ventes de cottages à Hull (n = %d)" % len(hull))44    ax.plot(S, y_ref, color=palette.cat[0], lw=2.0, label="Semi-log (référence)")45    ax.plot(S, y_lin, color=palette.cat[1], lw=1.6, ls="--", label="Linéaire")46    ax.plot(S, y_log, color=palette.cat[2], lw=1.6, ls="-.", label="Log-log")47    ax.plot(S, y_bc, color=palette.cat[3], lw=1.6, ls=":", label=f"Box–Cox (λ = {lam:.2f})")48    ax.axvline(s["superficie_habitable_m2"], color=palette.ink2, lw=0.8, ls=":")49    ax.text(s["superficie_habitable_m2"] + 2, 330_000, "sujet (165 m²)", fontsize=7, color=palette.ink2)50    ax.set_xlabel("Superficie habitable (m²)"); ax.set_ylabel("Prix prédit (milliers de $)")51    ax.yaxis.set_major_formatter(FuncFormatter(_k)); ax.set_xlim(88, 262); ax.set_ylim(300_000, 950_000)52    ax.legend(loc="upper left", fontsize=7.5)53    return savefig(fig, "ch18/ch18-fig-formes.pdf")5455def fig_diagnostics(res):56    ref = res["modele"]; e = ref["residus"]; yhat = ref["X"].values @ ref["beta"]57    fig, (a1, a2) = plt.subplots(1, 2, figsize=(LARGEUR, 3.0))58    a1.scatter(yhat, e, s=8, color=palette.cat[0], alpha=0.45, lw=0)59    a1.axhline(0, color=palette.ink, lw=0.8)60    # médianes et écarts interquartiles par décile de valeurs ajustées61    dec = pd.qcut(yhat, 10, labels=False); df = pd.DataFrame({"y": yhat, "e": e, "d": dec}).groupby("d")62    a1.plot(df["y"].median(), df["e"].median(), color=palette.cat[3], marker="o", ms=4, lw=1.4, label="Médiane par décile")63    a1.plot(df["y"].median(), df["e"].quantile(0.9), color=palette.cat[3], lw=1.0, ls="--", label="Déciles 1 et 9")64    a1.plot(df["y"].median(), df["e"].quantile(0.1), color=palette.cat[3], lw=1.0, ls="--")65    a1.set_xlabel("ln(prix) ajusté"); a1.set_ylabel("Résidu (log)"); a1.set_title("Résidus et valeurs ajustées", fontsize=9)66    a1.legend(loc="upper left", fontsize=7)67    (osm, osr), (pente, ordo, r) = stats.probplot(e, dist="norm")68    a2.scatter(osm, osr, s=8, color=palette.cat[0], alpha=0.55, lw=0)69    xx = np.array([osm.min(), osm.max()]); a2.plot(xx, ordo + pente * xx, color=palette.ink, lw=0.9)70    a2.set_xlabel("Quantiles théoriques (loi normale)"); a2.set_ylabel("Quantiles des résidus"); a2.set_title("Diagramme quantile-quantile", fontsize=9)71    fig.tight_layout(w_pad=2.0)72    return savefig(fig, "ch18/ch18-fig-diagnostics.pdf")7374def fig_predit_observe(res):75    v = res["ventes"]; ref = res["modele"]; comps = res["comparables"]76    pred = res["variantes"]["semi-log (référence)"]["ajuste"]; P = v["prix_vente"].to_numpy(float); sig = ref["sigma"]77    fig, ax = plt.subplots(figsize=(LARGEUR, 4.1))78    ax.scatter(P, pred, s=9, color=palette.cat[0], alpha=0.4, lw=0, label="600 ventes")79    xx = np.linspace(200_000, 1_150_000, 10)80    ax.plot(xx, xx, color=palette.ink, lw=0.9, label="Prédit = observé")81    ax.plot(xx, xx * math.exp(2 * sig), color=palette.ink2, lw=0.8, ls="--", label="± 2 σ (± %.0f %%)" % (100 * (math.exp(2 * sig) - 1)))82    ax.plot(xx, xx * math.exp(-2 * sig), color=palette.ink2, lw=0.8, ls="--")83    pc = pd.Series(pred, index=v.index).loc[comps.index]84    ax.scatter(comps["prix_vente"], pc, s=46, facecolor="none", edgecolor=palette.cat[3], lw=1.4, label="Cinq comparables du chapitre 16", zorder=5)85    decal = {"G0574": (7, 4), "G0404": (-34, 6), "G0468": (7, 4), "G0477": (7, -3), "G0575": (7, -9)}86    for (i, c), p in zip(comps.iterrows(), pc):87        ax.annotate(c["id_vente"], (c["prix_vente"], p), xytext=decal[c["id_vente"]], textcoords="offset points", fontsize=7, color=palette.ink)88    ax.set_xlabel("Prix de vente observé (milliers de $)"); ax.set_ylabel("Prix prédit par le modèle (milliers de $)")89    ax.xaxis.set_major_formatter(FuncFormatter(_k)); ax.yaxis.set_major_formatter(FuncFormatter(_k))90    ax.set_xlim(200_000, 1_150_000); ax.set_ylim(200_000, 1_150_000); ax.set_aspect("equal")91    ax.legend(loc="upper left", fontsize=7.5)92    return savefig(fig, "ch18/ch18-fig-predit-observe.pdf")9394def fig_surface(res):95    s = res["sujet"]; ref = res["modele"]96    S = np.linspace(90, 260, 86); A = np.linspace(0, 100, 101); SS, AA = np.meshgrid(S, A)97    base = cas1.sujet_df(s); g = pd.concat([base] * SS.size, ignore_index=True)98    g["superficie_habitable_m2"] = SS.ravel(); g["age"] = AA.ravel()99    Z = cas1.predire(ref, g).reshape(SS.shape) / 1000100    fig, ax = plt.subplots(figsize=(LARGEUR, 3.9))101    cf = ax.contourf(SS, AA, Z, levels=18, cmap=palette.cmap_seq())102    cs = ax.contour(SS, AA, Z, levels=[400, 450, 500, 550, 600, 650, 700, 800, 900], colors=palette.ink, linewidths=0.6)103    ax.clabel(cs, fmt=lambda x: f"{x:.0f} k$", fontsize=7, colors=palette.ink)104    ax.plot(s["superficie_habitable_m2"], s["age"], marker="*", ms=11, color=palette.cat[3], mec=palette.ink, mew=0.6)105    ax.text(s["superficie_habitable_m2"] + 4, s["age"] + 3, "sujet : 165 m², 28 ans", fontsize=7.5, color=palette.ink)106    cb = fig.colorbar(cf, ax=ax, pad=0.02); cb.set_label("Prix prédit (milliers de $)")107    ax.set_xlabel("Superficie habitable (m²)"); ax.set_ylabel("Âge du bâtiment (années)"); ax.grid(False)108    return savefig(fig, "ch18/ch18-fig-surface.pdf")109110def fig_prix_implicites(res):111    pi = res["prix_implicites"].set_index("caracteristique")112    items = [("Secteur Aylmer plutôt que Hull", pi.loc["Secteur Aylmer vs Hull", "prix_implicite"], pi.loc["Secteur Aylmer vs Hull", "se"]),113             ("Qualité : un échelon de plus", pi.loc["Qualité (un échelon de plus)", "prix_implicite"], pi.loc["Qualité (un échelon de plus)", "se"]),114             ("Superficie habitable : +10 m²", 10 * pi.loc["Superficie habitable", "prix_implicite"], 10 * pi.loc["Superficie habitable", "se"]),115             ("Garage : une place de plus", pi.loc["Garage (une place de plus)", "prix_implicite"], pi.loc["Garage (une place de plus)", "se"]),116             ("Cottage à étages plutôt que plain-pied", pi.loc["Cottage à étages vs plain-pied", "prix_implicite"], pi.loc["Cottage à étages vs plain-pied", "se"]),117             ("Piscine", pi.loc["Piscine (ajout)", "prix_implicite"], pi.loc["Piscine (ajout)", "se"]),118             ("Sous-sol fini", pi.loc["Sous-sol fini (contribution : fini vs non fini)", "prix_implicite"], pi.loc["Sous-sol fini (contribution : fini vs non fini)", "se"]),119             ("Terrain : +100 m²", 100 * pi.loc["Superficie du terrain", "prix_implicite"], 100 * pi.loc["Superficie du terrain", "se"]),120             ("Conditions du marché : 12 mois", 12 * pi.loc["Conditions du marché (un mois)", "prix_implicite"], 12 * pi.loc["Conditions du marché (un mois)", "se"]),121             ("Secteur Gatineau plutôt que Hull", pi.loc["Secteur Gatineau vs Hull", "prix_implicite"], pi.loc["Secteur Gatineau vs Hull", "se"]),122             ("Âge : 10 années de plus (à 28 ans)", 10 * pi.loc["Âge (une année de plus, à 28 ans)", "prix_implicite"], 10 * pi.loc["Âge (une année de plus, à 28 ans)", "se"])]123    items = sorted(items, key=lambda t: t[1])124    fig, ax = plt.subplots(figsize=(LARGEUR, 3.8))125    y = np.arange(len(items)); vals = np.array([t[1] for t in items]); ses = np.array([t[2] for t in items])126    cols = [palette.cat[0] if x >= 0 else palette.cat[3] for x in vals]127    ax.barh(y, vals, color=cols, height=0.62, xerr=1.96 * ses, error_kw={"ecolor": palette.ink, "elinewidth": 0.8, "capsize": 2})128    ax.set_yticks(y); ax.set_yticklabels([t[0] for t in items], fontsize=8)129    ax.axvline(0, color=palette.ink, lw=0.8)130    for yi, x, se in zip(y, vals, ses):131        ax.text(x + (1.96 * se + 2000 if x >= 0 else -(1.96 * se + 2000)), yi, f"{x/1000:+.1f} k$".replace(".", ","), va="center", ha="left" if x >= 0 else "right", fontsize=7.5, color=palette.ink)132    ax.xaxis.set_major_formatter(FuncFormatter(_k)); ax.set_xlabel("Prix implicite au sujet (milliers de $) ; barres d'erreur : ± 1,96 erreur-type HC1")133    ax.set_xlim(-50_000, 105_000); ax.grid(axis="y", visible=False)134    return savefig(fig, "ch18/ch18-fig-prix-implicites.pdf")135136def fig_indice(res):137    idx = res["indice"]138    fig, ax = plt.subplots(figsize=(LARGEUR, 3.2))139    x = np.arange(len(idx))140    ax.fill_between(x, idx["bas"], idx["haut"], color=palette.cat[0], alpha=0.15, lw=0, label="IC à 95 % (HC1) des indicatrices")141    ax.plot(x, idx["indice_dummies"], color=palette.cat[0], marker="o", ms=4, lw=1.6, label="Indicatrices trimestrielles")142    ax.plot(x, idx["indice_tendance"], color=palette.cat[3], lw=1.8, ls="--", label="Tendance linéaire (0,36 %/mois, modèle de référence)")143    ax.set_xticks(x); ax.set_xticklabels(idx["trimestre"], rotation=45, ha="right", fontsize=7)144    ax.set_ylabel("Indice (2023T1 = 100)"); ax.set_ylim(96, 128)145    for xi, n in zip(x, idx["n_ventes"]):146        ax.text(xi, 97, str(n), ha="center", fontsize=6.5, color=palette.ink2)147    ax.text(x[0] - 0.55, 97, "n :", fontsize=6.5, color=palette.ink2, ha="right"); ax.set_xlim(-1.1, len(idx) - 0.5)148    ax.legend(loc="upper left", fontsize=7.5)149    return savefig(fig, "ch18/ch18-fig-indice.pdf")150151def main():152    setup()153    res = c18.main(verbose=False)154    for f in (fig_formes, fig_diagnostics, fig_predit_observe, fig_surface, fig_prix_implicites, fig_indice):155        print("OK", f(res))156157if __name__ == "__main__":158    main()159