#!/usr/bin/env python3
"""Chapitre 40 — Monte-Carlo et stress tests (DONNÉES SYNTHÉTIQUES : cas 3, 12 logements à Montréal ; cas 4, portefeuille).

Ce script expose la MÉTHODE que le chapitre 47 (§47.12) APPLIQUE : il réutilise, sans les redéfinir, les lois `MC_LOIS`,
la graine 20261031 et les fonctions vectorisées `tirer`, `projeter_vectorise`, `valeur_vectorisee`, `monte_carlo` et
`tornade` du script `ch47_dossier.py`, de sorte que les chiffres des deux chapitres coïncident exactement
(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).

Il y ajoute : moments théoriques des lois et paramétrage par quantiles (élicitation), bêta-PERT, corrélation r–R_T,
copule gaussienne / Cholesky, erreur Monte-Carlo (écart-type de la moyenne et des quantiles : théorie vs répétitions),
convergence des quantiles selon N, échantillonnage par hypercube latin (variance comparée à l'échantillonnage simple),
indices de Sobol de premier ordre et totaux (estimateur par échantillons appariés), corrélations de rang, mesures de
risque (VaR, CVaR de valeur, probabilités de bris de covenant, vérification analytique), stress tests déterministes
avec probabilité attachée, stress test inversé, agrégation d'un portefeuille de 48 immeubles avec corrélation régionale
(fonction réutilisable `simuler_portefeuille`, cas 4, chapitre 48) et les exercices du chapitre.

Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch40_montecarlo.py [--rapide]
"""
from __future__ import annotations
import sys, pathlib, argparse, math
import numpy as np, pandas as pd
from scipy import stats, optimize
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
from traite import cas3, cas4, finance as fin
import ch47_dossier as c47
from ch47_dossier import MC_LOIS, R_T_PLANCHER, tirer, projeter_vectorise, valeur_vectorisee, monte_carlo, tornade  # noqa: F401 (réexportés)

ROOT = pathlib.Path(__file__).resolve().parents[3]
GRAINE = c47.GRAINE          # 20261031
N_MC = 10_000
Z90 = stats.norm.ppf(0.95)   # 1,6449

# ----------------------------------------------------------------------------- 40.2 lois : moments théoriques, quantiles, élicitation

def moments_loi(typ: str, p: dict) -> dict:
    """Moyenne, écart-type, quantiles 5 % et 95 % théoriques d'une loi déclarée (mêmes paramétrages que MC_LOIS)."""
    if typ == "lognormale":
        mu, s = math.log(p["mediane"]), p["sigma_log"]
        return {"moyenne": math.exp(mu + s**2 / 2), "ecart_type": math.sqrt((math.exp(s**2) - 1) * math.exp(2 * mu + s**2)),
                "q05": math.exp(mu - Z90 * s), "q95": math.exp(mu + Z90 * s), "mediane": p["mediane"]}
    if typ == "normale":
        return {"moyenne": p["mu"], "ecart_type": p["sigma"], "q05": p["mu"] - Z90 * p["sigma"], "q95": p["mu"] + Z90 * p["sigma"], "mediane": p["mu"]}
    if typ == "triangulaire":
        a, m, b = p["min"], p["mode"], p["max"]
        var = (a * a + b * b + m * m - a * b - a * m - b * m) / 18
        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),
                "mediane": quantile_triangulaire(0.5, a, m, b)}
    raise ValueError(typ)

def quantile_triangulaire(p, a, m, b):
    Fm = (m - a) / (b - a)
    p = np.asarray(p, dtype=float)
    return np.where(p <= Fm, a + np.sqrt(p * (b - a) * (m - a)), b - np.sqrt((1 - p) * (b - a) * (b - m)))

def quantile_loi(typ: str, p: dict, u):
    """Fonction quantile F^{-1}(u) d'une loi déclarée (sert à l'hypercube latin et à la copule)."""
    u = np.asarray(u, dtype=float)
    if typ == "lognormale": return np.exp(math.log(p["mediane"]) + p["sigma_log"] * stats.norm.ppf(u))
    if typ == "normale": return p["mu"] + p["sigma"] * stats.norm.ppf(u)
    if typ == "triangulaire": return quantile_triangulaire(u, p["min"], p["mode"], p["max"])
    raise ValueError(typ)

