#!/usr/bin/env python3
"""Chapitre 36 — Données de panel et identification causale : choc simulé sur les séries de RNE du cas 4
(DONNÉES SYNTHÉTIQUES, graine 20261031).

Le panel de départ est celui de `traite.cas4.charger()` : 48 immeubles × 42 trimestres (2016T1-2026T2) de RNE trimestriel.
Ces séries sont prises comme résultats potentiels SANS traitement, Y(0). Le script y superpose un traitement déclaré :
  - cohorte A : 8 immeubles de la RMR de Montréal « à distance de marche » d'une ligne de transport collectif fictive mise
    en service au 3e trimestre de 2022 (période 26) ;
  - cohorte B : 5 immeubles de la RMR de Gatineau desservis par un second tronçon fictif au 3e trimestre de 2024 (période 34) ;
  - effet déclaré sur ln RNE, identique pour les deux cohortes, fonction du temps écoulé k depuis la mise en service :
    δ(k) = 0,01 (k + 1) pour k = 0, 1, 2, 3 (montée en régime au rythme du renouvellement des baux), puis
    δ(k) = 0,04 + 0,0025 (k − 3) pour k ≥ 4 (régime de +4 % atteint, puis valorisation lente du voisinage).
Les 35 autres immeubles ne sont jamais traités. Le lecteur connaît la vérité ; chaque estimateur lui est confronté.

Estimateurs : DiD 2 × 2 (cellules et différences longues), TWFE statique (erreurs-types classiques / HC1 / grappes),
étude d'événement, Callaway–Sant'Anna (contrôles jamais traités, même RMR), décomposition de Goodman-Bacon, contrôle
synthétique (cohorte A agrégée vs 35 donneurs), discontinuité spatiale (simulation transversale séparée), placebo temporel,
sensibilité (tendances parallèles violées, choix du groupe témoin, traitement de l'autocorrélation).
Contrôle croisé (--controle) : linearmodels.PanelOLS et statsmodels (grappes) reproduisent le TWFE et ses erreurs-types.
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch36_panel.py [--controle]"""
from __future__ import annotations
import sys, pathlib
import numpy as np, pandas as pd
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
from traite import cas4
from traite import hedonique as h
from traite import panel as pn

ROOT = pathlib.Path(__file__).resolve().parents[3]
GRAINE = 20261031
T_A, T_B = 26, 34                      # 2022T3 et 2024T3 (2016T1 = 0)
N_A, N_B = 8, 5                        # immeubles traités à Montréal et à Gatineau
RMR_A, RMR_B = "Montréal", "Gatineau"
FENETRE_AVANT = [22, 23, 24, 25]       # 2021T3-2022T2
FENETRE_APRES = [30, 31, 32, 33]       # 2023T3-2024T2 (deuxième année : régime atteint)
PLACEBO_T = 14                         # 2019T3
B_BOOT = 499

_CACHE: dict = {}

def effet_vrai(k) -> np.ndarray:
    """Effet déclaré δ(k) sur ln RNE selon le temps relatif k (0 avant la mise en service)."""
    k = np.asarray(k, float)
    return np.where(k < 0, 0.0, np.where(k <= 3, 0.01 * (k + 1), 0.04 + 0.0025 * (k - 3)))

# ----------------------------------------------------------------------------- panel simulé

def panel_base() -> pd.DataFrame:
    d = cas4.charger(); se = d["series"].copy(); pf = d["portefeuille"]
    trims = sorted(se["trimestre"].unique()); idx = {p: i for i, p in enumerate(trims)}
    se["t"] = se["trimestre"].map(idx).astype(int)
    se = se.merge(pf[["id_immeuble", "rmr", "nb_logements", "zone_inondable"]], on="id_immeuble")
    se["y0"] = np.log(se["rne_trimestriel"].astype(float))
    return se.sort_values(["id_immeuble", "t"]).reset_index(drop=True)

