Causal Impactとは？

(画像は、Geminiで作成されたものです)

Causal Impactの概要

Causal Impactは、2015年にGoogleのKay H. Brodersen氏らによって提案された統計的アプローチです。元々は “Inferring causal effects from observational time series data using Bayesian structural time-series models” という正式論文で発表された理論ですが、その後、同名のRパッケージ（ツール）の公開に合わせて、内容をより実践向けに落とし込んだ論文「Inferring causal impact using Bayesian structural time-series models」が公開されました。この手法の最大の特徴は、「反事実（Counterfactual）」という概念を時系列データに持ち込んだ点にあります。論文の中では、その核心について次のように述べられています。 “The basic idea is to generalize the widely used difference-in-differences approach to a state-space framework that allows us to infer the counterfactual market response that would have occurred in the absence of the intervention.” （意訳）基本的なアイデアは、広く使われている「差の差の分析（DID）」を状態空間モデルへと拡張し、もし施策（介入）を行っていなかったら発生していたであろう「反事実」の市場の反応を推論することである。
従来の「施策の前後で平均値を比べる」ような単純な方法では、季節による変動や、市場全体の自然な成長といった「外部要因」を施策の効果だと誤認してしまうリスク（選択バイアス）がありました。 CausalImpactは、施策の影響を受けていない他のデータ（例：競合の動きや他地域のデータなど）を組み込んだ「ベイズ構造時系列モデル（BSTS）」を構築します。これにより、「もし施策を行っていなかったら、今頃どうなっていたか」というパラレルワールドのタイムライン（反事実）を高精度に予測します。そして、「実際の観測値」と「反事実の予測値」のギャップを測ることで、外部要因のノイズを綺麗に取り除いた「純粋な施策の効果」を浮き彫りにするのです。

CausalImpactの基本ロジック：「反事実」の予測

CausalImpactがやっていることを一言で言えば、「もしこの施策をやっていなかったら、データはどうなっていたか？（反事実：Counterfactual）」というパラレルワールドの時系列データを予測し、実際の観測値と比較することです。

具体的な効果測定のメカニズムと、それを支える数理モデルの構造を紐解いていきましょう。施策効果を導く「引き算」のシンプルさCausalImpactにおいて、施策による因果効果（Causal Effect）は以下の非常にシンプルな引き算で定義されます。

\tau_t = y_t - \hat{y}_t

$\tau_t$ ：時点 $t$ における施策効果（Causal Effect）
$y_t$ ：時点 $t$ での実際の観測値（施策を行った現実のデータ）
$\hat{y}_t$ ：時点 $t$ での予測値（反事実）（もし施策を行っていなかった場合のデータ）

現実のデータ（ $y_t$ ）は手元にありますから、私たちが解くべき問題は「どうやって施策がなかった場合の予測値（ $\hat{y}_t$ ）を高精度に叩き出すか」の1点に絞られます。

反事実を予測する「ベイズ構造時系列モデル（BSTS）」の数理

Googleの論文で提案されたCausalImpactのコアは、「ベイズ構造時系列モデル（BSTS: Bayesian Structural Time Series）」という柔軟な時系列モデルです。本記事のスクラッチ実装でもベースとなる、最も標準的なモデル（ローカルレベルモデル＋共変量）の数式は以下の2つのパーツ（方程式）で構成されています。

① 観測方程式（Observation Equation）
手元で見えているデータが、どのような要素で構成されているかを表す式です。

y_t = \mu_t + \beta X_t + \epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0, \sigma_y^2)

$y_t$ ：ターゲット変数（自社の売上やコンバージョン数など）
$\mu_t$ ：その時々のベースとなるトレンド（潜在状態）
$X_t$ ：共変量（コントロール変数）（施策の影響を受けていない、競合のデータや検索トレンドなど）
$\beta$ ：共変量がターゲットに与える影響の強さ（係数）
$\epsilon_t$ ：その日の偶然のばらつきを表す観測ノイズ（平均0、分散 $\sigma_y^2$ の正規分布に従う）

② 状態方程式（System Equation）
目に見えない「ベースのトレンド（ $\mu_t$ ）」が、時間の経過とともにどう変化していくかを表す式です。

\mu_t = \mu_{t-1} + \eta_t, \quad \eta_t \sim \mathcal{N}(0, \sigma_\mu^2)

$\mu_t$ は、1つ前の時点のトレンド $\mu_{t-1}$ に、その時々のランダムな変動 $\eta_t$ （システムノイズ）が加わることで日々緩やかに変化していきます（これをランダムウォークと呼びます）。

この数式がなぜ「施策効果」を炙り出せるのか？
CausalImpactの予測マジックは、データを「施策前の期間（Pre-period）」と「施策後の期間（Post-period）」に分けて処理することで成立します。

【施策前期間（学習）】ターゲット(y) と 共変量(X) の関係性（βやノイズ）をじっくり学習
                                ↓ 施策実行  
【施策後期間（予測・比較）】施策後の共変量(X)の動きをベースに「もし施策がなかったら(ŷ)」を予測

「通常時の関係性」を学習する（施策前期間）：
まだ施策を行っていない綺麗な期間のデータを使って、「共変量 $X$ が1動くとき、ターゲット $y$ はどれくらい動くか（ $\beta$ ）」「ベースのトレンドは普段どれくらい上下ブレするのか（ $\sigma_\mu, \sigma_y$ ）」というデータの構造（法則）をMCMCサンプリング等によって学習します。
「もしも」の未来をシミュレーションする（施策後期間）：
施策が始まった後の期間に対して、「施策の影響を受けていない共変量 $X_t$ 」のデータだけをモデルに投入します。モデルはステップ1で学習した「通常時の法則」に従って、「もし施策がなかったら、この $X_t$ の動きからして、売上 $y_t$ はこれくらいになっているはずだ」という予測値 $\hat{y}_t$ （反事実）を叩き出します。
差分を因果効果とする：
こうして作った「もしもの予測値（青い点線）」と「現実の売上（黒い実線）」を引くことで、曜日要因や市場全体のトレンドといったノイズを完璧に相殺した「純粋な施策のインパクト（ $\tau_t$ ）」が手に入るのです。