def parametres_depuis_quantiles(q05: float, q95: float, loi: str = "normale") -> dict:
    """Élicitation : un expert donne les quantiles 5 % et 95 % ; on en déduit les paramètres (normale ou log-normale)."""
    if loi == "normale": return {"mu": (q05 + q95) / 2, "sigma": (q95 - q05) / (2 * Z90)}
    if loi == "lognormale": return {"mediane": math.sqrt(q05 * q95), "sigma_log": math.log(q95 / q05) / (2 * Z90)}
    raise ValueError(loi)

def pert(a: float, m: float, b: float, lam: float = 4.0) -> dict:
    """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)."""
    al = 1 + lam * (m - a) / (b - a); be = 1 + lam * (b - m) / (b - a)
    moy = (a + lam * m + b) / (lam + 2)
    var = (moy - a) * (b - moy) / (lam + 3)
    d = stats.beta(al, be, loc=a, scale=b - a)
    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}

def tableau_lois(x: dict | None = None) -> pd.DataFrame:
    """Moments théoriques de chaque loi de MC_LOIS et, si des tirages sont fournis, leurs statistiques empiriques."""
    rows = []
    for nom, (typ, p) in MC_LOIS.items():
        m = moments_loi(typ, p); row = {"intrant": nom, "loi": typ, **{k: m[k] for k in ["moyenne", "ecart_type", "q05", "q95"]}}
        if x is not None:
            v = x[nom] if nom != "facteur_capex" else x[nom].ravel()
            row.update({"moyenne_tiree": float(v.mean()), "q05_tire": float(np.quantile(v, 0.05)), "q95_tire": float(np.quantile(v, 0.95))})
        rows.append(row)
    return pd.DataFrame(rows)

# ----------------------------------------------------------------------------- 40.2 corrélations : Cholesky, copule gaussienne, lien r–R_T

def cholesky_2x2(s1: float, s2: float, rho: float) -> np.ndarray:
    """Facteur de Cholesky L de la matrice de covariance [[s1², rho s1 s2], [rho s1 s2, s2²]] : Sigma = L L'."""
    return np.array([[s1, 0.0], [rho * s2, s2 * math.sqrt(1 - rho**2)]])

def copule_gaussienne(n: int, graine: int, corr: np.ndarray, lois: list[tuple[str, dict]]) -> np.ndarray:
    """Tirages corrélés par copule gaussienne : Z = L eps, U = Phi(Z), X_k = F_k^{-1}(U_k). Retourne (n, d)."""
    rng = np.random.default_rng(graine)
    L = np.linalg.cholesky(corr); Z = rng.standard_normal((n, len(lois))) @ L.T
    U = stats.norm.cdf(Z)
    return np.column_stack([quantile_loi(t, p, U[:, k]) for k, (t, p) in enumerate(lois)])

def spearman_depuis_pearson_gaussien(rho: float) -> float:
    """Corrélation de rang induite par une copule gaussienne de corrélation rho : rho_S = (6/pi) arcsin(rho/2)."""
    return 6 / math.pi * math.asin(rho / 2)

def correlation_r_RT(x: dict | None = None) -> dict:
    """Corrélation théorique entre r et R_T = r − delta (sans plancher) et, si tirages fournis, corrélation empirique."""
    sr = MC_LOIS["taux_actualisation"][1]["sigma"]; sd = MC_LOIS["ecart_terminal"][1]["sigma"]
    out = {"rho_theorique": sr / math.sqrt(sr**2 + sd**2), "sigma_RT_theorique": math.sqrt(sr**2 + sd**2)}
    if x is not None:
        out["rho_empirique"] = float(np.corrcoef(x["taux_actualisation"], x["taux_terminal"])[0, 1])
        out["spearman_empirique"] = float(stats.spearmanr(x["taux_actualisation"], x["taux_terminal"])[0])
        out["part_plancher"] = float((x["taux_actualisation"] - x["ecart_terminal"] < R_T_PLANCHER).mean())
        out["sigma_RT_empirique"] = float(x["taux_terminal"].std(ddof=1))
    out["spearman_theorique"] = spearman_depuis_pearson_gaussien(out["rho_theorique"])
    return out

# ----------------------------------------------------------------------------- 40.3 simulation : évaluation vectorisée, erreur Monte-Carlo, convergence

def evaluer(r: dict, x: dict) -> dict:
    """Valeur DCF et RNE pour un dictionnaire de tirages (même structure que `tirer`)."""
    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"])
    val = valeur_vectorisee(proj, x["taux_actualisation"], x["taux_terminal"], r["hyp"]["couts_vente"])
    return {"valeur": val["valeur"], "rne": proj["rne"], "rcd1": proj["rne"][:, 0] / r["financement"]["service_dette"]}

