#!/usr/bin/env python3
"""Chapitre 43 — Risques environnementaux et sinistres : DONNÉES SYNTHÉTIQUES et calculs de l'exemple intégral.

Quatre blocs, tous reproductibles (graine 20261031) :

1. Cas 4 (portefeuille, `traite.cas4`) : lecture du coefficient `inondable` et `sol` de `effets_regionaux`
   (prime de risque dans le log du taux global ; loi déclarée du générateur : +5 % du taux en zone inondable,
   +3 % si contamination), part de la valeur du portefeuille située en zone inondable.
2. Cas 4, doubles différences : on injecte sur les séries trimestrielles de RNE des trois immeubles
   `zone_inondable == 1` un choc d'inondation déclaré au printemps 2019 (chute transitoire puis perte
   permanente de 5 % du RNE), puis on l'estime par DiD 2 × 2 (naïf et par strate de RMR), par étude
   d'événement à effets fixes immeuble et RMR × trimestre (`traite.panel`) et par inférence de
   randomisation exacte (toutes les affectations de 3 immeubles traités parmi 48).
3. Municipalité riveraine fictive « Rivière-Fictive » : 625 lots dont 175 en zone inondable (classes très
   élevée / élevée / modérée / faible) le long d'une rivière et 75 lots riverains d'un lac non inondable, ventes
   2016-2026 avec décote déclarée par classe, stigmate après une inondation le 1er mai 2019 et aménité « vue sur
   l'eau » (rivière et lac) qui brouille la mesure si on l'omet. Quatre méthodes de mesure :
   ventes appariées, hédonique à indicatrice (avec et sans contrôle de l'aménité), hédonique par classe avec
   interaction classe × après, doubles différences (cellules 2 × 2 et régression).
4. Modèles de valeur : terrain contaminé (coût de réhabilitation + délai + stigmate), immeuble sinistré
   (remise en état + perte d'usage + stigmate résiduel ; rôle de l'assurance), désuétude vs coût de correction.

Écrit data/synthetic/ch43_riverain_ventes.csv et ch43_riverain_lots.csv s'ils n'existent pas.
Usage : python3 analysis/python/cases/ch43_risques.py
"""
from __future__ import annotations
import sys, pathlib, math, itertools
from datetime import date, timedelta
import numpy as np
import pandas as pd

ROOT = pathlib.Path(__file__).resolve().parents[3]
sys.path.insert(0, str(ROOT / "analysis/python"))
from traite import cas4, comparables as cp, hedonique as h, panel  # noqa: E402

SEED = 20261031
D = ROOT / "data/synthetic"

# =============================================================================== 1. cas 4 : prime de risque dans le taux

def prime_taux_cas4() -> dict:
    """Coefficients `inondable` et `sol` de cas4.effets_regionaux : effet en % du taux, IC à 95 % (HC1),
    comparaison à la loi déclarée (+5 % / +3 %) et part de la valeur en zone inondable."""
    d = cas4.charger(); pf = d["portefeuille"]
    eff = cas4.effets_regionaux(pf); som = cas4.sommaire(pf)
    out = {"n": len(pf), "n_inondable": int(pf["zone_inondable"].sum()), "n_sol": int(pf["contamination_sol"].sum()),
           "part_valeur_inondable": som["total"]["part_zone_inondable"], "valeur_totale": som["total"]["valeur"]}
    for v, declare in [("inondable", 0.05), ("sol", 0.03)]:
        b, se = float(eff.loc[v, "coef"]), float(eff.loc[v, "se_hc1"])
        out[v] = {"coef": b, "se": se, "effet_pct": math.exp(b) - 1, "ic_bas_pct": math.exp(b - 1.96 * se) - 1,
                  "ic_haut_pct": math.exp(b + 1.96 * se) - 1, "declare_pct": declare, "declare_log": math.log(1 + declare),
                  "t_vs_declare": (b - math.log(1 + declare)) / se,
                  # correction de Kennedy (variance non négligeable) : exp(b − V/2) − 1
                  "effet_kennedy_pct": math.exp(b - 0.5 * se ** 2) - 1}
    # décote de valeur impliquée par la prime de taux, à RNE donné : V = RNE/R → ΔlnV = −Δln R
    out["decote_valeur_inondable_pct"] = math.exp(-out["inondable"]["coef"]) - 1
    out["decote_valeur_declaree_pct"] = 1 / 1.05 - 1
    out["immeubles_inondables"] = pf.loc[pf["zone_inondable"] == 1, ["id_immeuble", "rmr", "nb_logements", "valeur_2026", "rne_2026", "taux_global"]].reset_index(drop=True)
    out["immeubles_sol"] = pf.loc[pf["contamination_sol"] == 1, ["id_immeuble", "rmr", "nb_logements", "valeur_2026", "rne_2026", "taux_global"]].reset_index(drop=True)
    out["portefeuille"] = pf
    return out

# =============================================================================== 2. cas 4 : DiD sur les séries de RNE

