spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Figures du chapitre 36 — données de panel et identification causale (cas 4, DONNÉES SYNTHÉTIQUES, graine 20261031).3Sorties (assets/figures/ch36/) :4  ch36-fig-trajectoires.pdf    ln RNE moyen (cumul depuis 2016T1) des cohortes traitées et des témoins de la même RMR5  ch36-fig-event-study.pdf     étude d'événement TWFE (RMR × trimestre) et agrégation dynamique de Callaway–Sant'Anna, avec la vérité6  ch36-fig-bacon.pdf           décomposition de Goodman-Bacon : cas 4 et laboratoire à effets croissants7  ch36-fig-synthetique.pdf     contrôle synthétique de la cohorte A : trajectoires, écarts et placebos, poids des donneurs8  ch36-fig-discontinuite.pdf   discontinuité spatiale à une frontière de zone (simulation transversale)9  ch36-fig-sensibilite.pdf     étude d'événement lorsque les tendances parallèles sont violées (+1 pt/an pour les traités)10Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch36_figures.py"""11import sys, pathlib12ROOT = pathlib.Path(__file__).resolve().parents[3]13sys.path.insert(0, str(ROOT / "analysis/python")); sys.path.insert(0, str(ROOT / "analysis/python/cases"))14import numpy as np, pandas as pd15import matplotlib.pyplot as plt16from traite.style import setup, palette, savefig, LARGEUR17import ch36_panel as c361819TRIMS = pd.period_range("2016Q1", "2026Q2", freq="Q")20def _etiq(ax, pas=8):21    ax.set_xticks(range(0, len(TRIMS), pas)); ax.set_xticklabels([str(TRIMS[i]).replace("Q", "T") for i in range(0, len(TRIMS), pas)], fontsize=7.5)2223def fig_trajectoires(res):24    se = res["panel"]25    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.0), sharey=True)26    for ax, (rmr, coh, T0, titre) in zip(axes, [(c36.RMR_A, "A", c36.T_A, "Montréal : cohorte A (2022T3)"), (c36.RMR_B, "B", c36.T_B, "Gatineau : cohorte B (2024T3)")]):27        for i, (lab, m) in enumerate([("Immeubles traités", (se["cohorte"] == coh)), ("Témoins jamais traités, même RMR", (se["cohorte"] == "jamais") & (se["rmr"] == rmr))]):28            W = se[m].pivot(index="id_immeuble", columns="t", values="y")29            traj = (W.sub(W[0], axis=0)).mean(axis=0)30            ax.plot(traj.index, traj.values, color=palette.cat[i], lw=1.7, label=f"{lab} (n = {W.shape[0]})")31            # contrefactuel déclaré : trajectoire traitée moins l'effet vrai moyen32            if i == 0:33                d = se[m].groupby("t")["delta_vrai"].mean()34                ax.plot(traj.index, traj.values - d.values, color=palette.cat[0], lw=1.0, ls=":", label="Traités sans traitement, Y(0) (connu ici)")35        ax.axvline(T0, color=palette.ink2, lw=0.9, ls="--"); ax.text(T0 + 0.4, 0.0, "mise en service", fontsize=7, color=palette.ink2, rotation=90, va="bottom")36        ax.set_title(titre, fontsize=8.5); _etiq(ax); ax.grid(True, lw=0.4)37    axes[0].set_ylabel("ln RNE moyen, écart à 2016T1"); axes[0].legend(fontsize=6.8, loc="upper left")38    savefig(fig, "ch36/ch36-fig-trajectoires.pdf")3940def fig_event_study(res):41    es = res["es"]; cs = res["cs"]["dynamique"].loc[-8:15]42    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.1), sharey=True)43    ax = axes[0]44    ax.fill_between(es.index, es["ic_bas"], es["ic_haut"], color=palette.cat[0], alpha=0.18, lw=0)45    ax.plot(es.index, es["coef"], "o-", color=palette.cat[0], ms=3.8, lw=1.4, label="Coefficient estimé (IC à 95 %, grappes)")46    ax.step(es.index, es["vrai"], where="mid", color=palette.ink, lw=1.1, ls="--", label="Effet déclaré δ(k)")47    ax.set_title("Étude d'événement TWFE (RMR × trimestre), k de −8 à +8 regroupés", fontsize=8); ax.set_xlabel("Trimestres depuis la mise en service, k")48    ax.set_ylabel("Effet sur ln RNE (points de log)")49    ax = axes[1]50    ax.fill_between(cs.index, cs["ic_bas"], cs["ic_haut"], color=palette.cat[1], alpha=0.18, lw=0)51    ax.plot(cs.index, cs["att"], "s-", color=palette.cat[1], ms=3.4, lw=1.4, label="ATT(k) de Callaway–Sant'Anna (IC bootstrap)")52    kk = np.arange(-8, 16); ax.step(kk, c36.effet_vrai(kk), where="mid", color=palette.ink, lw=1.1, ls="--", label="Effet déclaré δ(k)")53    ax.set_title("Agrégation dynamique de Callaway–Sant'Anna, témoins de la même RMR", fontsize=8); ax.set_xlabel("Trimestres depuis la mise en service, k")54    for ax in axes:55        ax.axhline(0, color=palette.ink2, lw=0.8); ax.axvline(-0.5, color=palette.ink2, lw=0.8, ls=":"); ax.legend(fontsize=6.8, loc="upper left"); ax.grid(True, lw=0.4)56    savefig(fig, "ch36/ch36-fig-event-study.pdf")5758def fig_bacon(res):59    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.0))60    marqueurs = {"traités vs jamais traités": "o", "précoce vs tardive (pas encore traitée)": "^", "tardive vs précoce (déjà traitée)": "D"}61    for ax, (b, vrai, titre) in zip(axes, [(res["bacon"], res["verite"]["att_global"], "Cas 4 : cohortes A et B, 35 témoins"), (res["labo"]["bacon"], res["labo"]["verite"], "Laboratoire : effets croissants, 10 témoins")]):62        t = b["table"]63        for j, (typ, mk) in enumerate(marqueurs.items()):64            s = t[t["type"] == typ]65            ax.scatter(s["poids"], s["beta"], marker=mk, s=55, color=palette.cat[j], edgecolor=palette.ink, lw=0.5, label=typ, zorder=5)66        ax.axhline(b["beta_twfe"], color=palette.ink, lw=1.1, label=f"TWFE = Σ poids × 2×2 = {b['beta_twfe']:.3f}")67        ax.axhline(vrai, color=palette.ink2, lw=1.0, ls="--", label=f"Effet moyen vrai = {vrai:.3f}")68        ax.axhline(0, color=palette.grid, lw=0.8)69        ax.set_xlabel("Poids de la comparaison 2 × 2"); ax.set_title(titre, fontsize=8.5); ax.grid(True, lw=0.4)70        ax.legend(fontsize=6.3, loc="best")71    axes[0].set_ylabel("Estimation 2 × 2 (points de log)")72    savefig(fig, "ch36/ch36-fig-bacon.pdf")7374def fig_synthetique(res):75    s = res["sc"]; T0 = s["T0"]; t = np.arange(len(s["traite"]))76    fig, axes = plt.subplots(1, 3, figsize=(LARGEUR, 2.9), gridspec_kw={"width_ratios": [1.3, 1.3, 1.0]})77    ax = axes[0]78    ax.plot(t, s["traite"], color=palette.cat[0], lw=1.7, label="Cohorte A (moyenne des 8 immeubles)")79    ax.plot(t, s["synthetique"], color=palette.cat[1], lw=1.5, ls="--", label="Contrôle synthétique")80    ax.plot(t, s["traite"] - np.concatenate([np.zeros(T0), c36.effet_vrai(np.arange(len(t) - T0))]), color=palette.cat[0], lw=0.9, ls=":", label="Y(0) déclaré")81    ax.axvline(T0 - 0.5, color=palette.ink2, lw=0.9, ls="--"); ax.set_title("Trajectoires (ln RNE centré sur la moyenne pré)", fontsize=8); _etiq(ax, 10); ax.legend(fontsize=6.3, loc="upper left"); ax.grid(True, lw=0.4)82    ax = axes[1]83    W = res["panel"].pivot(index="id_immeuble", columns="t", values="y"); coh = res["panel"].groupby("id_immeuble")["cohorte"].first()84    donneurs = W.loc[coh[coh == "jamais"].index]; Y0 = (donneurs.values - donneurs.values[:, :T0].mean(axis=1, keepdims=True)).T85    from traite import panel as pn86    for j in range(Y0.shape[1]):87        r = pn.controle_synthetique(np.delete(Y0, j, axis=1), Y0[:, j], T0)88        ax.plot(t, r["ecarts"], color=palette.grid, lw=0.7, alpha=0.9, label="Placebos (35 donneurs)" if j == 0 else None)89    ax.plot(t, s["ecarts"], color=palette.cat[0], lw=1.8, label="Cohorte A")90    ax.plot(t, np.concatenate([np.zeros(T0), c36.effet_vrai(np.arange(len(t) - T0))]), color=palette.ink, lw=1.0, ls="--", label="Effet déclaré")91    ax.axvline(T0 - 0.5, color=palette.ink2, lw=0.9, ls="--"); ax.axhline(0, color=palette.ink2, lw=0.7)92    ax.set_title("Écarts traité − synthétique et placebos", fontsize=8); _etiq(ax, 10); ax.legend(fontsize=6.3, loc="upper left"); ax.grid(True, lw=0.4); ax.set_ylim(-0.12, 0.12)93    ax = axes[2]94    w = s["poids"].sort_values(ascending=True); w = w[w > 0.005]95    ax.barh(range(len(w)), w.values, color=palette.cat[1], height=0.7)96    ax.set_yticks(range(len(w))); ax.set_yticklabels(w.index, fontsize=7); ax.set_xlabel("Poids"); ax.set_title(f"Poids > 0,5 % ({len(w)} donneurs)", fontsize=8); ax.grid(True, axis="x", lw=0.4)97    savefig(fig, "ch36/ch36-fig-synthetique.pdf")9899def fig_discontinuite(res):100    r = res["rdd"]; x, y = r["x"], r["y"]101    fig, ax = plt.subplots(figsize=(LARGEUR, 3.0))102    ax.scatter(x, y, s=6, color=palette.grid, lw=0, label="Immeubles simulés (n = 400)")103    bords = np.arange(-2, 2.01, 0.2); centres = (bords[:-1] + bords[1:]) / 2104    moy = [y[(x >= a) & (x < b)].mean() for a, b in zip(bords[:-1], bords[1:])]105    ax.plot(centres, moy, "o", color=palette.cat[0], ms=4.5, label="Moyenne par tranche de 200 m")106    e = r["estimations"][1.0]; b = e["beta"]107    xg = np.linspace(-1, 0, 50); ax.plot(xg, b[0] + b[2] * xg, color=palette.cat[1], lw=1.8, label="Ajustement local linéaire (h = 1 km)")108    xd = np.linspace(0, 1, 50); ax.plot(xd, b[0] + b[1] + (b[2] + b[3]) * xd, color=palette.cat[1], lw=1.8)109    ax.axvline(0, color=palette.ink2, lw=0.9, ls="--"); ax.text(0.03, y.min() + 0.01, "frontière de zone", fontsize=7.5, color=palette.ink2)110    ax.annotate("", xy=(0, b[0] + b[1]), xytext=(0, b[0]), arrowprops=dict(arrowstyle="<->", color=palette.ink, lw=1.0))111    ax.text(0.05, b[0] + b[1] / 2, f"saut estimé τ̂ = {e['tau']:.3f}\n(vrai 0,050)", fontsize=7.5, color=palette.ink, va="center")112    ax.set_xlabel("Distance signée à la frontière (km ; positif = zone de densité supérieure)"); ax.set_ylabel("ln de la valeur simulée")113    ax.legend(fontsize=7, loc="upper left"); ax.grid(True, lw=0.4)114    savefig(fig, "ch36/ch36-fig-discontinuite.pdf")115116def fig_sensibilite(res):117    es0 = res["es"]; esv = res["sens"]["event_study"]118    fig, ax = plt.subplots(figsize=(LARGEUR, 3.0))119    ax.fill_between(esv.index, esv["ic_bas"], esv["ic_haut"], color=palette.cat[1], alpha=0.15, lw=0)120    ax.plot(es0.index, es0["coef"], "o-", color=palette.cat[0], ms=3.6, lw=1.3, label="Tendances parallèles respectées (générateur)")121    ax.plot(esv.index, esv["coef"], "s-", color=palette.cat[1], ms=3.6, lw=1.3, label="Traités avec +1 point de croissance par an (IC à 95 %)")122    ax.step(es0.index, es0["vrai"], where="mid", color=palette.ink, lw=1.1, ls="--", label="Effet déclaré δ(k)")123    pre = esv.loc[-7:-2]; p = np.polyfit(pre.index.values.astype(float), pre["coef"].values, 1)124    xx = np.linspace(-8, -1, 20); ax.plot(xx, np.polyval(p, xx), color=palette.cat[1], lw=0.9, ls=":")125    ax.text(-7.9, np.polyval(p, -7.9) + 0.012, f"pente pré-traitement\n{p[0]:.4f} par trimestre", fontsize=7, color=palette.ink)126    ax.axhline(0, color=palette.ink2, lw=0.8); ax.axvline(-0.5, color=palette.ink2, lw=0.8, ls=":")127    ax.set_xlabel("Trimestres depuis la mise en service, k"); ax.set_ylabel("Effet sur ln RNE (points de log)"); ax.legend(fontsize=7, loc="upper left"); ax.grid(True, lw=0.4)128    savefig(fig, "ch36/ch36-fig-sensibilite.pdf")129130def main():131    setup(); res = c36.main(verbose=False)132    for f in (fig_trajectoires, fig_event_study, fig_bacon, fig_synthetique, fig_discontinuite, fig_sensibilite):133        f(res); print("ok", f.__name__)134135if __name__ == "__main__":136    main()137