def erreur_type_quantile(V: np.ndarray, p: float) -> dict:
    """Écart-type asymptotique d'un quantile empirique : sqrt(p(1-p)/N) / f(q_p), densité f estimée par noyau gaussien."""
    N = len(V); q = float(np.quantile(V, p)); f = float(stats.gaussian_kde(V)(q)[0])
    return {"p": p, "quantile": q, "densite": f, "se_theorique": math.sqrt(p * (1 - p) / N) / f}

def erreurs_theoriques(mc: dict) -> pd.DataFrame:
    V = mc["valeur"]; N = len(V); rows = [{"statistique": "moyenne", "estimation": float(V.mean()), "se_theorique": float(V.std(ddof=1) / math.sqrt(N))}]
    for p, nom in [(0.05, "q05"), (0.5, "mediane"), (0.95, "q95")]:
        e = erreur_type_quantile(V, p); rows.append({"statistique": nom, "estimation": e["quantile"], "se_theorique": e["se_theorique"], "densite": e["densite"]})
    rcd1 = mc["rcd"][:, 0]; p = float((rcd1 < 1.2).mean())
    rows.append({"statistique": "P(RCD1<1,20)", "estimation": p, "se_theorique": math.sqrt(p * (1 - p) / N)})
    return pd.DataFrame(rows)

def repetitions(r: dict, n: int, k: int, graine0: int = 1, lhs: bool = False) -> pd.DataFrame:
    """k simulations indépendantes de n tirages (graines graine0, graine0+1, …) ; statistiques de chaque répétition."""
    rows = []
    nb_capex = len(r["hyp"]["capex"])
    for g in range(graine0, graine0 + k):
        x = tirer_lhs(n, g, nb_capex) if lhs else tirer(n, g, nb_capex)
        e = evaluer(r, x); V = e["valeur"]
        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),
                     "p_rcd1_120": (e["rcd1"] < 1.2).mean(), "p_V_2M": (V < 2e6).mean()})
    return pd.DataFrame(rows)

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:
    """Quantiles cumulés d'une seule suite de N_max tirages (graine), lus après N tirages (N dans `points`)."""
    x = tirer(N_max, graine, len(r["hyp"]["capex"])); V = evaluer(r, x)["valeur"]
    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]
    return pd.DataFrame(rows)

def se_empirique_selon_N(r: dict, Ns=(250, 500, 1000, 2000, 4000, 8000), k: int = 40) -> pd.DataFrame:
    """É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)."""
    rows = []
    for N in Ns:
        rep = repetitions(r, N, k, graine0=100)
        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)})
    d = pd.DataFrame(rows)
    b = np.polyfit(np.log(d["N"]), np.log(d["se_mediane"]), 1)[0]
    d.attrs["pente_log_mediane"] = float(b)
    return d

def n_requis(se_unitaire: float, erreur_cible: float, confiance: float = 0.95) -> int:
    """N tel que z_(1-a/2) · se_unitaire / sqrt(N) <= erreur_cible, où se_unitaire est l'écart-type « pour N = 1 »."""
    z = stats.norm.ppf(0.5 + confiance / 2)
    return int(math.ceil((z * se_unitaire / erreur_cible) ** 2))

# ----------------------------------------------------------------------------- 40.4 hypercube latin

def hypercube_latin(n: int, d: int, rng: np.random.Generator) -> np.ndarray:
    """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."""
    U = np.empty((n, d))
    for k in range(d):
        U[:, k] = (rng.permutation(n) + rng.random(n)) / n
    return U

def tirer_lhs(n: int, graine: int, nb_capex: int, lois=MC_LOIS) -> dict:
    """Même structure que `ch47_dossier.tirer` (mêmes lois, même ordre, mêmes bornes) mais par hypercube latin et fonction quantile."""
    rng = np.random.default_rng(graine)
    noms = ["inoccupation_pertes", "croissance_loyers_marche", "rotation_annuelle", "croissance_depenses", "taux_actualisation", "ecart_terminal"]
    U = hypercube_latin(n, len(noms) + nb_capex, rng)
    x = {nom: quantile_loi(*lois[nom], U[:, k]) for k, nom in enumerate(noms)}
    x["inoccupation_pertes"] = np.clip(x["inoccupation_pertes"], 0.0, 0.30); x["taux_actualisation"] = np.clip(x["taux_actualisation"], 0.04, None)
    typ, p = lois["facteur_capex"]
    x["facteur_capex"] = np.column_stack([quantile_loi(typ, p, U[:, len(noms) + j]) for j in range(nb_capex)])
    x["taux_terminal"] = np.maximum(x["taux_actualisation"] - x["ecart_terminal"], R_T_PLANCHER)
    return x