検証用ダミーデータの生成コード

今回作成するデータは、先ほど解説したベイズ構造時系列モデル（BSTS）の構造をシンプルにした、以下の数式（データ生成プロセス）に基づいています。
施策前（ $t < 70$ ）のデータ生成式

\begin{aligned} X_t &\sim \mathcal{N}(10, 2^2) \\ y_t &= 2.5 X_t + \epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0, 1.5^2) \end{aligned}

ここでは、共変量 $X_t$ は平均10、標準偏差2の正規分布から独立に生成されるデータ（市場全体のトレンドなど）とします。ターゲットである自社の売上 $y_t$ は、この $X_t$ の値を 2.5倍したものに、標準偏差1.5のランダムな観測ノイズ $\epsilon_t$ が加わったものです。

施策後（ $t \ge 70$ ）のデータ生成式
70日目に施策（介入）が行われると、数式は以下のように変化します。

\begin{aligned} X_t &\sim \mathcal{N}(10, 2^2) \\ y_t &= 2.5 X_t + \text{True Effect} + \epsilon_t \\ &= 2.5 X_t + 10 + \epsilon_t \end{aligned}

ポイントは、共変量 $X_t$ の生成プロセスは施策前後で一切変わっていない（独立している）のに対し、ターゲット変数 $y_t$ にだけ一律で $+10$ の底上げ（True Effect）が加わっている点です。スクラッチ実装のゴールは、モデルに「 $y$ と $X$ のデータ」だけを渡し、この数式の裏に隠された「傾き $2.5$ 」と「施策効果 $+10$ 」という正解の数値を、MCMCによって正確に逆算（推論）できるかどうかを検証することになります。以下のコードでは、CausalImpactの前提条件である「施策の影響を受けないが、ターゲットと強い相関を持つ共変量（ $X$ ）」と、「施策によって途中でトレンドが跳ね上がるターゲット（ $y$ ）」をシミュレーションしています。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib

# 1. 乱数シードの固定（読者が同じ結果を再現できるようにするため）
np.random.seed(123)

# 2. 期間パラメータの設定
n_periods = 100          # 全期間（100日間）
intervention_point = 70  # 施策実行日（70日目）

# 3. 共変量（コントロール変数）X の作成
# 施策の影響を一切受けない、市場のトレンドや競合の動きなどを想定
X = np.random.normal(loc=10, scale=2, size=n_periods)

# 4. ターゲット変数 y のベース（施策がない場合の本来の推移）を作成
# X と強い相関を持たせつつ、ランダムなノイズ（観測誤差）を加える
y_base = 2.5 * X + np.random.normal(loc=0, scale=1.5, size=n_periods)
y = y_base.copy()

# 5. 施策効果（介入効果）の注入
# 70日目以降のターゲットデータにのみ、一律で「+10」の効果を加える
true_effect = 10
y[intervention_point:] += true_effect

# DataFrameにまとめる
df = pd.DataFrame({
    'y': y,
    'X': X
})

# 6. 生成したデータの可視化
plt.figure(figsize=(10, 5))
plt.plot(df.index, df['y'], label='Target (y: 自社の売上)', color='black', linewidth=2)
plt.plot(df.index, df['X'], label='Covariate (X: 競合・市場トレンド)', color='gray', linestyle='--')
plt.axvline(x=intervention_point, color='red', linestyle=':', linewidth=2, label='Intervention (施策実行日)')
plt.title('Generated Dummy Data for CausalImpact')
plt.xlabel('Time (Days)')
plt.ylabel('Value')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

ここでは、検証用に用いるダミーデータの生成と可視化を行っています。まず、np.random.seed(123) で乱数のシード値を固定し、実行のたびに同じデータが再現されるように設定します。全期間 n_periods = 100（100日間）とし、70日目（intervention_point = 70）を施策の実行日として定義します。共変量 X は、np.random.normal(loc=10, scale=2, size=n_periods) を用いて平均10、標準偏差2の正規分布から生成され、施策の影響を受けない市場トレンドなどを表現しています。ターゲット変数である自社の売上 y の本来の推移は、y_base = 2.5 * X + np.random.normal(loc=0, scale=1.5, size=n_periods) によって、共変量 X を2.5倍した値に、標準偏差1.5の観測誤差を加えることで構築しています。そして施策効果を反映するため、施策日である70日目以降のターゲットデータに対してのみ、y[intervention_point:] += true_effect で一律に +10 の底上げ（True Effect）を注入しています。最後に、生成されたターゲット変数 y と共変量 X を pd.DataFrame に格納し、matplotlib および日本語表記をサポートする japanize_matplotlib を用いて、施策実行日の位置（赤の点線）とともに可視化しています。

上記のコードを実行すると以下の結果が得られます。生成データ

PyMCによるモデル定義とMCMCサンプリング

PyMCを使ってベイズ構造時系列モデル（BSTS）を構築します。

ここでのポイントは「施策前（Pre-period）のデータのみをモデルに与える」ということです。施策の影響を受けていない純粋な期間のデータから、データ生成の「法則」をMCMCに学習させます。

モデル定義とサンプリングのコード

以下のコードを実行して、モデルの定義とMCMCサンプリング（パラメータの探索）を行います。

import pymc as pm
import numpy as np

# 1. データを施策前の期間（Pre-period）に絞り込む
y_pre = df['y'].iloc[:intervention_point].values
X_pre = df['X'].iloc[:intervention_point].values
T = len(y_pre)