CHOC_2019 = {"2019Q2": -0.25, "2019Q3": -0.12, "2019Q4": -0.06}   # ln RNE, transitoire (perte d'usage, remise en état)
CHOC_PERMANENT = -0.05                                             # ln RNE à partir de 2020Q1 (loyers, inoccupation, assurance)
T_EVENEMENT = "2019Q2"

def panel_choc(seed: int = SEED) -> pd.DataFrame:
    """Panel long (48 immeubles × 42 trimestres) de ln RNE avec le choc déclaré injecté sur les immeubles
    `zone_inondable == 1` ; colonnes : id, rmr, t (0..41), trimestre, traite, post, ln_rne, ln_rne_sans_choc, g."""
    d = cas4.charger(); pf, se = d["portefeuille"], d["series"].copy()
    se["trimestre"] = se["trimestre"].astype(str)
    se = se.merge(pf[["id_immeuble", "rmr", "zone_inondable"]], on="id_immeuble")
    trims = sorted(se["trimestre"].unique()); idx = {q: i for i, q in enumerate(trims)}
    se["t"] = se["trimestre"].map(idx); t_ev = idx[T_EVENEMENT]
    se["traite"] = se["zone_inondable"].astype(int); se["post"] = (se["t"] >= t_ev).astype(int)
    se["ln_rne_sans_choc"] = np.log(se["rne_trimestriel"].astype(float))
    choc = np.zeros(len(se))
    for q, c in CHOC_2019.items():
        choc[(se["trimestre"] == q).values & (se["traite"] == 1).values] = c
    choc[(se["t"] >= idx["2020Q1"]).values & (se["traite"] == 1).values] = CHOC_PERMANENT
    se["choc"] = choc; se["ln_rne"] = se["ln_rne_sans_choc"] + choc
    se["g"] = np.where(se["traite"] == 1, float(t_ev), np.nan)
    se["unite"] = se["id_immeuble"]
    return se.sort_values(["id_immeuble", "t"]).reset_index(drop=True)

def did_cas4(pan: pd.DataFrame) -> dict:
    trims = sorted(pan["trimestre"].unique()); idx = {q: i for i, q in enumerate(trims)}
    avant = list(range(idx["2017Q1"], idx["2019Q1"] + 1)); apres = list(range(idx["2020Q1"], idx["2026Q2"] + 1))
    r = panel.did_2x2(pan, "ln_rne", "unite", "t", "traite", avant, apres)
    # version par strate : Δ_i centré sur la moyenne des témoins de la même RMR
    dl = r["deltas"].merge(pan[["unite", "rmr"]].drop_duplicates(), on="unite")
    m_rmr = dl[dl["traite"] == 0].groupby("rmr")["dy"].mean()
    dl["dy_strate"] = dl["dy"] - dl["rmr"].map(m_rmr)
    est_strate = float(dl.loc[dl["traite"] == 1, "dy_strate"].mean())
    # inférence de randomisation exacte : toutes les affectations de 3 traités parmi 48 (C(48,3) = 17 296)
    dy = dl["dy"].to_numpy(); ids = np.arange(len(dy)); n_tr = int(dl["traite"].sum())
    tot = dy.sum(); n = len(dy)
    combos = np.array(list(itertools.combinations(ids, n_tr)))
    s_tr = dy[combos].sum(axis=1)
    est_all = s_tr / n_tr - (tot - s_tr) / (n - n_tr)
    obs = r["estimateur"]
    p_ri = float(np.mean(np.abs(est_all) >= abs(obs) - 1e-12))
    # idem pour l'estimateur par strate
    dys = dl["dy_strate"].to_numpy()  # témoins centrés par RMR ; pour les placebos on recentre à chaque tirage
    rmr_codes = dl["rmr"].astype("category").cat.codes.to_numpy()
    est_strate_all = np.empty(len(combos))
    for k, c in enumerate(combos):
        mask = np.zeros(n, bool); mask[c] = True
        # moyenne des témoins par RMR
        num = np.bincount(rmr_codes[~mask], weights=dy[~mask], minlength=rmr_codes.max() + 1)
        den = np.bincount(rmr_codes[~mask], minlength=rmr_codes.max() + 1)
        mr = np.where(den > 0, num / np.maximum(den, 1), 0.0)
        est_strate_all[k] = float(np.mean(dy[mask] - mr[rmr_codes[mask]]))
    p_ri_strate = float(np.mean(np.abs(est_strate_all) >= abs(est_strate) - 1e-12))
    # étude d'événement : effets fixes immeuble + RMR × trimestre, fenêtre −8..+12, référence −1
    es = panel.event_study(pan, "ln_rne", "unite", "t", "g", fenetre=(-8, 12), ref=-1, fe_supp="rmr")
    # TWFE statique avec RMR × trimestre
    pan2 = pan.copy(); pan2["traite_post"] = pan2["traite"] * pan2["post"]
    pan2["traite_permanent"] = pan2["traite"] * (pan2["t"] >= idx["2020Q1"]).astype(int)
    pan2["traite_transitoire"] = pan2["traite"] * pan2["trimestre"].isin(list(CHOC_2019)).astype(int)
    tw = panel.twfe(pan2, "ln_rne", ["traite_transitoire", "traite_permanent"], "unite", "t", fe_supp="rmr")
    # lien avec la valeur : Δln V = Δln RNE − Δln R
    pr = prime_taux_cas4()
    return {"cellules": r["cellules"], "estimateur_2x2": obs, "coef_reg": r["coef_reg"], "se_hc1": r["se_hc1"],
            "estimateur_strate": est_strate, "p_ri": p_ri, "p_ri_strate": p_ri_strate, "n_combinaisons": int(len(combos)),
            "rang_obs": int(np.sum(est_all <= obs)), "n_traites": r["n_traites"], "n_temoins": r["n_temoins"],
            "deltas": dl, "event_study": es, "twfe": tw["coef"], "placebo_dist": est_all, "placebo_dist_strate": est_strate_all,
            "verite_permanente": CHOC_PERMANENT, "verite_permanente_pct": math.exp(CHOC_PERMANENT) - 1,
            "verite_transitoire_moy": float(np.mean(list(CHOC_2019.values()))),
            "effet_valeur_estime_log": tw["coef"].loc["traite_permanent", "coef"] - pr["inondable"]["coef"],
            "effet_valeur_estime_pct": math.exp(tw["coef"].loc["traite_permanent", "coef"] - pr["inondable"]["coef"]) - 1,
            "effet_valeur_declare_log": CHOC_PERMANENT - math.log(1.05),
            "effet_valeur_declare_pct": math.exp(CHOC_PERMANENT - math.log(1.05)) - 1,
            "idx": idx, "trims": trims}