def choisir_traites(pf_ids: pd.DataFrame, graine: int = GRAINE) -> dict:
    rng = np.random.default_rng(graine)
    a = sorted(rng.choice(sorted(pf_ids.loc[pf_ids["rmr"] == RMR_A, "id_immeuble"]), N_A, replace=False))
    b = sorted(rng.choice(sorted(pf_ids.loc[pf_ids["rmr"] == RMR_B, "id_immeuble"]), N_B, replace=False))
    return {"A": list(a), "B": list(b)}

def simuler(tendance_diff_annuelle: float = 0.0) -> pd.DataFrame:
    """Panel long avec traitement échelonné. `tendance_diff_annuelle` > 0 ajoute aux immeubles traités une croissance
    supplémentaire (par an, en log) sur toute la période : violation déclarée des tendances parallèles."""
    se = panel_base(); ids = se.groupby("id_immeuble")[["rmr"]].first().reset_index()
    tr = choisir_traites(ids)
    se["g"] = np.nan
    se.loc[se["id_immeuble"].isin(tr["A"]), "g"] = T_A
    se.loc[se["id_immeuble"].isin(tr["B"]), "g"] = T_B
    se["cohorte"] = np.where(se["id_immeuble"].isin(tr["A"]), "A", np.where(se["id_immeuble"].isin(tr["B"]), "B", "jamais"))
    se["k"] = (se["t"] - se["g"]).where(se["g"].notna())
    se["D"] = ((se["k"] >= 0) & se["g"].notna()).astype(float)
    se["delta_vrai"] = np.where(se["g"].notna(), effet_vrai(se["k"].fillna(-1)), 0.0)
    if tendance_diff_annuelle:
        se.loc[se["g"].notna(), "y0"] = se.loc[se["g"].notna(), "y0"] + tendance_diff_annuelle * se.loc[se["g"].notna(), "t"] / 4.0
    se["y"] = se["y0"] + se["delta_vrai"]
    se["traite_A"] = se["id_immeuble"].isin(tr["A"]).astype(int)
    return se

def verite(se: pd.DataFrame) -> dict:
    post = se[se["D"] == 1]
    return {"att_global": float(post["delta_vrai"].mean()),
            "att_A": float(post.loc[post["cohorte"] == "A", "delta_vrai"].mean()),
            "att_B": float(post.loc[post["cohorte"] == "B", "delta_vrai"].mean()),
            "att_fenetre_2x2": float(effet_vrai(np.array(FENETRE_APRES) - T_A).mean()),
            "profil": {int(k): float(effet_vrai(k)) for k in range(0, 16)}}

# ----------------------------------------------------------------------------- estimations

def did_deux_periodes(se: pd.DataFrame, temoins: str = "tous") -> dict:
    """2 × 2 : cohorte A vs jamais traités (tous, ou de la même RMR)."""
    d = se[se["cohorte"] != "B"].copy()
    if temoins == "meme_rmr":
        d = d[d["rmr"] == RMR_A]
    return pn.did_2x2(d, "y", "id_immeuble", "t", "traite_A", FENETRE_AVANT, FENETRE_APRES)

def placebo_temporel(se: pd.DataFrame) -> dict:
    d = se[(se["cohorte"] != "B") & (se["rmr"] == RMR_A)].copy()
    avant = [PLACEBO_T - 4, PLACEBO_T - 3, PLACEBO_T - 2, PLACEBO_T - 1]; apres = [PLACEBO_T, PLACEBO_T + 1, PLACEBO_T + 2, PLACEBO_T + 3]
    return pn.did_2x2(d, "y", "id_immeuble", "t", "traite_A", avant, apres)

def twfe_statique(se: pd.DataFrame) -> dict:
    commun = pn.twfe(se, "y", ["D"], "id_immeuble", "t")
    regional = pn.twfe(se, "y", ["D"], "id_immeuble", "t", fe_supp="rmr")
    # équivalence within : D doublement centré (cylindré) → même coefficient que les indicatrices
    yw, Xw = pn.within_deux_voies(se["y"].values, se[["D"]].values, se["id_immeuble"].values, se["t"].values)
    bw = float(np.linalg.lstsq(Xw, yw, rcond=None)[0][0])
    return {"commun": commun, "regional": regional, "within_deux_voies": bw}

