spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Chapitre 36 — Données de panel et identification causale : choc simulé sur les séries de RNE du cas 43(DONNÉES SYNTHÉTIQUES, graine 20261031).45Le panel de départ est celui de `traite.cas4.charger()` : 48 immeubles × 42 trimestres (2016T1-2026T2) de RNE trimestriel.6Ces séries sont prises comme résultats potentiels SANS traitement, Y(0). Le script y superpose un traitement déclaré :7  - cohorte A : 8 immeubles de la RMR de Montréal « à distance de marche » d'une ligne de transport collectif fictive mise8    en service au 3e trimestre de 2022 (période 26) ;9  - cohorte B : 5 immeubles de la RMR de Gatineau desservis par un second tronçon fictif au 3e trimestre de 2024 (période 34) ;10  - effet déclaré sur ln RNE, identique pour les deux cohortes, fonction du temps écoulé k depuis la mise en service :11    δ(k) = 0,01 (k + 1) pour k = 0, 1, 2, 3 (montée en régime au rythme du renouvellement des baux), puis12    δ(k) = 0,04 + 0,0025 (k − 3) pour k ≥ 4 (régime de +4 % atteint, puis valorisation lente du voisinage).13Les 35 autres immeubles ne sont jamais traités. Le lecteur connaît la vérité ; chaque estimateur lui est confronté.1415Estimateurs : DiD 2 × 2 (cellules et différences longues), TWFE statique (erreurs-types classiques / HC1 / grappes),16étude d'événement, Callaway–Sant'Anna (contrôles jamais traités, même RMR), décomposition de Goodman-Bacon, contrôle17synthétique (cohorte A agrégée vs 35 donneurs), discontinuité spatiale (simulation transversale séparée), placebo temporel,18sensibilité (tendances parallèles violées, choix du groupe témoin, traitement de l'autocorrélation).19Contrôle croisé (--controle) : linearmodels.PanelOLS et statsmodels (grappes) reproduisent le TWFE et ses erreurs-types.20Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch36_panel.py [--controle]"""21from __future__ import annotations22import sys, pathlib23import numpy as np, pandas as pd24sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))25from traite import cas426from traite import hedonique as h27from traite import panel as pn2829ROOT = pathlib.Path(__file__).resolve().parents[3]30GRAINE = 2026103131T_A, T_B = 26, 34                      # 2022T3 et 2024T3 (2016T1 = 0)32N_A, N_B = 8, 5                        # immeubles traités à Montréal et à Gatineau33RMR_A, RMR_B = "Montréal", "Gatineau"34FENETRE_AVANT = [22, 23, 24, 25]       # 2021T3-2022T235FENETRE_APRES = [30, 31, 32, 33]       # 2023T3-2024T2 (deuxième année : régime atteint)36PLACEBO_T = 14                         # 2019T337B_BOOT = 4993839_CACHE: dict = {}4041def effet_vrai(k) -> np.ndarray:42    """Effet déclaré δ(k) sur ln RNE selon le temps relatif k (0 avant la mise en service)."""43    k = np.asarray(k, float)44    return np.where(k < 0, 0.0, np.where(k <= 3, 0.01 * (k + 1), 0.04 + 0.0025 * (k - 3)))4546# ----------------------------------------------------------------------------- panel simulé4748def panel_base() -> pd.DataFrame:49    d = cas4.charger(); se = d["series"].copy(); pf = d["portefeuille"]50    trims = sorted(se["trimestre"].unique()); idx = {p: i for i, p in enumerate(trims)}51    se["t"] = se["trimestre"].map(idx).astype(int)52    se = se.merge(pf[["id_immeuble", "rmr", "nb_logements", "zone_inondable"]], on="id_immeuble")53    se["y0"] = np.log(se["rne_trimestriel"].astype(float))54    return se.sort_values(["id_immeuble", "t"]).reset_index(drop=True)5556def choisir_traites(pf_ids: pd.DataFrame, graine: int = GRAINE) -> dict:57    rng = np.random.default_rng(graine)58    a = sorted(rng.choice(sorted(pf_ids.loc[pf_ids["rmr"] == RMR_A, "id_immeuble"]), N_A, replace=False))59    b = sorted(rng.choice(sorted(pf_ids.loc[pf_ids["rmr"] == RMR_B, "id_immeuble"]), N_B, replace=False))60    return {"A": list(a), "B": list(b)}6162def simuler(tendance_diff_annuelle: float = 0.0) -> pd.DataFrame:63    """Panel long avec traitement échelonné. `tendance_diff_annuelle` > 0 ajoute aux immeubles traités une croissance64    supplémentaire (par an, en log) sur toute la période : violation déclarée des tendances parallèles."""65    se = panel_base(); ids = se.groupby("id_immeuble")[["rmr"]].first().reset_index()66    tr = choisir_traites(ids)67    se["g"] = np.nan68    se.loc[se["id_immeuble"].isin(tr["A"]), "g"] = T_A69    se.loc[se["id_immeuble"].isin(tr["B"]), "g"] = T_B70    se["cohorte"] = np.where(se["id_immeuble"].isin(tr["A"]), "A", np.where(se["id_immeuble"].isin(tr["B"]), "B", "jamais"))71    se["k"] = (se["t"] - se["g"]).where(se["g"].notna())72    se["D"] = ((se["k"] >= 0) & se["g"].notna()).astype(float)73    se["delta_vrai"] = np.where(se["g"].notna(), effet_vrai(se["k"].fillna(-1)), 0.0)74    if tendance_diff_annuelle:75        se.loc[se["g"].notna(), "y0"] = se.loc[se["g"].notna(), "y0"] + tendance_diff_annuelle * se.loc[se["g"].notna(), "t"] / 4.076    se["y"] = se["y0"] + se["delta_vrai"]77    se["traite_A"] = se["id_immeuble"].isin(tr["A"]).astype(int)78    return se7980def verite(se: pd.DataFrame) -> dict:81    post = se[se["D"] == 1]82    return {"att_global": float(post["delta_vrai"].mean()),83            "att_A": float(post.loc[post["cohorte"] == "A", "delta_vrai"].mean()),84            "att_B": float(post.loc[post["cohorte"] == "B", "delta_vrai"].mean()),85            "att_fenetre_2x2": float(effet_vrai(np.array(FENETRE_APRES) - T_A).mean()),86            "profil": {int(k): float(effet_vrai(k)) for k in range(0, 16)}}8788# ----------------------------------------------------------------------------- estimations8990def did_deux_periodes(se: pd.DataFrame, temoins: str = "tous") -> dict:91    """2 × 2 : cohorte A vs jamais traités (tous, ou de la même RMR)."""92    d = se[se["cohorte"] != "B"].copy()93    if temoins == "meme_rmr":94        d = d[d["rmr"] == RMR_A]95    return pn.did_2x2(d, "y", "id_immeuble", "t", "traite_A", FENETRE_AVANT, FENETRE_APRES)9697def placebo_temporel(se: pd.DataFrame) -> dict:98    d = se[(se["cohorte"] != "B") & (se["rmr"] == RMR_A)].copy()99    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]100    return pn.did_2x2(d, "y", "id_immeuble", "t", "traite_A", avant, apres)101102def twfe_statique(se: pd.DataFrame) -> dict:103    commun = pn.twfe(se, "y", ["D"], "id_immeuble", "t")104    regional = pn.twfe(se, "y", ["D"], "id_immeuble", "t", fe_supp="rmr")105    # équivalence within : D doublement centré (cylindré) → même coefficient que les indicatrices106    yw, Xw = pn.within_deux_voies(se["y"].values, se[["D"]].values, se["id_immeuble"].values, se["t"].values)107    bw = float(np.linalg.lstsq(Xw, yw, rcond=None)[0][0])108    return {"commun": commun, "regional": regional, "within_deux_voies": bw}109110def etude_evenement(se: pd.DataFrame, fe_supp: str | None = "rmr") -> pd.DataFrame:111    es = pn.event_study(se, "y", "id_immeuble", "t", "g", fenetre=(-8, 8), ref=-1, fe_supp=fe_supp)112    # vérité par k : moyenne de δ sur les observations traitées regroupées dans chaque indicatrice113    tr = se[se["g"].notna()].copy(); tr["kb"] = tr["k"].clip(-8, 8)114    v = tr.groupby("kb")["delta_vrai"].mean()115    es["vrai"] = [float(v.get(k, 0.0)) if k >= 0 else 0.0 for k in es.index]116    return es117118def cs(se: pd.DataFrame, strate: str | None = "rmr", B: int = B_BOOT) -> dict:119    return pn.callaway_santanna(se, "y", "id_immeuble", "t", "g", strate=strate, B=B)120121def synthetique(se: pd.DataFrame) -> dict:122    """Unité traitée = moyenne des ln RNE de la cohorte A, centrée sur sa moyenne pré-traitement ; donneurs = 35 jamais traités,123    chacun centré de même (contrôle synthétique « en écarts à la moyenne pré-traitement »)."""124    W = se.pivot(index="id_immeuble", columns="t", values="y")125    coh = se.groupby("id_immeuble")["cohorte"].first()126    A = W.loc[coh[coh == "A"].index].mean(axis=0).values127    donneurs = W.loc[coh[coh == "jamais"].index]128    T0 = T_A129    A_c = A - A[:T0].mean()130    Y0 = (donneurs.values - donneurs.values[:, :T0].mean(axis=1, keepdims=True)).T   # (T × J)131    r = pn.controle_synthetique(Y0, A_c, T0)132    r["poids"] = pd.Series(r["poids"], index=donneurs.index)133    r["traite"] = A_c; r["T0"] = T0134    r["vrai_post"] = float(effet_vrai(np.arange(0, W.shape[1] - T0)).mean())135    r["placebos"] = pn.placebos_synthetiques(Y0, T0); r["placebos"]["donneur"] = donneurs.index[r["placebos"]["donneur"].values]136    r["ratio_traite"] = r["rmspe_post"] / r["rmspe_pre"]137    r["rang_ratio"] = int((r["placebos"]["ratio"] >= r["ratio_traite"]).sum() + 1)138    return r139140def discontinuite(graine: int = GRAINE, n: int = 400, saut: float = 0.05) -> dict:141    """Simulation transversale (fictive) : n immeubles répartis de part et d'autre d'une frontière de zone (x = distance signée en km,142    positive du côté de la zone dense), ln V = 12,9 + 0,03 x + saut·1[x ≥ 0] + ε, ε ~ N(0 ; 0,08²)."""143    rng = np.random.default_rng(graine + 7)144    x = rng.uniform(-2, 2, n); y = 12.9 + 0.03 * x + saut * (x >= 0) + rng.normal(0, 0.08, n)145    res = {hb: pn.rdd_local_lineaire(x, y, hb) for hb in (0.5, 1.0, 2.0)}146    naif = float(y[x >= 0].mean() - y[x < 0].mean())147    return {"x": x, "y": y, "saut_vrai": saut, "estimations": res, "naif": naif}148149def sensibilite(se: pd.DataFrame, tendance: float = 0.01) -> dict:150    """Tendances parallèles violées : +`tendance` par an (log) pour les immeubles traités depuis 2016."""151    sv = simuler(tendance_diff_annuelle=tendance)152    did_v = pn.did_2x2(sv[(sv["cohorte"] != "B") & (sv["rmr"] == RMR_A)], "y", "id_immeuble", "t", "traite_A", FENETRE_AVANT, FENETRE_APRES)153    tw_v = pn.twfe(sv, "y", ["D"], "id_immeuble", "t", fe_supp="rmr")154    es_v = pn.event_study(sv, "y", "id_immeuble", "t", "g", fenetre=(-8, 8), ref=-1, fe_supp="rmr")155    # pente des coefficients pré-traitement (k = −8 … −2) : diagnostic156    pre = es_v.loc[[k for k in es_v.index if -8 < k < -1]]157    pente_pre = float(np.polyfit(pre.index.values.astype(float), pre["coef"].values, 1)[0])158    # correction par tendances linéaires propres aux unités traitées (groupe × tendance)159    sv2 = sv.copy(); sv2["tend_tr"] = sv2["g"].notna().astype(float) * sv2["t"] / 4.0160    tw_tend = pn.twfe(sv2, "y", ["D", "tend_tr"], "id_immeuble", "t", fe_supp="rmr")161    # correction par extrapolation de la tendance pré-traitement (2 × 2 détendancé) : Δy_i − pente_i × écart temporel162    d = sv[(sv["cohorte"] != "B") & (sv["rmr"] == RMR_A) & (sv["t"] < T_A)]163    pentes = d.groupby("id_immeuble").apply(lambda g: np.polyfit(g["t"] / 4.0, g["y"], 1)[0], include_groups=False)164    deltas = did_v["deltas"].set_index("id_immeuble")165    ecart_ans = (np.mean(FENETRE_APRES) - np.mean(FENETRE_AVANT)) / 4.0166    deltas["dy_det"] = deltas["dy"] - pentes.reindex(deltas.index).values * ecart_ans167    did_det = float(deltas.loc[deltas["traite_A"] == 1, "dy_det"].mean() - deltas.loc[deltas["traite_A"] == 0, "dy_det"].mean())168    return {"tendance": tendance, "did_2x2": did_v, "twfe_regional": tw_v, "event_study": es_v, "pente_pre": pente_pre,169            "twfe_tendance": tw_tend, "did_detendance": did_det, "ecart_ans": ecart_ans}170171def groupes_temoins(se: pd.DataFrame) -> pd.DataFrame:172    """Sensibilité au choix des témoins pour la cohorte A, 2 × 2 sur les fenêtres de référence."""173    rows = []174    for nom, d in [("jamais traités, toutes RMR", se[se["cohorte"] != "B"]),175                   ("jamais traités, Montréal", se[(se["cohorte"] != "B") & (se["rmr"] == RMR_A)]),176                   ("jamais traités, hors Montréal", se[(se["cohorte"] != "B") & ((se["rmr"] != RMR_A) | (se["traite_A"] == 1))]),177                   ("pas encore traités (cohorte B)", se[(se["cohorte"] == "B") | (se["traite_A"] == 1)])]:178        r = pn.did_2x2(d, "y", "id_immeuble", "t", "traite_A", FENETRE_AVANT, FENETRE_APRES)179        rows.append({"temoins": nom, "n_temoins": r["n_temoins"], "estimateur": r["estimateur"], "se_hc1": r["se_hc1"]})180    return pd.DataFrame(rows)181182# ----------------------------------------------------------------------------- laboratoire : biais du TWFE échelonné183184def 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,185                pente_effet: float = 0.005, sigma: float = 0.02) -> dict:186    """Panel artificiel à tendance commune (α_i + λ_t + ε, sans tendance régionale) : deux cohortes (g_k, g_l) et n_u unités187    jamais traitées ; effet croissant δ(k) = pente_effet·(k + 1). Isole le biais du TWFE statique dû aux comparaisons188    « tardive vs précoce déjà traitée » (Goodman-Bacon 2021) et montre que Callaway–Sant'Anna le corrige."""189    rng = np.random.default_rng(graine + 36)190    N = n_u + n_k + n_l191    alpha = rng.normal(11.0, 0.3, N); lam = np.cumsum(rng.normal(0.007, 0.004, T))192    g = np.array([np.nan] * n_u + [g_k] * n_k + [g_l] * n_l, float)193    rows = []194    for i in range(N):195        for t in range(T):196            k = t - g[i] if not np.isnan(g[i]) else np.nan197            delta = pente_effet * (k + 1) if (not np.isnan(k) and k >= 0) else 0.0198            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)})199    d = pd.DataFrame(rows); d["D"] = ((d["t"] >= d["g"]) & d["g"].notna()).astype(float)200    tw = pn.twfe(d, "y", ["D"], "unite", "t")201    bac = pn.bacon(d, "y", "unite", "t", "g")202    csl = pn.callaway_santanna(d, "y", "unite", "t", "g", strate=None, controle="jamais", B=199, graine=graine)203    verite = float(d.loc[d["D"] == 1, "delta_vrai"].mean())204    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,205            "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())}206207# ----------------------------------------------------------------------------- solutions des exercices (petits jeux)208209def exercices() -> dict:210    out = {}211    # ex. 1 : 2 × 2 à la main212    c = {"tr_av": 11.52, "tr_ap": 11.64, "te_av": 11.40, "te_ap": 11.47}213    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)}214    # ex. 2 : within sur trois immeubles × deux périodes (y = ln RNE, x = indicatrice de rénovation)215    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])216    yw, Xw = pn.within(y, x[:, None], u)217    b_within = float((Xw[:, 0] @ yw) / (Xw[:, 0] @ Xw[:, 0]))218    b_pooled = float(np.polyfit(x, y, 1)[0])219    out["ex2"] = {"yw": yw, "xw": Xw[:, 0], "b_within": b_within, "b_pooled": b_pooled}220    # ex. 3 : rapport des erreurs-types (formule de Moulton approchée) : 1 + (m − 1) ρ221    out["ex3"] = {(m, rho): float(np.sqrt(1 + (m - 1) * rho)) for m, rho in [(42, 0.5), (42, 0.9), (8, 0.5)]}222    # ex. 5 : poids de Goodman-Bacon pour un plan à deux cohortes sans jamais traités, T = 20223    T = 20; gk, gl = 5, 15; nk = nl = 0.5224    Dk = (T - gk) / T; Dl = (T - gl) / T225    nkl = nk / (nk + nl)226    # variance de D doublement centré, calculée exactement sur le panel227    D = np.zeros((2, T)); D[0, gk:] = 1; D[1, gl:] = 1228    Dt = D - D.mean(1, keepdims=True) - D.mean(0, keepdims=True) + D.mean(); VD = float((Dt ** 2).mean())229    s_kl = ((nk + nl) * (1 - Dl)) ** 2 * nkl * (1 - nkl) * ((Dk - Dl) / (1 - Dl)) * ((1 - Dk) / (1 - Dl)) / VD230    s_lk = ((nk + nl) * Dk) ** 2 * nkl * (1 - nkl) * (Dl / Dk) * ((Dk - Dl) / Dk) / VD231    out["ex5"] = {"Dk": Dk, "Dl": Dl, "VD": VD, "s_kl": s_kl, "s_lk": s_lk, "somme": s_kl + s_lk}232    # ex. 6 : Callaway–Sant'Anna à la main (cohorte g = 3, base t = 2, périodes 1..4, deux traités, trois témoins)233    tr = np.array([[10.00, 10.05, 10.12, 10.20], [9.50, 9.56, 9.66, 9.72]])234    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]])235    att = {t: float((tr[:, t] - tr[:, 1]).mean() - (te[:, t] - te[:, 1]).mean()) for t in (0, 2, 3)}236    out["ex6"] = {"att": att, "moyenne_post": float(np.mean([att[2], att[3]]))}237    return out238239def exercice_7(se: pd.DataFrame) -> pd.DataFrame:240    """Contrôle synthétique de la cohorte A selon le bassin de donneurs, et placebo temporel (mise en service fictive 2020T3)."""241    W = se.pivot(index="id_immeuble", columns="t", values="y"); coh = se.groupby("id_immeuble")["cohorte"].first()242    rmr = se.groupby("id_immeuble")["rmr"].first()243    A = W.loc[coh[coh == "A"].index].mean(axis=0).values244    rows = []245    for nom, ids, T0 in [("35 donneurs, six RMR", coh[coh == "jamais"].index, T_A),246                         ("9 donneurs, Montréal", coh[(coh == "jamais") & (rmr == RMR_A)].index, T_A),247                         ("26 donneurs, hors Montréal", coh[(coh == "jamais") & (rmr != RMR_A)].index, T_A),248                         ("placebo 2020T3, 35 donneurs", coh[coh == "jamais"].index, 18)]:249        D = W.loc[ids]; Y0 = (D.values - D.values[:, :T0].mean(axis=1, keepdims=True)).T; A_c = A - A[:T0].mean()250        r = pn.controle_synthetique(Y0, A_c, T0)251        fin = T_A if T0 < T_A else W.shape[1]          # placebo : écart moyen jusqu'à la vraie mise en service252        effet = float(r["ecarts"][T0:fin].mean())253        w = pd.Series(r["poids"], index=ids); top = ", ".join(f"{i} {v:.2f}" for i, v in w.sort_values(ascending=False).head(3).items())254        rows.append({"variante": nom, "n_donneurs": len(ids), "rmspe_pre": r["rmspe_pre"], "effet_post": effet, "n_poids_sup_1pct": int((w > 0.01).sum()), "trois_premiers": top})255    return pd.DataFrame(rows)256257# ----------------------------------------------------------------------------- orchestration258259def main(verbose: bool = True) -> dict:260    if "res" in _CACHE: return _CACHE["res"]261    se = simuler(); v = verite(se)262    res = {"panel": se, "verite": v, "traites": choisir_traites(se.groupby("id_immeuble")[["rmr"]].first().reset_index())}263    res["did_tous"] = did_deux_periodes(se, "tous"); res["did_rmr"] = did_deux_periodes(se, "meme_rmr")264    res["placebo"] = placebo_temporel(se)265    res["twfe"] = twfe_statique(se)266    res["es"] = etude_evenement(se, "rmr"); res["es_commun"] = etude_evenement(se, None)267    res["cs"] = cs(se, "rmr"); res["cs_tous"] = cs(se, None, B=99)268    res["bacon"] = pn.bacon(se, "y", "id_immeuble", "t", "g")269    res["sc"] = synthetique(se)270    res["rdd"] = discontinuite()271    res["sens"] = sensibilite(se); res["temoins"] = groupes_temoins(se)272    res["labo"] = laboratoire()273    res["ex"] = exercices(); res["ex7"] = exercice_7(se)274    _CACHE["res"] = res275    if verbose: imprimer(res)276    return res277278def imprimer(res: dict) -> None:279    v = res["verite"]; tr = res["traites"]280    print("== Ch. 36 — choc simulé sur le cas 4 (données synthétiques) ==")281    print(f"Cohorte A ({RMR_A}, 2022T3) : {', '.join(tr['A'])}\nCohorte B ({RMR_B}, 2024T3) : {', '.join(tr['B'])}")282    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}")283    for nom in ("did_tous", "did_rmr", "placebo"):284        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']}")285    tw = res["twfe"]286    for nom in ("commun", "regional"):287        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']})")288    print(f"within deux voies : {tw['within_deux_voies']:.4f}")289    print("Event study (rmr × t) :\n", res["es"][["coef", "se_grappes", "ic_bas", "ic_haut", "vrai"]].round(4))290    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)}")291    c2 = res["cs_tous"]; print(f"CS (tous témoins) : global {c2['globale']:.4f}\n{c2['par_cohorte'].round(4)}")292    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}")293    s = res["sc"]; w = s["poids"][s["poids"] > 0.005].sort_values(ascending=False)294    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)}")295    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))296    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}")297    print(sv["event_study"][["coef", "se_grappes"]].round(4).T)298    print("Témoins :\n", res["temoins"].round(4))299    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}")300    print("Exercice 7 :\n", res["ex7"].round(4).to_string())301    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()})302303def controle_croise(res: dict) -> None:304    import statsmodels.api as sm305    from linearmodels.panel import PanelOLS306    se = res["panel"].copy()307    d = se.set_index(["id_immeuble", "t"])308    m = PanelOLS(d["y"], d[["D"]], entity_effects=True, time_effects=True).fit(cov_type="clustered", cluster_entity=True)309    mine = res["twfe"]["commun"]["coef"].loc["D"]310    print(f"linearmodels PanelOLS : D = {m.params['D']:.6f} se(grappes) {m.std_errors['D']:.6f} | maison {mine['coef']:.6f} / {mine['se_grappes']:.6f}")311    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)312    o = sm.OLS(se["y"].values, X.values).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(se["id_immeuble"])[0]})313    print(f"statsmodels OLS grappes : D = {o.params[1]:.6f} se {o.bse[1]:.6f}")314    o2 = sm.OLS(se["y"].values, X.values).fit(cov_type="HC1")315    print(f"statsmodels HC1 : se {o2.bse[1]:.6f} | maison {mine['se_hc1']:.6f}")316317if __name__ == "__main__":318    r = main(verbose=True)319    if "--controle" in sys.argv: controle_croise(r)320