# =============================================================================== 3. Rivière-Fictive : ventes riveraines

NCOL, NLIG, PAS = 25, 25, 40.0            # 625 lots de 40 m × 40 m ; la rivière longe le bord y = 0, un lac (non inondable) le bord x = 24
COLS_LAC = {22, 23, 24}                    # lots riverains du lac : aménité « vue sur l'eau » sans exposition à l'inondation
CLASSES = ["très élevée", "élevée", "modérée", "faible", "hors zone"]
LIGNES_CLASSE = {0: "très élevée", 1: "élevée", 2: "élevée", 3: "modérée", 4: "modérée", 5: "faible", 6: "faible"}
DECOTE_AVANT = {"très élevée": -0.10, "élevée": -0.07, "modérée": -0.04, "faible": -0.015, "hors zone": 0.0}  # ln P, hypothèses
STIGMATE = 0.5                             # après l'inondation du 1er mai 2019, la décote est multipliée par (1 + STIGMATE)
AMENITE_RIVE = 0.04                        # ln P, lots à moins de 120 m de la rivière OU riverains du lac (vue, accès) ; indépendant du risque
T_INOND = date(2019, 5, 1)
DATE_DEBUT, DATE_FIN, DATE_REF = date(2016, 1, 1), date(2026, 8, 31), date(2021, 1, 1)
PRIX_BASE, TAUX_MARCHE, SIGMA = 390_000.0, 0.035, 0.07
BETA = {"ln_hab": 0.55, "ln_terr": 0.10, "age": -0.005, "qualite": 0.06, "garage": 0.035}

def generer_riverain(seed: int = SEED) -> tuple[pd.DataFrame, pd.DataFrame]:
    rng = np.random.default_rng(seed)
    lots = []
    for y in range(NLIG):
        for x in range(NCOL):
            cl = LIGNES_CLASSE.get(y, "hors zone")
            dist = (y + 0.5) * PAS + rng.uniform(-8, 8)
            lots.append({"id_lot": f"RF{y:02d}{x:02d}", "x": x, "y": y, "x_m": (x + 0.5) * PAS, "y_m": (y + 0.5) * PAS,
                         "dist_riviere_m": round(dist, 1), "classe": cl, "expose": int(cl != "hors zone"), "rive": int(dist < 120 or x in COLS_LAC)})
    lots = pd.DataFrame(lots)
    # ventes : chaque lot se vend en moyenne 1,6 fois sur la période (Poisson), dates uniformes
    ventes = []
    n_jours = (DATE_FIN - DATE_DEBUT).days
    k = 0
    for _, L in lots.iterrows():
        nv = rng.poisson(1.6)
        hab = float(np.exp(rng.normal(np.log(140), 0.22))); terr = float(np.clip(rng.normal(560, 110), 300, 1200))
        annee = int(np.clip(rng.normal(1985, 18), 1930, 2016)); qual = int(np.clip(round(rng.normal(3, 0.8)), 1, 5)); gar = int(rng.random() < 0.55)
        for _ in range(nv):
            dv = DATE_DEBUT + timedelta(days=int(rng.integers(0, n_jours + 1)))
            age = dv.year - annee; post = int(dv >= T_INOND)
            decote = DECOTE_AVANT[L["classe"]] * (1 + STIGMATE * post)
            lnp = (math.log(PRIX_BASE) + BETA["ln_hab"] * math.log(hab / 140) + BETA["ln_terr"] * math.log(terr / 560) + BETA["age"] * age
                   + BETA["qualite"] * (qual - 3) + BETA["garage"] * gar + AMENITE_RIVE * L["rive"] + decote
                   + math.log(1 + TAUX_MARCHE) * ((dv - DATE_REF).days / 365.25) + rng.normal(0, SIGMA))
            k += 1
            ventes.append({"id_vente": f"V{k:04d}", "id_lot": L["id_lot"], "date_vente": dv.isoformat(), "prix_vente": round(math.exp(lnp) / 500) * 500,
                           "sup_habitable_m2": round(hab, 1), "sup_terrain_m2": round(terr), "annee_construction": annee, "age": age, "qualite": qual,
                           "garage": gar, "classe": L["classe"], "expose": L["expose"], "rive": L["rive"], "dist_riviere_m": L["dist_riviere_m"],
                           "post": post, "decote_declaree": round(decote, 5)})
    ventes = pd.DataFrame(ventes).sort_values("date_vente").reset_index(drop=True)
    return lots, ventes