def etude_evenement(se: pd.DataFrame, fe_supp: str | None = "rmr") -> pd.DataFrame:
    es = pn.event_study(se, "y", "id_immeuble", "t", "g", fenetre=(-8, 8), ref=-1, fe_supp=fe_supp)
    # vérité par k : moyenne de δ sur les observations traitées regroupées dans chaque indicatrice
    tr = se[se["g"].notna()].copy(); tr["kb"] = tr["k"].clip(-8, 8)
    v = tr.groupby("kb")["delta_vrai"].mean()
    es["vrai"] = [float(v.get(k, 0.0)) if k >= 0 else 0.0 for k in es.index]
    return es

def cs(se: pd.DataFrame, strate: str | None = "rmr", B: int = B_BOOT) -> dict:
    return pn.callaway_santanna(se, "y", "id_immeuble", "t", "g", strate=strate, B=B)

def synthetique(se: pd.DataFrame) -> dict:
    """Unité traitée = moyenne des ln RNE de la cohorte A, centrée sur sa moyenne pré-traitement ; donneurs = 35 jamais traités,
    chacun centré de même (contrôle synthétique « en écarts à la moyenne pré-traitement »)."""
    W = se.pivot(index="id_immeuble", columns="t", values="y")
    coh = se.groupby("id_immeuble")["cohorte"].first()
    A = W.loc[coh[coh == "A"].index].mean(axis=0).values
    donneurs = W.loc[coh[coh == "jamais"].index]
    T0 = T_A
    A_c = A - A[:T0].mean()
    Y0 = (donneurs.values - donneurs.values[:, :T0].mean(axis=1, keepdims=True)).T   # (T × J)
    r = pn.controle_synthetique(Y0, A_c, T0)
    r["poids"] = pd.Series(r["poids"], index=donneurs.index)
    r["traite"] = A_c; r["T0"] = T0
    r["vrai_post"] = float(effet_vrai(np.arange(0, W.shape[1] - T0)).mean())
    r["placebos"] = pn.placebos_synthetiques(Y0, T0); r["placebos"]["donneur"] = donneurs.index[r["placebos"]["donneur"].values]
    r["ratio_traite"] = r["rmspe_post"] / r["rmspe_pre"]
    r["rang_ratio"] = int((r["placebos"]["ratio"] >= r["ratio_traite"]).sum() + 1)
    return r

def discontinuite(graine: int = GRAINE, n: int = 400, saut: float = 0.05) -> dict:
    """Simulation transversale (fictive) : n immeubles répartis de part et d'autre d'une frontière de zone (x = distance signée en km,
    positive du côté de la zone dense), ln V = 12,9 + 0,03 x + saut·1[x ≥ 0] + ε, ε ~ N(0 ; 0,08²)."""
    rng = np.random.default_rng(graine + 7)
    x = rng.uniform(-2, 2, n); y = 12.9 + 0.03 * x + saut * (x >= 0) + rng.normal(0, 0.08, n)
    res = {hb: pn.rdd_local_lineaire(x, y, hb) for hb in (0.5, 1.0, 2.0)}
    naif = float(y[x >= 0].mean() - y[x < 0].mean())
    return {"x": x, "y": y, "saut_vrai": saut, "estimations": res, "naif": naif}

