spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Chapitre 40 — figures (DONNÉES SYNTHÉTIQUES : cas 3 et cas 4, graine 20261031).3Produit assets/figures/ch40/ch40-fig-{distribution,convergence,nuage,sobol,rcd,lhs,portefeuille}.pdf4Usage (racine) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch40_figures.py [--rapide]"""5from __future__ import annotations6import sys, pathlib7import numpy as np, pandas as pd8import matplotlib.pyplot as plt9from matplotlib.ticker import FuncFormatter, PercentFormatter10sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))11sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1] / "cases"))12from traite.style import setup, palette, savefig, LARGEUR, fmt_dollars13from traite import finance as fin14import ch40_montecarlo as m1516ETIQ = {"inoccupation_pertes": "Inoccupation et pertes", "croissance_loyers_marche": "Croissance loyers de marché", "rotation_annuelle": "Rotation annuelle",17        "croissance_depenses": "Croissance des dépenses", "taux_actualisation": "Taux d'actualisation r", "ecart_terminal": "Écart terminal r − R_T", "facteur_capex": "Facteur CAPEX", "taux_terminal": "Taux terminal R_T"}18MS = FuncFormatter(lambda x, p: f"{x/1e6:.1f}".replace(".", ",") + " M$")19fr = lambda v, d=2: f"{v:.{d}f}".replace(".", ",")2021def fig_distribution(D):22    V = D["mc"]["valeur"]; st = D["mc"]["stats"]23    fig, (a1, a2) = plt.subplots(1, 2, figsize=(LARGEUR, 2.9), gridspec_kw={"wspace": 0.28})24    a1.hist(V / 1e6, bins=70, range=(1.2, 3.8), density=True, color=palette.cat[0], alpha=0.85, edgecolor="white", linewidth=0.3)25    ymax = a1.get_ylim()[1]26    for q, lab in [(st["q05"], "q₅"), (st["mediane"], "médiane"), (st["q95"], "q₉₅")]:27        a1.axvline(q / 1e6, color=palette.ink, lw=0.9, ls="--"); a1.text(q / 1e6, ymax * 0.98, f"{lab}\n{fr(q/1e6)}", ha="center", va="top", fontsize=7)28    a1.axvline(D["dcf"]["valeur"] / 1e6, color=palette.cat[1], lw=1.5); a1.text(D["dcf"]["valeur"] / 1e6 + 0.3, ymax * 0.62, "DCF déterministe\n" + fr(D["dcf"]["valeur"] / 1e6), fontsize=7, color=palette.ink)29    a1.set_xlabel("Valeur simulée (millions de dollars)"); a1.set_ylabel("Densité"); a1.set_title("Histogramme (10 000 tirages)", loc="left")30    Vs = np.sort(V); F = np.arange(1, len(Vs) + 1) / len(Vs)31    a2.plot(Vs / 1e6, F, color=palette.cat[0], lw=1.6)32    q05 = st["q05"]; a2.fill_between(Vs[Vs <= q05] / 1e6, 0, F[Vs <= q05], color=palette.cat[2], alpha=0.35, lw=0)33    a2.axhline(0.05, color=palette.ink2, lw=0.7, ls=":"); a2.axvline(q05 / 1e6, color=palette.ink2, lw=0.7, ls=":")34    a2.annotate(f"q₅ = {fr(q05/1e6)} M$\nVaR₅ = {fmt_dollars(D['dcf']['valeur']-q05)}", xy=(q05 / 1e6, 0.05), xytext=(1.25, 0.42), fontsize=7, arrowprops=dict(arrowstyle="-", color=palette.ink2, lw=0.6))35    a2.axvline(D["valeur_capitalisation"] / 1e6, color=palette.cat[1], lw=1.2); a2.text(D["valeur_capitalisation"] / 1e6 + 0.03, 0.9, f"capitalisation\ndirecte : F = {fr(D['mc']['probabilites']['V_sous_cap'],3)}", fontsize=7)36    a2.set_xlim(1.2, 3.8); a2.set_ylim(0, 1); a2.set_xlabel("Valeur simulée (millions de dollars)"); a2.set_ylabel("Fonction de répartition empirique"); a2.set_title("Répartition et VaR de valeur", loc="left")37    savefig(fig, "ch40/ch40-fig-distribution.pdf")3839def fig_convergence(D):40    c = D["convergence"]; s = D["se_N"]41    fig, (a1, a2) = plt.subplots(1, 2, figsize=(LARGEUR, 2.9), gridspec_kw={"wspace": 0.3})42    for col, lab, k in [("q95", "quantile 95 %", 2), ("mediane", "médiane", 0), ("q05", "quantile 5 %", 1)]:43        a1.plot(c["N"], c[col] / 1e6, marker="o", ms=3, color=palette.cat[k], label=lab)44        a1.axhline(c[col].iloc[-1] / 1e6, color=palette.cat[k], lw=0.6, ls=":")45    st = D["mc"]["stats"]46    for col, k in [("q95", 2), ("mediane", 0), ("q05", 1)]:47        a1.scatter([10_000], [st[col] / 1e6], marker="D", s=22, color=palette.cat[k], edgecolor=palette.ink, linewidth=0.5, zorder=5)48    a1.text(10_000, st["mediane"] / 1e6 - 0.07, "graine\n20261031", ha="center", va="top", fontsize=6.8, color=palette.ink2)49    a1.set_xscale("log"); a1.set_xlabel("Nombre de tirages N (échelle log)"); a1.set_ylabel("Quantile cumulé (M$)"); a1.legend(loc="center right", fontsize=7); a1.set_title("Quantiles cumulés (suite de 100 000 tirages)", loc="left")50    a1.yaxis.set_major_formatter(FuncFormatter(lambda x, p: fr(x, 1)))51    for col, lab, k in [("se_mediane", "médiane", 0), ("se_q05", "quantile 5 %", 1), ("se_q95", "quantile 95 %", 2), ("se_moyenne", "moyenne", 3)]:52        a2.plot(s["N"], s[col], marker="o", ms=3, color=palette.cat[k], label=lab)53    ref = s["se_mediane"].iloc[0] * np.sqrt(s["N"].iloc[0] / s["N"])54    a2.plot(s["N"], ref, color=palette.ink2, lw=0.8, ls="--", label="pente −1/2")55    a2.set_xscale("log"); a2.set_yscale("log"); a2.set_xlabel("Nombre de tirages N (échelle log)"); a2.set_ylabel("Écart-type empirique de l'estimateur ($)"); a2.legend(fontsize=7, ncol=2); a2.set_title("Erreur Monte-Carlo (40 répétitions par N)", loc="left")56    a2.yaxis.set_major_formatter(FuncFormatter(lambda x, p: f"{x:,.0f}".replace(",", " ")))57    savefig(fig, "ch40/ch40-fig-convergence.pdf")5859def fig_nuage(D):60    x, V = D["mc"]["tirages"], D["mc"]["valeur"]; rng = np.random.default_rng(1); idx = rng.choice(len(V), 3000, replace=False)61    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 2.9), gridspec_kw={"wspace": 0.25})62    for ax, nom, xl in zip(axes, ["taux_actualisation", "taux_terminal"], ["Taux d'actualisation r (%)", "Taux terminal R_T (%)"]):63        xv = x[nom]; ax.scatter(xv[idx] * 100, V[idx] / 1e6, s=4, color=palette.cat[0], alpha=0.35, lw=0)64        bins = np.quantile(xv, np.linspace(0, 1, 21)); mids, meds = [], []65        for lo, hi in zip(bins[:-1], bins[1:]):66            sel = (xv >= lo) & (xv < hi); mids.append(xv[sel].mean() * 100); meds.append(np.median(V[sel]) / 1e6)67        ax.plot(mids, meds, color=palette.cat[1], lw=1.8, label="médiane conditionnelle")68        rho = D["rangs"].set_index("intrant").loc[nom, "spearman"]69        ax.text(0.03, 0.05, f"ρ de Spearman = {fr(rho)}", transform=ax.transAxes, fontsize=7.5)70        ax.set_xlabel(xl); ax.set_ylabel("Valeur simulée (M$)"); ax.yaxis.set_major_formatter(FuncFormatter(lambda v, p: fr(v, 1)))71    axes[0].legend(loc="upper right", fontsize=7); axes[0].set_title("Valeur et taux d'actualisation", loc="left"); axes[1].set_title("Valeur et taux terminal", loc="left")72    savefig(fig, "ch40/ch40-fig-nuage.pdf")7374def fig_sobol(D):75    so = D["sobol"].set_index("intrant"); rg = D["rangs"].set_index("intrant"); to = D["tornade"].set_index("intrant")76    ordre = so["ST"].sort_values().index.tolist(); y = np.arange(len(ordre)); h = 0.2677    part_tornade = (to["amplitude"] ** 2 / (to["amplitude"] ** 2).sum())78    fig, ax = plt.subplots(figsize=(LARGEUR, 3.0))79    ax.barh(y + h, so.loc[ordre, "ST"], h, color=palette.cat[0], label="indice total S_T (Sobol)")80    ax.barh(y, so.loc[ordre, "S1"], h, color=palette.cat[1], label="indice de premier ordre S₁ (Sobol)")81    ax.barh(y - h, rg.loc[ordre, "spearman2"], h, color=palette.cat[2], label="ρ² de Spearman (corrélation de rang)")82    ax.scatter(part_tornade.loc[ordre], y, marker="|", s=120, color=palette.ink, zorder=4, label="part de la tornade (amplitude² normalisée)")83    ax.set_yticks(y); ax.set_yticklabels([ETIQ[n] for n in ordre]); ax.set_xlabel("Part de la variance de la valeur expliquée"); ax.xaxis.set_major_formatter(PercentFormatter(1.0, decimals=0))84    ax.legend(loc="lower right", fontsize=7); ax.set_xlim(0, 0.62)85    for yi, n in zip(y, ordre): ax.text(so.loc[n, "ST"] + 0.008, yi + h, fr(so.loc[n, "ST"]), va="center", fontsize=6.8, color=palette.ink)86    savefig(fig, "ch40/ch40-fig-sobol.pdf")8788def fig_rcd(D):89    mc = D["mc"]; rne = mc["rne"]; f = D["financement"]; h = f["hypotheque"]; terme = D["hyp"]["financement"]["terme_ans"]90    h6 = fin.Hypotheque(h.solde(terme * 12), h.taux_nominal + 0.02, h.amortissement_mois - terme * 12)91    sd = np.where(np.arange(1, 11) <= terme, f["service_dette"], h6.service_dette_annuel)92    rcd = rne[:, :10] / sd; annees = np.arange(1, 11)93    fig, ax = plt.subplots(figsize=(LARGEUR, 3.0))94    for lo, hi, a in [(0.05, 0.95, 0.18), (0.25, 0.75, 0.32)]:95        ax.fill_between(annees, np.quantile(rcd, lo, axis=0), np.quantile(rcd, hi, axis=0), color=palette.cat[0], alpha=a, lw=0, label=f"bande {round((hi-lo)*100)} %")96    ax.plot(annees, np.median(rcd, axis=0), color=palette.cat[0], lw=1.8, label="médiane")97    rng = np.random.default_rng(3)98    for i in rng.choice(len(rcd), 12, replace=False): ax.plot(annees, rcd[i], color=palette.ink2, lw=0.5, alpha=0.6)99    ax.axhline(1.20, color=palette.cat[1], lw=1.0, ls="--"); ax.text(10.1, 1.20, "covenant 1,20", va="center", fontsize=7, color=palette.ink)100    ax.axhline(1.00, color=palette.cat[3], lw=1.0, ls="--"); ax.text(10.1, 1.00, "défaut 1,00", va="center", fontsize=7, color=palette.ink)101    ax.set_ylim(0.9, 2.15); ax.axvline(terme + 0.5, color=palette.ink2, lw=0.7, ls=":"); ax.text(terme + 0.6, 2.12, "renouvellement\nà +200 pb", fontsize=7, va="top", color=palette.ink2)102    ax.set_xlabel("Année d'exploitation"); ax.set_ylabel("Ratio de couverture de la dette"); ax.set_xticks(annees); ax.set_xlim(0.8, 11.2); ax.legend(loc="upper left", fontsize=7)103    ax.yaxis.set_major_formatter(FuncFormatter(lambda v, p: fr(v, 2)))104    savefig(fig, "ch40/ch40-fig-rcd.pdf")105106def fig_lhs(D):107    rng = np.random.default_rng(5); n = 12108    fig, (a1, a2, a3) = plt.subplots(1, 3, figsize=(LARGEUR, 2.6), gridspec_kw={"wspace": 0.85, "width_ratios": [0.85, 0.85, 1.7]})109    U1 = rng.random((n, 2)); U2 = m.hypercube_latin(n, 2, rng)110    for ax, U, t in [(a1, U1, "Simple"), (a2, U2, "Hypercube latin")]:111        for g in np.arange(0, 1.0001, 1 / n): ax.axhline(g, color=palette.grid, lw=0.5); ax.axvline(g, color=palette.grid, lw=0.5)112        ax.scatter(U[:, 0], U[:, 1], s=18, color=palette.cat[0], zorder=3); ax.set_xlim(0, 1); ax.set_ylim(0, 1); ax.set_aspect("equal"); ax.set_title(t + f", n = {n}", loc="left", fontsize=8)113        ax.set_xticks([0, 0.5, 1]); ax.set_yticks([0, 0.5, 1]); ax.grid(False); ax.set_xlabel("u₁"); ax.set_ylabel("u₂")114    l = D["lhs"].set_index("statistique"); noms = {"moyenne": "moyenne", "ecart_type": "écart-type", "mediane": "médiane", "q05": "q₅", "q95": "q₉₅", "p_V_2M": "P(V<2 M$)"}115    ordre = ["moyenne", "ecart_type", "mediane", "p_V_2M", "q95", "q05"]; vals = [min(l.loc[k, "ratio_variance"], 300) for k in ordre]116    a3.barh(range(len(ordre)), vals, color=palette.cat[1]); a3.set_yticks(range(len(ordre))); a3.set_yticklabels([noms[k] for k in ordre]); a3.set_xscale("log")117    a3.axvline(1, color=palette.ink, lw=0.8); a3.set_xlabel("Variance simple / variance hypercube latin (log)")118    for i, k in enumerate(ordre): a3.text(vals[i] * 1.12, i, f"× {l.loc[k, 'ratio_variance']:.1f}".replace(".", ","), va="center", fontsize=7)119    a3.set_xlim(0.5, 1000); a3.set_title("Gain de variance, N = 1 000, 100 répétitions", loc="left", fontsize=8)120    savefig(fig, "ch40/ch40-fig-lhs.pdf")121122def fig_portefeuille(D):123    pf, cr = D["cas4"]["portefeuille"], D["cas4_croissance"]124    fig, ax = plt.subplots(figsize=(LARGEUR, 2.9))125    for (nom, rp, rr), k in zip([("indépendance", 0.0, 0.0), ("corrélation régionale (0,30 ; 0,60)", 0.30, 0.60), ("corrélation parfaite", 1.0, 1.0)], range(3)):126        VP = m.simuler_portefeuille(pf, 10_000, m.GRAINE, rho_rmr=rr, rho_prov=rp, croissance=cr)["valeur_portefeuille"]127        Vs = np.sort(VP); ax.plot(Vs / 1e6, np.arange(1, len(Vs) + 1) / len(Vs), color=palette.cat[k], lw=1.6, label=nom)128    det = D["portefeuille"]["stats"]["valeur_deterministe"]129    ax.axvline(det / 1e6, color=palette.ink, lw=0.9, ls="--"); ax.text(det / 1e6 + 14, 0.12, f"valeur déterministe\nà un an : {fr(det/1e6,1)} M$", fontsize=7)130    ax.axhline(0.05, color=palette.ink2, lw=0.6, ls=":")131    ax.set_xlabel("Valeur du portefeuille de 48 immeubles à un an (millions de dollars)"); ax.set_ylabel("Fonction de répartition empirique"); ax.legend(loc="upper left", fontsize=7)132    ax.xaxis.set_major_formatter(FuncFormatter(lambda v, p: f"{v:.0f}"))133    savefig(fig, "ch40/ch40-fig-portefeuille.pdf")134135if __name__ == "__main__":136    setup(); rapide = "--rapide" in sys.argv137    D = m.tout(rapide=rapide)138    for f in [fig_distribution, fig_convergence, fig_nuage, fig_sobol, fig_rcd, fig_lhs, fig_portefeuille]:139        f(D); print("OK", f.__name__)140