def comparer_lhs(r: dict, n: int = 1000, k: int = 100) -> pd.DataFrame:
    """Variance des estimateurs (moyenne, médiane, q05, q95, P(RCD1<1,20)) sur k répétitions : échantillonnage simple vs hypercube latin."""
    simple = repetitions(r, n, k, graine0=1000, lhs=False); lhs = repetitions(r, n, k, graine0=1000, lhs=True)
    rows = []
    for s in ["moyenne", "mediane", "q05", "q95", "ecart_type", "p_rcd1_120", "p_V_2M"]:
        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),
                     "ratio_variance": (simple[s].var(ddof=1) / lhs[s].var(ddof=1)) if lhs[s].var(ddof=1) > 0 else float("inf")})
    return pd.DataFrame(rows)

# ----------------------------------------------------------------------------- 40.5 mesures de risque, vérification analytique, stress tests

def mesures_risque(mc: dict, V_ref: float, alphas=(0.05, 0.10)) -> pd.DataFrame:
    """VaR_alpha = V_ref − q_alpha (perte par rapport à la valeur de référence) ; CVaR_alpha = V_ref − E[V | V <= q_alpha]."""
    V = mc["valeur"]; rows = []
    for a in alphas:
        q = float(np.quantile(V, a)); cvar_v = float(V[V <= q].mean())
        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})
    return pd.DataFrame(rows)

def cvar_normale(mu: float, sigma: float, alpha: float) -> float:
    """E[V | V <= q_alpha] pour V ~ N(mu, sigma²) : mu − sigma · phi(z_alpha)/alpha."""
    z = stats.norm.ppf(alpha); return mu - sigma * stats.norm.pdf(z) / alpha

def rcd1_analytique(r: dict) -> dict:
    """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 ;
    seuils v* (RCD1 = 1,20) et v_1 (RCD1 = 1,00) et probabilités exactes sous la loi log-normale déclarée."""
    b, hyp, f = r["base"], r["hyp"], r["financement"]
    rbp, A, RR, SD = b["rbp_contrat"], hyp["autres_revenus_annuels"], b["reserve"], f["service_dette"]
    de_hg = b["depenses_hors_reserve"] - 0.04 * b["rbe"]
    rcd = lambda v: (0.96 * (rbp * (1 - v) + A) - de_hg - RR) / SD
    v_seuil = lambda cible: 1 - ((cible * SD + de_hg + RR) / 0.96 - A) / rbp
    med, s = MC_LOIS["inoccupation_pertes"][1]["mediane"], MC_LOIS["inoccupation_pertes"][1]["sigma_log"]
    P = lambda v: float(1 - stats.norm.cdf(math.log(v / med) / s))
    v120, v100 = v_seuil(1.20), v_seuil(1.00)
    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,
            "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)}

def valeur_deterministe(r: dict, **par) -> float:
    """Valeur DCF pour un jeu d'intrants déterministes (défaut : hypothèses de base) ; clés = noms de MC_LOIS (+ taux_terminal direct)."""
    hyp = r["hyp"]
    base = {"inoccupation_pertes": hyp["taux_inoccupation_pertes"], "croissance_loyers_marche": hyp["croissance_loyers_marche"], "rotation_annuelle": hyp["rotation_annuelle"],
            "croissance_depenses": hyp["croissance_depenses"], "taux_actualisation": hyp["taux_actualisation"], "ecart_terminal": hyp["taux_actualisation"] - hyp["taux_terminal"], "facteur_capex": 1.0}
    base.update({k: v for k, v in par.items() if k != "taux_terminal"})
    R_T = par.get("taux_terminal", max(base["taux_actualisation"] - base["ecart_terminal"], R_T_PLANCHER))
    x = {k: np.array([base[k]]) for k in ["inoccupation_pertes", "croissance_loyers_marche", "rotation_annuelle", "croissance_depenses", "taux_actualisation"]}
    x["facteur_capex"] = np.full((1, len(hyp["capex"])), base["facteur_capex"]); x["taux_terminal"] = np.array([R_T])
    return float(evaluer(r, x)["valeur"][0])