def sensibilite(se: pd.DataFrame, tendance: float = 0.01) -> dict:
    """Tendances parallèles violées : +`tendance` par an (log) pour les immeubles traités depuis 2016."""
    sv = simuler(tendance_diff_annuelle=tendance)
    did_v = pn.did_2x2(sv[(sv["cohorte"] != "B") & (sv["rmr"] == RMR_A)], "y", "id_immeuble", "t", "traite_A", FENETRE_AVANT, FENETRE_APRES)
    tw_v = pn.twfe(sv, "y", ["D"], "id_immeuble", "t", fe_supp="rmr")
    es_v = pn.event_study(sv, "y", "id_immeuble", "t", "g", fenetre=(-8, 8), ref=-1, fe_supp="rmr")
    # pente des coefficients pré-traitement (k = −8 … −2) : diagnostic
    pre = es_v.loc[[k for k in es_v.index if -8 < k < -1]]
    pente_pre = float(np.polyfit(pre.index.values.astype(float), pre["coef"].values, 1)[0])
    # correction par tendances linéaires propres aux unités traitées (groupe × tendance)
    sv2 = sv.copy(); sv2["tend_tr"] = sv2["g"].notna().astype(float) * sv2["t"] / 4.0
    tw_tend = pn.twfe(sv2, "y", ["D", "tend_tr"], "id_immeuble", "t", fe_supp="rmr")
    # correction par extrapolation de la tendance pré-traitement (2 × 2 détendancé) : Δy_i − pente_i × écart temporel
    d = sv[(sv["cohorte"] != "B") & (sv["rmr"] == RMR_A) & (sv["t"] < T_A)]
    pentes = d.groupby("id_immeuble").apply(lambda g: np.polyfit(g["t"] / 4.0, g["y"], 1)[0], include_groups=False)
    deltas = did_v["deltas"].set_index("id_immeuble")
    ecart_ans = (np.mean(FENETRE_APRES) - np.mean(FENETRE_AVANT)) / 4.0
    deltas["dy_det"] = deltas["dy"] - pentes.reindex(deltas.index).values * ecart_ans
    did_det = float(deltas.loc[deltas["traite_A"] == 1, "dy_det"].mean() - deltas.loc[deltas["traite_A"] == 0, "dy_det"].mean())
    return {"tendance": tendance, "did_2x2": did_v, "twfe_regional": tw_v, "event_study": es_v, "pente_pre": pente_pre,
            "twfe_tendance": tw_tend, "did_detendance": did_det, "ecart_ans": ecart_ans}

def groupes_temoins(se: pd.DataFrame) -> pd.DataFrame:
    """Sensibilité au choix des témoins pour la cohorte A, 2 × 2 sur les fenêtres de référence."""
    rows = []
    for nom, d in [("jamais traités, toutes RMR", se[se["cohorte"] != "B"]),
                   ("jamais traités, Montréal", se[(se["cohorte"] != "B") & (se["rmr"] == RMR_A)]),
                   ("jamais traités, hors Montréal", se[(se["cohorte"] != "B") & ((se["rmr"] != RMR_A) | (se["traite_A"] == 1))]),
                   ("pas encore traités (cohorte B)", se[(se["cohorte"] == "B") | (se["traite_A"] == 1)])]:
        r = pn.did_2x2(d, "y", "id_immeuble", "t", "traite_A", FENETRE_AVANT, FENETRE_APRES)
        rows.append({"temoins": nom, "n_temoins": r["n_temoins"], "estimateur": r["estimateur"], "se_hc1": r["se_hc1"]})
    return pd.DataFrame(rows)

# ----------------------------------------------------------------------------- laboratoire : biais du TWFE échelonné

