spb/traite-evaluation-qc

1"""Données de panel et identification causale (chapitre 36) : transformation within, erreurs-types en grappes (CR1),2effets fixes bidirectionnels par indicatrices (via hedonique.mco), doubles différences 2 × 2, étude d'événement,3estimateur de Callaway–Sant'Anna (2021) par cohortes avec contrôles « jamais traités », décomposition de4Goodman-Bacon (2021), contrôle synthétique (Abadie–Diamond–Hainmueller 2010) et discontinuité locale linéaire.56Toutes les fonctions sont pures (numpy/pandas/scipy) et testées dans analysis/python/tests/test_ch36.py.7Conventions : un panel « long » est un DataFrame avec une colonne d'unité, une colonne de période entière (0, 1, …, T−1),8une variable dépendante et, pour les plans échelonnés, une colonne `g` = première période traitée (NaN si jamais traité)."""9from __future__ import annotations10import numpy as np, pandas as pd11from scipy import optimize12from . import hedonique as h1314# ----------------------------------------------------------------------------- MCO en grappes et within1516def mco_grappes(X: np.ndarray, y: np.ndarray, grappes: np.ndarray) -> dict:17    """MCO avec erreurs-types classiques, HC1 et robustes aux grappes (CR1 : facteur G/(G−1)·(N−1)/(N−K)).18    X contient la constante (ou les indicatrices qui en tiennent lieu)."""19    X = np.asarray(X, float); y = np.asarray(y, float); grappes = np.asarray(grappes)20    base = h.mco(X, y)21    n, k = X.shape; e = base["residus"]22    bread = np.linalg.pinv(X.T @ X)23    ids = np.unique(grappes); G = len(ids)24    meat = np.zeros((k, k))25    for gid in ids:26        m = grappes == gid27        s = X[m].T @ e[m]28        meat += np.outer(s, s)29    c = G / (G - 1) * (n - 1) / (n - k)30    V = c * bread @ meat @ bread31    out = dict(base); out.update({"se_grappes": np.sqrt(np.diag(V)), "V_grappes": V, "G": G, "correction_cr1": c})32    return out3334def within(y: np.ndarray, X: np.ndarray, unites: np.ndarray) -> tuple[np.ndarray, np.ndarray]:35    """Transformation within (écarts aux moyennes d'unité) de y et de chaque colonne de X."""36    df = pd.DataFrame(np.column_stack([y, X])); g = df.groupby(np.asarray(unites))37    Z = (df - g.transform("mean")).to_numpy()38    return Z[:, 0], Z[:, 1:]3940def within_deux_voies(y: np.ndarray, X: np.ndarray, unites: np.ndarray, temps: np.ndarray) -> tuple[np.ndarray, np.ndarray]:41    """Transformation within bidirectionnelle (panel cylindré) : z − z̄_i − z̄_t + z̄."""42    df = pd.DataFrame(np.column_stack([y, X])); u = np.asarray(unites); t = np.asarray(temps)43    Z = (df - df.groupby(u).transform("mean") - df.groupby(t).transform("mean") + df.mean()).to_numpy()44    return Z[:, 0], Z[:, 1:]4546def estimateur_within(y, X, unites, grappes=None) -> dict:47    """Estimateur within (effets fixes d'unité) : MCO sur les données transformées, ddl corrigés du nombre d'unités,48    erreurs-types en grappes par unité (par défaut)."""49    yw, Xw = within(y, X, unites)50    grappes = np.asarray(unites) if grappes is None else np.asarray(grappes)51    r = mco_grappes(Xw, yw, grappes)52    n, k = Xw.shape; N = len(np.unique(unites))53    # correction des degrés de liberté pour les erreurs-types classiques (N constantes éliminées)54    s2 = r["residus"] @ r["residus"] / (n - k - N)55    r["se"] = np.sqrt(np.diag(s2 * np.linalg.pinv(Xw.T @ Xw))); r["sigma2"] = s256    return r5758def indicatrices(codes, omettre_premiere: bool = True, prefixe: str = "") -> pd.DataFrame:59    """Indicatrices (0/1) d'un vecteur de codes, dans l'ordre trié ; la première modalité est omise par défaut."""60    s = pd.Series(np.asarray(codes))61    D = pd.get_dummies(s, prefix=prefixe, prefix_sep="" if prefixe == "" else ":", dtype=float)62    if omettre_premiere: D = D.iloc[:, 1:]63    return D.reset_index(drop=True)6465def twfe(df: pd.DataFrame, y: str, regresseurs: list[str], unite: str, temps: str, fe_supp: str | None = None,66         grappes: str | None = None) -> dict:67    """Effets fixes bidirectionnels par indicatrices : y = Σ β_r x_r + α_i + λ_t (+ θ_{s,t} si fe_supp) + ε.68    fe_supp : colonne de strate (ex. RMR) → indicatrices strate × période (période 0 omise dans chaque strate), qui69    remplacent les effets de période communs. Erreurs-types en grappes par `grappes` (défaut : l'unité)."""70    d = df.reset_index(drop=True)71    parts = [pd.DataFrame({"const": np.ones(len(d))}), d[regresseurs].astype(float).reset_index(drop=True),72             indicatrices(d[unite], prefixe="u")]73    if fe_supp is None:74        parts.append(indicatrices(d[temps], prefixe="t"))75    else:76        st = d[fe_supp].astype(str) + "@" + d[temps].astype(int).astype(str)77        D = indicatrices(st, omettre_premiere=False, prefixe="st")78        t0 = d[temps].min()79        D = D[[c for c in D.columns if not c.endswith(f"@{int(t0)}")]]80        parts.append(D.reset_index(drop=True))81    X = pd.concat(parts, axis=1)82    cl = d[unite].values if grappes is None else d[grappes].values83    r = mco_grappes(X.values, d[y].values, cl)84    idx = {nom: i for i, nom in enumerate(X.columns)}85    coef = pd.DataFrame({"coef": [r["beta"][idx[v]] for v in regresseurs], "se": [r["se"][idx[v]] for v in regresseurs],86                         "se_hc1": [r["se_hc1"][idx[v]] for v in regresseurs], "se_grappes": [r["se_grappes"][idx[v]] for v in regresseurs]},87                        index=regresseurs)88    return {"coef": coef, "n": r["n"], "k": r["k"], "G": r["G"], "sigma2": r["sigma2"], "r2": r["r2"], "residus": r["residus"], "colonnes": list(X.columns), "beta": r["beta"]}8990# ----------------------------------------------------------------------------- doubles différences 2 × 29192def did_2x2(df: pd.DataFrame, y: str, unite: str, temps: str, traite: str, avant, apres) -> dict:93    """DiD à deux périodes : quatre moyennes de cellule (fenêtres `avant` et `apres` = listes de périodes), estimateur94    (ȳ_tr,ap − ȳ_tr,av) − (ȳ_té,ap − ȳ_té,av), et version agrégée par unité (différence longue Δy_i régressée sur95    l'indicatrice de traitement, erreurs-types HC1 : robustes à l'autocorrélation intra-unité)."""96    d = df[df[temps].isin(list(avant) + list(apres))].copy()97    d["periode"] = np.where(d[temps].isin(list(apres)), 1, 0)98    cell = d.groupby([traite, "periode"])[y].mean().unstack()99    est = (cell.loc[1, 1] - cell.loc[1, 0]) - (cell.loc[0, 1] - cell.loc[0, 0])100    # agrégation par unité : Δy_i = ȳ_i,après − ȳ_i,avant101    u = d.groupby([unite, traite, "periode"])[y].mean().unstack("periode")102    delta = (u[1] - u[0]).reset_index(); delta.columns = [unite, traite, "dy"]103    X = np.column_stack([np.ones(len(delta)), delta[traite].values.astype(float)])104    r = h.mco(X, delta["dy"].values)105    return {"cellules": cell, "estimateur": float(est), "coef_reg": float(r["beta"][1]), "se_hc1": float(r["se_hc1"][1]),106            "se_classique": float(r["se"][1]), "n_traites": int(delta[traite].sum()), "n_temoins": int((delta[traite] == 0).sum()),107            "deltas": delta}108109# ----------------------------------------------------------------------------- étude d'événement110111def event_study(df: pd.DataFrame, y: str, unite: str, temps: str, g: str, fenetre: tuple[int, int] = (-8, 8), ref: int = -1,112                fe_supp: str | None = None, cohortes=None) -> pd.DataFrame:113    """Étude d'événement par effets fixes bidirectionnels : indicatrices du temps relatif k = t − g (bornes de `fenetre`114    regroupées : « ≤ kmin » et « ≥ kmax »), période de référence `ref` omise, unités jamais traitées (g NaN) à zéro.115    `cohortes` : restreint les unités traitées à ces valeurs de g (les autres cohortes sont exclues du panel)."""116    d = df.copy()117    if cohortes is not None:118        d = d[d[g].isna() | d[g].isin(list(cohortes))].copy()119    kmin, kmax = fenetre120    k = d[temps] - d[g]121    k = k.where(d[g].notna())122    kb = k.clip(lower=kmin, upper=kmax)123    noms = []124    for j in range(kmin, kmax + 1):125        if j == ref: continue126        nom = f"k{j:+d}"; d[nom] = ((kb == j) & d[g].notna()).astype(float); noms.append(nom)127    r = twfe(d, y, noms, unite, temps, fe_supp=fe_supp)128    out = r["coef"].copy(); out["k"] = [int(n[1:]) for n in out.index]129    out = pd.concat([out, pd.DataFrame({"coef": [0.0], "se": [0.0], "se_hc1": [0.0], "se_grappes": [0.0], "k": [ref]}, index=[f"k{ref:+d}"])]).sort_values("k")130    out["ic_bas"] = out["coef"] - 1.96 * out["se_grappes"]; out["ic_haut"] = out["coef"] + 1.96 * out["se_grappes"]131    return out.set_index("k")132133# ----------------------------------------------------------------------------- Callaway–Sant'Anna (simplifié)134135def callaway_santanna(df: pd.DataFrame, y: str, unite: str, temps: str, g: str, strate: str | None = None,136                      controle: str = "jamais", B: int = 999, graine: int = 20261031) -> dict:137    """ATT(g, t) = E[y_t − y_{g−1} | G = g] − E[y_t − y_{g−1} | C], période de base universelle g − 1, contrôles138    « jamais traités » (controle='jamais') ou « pas encore traités à t » (controle='pas_encore') ; si `strate` est donnée,139    les contrôles sont restreints à la strate de la cohorte (ex. même RMR). Agrégations : dynamique (par k), par cohorte,140    globale (pondérée par la taille des cohortes). Erreurs-types : bootstrap par blocs d'unités (B tirages)."""141    W = df.pivot(index=unite, columns=temps, values=y)142    G = df.groupby(unite)[g].first().reindex(W.index)143    S = df.groupby(unite)[strate].first().reindex(W.index) if strate else None144    periodes = list(W.columns); cohortes = sorted(G.dropna().unique())145146    def _att_table(W, G, S):147        rows = []148        for gg in cohortes:149            tr = G.index[G == gg]150            base = int(gg) - 1151            for t in periodes:152                if controle == "jamais":153                    ctl = G.index[G.isna()]154                else:155                    ctl = G.index[G.isna() | (G > max(t, gg))]156                if S is not None:157                    s_g = S.loc[tr].mode().iloc[0]; ctl = [c for c in ctl if S.loc[c] == s_g]158                dy_tr = W.loc[tr, t] - W.loc[tr, base]; dy_c = W.loc[ctl, t] - W.loc[ctl, base]159                att = dy_tr.mean() - dy_c.mean()160                se = np.sqrt(dy_tr.var(ddof=1) / len(dy_tr) + dy_c.var(ddof=1) / len(dy_c)) if len(dy_tr) > 1 and len(dy_c) > 1 else np.nan161                rows.append({"g": int(gg), "t": int(t), "k": int(t) - int(gg), "att": float(att), "se_analytique": float(se), "n_g": len(tr), "n_c": len(ctl)})162        return pd.DataFrame(rows)163164    def _agreger(tab):165        post = tab[tab["k"] >= 0]166        par_cohorte = post.groupby("g").apply(lambda d: d["att"].mean(), include_groups=False)167        taille = post.groupby("g")["n_g"].first()168        globale = float((par_cohorte * taille).sum() / taille.sum())169        dyn = tab.groupby("k").apply(lambda d: float((d["att"] * d["n_g"]).sum() / d["n_g"].sum()), include_groups=False)170        return par_cohorte, globale, dyn, taille171172    tab = _att_table(W, G, S)173    par_cohorte, globale, dyn, taille = _agreger(tab)174    rng = np.random.default_rng(graine)175    boot_c, boot_g, boot_d = [], [], []176    unites = np.array(W.index)177    for _ in range(B):178        # bootstrap stratifié : chaque cohorte et chaque strate de contrôles (ex. RMR) est rééchantillonnée séparément,179        # de sorte que le plan (tailles des groupes, présence de témoins dans chaque strate) est conservé180        idx = []181        for gg in cohortes:182            u = unites[(G == gg).values]; idx += list(rng.choice(u, size=len(u), replace=True))183        strates_ctl = S.loc[G.isna()].unique() if S is not None else [None]184        for st in strates_ctl:185            m = G.isna().values if st is None else (G.isna() & (S == st)).values186            u = unites[m]; idx += list(rng.choice(u, size=len(u), replace=True))187        Wb = W.loc[idx].reset_index(drop=True); Gb = G.loc[idx].reset_index(drop=True); Sb = S.loc[idx].reset_index(drop=True) if S is not None else None188        Wb.index = range(len(Wb)); Gb.index = range(len(Gb))189        if Sb is not None: Sb.index = range(len(Sb))190        tb = _att_table(Wb, Gb, Sb); pc, gl, dn, _ = _agreger(tb)191        boot_c.append(pc.reindex(par_cohorte.index).values); boot_g.append(gl); boot_d.append(dn.reindex(dyn.index).values)192    boot_c = np.array(boot_c); boot_d = np.array(boot_d)193    cohortes_df = pd.DataFrame({"att": par_cohorte, "se_boot": boot_c.std(axis=0, ddof=1), "n_g": taille})194    dyn_df = pd.DataFrame({"att": dyn, "se_boot": np.nanstd(boot_d, axis=0, ddof=1)})195    dyn_df["ic_bas"] = dyn_df["att"] - 1.96 * dyn_df["se_boot"]; dyn_df["ic_haut"] = dyn_df["att"] + 1.96 * dyn_df["se_boot"]196    return {"att_gt": tab, "par_cohorte": cohortes_df, "globale": globale, "globale_se_boot": float(np.std(boot_g, ddof=1)), "dynamique": dyn_df}197198# ----------------------------------------------------------------------------- décomposition de Goodman-Bacon199200def bacon(df: pd.DataFrame, y: str, unite: str, temps: str, g: str) -> dict:201    """Décomposition de Goodman-Bacon (2021) de l'estimateur TWFE statique (panel cylindré, sans covariable) en202    comparaisons 2 × 2 : chaque cohorte vs jamais traités ; cohorte précoce (traitée) vs cohorte tardive (témoin, avant203    son traitement) ; cohorte tardive (traitée) vs cohorte précoce (témoin déjà traité). Vérifie Σ s = 1 et204    Σ s·β_2x2 = β_TWFE."""205    W = df.pivot(index=unite, columns=temps, values=y)206    G = df.groupby(unite)[g].first().reindex(W.index)207    T = W.shape[1]; periodes = np.array(W.columns); N = W.shape[0]208    cohortes = sorted(G.dropna().unique())209    n = {gg: (G == gg).sum() / N for gg in cohortes}; n_U = G.isna().sum() / N210    Dbar = {gg: (periodes >= gg).mean() for gg in cohortes}211    # variance de D doublement centré212    d = df.copy(); d["D"] = ((d[temps] >= d[g]) & d[g].notna()).astype(float)213    Dt = d["D"] - d.groupby(unite)["D"].transform("mean") - d.groupby(temps)["D"].transform("mean") + d["D"].mean()214    VD = float((Dt ** 2).mean())215216    def m(ids, mask_t):217        return W.loc[ids, periodes[mask_t]].values.mean()218    rows = []219    U = G.index[G.isna()]220    for gg in cohortes:221        tr = G.index[G == gg]222        if len(U) > 0:223            b = (m(tr, periodes >= gg) - m(tr, periodes < gg)) - (m(U, periodes >= gg) - m(U, periodes < gg))224            nkU = n[gg] / (n[gg] + n_U)225            s = (n[gg] + n_U) ** 2 * nkU * (1 - nkU) * Dbar[gg] * (1 - Dbar[gg]) / VD226            rows.append({"traites": int(gg), "temoins": "jamais", "type": "traités vs jamais traités", "beta": float(b), "poids": float(s)})227    for i, gk in enumerate(cohortes):228        for gl in cohortes[i + 1:]:229            k_ids = G.index[G == gk]; l_ids = G.index[G == gl]230            nkl = n[gk] / (n[gk] + n[gl])231            # précoce (k) traité, tardif (l) témoin : fenêtre t < gl232            b_kl = (m(k_ids, (periodes >= gk) & (periodes < gl)) - m(k_ids, periodes < gk)) - (m(l_ids, (periodes >= gk) & (periodes < gl)) - m(l_ids, periodes < gk))233            s_kl = ((n[gk] + n[gl]) * (1 - Dbar[gl])) ** 2 * nkl * (1 - nkl) * ((Dbar[gk] - Dbar[gl]) / (1 - Dbar[gl])) * ((1 - Dbar[gk]) / (1 - Dbar[gl])) / VD234            rows.append({"traites": int(gk), "temoins": int(gl), "type": "précoce vs tardive (pas encore traitée)", "beta": float(b_kl), "poids": float(s_kl)})235            # tardif (l) traité, précoce (k) témoin déjà traité : fenêtre t ≥ gk236            b_lk = (m(l_ids, periodes >= gl) - m(l_ids, (periodes >= gk) & (periodes < gl))) - (m(k_ids, periodes >= gl) - m(k_ids, (periodes >= gk) & (periodes < gl)))237            s_lk = ((n[gk] + n[gl]) * Dbar[gk]) ** 2 * nkl * (1 - nkl) * (Dbar[gl] / Dbar[gk]) * ((Dbar[gk] - Dbar[gl]) / Dbar[gk]) / VD238            rows.append({"traites": int(gl), "temoins": int(gk), "type": "tardive vs précoce (déjà traitée)", "beta": float(b_lk), "poids": float(s_lk)})239    tab = pd.DataFrame(rows)240    tw = twfe(d, y, ["D"], unite, temps)241    return {"table": tab, "somme_poids": float(tab["poids"].sum()), "beta_pondere": float((tab["beta"] * tab["poids"]).sum()),242            "beta_twfe": float(tw["coef"].loc["D", "coef"]), "var_D": VD}243244# ----------------------------------------------------------------------------- contrôle synthétique245246def controle_synthetique(Y0: np.ndarray, y1: np.ndarray, T0: int) -> dict:247    """Poids w ≥ 0, Σ w = 1, minimisant Σ_{t<T0} (y1_t − Y0_t w)² (Abadie, Diamond et Hainmueller 2010, sans prédicteurs248    auxiliaires). Y0 : (T × J) donneurs ; y1 : (T,) unité traitée. Retourne poids, trajectoire synthétique, écarts,249    RMSPE pré et post."""250    Y0 = np.asarray(Y0, float); y1 = np.asarray(y1, float); J = Y0.shape[1]251    A = Y0[:T0]; b = y1[:T0]252    def f(w): r = A @ w - b; return float(r @ r)253    def grad(w): return 2 * A.T @ (A @ w - b)254    cons = ({"type": "eq", "fun": lambda w: w.sum() - 1.0, "jac": lambda w: np.ones(J)},)255    w0 = np.ones(J) / J256    res = optimize.minimize(f, w0, jac=grad, bounds=[(0, 1)] * J, constraints=cons, method="SLSQP", options={"maxiter": 2000, "ftol": 1e-14})257    w = np.clip(res.x, 0, None); w = w / w.sum()258    synth = Y0 @ w; gap = y1 - synth259    return {"poids": w, "synthetique": synth, "ecarts": gap, "rmspe_pre": float(np.sqrt(np.mean(gap[:T0] ** 2))),260            "rmspe_post": float(np.sqrt(np.mean(gap[T0:] ** 2))), "effet_moyen_post": float(gap[T0:].mean()), "convergence": bool(res.success)}261262def placebos_synthetiques(Y0: np.ndarray, T0: int) -> pd.DataFrame:263    """Placebos « dans l'espace » : chaque donneur devient à son tour l'unité traitée, les autres forment le bassin."""264    Y0 = np.asarray(Y0, float); J = Y0.shape[1]; rows = []265    for j in range(J):266        r = controle_synthetique(np.delete(Y0, j, axis=1), Y0[:, j], T0)267        rows.append({"donneur": j, "rmspe_pre": r["rmspe_pre"], "rmspe_post": r["rmspe_post"], "ratio": r["rmspe_post"] / r["rmspe_pre"], "effet_moyen_post": r["effet_moyen_post"]})268    return pd.DataFrame(rows)269270# ----------------------------------------------------------------------------- discontinuité (locale linéaire)271272def rdd_local_lineaire(x: np.ndarray, y: np.ndarray, h_bande: float, c: float = 0.0) -> dict:273    """Saut à la frontière c par régression locale linéaire à noyau triangulaire de bande h : y = a + τ·1[x ≥ c] + b(x−c)274    + d·1[x ≥ c](x−c), pondérée par K((x−c)/h). Erreur-type HC1 de τ."""275    x = np.asarray(x, float); y = np.asarray(y, float); u = (x - c) / h_bande276    w = np.clip(1 - np.abs(u), 0, None); m = w > 0277    D = (x[m] >= c).astype(float); xc = x[m] - c278    X = np.column_stack([np.ones(m.sum()), D, xc, D * xc]); sw = np.sqrt(w[m])279    r = h.mco(X * sw[:, None], y[m] * sw)280    return {"tau": float(r["beta"][1]), "se_hc1": float(r["se_hc1"][1]), "n_eff": int(m.sum()), "beta": r["beta"], "bande": h_bande}281