def charger_ou_generer(ecrire: bool = True) -> tuple[pd.DataFrame, pd.DataFrame]:
    fl, fv = D / "ch43_riverain_lots.csv", D / "ch43_riverain_ventes.csv"
    lots, ventes = generer_riverain()
    if ecrire and not (fl.exists() and fv.exists()):
        lots.to_csv(fl, index=False); ventes.to_csv(fv, index=False)
    return lots, ventes

def _X(v: pd.DataFrame, colonnes_sup: list[pd.Series]) -> np.ndarray:
    t = (pd.to_datetime(v["date_vente"]) - pd.Timestamp(DATE_REF)).dt.days / 365.25
    base = [np.ones(len(v)), np.log(v["sup_habitable_m2"] / 140), np.log(v["sup_terrain_m2"] / 560), v["age"].astype(float),
            (v["qualite"] - 3).astype(float), v["garage"].astype(float), t]
    return np.column_stack(base + [np.asarray(c, float) for c in colonnes_sup])

def apparier(ventes: pd.DataFrame, fenetre_mois: float = 12.0) -> pd.DataFrame:
    """Ventes appariées après l'inondation : chaque vente exposée est appariée à la vente hors zone la plus proche
    (distance d'attributs : superficies, âge, qualité, garage, date ≤ 12 mois) ; témoin ajusté du temps et des écarts
    physiques par les coefficients hédoniques estimés sur les témoins (ch. 17)."""
    post = ventes[ventes["post"] == 1].copy(); post["d"] = pd.to_datetime(post["date_vente"])
    ex = post[post["expose"] == 1]; te = post[post["expose"] == 0].copy()
    # coefficients physiques sur les témoins seulement (hors zone)
    r = h.mco(_X(te, []), np.log(te["prix_vente"].values)); b = r["beta"]
    ech = {"sup_habitable_m2": 25, "sup_terrain_m2": 150, "age": 10, "qualite": 1, "garage": 1}
    rows = []; pris = set()
    for _, s in ex.iterrows():
        cand = te[(te["d"] - s["d"]).abs() <= pd.Timedelta(days=int(fenetre_mois * 30.4375))]
        cand = cand[~cand["id_vente"].isin(pris)]
        if cand.empty: continue
        dist = sum(((cand[k] - s[k]) / e) ** 2 for k, e in ech.items()) + ((cand["d"] - s["d"]).dt.days / 182.6) ** 2
        tem = cand.loc[dist.idxmin()]; pris.add(tem["id_vente"])
        mois = (s["d"] - tem["d"]).days / 30.4375
        aj_temps = cp.ajustement_temps_taux(tem["prix_vente"], TAUX_MARCHE, mois)
        fact = math.exp(b[1] * math.log(s["sup_habitable_m2"] / tem["sup_habitable_m2"]) + b[2] * math.log(s["sup_terrain_m2"] / tem["sup_terrain_m2"])
                        + b[3] * (s["age"] - tem["age"]) + b[4] * (s["qualite"] - tem["qualite"]) + b[5] * (s["garage"] - tem["garage"]))
        p_aj = (tem["prix_vente"] + aj_temps) * fact
        rows.append({"id_expose": s["id_vente"], "id_temoin": tem["id_vente"], "classe": s["classe"], "rive": s["rive"], "prix_expose": s["prix_vente"],
                     "prix_temoin": tem["prix_vente"], "ajust_temps": aj_temps, "facteur_physique": fact, "prix_temoin_ajuste": p_aj,
                     "difference": s["prix_vente"] - p_aj, "difference_pct": s["prix_vente"] / p_aj - 1, "decote_declaree": s["decote_declaree"]})
    return pd.DataFrame(rows)

def resume_paires(p: pd.DataFrame) -> dict:
    r = cp.ventes_appariees(p["prix_expose"].tolist(), p["prix_temoin_ajuste"].tolist())
    pct = p["difference_pct"]; n = len(p)
    out = {"n": n, "moyenne": r["moyenne"], "mediane": r["mediane"], "ecart_type": r["ecart_type"], "erreur_type": r["ecart_type"] / math.sqrt(n),
           "moyenne_pct": float(pct.mean()), "mediane_pct": float(pct.median()), "et_pct": float(pct.std(ddof=1) / math.sqrt(n)),
           "verite_pct": float(np.exp(p["decote_declaree"]).mean() - 1), "verite_log": float(p["decote_declaree"].mean())}
    out["ic95_pct"] = (out["moyenne_pct"] - 1.96 * out["et_pct"], out["moyenne_pct"] + 1.96 * out["et_pct"])
    out["par_classe"] = p.groupby("classe")["difference_pct"].agg(["count", "mean"]).rename(columns={"count": "n", "mean": "moyenne_pct"})
    out["par_rive"] = p.groupby("rive")["difference_pct"].agg(["count", "mean"]).rename(columns={"count": "n", "mean": "moyenne_pct"})
    return out