def laboratoire(graine: int = GRAINE, n_u: int = 10, n_k: int = 20, n_l: int = 20, T: int = 42, g_k: int = 10, g_l: int = 26,
                pente_effet: float = 0.005, sigma: float = 0.02) -> dict:
    """Panel artificiel à tendance commune (α_i + λ_t + ε, sans tendance régionale) : deux cohortes (g_k, g_l) et n_u unités
    jamais traitées ; effet croissant δ(k) = pente_effet·(k + 1). Isole le biais du TWFE statique dû aux comparaisons
    « tardive vs précoce déjà traitée » (Goodman-Bacon 2021) et montre que Callaway–Sant'Anna le corrige."""
    rng = np.random.default_rng(graine + 36)
    N = n_u + n_k + n_l
    alpha = rng.normal(11.0, 0.3, N); lam = np.cumsum(rng.normal(0.007, 0.004, T))
    g = np.array([np.nan] * n_u + [g_k] * n_k + [g_l] * n_l, float)
    rows = []
    for i in range(N):
        for t in range(T):
            k = t - g[i] if not np.isnan(g[i]) else np.nan
            delta = pente_effet * (k + 1) if (not np.isnan(k) and k >= 0) else 0.0
            rows.append({"unite": f"L{i:02d}", "t": t, "g": g[i], "delta_vrai": delta, "y": alpha[i] + lam[t] + delta + rng.normal(0, sigma)})
    d = pd.DataFrame(rows); d["D"] = ((d["t"] >= d["g"]) & d["g"].notna()).astype(float)
    tw = pn.twfe(d, "y", ["D"], "unite", "t")
    bac = pn.bacon(d, "y", "unite", "t", "g")
    csl = pn.callaway_santanna(d, "y", "unite", "t", "g", strate=None, controle="jamais", B=199, graine=graine)
    verite = float(d.loc[d["D"] == 1, "delta_vrai"].mean())
    return {"panel": d, "twfe": float(tw["coef"].loc["D", "coef"]), "se_grappes": float(tw["coef"].loc["D", "se_grappes"]), "bacon": bac, "cs": csl, "verite": verite,
            "verite_k": float(d.loc[(d["D"] == 1) & (d["g"] == g_k), "delta_vrai"].mean()), "verite_l": float(d.loc[(d["D"] == 1) & (d["g"] == g_l), "delta_vrai"].mean())}

# ----------------------------------------------------------------------------- solutions des exercices (petits jeux)

def exercices() -> dict:
    out = {}
    # ex. 1 : 2 × 2 à la main
    c = {"tr_av": 11.52, "tr_ap": 11.64, "te_av": 11.40, "te_ap": 11.47}
    out["ex1"] = {"did": (c["tr_ap"] - c["tr_av"]) - (c["te_ap"] - c["te_av"]), "pct": 100 * (np.exp((c["tr_ap"] - c["tr_av"]) - (c["te_ap"] - c["te_av"])) - 1)}
    # ex. 2 : within sur trois immeubles × deux périodes (y = ln RNE, x = indicatrice de rénovation)
    y = np.array([11.0, 11.3, 11.6, 11.7, 10.8, 10.9]); x = np.array([0, 1, 0, 1, 0, 0], float); u = np.array([1, 1, 2, 2, 3, 3])
    yw, Xw = pn.within(y, x[:, None], u)
    b_within = float((Xw[:, 0] @ yw) / (Xw[:, 0] @ Xw[:, 0]))
    b_pooled = float(np.polyfit(x, y, 1)[0])
    out["ex2"] = {"yw": yw, "xw": Xw[:, 0], "b_within": b_within, "b_pooled": b_pooled}
    # ex. 3 : rapport des erreurs-types (formule de Moulton approchée) : 1 + (m − 1) ρ
    out["ex3"] = {(m, rho): float(np.sqrt(1 + (m - 1) * rho)) for m, rho in [(42, 0.5), (42, 0.9), (8, 0.5)]}
    # ex. 5 : poids de Goodman-Bacon pour un plan à deux cohortes sans jamais traités, T = 20
    T = 20; gk, gl = 5, 15; nk = nl = 0.5
    Dk = (T - gk) / T; Dl = (T - gl) / T
    nkl = nk / (nk + nl)
    # variance de D doublement centré, calculée exactement sur le panel
    D = np.zeros((2, T)); D[0, gk:] = 1; D[1, gl:] = 1
    Dt = D - D.mean(1, keepdims=True) - D.mean(0, keepdims=True) + D.mean(); VD = float((Dt ** 2).mean())
    s_kl = ((nk + nl) * (1 - Dl)) ** 2 * nkl * (1 - nkl) * ((Dk - Dl) / (1 - Dl)) * ((1 - Dk) / (1 - Dl)) / VD
    s_lk = ((nk + nl) * Dk) ** 2 * nkl * (1 - nkl) * (Dl / Dk) * ((Dk - Dl) / Dk) / VD
    out["ex5"] = {"Dk": Dk, "Dl": Dl, "VD": VD, "s_kl": s_kl, "s_lk": s_lk, "somme": s_kl + s_lk}
    # ex. 6 : Callaway–Sant'Anna à la main (cohorte g = 3, base t = 2, périodes 1..4, deux traités, trois témoins)
    tr = np.array([[10.00, 10.05, 10.12, 10.20], [9.50, 9.56, 9.66, 9.72]])
    te = np.array([[10.20, 10.24, 10.27, 10.31], [9.80, 9.85, 9.89, 9.92], [9.90, 9.93, 9.98, 10.01]])
    att = {t: float((tr[:, t] - tr[:, 1]).mean() - (te[:, t] - te[:, 1]).mean()) for t in (0, 2, 3)}
    out["ex6"] = {"att": att, "moyenne_post": float(np.mean([att[2], att[3]]))}
    return out

