🙎‍♂️

合成コントロール法は非定常データに弱い？

に公開
統計学
因果推論
効果検証
合成コントロール法
tech

本記事では以下の論文を基に議論していきます．
https://doi.org/10.1080/07350015.2020.1799814

実務で因果推論を行う際に，データの定常性を確認したことはありますか？

観測できない事象に対して因果を明らかにしようとする試みは，事業に携わる人間にとって魅力的であり，誰もが関心を抱くテーマである．
しかし，その真の効果を直接確認することは不可能であるため，因果推論の妥当性を評価することは極めて困難である．

もちろん，ドメイン知識や仮説に基づいて一定の判断を下すことは可能であるが，その判断には再現性が乏しいという問題が残る．

近年，生成AIの発展により，因果推論を伴う効果検証のプロセスは次第に容易になりつつある．

こうした状況の中で，我々データアナリストやデータサイエンティストが価値を創出するうえで重要な要素の一つは，分析結果の妥当性と再現性をいかに担保するかであると言えるだろう．
少なくとも私はそう思う．
(と言いつつ筆者はしがない大学院生なのですが)

本記事の要点

・合成コントロール法を含む多くの反実仮想推定法は，時系列の非定常性を無視している
・特にデータが単位根過程であり共和分関係にない場合，推定される反実仮想は意味をなさない
・因果推論を行う前にデータの定常性と共和分関係を確認するべき

単位根過程であり共和分関係にない場合どのような推定をするか

データ分析で因果推論に取り組む際，合成コントロール法（以下 SCM）は非常に強力なツールです．

しかし，Masini & Medeiros (2022) はこの推定が非定常なデータ（単位根を含む共和分関係にない時系列）に対して極めて危険になりうることを明らかにしています．

前提知識として見せかけの回帰，単位根過程と共和分関係について理解しておくとよいです．以下の記事がより詳細です．
https://tjo.hatenablog.com/entry/2013/04/23/190417
まず，各系列が単位根を持つ（非定常な）ランダムウォーク型の過程だとします．

Y_{it} = Y_{i,t-1} + u_{it}, \quad u_{it} \sim \text{i.i.d.}(0, \sigma_i^2)

ここでi = 1は処置群，i = 2,..,J +1はコントロール群です．
この時，各系列は時間とともに累積ノイズとなります．

Y_{it} = Y_{i0} + \sum_{s=1}^{t} u_{is}

SCMは、介入前のデータから最適な重み
w^{*} = (w_2, \ldots, w_{J+1}) を選び，処置群の系列を次のような線形結合で近似します．

\widehat{Y}_{1t} = \sum_{j=2}^{J+1} w_j^{*} Y_{jt}

ここで「良い重み」とは、介入前の期間 t \leq T_0 における誤差

e_t = Y_{1t} - \widehat{Y}_{1t}

を最小化するように選ばれます．

ここで，共和分関係がないとは各 Y_{it} が独立したランダムウォークであり，どの線形結合をとっても定常な系列が得られないことです．

つまり，任意の重みベクトル w に対して，以下のZ_tが定常でない場合です．

Z_t = Y_{1t} - w'Y_t

この時，推定誤差e_t は次のように展開されます．

\begin{aligned} e_t &= (Y_{1,t-1} + u_{1t}) - \sum_j w_j^{*}(Y_{j,t-1} + u_{jt}) \\[-4pt] &= \underbrace{(Y_{1,t-1} - \sum_j w_j^{*} Y_{j,t-1})}_{=\, e_{t-1}} + \left(u_{1t} - \sum_j w_j^{*} u_{jt}\right) \\[-4pt] &= e_{t-1} + \left(u_{1t} - \sum_j w_j^{*} u_{jt}\right) \end{aligned}

したがって，e_t自体も単位根過程になります．

誤差は時間とともにランダムウォーク的に拡散し，介入がなくても効果があるように見えるという見せかけの効果が生じます．

これは，処置群とコントロール群の系列が独立した確率トレンドを持つため，偶然のトレンドの一致が介入効果に誤認されることを意味します．

つまり，何でもかんでも時系列データに対して合成コントロール法を使うと痛い目を見るということですね．

シミュレーションで確認する

# === 共和分ドナー 0 vs 30 比較（tau=0, 一括実行セル） ===
import numpy as np
import matplotlib.pyplot as plt
from dataclasses import dataclass
from typing import Optional, Tuple, Dict

# --- 基本関数群（SC推定 + 射影勾配法） ---

def project_to_simplex(v, z=1.0):
    u = np.sort(v)[::-1]
    cssv = np.cumsum(u) - z
    rho = np.nonzero(u - cssv / (np.arange(v.size) + 1) > 0)[0]
    theta = cssv[rho[-1]] / (rho[-1] + 1.0) if rho.size else 0.0
    return np.maximum(v - theta, 0.0)