def hedoniques(ventes: pd.DataFrame) -> dict:
    """(a) indicatrice `exposé` sans contrôle de l'aménité ; (b) avec `rive` ; (c) par classe avec interactions classe × après ;
    (d) DiD par régression exposé × après (avec rive)."""
    v = ventes; y = np.log(v["prix_vente"].values); out = {}
    ra = h.mco(_X(v, [v["expose"]]), y)
    out["a_sans_rive"] = {"coef": ra["beta"][7], "se": ra["se_hc1"][7], "r2": ra["r2"], "sigma": math.sqrt(ra["sigma2"])}
    rb = h.mco(_X(v, [v["expose"], v["rive"]]), y)
    out["b_avec_rive"] = {"coef": rb["beta"][7], "se": rb["se_hc1"][7], "coef_rive": rb["beta"][8], "se_rive": rb["se_hc1"][8], "r2": rb["r2"], "sigma": math.sqrt(rb["sigma2"])}
    cls = ["très élevée", "élevée", "modérée", "faible"]
    cols = [(v["classe"] == c).astype(float) for c in cls] + [(v["classe"] == c).astype(float) * v["post"] for c in cls] + [v["rive"]]
    rc = h.mco(_X(v, cols), y)
    noms = [f"{c}" for c in cls] + [f"{c} × après" for c in cls] + ["rive"]
    out["c_classes"] = pd.DataFrame({"coef": rc["beta"][7:], "se": rc["se_hc1"][7:]}, index=noms)
    out["c_classes"]["effet_pct"] = np.exp(out["c_classes"]["coef"]) - 1
    out["c_classes"]["verite"] = [DECOTE_AVANT[c] for c in cls] + [DECOTE_AVANT[c] * STIGMATE for c in cls] + [AMENITE_RIVE]
    out["c_classes"]["r2"] = rc["r2"]
    rd = h.mco(_X(v, [v["expose"], v["expose"] * v["post"], v["rive"]]), y)   # `post` seul est colinéaire avec la tendance ? non : on l'ajoute
    rd = h.mco(_X(v, [v["expose"], v["post"], v["expose"] * v["post"], v["rive"]]), y)
    out["d_did_reg"] = {"coef_expose": rd["beta"][7], "se_expose": rd["se_hc1"][7], "coef_post": rd["beta"][8], "coef_did": rd["beta"][9], "se_did": rd["se_hc1"][9],
                        "coef_rive": rd["beta"][10], "r2": rd["r2"], "sigma": math.sqrt(rd["sigma2"])}
    # cellules 2 × 2 sur ln P (sans contrôle)
    m = v.assign(lnp=y).groupby(["expose", "post"])["lnp"].mean()
    out["d_did_2x2"] = {"moyennes": m, "estimateur": float((m[(1, 1)] - m[(1, 0)]) - (m[(0, 1)] - m[(0, 0)])),
                        "effectifs": v.groupby(["expose", "post"]).size().to_dict()}
    # vérités : décote moyenne (log) des ventes exposées avant / après, stigmate moyen
    ex = v[v["expose"] == 1]
    out["verite"] = {"decote_avant_log": float(ex.loc[ex["post"] == 0, "decote_declaree"].mean()), "decote_apres_log": float(ex.loc[ex["post"] == 1, "decote_declaree"].mean()),
                     "decote_toutes_log": float(ex["decote_declaree"].mean())}
    out["verite"]["stigmate_log"] = out["verite"]["decote_apres_log"] - out["verite"]["decote_avant_log"]
    out["effectifs"] = {"ventes": len(v), "exposees": int(v["expose"].sum()), "rive": int(v["rive"].sum()), "apres": int(v["post"].sum())}
    return out