# 2. ベイズ構造時系列モデル（BSTS）の定義
with pm.Model() as bsts_model:
    
    # --- 事前分布（Priors）の設定 ---
    # 観測ノイズ（データ自体のブレ）の標準偏差
    sigma_obs = pm.HalfNormal("sigma_obs", sigma=5)
    
    # トレンドノイズ（日々のトレンドの変動幅）の標準偏差
    # 時系列が急激に変化しすぎないよう、ややタイトな範囲を想定
    sigma_level = pm.HalfNormal("sigma_level", sigma=2)
    
    # 共変量 X の影響度（回帰係数 beta）
    beta = pm.Normal("beta", mu=0, sigma=5)
    
    # --- 状態空間の定義（非中心化パラメータ化） ---
    # サンプラーが「漏斗状の谷（ファンネル）」にハマって発散するのを防ぐため、
    # 状態方程式のランダムウォークを分解して記述します。
    
    # ① トレンドの初期値（スタート地点）
    mu_init = pm.Normal("mu_init", mu=y_pre[0], sigma=5)
    
    # ② 各時点におけるトレンドの「無次元化された変化量」
    # サンプラーが探索しやすい標準正規分布 N(0, 1) から独立にサンプリング
    mu_innovation = pm.Normal("mu_innovation", mu=0, sigma=1, shape=T-1)
    
    # ③ 変化量を sigma_level でスケールさせ、累積和（cumsum）をとってトレンドを構築
    mu_trend = pm.math.concatenate([
        [mu_init], 
        mu_init + pm.math.cumsum(mu_innovation * sigma_level)
    ])
    
    # --- 観測方程式（データ生成プロセス） ---
    # ターゲット変数の予測平均値 ＝ トレンド ＋ 共変量による効果
    nu = mu_trend + beta * X_pre
    
    # 尤度（実際の観測データ y_pre をモデルに適合）
    obs = pm.Normal("obs", mu=nu, sigma=sigma_obs, observed=y_pre)
    
    # --- MCMCサンプリングの実行 ---
    # 時系列モデルの複雑な幾何学的形状を確実に、かつ細かく探索させるため、
    # tune（慣らし運転）を1500回に増やし、target_accept（採択率の目標）を0.99と高く設定します
    trace = pm.sample(draws=1000, tune=1500, target_accept=0.99, random_seed=42, progressbar=True)

# パラメータの推定結果の要約を表示（MCMCが正解を捉えているか確認）
pm.summary(trace, var_names=["beta", "sigma_obs", "sigma_level"])

ここでは、PyMCを用いてベイズ構造時系列モデル（BSTS）を定義し、施策前期間（Pre-period）のデータだけを用いてMCMCサンプリングによるパラメータ推定を実行しています。まず、推定対象のデータとして、施策開始前の70日分（y_pre および X_pre）を抽出します。 with pm.Model() as bsts_model でPyMCのモデルブロックを開始し、観測ノイズの標準偏差 sigma_obs およびトレンドノイズの標準偏差 sigma_level に非負の事前分布である半正規分布 pm.HalfNormal を設定します。共変量の影響度を表す回帰係数 beta には、平均0、標準偏差5の正規分布 pm.Normal を設定しています。モデルの核となる「見えないトレンド mu_trend」の表現には、サンプリングを安定させるための「非中心化パラメータ化」を採用しています。初期値 mu_init に加えて、各時点の標準的なノイズ成分 mu_innovation を標準正規分布 pm.Normal(mu=0, sigma=1) から抽出し、それに sigma_level を乗じた累積和を pm.math.cumsum で計算することで、トレンドのランダムウォークを構成します。観測方程式として、ターゲット変数の平均 nu をトレンド mu_trend と共変量効果 beta * X_pre の和で表し、実際の観測データ y_pre が平均 nu、標準偏差 sigma_obs の正規分布に従う尤度 pm.Normal("obs", ...) として記述します。最後に、pm.sample を呼び出し、MCMCによるサンプリングを実行します。サンプリング空間の複雑な幾何学的形状に対応するため、tune=1500 でバーンイン期間（慣らし運転）を長めに確保し、target_accept=0.99 とすることで、サンプラーが谷間（ファンネル）に捕まることによるサンプリングの「発散（Divergences）」を完全に回避します。サンプリング完了後、pm.summary により各パラメータの推定要約（平均値、標準偏差、信用区間など）を表示させています。

① 状態空間の定義（非中心化パラメータ化）によるトレンドの表現
状態空間モデルの肝は、「時間とともに変化する見えないトレンド $\mu_t$ 」をどうモデル化するかです。このコードでは、時系列データが持つ「昨日の売上ベースから、今日も少しランダムに変動する（ランダムウォーク）」という性質をそのまま数式通りに記述するのではなく、「非中心化パラメータ化（Non-centered Parameterization）」という強力なテクニックを用いています。愚直に書いた場合の数式（中心化モデル）本来、ランダムウォークは以下のように「昨日の状態」をベースに次の状態の確率分布を定義します。

\mu_t = \mu_{t-1} + \eta_t, \quad \eta_t \sim \mathcal{N}(0, \sigma_\mu^2)

しかし、この数式の通りにモデル化してしまうと、各時点の $\mu_t$ が前後のパラメータと強力に結びついてしまいます。その結果、サンプリング空間の幾何学的な形状が「漏斗（ファンネル）状」の狭い谷になり、PyMCの標準アルゴリズムであるNUTS（No-U-Turn Sampler）がその谷間に入り込んで動けなくなり、サンプリングが発散（Divergences）するエラーを引き起こしやすくなります。非中心化パラメータ化の数式（実際にコード化したもの）そこで、この依存関係を断ち切るために、数式を以下のように分解して再構築します。

\begin{aligned} \tilde{\eta}_t &\sim \mathcal{N}(0, 1^2) \\ \mu_t &= \mu_0 + \sum_{i=1}^{t} (\sigma_\mu \times \tilde{\eta}_i) \end{aligned}