# ----------------------------------------------------------------------------- orchestration

def main(verbose: bool = True) -> dict:
    if "res" in _CACHE: return _CACHE["res"]
    se = simuler(); v = verite(se)
    res = {"panel": se, "verite": v, "traites": choisir_traites(se.groupby("id_immeuble")[["rmr"]].first().reset_index())}
    res["did_tous"] = did_deux_periodes(se, "tous"); res["did_rmr"] = did_deux_periodes(se, "meme_rmr")
    res["placebo"] = placebo_temporel(se)
    res["twfe"] = twfe_statique(se)
    res["es"] = etude_evenement(se, "rmr"); res["es_commun"] = etude_evenement(se, None)
    res["cs"] = cs(se, "rmr"); res["cs_tous"] = cs(se, None, B=99)
    res["bacon"] = pn.bacon(se, "y", "id_immeuble", "t", "g")
    res["sc"] = synthetique(se)
    res["rdd"] = discontinuite()
    res["sens"] = sensibilite(se); res["temoins"] = groupes_temoins(se)
    res["labo"] = laboratoire()
    res["ex"] = exercices()
    _CACHE["res"] = res
    if verbose: imprimer(res)
    return res

def imprimer(res: dict) -> None:
    v = res["verite"]; tr = res["traites"]
    print("== Ch. 36 — choc simulé sur le cas 4 (données synthétiques) ==")
    print(f"Cohorte A ({RMR_A}, 2022T3) : {', '.join(tr['A'])}\nCohorte B ({RMR_B}, 2024T3) : {', '.join(tr['B'])}")
    print(f"Vérité : ATT global {v['att_global']:.4f} ; A {v['att_A']:.4f} ; B {v['att_B']:.4f} ; fenêtre 2x2 {v['att_fenetre_2x2']:.4f}")
    for nom in ("did_tous", "did_rmr", "placebo"):
        r = res[nom]; print(f"{nom}: cellules\n{r['cellules'].round(4)}\n  estimateur {r['estimateur']:.4f} ; reg {r['coef_reg']:.4f} (HC1 {r['se_hc1']:.4f}, classique {r['se_classique']:.4f}) ; n_tr {r['n_traites']} n_te {r['n_temoins']}")
    tw = res["twfe"]
    for nom in ("commun", "regional"):
        c = tw[nom]["coef"].loc["D"]; print(f"TWFE {nom}: {c['coef']:.4f} se {c['se']:.4f} hc1 {c['se_hc1']:.4f} grappes {c['se_grappes']:.4f} (k={tw[nom]['k']})")
    print(f"within deux voies : {tw['within_deux_voies']:.4f}")
    print("Event study (rmr × t) :\n", res["es"][["coef", "se_grappes", "ic_bas", "ic_haut", "vrai"]].round(4))
    c = res["cs"]; print(f"CS (même RMR) : global {c['globale']:.4f} (se {c['globale_se_boot']:.4f})\n{c['par_cohorte'].round(4)}\n{c['dynamique'].round(4)}")
    c2 = res["cs_tous"]; print(f"CS (tous témoins) : global {c2['globale']:.4f}\n{c2['par_cohorte'].round(4)}")
    b = res["bacon"]; print("Bacon :\n", b["table"].round(4), f"\n somme poids {b['somme_poids']:.6f} ; pondéré {b['beta_pondere']:.6f} ; TWFE {b['beta_twfe']:.6f}")
    s = res["sc"]; w = s["poids"][s["poids"] > 0.005].sort_values(ascending=False)
    print(f"Contrôle synthétique : effet moyen post {s['effet_moyen_post']:.4f} (vrai {s['vrai_post']:.4f}) ; RMSPE pré {s['rmspe_pre']:.4f} post {s['rmspe_post']:.4f} ratio {s['ratio_traite']:.2f} rang {s['rang_ratio']}/36\n poids :\n{w.round(3)}")
    r = res["rdd"]; print("RDD :", {hb: (round(e['tau'], 4), round(e['se_hc1'], 4), e['n_eff']) for hb, e in r["estimations"].items()}, "naïf", round(r["naif"], 4))
    sv = res["sens"]; print(f"Sensibilité (tendance +{sv['tendance']}/an) : 2x2 {sv['did_2x2']['estimateur']:.4f} ; TWFE {sv['twfe_regional']['coef'].loc['D','coef']:.4f} ; pente pré {sv['pente_pre']:.4f} ; TWFE+tend {sv['twfe_tendance']['coef'].loc['D','coef']:.4f} ; 2x2 détendancé {sv['did_detendance']:.4f}")
    print(sv["event_study"][["coef", "se_grappes"]].round(4).T)
    print("Témoins :\n", res["temoins"].round(4))
    L = res["labo"]; print(f"Laboratoire : TWFE {L['twfe']:.4f} (se {L['se_grappes']:.4f}) ; vérité {L['verite']:.4f} (k {L['verite_k']:.4f}, l {L['verite_l']:.4f}) ; CS {L['cs']['globale']:.4f} (se {L['cs']['globale_se_boot']:.4f})\n{L['cs']['par_cohorte'].round(4)}\n{L['bacon']['table'].round(4)}\n pondéré {L['bacon']['beta_pondere']:.6f} TWFE {L['bacon']['beta_twfe']:.6f}")
    print("Exercices :", {k: (vv if not isinstance(vv, dict) else {kk: (np.round(x, 4) if isinstance(x, (float, np.ndarray)) else x) for kk, x in vv.items()}) for k, vv in res["ex"].items()})