def stress_deterministes(r: dict, mc: dict) -> pd.DataFrame:
    """Scénarios déterministes classiques et probabilité, dans la simulation, d'un résultat au moins aussi défavorable."""
    hyp = r["hyp"]; V = mc["valeur"]; q = lambda nom, p: float(quantile_loi(*MC_LOIS[nom], p))
    sc = [("Base (hypothèses du dossier)", {}),
          ("Hausse des taux : r + 75 pb, R_T + 50 pb (ch. 30)", {"taux_actualisation": 0.08, "taux_terminal": 0.055}),
          ("Marché détendu : inoccupation 6 %", {"inoccupation_pertes": 0.06}),
          ("Croissance faible : loyers de marché 2,5 %, dépenses 4 %", {"croissance_loyers_marche": 0.025, "croissance_depenses": 0.04}),
          ("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)}),
          ("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),
                                                                   "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)})]
    rows = []
    for nom, par in sc:
        v = valeur_deterministe(r, **par)
        rows.append({"scenario": nom, "valeur": v, "ecart_pct": v / r["dcf"]["valeur"] - 1, "P_V_inferieure": float((V <= v).mean())})
    return pd.DataFrame(rows)

def stress_inverse(r: dict) -> dict:
    """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)."""
    an = rcd1_analytique(r)
    f_r = lambda rr: valeur_deterministe(r, taux_actualisation=rr, taux_terminal=max(rr - 0.0225, R_T_PLANCHER)) - 2_000_000
    r_2M = optimize.brentq(f_r, 0.05, 0.12)
    f_RT = lambda RT: valeur_deterministe(r, taux_terminal=RT) - r["valeur_capitalisation"]
    RT_cap = optimize.brentq(f_RT, 0.035, 0.09)
    f_g = lambda g: valeur_deterministe(r, croissance_loyers_marche=g) - r["valeur_capitalisation"]
    g_cap = optimize.brentq(f_g, -0.05, 0.035)
    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}

# ----------------------------------------------------------------------------- 40.6 sensibilité globale : Sobol (échantillons appariés) et corrélations de rang

NOMS_SOBOL = ["inoccupation_pertes", "croissance_loyers_marche", "rotation_annuelle", "croissance_depenses", "taux_actualisation", "ecart_terminal", "facteur_capex"]

def _completer(x: dict) -> dict:
    x = dict(x); x["taux_terminal"] = np.maximum(x["taux_actualisation"] - x["ecart_terminal"], R_T_PLANCHER); return x

def sobol(r: dict, n: int = 2 * N_MC, graine: int = GRAINE) -> pd.DataFrame:
    """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
    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.
    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)."""
    nb = len(r["hyp"]["capex"])
    A = tirer(n, graine, nb); B = tirer(n, graine + 1, nb)
    fA = evaluer(r, A)["valeur"]; fB = evaluer(r, B)["valeur"]
    tous = np.concatenate([fA, fB]); var = tous.var(ddof=1); mu = tous.mean(); rows = []
    for nom in NOMS_SOBOL:
        AB = dict(A); AB[nom] = B[nom]; fAB = evaluer(r, _completer(AB))["valeur"]
        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éduite
        rows.append({"intrant": nom, "S1": S1, "ST": ST})
    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))
    return d

def correlations_rang(mc: dict) -> pd.DataFrame:
    """Corrélations de Spearman et de Pearson entre chaque intrant et la valeur ; rho² comme approximation de S1 pour un modèle quasi additif."""
    x, V = mc["tirages"], mc["valeur"]; rows = []
    for nom in NOMS_SOBOL + ["taux_terminal"]:
        v = x[nom].mean(1) if nom == "facteur_capex" else x[nom]
        rows.append({"intrant": nom, "spearman": float(stats.spearmanr(v, V)[0]), "pearson": float(np.corrcoef(v, V)[0, 1])})
    d = pd.DataFrame(rows); d["spearman2"] = d["spearman"] ** 2; d["pearson2"] = d["pearson"] ** 2
    return d

def regression_standardisee(mc: dict) -> dict:
    """Coefficients standardisés d'une régression linéaire de V sur les sept intrants : R² mesure la part « additive » de la variance."""
    x, V = mc["tirages"], mc["valeur"]
    X = np.column_stack([x[n].mean(1) if n == "facteur_capex" else x[n] for n in NOMS_SOBOL])
    Xs = (X - X.mean(0)) / X.std(0, ddof=1); Vs = (V - V.mean()) / V.std(ddof=1)
    beta, *_ = np.linalg.lstsq(np.column_stack([np.ones(len(V)), Xs]), Vs, rcond=None)
    pred = np.column_stack([np.ones(len(V)), Xs]) @ beta
    return {"beta": dict(zip(NOMS_SOBOL, beta[1:])), "R2": float(1 - ((Vs - pred) ** 2).sum() / (Vs ** 2).sum())}