$\tilde{\eta}_t$ （mu_innovation）：各時点の「ブレ」を、前後の状態から完全に独立した、サンプラーが最も探索しやすい広々とした形状の「標準正規分布 $\mathcal{N}(0, 1)$ 」からサンプリングします。累積和（ $\sum$ / pm.math.cumsum）：巡回しやすい空間で集めた綺麗なブレの塊に対して、後からトレンドの変動の激しさ $\sigma_\mu$ （sigma_level）を掛け算し、初期値 $\mu_0$ （mu_init）から順に足し上げていくことで、数学的に全く等価なランダムウォーク構造を作ります。コードを見ると、この数式がそのまま再現されていることが分かります。

コードを見ると、この数式がそのまま再現されていることが分かります。
数式： $\tilde{\eta}_t \sim \mathcal{N}(0, 1)$

mu_innovation = pm.Normal("mu_innovation", mu=0, sigma=1, shape=T-1)

数式： $\mu_0 + \sum(\sigma_\mu * \tilde{\eta}_i)$

mu_trend = pm.math.concatenate([
    [mu_init], 
    mu_init + pm.math.cumsum(mu_innovation * sigma_level)
])

② なぜパラメータを「分布」で持つのか？
従来の機械学習のように「傾き $\beta = 2.5$ 」と一つの値を決める（点推定する）だけでは、データが持つ本来の「不確実性」を評価できません。MCMCは「 $\beta$ は 2.48〜2.52 の間にありそうだ」という確率分布（事後分布）を丸ごと出力してくれます。今回のコードで、サンプラーの設定を target_accept=0.99 に引き上げ、慣らし運転（tune=1500）を増やしたのも、この「分布の形」を極めて正確に描き出すためです。

非中心化パラメータ化によってNUTSのトラップを回避し、さらにサンプラーの歩幅を細かく（慎重に）調整したことで、あり得るパターンのバリエーション（有効サンプルサイズ）が十分に確保されます。これにより、パラメータの「ブレの幅」が正しく保証され、CausalImpactの最大の強みである「パラレルワールド（反事実）の予測における95%信用区間」の『帯の幅』を、一切の妥協なく数学的に正しく計算できるようになるのです。

上記のコードを実行すると以下の結果が得られます。実行結果施策前 MCMCサンプリングの結果得られた要約テーブル（pm.summary の出力）から、モデルがデータ生成プロセスの真の構造を極めて正確に捉えられていることが確認できます。

回帰係数 beta の推定値は平均値（mean）が 2.416 となっており、データ生成式で設定した真の値である 2.5 に非常に近い値を推定できています。94%信用区間（hdi_3% 〜 hdi_97%）も [2.266, 2.564] となっており、真の値をしっかりとその範囲内に内包しています。
観測ノイズ sigma_obs の推定平均値は 1.506 であり、これもデータ生成時にノイズとして加えた真の標準偏差 1.5 とほぼ一致しています。
トレンドノイズ sigma_level は平均値が 0.092 と小さく推定されています。これは、ダミーデータ生成の際、トレンドのランダムウォークによる変化（システムノイズ）を事実上ゼロ（y_base は共変量と観測ノイズのみで決定）として生成したため、モデルが「ベーストレンドはほぼ時間経過で変化しない（ほぼ一定である）」と正しく見抜いたことを意味します。
すべての主要パラメータについて、収束判定指標である r_hat（Gelman-Rubin統計量）が収束の目安である 1.0（厳密には1.05以下）となっており、MCMCサンプリングのチェーンが十分に収束し、信頼性の高い事後分布が得られていることが保証されています。

施策後の反事実シミュレーションと因果効果の可視化

MCMCサンプリングによって、「通常時におけるデータ生成の法則（ $\beta, \sigma_\mu, \sigma_y$ ）」の事後分布を手に入れることができました。

最後のステップでは、この数千パターンのパラメータの組み合わせを使って、施策後の期間（70日目〜100日目）に対して「もし施策がなかったら（反事実）」のシミュレーションを数千通り行います。これにより、予測値の平均だけでなく、「95%信用区間（予測のばらつきの幅）」を数学的に正しく導き出せます。

反事実のシミュレーションと効果測定のコード

MCMCサンプリングによって無事にパラメータの事後分布が得られたら、次はいよいよ施策後期間（ $t \ge 70$ ）における「もし施策がなかったら（反事実）」のシミュレーションを行います。

この未来のタイムライン（パラレルワールド）を1日ずつ描き出す予測ロジックは、最初に定義した状態空間モデルの数式そのものです。
反事実予測のシミュレーション数式施策後の各時点 $t$ において、モデルは以下の2つのステップを1日ずつ未来に向かって繰り返します。

\begin{aligned} \mu_t^{(i)} &= \mu_{t-1}^{(i)} + \eta_t^{(i)}, \quad \eta_t^{(i)} \sim \mathcal{N}\left(0, (\sigma_\mu^{(i)})^2\right) \\ \hat{y}_t^{(i)} &= \mu_t^{(i)} + \beta^{(i)} X_t + \epsilon_t^{(i)}, \quad \epsilon_t^{(i)} \sim \mathcal{N}\left(0, (\sigma_y^{(i)})^2\right) \end{aligned}

$(i)$ は、MCMCで得られた数千個あるパラメータサンプルのうちの「 $i$ 番目のパターン」であることを意味しています。
ステップ1（状態の更新）：昨日のトレンド $\mu_{t-1}$ に、学習されたトレンドの変動幅 $\sigma_\mu$ を持つシステムノイズ $\eta_t$ を足し算して、今日の潜在トレンド $\mu_t$ をシミュレーションします。
ステップ2（観測値の生成）：今日の潜在トレンド $\mu_t$ に、学習された係数 $\beta$ と現実の（施策の影響を受けていない）共変量 $X_t$ を掛け合わせたものを足し、最後に観測ノイズ $\epsilon_t$ を乗せて、施策がなかった場合の予測値 $\hat{y}_t$ （反事実）を1つ生成します。