def controle_croise(res: dict) -> None:
    import statsmodels.api as sm
    from linearmodels.panel import PanelOLS
    se = res["panel"].copy()
    d = se.set_index(["id_immeuble", "t"])
    m = PanelOLS(d["y"], d[["D"]], entity_effects=True, time_effects=True).fit(cov_type="clustered", cluster_entity=True)
    mine = res["twfe"]["commun"]["coef"].loc["D"]
    print(f"linearmodels PanelOLS : D = {m.params['D']:.6f} se(grappes) {m.std_errors['D']:.6f} | maison {mine['coef']:.6f} / {mine['se_grappes']:.6f}")
    X = pd.concat([pd.DataFrame({"const": 1.0, "D": se["D"]}), pn.indicatrices(se["id_immeuble"], prefixe="u"), pn.indicatrices(se["t"], prefixe="t")], axis=1)
    o = sm.OLS(se["y"].values, X.values).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(se["id_immeuble"])[0]})
    print(f"statsmodels OLS grappes : D = {o.params[1]:.6f} se {o.bse[1]:.6f}")
    o2 = sm.OLS(se["y"].values, X.values).fit(cov_type="HC1")
    print(f"statsmodels HC1 : se {o2.bse[1]:.6f} | maison {mine['se_hc1']:.6f}")

if __name__ == "__main__":
    r = main(verbose=True)
    if "--controle" in sys.argv: controle_croise(r)