def effets_par_methode(ventes: pd.DataFrame, paires: pd.DataFrame) -> pd.DataFrame:
    """Tableau des estimations de la décote (ln P) des ventes exposées après l'inondation, par méthode, avec IC à 95 %."""
    he = hedoniques(ventes); rp = resume_paires(paires)
    lignes = []
    lignes.append({"methode": "Ventes appariées (après 2019, sans contrôle de l'aménité)", "coef": math.log(1 + rp["moyenne_pct"]), "se": rp["et_pct"] / (1 + rp["moyenne_pct"]), "cible": "décote après"})
    lignes.append({"methode": "Hédonique, indicatrice exposé (sans « rive »)", "coef": he["a_sans_rive"]["coef"], "se": he["a_sans_rive"]["se"], "cible": "décote moyenne"})
    lignes.append({"methode": "Hédonique, indicatrice exposé + « rive »", "coef": he["b_avec_rive"]["coef"], "se": he["b_avec_rive"]["se"], "cible": "décote moyenne"})
    d = he["d_did_reg"]
    lignes.append({"methode": "DiD régression : exposé (avant)", "coef": d["coef_expose"], "se": d["se_expose"], "cible": "décote avant"})
    lignes.append({"methode": "DiD régression : exposé × après (stigmate)", "coef": d["coef_did"], "se": d["se_did"], "cible": "stigmate"})
    lignes.append({"methode": "DiD régression : décote après = somme", "coef": d["coef_expose"] + d["coef_did"], "se": float("nan"), "cible": "décote après"})
    lignes.append({"methode": "DiD 2 × 2 sur les moyennes de ln P (stigmate)", "coef": he["d_did_2x2"]["estimateur"], "se": float("nan"), "cible": "stigmate"})
    t = pd.DataFrame(lignes)
    ver = he["verite"]
    t["verite"] = t["cible"].map({"décote après": ver["decote_apres_log"], "décote moyenne": ver["decote_toutes_log"], "décote avant": ver["decote_avant_log"], "stigmate": ver["stigmate_log"]})
    t["effet_pct"] = np.exp(t["coef"]) - 1; t["ic_bas"] = t["coef"] - 1.96 * t["se"]; t["ic_haut"] = t["coef"] + 1.96 * t["se"]
    return t

# =============================================================================== 4. modèles de valeur

def valeur_terrain_contamine(v_propre: float, cout_rehab: float, delai_mois: float, r: float, stigmate: float, contingence: float = 0.0, perte_revenu_mensuelle: float = 0.0) -> dict:
    """V_c = V_p (1 − s)(1 + r)^(−T) − C(1 + k)(1 + r)^(−T/2) − Σ pertes de revenu actualisées pendant les travaux.
    T en années ; C payé à mi-parcours ; pertes mensuelles actualisées au taux périodique équivalent."""
    T = delai_mois / 12.0; rm = (1 + r) ** (1 / 12) - 1
    va_propre = v_propre * (1 - stigmate) / (1 + r) ** T
    va_cout = cout_rehab * (1 + contingence) / (1 + r) ** (T / 2)
    va_perte = sum(perte_revenu_mensuelle / (1 + rm) ** m for m in range(1, int(round(delai_mois)) + 1))
    v = va_propre - va_cout - va_perte
    return {"valeur": v, "va_propre_stigmatisee": va_propre, "va_cout": va_cout, "va_perte_revenu": va_perte, "cout_delai": v_propre * (1 - stigmate) - va_propre,
            "cout_stigmate": v_propre * stigmate, "ecart_total": v_propre - v, "ecart_pct": v / v_propre - 1}

def valeur_apres_sinistre(v_avant: float, cout_remise: float, loyer_mensuel_perdu: float, mois_perte: float, r: float, stigmate: float, indemnite: float = 0.0) -> dict:
    """Valeur de l'immeuble sinistré non réparé : V_avant − coût de remise en état − VA de la perte d'usage − stigmate résiduel × V_avant ;
    position du propriétaire assuré : valeur + indemnité (l'indemnité couvre des coûts, jamais le stigmate)."""
    rm = (1 + r) ** (1 / 12) - 1
    va_perte = sum(loyer_mensuel_perdu / (1 + rm) ** m for m in range(1, int(round(mois_perte)) + 1))
    stig = stigmate * v_avant
    v = v_avant - cout_remise - va_perte - stig
    return {"valeur": v, "va_perte_usage": va_perte, "stigmate_dollars": stig, "perte_totale": v_avant - v, "perte_pct": v / v_avant - 1,
            "position_assure": v + indemnite, "perte_nette_assure": v_avant - (v + indemnite), "part_non_assuree": (v_avant - (v + indemnite)) / (v_avant - v)}

def desuetude_vs_correction(cout_correction: float, perte_usage: float, perte_rne_annuelle: float, R: float, curable: bool) -> dict:
    """Pathologie curable : D = coût de correction + perte d'usage ; incurable : D = perte de RNE / R (capitalisation) — ch. 24."""
    d_cur = cout_correction + perte_usage; d_inc = perte_rne_annuelle / R
    return {"curable": d_cur, "incurable": d_inc, "retenu": d_cur if curable else d_inc, "regle": "curable si coût de correction ≤ gain de valeur"}