これをすべてのMCMCサンプル（数千パターン）× 施策後の期間（30日間）の分だけループ処理します。
実際のコードを見ると、この2つの数式ステップが完全にそのまま移植されていることが分かります。

数式： $\mu_t = \mu_{t-1} + \eta_t, \eta_t ~ N(0, \sigma_mu^2)$
① 状態の更新：前日のトレンドにノイズを加えて今日のトレンドを生成

current_mu = np.random.normal(current_mu, sigma_level_samples[i])

数式： $\hat{y}_t = \mu_t + \beta X_t + \epsilon_t, \epsilon_t ~ N(0, \sigma_y^2)$
② 観測値の生成：今日のトレンド＋今日の共変量効果＋観測ノイズ

y_post_pred_samples[i, t] = np.random.normal(
    current_mu + beta_samples[i] * X_post[t], 
    sigma_obs_samples[i]
)

この「数式通りのシミュレーションを何千回も泥臭く繰り返す」ことによって、あり得る未来のバリエーションがすべて配列（y_post_pred_samples）に蓄積されます。あとはその平均値やパーセンタイル（2.5%点、97.5%点）を間引くだけで、綺麗な「95%信用区間」が計算できるわけです。

それでは、このシミュレーションから効果測定、そしてグラフ描画までを行う全体のコードを見てみましょう。

import numpy as np
import matplotlib.pyplot as plt

# 1. MCMCの事後分布から各パラメータの全サンプルを抽出
post = trace.posterior
beta_samples = post["beta"].values.flatten()
sigma_obs_samples = post["sigma_obs"].values.flatten()
sigma_level_samples = post["sigma_level"].values.flatten()

# 施策直前（Pre期間の最後）のトレンド状態を抽出
mu_last_samples = post["mu_init"].values.flatten() + post["mu_innovation"].values.sum(axis=-1).flatten()

n_samples = len(beta_samples)
X_post = df['X'].iloc[intervention_point:].values
y_post_actual = df['y'].iloc[intervention_point:].values
n_post = len(X_post)

# 反事実の予測結果を格納する配列（サンプル数 × 施策後の期間）
y_post_pred_samples = np.zeros((n_samples, n_post))

# 2. 全てのMCMCサンプルに対して、施策後の未来（パラレルワールド）を1日ずつシミュレーション
for i in range(n_samples):
    current_mu = mu_last_samples[i]
    
    for t in range(n_post):
        # ① 状態の更新：前日のトレンドにノイズを加えて今日のトレンドを生成
        current_mu = np.random.normal(current_mu, sigma_level_samples[i])
        
        # ② 観測値の生成：今日のトレンド ＋ 今日の共変量効果 ＋ 観測ノイズ
        y_post_pred_samples[i, t] = np.random.normal(
            current_mu + beta_samples[i] * X_post[t], 
            sigma_obs_samples[i]
        )

# 3. 予測の平均値と95%信用区間（下側2.5%, 上側97.5%）を算出
y_post_mean = y_post_pred_samples.mean(axis=0)
y_post_lower, y_post_upper = np.percentile(y_post_pred_samples, [2.5, 97.5], axis=0)

# 4. 点推定効果の算出（実際の値 - 反事実のシミュレーション値）
point_effect_samples = y_post_actual - y_post_pred_samples
point_effect_mean = point_effect_samples.mean(axis=0)
point_effect_lower, point_effect_upper = np.percentile(point_effect_samples, [2.5, 97.5], axis=0)

# 5. 累積効果の算出（時間経過とともに効果を足し上げる）
cumulative_effect_samples = point_effect_samples.cumsum(axis=1)
cumulative_effect_mean = cumulative_effect_samples.mean(axis=0)
cumulative_effect_lower, cumulative_effect_upper = np.percentile(cumulative_effect_samples, [2.5, 97.5], axis=0)

# ==========================================
# 6. オリジナルを再現した3部構成グラフの描画
# ==========================================
time_post = np.arange(intervention_point, n_periods)
time_all = np.arange(n_periods)

fig, axes = plt.subplots(3, 1, figsize=(11, 12), sharex=True)

# グラフ1: 実際の観測値と反事実の予測（95%信用区間の帯付き）
axes[0].plot(time_all, df['y'], label='Actual (施策ありの現実)', color='black', linewidth=1.5)
axes[0].plot(time_post, y_post_mean, label='Counterfactual (予測: 施策なし)', color='blue', linestyle='--')
axes[0].fill_between(time_post, y_post_lower, y_post_upper, color='blue', alpha=0.15, label='95% Credible Interval')
axes[0].axvline(x=intervention_point, color='red', linestyle=':', linewidth=2, label='Intervention (施策日)')
axes[0].set_title('1. 実際の値と反事実の予測比較')
axes[0].set_ylabel('Target Value')
axes[0].legend(loc='upper left')
axes[0].grid(True, alpha=0.3)

# グラフ2: 各時点における施策効果（Point Effect）
axes[1].plot(time_post, point_effect_mean, color='blue', label='Estimated Effect (推定効果)')
axes[1].fill_between(time_post, point_effect_lower, point_effect_upper, color='blue', alpha=0.15)
axes[1].axhline(y=0, color='black', linestyle='-')
axes[1].axhline(y=true_effect, color='green', linestyle=':', linewidth=2, label='True Effect (+10)')
axes[1].axvline(x=intervention_point, color='red', linestyle=':', linewidth=2)
axes[1].set_title('2. 各時点における施策効果 (Point Effect)')
axes[1].set_ylabel('Effect Size')
axes[1].legend(loc='upper left')
axes[1].grid(True, alpha=0.3)