# ----------------------------------------------------------------------------- 40.8 agrégation : portefeuille de 48 immeubles avec corrélation régionale (cas 4 → ch. 48)

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,
                         croissance: pd.DataFrame | None = None, sigma_g_defaut: float = 0.03, choc_taux_renouvellement: float = 0.02, horizon_renouvellement: int = 2028,
                         amortissement_mois: int = 480, seuil_rcd: float = 1.20) -> dict:
    """Simulation à un an de la valeur et de la couverture de la dette d'un portefeuille d'immeubles à revenus avec corrélation régionale.

    Modèle (hypothèses pédagogiques déclarées au ch. 40, §40.8) : pour l'immeuble i de la RMR k,
        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,
        R_i'      = R_i + sigma_R · u_i,                u_i de même structure factorielle (facteurs indépendants de ceux du RNE),
        V_i'      = RNE_i' / R_i',
    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
    de la dette étant celui de la dette `dette` au taux `taux_dette` (composition semestrielle, amortissement `amortissement_mois`) ou, si
    `echeance_dette` <= horizon_renouvellement, au taux majoré de `choc_taux_renouvellement` (renouvellement sous stress).
    `croissance` : DataFrame (id_immeuble, croissance_annuelle, volatilite) de `traite.cas4.croissance_rne` ; sinon g = 0, sigma_g = sigma_g_defaut.
    Retourne les valeurs simulées (n, N_immeubles), la valeur du portefeuille, les RCD, et un dictionnaire de statistiques.
    Réutilisable par le chapitre 48 : `from ch40_montecarlo import simuler_portefeuille`."""
    assert 0 <= rho_prov <= rho_rmr <= 1
    rng = np.random.default_rng(graine)
    pf = pf.reset_index(drop=True); N = len(pf); rmrs = pf["rmr"].astype("category"); k = rmrs.cat.codes.values; K = len(rmrs.cat.categories)
    g = np.zeros(N); sg = np.full(N, sigma_g_defaut)
    if croissance is not None:
        c = croissance.set_index("id_immeuble").reindex(pf["id_immeuble"]); g = c["croissance_annuelle"].values; sg = c["volatilite"].values
    def facteur():
        Z0 = rng.standard_normal((n, 1)); Zk = rng.standard_normal((n, K)); eps = rng.standard_normal((n, N))
        return math.sqrt(rho_prov) * Z0 + math.sqrt(max(rho_rmr - rho_prov, 0.0)) * Zk[:, k] + math.sqrt(1 - rho_rmr) * eps
    e = facteur(); u = facteur()
    rne = pf["rne_2026"].values * np.exp(g + sg * e)
    R = np.maximum(pf["taux_global"].values + sigma_R * u, 0.02)
    V = rne / R
    # service de la dette : amortissement résiduel non observé → amortissement `amortissement_mois` (hypothèse), renouvellement sous stress si échéance proche
    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
                   for d, t, ech in zip(pf["dette"], pf["taux_dette"], pf["echeance_dette"])])
    rcd = rne / sd
    sd0 = np.array([fin.Hypotheque(float(d), float(t), amortissement_mois).service_dette_annuel for d, t in zip(pf["dette"], pf["taux_dette"])])
    rcd_det = pf["rne_2026"].values / sd0
    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())
    q = lambda a, p: float(np.quantile(a, p))
    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)),
              "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()),
              "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()),
              "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()),
              "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()),
              "nb_renouvellements_stress": int((pf["echeance_dette"] <= horizon_renouvellement).sum()),
              "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()),
              "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)),
              "P_plus_de_10_sous_1": float(((rcd < 1.0).sum(1) > 10).mean()), "rpv_portefeuille": float(pf["dette"].sum() / V0)}
    stats_["ratio_diversification"] = stats_["ecart_type"] / stats_["somme_ecarts_types"]
    return {"valeurs": V, "valeur_portefeuille": VP, "rne": rne, "rcd": rcd, "service_dette": sd, "stats": stats_,
            "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)]})}

def comparer_structures_correlation(pf: pd.DataFrame, croissance: pd.DataFrame | None, n: int = N_MC, **kw) -> pd.DataFrame:
    """Trois structures : indépendance (0, 0), corrélation régionale (0,30 ; 0,60), corrélation parfaite (1, 1). kw → simuler_portefeuille
    (p. ex. croissance=None, sigma_g_defaut=0.10 pour un choc de revenu de type récession)."""
    rows = []
    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)]:
        s = simuler_portefeuille(pf, n, GRAINE, rho_rmr=rr, rho_prov=rp, croissance=croissance, **kw)["stats"]
        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"]}})
    return pd.DataFrame(rows)

