spb/traite-evaluation-qc
Public
1#!/usr/bin/env python32"""Chapitre 40 — Monte-Carlo et stress tests (DONNÉES SYNTHÉTIQUES : cas 3, 12 logements à Montréal ; cas 4, portefeuille).34Ce script expose la MÉTHODE que le chapitre 47 (§47.12) APPLIQUE : il réutilise, sans les redéfinir, les lois `MC_LOIS`,5la graine 20261031 et les fonctions vectorisées `tirer`, `projeter_vectorise`, `valeur_vectorisee`, `monte_carlo` et6`tornade` du script `ch47_dossier.py`, de sorte que les chiffres des deux chapitres coïncident exactement7(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).89Il y ajoute : moments théoriques des lois et paramétrage par quantiles (élicitation), bêta-PERT, corrélation r–R_T,10copule gaussienne / Cholesky, erreur Monte-Carlo (écart-type de la moyenne et des quantiles : théorie vs répétitions),11convergence des quantiles selon N, échantillonnage par hypercube latin (variance comparée à l'échantillonnage simple),12indices de Sobol de premier ordre et totaux (estimateur par échantillons appariés), corrélations de rang, mesures de13risque (VaR, CVaR de valeur, probabilités de bris de covenant, vérification analytique), stress tests déterministes14avec probabilité attachée, stress test inversé, agrégation d'un portefeuille de 48 immeubles avec corrélation régionale15(fonction réutilisable `simuler_portefeuille`, cas 4, chapitre 48) et les exercices du chapitre.1617Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch40_montecarlo.py [--rapide]18"""19from __future__ import annotations20import sys, pathlib, argparse, math21import numpy as np, pandas as pd22from scipy import stats, optimize23sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))24sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))25from traite import cas3, cas4, finance as fin26import ch47_dossier as c4727from ch47_dossier import MC_LOIS, R_T_PLANCHER, tirer, projeter_vectorise, valeur_vectorisee, monte_carlo, tornade # noqa: F401 (réexportés)2829ROOT = pathlib.Path(__file__).resolve().parents[3]30GRAINE = c47.GRAINE # 2026103131N_MC = 10_00032Z90 = stats.norm.ppf(0.95) # 1,64493334# ----------------------------------------------------------------------------- 40.2 lois : moments théoriques, quantiles, élicitation3536def 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) / 1847 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)5051def 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)))5556def 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)6364def 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)6970def 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}7778def 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)8889# ----------------------------------------------------------------------------- 40.2 corrélations : Cholesky, copule gaussienne, lien r–R_T9091def 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)]])9495def 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.T99 U = stats.norm.cdf(Z)100 return np.column_stack([quantile_loi(t, p, U[:, k]) for k, (t, p) in enumerate(lois)])101102def 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)105106def 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 out117118# ----------------------------------------------------------------------------- 40.3 simulation : évaluation vectorisée, erreur Monte-Carlo, convergence119120def 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"]}125126def 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}130131def 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)138139def 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)149150def 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)155156def 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 d166167def 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))171172# ----------------------------------------------------------------------------- 40.4 hypercube latin173174def 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)) / n179 return U180181def 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 x192193def 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)) if lhs[s].var(ddof=1) > 0 else float("inf")})200 return pd.DataFrame(rows)201202# ----------------------------------------------------------------------------- 40.5 mesures de risque, vérification analytique, stress tests203204def 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)211212def 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) / alpha215216def 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) / SD223 v_seuil = lambda cible: 1 - ((cible * SD + de_hg + RR) / 0.96 - A) / rbp224 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)}229230def 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])240241def 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)256257def 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_000261 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}267268# ----------------------------------------------------------------------------- 40.6 sensibilité globale : Sobol (échantillons appariés) et corrélations de rang269270NOMS_SOBOL = ["inoccupation_pertes", "croissance_loyers_marche", "rotation_annuelle", "croissance_depenses", "taux_actualisation", "ecart_terminal", "facteur_capex"]271272def _completer(x: dict) -> dict:273 x = dict(x); x["taux_terminal"] = np.maximum(x["taux_actualisation"] - x["ecart_terminal"], R_T_PLANCHER); return x274275def sobol(r: dict, n: int = 2 * 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 tirages277 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 tous = np.concatenate([fA, fB]); var = tous.var(ddof=1); mu = tous.mean(); 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 - mu) * (fAB - fA)) / var); ST = float(0.5 * np.mean((fA - fAB) ** 2) / var) # centrage : même espérance, variance réduite286 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 d289290def 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"] ** 2297 return d298299def 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]) @ beta306 return {"beta": dict(zip(NOMS_SOBOL, beta[1:])), "R2": float(1 - ((Vs - pred) ** 2).sum() / (Vs ** 2).sum())}307308# ----------------------------------------------------------------------------- 40.8 agrégation : portefeuille de 48 immeubles avec corrélation régionale (cas 4 → ch. 48)309310def 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 = 480, 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.314315 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 service320 de la dette étant celui de la dette `dette` au taux `taux_dette` (composition semestrielle, amortissement `amortissement_mois`) ou, si321 `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 <= 1326 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"].values331 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) * eps334 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 / R338 # service de la dette : amortissement résiduel non observé → amortissement `amortissement_mois` (hypothèse), renouvellement sous stress si échéance proche339 sd = np.array([fin.Hypotheque(float(d), float(t) + (choc_taux_renouvellement if int(ech) <= horizon_renouvellement else 0.0), amortissement_mois).service_dette_annuel340 for d, t, ech in zip(pf["dette"], pf["taux_dette"], pf["echeance_dette"])])341 rcd = rne / sd342 sd0 = np.array([fin.Hypotheque(float(d), float(t), amortissement_mois).service_dette_annuel for d, t in zip(pf["dette"], pf["taux_dette"])])343 rcd_det = pf["rne_2026"].values / sd0344 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())345 q = lambda a, p: float(np.quantile(a, p))346 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)),347 "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()),348 "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()),349 "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()),350 "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()),351 "nb_renouvellements_stress": int((pf["echeance_dette"] <= horizon_renouvellement).sum()),352 "nb_rcd_det_sous_seuil": int((rcd_det < seuil_rcd).sum()), "nb_rcd_det_sous_1": int((rcd_det < 1.0).sum()), "rcd_portefeuille_det": float(pf["rne_2026"].sum() / sd0.sum()),353 "nb_rcd_sous_seuil_q95": float(np.quantile((rcd < seuil_rcd).sum(1), 0.95)), "nb_rcd_sous_1_q95": float(np.quantile((rcd < 1.0).sum(1), 0.95)),354 "P_plus_de_10_sous_1": float(((rcd < 1.0).sum(1) > 10).mean()), "rpv_portefeuille": float(pf["dette"].sum() / V0)}355 stats_["ratio_diversification"] = stats_["ecart_type"] / stats_["somme_ecarts_types"]356 return {"valeurs": V, "valeur_portefeuille": VP, "rne": rne, "rcd": rcd, "service_dette": sd, "stats": stats_,357 "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)]})}358359def comparer_structures_correlation(pf: pd.DataFrame, croissance: pd.DataFrame | None, n: int = N_MC, **kw) -> pd.DataFrame:360 """Trois structures : indépendance (0, 0), corrélation régionale (0,30 ; 0,60), corrélation parfaite (1, 1). kw → simuler_portefeuille361 (p. ex. croissance=None, sigma_g_defaut=0.10 pour un choc de revenu de type récession)."""362 rows = []363 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)]:364 s = simuler_portefeuille(pf, n, GRAINE, rho_rmr=rr, rho_prov=rp, croissance=croissance, **kw)["stats"]365 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_seuil_q95", "nb_rcd_sous_1_moyen", "nb_rcd_sous_1_q95", "P_plus_de_10_sous_1", "rcd_portefeuille_mediane"]}})366 return pd.DataFrame(rows)367368def variance_somme_correlee(sigmas: np.ndarray, rho: float) -> float:369 """É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)."""370 s2 = (sigmas ** 2).sum(); tot = sigmas.sum() ** 2371 return math.sqrt(s2 + rho * (tot - s2))372373# ----------------------------------------------------------------------------- exercices374375def exercices(r: dict, mc: dict) -> dict:376 out = {}377 # Ex. 1 : élicitation — l'expert donne q05 = 2,4 %, q95 = 6,3 % d'inoccupation → log-normale ; moyenne, P(v > 5 %)378 p = parametres_depuis_quantiles(0.024, 0.063, "lognormale"); m = moments_loi("lognormale", p)379 out["ex1"] = {**p, **m, "P_sup_5pct": float(1 - stats.norm.cdf(math.log(0.05 / p["mediane"]) / p["sigma_log"]))}380 # Ex. 2 : PERT vs triangulaire pour la rotation (0,08 ; 0,15 ; 0,25)381 t = moments_loi("triangulaire", MC_LOIS["rotation_annuelle"][1]); pe = pert(0.08, 0.15, 0.25)382 a_, m_, b_ = 0.08, 0.15, 0.25383 out["ex2"] = {"tri": t, "pert": {k: pe[k] for k in ["alpha", "beta", "moyenne", "ecart_type", "q05", "q95"]},384 "P_tri_sup_20": (b_ - 0.20) ** 2 / ((b_ - a_) * (b_ - m_)), "P_pert_sup_20": float(pe["loi"].sf(0.20))}385 # Ex. 3 : nombre de tirages pour une médiane à ±5 000 $ (95 %) et pour P(RCD<1,20) à ±1 point386 V = mc["valeur"]; e = erreur_type_quantile(V, 0.5)387 se1_med = e["se_theorique"] * math.sqrt(len(V)) # écart-type « pour N = 1 »388 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),389 "N_proba_1pt": n_requis(0.5, 0.01), "N_proba_0.5pt": n_requis(0.5, 0.005)}390 # Ex. 4 : Cholesky 2×2 pour (r, R_T) avec rho = 0,78 et sigma_RT = 0,0064 ; vérification Sigma = L L' ; Spearman induit391 L = cholesky_2x2(0.005, 0.0064, 0.78); eps = np.array([1.2, -0.5]); z = L @ eps392 out["ex4"] = {"L": L, "Sigma": L @ L.T, "spearman": spearman_depuis_pearson_gaussien(0.78), "r_tire": 0.0725 + z[0], "RT_tire": 0.05 + z[1]}393 # Ex. 5 : VaR/CVaR 10 % de valeur et loi normale de même moyenne/écart-type ; comparaison394 st = mc["stats"]; mr = mesures_risque(mc, r["dcf"]["valeur"], alphas=(0.10,)).iloc[0].to_dict()395 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)}396 # Ex. 6 : simulation avec sigma_r doublé (100 pb) : écart-type, IC 90 %, P(V<2 M$)397 lois = dict(MC_LOIS); lois["taux_actualisation"] = ("normale", {"mu": 0.0725, "sigma": 0.010})398 x = tirer(N_MC, GRAINE, len(r["hyp"]["capex"]), lois); ev = evaluer(r, x); Vx = ev["valeur"]399 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()),400 "P_V_2M_base": float((V < 2e6).mean()), "cv": float(Vx.std(ddof=1) / Vx.mean()), "ratio_ecart_type": float(Vx.std(ddof=1) / V.std(ddof=1))}401 S_r = 0.528; out["ex6"]["S_r_approx"] = 4 * S_r / (1 + 3 * S_r) # modèle additif : la variance due à r quadruple402 return out403404# ----------------------------------------------------------------------------- tout405406def tout(rapide: bool = False) -> dict:407 r = cas3.tout()408 mc = monte_carlo(r, N_MC, GRAINE); torn = tornade(r)409 x = mc["tirages"]410 D = {**r, "mc": mc, "tornade": torn, "lois": tableau_lois(x), "corr_rRT": correlation_r_RT(x), "erreurs": erreurs_theoriques(mc),411 "risque": mesures_risque(mc, r["dcf"]["valeur"]), "rcd_analytique": rcd1_analytique(r), "stress": stress_deterministes(r, mc), "inverse": stress_inverse(r),412 "sobol": sobol(r), "rangs": correlations_rang(mc), "reg": regression_standardisee(mc), "pert_rotation": pert(0.08, 0.15, 0.25)}413 k_rep = 20 if rapide else 100; k_se = 10 if rapide else 40; k_med = 40 if rapide else 300414 D["repetitions"] = repetitions(r, N_MC, k_med, graine0=1)415 rep = D["repetitions"]; D["graine_traite"] = {"mediane": mc["stats"]["mediane"], "moyenne_des_medianes": float(rep["mediane"].mean()), "se_mediane_empirique": float(rep["mediane"].std(ddof=1)),416 "z_mediane": float((mc["stats"]["mediane"] - rep["mediane"].mean()) / rep["mediane"].std(ddof=1)),417 "z_q95": float((mc["stats"]["q95"] - rep["q95"].mean()) / rep["q95"].std(ddof=1)), "z_q05": float((mc["stats"]["q05"] - rep["q05"].mean()) / rep["q05"].std(ddof=1)),418 "rang_mediane": int((rep["mediane"] < mc["stats"]["mediane"]).sum()), "k": len(rep),419 "moyenne_tiree_r": float(x["taux_actualisation"].mean()), "moyenne_tiree_gm": float(x["croissance_loyers_marche"].mean())}420 D["convergence"] = convergence(r, N_max=20_000 if rapide else 100_000)421 D["se_N"] = se_empirique_selon_N(r, k=k_se)422 D["lhs"] = comparer_lhs(r, n=1000, k=k_rep)423 c4 = cas4.charger(); cr = cas4.croissance_rne(c4["series"])424 D["portefeuille"] = simuler_portefeuille(c4["portefeuille"], N_MC, GRAINE, croissance=cr); D["portefeuille_structures"] = comparer_structures_correlation(c4["portefeuille"], cr)425 D["portefeuille_recession"] = comparer_structures_correlation(c4["portefeuille"], None, sigma_g_defaut=0.10) # choc de revenu sigma_g = 10 %, g = 0426 D["cas4"] = c4; D["cas4_croissance"] = cr427 D["exercices"] = exercices(r, mc)428 return D429430def _d(x): return f"{x:,.0f}".replace(",", " ") + " $"431def _p(x, d=2): return f"{100*x:.{d}f}".replace(".", ",") + " %"432433def imprimer(D: dict) -> None:434 pd.set_option("display.width", 200); pd.set_option("display.max_columns", 20); pd.set_option("display.float_format", lambda v: f"{v:,.6g}")435 mc, st, pr = D["mc"], D["mc"]["stats"], D["mc"]["probabilites"]436 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"]))437 print("\n=== 40.2 Lois déclarées : moments théoriques et tirés\n", D["lois"].to_string(index=False))438 print("\nPERT rotation (0,08 ; 0,15 ; 0,25) :", {k: round(v, 4) for k, v in D["pert_rotation"].items() if k != "loi"})439 print("Corrélation r–R_T :", {k: round(v, 4) for k, v in D["corr_rRT"].items()})440 print("\n=== 40.3 Monte-Carlo (N = 10 000, graine 20261031) — identique au ch. 47")441 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'])}")442 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)}")443 print(f"RCD1 q05 {mc['rcd1_q05']:.3f} | RCD6 choc médiane {mc['rcd6_choc200_mediane']:.3f} q05 {mc['rcd6_choc200_q05']:.3f}")444 print("\nErreur Monte-Carlo théorique (N = 10 000) :\n", D["erreurs"].to_string(index=False))445 rep = D["repetitions"]446 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}")447 print(f" moyennes des 300 : moyenne {rep['moyenne'].mean():,.0f} | q05 {rep['q05'].mean():,.0f} | q95 {rep['q95'].mean():,.0f}")448 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}")449 print(" graine du traité (20261031) :", {k: (round(v, 4) if isinstance(v, float) else v) for k, v in D["graine_traite"].items()})450 print("\nConvergence (une suite, graine 20261031) :\n", D["convergence"].to_string(index=False))451 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))452 print("\n=== 40.4 Hypercube latin vs simple (N = 1 000, répétitions) :\n", D["lhs"].to_string(index=False))453 print("\n=== 40.5 Mesures de risque (référence : valeur DCF déterministe)\n", D["risque"].to_string(index=False))454 print("RCD1 analytique :", {k: round(v, 5) for k, v in D["rcd_analytique"].items()})455 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"])))456 print("\nStress tests déterministes :\n", D["stress"].to_string(index=False))457 print("Stress inversé :", {k: round(v, 5) for k, v in D["inverse"].items()})458 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))459 print("\nCorrélations de rang / Pearson :\n", D["rangs"].to_string(index=False))460 print("Régression standardisée : R² = %.4f ; bêtas =" % D["reg"]["R2"], {k: round(v, 3) for k, v in D["reg"]["beta"].items()})461 print("\nTornade (ch. 47) :\n", D["tornade"].to_string(index=False))462 print("\n=== 40.8 Portefeuille (cas 4, 48 immeubles) — structure régionale (0,30 ; 0,60)")463 print({k: (round(v, 4) if isinstance(v, float) else v) for k, v in D["portefeuille"]["stats"].items()})464 print(D["portefeuille"]["par_rmr"].to_string(index=False))465 print("\nTrois structures de corrélation :\n", D["portefeuille_structures"].to_string(index=False))466 print("\nTrois structures, choc de revenu sigma_g = 10 % (g = 0) :\n", D["portefeuille_recession"].to_string(index=False))467 print("\n=== Exercices")468 for k, v in D["exercices"].items(): print(k, v)469470if __name__ == "__main__":471 ap = argparse.ArgumentParser(); ap.add_argument("--rapide", action="store_true"); a = ap.parse_args()472 imprimer(tout(rapide=a.rapide))473