# グラフ3: 累積の施策効果（Cumulative Effect）
axes[2].plot(time_post, cumulative_effect_mean, color='blue', label='Cumulative Effect (累積効果)')
axes[2].fill_between(time_post, cumulative_effect_lower, cumulative_effect_upper, color='blue', alpha=0.15)
axes[2].axvline(x=intervention_point, color='red', linestyle=':', linewidth=2)
axes[2].set_title('3. 累積の施策効果 (Cumulative Effect)')
axes[2].set_ylabel('Cumulative Size')
axes[2].set_xlabel('Time (Days)')
axes[2].legend(loc='upper left')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

ここでは、MCMCサンプリングによって得られたパラメータの事後分布を用いて、施策後の期間（70日目〜100日目）における「もし施策を行っていなかったら（反事実）」のシミュレーションと、因果効果の算出・可視化を行っています。

パラメータサンプルの抽出: MCMCトレースからパラメータの事後分布サンプルをフラットな配列として取り出します。回帰係数 beta_samples、観測ノイズの標準偏差 sigma_obs_samples、トレンドノイズの標準偏差 sigma_level_samples を取得します。施策直前（69日目）における潜在トレンドの最終状態 mu_last_samples は、トレンドの初期値サンプルの配列と、それまでのトレンド変化量 mu_innovation の総和を足し合わせることで、各サンプリングパスごとに計算します。
反事実シミュレーションのループ: MCMCサンプルの総数分だけループ（for i in range(n_samples)）を回し、各サンプルにおける「施策がなかった場合の未来」を1日ずつ時系列順にシミュレーションします。各日のステップにおいて、まず np.random.normal(current_mu, sigma_level_samples[i]) によりトレンドをランダムウォークさせて更新（状態の更新）します。その後、更新されたトレンドに、そのサンプルの回帰係数 beta_samples[i] と施策後期間の実際の共変量 X_post[t] を掛け合わせたものを加え、さらに観測誤差 sigma_obs_samples[i] のノイズを上乗せした正規分布から値を生成します。これにより、ノイズも含めた「あり得たかもしれない観測値（反事実）」を1日ずつ、すべてのサンプル分蓄積します。
効果量の算出:
- 点推定効果（Point Effect）: 実際の観測値 y_post_actual から、シミュレーションした反事実の予測サンプルを差し引く（y_post_actual - y_post_pred_samples）ことで、各時点における施策効果を計算します。
- 累積効果（Cumulative Effect）: 点推定効果を時間経過とともに累積加算（cumsum(axis=1)）し、施策開始から現在までの効果の積み上げ量を算出します。
95%信用区間の計算とプロット: 各指標に対して np.percentile(..., [2.5, 97.5]) を用いることで、MCMCサンプルに基づく「予測の95%信用区間（下限と上限）」を算出します。最後に、matplotlib で3行1列のグラフを作成し、「実際の値と予測の比較」「各時点の施策効果（Point Effect）」「累積効果（Cumulative Effect）」をプロットしています。fill_between で描画された青い半透明の領域（信用区間）が、未来に進むにつれて不確実性の高まり（ランダムウォークの拡散）によって広がっていく様子が確認できます。

実行結果上記のコードを実行すると以下の結果が得られます。実行結果施策後

① 「正解」を完璧に見抜いたグラフ2
グラフ2の点線で示した True Effect (+10) と、MCMCシミュレーションが叩き出した Estimated Effect（青線）を見比べてみてください。モデルには「＋10の施策を行った」という情報は一切与えていないにもかかわらず、見事に＋10のラインを中心とした効果を逆算できていることが視覚的に分かります。

② 不確実性を表す「青い帯（95%信用区間）」の意味
グラフが未来（右側）に進むにつれて、青い帯の幅がじわじわと広がっていることに注目してください。これは、時系列モデルのランダムウォークが「未来に行けば行くほど、予測の不確実性が増してブレが大きくなる」という時間の性質を数理的に正しく表現できている証拠です。

単なる前後の平均値の比較では、「その効果が偶然なのか、それとも本当に施策のおかげなのか」の区別がつきません。しかし、このスクラッチ実装によって、「95%の確率の帯がゼロのライン（黒い横線）を完全に上回っているから、この施策は統計的に有意に効果があった」と自信を持って言えるロジックが完成しました。

causalimpact ライブラリでの実装

ここまでCausalImpactの数理ロジックとMCMCによるスクラッチ実装を解説してきましたが、実務の現場で毎回モデルを1から組むのは大変です。

実は、PythonにはGoogleのオリジナル（R言語版）を移植した causalimpact という強力なライブラリが存在します。これを使うと、私たちが裏側で実装した複雑な状態空間モデルの構築や反事実のシミュレーション、さらには効果の統計的有意差のレポート出力までを、わずか数行で実行してくれます。

ライブラリのインストール

事前に以下のコマンドでパッケージをインストールしておきます。

pip install pycausalimpact

ライブラリを使った効果測定コード

上記で作ったダミーデータ（df）をそのまま使って、ライブラリで効果測定を行ってみましょう。

from causalimpact import CausalImpact
import matplotlib.pyplot as plt

# 1. 施策前（Pre）と施策後（Post）の期間を指定し、モデルを学習させる
pre_period = [0, intervention_point - 1]
post_period = [intervention_point, n_periods - 1]
ci = CausalImpact(df, pre_period, post_period)

# ==========================================================
# 💡 Tips: y軸のスケールを最適化して見やすいグラフにするための処理
# ==========================================================
# ci.plot() は実行と同時にグラフを閉じてしまう仕様があるため、
# 一時的に plt.show() を「何もしない関数」にすり替えてミュートします。
original_show = plt.show
plt.show = lambda: None

# 2. グラフを生成（ミュート中のため、画面には出力されず裏側に保持されます）
ci.plot()

