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# ----------------------------------------------------------------------------- solutions des exercices (petits jeux)183184def exercices() -> dict:185    out = {}186    # ex. 1 : 2 × 2 à la main187    c = {"tr_av": 11.52, "tr_ap": 11.64, "te_av": 11.40, "te_ap": 11.47}188    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)}189    # ex. 2 : within sur trois immeubles × deux périodes (y = ln RNE, x = indicatrice de rénovation)190    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])191    yw, Xw = pn.within(y, x[:, None], u)192    b_within = float((Xw[:, 0] @ yw) / (Xw[:, 0] @ Xw[:, 0]))193    b_pooled = float(np.polyfit(x, y, 1)[0])194    out["ex2"] = {"yw": yw, "xw": Xw[:, 0], "b_within": b_within, "b_pooled": b_pooled}195    # ex. 3 : rapport des erreurs-types (formule de Moulton approchée) : 1 + (m − 1) ρ196    out["ex3"] = {m: float(np.sqrt(1 + (m - 1) * rho)) for m, rho in [(42, 0.5), (42, 0.9), (8, 0.5)]}197    # ex. 5 : poids de Goodman-Bacon pour un plan à deux cohortes sans jamais traités, T = 20198    T = 20; gk, gl = 5, 15; nk = nl = 0.5199    Dk = (T - gk) / T; Dl = (T - gl) / T200    nkl = nk / (nk + nl)201    # variance de D doublement centré, calculée exactement sur le panel202    D = np.zeros((2, T)); D[0, gk:] = 1; D[1, gl:] = 1203    Dt = D - D.mean(1, keepdims=True) - D.mean(0, keepdims=True) + D.mean(); VD = float((Dt ** 2).mean())204    s_kl = ((nk + nl) * (1 - Dl)) ** 2 * nkl * (1 - nkl) * ((Dk - Dl) / (1 - Dl)) * ((1 - Dk) / (1 - Dl)) / VD205    s_lk = ((nk + nl) * Dk) ** 2 * nkl * (1 - nkl) * (Dl / Dk) * ((Dk - Dl) / Dk) / VD206    out["ex5"] = {"Dk": Dk, "Dl": Dl, "VD": VD, "s_kl": s_kl, "s_lk": s_lk, "somme": s_kl + s_lk}207    # ex. 6 : Callaway–Sant'Anna à la main (cohorte g = 3, base t = 2, périodes 1..4, deux traités, trois témoins)208    tr = np.array([[10.00, 10.05, 10.12, 10.20], [9.50, 9.56, 9.66, 9.72]])209    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]])210    att = {t: float((tr[:, t] - tr[:, 1]).mean() - (te[:, t] - te[:, 1]).mean()) for t in (0, 2, 3)}211    out["ex6"] = {"att": att, "moyenne_post": float(np.mean([att[2], att[3]]))}212    return out213214# ----------------------------------------------------------------------------- orchestration215216def main(verbose: bool = True) -> dict:217    if "res" in _CACHE: return _CACHE["res"]218    se = simuler(); v = verite(se)219    res = {"panel": se, "verite": v, "traites": choisir_traites(se.groupby("id_immeuble")[["rmr"]].first().reset_index())}220    res["did_tous"] = did_deux_periodes(se, "tous"); res["did_rmr"] = did_deux_periodes(se, "meme_rmr")221    res["placebo"] = placebo_temporel(se)222    res["twfe"] = twfe_statique(se)223    res["es"] = etude_evenement(se, "rmr"); res["es_commun"] = etude_evenement(se, None)224    res["cs"] = cs(se, "rmr"); res["cs_tous"] = cs(se, None, B=99)225    res["bacon"] = pn.bacon(se, "y", "id_immeuble", "t", "g")226    res["sc"] = synthetique(se)227    res["rdd"] = discontinuite()228    res["sens"] = sensibilite(se); res["temoins"] = groupes_temoins(se)229    res["ex"] = exercices()230    _CACHE["res"] = res231    if verbose: imprimer(res)232    return res233234def imprimer(res: dict) -> None:235    v = res["verite"]; tr = res["traites"]236    print("== Ch. 36 — choc simulé sur le cas 4 (données synthétiques) ==")237    print(f"Cohorte A ({RMR_A}, 2022T3) : {', '.join(tr['A'])}\nCohorte B ({RMR_B}, 2024T3) : {', '.join(tr['B'])}")238    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}")239    for nom in ("did_tous", "did_rmr", "placebo"):240        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']}")241    tw = res["twfe"]242    for nom in ("commun", "regional"):243        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']})")244    print(f"within deux voies : {tw['within_deux_voies']:.4f}")245    print("Event study (rmr × t) :\n", res["es"][["coef", "se_grappes", "ic_bas", "ic_haut", "vrai"]].round(4))246    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)}")247    c2 = res["cs_tous"]; print(f"CS (tous témoins) : global {c2['globale']:.4f}\n{c2['par_cohorte'].round(4)}")248    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}")249    s = res["sc"]; w = s["poids"][s["poids"] > 0.005].sort_values(ascending=False)250    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)}")251    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))252    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}")253    print(sv["event_study"][["coef", "se_grappes"]].round(4).T)254    print("Témoins :\n", res["temoins"].round(4))255    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()})256257def controle_croise(res: dict) -> None:258    import statsmodels.api as sm259    from linearmodels.panel import PanelOLS260    se = res["panel"].copy()261    d = se.set_index(["id_immeuble", "t"])262    m = PanelOLS(d["y"], d[["D"]], entity_effects=True, time_effects=True).fit(cov_type="clustered", cluster_entity=True)263    mine = res["twfe"]["commun"]["coef"].loc["D"]264    print(f"linearmodels PanelOLS : D = {m.params['D']:.6f} se(grappes) {m.std_errors['D']:.6f} | maison {mine['coef']:.6f} / {mine['se_grappes']:.6f}")265    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)266    o = sm.OLS(se["y"].values, X.values).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(se["id_immeuble"])[0]})267    print(f"statsmodels OLS grappes : D = {o.params[1]:.6f} se {o.bse[1]:.6f}")268    o2 = sm.OLS(se["y"].values, X.values).fit(cov_type="HC1")269    print(f"statsmodels HC1 : se {o2.bse[1]:.6f} | maison {mine['se_hc1']:.6f}")270271if __name__ == "__main__":272    r = main(verbose=True)273    if "--controle" in sys.argv: controle_croise(r)274