def pgd_quadratic(X, y, lam=0.0, max_iter=5000, tol=1e-8, step=None):
    n = X.shape[1]
    w = np.full(n, 1/n)
    if step is None:
        v = np.random.default_rng(0).normal(size=n)
        v /= np.linalg.norm(v)
        for _ in range(30):
            v = X.T @ (X @ v)
            v /= np.linalg.norm(v)
        L = v @ (X.T @ (X @ v))
        step = 1.0 / (L + lam + 1e-8)
    for _ in range(max_iter):
        grad = X.T @ (X @ w - y) + lam * w
        w_new = project_to_simplex(w - step * grad)
        if np.linalg.norm(w_new - w) < tol:
            break
        w = w_new
    return w

def fit_sc(y1, Y0, T0, lam=0.0):
    X, y = Y0[:T0], y1[:T0]
    try:
        w_ls = np.linalg.lstsq(X, y, rcond=None)[0]
    except np.linalg.LinAlgError:
        w_ls = np.zeros(X.shape[1])
    w0 = np.clip(w_ls, 0, None)
    if w0.sum() == 0:
        w0[:] = 1.0 / len(w0)
    w0 /= w0.sum()
    w = pgd_quadratic(X, y, lam=lam)
    return Y0 @ w, w

def sc_resample_test(y1, Y0, T0, w, B=500, seed=42):
    rng = np.random.default_rng(seed)
    T, T2 = y1.shape[0], y1.shape[0] - T0
    yhat = Y0 @ w
    obs = (y1[T0:] - yhat[T0:]).mean()
    starts = rng.integers(0, T0 - T2 + 1, size=B)
    pseudo = [(y1[s:s+T2] - (Y0[s:s+T2] @ w)).mean() for s in starts]
    p = (np.abs(pseudo) >= abs(obs)).mean()
    return {"ate": obs, "p_value": p}

@dataclass
class MixConfig:
    T: int = 80
    T0: int = 40
    n_donors: int = 30
    k_cointeg: int = 10
    tau: float = 0.0
    seed: Optional[int] = 10
    sigma_eps: float = 1.0
    sigma_u: float = 1.0

def simulate_mixed_panel(cfg: MixConfig):
    rng = np.random.default_rng(cfg.seed)
    s = rng.normal(0, cfg.sigma_eps, size=cfg.T).cumsum()
    y1 = s + rng.normal(0, cfg.sigma_u, size=cfg.T)
    loads = rng.normal(0.8, 0.15, size=cfg.k_cointeg)
    Yc = s[:, None] * loads[None, :] + rng.normal(0, cfg.sigma_u, size=(cfg.T, cfg.k_cointeg))
    Ys = rng.normal(0, cfg.sigma_eps, size=(cfg.T, cfg.n_donors - cfg.k_cointeg)).cumsum(axis=0)
    Y0 = np.concatenate([Yc, Ys], axis=1) if cfg.k_cointeg < cfg.n_donors else Yc
    y1 = y1.copy(); y1[cfg.T0:] += cfg.tau
    return y1, Y0

# --- 比較実験：k=0 vs k=30 ---
configs = [
    MixConfig(k_cointeg=0, tau=0.0, seed=10),
    MixConfig(k_cointeg=30, tau=0.0, seed=10)
]

plt.figure(figsize=(8,5))
colors = ["red", "blue"]
labels = ["All spurious (k=0)", "All cointegrated (k=30)"]

for cfg, c, lab in zip(configs, colors, labels):
    y1, Y0 = simulate_mixed_panel(cfg)
    yhat, w = fit_sc(y1, Y0, cfg.T0)
    res = sc_resample_test(y1, Y0, cfg.T0, w, B=500)
    plt.plot(yhat, color=c, label=f"{lab}, p={res['p_value']:.3f}")
    plt.axvline(cfg.T0, linestyle="--", color="gray", alpha=0.6)

plt.plot(y1, color="black", alpha=0.3, label="Treated (true, tau=0)")
plt.title("SC fit comparison: All spurious vs All cointegrated (tau=0)")
plt.xlabel("Time")
plt.ylabel("Value")
plt.legend()
plt.tight_layout()
plt.show()

シミュレーションの設定は以下です．
時点数：T = 80
介入前期間：T_0 = 40
ドナーユニット：n_{donors} = 30
真の介入効果：tau = 0

以上の結果より，介入効果が0であっても単位根過程で共和分関係がない場合，あたかも介入効果があるように見誤る可能性があることがわかります．

まとめ

本記事では，合成コントロールを用いるときに暗黙的に定常性を仮定しており，ある状況下で著しく介入効果の信頼性が下がることを解説しました．
因果推論って難しいですね．

Discussion