# 3. 裏側に保持されているグラフのパーツ（FigureとAxes）を取得
fig = plt.gcf()
axes = fig.get_axes()

# 4. 0日目の散漫初期化による極端なブレを画面外にカットし、適切な表示範囲に調整
y_min, y_max = df['y'].min(), df['y'].max()
margin = (y_max - y_min) * 0.2

# 1枚目（実際の値と予測）を現実的なデータ幅に設定
axes[0].set_ylim(y_min - margin, y_max + margin) 
# 2枚目（各時点の施策効果）を、想定される効果量の周辺にクローズアップ
axes[1].set_ylim(-20, 30)                        

# 5. 本来の表示機能を元に戻し、スケール調整済みの綺麗なグラフを出力
plt.show = original_show
plt.show()
# ==========================================================

# 6. 効果測定の結果（統計的な有意差など）をまとめたテキストレポートの出力
print(ci.summary())

ここでは、Pythonの pycausalimpact ライブラリを用いて、同様の因果効果測定をより簡潔なコードで実装しています。

モデルの初期化と学習: CausalImpact(df, pre_period, post_period) クラスをインスタンス化するだけで、ライブラリ内部で自動的にベイズ構造時系列モデルが構築・学習されます。引数には、元のデータフレーム df、施策前の期間インデックスを示すリスト pre_period（[0, 69]）、施策後の期間インデックス post_period（[70, 99]）を渡します。デフォルトでは、ターゲット変数（データフレームの最初の列、ここでは y）以外の列（ここでは X）が自動的に共変量（コントロール変数）として認識されます。
グラフの描画と表示範囲のカスタマイズ: 通常は ci.plot() を実行するだけで、CausalImpactが持つ「3部構成の可視化グラフ」が自動生成されます。しかし、ライブラリの仕様として、モデル構築の初期段階（特に0日目付近）における事前分布の初期化ブレが大きく、そのまま描画するとグラフ全体の縦軸スケールが引っ張られて非常に見づらくなってしまいます。このコードでは、plt.show を一時的にダミー関数に置き換えて描画処理を裏側に閉じ込めた上で、fig.get_axes() から各サブプロットの軸オブジェクトを取得しています。ターゲット変数の実データ幅に合わせて、1枚目のプロットのY軸表示範囲（axes[0].set_ylim）や、2枚目のプロットのY軸表示範囲（axes[1].set_ylim(-20, 30)）を適切にクローズアップすることで、初期値の大きなブレによる描画崩れをカットした非常に見やすいグラフを再構築しています。
レポートの出力: 最後に print(ci.summary()) を実行することで、推定された施策効果の平均値や累積値、それぞれの標準偏差や95%信用区間、そして効果の統計的有意性に関する要約データをテキスト形式で出力しています。

実行結果上記のコードを実行すると以下の結果が得られます。ライブラリ実行結果

Posterior Inference {Causal Impact}
                          Average            Cumulative
Actual                    34.2               1026.04
Prediction (s.d.)         24.61 (0.42)       738.42 (12.68)
95% CI                    [23.77, 25.42]     [712.98, 762.69]

Absolute effect (s.d.)    9.59 (0.42)        287.62 (12.68)
95% CI                    [8.78, 10.44]      [263.34, 313.06]

Relative effect (s.d.)    38.95% (1.72%)     38.95% (1.72%)
95% CI                    [35.66%, 42.4%]    [35.66%, 42.4%]

Posterior tail-area probability p: 0.0
Posterior prob. of a causal effect: 100.0%

For more details run the command: print(impact.summary('report'))

出力されたテキストレポート（ci.summary()）は、ベイズ統計に基づく推定結果を「平均値（Average）」と「累積値（Cumulative）」の2つの側面から要約して示しています。

観測値と予測値（Actual & Prediction）: 施策後期間における実際の観測値（Actual）の平均値は 34.2 であったのに対し、施策が行われなかった場合の予測値（Prediction）は平均 24.61 となり、95%信用区間（95% CI）は [23.77, 25.42] であると推定されています。
絶対効果（Absolute effect）: 実際の値と予測値の差分である「純粋な施策効果」は、1日あたり平均 9.59（95%信用区間：[8.78, 10.44]）と推定されました。データ生成時に私たちが設定した真の効果（True Effect = 10）を、極めて狭い信用区間の幅の中に高精度で捉えられていることがわかります。また、30日間の累積効果（Cumulative）は計 287.62（95%信用区間：[263.34, 313.06]）の売上底上げに寄与したことを示しています。
相対効果（Relative effect）: 施策が行われなかった場合のベースライン予測と比較して、施策によって売上が +38.95%（95%信用区間：[35.66%, 42.4%]）増加したことを示しており、施策のインパクトを直感的な割合で評価できます。
有意性の判定（Posterior tail-area probability & prob. of a causal effect）: ベイズ的な片側p値に相当する Posterior tail-area probability p が 0.0 となっており、施策による効果が偶然得られる確率は極めて低いことを意味します。また、因果効果が存在する確率（Posterior prob. of a causal effect）は 100.0% と判定されており、この施策に統計的に有意な効果があったと自信を持って結論付けることができます。

詳細な文章レポートの自動生成

さらに、このライブラリの非常に面白い機能として、ci.summary(output='report') を実行すると、「この施策には統計的な有意差があったのか、累積でどれくらいのインパクトがあったのか」を人間が読める英語の文章で自動生成してくれます。

ci.summary(output='report')を実行した結果

Analysis report {CausalImpact}


During the post-intervention period, the response variable had
an average value of approx. 34.2. By contrast, in the absence of an
intervention, we would have expected an average response of 24.61.
The 95% interval of this counterfactual prediction is [23.85, 25.4].
Subtracting this prediction from the observed response yields
an estimate of the causal effect the intervention had on the
response variable. This effect is 9.59 with a 95% interval of
[8.8, 10.36]. For a discussion of the significance of this effect,
see below.