def exemple_contamination() -> dict:
    """Hypothèses déclarées (fictives) : terrain d'un ancien garage avec immeuble de 8 logements à construire ou existant ;
    valeur propre 2 100 000 $ ; phase I 3 500 $, phase II 18 000 $, réhabilitation 185 000 $ (+ 25 % contingence) ; 15 mois ; 6,5 % ; stigmate 5 %."""
    hyp = {"v_propre": 2_100_000.0, "phase1": 3_500.0, "phase2": 18_000.0, "rehab": 185_000.0, "contingence": 0.25, "delai_mois": 15, "r": 0.065, "stigmate": 0.05, "perte_mensuelle": 0.0}
    c_total = hyp["phase1"] + hyp["phase2"] + hyp["rehab"]
    base = valeur_terrain_contamine(hyp["v_propre"], c_total, hyp["delai_mois"], hyp["r"], hyp["stigmate"], hyp["contingence"], hyp["perte_mensuelle"])
    sens = {}
    for s in [0.0, 0.05, 0.10, 0.15]:
        for T in [6, 15, 24]:
            sens[(s, T)] = valeur_terrain_contamine(hyp["v_propre"], c_total, T, hyp["r"], s, hyp["contingence"])["valeur"]
    sens_cout = {k: valeur_terrain_contamine(hyp["v_propre"], c_total * k, hyp["delai_mois"], hyp["r"], hyp["stigmate"], hyp["contingence"])["valeur"] for k in [0.7, 1.0, 1.5, 2.0]}
    return {"hyp": hyp, "cout_total": c_total, "base": base, "sensibilite": sens, "sensibilite_cout": sens_cout}

def exemple_sinistre() -> dict:
    """Immeuble P41 du cas 4 (Montréal, zone inondable) : inondation fictive ; hypothèses déclarées."""
    pf = cas4.charger()["portefeuille"]; p = pf[pf["id_immeuble"] == "P41"].iloc[0]
    hyp = {"v_avant": float(p["valeur_2026"]), "rne": float(p["rne_2026"]), "nb_logements": int(p["nb_logements"]), "loyer": float(p["loyer_mensuel_moyen"]),
           "cout_remise": 240_000.0, "logements_touches": 4, "mois_perte": 5, "r": 0.065, "stigmate": 0.05, "franchise": 10_000.0, "part_couverte": 0.85}
    loyer_perdu = hyp["loyer"] * hyp["logements_touches"]
    indem = max(0.0, hyp["part_couverte"] * hyp["cout_remise"] - hyp["franchise"])
    base = valeur_apres_sinistre(hyp["v_avant"], hyp["cout_remise"], loyer_perdu, hyp["mois_perte"], hyp["r"], hyp["stigmate"], indem)
    sens = {s: valeur_apres_sinistre(hyp["v_avant"], hyp["cout_remise"], loyer_perdu, hyp["mois_perte"], hyp["r"], s, indem)["valeur"] for s in [0.0, 0.025, 0.05, 0.10, 0.15]}
    sans_assurance = valeur_apres_sinistre(hyp["v_avant"], hyp["cout_remise"], loyer_perdu, hyp["mois_perte"], hyp["r"], hyp["stigmate"], 0.0)
    return {"hyp": hyp, "loyer_perdu_mensuel": loyer_perdu, "indemnite": indem, "base": base, "sensibilite_stigmate": sens, "sans_assurance": sans_assurance, "immeuble": p}

def exemple_pyrite() -> dict:
    """Pyrite (renvoi ch. 11 et 24) : excavation 25 000 $, perte d'usage 2 000 $ ; variante incurable : perte de RNE 1 800 $/an capitalisée à 5,5 %."""
    return {"curable": desuetude_vs_correction(25_000, 2_000, 0, 0.055, True), "incurable": desuetude_vs_correction(0, 0, 1_800, 0.055, False)}

# =============================================================================== main