def variance_somme_correlee(sigmas: np.ndarray, rho: float) -> float:
    """É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)."""
    s2 = (sigmas ** 2).sum(); tot = sigmas.sum() ** 2
    return math.sqrt(s2 + rho * (tot - s2))

# ----------------------------------------------------------------------------- exercices

def exercices(r: dict, mc: dict) -> dict:
    out = {}
    # Ex. 1 : élicitation — l'expert donne q05 = 2,4 %, q95 = 6,3 % d'inoccupation → log-normale ; moyenne, P(v > 5 %)
    p = parametres_depuis_quantiles(0.024, 0.063, "lognormale"); m = moments_loi("lognormale", p)
    out["ex1"] = {**p, **m, "P_sup_5pct": float(1 - stats.norm.cdf(math.log(0.05 / p["mediane"]) / p["sigma_log"]))}
    # Ex. 2 : PERT vs triangulaire pour la rotation (0,08 ; 0,15 ; 0,25)
    t = moments_loi("triangulaire", MC_LOIS["rotation_annuelle"][1]); pe = pert(0.08, 0.15, 0.25)
    a_, m_, b_ = 0.08, 0.15, 0.25
    out["ex2"] = {"tri": t, "pert": {k: pe[k] for k in ["alpha", "beta", "moyenne", "ecart_type", "q05", "q95"]},
                  "P_tri_sup_20": (b_ - 0.20) ** 2 / ((b_ - a_) * (b_ - m_)), "P_pert_sup_20": float(pe["loi"].sf(0.20))}
    # Ex. 3 : nombre de tirages pour une médiane à ±5 000 $ (95 %) et pour P(RCD<1,20) à ±1 point
    V = mc["valeur"]; e = erreur_type_quantile(V, 0.5)
    se1_med = e["se_theorique"] * math.sqrt(len(V))          # écart-type « pour N = 1 »
    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),
                  "N_proba_1pt": n_requis(0.5, 0.01), "N_proba_0.5pt": n_requis(0.5, 0.005)}
    # 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
    L = cholesky_2x2(0.005, 0.0064, 0.78); eps = np.array([1.2, -0.5]); z = L @ eps
    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]}
    # Ex. 5 : VaR/CVaR 10 % de valeur et loi normale de même moyenne/écart-type ; comparaison
    st = mc["stats"]; mr = mesures_risque(mc, r["dcf"]["valeur"], alphas=(0.10,)).iloc[0].to_dict()
    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)}
    # Ex. 6 : simulation avec sigma_r doublé (100 pb) : écart-type, IC 90 %, P(V<2 M$)
    lois = dict(MC_LOIS); lois["taux_actualisation"] = ("normale", {"mu": 0.0725, "sigma": 0.010})
    x = tirer(N_MC, GRAINE, len(r["hyp"]["capex"]), lois); ev = evaluer(r, x); Vx = ev["valeur"]
    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()),
                  "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))}
    S_r = 0.528; out["ex6"]["S_r_approx"] = 4 * S_r / (1 + 3 * S_r)   # modèle additif : la variance due à r quadruple
    return out

# ----------------------------------------------------------------------------- tout

def tout(rapide: bool = False) -> dict:
    r = cas3.tout()
    mc = monte_carlo(r, N_MC, GRAINE); torn = tornade(r)
    x = mc["tirages"]
    D = {**r, "mc": mc, "tornade": torn, "lois": tableau_lois(x), "corr_rRT": correlation_r_RT(x), "erreurs": erreurs_theoriques(mc),
         "risque": mesures_risque(mc, r["dcf"]["valeur"]), "rcd_analytique": rcd1_analytique(r), "stress": stress_deterministes(r, mc), "inverse": stress_inverse(r),
         "sobol": sobol(r), "rangs": correlations_rang(mc), "reg": regression_standardisee(mc), "pert_rotation": pert(0.08, 0.15, 0.25)}
    k_rep = 20 if rapide else 100; k_se = 10 if rapide else 40; k_med = 40 if rapide else 300
    D["repetitions"] = repetitions(r, N_MC, k_med, graine0=1)
    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)),
                                                  "z_mediane": float((mc["stats"]["mediane"] - rep["mediane"].mean()) / rep["mediane"].std(ddof=1)),
                                                  "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)),
                                                  "rang_mediane": int((rep["mediane"] < mc["stats"]["mediane"]).sum()), "k": len(rep),
                                                  "moyenne_tiree_r": float(x["taux_actualisation"].mean()), "moyenne_tiree_gm": float(x["croissance_loyers_marche"].mean())}
    D["convergence"] = convergence(r, N_max=20_000 if rapide else 100_000)
    D["se_N"] = se_empirique_selon_N(r, k=k_se)
    D["lhs"] = comparer_lhs(r, n=1000, k=k_rep)
    c4 = cas4.charger(); cr = cas4.croissance_rne(c4["series"])
    D["portefeuille"] = simuler_portefeuille(c4["portefeuille"], N_MC, GRAINE, croissance=cr); D["portefeuille_structures"] = comparer_structures_correlation(c4["portefeuille"], cr)
    D["portefeuille_recession"] = comparer_structures_correlation(c4["portefeuille"], None, sigma_g_defaut=0.10)   # choc de revenu sigma_g = 10 %, g = 0
    D["cas4"] = c4; D["cas4_croissance"] = cr
    D["exercices"] = exercices(r, mc)
    return D

