Ch. 30 : réserve 675 $/logement alignée sur le code et les ch. 28/29/47
16 changed files +648 −5
modified
analysis/python/cases/ch20_reconciliation.py
+3 −1
@@ -196,7 +196,8 @@ def cout(b: dict) -> dict: | |||||
| 196 | def bootstrap_cout(b: dict, bs: dict, co: dict, B: int = B_COUT) -> dict: | 196 | def bootstrap_cout(b: dict, bs: dict, co: dict, B: int = B_COUT) -> dict: | ||
| 197 | """Bootstrap de l'indication par le coût du ch. 24 : mêmes rééchantillons des 600 ventes que le modèle (idxB) ; pour chaque | 197 | """Bootstrap de l'indication par le coût du ch. 24 : mêmes rééchantillons des 600 ventes que le modèle (idxB) ; pour chaque | ||
| 198 | réplication, prix ramenés à la date par la tendance du modèle rééstimé, terrain extrait des ventes récentes de Hull, | 198 | réplication, prix ramenés à la date par la tendance du modèle rééstimé, terrain extrait des ventes récentes de Hull, | ||
| 199 | − dépréciation extraite à 28 ans (ajustement quadratique), désuétude fonctionnelle fixe.""" | 199 | + dépréciation extraite à 28 ans (ajustement quadratique), désuétude fonctionnelle fixe. Les ventes d'un secteur sans vente | ||
| 200 | + récente dans la réplication (terrain non extractible) sont exclues de l'ajustement du profil.""" | ||||
| 200 | import ch24_depreciation as c24 | 201 | import ch24_depreciation as c24 | ||
| 201 | v = b["ventes"]; CN = co["CN"]; fonc = co["D_fonc"]; i_mois = b["cols"].index("mois") | 202 | v = b["ventes"]; CN = co["CN"]; fonc = co["D_fonc"]; i_mois = b["cols"].index("mois") | ||
| 202 | VK = np.empty(B); VL = np.empty(B); T28 = np.empty(B) | 203 | VK = np.empty(B); VL = np.empty(B); T28 = np.empty(B) | ||
@@ -204,6 +205,7 @@ def bootstrap_cout(b: dict, bs: dict, co: dict, B: int = B_COUT) -> dict: | |||||
| 204 | vb = v.iloc[bs["idxB"][r]].reset_index(drop=True) | 205 | vb = v.iloc[bs["idxB"][r]].reset_index(drop=True) | ||
| 205 | mb = {"croissance_mensuelle": float(np.exp(bs["betas"][r][i_mois]) - 1.0)} | 206 | mb = {"croissance_mensuelle": float(np.exp(bs["betas"][r][i_mois]) - 1.0)} | ||
| 206 | w, vl = c24.extraction(vb, mb, CN) | 207 | w, vl = c24.extraction(vb, mb, CN) | ||
| 208 | + w = w.dropna(subset=["taux"]) # secteurs sans vente récente dans la réplication : pas de terrain, exclus du profil | ||||
| 207 | q = c24.ajuster_profil(w)["quadratique"]; t28 = float(c24.taux_quadratique(q, AGE_EFFECTIF)) | 209 | q = c24.ajuster_profil(w)["quadratique"]; t28 = float(c24.taux_quadratique(q, AGE_EFFECTIF)) | ||
| 208 | VL[r] = float(vl.loc["Hull", "vl_mediane"]); T28[r] = t28 | 210 | VL[r] = float(vl.loc["Hull", "vl_mediane"]); T28[r] = t28 | ||
| 209 | VK[r] = ct.valeur_par_cout(VL[r], CN, t28 * CN, fonc)["valeur"] | 211 | VK[r] = ct.valeur_par_cout(VL[r], CN, t28 * CN, fonc)["valeur"] | ||
added
analysis/python/cases/ch40_montecarlo.py
+453 −0
@@ -0,0 +1,453 @@ | |||||
| 1 | +#!/usr/bin/env python3 | ||||
| 2 | +"""Chapitre 40 — Monte-Carlo et stress tests (DONNÉES SYNTHÉTIQUES : cas 3, 12 logements à Montréal ; cas 4, portefeuille). | ||||
| 3 | + | ||||
| 4 | +Ce script expose la MÉTHODE que le chapitre 47 (§47.12) APPLIQUE : il réutilise, sans les redéfinir, les lois `MC_LOIS`, | ||||
| 5 | +la graine 20261031 et les fonctions vectorisées `tirer`, `projeter_vectorise`, `valeur_vectorisee`, `monte_carlo` et | ||||
| 6 | +`tornade` du script `ch47_dossier.py`, de sorte que les chiffres des deux chapitres coïncident exactement | ||||
| 7 | +(médiane 2 191 608 $, IC 90 % [1 735 194 ; 2 876 902], P(RCD_1 < 1,20) = 49,4 %, P(RCD_1 < 1) = 0). | ||||
| 8 | + | ||||
| 9 | +Il y ajoute : moments théoriques des lois et paramétrage par quantiles (élicitation), bêta-PERT, corrélation r–R_T, | ||||
| 10 | +copule gaussienne / Cholesky, erreur Monte-Carlo (écart-type de la moyenne et des quantiles : théorie vs répétitions), | ||||
| 11 | +convergence des quantiles selon N, échantillonnage par hypercube latin (variance comparée à l'échantillonnage simple), | ||||
| 12 | +indices de Sobol de premier ordre et totaux (estimateur par échantillons appariés), corrélations de rang, mesures de | ||||
| 13 | +risque (VaR, CVaR de valeur, probabilités de bris de covenant, vérification analytique), stress tests déterministes | ||||
| 14 | +avec probabilité attachée, stress test inversé, agrégation d'un portefeuille de 48 immeubles avec corrélation régionale | ||||
| 15 | +(fonction réutilisable `simuler_portefeuille`, cas 4, chapitre 48) et les exercices du chapitre. | ||||
| 16 | + | ||||
| 17 | +Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch40_montecarlo.py [--rapide] | ||||
| 18 | +""" | ||||
| 19 | +from __future__ import annotations | ||||
| 20 | +import sys, pathlib, argparse, math | ||||
| 21 | +import numpy as np, pandas as pd | ||||
| 22 | +from scipy import stats, optimize | ||||
| 23 | +sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1])) | ||||
| 24 | +sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent)) | ||||
| 25 | +from traite import cas3, cas4, finance as fin | ||||
| 26 | +import ch47_dossier as c47 | ||||
| 27 | +from ch47_dossier import MC_LOIS, R_T_PLANCHER, tirer, projeter_vectorise, valeur_vectorisee, monte_carlo, tornade # noqa: F401 (réexportés) | ||||
| 28 | + | ||||
| 29 | +ROOT = pathlib.Path(__file__).resolve().parents[3] | ||||
| 30 | +GRAINE = c47.GRAINE # 20261031 | ||||
| 31 | +N_MC = 10_000 | ||||
| 32 | +Z90 = stats.norm.ppf(0.95) # 1,6449 | ||||
| 33 | + | ||||
| 34 | +# ----------------------------------------------------------------------------- 40.2 lois : moments théoriques, quantiles, élicitation | ||||
| 35 | + | ||||
| 36 | +def moments_loi(typ: str, p: dict) -> dict: | ||||
| 37 | + """Moyenne, écart-type, quantiles 5 % et 95 % théoriques d'une loi déclarée (mêmes paramétrages que MC_LOIS).""" | ||||
| 38 | + if typ == "lognormale": | ||||
| 39 | + mu, s = math.log(p["mediane"]), p["sigma_log"] | ||||
| 40 | + return {"moyenne": math.exp(mu + s**2 / 2), "ecart_type": math.sqrt((math.exp(s**2) - 1) * math.exp(2 * mu + s**2)), | ||||
| 41 | + "q05": math.exp(mu - Z90 * s), "q95": math.exp(mu + Z90 * s), "mediane": p["mediane"]} | ||||
| 42 | + if typ == "normale": | ||||
| 43 | + return {"moyenne": p["mu"], "ecart_type": p["sigma"], "q05": p["mu"] - Z90 * p["sigma"], "q95": p["mu"] + Z90 * p["sigma"], "mediane": p["mu"]} | ||||
| 44 | + if typ == "triangulaire": | ||||
| 45 | + a, m, b = p["min"], p["mode"], p["max"] | ||||
| 46 | + var = (a * a + b * b + m * m - a * b - a * m - b * m) / 18 | ||||
| 47 | + return {"moyenne": (a + m + b) / 3, "ecart_type": math.sqrt(var), "q05": quantile_triangulaire(0.05, a, m, b), "q95": quantile_triangulaire(0.95, a, m, b), | ||||
| 48 | + "mediane": quantile_triangulaire(0.5, a, m, b)} | ||||
| 49 | + raise ValueError(typ) | ||||
| 50 | + | ||||
| 51 | +def quantile_triangulaire(p, a, m, b): | ||||
| 52 | + Fm = (m - a) / (b - a) | ||||
| 53 | + p = np.asarray(p, dtype=float) | ||||
| 54 | + return np.where(p <= Fm, a + np.sqrt(p * (b - a) * (m - a)), b - np.sqrt((1 - p) * (b - a) * (b - m))) | ||||
| 55 | + | ||||
| 56 | +def quantile_loi(typ: str, p: dict, u): | ||||
| 57 | + """Fonction quantile F^{-1}(u) d'une loi déclarée (sert à l'hypercube latin et à la copule).""" | ||||
| 58 | + u = np.asarray(u, dtype=float) | ||||
| 59 | + if typ == "lognormale": return np.exp(math.log(p["mediane"]) + p["sigma_log"] * stats.norm.ppf(u)) | ||||
| 60 | + if typ == "normale": return p["mu"] + p["sigma"] * stats.norm.ppf(u) | ||||
| 61 | + if typ == "triangulaire": return quantile_triangulaire(u, p["min"], p["mode"], p["max"]) | ||||
| 62 | + raise ValueError(typ) | ||||
| 63 | + | ||||
| 64 | +def parametres_depuis_quantiles(q05: float, q95: float, loi: str = "normale") -> dict: | ||||
| 65 | + """Élicitation : un expert donne les quantiles 5 % et 95 % ; on en déduit les paramètres (normale ou log-normale).""" | ||||
| 66 | + if loi == "normale": return {"mu": (q05 + q95) / 2, "sigma": (q95 - q05) / (2 * Z90)} | ||||
| 67 | + if loi == "lognormale": return {"mediane": math.sqrt(q05 * q95), "sigma_log": math.log(q95 / q05) / (2 * Z90)} | ||||
| 68 | + raise ValueError(loi) | ||||
| 69 | + | ||||
| 70 | +def pert(a: float, m: float, b: float, lam: float = 4.0) -> dict: | ||||
| 71 | + """Loi bêta-PERT sur [a, b] de mode m : alpha = 1 + lam (m-a)/(b-a), beta = 1 + lam (b-m)/(b-a) ; moyenne (a + lam m + b)/(lam + 2).""" | ||||
| 72 | + al = 1 + lam * (m - a) / (b - a); be = 1 + lam * (b - m) / (b - a) | ||||
| 73 | + moy = (a + lam * m + b) / (lam + 2) | ||||
| 74 | + var = (moy - a) * (b - moy) / (lam + 3) | ||||
| 75 | + d = stats.beta(al, be, loc=a, scale=b - a) | ||||
| 76 | + return {"alpha": al, "beta": be, "moyenne": moy, "ecart_type": math.sqrt(var), "q05": float(d.ppf(0.05)), "q95": float(d.ppf(0.95)), "loi": d} | ||||
| 77 | + | ||||
| 78 | +def tableau_lois(x: dict | None = None) -> pd.DataFrame: | ||||
| 79 | + """Moments théoriques de chaque loi de MC_LOIS et, si des tirages sont fournis, leurs statistiques empiriques.""" | ||||
| 80 | + rows = [] | ||||
| 81 | + for nom, (typ, p) in MC_LOIS.items(): | ||||
| 82 | + m = moments_loi(typ, p); row = {"intrant": nom, "loi": typ, **{k: m[k] for k in ["moyenne", "ecart_type", "q05", "q95"]}} | ||||
| 83 | + if x is not None: | ||||
| 84 | + v = x[nom] if nom != "facteur_capex" else x[nom].ravel() | ||||
| 85 | + row.update({"moyenne_tiree": float(v.mean()), "q05_tire": float(np.quantile(v, 0.05)), "q95_tire": float(np.quantile(v, 0.95))}) | ||||
| 86 | + rows.append(row) | ||||
| 87 | + return pd.DataFrame(rows) | ||||
| 88 | + | ||||
| 89 | +# ----------------------------------------------------------------------------- 40.2 corrélations : Cholesky, copule gaussienne, lien r–R_T | ||||
| 90 | + | ||||
| 91 | +def cholesky_2x2(s1: float, s2: float, rho: float) -> np.ndarray: | ||||
| 92 | + """Facteur de Cholesky L de la matrice de covariance [[s1², rho s1 s2], [rho s1 s2, s2²]] : Sigma = L L'.""" | ||||
| 93 | + return np.array([[s1, 0.0], [rho * s2, s2 * math.sqrt(1 - rho**2)]]) | ||||
| 94 | + | ||||
| 95 | +def copule_gaussienne(n: int, graine: int, corr: np.ndarray, lois: list[tuple[str, dict]]) -> np.ndarray: | ||||
| 96 | + """Tirages corrélés par copule gaussienne : Z = L eps, U = Phi(Z), X_k = F_k^{-1}(U_k). Retourne (n, d).""" | ||||
| 97 | + rng = np.random.default_rng(graine) | ||||
| 98 | + L = np.linalg.cholesky(corr); Z = rng.standard_normal((n, len(lois))) @ L.T | ||||
| 99 | + U = stats.norm.cdf(Z) | ||||
| 100 | + return np.column_stack([quantile_loi(t, p, U[:, k]) for k, (t, p) in enumerate(lois)]) | ||||
| 101 | + | ||||
| 102 | +def spearman_depuis_pearson_gaussien(rho: float) -> float: | ||||
| 103 | + """Corrélation de rang induite par une copule gaussienne de corrélation rho : rho_S = (6/pi) arcsin(rho/2).""" | ||||
| 104 | + return 6 / math.pi * math.asin(rho / 2) | ||||
| 105 | + | ||||
| 106 | +def correlation_r_RT(x: dict | None = None) -> dict: | ||||
| 107 | + """Corrélation théorique entre r et R_T = r − delta (sans plancher) et, si tirages fournis, corrélation empirique.""" | ||||
| 108 | + sr = MC_LOIS["taux_actualisation"][1]["sigma"]; sd = MC_LOIS["ecart_terminal"][1]["sigma"] | ||||
| 109 | + out = {"rho_theorique": sr / math.sqrt(sr**2 + sd**2), "sigma_RT_theorique": math.sqrt(sr**2 + sd**2)} | ||||
| 110 | + if x is not None: | ||||
| 111 | + out["rho_empirique"] = float(np.corrcoef(x["taux_actualisation"], x["taux_terminal"])[0, 1]) | ||||
| 112 | + out["spearman_empirique"] = float(stats.spearmanr(x["taux_actualisation"], x["taux_terminal"])[0]) | ||||
| 113 | + out["part_plancher"] = float((x["taux_actualisation"] - x["ecart_terminal"] < R_T_PLANCHER).mean()) | ||||
| 114 | + out["sigma_RT_empirique"] = float(x["taux_terminal"].std(ddof=1)) | ||||
| 115 | + out["spearman_theorique"] = spearman_depuis_pearson_gaussien(out["rho_theorique"]) | ||||
| 116 | + return out | ||||
| 117 | + | ||||
| 118 | +# ----------------------------------------------------------------------------- 40.3 simulation : évaluation vectorisée, erreur Monte-Carlo, convergence | ||||
| 119 | + | ||||
| 120 | +def evaluer(r: dict, x: dict) -> dict: | ||||
| 121 | + """Valeur DCF et RNE pour un dictionnaire de tirages (même structure que `tirer`).""" | ||||
| 122 | + proj = projeter_vectorise(r["etat"], r["hyp"], r["base"], x["inoccupation_pertes"], x["croissance_loyers_marche"], x["rotation_annuelle"], x["croissance_depenses"], x["facteur_capex"]) | ||||
| 123 | + val = valeur_vectorisee(proj, x["taux_actualisation"], x["taux_terminal"], r["hyp"]["couts_vente"]) | ||||
| 124 | + return {"valeur": val["valeur"], "rne": proj["rne"], "rcd1": proj["rne"][:, 0] / r["financement"]["service_dette"]} | ||||
| 125 | + | ||||
| 126 | +def erreur_type_quantile(V: np.ndarray, p: float) -> dict: | ||||
| 127 | + """Écart-type asymptotique d'un quantile empirique : sqrt(p(1-p)/N) / f(q_p), densité f estimée par noyau gaussien.""" | ||||
| 128 | + N = len(V); q = float(np.quantile(V, p)); f = float(stats.gaussian_kde(V)(q)[0]) | ||||
| 129 | + return {"p": p, "quantile": q, "densite": f, "se_theorique": math.sqrt(p * (1 - p) / N) / f} | ||||
| 130 | + | ||||
| 131 | +def erreurs_theoriques(mc: dict) -> pd.DataFrame: | ||||
| 132 | + V = mc["valeur"]; N = len(V); rows = [{"statistique": "moyenne", "estimation": float(V.mean()), "se_theorique": float(V.std(ddof=1) / math.sqrt(N))}] | ||||
| 133 | + for p, nom in [(0.05, "q05"), (0.5, "mediane"), (0.95, "q95")]: | ||||
| 134 | + e = erreur_type_quantile(V, p); rows.append({"statistique": nom, "estimation": e["quantile"], "se_theorique": e["se_theorique"], "densite": e["densite"]}) | ||||
| 135 | + rcd1 = mc["rcd"][:, 0]; p = float((rcd1 < 1.2).mean()) | ||||
| 136 | + rows.append({"statistique": "P(RCD1<1,20)", "estimation": p, "se_theorique": math.sqrt(p * (1 - p) / N)}) | ||||
| 137 | + return pd.DataFrame(rows) | ||||
| 138 | + | ||||
| 139 | +def repetitions(r: dict, n: int, k: int, graine0: int = 1, lhs: bool = False) -> pd.DataFrame: | ||||
| 140 | + """k simulations indépendantes de n tirages (graines graine0, graine0+1, …) ; statistiques de chaque répétition.""" | ||||
| 141 | + rows = [] | ||||
| 142 | + nb_capex = len(r["hyp"]["capex"]) | ||||
| 143 | + for g in range(graine0, graine0 + k): | ||||
| 144 | + x = tirer_lhs(n, g, nb_capex) if lhs else tirer(n, g, nb_capex) | ||||
| 145 | + e = evaluer(r, x); V = e["valeur"] | ||||
| 146 | + rows.append({"graine": g, "moyenne": V.mean(), "q05": np.quantile(V, 0.05), "mediane": np.median(V), "q95": np.quantile(V, 0.95), "ecart_type": V.std(ddof=1), | ||||
| 147 | + "p_rcd1_120": (e["rcd1"] < 1.2).mean(), "p_V_2M": (V < 2e6).mean()}) | ||||
| 148 | + return pd.DataFrame(rows) | ||||
| 149 | + | ||||
| 150 | +def convergence(r: dict, N_max: int = 100_000, graine: int = GRAINE, points: tuple = (100, 200, 500, 1000, 2000, 5000, 10_000, 20_000, 50_000, 100_000)) -> pd.DataFrame: | ||||
| 151 | + """Quantiles cumulés d'une seule suite de N_max tirages (graine), lus après N tirages (N dans `points`).""" | ||||
| 152 | + x = tirer(N_max, graine, len(r["hyp"]["capex"])); V = evaluer(r, x)["valeur"] | ||||
| 153 | + rows = [{"N": N, "q05": np.quantile(V[:N], 0.05), "mediane": np.median(V[:N]), "q95": np.quantile(V[:N], 0.95), "moyenne": V[:N].mean()} for N in points if N <= N_max] | ||||
| 154 | + return pd.DataFrame(rows) | ||||
| 155 | + | ||||
| 156 | +def se_empirique_selon_N(r: dict, Ns=(250, 500, 1000, 2000, 4000, 8000), k: int = 40) -> pd.DataFrame: | ||||
| 157 | + """Écart-type empirique (sur k répétitions) de la médiane, de q05 et de la moyenne selon N : doit décroître en N^(-1/2).""" | ||||
| 158 | + rows = [] | ||||
| 159 | + for N in Ns: | ||||
| 160 | + rep = repetitions(r, N, k, graine0=100) | ||||
| 161 | + rows.append({"N": N, "se_mediane": rep["mediane"].std(ddof=1), "se_q05": rep["q05"].std(ddof=1), "se_moyenne": rep["moyenne"].std(ddof=1), "se_q95": rep["q95"].std(ddof=1)}) | ||||
| 162 | + d = pd.DataFrame(rows) | ||||
| 163 | + b = np.polyfit(np.log(d["N"]), np.log(d["se_mediane"]), 1)[0] | ||||
| 164 | + d.attrs["pente_log_mediane"] = float(b) | ||||
| 165 | + return d | ||||
| 166 | + | ||||
| 167 | +def n_requis(se_unitaire: float, erreur_cible: float, confiance: float = 0.95) -> int: | ||||
| 168 | + """N tel que z_(1-a/2) · se_unitaire / sqrt(N) <= erreur_cible, où se_unitaire est l'écart-type « pour N = 1 ».""" | ||||
| 169 | + z = stats.norm.ppf(0.5 + confiance / 2) | ||||
| 170 | + return int(math.ceil((z * se_unitaire / erreur_cible) ** 2)) | ||||
| 171 | + | ||||
| 172 | +# ----------------------------------------------------------------------------- 40.4 hypercube latin | ||||
| 173 | + | ||||
| 174 | +def hypercube_latin(n: int, d: int, rng: np.random.Generator) -> np.ndarray: | ||||
| 175 | + """n points dans [0,1]^d : pour chaque dimension, une permutation des n strates [j/n, (j+1)/n) et un point uniforme dans chaque strate.""" | ||||
| 176 | + U = np.empty((n, d)) | ||||
| 177 | + for k in range(d): | ||||
| 178 | + U[:, k] = (rng.permutation(n) + rng.random(n)) / n | ||||
| 179 | + return U | ||||
| 180 | + | ||||
| 181 | +def tirer_lhs(n: int, graine: int, nb_capex: int, lois=MC_LOIS) -> dict: | ||||
| 182 | + """Même structure que `ch47_dossier.tirer` (mêmes lois, même ordre, mêmes bornes) mais par hypercube latin et fonction quantile.""" | ||||
| 183 | + rng = np.random.default_rng(graine) | ||||
| 184 | + noms = ["inoccupation_pertes", "croissance_loyers_marche", "rotation_annuelle", "croissance_depenses", "taux_actualisation", "ecart_terminal"] | ||||
| 185 | + U = hypercube_latin(n, len(noms) + nb_capex, rng) | ||||
| 186 | + x = {nom: quantile_loi(*lois[nom], U[:, k]) for k, nom in enumerate(noms)} | ||||
| 187 | + x["inoccupation_pertes"] = np.clip(x["inoccupation_pertes"], 0.0, 0.30); x["taux_actualisation"] = np.clip(x["taux_actualisation"], 0.04, None) | ||||
| 188 | + typ, p = lois["facteur_capex"] | ||||
| 189 | + x["facteur_capex"] = np.column_stack([quantile_loi(typ, p, U[:, len(noms) + j]) for j in range(nb_capex)]) | ||||
| 190 | + x["taux_terminal"] = np.maximum(x["taux_actualisation"] - x["ecart_terminal"], R_T_PLANCHER) | ||||
| 191 | + return x | ||||
| 192 | + | ||||
| 193 | +def comparer_lhs(r: dict, n: int = 1000, k: int = 100) -> pd.DataFrame: | ||||
| 194 | + """Variance des estimateurs (moyenne, médiane, q05, q95, P(RCD1<1,20)) sur k répétitions : échantillonnage simple vs hypercube latin.""" | ||||
| 195 | + simple = repetitions(r, n, k, graine0=1000, lhs=False); lhs = repetitions(r, n, k, graine0=1000, lhs=True) | ||||
| 196 | + rows = [] | ||||
| 197 | + for s in ["moyenne", "mediane", "q05", "q95", "ecart_type", "p_rcd1_120", "p_V_2M"]: | ||||
| 198 | + rows.append({"statistique": s, "moy_simple": simple[s].mean(), "moy_lhs": lhs[s].mean(), "se_simple": simple[s].std(ddof=1), "se_lhs": lhs[s].std(ddof=1), | ||||
| 199 | + "ratio_variance": simple[s].var(ddof=1) / lhs[s].var(ddof=1)}) | ||||
| 200 | + return pd.DataFrame(rows) | ||||
| 201 | + | ||||
| 202 | +# ----------------------------------------------------------------------------- 40.5 mesures de risque, vérification analytique, stress tests | ||||
| 203 | + | ||||
| 204 | +def mesures_risque(mc: dict, V_ref: float, alphas=(0.05, 0.10)) -> pd.DataFrame: | ||||
| 205 | + """VaR_alpha = V_ref − q_alpha (perte par rapport à la valeur de référence) ; CVaR_alpha = V_ref − E[V | V <= q_alpha].""" | ||||
| 206 | + V = mc["valeur"]; rows = [] | ||||
| 207 | + for a in alphas: | ||||
| 208 | + q = float(np.quantile(V, a)); cvar_v = float(V[V <= q].mean()) | ||||
| 209 | + rows.append({"alpha": a, "quantile": q, "VaR": V_ref - q, "VaR_pct": (V_ref - q) / V_ref, "esperance_queue": cvar_v, "CVaR": V_ref - cvar_v, "CVaR_pct": (V_ref - cvar_v) / V_ref}) | ||||
| 210 | + return pd.DataFrame(rows) | ||||
| 211 | + | ||||
| 212 | +def cvar_normale(mu: float, sigma: float, alpha: float) -> float: | ||||
| 213 | + """E[V | V <= q_alpha] pour V ~ N(mu, sigma²) : mu − sigma · phi(z_alpha)/alpha.""" | ||||
| 214 | + z = stats.norm.ppf(alpha); return mu - sigma * stats.norm.pdf(z) / alpha | ||||
| 215 | + | ||||
| 216 | +def rcd1_analytique(r: dict) -> dict: | ||||
| 217 | + """RCD de première année comme fonction affine de l'inoccupation v : RCD1(v) = [0,96 (RBP(1−v) + A) − DE_hors_gestion − RR] / SD ; | ||||
| 218 | + seuils v* (RCD1 = 1,20) et v_1 (RCD1 = 1,00) et probabilités exactes sous la loi log-normale déclarée.""" | ||||
| 219 | + b, hyp, f = r["base"], r["hyp"], r["financement"] | ||||
| 220 | + rbp, A, RR, SD = b["rbp_contrat"], hyp["autres_revenus_annuels"], b["reserve"], f["service_dette"] | ||||
| 221 | + de_hg = b["depenses_hors_reserve"] - 0.04 * b["rbe"] | ||||
| 222 | + rcd = lambda v: (0.96 * (rbp * (1 - v) + A) - de_hg - RR) / SD | ||||
| 223 | + v_seuil = lambda cible: 1 - ((cible * SD + de_hg + RR) / 0.96 - A) / rbp | ||||
| 224 | + med, s = MC_LOIS["inoccupation_pertes"][1]["mediane"], MC_LOIS["inoccupation_pertes"][1]["sigma_log"] | ||||
| 225 | + P = lambda v: float(1 - stats.norm.cdf(math.log(v / med) / s)) | ||||
| 226 | + v120, v100 = v_seuil(1.20), v_seuil(1.00) | ||||
| 227 | + return {"pente_par_point": (rcd(0.01) - rcd(0.0)), "rcd1_base": rcd(hyp["taux_inoccupation_pertes"]), "v_seuil_120": v120, "v_seuil_100": v100, | ||||
| 228 | + "P_rcd1_sous_120_exacte": P(v120), "P_rcd1_sous_100_exacte": P(v100), "attendu_sous_100_sur_10000": 10_000 * P(v100), "rcd1_v30": rcd(0.30)} | ||||
| 229 | + | ||||
| 230 | +def valeur_deterministe(r: dict, **par) -> float: | ||||
| 231 | + """Valeur DCF pour un jeu d'intrants déterministes (défaut : hypothèses de base) ; clés = noms de MC_LOIS (+ taux_terminal direct).""" | ||||
| 232 | + hyp = r["hyp"] | ||||
| 233 | + base = {"inoccupation_pertes": hyp["taux_inoccupation_pertes"], "croissance_loyers_marche": hyp["croissance_loyers_marche"], "rotation_annuelle": hyp["rotation_annuelle"], | ||||
| 234 | + "croissance_depenses": hyp["croissance_depenses"], "taux_actualisation": hyp["taux_actualisation"], "ecart_terminal": hyp["taux_actualisation"] - hyp["taux_terminal"], "facteur_capex": 1.0} | ||||
| 235 | + base.update({k: v for k, v in par.items() if k != "taux_terminal"}) | ||||
| 236 | + R_T = par.get("taux_terminal", max(base["taux_actualisation"] - base["ecart_terminal"], R_T_PLANCHER)) | ||||
| 237 | + x = {k: np.array([base[k]]) for k in ["inoccupation_pertes", "croissance_loyers_marche", "rotation_annuelle", "croissance_depenses", "taux_actualisation"]} | ||||
| 238 | + x["facteur_capex"] = np.full((1, len(hyp["capex"])), base["facteur_capex"]); x["taux_terminal"] = np.array([R_T]) | ||||
| 239 | + return float(evaluer(r, x)["valeur"][0]) | ||||
| 240 | + | ||||
| 241 | +def stress_deterministes(r: dict, mc: dict) -> pd.DataFrame: | ||||
| 242 | + """Scénarios déterministes classiques et probabilité, dans la simulation, d'un résultat au moins aussi défavorable.""" | ||||
| 243 | + hyp = r["hyp"]; V = mc["valeur"]; q = lambda nom, p: float(quantile_loi(*MC_LOIS[nom], p)) | ||||
| 244 | + sc = [("Base (hypothèses du dossier)", {}), | ||||
| 245 | + ("Hausse des taux : r + 75 pb, R_T + 50 pb (ch. 30)", {"taux_actualisation": 0.08, "taux_terminal": 0.055}), | ||||
| 246 | + ("Marché détendu : inoccupation 6 %", {"inoccupation_pertes": 0.06}), | ||||
| 247 | + ("Croissance faible : loyers de marché 2,5 %, dépenses 4 %", {"croissance_loyers_marche": 0.025, "croissance_depenses": 0.04}), | ||||
| 248 | + ("Trois intrants à leur quantile défavorable 90 %", {"taux_actualisation": q("taux_actualisation", 0.90), "croissance_loyers_marche": q("croissance_loyers_marche", 0.10), "inoccupation_pertes": q("inoccupation_pertes", 0.90)}), | ||||
| 249 | + ("Tous les intrants à leur quantile défavorable 90 %", {"taux_actualisation": q("taux_actualisation", 0.90), "ecart_terminal": q("ecart_terminal", 0.10), "croissance_loyers_marche": q("croissance_loyers_marche", 0.10), | ||||
| 250 | + "inoccupation_pertes": q("inoccupation_pertes", 0.90), "croissance_depenses": q("croissance_depenses", 0.90), "rotation_annuelle": q("rotation_annuelle", 0.10), "facteur_capex": q("facteur_capex", 0.90)})] | ||||
| 251 | + rows = [] | ||||
| 252 | + for nom, par in sc: | ||||
| 253 | + v = valeur_deterministe(r, **par) | ||||
| 254 | + rows.append({"scenario": nom, "valeur": v, "ecart_pct": v / r["dcf"]["valeur"] - 1, "P_V_inferieure": float((V <= v).mean())}) | ||||
| 255 | + return pd.DataFrame(rows) | ||||
| 256 | + | ||||
| 257 | +def stress_inverse(r: dict) -> dict: | ||||
| 258 | + """Stress test inversé : quel intrant, seul, ramène la valeur à un seuil donné ? (inoccupation pour RCD1 = 1 ; r pour V = 2 M$ ; R_T pour V = V_cap).""" | ||||
| 259 | + an = rcd1_analytique(r) | ||||
| 260 | + f_r = lambda rr: valeur_deterministe(r, taux_actualisation=rr, taux_terminal=max(rr - 0.0225, R_T_PLANCHER)) - 2_000_000 | ||||
| 261 | + r_2M = optimize.brentq(f_r, 0.05, 0.12) | ||||
| 262 | + f_RT = lambda RT: valeur_deterministe(r, taux_terminal=RT) - r["valeur_capitalisation"] | ||||
| 263 | + RT_cap = optimize.brentq(f_RT, 0.035, 0.09) | ||||
| 264 | + f_g = lambda g: valeur_deterministe(r, croissance_loyers_marche=g) - r["valeur_capitalisation"] | ||||
| 265 | + g_cap = optimize.brentq(f_g, -0.05, 0.035) | ||||
| 266 | + return {"v_rcd1_egal_1": an["v_seuil_100"], "r_valeur_2M": r_2M, "RT_valeur_cap": RT_cap, "g_marche_valeur_cap": g_cap} | ||||
| 267 | + | ||||
| 268 | +# ----------------------------------------------------------------------------- 40.6 sensibilité globale : Sobol (échantillons appariés) et corrélations de rang | ||||
| 269 | + | ||||
| 270 | +NOMS_SOBOL = ["inoccupation_pertes", "croissance_loyers_marche", "rotation_annuelle", "croissance_depenses", "taux_actualisation", "ecart_terminal", "facteur_capex"] | ||||
| 271 | + | ||||
| 272 | +def _completer(x: dict) -> dict: | ||||
| 273 | + x = dict(x); x["taux_terminal"] = np.maximum(x["taux_actualisation"] - x["ecart_terminal"], R_T_PLANCHER); return x | ||||
| 274 | + | ||||
| 275 | +def sobol(r: dict, n: int = N_MC, graine: int = GRAINE) -> pd.DataFrame: | ||||
| 276 | + """Indices de Sobol de premier ordre S_i et totaux S_Ti par échantillons appariés : deux matrices indépendantes A et B (n tirages | ||||
| 277 | + chacune), puis A_B^(i) = A où la seule colonne i vient de B. S_i ≈ (1/n) Σ f(B)·[f(A_B^(i)) − f(A)] / Var ; S_Ti ≈ (1/2n) Σ [f(A) − f(A_B^(i))]² / Var. | ||||
| 278 | + Les intrants de base (r et delta, non r et R_T) sont indépendants, condition de validité de la décomposition. Le facteur CAPEX est traité en groupe (3 facteurs).""" | ||||
| 279 | + nb = len(r["hyp"]["capex"]) | ||||
| 280 | + A = tirer(n, graine, nb); B = tirer(n, graine + 1, nb) | ||||
| 281 | + fA = evaluer(r, A)["valeur"]; fB = evaluer(r, B)["valeur"] | ||||
| 282 | + var = np.concatenate([fA, fB]).var(ddof=1); rows = [] | ||||
| 283 | + for nom in NOMS_SOBOL: | ||||
| 284 | + AB = dict(A); AB[nom] = B[nom]; fAB = evaluer(r, _completer(AB))["valeur"] | ||||
| 285 | + S1 = float(np.mean(fB * (fAB - fA)) / var); ST = float(0.5 * np.mean((fA - fAB) ** 2) / var) | ||||
| 286 | + rows.append({"intrant": nom, "S1": S1, "ST": ST}) | ||||
| 287 | + d = pd.DataFrame(rows); d.attrs["somme_S1"] = float(d["S1"].sum()); d.attrs["variance"] = float(var); d.attrs["n_evaluations"] = n * (2 + len(NOMS_SOBOL)) | ||||
| 288 | + return d | ||||
| 289 | + | ||||
| 290 | +def correlations_rang(mc: dict) -> pd.DataFrame: | ||||
| 291 | + """Corrélations de Spearman et de Pearson entre chaque intrant et la valeur ; rho² comme approximation de S1 pour un modèle quasi additif.""" | ||||
| 292 | + x, V = mc["tirages"], mc["valeur"]; rows = [] | ||||
| 293 | + for nom in NOMS_SOBOL + ["taux_terminal"]: | ||||
| 294 | + v = x[nom].mean(1) if nom == "facteur_capex" else x[nom] | ||||
| 295 | + rows.append({"intrant": nom, "spearman": float(stats.spearmanr(v, V)[0]), "pearson": float(np.corrcoef(v, V)[0, 1])}) | ||||
| 296 | + d = pd.DataFrame(rows); d["spearman2"] = d["spearman"] ** 2; d["pearson2"] = d["pearson"] ** 2 | ||||
| 297 | + return d | ||||
| 298 | + | ||||
| 299 | +def regression_standardisee(mc: dict) -> dict: | ||||
| 300 | + """Coefficients standardisés d'une régression linéaire de V sur les sept intrants : R² mesure la part « additive » de la variance.""" | ||||
| 301 | + x, V = mc["tirages"], mc["valeur"] | ||||
| 302 | + X = np.column_stack([x[n].mean(1) if n == "facteur_capex" else x[n] for n in NOMS_SOBOL]) | ||||
| 303 | + Xs = (X - X.mean(0)) / X.std(0, ddof=1); Vs = (V - V.mean()) / V.std(ddof=1) | ||||
| 304 | + beta, *_ = np.linalg.lstsq(np.column_stack([np.ones(len(V)), Xs]), Vs, rcond=None) | ||||
| 305 | + pred = np.column_stack([np.ones(len(V)), Xs]) @ beta | ||||
| 306 | + return {"beta": dict(zip(NOMS_SOBOL, beta[1:])), "R2": float(1 - ((Vs - pred) ** 2).sum() / (Vs ** 2).sum())} | ||||
| 307 | + | ||||
| 308 | +# ----------------------------------------------------------------------------- 40.8 agrégation : portefeuille de 48 immeubles avec corrélation régionale (cas 4 → ch. 48) | ||||
| 309 | + | ||||
| 310 | +def simuler_portefeuille(pf: pd.DataFrame, n: int = N_MC, graine: int = GRAINE, sigma_R: float = 0.005, rho_rmr: float = 0.60, rho_prov: float = 0.30, | ||||
| 311 | + croissance: pd.DataFrame | None = None, sigma_g_defaut: float = 0.03, choc_taux_renouvellement: float = 0.02, horizon_renouvellement: int = 2028, | ||||
| 312 | + amortissement_mois: int = 300, seuil_rcd: float = 1.20) -> dict: | ||||
| 313 | + """Simulation à un an de la valeur et de la couverture de la dette d'un portefeuille d'immeubles à revenus avec corrélation régionale. | ||||
| 314 | + | ||||
| 315 | + Modèle (hypothèses pédagogiques déclarées au ch. 40, §40.8) : pour l'immeuble i de la RMR k, | ||||
| 316 | + ln RNE_i' = ln RNE_i + g_i + sigma_g,i · e_i, e_i = sqrt(rho_prov) Z_0 + sqrt(rho_rmr − rho_prov) Z_k + sqrt(1 − rho_rmr) eps_i, | ||||
| 317 | + R_i' = R_i + sigma_R · u_i, u_i de même structure factorielle (facteurs indépendants de ceux du RNE), | ||||
| 318 | + V_i' = RNE_i' / R_i', | ||||
| 319 | + de sorte que Corr(e_i, e_j) = rho_rmr si même RMR, rho_prov sinon (0 <= rho_prov <= rho_rmr <= 1). Le RCD est RNE_i' / SD_i, le service | ||||
| 320 | + de la dette étant celui de la dette `dette` au taux `taux_dette` (composition semestrielle, amortissement `amortissement_mois`) ou, si | ||||
| 321 | + `echeance_dette` <= horizon_renouvellement, au taux majoré de `choc_taux_renouvellement` (renouvellement sous stress). | ||||
| 322 | + `croissance` : DataFrame (id_immeuble, croissance_annuelle, volatilite) de `traite.cas4.croissance_rne` ; sinon g = 0, sigma_g = sigma_g_defaut. | ||||
| 323 | + Retourne les valeurs simulées (n, N_immeubles), la valeur du portefeuille, les RCD, et un dictionnaire de statistiques. | ||||
| 324 | + Réutilisable par le chapitre 48 : `from ch40_montecarlo import simuler_portefeuille`.""" | ||||
| 325 | + assert 0 <= rho_prov <= rho_rmr <= 1 | ||||
| 326 | + rng = np.random.default_rng(graine) | ||||
| 327 | + pf = pf.reset_index(drop=True); N = len(pf); rmrs = pf["rmr"].astype("category"); k = rmrs.cat.codes.values; K = len(rmrs.cat.categories) | ||||
| 328 | + g = np.zeros(N); sg = np.full(N, sigma_g_defaut) | ||||
| 329 | + if croissance is not None: | ||||
| 330 | + c = croissance.set_index("id_immeuble").reindex(pf["id_immeuble"]); g = c["croissance_annuelle"].values; sg = c["volatilite"].values | ||||
| 331 | + def facteur(): | ||||
| 332 | + Z0 = rng.standard_normal((n, 1)); Zk = rng.standard_normal((n, K)); eps = rng.standard_normal((n, N)) | ||||
| 333 | + return math.sqrt(rho_prov) * Z0 + math.sqrt(max(rho_rmr - rho_prov, 0.0)) * Zk[:, k] + math.sqrt(1 - rho_rmr) * eps | ||||
| 334 | + e = facteur(); u = facteur() | ||||
| 335 | + rne = pf["rne_2026"].values * np.exp(g + sg * e) | ||||
| 336 | + R = np.maximum(pf["taux_global"].values + sigma_R * u, 0.02) | ||||
| 337 | + V = rne / R | ||||
| 338 | + # service de la dette : amortissement résiduel non observé → amortissement `amortissement_mois` (hypothèse), renouvellement sous stress si échéance proche | ||||
| 339 | + sd = np.array([fin.Hypotheque(float(d), float(t) + (choc_taux_renouvellement if int(ech) <= horizon_renouvellement else 0.0), amortissement_mois).service_dette_annuel | ||||
| 340 | + for d, t, ech in zip(pf["dette"], pf["taux_dette"], pf["echeance_dette"])]) | ||||
| 341 | + rcd = rne / sd | ||||
| 342 | + VP = V.sum(1); V0 = float(pf["valeur_2026"].sum()); VP_det = float((pf["rne_2026"].values * np.exp(g) / pf["taux_global"].values).sum()) | ||||
| 343 | + q = lambda a, p: float(np.quantile(a, p)) | ||||
| 344 | + stats_ = {"n": n, "N_immeubles": N, "valeur_2026": V0, "valeur_deterministe": VP_det, "mediane": q(VP, 0.5), "moyenne": float(VP.mean()), "ecart_type": float(VP.std(ddof=1)), | ||||
| 345 | + "q05": q(VP, 0.05), "q95": q(VP, 0.95), "VaR05": VP_det - q(VP, 0.05), "CVaR05": VP_det - float(VP[VP <= q(VP, 0.05)].mean()), | ||||
| 346 | + "somme_ecarts_types": float(V.std(0, ddof=1).sum()), "nb_rcd_sous_seuil_moyen": float((rcd < seuil_rcd).sum(1).mean()), "nb_rcd_sous_1_moyen": float((rcd < 1.0).sum(1).mean()), | ||||
| 347 | + "P_au_moins_un_rcd_sous_1": float(((rcd < 1.0).any(1)).mean()), "P_au_moins_un_rcd_sous_seuil": float(((rcd < seuil_rcd).any(1)).mean()), | ||||
| 348 | + "rcd_portefeuille_mediane": q(rne.sum(1) / sd.sum(), 0.5), "P_rcd_portefeuille_sous_seuil": float((rne.sum(1) / sd.sum() < seuil_rcd).mean()), | ||||
| 349 | + "nb_renouvellements_stress": int((pf["echeance_dette"] <= horizon_renouvellement).sum())} | ||||
| 350 | + stats_["ratio_diversification"] = stats_["ecart_type"] / stats_["somme_ecarts_types"] | ||||
| 351 | + return {"valeurs": V, "valeur_portefeuille": VP, "rne": rne, "rcd": rcd, "service_dette": sd, "stats": stats_, | ||||
| 352 | + "par_rmr": pd.DataFrame({"rmr": rmrs.cat.categories, "ecart_type_valeur": [V[:, k == j].sum(1).std(ddof=1) for j in range(K)], "valeur_2026": [pf.loc[k == j, "valeur_2026"].sum() for j in range(K)]})} | ||||
| 353 | + | ||||
| 354 | +def comparer_structures_correlation(pf: pd.DataFrame, croissance: pd.DataFrame, n: int = N_MC) -> pd.DataFrame: | ||||
| 355 | + """Trois structures : indépendance (0, 0), corrélation régionale (0,30 ; 0,60), corrélation parfaite (1, 1).""" | ||||
| 356 | + rows = [] | ||||
| 357 | + for nom, rp, rr in [("Indépendance", 0.0, 0.0), ("Régionale (0,30 ; 0,60)", 0.30, 0.60), ("Parfaite", 1.0, 1.0)]: | ||||
| 358 | + s = simuler_portefeuille(pf, n, GRAINE, rho_rmr=rr, rho_prov=rp, croissance=croissance)["stats"] | ||||
| 359 | + rows.append({"structure": nom, "rho_prov": rp, "rho_rmr": rr, **{k: s[k] for k in ["mediane", "ecart_type", "q05", "VaR05", "CVaR05", "ratio_diversification", "nb_rcd_sous_seuil_moyen", "nb_rcd_sous_1_moyen", "P_au_moins_un_rcd_sous_1", "P_rcd_portefeuille_sous_seuil"]}}) | ||||
| 360 | + return pd.DataFrame(rows) | ||||
| 361 | + | ||||
| 362 | +def variance_somme_correlee(sigmas: np.ndarray, rho: float) -> float: | ||||
| 363 | + """Écart-type de la somme de N variables d'écarts-types sigma_i et de corrélation commune rho : sqrt(Σσ_i² + rho Σ_{i≠j} σ_i σ_j).""" | ||||
| 364 | + s2 = (sigmas ** 2).sum(); tot = sigmas.sum() ** 2 | ||||
| 365 | + return math.sqrt(s2 + rho * (tot - s2)) | ||||
| 366 | + | ||||
| 367 | +# ----------------------------------------------------------------------------- exercices | ||||
| 368 | + | ||||
| 369 | +def exercices(r: dict, mc: dict) -> dict: | ||||
| 370 | + out = {} | ||||
| 371 | + # Ex. 1 : élicitation — l'expert donne q05 = 2,4 %, q95 = 6,3 % d'inoccupation → log-normale ; moyenne, P(v > 5 %) | ||||
| 372 | + p = parametres_depuis_quantiles(0.024, 0.063, "lognormale"); m = moments_loi("lognormale", p) | ||||
| 373 | + out["ex1"] = {**p, **m, "P_sup_5pct": float(1 - stats.norm.cdf(math.log(0.05 / p["mediane"]) / p["sigma_log"]))} | ||||
| 374 | + # Ex. 2 : PERT vs triangulaire pour la rotation (0,08 ; 0,15 ; 0,25) | ||||
| 375 | + t = moments_loi("triangulaire", MC_LOIS["rotation_annuelle"][1]); pe = pert(0.08, 0.15, 0.25) | ||||
| 376 | + out["ex2"] = {"tri": t, "pert": {k: pe[k] for k in ["alpha", "beta", "moyenne", "ecart_type", "q05", "q95"]}} | ||||
| 377 | + # Ex. 3 : nombre de tirages pour une médiane à ±5 000 $ (95 %) et pour P(RCD<1,20) à ±1 point | ||||
| 378 | + V = mc["valeur"]; e = erreur_type_quantile(V, 0.5) | ||||
| 379 | + se1_med = e["se_theorique"] * math.sqrt(len(V)) # écart-type « pour N = 1 » | ||||
| 380 | + out["ex3"] = {"se_mediane_10000": e["se_theorique"], "densite": e["densite"], "N_mediane_5000": n_requis(se1_med, 5000), "N_mediane_1000": n_requis(se1_med, 1000), | ||||
| 381 | + "N_proba_1pt": n_requis(0.5, 0.01), "N_proba_0.5pt": n_requis(0.5, 0.005)} | ||||
| 382 | + # Ex. 4 : Cholesky 2×2 pour (r, R_T) avec rho = 0,78 et sigma_RT = 0,0064 ; vérification Sigma = L L' ; Spearman induit | ||||
| 383 | + L = cholesky_2x2(0.005, 0.0064, 0.78); out["ex4"] = {"L": L, "Sigma": L @ L.T, "spearman": spearman_depuis_pearson_gaussien(0.78)} | ||||
| 384 | + # Ex. 5 : VaR/CVaR 10 % de valeur et loi normale de même moyenne/écart-type ; comparaison | ||||
| 385 | + st = mc["stats"]; mr = mesures_risque(mc, r["dcf"]["valeur"], alphas=(0.10,)).iloc[0].to_dict() | ||||
| 386 | + out["ex5"] = {**mr, "q10_normale": stats.norm.ppf(0.10, st["moyenne"], st["ecart_type"]), "cvar_normale": cvar_normale(st["moyenne"], st["ecart_type"], 0.10)} | ||||
| 387 | + # Ex. 6 : simulation avec sigma_r doublé (100 pb) : écart-type, IC 90 %, P(V<2 M$) | ||||
| 388 | + lois = dict(MC_LOIS); lois["taux_actualisation"] = ("normale", {"mu": 0.0725, "sigma": 0.010}) | ||||
| 389 | + x = tirer(N_MC, GRAINE, len(r["hyp"]["capex"]), lois); ev = evaluer(r, x); Vx = ev["valeur"] | ||||
| 390 | + out["ex6"] = {"mediane": float(np.median(Vx)), "ecart_type": float(Vx.std(ddof=1)), "q05": float(np.quantile(Vx, 0.05)), "q95": float(np.quantile(Vx, 0.95)), "P_V_2M": float((Vx < 2e6).mean()), | ||||
| 391 | + "P_V_2M_base": float((V < 2e6).mean()), "cv": float(Vx.std(ddof=1) / Vx.mean())} | ||||
| 392 | + return out | ||||
| 393 | + | ||||
| 394 | +# ----------------------------------------------------------------------------- tout | ||||
| 395 | + | ||||
| 396 | +def tout(rapide: bool = False) -> dict: | ||||
| 397 | + r = cas3.tout() | ||||
| 398 | + mc = monte_carlo(r, N_MC, GRAINE); torn = tornade(r) | ||||
| 399 | + x = mc["tirages"] | ||||
| 400 | + D = {**r, "mc": mc, "tornade": torn, "lois": tableau_lois(x), "corr_rRT": correlation_r_RT(x), "erreurs": erreurs_theoriques(mc), | ||||
| 401 | + "risque": mesures_risque(mc, r["dcf"]["valeur"]), "rcd_analytique": rcd1_analytique(r), "stress": stress_deterministes(r, mc), "inverse": stress_inverse(r), | ||||
| 402 | + "sobol": sobol(r), "rangs": correlations_rang(mc), "reg": regression_standardisee(mc), "pert_rotation": pert(0.08, 0.15, 0.25)} | ||||
| 403 | + k_rep = 20 if rapide else 100; k_se = 10 if rapide else 40 | ||||
| 404 | + D["repetitions"] = repetitions(r, N_MC, k_rep, graine0=1) | ||||
| 405 | + D["convergence"] = convergence(r, N_max=20_000 if rapide else 100_000) | ||||
| 406 | + D["se_N"] = se_empirique_selon_N(r, k=k_se) | ||||
| 407 | + D["lhs"] = comparer_lhs(r, n=1000, k=k_rep) | ||||
| 408 | + c4 = cas4.charger(); cr = cas4.croissance_rne(c4["series"]) | ||||
| 409 | + D["portefeuille"] = simuler_portefeuille(c4["portefeuille"], N_MC, GRAINE, croissance=cr); D["portefeuille_structures"] = comparer_structures_correlation(c4["portefeuille"], cr) | ||||
| 410 | + D["cas4"] = c4; D["cas4_croissance"] = cr | ||||
| 411 | + D["exercices"] = exercices(r, mc) | ||||
| 412 | + return D | ||||
| 413 | + | ||||
| 414 | +def _d(x): return f"{x:,.0f}".replace(",", " ") + " $" | ||||
| 415 | +def _p(x, d=2): return f"{100*x:.{d}f}".replace(".", ",") + " %" | ||||
| 416 | + | ||||
| 417 | +def imprimer(D: dict) -> None: | ||||
| 418 | + pd.set_option("display.width", 200); pd.set_option("display.max_columns", 20); pd.set_option("display.float_format", lambda v: f"{v:,.6g}") | ||||
| 419 | + mc, st, pr = D["mc"], D["mc"]["stats"], D["mc"]["probabilites"] | ||||
| 420 | + print("=== Cas 3 : DCF déterministe", _d(D["dcf"]["valeur"]), "| capitalisation", _d(D["valeur_capitalisation"]), "| prêt", _d(D["financement"]["pret"]), "| SD", _d(D["financement"]["service_dette"])) | ||||
| 421 | + print("\n=== 40.2 Lois déclarées : moments théoriques et tirés\n", D["lois"].to_string(index=False)) | ||||
| 422 | + print("\nPERT rotation (0,08 ; 0,15 ; 0,25) :", {k: round(v, 4) for k, v in D["pert_rotation"].items() if k != "loi"}) | ||||
| 423 | + print("Corrélation r–R_T :", {k: round(v, 4) for k, v in D["corr_rRT"].items()}) | ||||
| 424 | + print("\n=== 40.3 Monte-Carlo (N = 10 000, graine 20261031) — identique au ch. 47") | ||||
| 425 | + print(f"moyenne {_d(st['moyenne'])} | écart-type {_d(st['ecart_type'])} | CV {_p(st['cv'],1)} | médiane {_d(st['mediane'])} | q05 {_d(st['q05'])} | q25 {_d(st['q25'])} | q75 {_d(st['q75'])} | q95 {_d(st['q95'])}") | ||||
| 426 | + print(f"P(V<cap) {_p(pr['V_sous_cap'],1)} | P(V<2M) {_p(pr['V_sous_2M'],1)} | P(V<rôle) {_p(pr['V_sous_role'],1)} | P(RCD1<1,20) {_p(pr['rcd1_sous_1.20'],1)} | P(RCD1<1) {pr['rcd1_sous_1.00']} | P(RCD6 +200pb <1,20) {_p(pr['rcd6_choc200_sous_1.20'],1)}") | ||||
| 427 | + print(f"RCD1 q05 {mc['rcd1_q05']:.3f} | RCD6 choc médiane {mc['rcd6_choc200_mediane']:.3f} q05 {mc['rcd6_choc200_q05']:.3f}") | ||||
| 428 | + print("\nErreur Monte-Carlo théorique (N = 10 000) :\n", D["erreurs"].to_string(index=False)) | ||||
| 429 | + rep = D["repetitions"] | ||||
| 430 | + print(f"\nRépétitions ({len(rep)} graines, N = 10 000) : écart-type empirique — moyenne {rep['moyenne'].std(ddof=1):,.0f} | médiane {rep['mediane'].std(ddof=1):,.0f} | q05 {rep['q05'].std(ddof=1):,.0f} | q95 {rep['q95'].std(ddof=1):,.0f} | P(RCD1<1,20) {rep['p_rcd1_120'].std(ddof=1):.4f}") | ||||
| 431 | + print(f" moyenne des médianes {rep['mediane'].mean():,.0f} ; min {rep['mediane'].min():,.0f} ; max {rep['mediane'].max():,.0f} ; P(RCD1<1,20) moyenne {rep['p_rcd1_120'].mean():.4f}") | ||||
| 432 | + print("\nConvergence (une suite, graine 20261031) :\n", D["convergence"].to_string(index=False)) | ||||
| 433 | + print(f"\nÉcart-type empirique selon N (pente log-log médiane = {D['se_N'].attrs['pente_log_mediane']:.3f}) :\n", D["se_N"].to_string(index=False)) | ||||
| 434 | + print("\n=== 40.4 Hypercube latin vs simple (N = 1 000, répétitions) :\n", D["lhs"].to_string(index=False)) | ||||
| 435 | + print("\n=== 40.5 Mesures de risque (référence : valeur DCF déterministe)\n", D["risque"].to_string(index=False)) | ||||
| 436 | + print("RCD1 analytique :", {k: round(v, 5) for k, v in D["rcd_analytique"].items()}) | ||||
| 437 | + print("CVaR 5 % si V était normale :", _d(cvar_normale(st["moyenne"], st["ecart_type"], 0.05)), "| q05 normal :", _d(stats.norm.ppf(0.05, st["moyenne"], st["ecart_type"]))) | ||||
| 438 | + print("\nStress tests déterministes :\n", D["stress"].to_string(index=False)) | ||||
| 439 | + print("Stress inversé :", {k: round(v, 5) for k, v in D["inverse"].items()}) | ||||
| 440 | + print("\n=== 40.6 Sobol (S1, ST ; somme S1 = %.3f ; %d évaluations)\n" % (D["sobol"].attrs["somme_S1"], D["sobol"].attrs["n_evaluations"]), D["sobol"].to_string(index=False)) | ||||
| 441 | + print("\nCorrélations de rang / Pearson :\n", D["rangs"].to_string(index=False)) | ||||
| 442 | + print("Régression standardisée : R² = %.4f ; bêtas =" % D["reg"]["R2"], {k: round(v, 3) for k, v in D["reg"]["beta"].items()}) | ||||
| 443 | + print("\nTornade (ch. 47) :\n", D["tornade"].to_string(index=False)) | ||||
| 444 | + print("\n=== 40.8 Portefeuille (cas 4, 48 immeubles) — structure régionale (0,30 ; 0,60)") | ||||
| 445 | + print({k: (round(v, 4) if isinstance(v, float) else v) for k, v in D["portefeuille"]["stats"].items()}) | ||||
| 446 | + print(D["portefeuille"]["par_rmr"].to_string(index=False)) | ||||
| 447 | + print("\nTrois structures de corrélation :\n", D["portefeuille_structures"].to_string(index=False)) | ||||
| 448 | + print("\n=== Exercices") | ||||
| 449 | + for k, v in D["exercices"].items(): print(k, v) | ||||
| 450 | + | ||||
| 451 | +if __name__ == "__main__": | ||||
| 452 | + ap = argparse.ArgumentParser(); ap.add_argument("--rapide", action="store_true"); a = ap.parse_args() | ||||
| 453 | + imprimer(tout(rapide=a.rapide)) | ||||
added
analysis/python/figures/ch38_figures.py
+187 −0
@@ -0,0 +1,187 @@ | |||||
| 1 | +#!/usr/bin/env python3 | ||||
| 2 | +"""Figures du chapitre 38 — SIG et risques géospatiaux (cas 1 Gatineau et cas 4 portefeuille, DONNÉES SYNTHÉTIQUES, graine 20261031). | ||||
| 3 | +Cartes geopandas/Matplotlib sans fond de tuiles (contextily absent) ; coordonnées projetées NAD83 / MTM zone 9 (cas 1) et | ||||
| 4 | +NAD83 / Québec Lambert (cas 4), exprimées en kilomètres ; toutes les géométries « zones », « corridor » et « pôles » sont FICTIVES. | ||||
| 5 | +Sorties (assets/figures/ch38/) : | ||||
| 6 | + ch38-fig-carte-cas1.pdf 600 ventes, zones d'exposition fictives, corridor fictif et ses tampons, pôles fictifs, sujet | ||||
| 7 | + ch38-fig-choroplethe.pdf médiane du prix au m² par cellule hexagonale selon trois classifications + résidus moyens (divergent) | ||||
| 8 | + ch38-fig-faux-positifs.pdf distribution des t de polygones aléatoires et part de « significatifs » selon le rayon et la spécification | ||||
| 9 | + ch38-fig-tampon.pdf coefficient de l'indicatrice de corridor selon le rayon du tampon (IC à 95 %) | ||||
| 10 | + ch38-fig-cas4.pdf 48 immeubles du cas 4 par RMR et part de la valeur en zone inondable | ||||
| 11 | +Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch38_figures.py""" | ||||
| 12 | +import sys, pathlib, warnings | ||||
| 13 | +ROOT = pathlib.Path(__file__).resolve().parents[3] | ||||
| 14 | +sys.path.insert(0, str(ROOT / "analysis/python")); sys.path.insert(0, str(ROOT / "analysis/python/cases")) | ||||
| 15 | +import numpy as np | ||||
| 16 | +import matplotlib.pyplot as plt | ||||
| 17 | +from matplotlib.colors import ListedColormap, BoundaryNorm | ||||
| 18 | +from matplotlib.lines import Line2D | ||||
| 19 | +from matplotlib.patches import Patch | ||||
| 20 | +from scipy.stats import norm | ||||
| 21 | +from traite.style import setup, palette, savefig, LARGEUR, DEMI, fmt_milliers | ||||
| 22 | +from traite import geo | ||||
| 23 | +import ch38_sig as c38 | ||||
| 24 | + | ||||
| 25 | +SECTEURS = ["Hull", "Gatineau", "Aylmer", "Buckingham", "Masson-Angers"] | ||||
| 26 | + | ||||
| 27 | + | ||||
| 28 | +def _km(ax, xy, marge=1.5): | ||||
| 29 | + ax.set_aspect("equal"); ax.set_xlim(xy[:, 0].min() - marge, xy[:, 0].max() + marge); ax.set_ylim(xy[:, 1].min() - marge, xy[:, 1].max() + marge) | ||||
| 30 | + ax.grid(True, lw=0.4) | ||||
| 31 | + | ||||
| 32 | + | ||||
| 33 | +def _poly_km(ax, geom, **kw): | ||||
| 34 | + """Trace un (multi)polygone shapely en km (coordonnées MTM divisées par 1000).""" | ||||
| 35 | + polys = getattr(geom, "geoms", [geom]) | ||||
| 36 | + for p in polys: | ||||
| 37 | + ext = np.asarray(p.exterior.coords) / 1000 | ||||
| 38 | + ax.fill(ext[:, 0], ext[:, 1], **kw) | ||||
| 39 | + for hole in p.interiors: | ||||
| 40 | + hh = np.asarray(hole.coords) / 1000; ax.fill(hh[:, 0], hh[:, 1], color=palette.surface, lw=0) | ||||
| 41 | + | ||||
| 42 | + | ||||
| 43 | +def fig_carte_cas1(r): | ||||
| 44 | + b, gf, v = r["base"], r["geom"], r["ventes"]; xy = b["xy_m"] / 1000; cs = np.array([b["xs"], b["ys"]]) / 1000 | ||||
| 45 | + fig, ax = plt.subplots(figsize=(LARGEUR, 4.3)) | ||||
| 46 | + zA, zB = gf["zones"].geometry.iloc[0], gf["zones"].geometry.iloc[1] | ||||
| 47 | + _poly_km(ax, zB, color=palette.seq[1], alpha=0.55, lw=0, zorder=1); _poly_km(ax, zA, color=palette.seq[5], alpha=0.6, lw=0, zorder=1) | ||||
| 48 | + riv = np.asarray(gf["cours_eau"].coords) / 1000; ax.plot(riv[:, 0], riv[:, 1], color=palette.seq[8], lw=1.1, zorder=2) | ||||
| 49 | + t300 = gf["tampons"].set_index("rayon_m").geometry.loc[c38.RAYON_CORRIDOR_REF_M] | ||||
| 50 | + _poly_km(ax, t300, color=palette.cat[1], alpha=0.3, lw=0, zorder=1) | ||||
| 51 | + cor = np.asarray(gf["corridor"].coords) / 1000; ax.plot(cor[:, 0], cor[:, 1], color=palette.cat[1], lw=1.4, zorder=2) | ||||
| 52 | + hors = (v["classe_expo"] == "hors zone").values | ||||
| 53 | + ax.scatter(xy[hors, 0], xy[hors, 1], s=6, color=palette.ink2, alpha=0.45, lw=0, zorder=3) | ||||
| 54 | + mA = (v["classe_expo"] == "A").values; mB = (v["classe_expo"] == "B").values | ||||
| 55 | + ax.scatter(xy[mB, 0], xy[mB, 1], s=14, color=palette.seq[4], lw=0, zorder=4); ax.scatter(xy[mA, 0], xy[mA, 1], s=16, color=palette.seq[9], lw=0, zorder=4) | ||||
| 56 | + pc = np.column_stack([gf["poles"].geometry.x.values, gf["poles"].geometry.y.values]) / 1000 | ||||
| 57 | + ax.scatter(pc[:, 0], pc[:, 1], marker="^", s=18 + gf["poles"]["masse"].values * 1.2, color=palette.cat[2], edgecolor=palette.ink, lw=0.6, zorder=5) | ||||
| 58 | + for (x, y), nom in zip(pc, gf["poles"]["pole"]): | ||||
| 59 | + ax.annotate(nom, (x, y), xytext=(5, 4), textcoords="offset points", fontsize=7, color=palette.ink2) | ||||
| 60 | + ax.scatter(cs[0], cs[1], marker="*", s=130, color=palette.ink, zorder=6) | ||||
| 61 | + for s in SECTEURS: | ||||
| 62 | + m = (v["secteur"] == s).values; ax.text(np.median(xy[m, 0]), np.median(xy[m, 1]) + 3.2, s, fontsize=7, color=palette.ink2, ha="center", | ||||
| 63 | + bbox=dict(boxstyle="round,pad=0.15", fc=palette.surface, ec="none", alpha=0.8)) | ||||
| 64 | + _km(ax, xy); ax.set_xlabel("Est, NAD83 / MTM zone 9 (km) — coordonnées synthétiques"); ax.set_ylabel("Nord (km)") | ||||
| 65 | + leg = [Patch(color=palette.seq[5], alpha=0.6, label="Zone d'exposition fictive, classe A (≤ 350 m)"), Patch(color=palette.seq[1], alpha=0.55, label="Classe B (350–900 m)"), | ||||
| 66 | + Line2D([], [], color=palette.cat[1], lw=1.4, label="Corridor de transport fictif et tampon de 300 m"), | ||||
| 67 | + Line2D([], [], marker="^", color="none", markerfacecolor=palette.cat[2], markeredgecolor=palette.ink, markersize=7, label="Pôles fictifs (taille ∝ masse)"), | ||||
| 68 | + Line2D([], [], marker="o", color="none", markerfacecolor=palette.seq[9], markersize=5, label="Vente en classe A"), | ||||
| 69 | + Line2D([], [], marker="o", color="none", markerfacecolor=palette.seq[4], markersize=5, label="Vente en classe B"), | ||||
| 70 | + Line2D([], [], marker="o", color="none", markerfacecolor=palette.ink2, markersize=4, alpha=0.6, label="Autre vente"), | ||||
| 71 | + Line2D([], [], marker="*", color="none", markerfacecolor=palette.ink, markersize=10, label="Sujet")] | ||||
| 72 | + ax.legend(handles=leg, fontsize=6.8, loc="upper left", ncol=2, framealpha=0.9, frameon=True, edgecolor="none", facecolor=palette.surface) | ||||
| 73 | + savefig(fig, "ch38/ch38-fig-carte-cas1.pdf") | ||||
| 74 | + | ||||
| 75 | + | ||||
| 76 | +def _carte_classes(ax, cells, z, bornes, cmap_colors, titre, fmt=lambda x: fmt_milliers(x) + " $"): | ||||
| 77 | + cmap = ListedColormap(cmap_colors); normc = BoundaryNorm(bornes, cmap.N) | ||||
| 78 | + for geom, val in zip(cells.geometry, z): | ||||
| 79 | + ext = np.asarray(geom.exterior.coords) / 1000 | ||||
| 80 | + ax.fill(ext[:, 0], ext[:, 1], color=cmap(normc(val)), lw=0.3, ec=palette.surface) | ||||
| 81 | + ax.set_title(titre, fontsize=8.5) | ||||
| 82 | + handles = [Patch(color=cmap(i), label=f"{fmt(bornes[i])} – {fmt(bornes[i + 1])}") for i in range(len(bornes) - 1)] | ||||
| 83 | + ax.legend(handles=handles, fontsize=6, loc="lower right", frameon=True, framealpha=0.9, edgecolor="none", facecolor=palette.surface, title_fontsize=6.5) | ||||
| 84 | + | ||||
| 85 | + | ||||
| 86 | +def fig_choroplethe(r): | ||||
| 87 | + g = r["grille"]; ok = g["cartographiees"]; occ = g["occupees"]; b = r["base"]; xy = b["xy_m"] / 1000 | ||||
| 88 | + z = ok["prix_m2_med"].values; k = c38.K_CLASSES | ||||
| 89 | + seq = [palette.seq[i] for i in np.linspace(0, 9, k).round().astype(int)] | ||||
| 90 | + fig, axes = plt.subplots(2, 2, figsize=(LARGEUR, 6.4), sharex=True, sharey=True) | ||||
| 91 | + titres = {"intervalles égaux": "(a) Intervalles égaux", "quantiles": "(b) Quantiles", "Jenks": "(c) Ruptures naturelles (Fisher–Jenks)"} | ||||
| 92 | + for ax, (nom, cl) in zip(axes.flat[:3], g["classif"].items()): | ||||
| 93 | + for geom in occ.loc[occ["n"] < c38.MIN_VENTES_CELLULE].geometry: | ||||
| 94 | + ext = np.asarray(geom.exterior.coords) / 1000; ax.fill(ext[:, 0], ext[:, 1], facecolor="none", edgecolor=palette.ink2, hatch="////", lw=0.3) | ||||
| 95 | + _carte_classes(ax, ok, z, cl["bornes"], seq, f"{titres[nom]} — GVF = {cl['gvf']:.3f}") | ||||
| 96 | + ax.scatter(xy[:, 0], xy[:, 1], s=1.2, color=palette.ink, alpha=0.25, lw=0) | ||||
| 97 | + ax = axes.flat[3]; res = ok["residu_moy"].values; vmax = float(np.abs(res).max()); bornes = np.linspace(-vmax, vmax, 7) | ||||
| 98 | + div = [palette.div["neg"], "#7fa8e0", "#c5d4ea", "#ead0cf", "#eb9a99", palette.div["pos"]] | ||||
| 99 | + for geom in occ.loc[occ["n"] < c38.MIN_VENTES_CELLULE].geometry: | ||||
| 100 | + ext = np.asarray(geom.exterior.coords) / 1000; ax.fill(ext[:, 0], ext[:, 1], facecolor="none", edgecolor=palette.ink2, hatch="////", lw=0.3) | ||||
| 101 | + _carte_classes(ax, ok, res, bornes, div, "(d) Résidu moyen du modèle de référence (divergent)", fmt=lambda x: f"{x:+.3f}".replace(".", ",")) | ||||
| 102 | + ax.scatter(xy[:, 0], xy[:, 1], s=1.2, color=palette.ink, alpha=0.25, lw=0) | ||||
| 103 | + for ax in axes.flat: | ||||
| 104 | + _km(ax, xy, marge=1.0) | ||||
| 105 | + for ax in axes[1]: | ||||
| 106 | + ax.set_xlabel("Est, MTM 9 (km)") | ||||
| 107 | + for ax in axes[:, 0]: | ||||
| 108 | + ax.set_ylabel("Nord (km)") | ||||
| 109 | + fig.text(0.5, 0.005, f"Cellules hexagonales de {g['aire_cellule_km2']:.2f} km² ; hachures : cellules de moins de {c38.MIN_VENTES_CELLULE} ventes (non cartographiées) ; points : ventes", ha="center", fontsize=7, color=palette.ink2) | ||||
| 110 | + fig.subplots_adjust(hspace=0.18, wspace=0.08) | ||||
| 111 | + savefig(fig, "ch38/ch38-fig-choroplethe.pdf") | ||||
| 112 | + | ||||
| 113 | + | ||||
| 114 | +def fig_faux_positifs(r): | ||||
| 115 | + fp = r["faux_positifs"]; tab = fp["tableau"]; T = fp["t"] | ||||
| 116 | + fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.1)) | ||||
| 117 | + ax = axes[0]; bins = np.linspace(-8, 8, 49); x = np.linspace(-8, 8, 400) | ||||
| 118 | + for i, (spec, col) in enumerate([("référence", palette.cat[0]), ("sans secteur", palette.cat[1])]): | ||||
| 119 | + t = np.array(T[spec][1.0])[:, 0] | ||||
| 120 | + ax.hist(t, bins=bins, density=True, histtype="stepfilled", alpha=0.45, color=col, lw=0, label=f"{spec} (écart-type {t.std():.2f})") | ||||
| 121 | + ax.plot(x, norm.pdf(x), color=palette.ink, lw=1.0, ls="--", label="Loi normale centrée réduite") | ||||
| 122 | + for s in (-1.96, 1.96): | ||||
| 123 | + ax.axvline(s, color=palette.ink2, lw=0.7, ls=":") | ||||
| 124 | + ax.set_xlabel("Statistique $t$ de l'indicatrice d'un disque aléatoire de 1 km"); ax.set_ylabel("Densité"); ax.set_title("(a) Distribution des $t$ (346 disques)", fontsize=8.5) | ||||
| 125 | + ax.legend(fontsize=6.8, loc="upper left") | ||||
| 126 | + ax = axes[1] | ||||
| 127 | + for i, spec in enumerate(["référence", "sans secteur", "sans secteur ni gradient"]): | ||||
| 128 | + w = tab[(tab["spec"] == spec) & (tab["erreur_type"] == "classique")].sort_values("rayon_km") | ||||
| 129 | + ax.plot(w["rayon_km"], w["part_5pct"] * 100, marker="o", color=palette.cat[i], label=spec) | ||||
| 130 | + ax.axhline(5, color=palette.ink2, lw=0.8, ls="--"); ax.text(2.0, 6.0, "seuil nominal 5 %", fontsize=7, color=palette.ink2, ha="right") | ||||
| 131 | + ax.set_xlabel("Rayon du disque aléatoire (km)"); ax.set_ylabel("Part des coefficients « significatifs » à 5 % (%)"); ax.set_ylim(0, 85) | ||||
| 132 | + ax.set_xticks([0.5, 1.0, 2.0]); ax.set_title("(b) Faux positifs selon la spécification", fontsize=8.5); ax.legend(fontsize=6.8, loc="upper left") | ||||
| 133 | + fig.tight_layout(); savefig(fig, "ch38/ch38-fig-faux-positifs.pdf") | ||||
| 134 | + | ||||
| 135 | + | ||||
| 136 | +def fig_tampon(r): | ||||
| 137 | + tab = r["tampon"] | ||||
| 138 | + fig, ax = plt.subplots(figsize=(LARGEUR, 3.0)) | ||||
| 139 | + for i, (spec, dx) in enumerate([("référence", -0.012), ("sans secteur", 0.012)]): | ||||
| 140 | + w = tab[tab["spec"] == spec].sort_values("rayon_m"); x = np.log10(w["rayon_m"]) + dx | ||||
| 141 | + ax.errorbar(x, w["coef"] * 100, yerr=1.96 * w["se"] * 100, fmt="o-", color=palette.cat[i], capsize=2, lw=1.2, ms=4, label=spec) | ||||
| 142 | + ax.axhline(0, color=palette.ink, lw=0.8) | ||||
| 143 | + w = tab[tab["spec"] == "référence"].sort_values("rayon_m") | ||||
| 144 | + for xr, n in zip(w["rayon_m"], w["n_dans_tampon"]): | ||||
| 145 | + ax.text(np.log10(xr), 13.5, f"n = {n}", fontsize=6.5, ha="center", color=palette.ink2) | ||||
| 146 | + ax.set_xticks(np.log10(w["rayon_m"])); ax.set_xticklabels([str(int(v)) for v in w["rayon_m"]]) | ||||
| 147 | + ax.set_xlabel("Rayon du tampon autour du corridor fictif (m, échelle logarithmique)"); ax.set_ylabel("Coefficient de l'indicatrice (%) et IC à 95 %") | ||||
| 148 | + ax.set_ylim(-20, 16); ax.legend(fontsize=7.5, loc="lower right") | ||||
| 149 | + savefig(fig, "ch38/ch38-fig-tampon.pdf") | ||||
| 150 | + | ||||
| 151 | + | ||||
| 152 | +def fig_cas4(r): | ||||
| 153 | + c4 = r["cas4"]; g = c4["gdf"]; par = c4["par_rmr"] | ||||
| 154 | + fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.6), gridspec_kw={"width_ratios": [1.55, 1]}) | ||||
| 155 | + ax = axes[0]; x = g.geometry.x.values / 1000; y = g.geometry.y.values / 1000; val = g["valeur_2026"].values | ||||
| 156 | + s = 12 + 110 * val / val.max() | ||||
| 157 | + ax.scatter(x, y, s=s, color=palette.cat[0], alpha=0.55, lw=0.4, edgecolor=palette.ink, zorder=3) | ||||
| 158 | + m = g["zone_inondable"].values == 1 | ||||
| 159 | + ax.scatter(x[m], y[m], s=s[m] + 60, facecolor="none", edgecolor=palette.cat[5], lw=1.6, zorder=4) | ||||
| 160 | + m2 = g["contamination_sol"].values == 1 | ||||
| 161 | + ax.scatter(x[m2], y[m2], marker="s", s=s[m2] + 70, facecolor="none", edgecolor=palette.cat[3], lw=1.3, zorder=4) | ||||
| 162 | + for rmr, w in g.groupby("rmr"): | ||||
| 163 | + cx, cy = w.geometry.x.mean() / 1000, w.geometry.y.mean() / 1000 | ||||
| 164 | + ax.annotate(f"{rmr} ({len(w)})", (cx, cy), xytext=(0, 16 if rmr != "Montréal" else -22), textcoords="offset points", fontsize=7, ha="center", color=palette.ink2) | ||||
| 165 | + ax.set_aspect("equal"); ax.grid(True, lw=0.4); ax.set_xlabel("Est, NAD83 / Québec Lambert (km)"); ax.set_ylabel("Nord (km)") | ||||
| 166 | + ax.set_title("(a) 48 immeubles, taille ∝ valeur 2026", fontsize=8.5) | ||||
| 167 | + leg = [Line2D([], [], marker="o", color="none", markerfacecolor=palette.cat[0], alpha=0.6, markersize=7, label="Immeuble"), | ||||
| 168 | + Line2D([], [], marker="o", color="none", markerfacecolor="none", markeredgecolor=palette.cat[5], markeredgewidth=1.6, markersize=9, label="Zone inondable (3)"), | ||||
| 169 | + Line2D([], [], marker="s", color="none", markerfacecolor="none", markeredgecolor=palette.cat[3], markeredgewidth=1.3, markersize=9, label="Contamination du sol (2)")] | ||||
| 170 | + ax.legend(handles=leg, fontsize=6.8, loc="lower right") | ||||
| 171 | + ax = axes[1]; ordre = par.index[::-1]; yy = np.arange(len(ordre)) | ||||
| 172 | + ax.barh(yy, par.loc[ordre, "part_inond_rmr"] * 100, color=palette.cat[5], alpha=0.85, height=0.6) | ||||
| 173 | + for yi, rmr in zip(yy, ordre): | ||||
| 174 | + p = par.loc[rmr, "part_inond_rmr"] * 100; ax.text(p + 0.8, yi, f"{p:.1f} %".replace(".", ","), va="center", fontsize=7) | ||||
| 175 | + ax.set_yticks(yy); ax.set_yticklabels(ordre, fontsize=7.5); ax.set_xlim(0, 36); ax.set_xlabel("Part de la valeur de la RMR en zone inondable (%)") | ||||
| 176 | + ax.set_title(f"(b) Ensemble : {c4['total']['part_zone_inondable'] * 100:.1f} % de la valeur".replace(".", ","), fontsize=8.5) | ||||
| 177 | + fig.tight_layout(); savefig(fig, "ch38/ch38-fig-cas4.pdf") | ||||
| 178 | + | ||||
| 179 | + | ||||
| 180 | +def main(): | ||||
| 181 | + warnings.filterwarnings("ignore"); setup(); r = c38.tout() | ||||
| 182 | + fig_carte_cas1(r); fig_choroplethe(r); fig_faux_positifs(r); fig_tampon(r); fig_cas4(r) | ||||
| 183 | + print("figures ch38 écrites dans", ROOT / "assets/figures/ch38") | ||||
| 184 | + | ||||
| 185 | + | ||||
| 186 | +if __name__ == "__main__": | ||||
| 187 | + main() | ||||
modified
assets/figures/ch36/ch36-fig-bacon.pdf
+0 −0
Binary file not shown.
modified
assets/figures/ch36/ch36-fig-discontinuite.pdf
+0 −0
Binary file not shown.
modified
assets/figures/ch36/ch36-fig-event-study.pdf
+0 −0
Binary file not shown.
modified
assets/figures/ch36/ch36-fig-sensibilite.pdf
+0 −0
Binary file not shown.
modified
assets/figures/ch36/ch36-fig-synthetique.pdf
+0 −0
Binary file not shown.
modified
assets/figures/ch36/ch36-fig-trajectoires.pdf
+0 −0
Binary file not shown.
added
assets/figures/ch38/ch38-fig-carte-cas1.pdf
+0 −0
Binary file not shown.
added
assets/figures/ch38/ch38-fig-cas4.pdf
+0 −0
Binary file not shown.
added
assets/figures/ch38/ch38-fig-choroplethe.pdf
+0 −0
Binary file not shown.
added
assets/figures/ch38/ch38-fig-faux-positifs.pdf
+0 −0
Binary file not shown.
added
assets/figures/ch38/ch38-fig-tampon.pdf
+0 −0
Binary file not shown.
modified
manuscript/chapters/ch30.tex
+4 −4
@@ -114,9 +114,9 @@ Le cas 3 retient $v = \pct{3.9}$. Ce taux est supérieur à l'inoccupation de \p | |||||
| 114 | 114 | ||||
| 115 | \subsection{Dépenses, réserve de remplacement et dépenses en capital} | 115 | \subsection{Dépenses, réserve de remplacement et dépenses en capital} | ||
| 116 | 116 | ||||
| 117 | −Les dépenses d'exploitation normalisées de l'année 1 proviennent du \cref{ch:28} : dernière année observée indexée, taxes recalculées sur la valeur au rôle et les taux de 2026, gestion imputée à \pct{4} du revenu brut effectif, réserve de remplacement de \dollars{650} par logement. Pour les années suivantes, les dépenses hors réserve croissent au taux $g_d = \pct{3}$, proche de l'inflation observée ; la réserve, qui est une convention de calcul et non une facture, est maintenue constante ; les autres revenus croissent comme les dépenses. | 117 | +Les dépenses d'exploitation normalisées de l'année 1 proviennent du \cref{ch:28} : dernière année observée indexée, taxes recalculées sur la valeur au rôle et les taux de 2026, gestion imputée à \pct{4} du revenu brut effectif, réserve de remplacement de \dollars{675} par logement (montant historique du dossier, \dollars{8103}). Pour les années suivantes, les dépenses hors réserve croissent au taux $g_d = \pct{3}$, proche de l'inflation observée ; la réserve, qui est une convention de calcul et non une facture, est maintenue constante ; les autres revenus croissent comme les dépenses. | ||
| 118 | 118 | ||||
| 119 | −La question délicate est celle des \motcle{dépenses en capital} (CAPEX). La convention de l'ouvrage, énoncée à l'\cref{app:conventions-calcul}, est d'inclure la réserve dans le revenu net d'exploitation servant à la capitalisation directe et de traiter explicitement, dans les flux actualisés, les dépenses en capital ponctuelles prévues : toiture, fenêtres, balcons, système de chauffage. La réserve couvre le renouvellement régulier des composantes à vie courte (électroménagers, revêtements, peinture) qu'un acheteur intègre dans son budget annuel ; les CAPEX couvrent les travaux lourds, discontinus, dont la date et le montant sont connus au moment de l'évaluation. Les deux ne doivent jamais couvrir les mêmes éléments : si la réserve est calibrée pour inclure la toiture, celle-ci ne doit pas réapparaître en CAPEX, et réciproquement. Le cas 3 programme \dollars{85000} en 2027 (maçonnerie et balcons), \dollars{40000} en 2029 (fenêtres) et \dollars{120000} en 2032 (toiture et chauffage) ; ces montants sont des hypothèses du dossier synthétique, et leur cohérence avec la réserve est assurée par construction, la réserve de \dollars{650} par logement ne visant que les composantes intérieures. Le revenu net et le flux avant dette de chaque année s'écrivent alors | 119 | +La question délicate est celle des \motcle{dépenses en capital} (CAPEX). La convention de l'ouvrage, énoncée à l'\cref{app:conventions-calcul}, est d'inclure la réserve dans le revenu net d'exploitation servant à la capitalisation directe et de traiter explicitement, dans les flux actualisés, les dépenses en capital ponctuelles prévues : toiture, fenêtres, balcons, système de chauffage. La réserve couvre le renouvellement régulier des composantes à vie courte (électroménagers, revêtements, peinture) qu'un acheteur intègre dans son budget annuel ; les CAPEX couvrent les travaux lourds, discontinus, dont la date et le montant sont connus au moment de l'évaluation. Les deux ne doivent jamais couvrir les mêmes éléments : si la réserve est calibrée pour inclure la toiture, celle-ci ne doit pas réapparaître en CAPEX, et réciproquement. Le cas 3 programme \dollars{85000} en 2027 (maçonnerie et balcons), \dollars{40000} en 2029 (fenêtres) et \dollars{120000} en 2032 (toiture et chauffage) ; ces montants sont des hypothèses du dossier synthétique, et leur cohérence avec la réserve est assurée par construction, la réserve de \dollars{675} par logement ne visant que les composantes intérieures. Le revenu net et le flux avant dette de chaque année s'écrivent alors | ||
| 120 | \begin{equation}\label{eq:30:flux} | 120 | \begin{equation}\label{eq:30:flux} | ||
| 121 | \NOI_t = RBE_t - DE_t - RR_t, \qquad CF_t = \NOI_t - CAPEX_t, | 121 | \NOI_t = RBE_t - DE_t - RR_t, \qquad CF_t = \NOI_t - CAPEX_t, | ||
| 122 | \end{equation} | 122 | \end{equation} | ||
@@ -147,7 +147,7 @@ Croissance des loyers de marché $g_m$ & \pct{3.5} par an & IPC Québec \pct{3.1 | |||||
| 147 | Rotation annuelle $\rho$ & \pct{15} des logements, les plus sous le marché d'abord & SCHL : rotation RMR Montréal \pct{9.4} (2025) \\ | 147 | Rotation annuelle $\rho$ & \pct{15} des logements, les plus sous le marché d'abord & SCHL : rotation RMR Montréal \pct{9.4} (2025) \\ | ||
| 148 | Inoccupation et pertes $v$ & \pct{3.9} du RBP & SCHL : inocc. \pct{2.9} (2025), prévision \pct{3.8}--\pct{4.8} \\ | 148 | Inoccupation et pertes $v$ & \pct{3.9} du RBP & SCHL : inocc. \pct{2.9} (2025), prévision \pct{3.8}--\pct{4.8} \\ | ||
| 149 | Croissance des dépenses $g_d$ & \pct{3.0} par an & IPC Québec \pct{3.1} \\ | 149 | Croissance des dépenses $g_d$ & \pct{3.0} par an & IPC Québec \pct{3.1} \\ | ||
| 150 | −Réserve de remplacement & \dollars{650}/logement, constante & convention de l'ouvrage (\cref{ch:28}) \\ | 150 | +Réserve de remplacement & \dollars{675}/logement (\dollars{8103}), constante & convention de l'ouvrage (\cref{ch:28}) \\ | ||
| 151 | CAPEX programmés & \dollars{85000} (2027), \dollars{40000} (2029), \dollars{120000} (2032) & hypothèse du dossier \\ | 151 | CAPEX programmés & \dollars{85000} (2027), \dollars{40000} (2029), \dollars{120000} (2032) & hypothèse du dossier \\ | ||
| 152 | Financement (\cref{ch:31}) & RPV max \pct{75}, RCD min $1{,}20$, \pct{4.12} nominal semestriel, 480 mois & BdC, fonds avancés assurés $\geq$ 5 ans, août 2026 \\ | 152 | Financement (\cref{ch:31}) & RPV max \pct{75}, RCD min $1{,}20$, \pct{4.12} nominal semestriel, 480 mois & BdC, fonds avancés assurés $\geq$ 5 ans, août 2026 \\ | ||
| 153 | \bottomrule | 153 | \bottomrule | ||
@@ -269,7 +269,7 @@ Tous les chiffres de cette section sont calculés par \texttt{analysis/python/ca | |||||
| 269 | 269 | ||||
| 270 | \subsection{Point de départ : la première année} | 270 | \subsection{Point de départ : la première année} | ||
| 271 | 271 | ||||
| 272 | −Le revenu brut potentiel contractuel des douze baux est de \dollars{165540} ; au loyer de marché, il serait de \dollars{186036}, soit un écart de \dollars{20496} (\pct{11.0}). La provision pour inoccupation et pertes de \pct{3.9} retranche \dollars{6456} ; les autres revenus (stationnement, buanderie) ajoutent \dollars{9360} ; le revenu brut effectif est de \dollars{168444}. Les dépenses d'exploitation normalisées hors réserve s'élèvent à \dollars{70385} et la réserve de remplacement à \dollars{8103} (douze fois \dollars{650}, arrondie dans le dossier), pour un ratio de dépenses de \pct{46.6} du revenu brut effectif ; le revenu net d'exploitation de la première année est de \dollars{89956}, soit \dollars{7496} par logement. Le taux global extrait des quarante ventes comparables a une médiane de \pct{4.415}, une moyenne de \pct{4.404}, une moyenne pondérée par le revenu net de \pct{4.332} et un intervalle interquartile de \pct{4.152} à \pct{4.638} ; la capitalisation directe au taux médian donne \dollars{2037503}. | 272 | +Le revenu brut potentiel contractuel des douze baux est de \dollars{165540} ; au loyer de marché, il serait de \dollars{186036}, soit un écart de \dollars{20496} (\pct{11.0}). La provision pour inoccupation et pertes de \pct{3.9} retranche \dollars{6456} ; les autres revenus (stationnement, buanderie) ajoutent \dollars{9360} ; le revenu brut effectif est de \dollars{168444}. Les dépenses d'exploitation normalisées hors réserve s'élèvent à \dollars{70385} et la réserve de remplacement à \dollars{8103} (montant historique 2025 du dossier, soit \dollars{675.25} par logement), pour un ratio de dépenses de \pct{46.6} du revenu brut effectif ; le revenu net d'exploitation de la première année est de \dollars{89956}, soit \dollars{7496} par logement. Le taux global extrait des quarante ventes comparables a une médiane de \pct{4.415}, une moyenne de \pct{4.404}, une moyenne pondérée par le revenu net de \pct{4.332} et un intervalle interquartile de \pct{4.152} à \pct{4.638} ; la capitalisation directe au taux médian donne \dollars{2037503}. | ||
| 273 | 273 | ||||
| 274 | \subsection{Projection des revenus et des dépenses} | 274 | \subsection{Projection des revenus et des dépenses} | ||
| 275 | 275 | ||||
modified
verification/review_log.md
+1 −0
@@ -8,3 +8,4 @@ | |||||
| 8 | - 2026-10-31 (ch. 3) : SRC-E-019 « RMR Montréal 2025 : 21 138 mises en chantier » ne couvre que 9 mois (34-10-0143 : juillet–septembre 2025 absents) → ne pas citer comme total annuel ; V39079 Valet commence en 2009 (SOURCES.md dit 2000). | 8 | - 2026-10-31 (ch. 3) : SRC-E-019 « RMR Montréal 2025 : 21 138 mises en chantier » ne couvre que 9 mois (34-10-0143 : juillet–septembre 2025 absents) → ne pas citer comme total annuel ; V39079 Valet commence en 2009 (SOURCES.md dit 2000). | ||
| 9 | - 2026-10-31 (ch. 15) : harmonisation à faire — ch. 4 « Sherbrooke 4,4 mois » vs 4,3 (ch. 3/15) ; rotation 15 % du cas 3 vs SCHL 8,5-9,4 % (déclaré). | 9 | - 2026-10-31 (ch. 15) : harmonisation à faire — ch. 4 « Sherbrooke 4,4 mois » vs 4,3 (ch. 3/15) ; rotation 15 % du cas 3 vs SCHL 8,5-9,4 % (déclaré). | ||
| 10 | - 2026-10-31 (ch. 32) : à trancher en harmonisation — taxes du dossier cas 2 (5 900 $) vs rôle × taux 2026 (5 438 $) ; prix cas 2 indexés ≈ +15 % vs médiane T3 2026 (normalisation au 2024-07-01). | 10 | - 2026-10-31 (ch. 32) : à trancher en harmonisation — taxes du dossier cas 2 (5 900 $) vs rôle × taux 2026 (5 438 $) ; prix cas 2 indexés ≈ +15 % vs médiane T3 2026 (normalisation au 2024-07-01). | ||
| 11 | +- 2026-10-31 (ch. 28) : harmonisation — ch. 30 écrit « réserve 650 $/logement » alors que code/ch. 28/29/47 = 675,25 $/log. (8 103 $) ; nature des CAPEX 2027/2029/2032 décrite différemment ch. 30 vs ch. 47. | ||||
| 11 | 12 | ||||