def main(ecrire: bool = True, verbeux: bool = True) -> dict:
    pr = prime_taux_cas4(); pan = panel_choc(); dd = did_cas4(pan)
    lots, ventes = charger_ou_generer(ecrire); paires = apparier(ventes); rp = resume_paires(paires); he = hedoniques(ventes); epm = effets_par_methode(ventes, paires)
    ct = exemple_contamination(); sn = exemple_sinistre(); py = exemple_pyrite()
    res = {"prime": pr, "panel": pan, "did": dd, "lots": lots, "ventes": ventes, "paires": paires, "resume_paires": rp, "hedoniques": he, "effets_methodes": epm,
           "contamination": ct, "sinistre": sn, "pyrite": py}
    if not verbeux: return res
    P = print
    P("=== 1. Cas 4 — prime de risque dans le taux global (effets_regionaux) ===")
    P(f"n = {pr['n']} immeubles ; inondables = {pr['n_inondable']} ; contaminés = {pr['n_sol']} ; part de la valeur en zone inondable = {pr['part_valeur_inondable']*100:.2f} %")
    for v in ["inondable", "sol"]:
        e = pr[v]
        P(f"{v:10s}: coef = {e['coef']:.4f} (se HC1 {e['se']:.4f}) → +{e['effet_pct']*100:.2f} % du taux [IC95 {e['ic_bas_pct']*100:.2f} ; {e['ic_haut_pct']*100:.2f}] ; Kennedy {e['effet_kennedy_pct']*100:.2f} % ; déclaré +{e['declare_pct']*100:.0f} % (t = {e['t_vs_declare']:.2f})")
    P(f"décote de valeur impliquée (à RNE donné) : {pr['decote_valeur_inondable_pct']*100:.2f} % (déclarée {pr['decote_valeur_declaree_pct']*100:.2f} %)")
    P(pr["immeubles_inondables"].to_string(index=False))
    P("\n=== 2. Cas 4 — DiD sur ln RNE, choc 2019 injecté ===")
    P(f"vérité : transitoire {CHOC_2019} ; permanent {CHOC_PERMANENT} ({dd['verite_permanente_pct']*100:.2f} %)")
    P("cellules (ln RNE moyen) :"); P(dd["cellules"].to_string())
    P(f"DiD 2×2 naïf = {dd['estimateur_2x2']:.4f} ; régression Δ (se HC1) = {dd['coef_reg']:.4f} ({dd['se_hc1']:.4f}) ; par strate RMR = {dd['estimateur_strate']:.4f}")
    P(f"inférence de randomisation exacte : {dd['n_combinaisons']} affectations ; p (naïf) = {dd['p_ri']:.4f} ; p (strate) = {dd['p_ri_strate']:.4f}")
    P("TWFE immeuble + RMR×trimestre :"); P(dd["twfe"].to_string())
    P("étude d'événement (k, coef, se grappes) :"); P(dd["event_study"][["coef", "se_grappes", "ic_bas", "ic_haut"]].round(4).to_string())
    P(f"effet sur la valeur (Δln RNE − Δln R) : estimé {dd['effet_valeur_estime_pct']*100:.2f} % ; déclaré {dd['effet_valeur_declare_pct']*100:.2f} %")
    P("\n=== 3. Rivière-Fictive — ventes riveraines ===")
    P(f"lots {len(lots)} ; exposés {int(lots['expose'].sum())} ; par classe {lots['classe'].value_counts().to_dict()}")
    P(f"ventes {he['effectifs']} ; vérités (ln) : avant {he['verite']['decote_avant_log']:.4f}, après {he['verite']['decote_apres_log']:.4f}, stigmate {he['verite']['stigmate_log']:.4f}")
    P(f"paires : n = {rp['n']} ; moyenne {rp['moyenne']:.0f} $ ({rp['moyenne_pct']*100:.2f} %) ; médiane {rp['mediane_pct']*100:.2f} % ; IC95 [{rp['ic95_pct'][0]*100:.2f} ; {rp['ic95_pct'][1]*100:.2f}] ; vérité {rp['verite_pct']*100:.2f} %")
    P(rp["par_classe"].to_string()); P(rp["par_rive"].to_string())
    P(f"hédonique (a) exposé sans rive : {he['a_sans_rive']['coef']:.4f} ({he['a_sans_rive']['se']:.4f}) ; (b) avec rive : {he['b_avec_rive']['coef']:.4f} ({he['b_avec_rive']['se']:.4f}), rive {he['b_avec_rive']['coef_rive']:.4f} ({he['b_avec_rive']['se_rive']:.4f}) ; σ = {he['b_avec_rive']['sigma']:.4f}, R² = {he['b_avec_rive']['r2']:.3f}")
    P("(c) par classe :"); P(he["c_classes"].round(4).to_string())
    d = he["d_did_reg"]; P(f"(d) DiD régression : exposé {d['coef_expose']:.4f} ({d['se_expose']:.4f}) ; exposé×après {d['coef_did']:.4f} ({d['se_did']:.4f}) ; post {d['coef_post']:.4f} ; rive {d['coef_rive']:.4f}")
    P(f"    DiD 2×2 : {he['d_did_2x2']['estimateur']:.4f} ; effectifs {he['d_did_2x2']['effectifs']}")
    P("effets par méthode :"); P(epm.round(4).to_string(index=False))
    P("\n=== 4. Modèles de valeur ===")
    b = ct["base"]; P(f"contamination : coût total {ct['cout_total']:.0f} $ ; V_c = {b['valeur']:.0f} $ ({b['ecart_pct']*100:.2f} %) ; VA propre stigmatisée {b['va_propre_stigmatisee']:.0f} ; VA coût {b['va_cout']:.0f} ; coût du délai {b['cout_delai']:.0f} ; stigmate {b['cout_stigmate']:.0f}")
    P("sensibilité (stigmate, mois) → valeur :"); P({k: round(v) for k, v in ct["sensibilite"].items()}); P("coût ×k :"); P({k: round(v) for k, v in ct["sensibilite_cout"].items()})
    s = sn["base"]; P(f"sinistre P41 : V_avant {sn['hyp']['v_avant']:.0f} ; loyer perdu {sn['loyer_perdu_mensuel']:.0f} $/mois × {sn['hyp']['mois_perte']} ; VA perte {s['va_perte_usage']:.0f} ; stigmate {s['stigmate_dollars']:.0f} ; V_après {s['valeur']:.0f} ({s['perte_pct']*100:.2f} %) ; indemnité {sn['indemnite']:.0f} ; position assuré {s['position_assure']:.0f} ; perte nette {s['perte_nette_assure']:.0f} ; part non assurée {s['part_non_assuree']*100:.1f} %")
    P("sensibilité stigmate :"); P({k: round(v) for k, v in sn["sensibilite_stigmate"].items()})
    P(f"pyrite : curable {py['curable']['retenu']:.0f} $ ; incurable {py['incurable']['retenu']:.0f} $")
    return res

if __name__ == "__main__":
    main()