def _d(x): return f"{x:,.0f}".replace(",", " ") + " $"
def _p(x, d=2): return f"{100*x:.{d}f}".replace(".", ",") + " %"

def imprimer(D: dict) -> None:
    pd.set_option("display.width", 200); pd.set_option("display.max_columns", 20); pd.set_option("display.float_format", lambda v: f"{v:,.6g}")
    mc, st, pr = D["mc"], D["mc"]["stats"], D["mc"]["probabilites"]
    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"]))
    print("\n=== 40.2 Lois déclarées : moments théoriques et tirés\n", D["lois"].to_string(index=False))
    print("\nPERT rotation (0,08 ; 0,15 ; 0,25) :", {k: round(v, 4) for k, v in D["pert_rotation"].items() if k != "loi"})
    print("Corrélation r–R_T :", {k: round(v, 4) for k, v in D["corr_rRT"].items()})
    print("\n=== 40.3 Monte-Carlo (N = 10 000, graine 20261031) — identique au ch. 47")
    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'])}")
    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)}")
    print(f"RCD1 q05 {mc['rcd1_q05']:.3f} | RCD6 choc médiane {mc['rcd6_choc200_mediane']:.3f} q05 {mc['rcd6_choc200_q05']:.3f}")
    print("\nErreur Monte-Carlo théorique (N = 10 000) :\n", D["erreurs"].to_string(index=False))
    rep = D["repetitions"]
    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}")
    print(f"   moyennes des 300 : moyenne {rep['moyenne'].mean():,.0f} | q05 {rep['q05'].mean():,.0f} | q95 {rep['q95'].mean():,.0f}")
    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}")
    print("   graine du traité (20261031) :", {k: (round(v, 4) if isinstance(v, float) else v) for k, v in D["graine_traite"].items()})
    print("\nConvergence (une suite, graine 20261031) :\n", D["convergence"].to_string(index=False))
    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))
    print("\n=== 40.4 Hypercube latin vs simple (N = 1 000, répétitions) :\n", D["lhs"].to_string(index=False))
    print("\n=== 40.5 Mesures de risque (référence : valeur DCF déterministe)\n", D["risque"].to_string(index=False))
    print("RCD1 analytique :", {k: round(v, 5) for k, v in D["rcd_analytique"].items()})
    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"])))
    print("\nStress tests déterministes :\n", D["stress"].to_string(index=False))
    print("Stress inversé :", {k: round(v, 5) for k, v in D["inverse"].items()})
    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))
    print("\nCorrélations de rang / Pearson :\n", D["rangs"].to_string(index=False))
    print("Régression standardisée : R² = %.4f ; bêtas =" % D["reg"]["R2"], {k: round(v, 3) for k, v in D["reg"]["beta"].items()})
    print("\nTornade (ch. 47) :\n", D["tornade"].to_string(index=False))
    print("\n=== 40.8 Portefeuille (cas 4, 48 immeubles) — structure régionale (0,30 ; 0,60)")
    print({k: (round(v, 4) if isinstance(v, float) else v) for k, v in D["portefeuille"]["stats"].items()})
    print(D["portefeuille"]["par_rmr"].to_string(index=False))
    print("\nTrois structures de corrélation :\n", D["portefeuille_structures"].to_string(index=False))
    print("\nTrois structures, choc de revenu sigma_g = 10 % (g = 0) :\n", D["portefeuille_recession"].to_string(index=False))
    print("\n=== Exercices")
    for k, v in D["exercices"].items(): print(k, v)

if __name__ == "__main__":
    ap = argparse.ArgumentParser(); ap.add_argument("--rapide", action="store_true"); a = ap.parse_args()
    imprimer(tout(rapide=a.rapide))