Summing up the individual data points during the post-intervention
period (which can only sometimes be meaningfully interpreted), the
response variable had an overall value of 1026.04.
By contrast, had the intervention not taken place, we would have expected
a sum of 738.42. The 95% interval of this prediction is [715.37, 761.89].


The above results are given in terms of absolute numbers. In relative
terms, the response variable showed an increase of +38.95%. The 95%
interval of this percentage is [35.77%, 42.07%].


This means that the positive effect observed during the intervention
period is statistically significant and unlikely to be due to random
fluctuations. It should be noted, however, that the question of whether
this increase also bears substantive significance can only be answered
by comparing the absolute effect (9.59) to the original goal
of the underlying intervention.


The probability of obtaining this effect by chance is very small
(Bayesian one-sided tail-area probability p = 0.0).
This means the causal effect can be considered statistically
significant.

(訳)
分析レポート {因果効果}

介入後期間において、応答変数の平均値は約34.2でした。一方、介入がなかった場合、平均応答値は24.61になると予測されました。

この反事実的予測値の95%信頼区間は[23.85, 25.4]です。

この予測値を観測された応答値から差し引くことで、介入が応答変数に及ぼした因果効果の推定値が得られます。この効果は9.59で、95%信頼区間は[8.8, 10.36]です。この効果の有意性については、以下を参照してください。

介入後期間における個々のデータポイントを合計すると（必ずしも意味のある解釈ができるとは限りませんが）、応答変数の全体値は1026.04となりました。

対照的に、介入が行われなかった場合、合計値は738.42になると予想されました。この予測値の95%信頼区間は[715.37, 761.89]です。

上記の結果は絶対値で示されています。相対値で見ると、応答変数は+38.95%増加しました。この割合の95%信頼区間は[35.77%, 42.07%]です。

これは、介入期間中に観察された正の効果が統計的に有意であり、ランダムな変動によるものではない可能性が高いことを意味します。ただし、この増加が実質的な意義を持つかどうかは、絶対値（9.59）を介入の当初の目標値と比較することによってのみ判断できることに留意する必要があります。

この効果が偶然に得られる確率は非常に小さいです（ベイズ片側裾面積確率p = 0.0）。

これは、因果関係が統計的に有意であると考えられることを意味します。

実務における注意点

「今回のスクラッチ実装では、ロジックの直感的な理解のためにサンプリング後の未来をPythonのforループで1日ずつシミュレーションしました。実務でデータ量が膨大な場合や、より複雑な季節性を扱う場合は、本家Googleのアルゴリズムのようにシミュレーション・スムーザー（Simulation Smoother）等を用いて状態そのものを一括でサンプリングするアプローチが取られます」
「実務のデータに適用する場合、曜日ごとの変動（7日周期）などの『季節性』を考慮する必要があります。数式としては観測方程式に $\gamma_t$ （季節成分）を追加し、状態方程式で『過去7日間の合計が0になる』といった制約を設けることで、曜日ノイズを綺麗に吸収した反事実予測が可能になります」
「実務でコントロール変数の選定や予測精度に自信が持てないときは、施策を実行していない過去の期間に対してあえてCausalImpactをかけてみる『プラセボテスト』が有効です。何もしていない期間で効果が『0（有意差なし）』と正しく出力されれば、そのモデルは通常時のトレンドを完璧に捉えられている（＝施策後の反事実予測も信頼できる）という強力な証明になります」

まとめ

本記事では、施策や介入による影響を「反事実（もし施策を行っていなかったら）」との比較によって高精度に評価するCausalImpact（因果効果測定） について解説しました。

記事を通じて、以下のステップを実践・学習しました。

CausalImpactの理論と数理の理解: ベイズ構造時系列モデル（BSTS）における「観測方程式」と「状態方程式」の役割、および「反事実」を用いた引き算による効果測定のロジックを学びました。
検証用データのシミュレーション: 施策日以降に一律で効果（+10）が加わるターゲット変数と、施策に影響されない共変量のダミーデータをPythonで生成しました。
PyMCによるスクラッチ実装: 状態の推移におけるサンプラーの「発散（Divergences）」を防ぐため、「非中心化パラメータ化（Non-centered Parameterization）」を施したモデルを構築し、MCMCによる安定したサンプリングを実証しました。
反事実のシミュレーションと可視化: MCMCの事後分布パラメータを泥臭くループ処理して施策なしの未来をシミュレーションし、3部構成（現実と予測の比較、点推定効果、累積効果）のグラフと「95%信用区間（予測の帯）」を再現しました。
pycausalimpact ライブラリの実用: 実務で強力なツールとなるライブラリ版の実装方法、初期化ブレを抑えるためのY軸のスケール調整テクニック、および自動生成される詳細なサマリーレポートの読み方を解説しました。

CausalImpactは、曜日要因や市場トレンドといった外部ノイズを綺麗に相殺し、予測が持つ「不確実性」までも信用区間という形で定量化できる極めて強力なフレームワークです。実務におけるマーケティング施策や新機能リリースの効果測定など、確度の高い意思決定が求められる局面でぜひ役立ててください。

本記事の文章・構成の一部に生成AIを使用しています。

Causal Impactの概要​

CausalImpactの基本ロジック：「反事実」の予測​

反事実を予測する「ベイズ構造時系列モデル（BSTS）」の数理​

検証用ダミーデータの生成コード​

PyMCによるモデル定義とMCMCサンプリング​

モデル定義とサンプリングのコード​

施策後の反事実シミュレーションと因果効果の可視化​

反事実のシミュレーションと効果測定のコード​

causalimpact ライブラリでの実装​

ライブラリのインストール​

ライブラリを使った効果測定コード​

詳細な文章レポートの自動生成​

実務における注意点​

まとめ​