@Maxwell_110 INFERRING CAUSAL IMPACT USING BAYESIAN STRUCTURAL TIME-SERIES MODELS Paper Summary Ver. 2024

01 Introduction 02 Bayesian structural time-series models 03 Application to simulated data 04 Application to empirical data

01 ◼ Causal Impact とは，Google による因果推論のフレームワーク package: "CausalImpact" ◼ A/B テストなどの RCT ができない環境下において，反実仮想的予測を行い介入効果（(A)TT: (Average) Treatment Effect on the Treated）を推定する ◼ 時系列ベイズ状態空間モデル： Simulation Smoother + Gibbs sampling にて，状態 および パラメータ をサンプリング ◼ 局所的トレンド・季節性・共変量などを組み込み可能 ◼ Spike and slab 事前分布による共変量の選択機能 ◼ Bayesian Model Averaging による過学習の抑制 ◼ 近年は公式以外の実装もある模様 ⇒ https://tjo.hatenablog.com/entry/2023/12/04/170000 Introduction Resources: https://research.google/pubs/pub41854/ https://youtu.be/GTgZfCltMm8 http://google.github.io/CausalImpact/ Summary

Causal Impact は 「反実仮想（counterfactual）の時系列を 扱えるように DID を一般化したもの」と謳っている Causal Effect in DID (traditional DID) Difference In Difference (DID) Causal Impact (CI) post-intervention pre-intervention post-intervention pre-intervention Causal Effect at any given time point Counterfactual time series Observed time series Treatment group Control group • Parallel trends assumption • Synthetic control (Abadie et al., 2010) 𝑌𝑖 = 𝛽0 + 𝛽1 𝑇𝑟𝑒𝑎𝑡𝑖 + 𝛽2 𝑃𝑜𝑠𝑡𝑖 + 𝜷𝟑 𝑇𝑟𝑒𝑎𝑡𝑖 × 𝑃𝑜𝑠𝑡𝑖 + 𝛽4 𝐶𝑜𝑣𝑎𝑟𝑖𝑎𝑡𝑒𝑠𝑖 + 𝑒𝑖 実際には，DID の中でも，特に Synthetic Control に inspire されているといって良さそう 論文中で使用される "control" という word は，普通の統制群の意味ではなく，この synthetic control の control の意に近いと思う

Synthetic Control ◼ Causal Impact と目的を同じくする手法で，Causal Impact は Synthetic Control に inspire されているといっても 良さそう ◼ A/B test 等の場合と異なり非介入群が存在しない場合，介入をうけていない reference data から「任意の時系列の 加重平均で求まる，介入群と似た性質をもつ時系列」を生成して比較に使用する手法 => Causal Impact と目的は同 じで，統制群がなくても impact を推定したい時に使う ◼ 例えば，Abadie and Gardeazabal, 2003 はスペインの Basque 州において発生したテロが経済に与えた影響を， synthetic control で推定した．その際に，2 つの別のスペインの地域の経済の時系列を組み合わせて，統制群を作 成した．他にも US の California における禁煙キャンペーンの影響を synthetic control で推定した例もある (Abadie et al., 2010) ◼ 擬似統制群の構築に使用可能な情報 1. 介入前の目的変数そのもの（トレンド・季節性など） 2. 他の時系列 - 介入の影響が無い - 介入前の目的変数に対する予測力がある 3. 事前分布に使用可能な知識（過去の研究結果など） これらは，Causal Impact でも同様に使用可能 Using units from the rest of the U.S. to make a synthetic control. A synthetic control matches to the interved time-seriese. Abadie et al., 2010 (Abadie et al., 2010; Abadie and Gardeazabal, 2003)

Limitations for Difference In Difference ◼ 通常の DID は，data が i.i.d. であることを前提とした静的な回帰モデルを 使用するため，推定した係数の標準誤差の過小評価に伴う Type I error が大 きい ◼ 大部分の DID は，介入前後の二時点のみの比較であり，介入効果がその後ど うなっていくかを分析しない（できない） ◼ DID の枠組み内で，統制群が存在しない場合に使用可能な "Synthetic Control" は，統制群の構成の仕方に数理的な制約がある．詳しくは以下参照． (Abadie et al., 2010; Abadie and Gardeazabal, 2003)

行列表示ではなく，季節性をもつ場合の式を陽に表すと・・・ 𝑦𝑡 = 𝜇𝑡 + 𝛾𝑡 + 𝛽𝑇 𝑥𝑡 + 𝜀𝑡 𝜇𝑡 = 𝜇𝑡−1 + 𝛿𝑡−1 + 𝜂𝜇, 𝑡−1 𝛿𝑡 = 𝛿𝑡−1 + 𝜂𝛿, 𝑡−1 𝛾𝑡 = − σ𝑠=1 𝑆−1 𝛾𝑡−𝑠 + 𝜂𝛾, 𝑡−1 状態空間モデルの一般式（行列表示） 観測方程式： 𝑦𝑡 = 𝑍𝑡 𝑇 𝛼𝑡 + 𝜀𝑡 𝜀𝑡 ~ 𝑁 0, 𝜎𝑡 2 状態方程式： 𝛼𝑡+1 = 𝑇𝑡 𝛼𝑡 + 𝑅𝑡 𝜂𝑡 𝜂𝑡 ~ 𝑁 0, 𝑄𝑡 Causal Impact の実体は状態空間モデル ◼ AR や ARIMA などの時系列モデルより表現力が高い （例えば AR や ARIMA は状態空間モデルの特殊形） ◼ 計算コストが低い bsts は simulation smoother (Durbin et al., 2002) + Gibbs sampling で推定 状態方程式 観測方程式 線形ガウス状態空間モデル 02 Bayesian structural time-series models 𝑦𝑡 : observed data 𝛼𝑡 : state vector 𝑍𝑡 : output vector 𝑇𝑡 : transition matrix 𝑅𝑡 : control matrix 𝜀𝑡 : scalar observation error 𝜎𝑡 2: noise variance 𝜂𝑡 : system error 𝑄𝑡 : state-diffusion matrix

Basic structural model Modular ・ 線形モデルであるため，コンポーネント 2 や 3 を加え たり外したりしてモデル構築が可能 Fully Bayesian Approach ・ 全コンポーネントをベイズの枠組みの中で計算可能 𝛾𝑡 = − σ𝑠=1 𝑆−1 𝛾𝑡−𝑠 + 𝜂𝛾, 𝑡 𝜂𝛾, 𝑡 ~ 𝑁 0, 𝜎𝛾 2 周期 S の季節性 春夏秋冬であれば S = 4 でいずれか一つの季節は reference（係数は 0） 𝛾𝑡+1 𝛾𝑡 𝛾𝑡−1 = −1 −1 −1 1 0 0 0 1 0 𝛾𝑡 𝛾𝑡−1 𝛾𝑡−2 + 𝜂𝛾, 𝑡 0 0 静的（static）・動的（dynamic）な係数をもつ共変量 「(spill-over 効果なども含めた) 介入の影響を受けない」共変量 𝑋𝑡 や係数 𝛽 が時間とともに変化する「動的係数 𝛽𝑡 」の設定が可能 𝛽𝑇𝑋𝑡 → 𝛽𝑡 𝑇𝑋𝑡 = σ 𝑗=1 𝐽 𝛽𝑗, 𝑡 𝑋𝑗, 𝑡 𝛽𝑗, 𝑡 = 𝛽𝑗, 𝑡−1 + 𝜂𝛽, 𝑗, 𝑡−1 𝜂𝛽, 𝑗, 𝑡 ~ 𝑁 0, 𝜎𝛽𝑗 2 Random walk → dynamic ローカル線形トレンド（確率的トレンド） トレンド成分 𝛿𝑡 が変化可能．例えば，最初は増加でその後，減少トレンドに転じるような場合も表現可能 𝜇𝑡 = 𝜇𝑡−1 + 𝛿𝑡−1 + 𝜂𝜇, 𝑡−1 𝜂𝜇, 𝑡 ~ 𝑁 0, 𝜎𝜇 2 𝛿𝑡 = 𝛿𝑡−1 + 𝜂𝛿, 𝑡−1 𝜂𝛿, 𝑡 ~ 𝑁 0, 𝜎𝛿 2 𝑦𝑡 = 𝜇𝑡 + 𝛾𝑡 + 𝛽𝑇𝑥𝑡 + 𝜀𝑡 1 2 3 1 2 3 Random walk →

1. 介入前において，共変量と介入群の関係性が安定 => 静的係数 => Forward-filtering と backward-sampling の枠組みの中で，spike-and-slab 事前分布（後述）を効率的に実装 できるため，変数選択による過学習の抑制が期待できる 2. 時間の経過に伴い共変量と介入群の線形関係が変化すると考えられる => 動的係数 => 共変量候補が多いと計算コストが大きくなってしまうが，いくつか代替案を述べている (Nakajima and West, 2013) => しかし，bsts の CRAN manual の bsts::add.dynamic.regression をみる限りでは計算コストの低減対策はさ れておらず，他の状態変数と同様の扱いの模様 介入効果を受けない共変量を多く用意すること自体が難しいので，計算コストはそこまで心配しなくて良い？ 共変量に対して静的係数と動的係数のどちらを選択すべきか？

Example: Local Linear Trend model with a static regression component 𝑦𝑡 = 𝜇𝑡 + 𝛽𝜌 𝑇𝑥𝑡 + 𝜀𝑡 𝜀𝑡 ~ 𝑁 0, 𝜎𝑦 2 𝜇𝑡 = 𝜇𝑡−1 + 𝛿𝑡−1 + 𝜂𝜇, 𝑡 𝜂𝜇, 𝑡 ~ 𝑁 0, 𝜎𝜇 2 𝛿𝑡 = 𝛿𝑡−1 + 𝜂𝛿, 𝑡 𝜂𝛿, 𝑡 ~ 𝑁 0, 𝜎𝛿 2 1 𝜎2 ~ 𝐺𝑎𝑚𝑚𝑎 𝜈 2 , 𝑠 2 共変量 𝛽 に関する項を除けば，MCMC でサンプリングするパラメータは観測・状態方程式内の各誤差項の 𝜎2 𝜎2 の事前分布は，典型的なものとして（正規分布に対して共役な）逆ガンマ分布を採用することが多い 𝜈 や 𝑠 の値の設定の仕方は論文を参照 (基本的には 𝑦𝑡 の不偏分散の値を利用して設定する) Brodersen et al., 2015 Figure. 2 b a a 状態方程式 (ローカル線形トレンド) b 観測方程式

共変量の選択機能をもつ spike-and-slab 事前分布 𝑝 𝛽, 𝜎𝜀 2, 𝜚 = 𝑝 𝛽𝜚 |𝜚, 𝜎𝜀 2 𝑝 𝜎𝜀 2|𝜚 𝑝 𝜚 Spike: 𝛽𝑗 が 0 になるかどうかを表す二項分布 𝑝 𝜚 = ς 𝑗=1 𝐽 𝜋 𝑗 𝜚𝑗 1 − 𝜋𝑗 1−𝜚𝑗 Indicator vector: 𝜚 = 𝜚1 , 𝜚2 , … , 𝜚𝐽 𝜚𝑗 = 1 if 𝛽𝑗 ≠ 0 𝜚𝑗 = 0 if 𝛽𝑗 = 0 Slab: 𝛽𝜚 は 𝜌𝑗 = 1 の 𝛽𝑗 の集合 𝛽𝜚 | 𝜎𝜀 2 ~ 𝑁 𝑏𝜚 , 𝜎𝜀 2 𝛴𝜚 1 𝜎𝜀 2 ~ 𝐺𝑎𝑚𝑚𝑎 𝜈𝜀 2 , 𝑠𝜀 2 観測方程式: 𝑦𝑡 = 𝜇𝑡 + 𝛾𝑡 + 𝛽𝑇𝑥𝑡 + 𝜀𝑡 共変量： 𝛽 = 𝛽1 , 𝛽2 , … , 𝛽𝐽 MCMC で 𝛽 の事後分布を sampling するにあたり， 以下の事前分布を設定 Spike-and-slab 事前分布 spike slab (Gaussian with large variance) spike は確率質量関数なので 実際の幅はゼロ 𝛽𝑗 spike slab 𝜋𝑗 や 𝛴𝜚 などの値の設定の 仕方は論文を参照

Inference 定常事後分布が 𝑝 𝜃, 𝛼 𝑦1:𝑛 ) で表される 𝜃, 𝛼 を以下の手続きで得る (a) と (b) を交互に繰り返し 𝜃, 𝛼 のサンプリングを繰り返す 𝜃: モデルパラメータ (𝛽, 𝜎𝑦 2, 𝜎𝜇 2, 𝜎𝛿 2, ...) 𝛼: 状態ベクトル (a) Data augmentation 𝑝 𝛼 𝑦1:𝑛 , 𝜃) から 𝛼 を simulation smoother でサンプリング（Durbin et al., 2002） (b) Parameter simulation 共変量 𝛽 と 𝜎𝑦 2 以外のパラメータ 𝜃′ を 𝑝 𝜃′ 𝑦1:𝑛 , 𝛼) からサンプリング (状態 𝛼 の分散 𝜎2 は事前分布の共役性から計算が容易) その後，各共変量 𝛽𝑗 に関する 各 𝜚𝑗 および 𝜎𝑦 2 を Gibbs サンプリングで得る 1. Posterior Simulation 2. Posterior Predictive Simulation (Bayesian Model Averaging) 介入後の反実仮想の観察量 𝑝 ෤ 𝑦𝑛+1:𝑚 | 𝑦1:𝑛 , 𝑥1:𝑛 を上記の 𝜃, 𝛼 のサンプルから計算し，それらの平均値をとる （介入後の全時点を同時分布としてサンプリングしている） 異なった「共変量の組み合わせ」や「パラメータ」を平均化している => Bayesian Model Averaging による sparsity や uncertainty の推論への反映 3. Evaluating Impact 最後に各サンプリングと各時点における Causal Impact 𝜙 𝑡 (𝜏) を 𝜙 𝑡 (𝜏) = 𝑦𝑡 − ෤ 𝑦 𝑡 (𝜏) で求める

"bsts" backend CausalImpact impact_analysis.R in "CausalImpact" RunWithData RunWithBstsModel https://github.com/cran/bsts/blob/master/src/bsts.cc with data arg? yes no ConstructModel Main components of "CausalImpact" package impact_model.R in "CausalImpact" CompilePosteriorInferences impact_inferece.R in "CausalImpact" Boom::SdPrior bsts::AddLocalLevel bsts::AddSeasonal bsts::bsts bsts::AddDynamicRegression Dynamic regression? yes no ※ The details of this flow are not shown here. ComputeResponseTrajectories GetPosteriorStateSamples ComputePointPredictions ComputeCumulativePredictions "bsts" functions

03 Application to simulated data 𝑦𝑡 = 𝛽𝑡, 1 𝑧𝑡, 1 + 𝛽𝑡, 2 𝑧𝑡, 2 + 𝜇𝑡 + 𝜀𝑡 𝛽𝑡 ~ 𝑁 𝛽𝑡−1 , 0.012 𝛽0 = 1 𝑧𝑡, 1 or 2 : 90 or 360 days の周期の sin 波 𝜇𝑡 ~ 𝑁 𝜇𝑡−1 , 0.12 𝜇0 = 0 𝜀𝑡 ~ 𝑁 0, 0.12 擬似データで Causal Impact の推定結果を検証 上記の時系列を生成し，介入後 (2014 - 01 以降) は， 介入効果として， 𝑦𝑡 に 1 + 𝑒 をかけ， 𝑒 を Causal Impact で推定 ◼ 効果の大きさ 0%, 0.1%, 1%, 10%, 25%, 50%, 100% それぞれに対して，256 回の simulation を行った ◼ 効果の検知の有無は，impact の事後分布の信用区間が zero を含むか含まないかで判定 ◼ 効果が 25% 以上だと，概ね 80% 以上の確率で検知可能であった (Fig 3 - b) ◼ 介入期間の長さを変化させても，効果の 95% 信用区間に真の効果の大きさ (10%) が含まれる割合は殆ど変化しなかった (Fig 3 - c) ◼ 但し，擬似データの生成モデルが介入後の途中で変化するような場合 (例えば，動的係数 𝛽 のノイズ成分の大きさが変化する)，推定精度は大 きく悪化する (論文中の Fig. 4) ◼ 良い感じに推定できているが，使用したモデルは生成モデルと同じ表現式を使用しているので，少し leaky な評価か？ (要は，真のモデルを知った上で推定しているが，それが分からない時でも同じような推定ができるのか？) Brodersen et al., 2015. Fig 3. Adequacy of posterior uncertainty

04 Application to empirical data 実データで Causal Impact の推定結果を検証 US における広告効果のデータ (Vaver and Koehler, 2012) を使って検証 ◼ 広告は，190 地域のうち 95 地域をランダムに選択し，6 週間実施 ◼ 介入効果は広告の click 数で測定 ◼ 広告介入した 95 地域の売り上げを合計して，一つの介入地域を生成 ◼ 10,000 MCMC samples ◼ 非介入 95 地域を静的係数とした local level model を採用 (季節性は非介入地域の時系列に含まれており，それらを共変量として使用しているため，組み込まず) (a) 広告介入前をよく予測できており，介入後では，counterfactual な推 定値と観測値との乖離がおおきくなっていき，明らかに広告効果があるこ とがみてとれる (b) 広告の効果を point-wise で表しており，3 週間付近でピークを向か えた後，段々と効果が薄れている (c) 広告の累積効果は 95% BCI (bayesian credible interval) でみても有 意であり，treatment-control comparison による推定累積効果 (Vaver and Koehler, 2012) と比較してもほぼ差がない Causal Impact による推定累積効果: 22% [13%, 30%], 88,400 clicks Vaver and Koehler による推定累積効果: 21% [19%, 22%], 84,700 clicks Brodersen et al., 2015. ← Model fit Prediction →

広告の対象となっていない地域 (controls) における効果の推定結果 ◼ もし，広告の効果が介入された地域にしかないのであれば，controls における効果 はないはず ◼ 広告元の業界に関連するキーワード検索 (Google Trends) を共変量として使用 有意な効果は一切みられず 推定累積効果: 2% [-6%, 10%]

1. 『効果検証入門 正しい比較のための因果推論／計量経済学の基礎』 安井翔太 著，技術評論社，2020 2. 『時系列分析と状態空間モデルの基礎』馬場真哉 著，プレアデス出版，2018 3. 『経済・ファイナンスデータの計量時系列分析』沖本竜義 著，朝倉書店，2010 4. George et al., 1997. Approaches for Bayesian variable selection. Statist.Sinica, 7, 339–374. 5. Durbin et al., 2002. A simple and efficient simulation smoother for state space time series analysis. Biometrika. 89:603–16. 6. Abadie et al., 2010. Synthetic control methods for comparative case studies: Estimating the effect of California’s tobacco control program. J. Amer. Statist. Assoc. 105 493–505. 7. Scott et al., 2014. Predicting the present with Bayesian structural time series. International Journal of Mathematical Modeling and Optimization 5 4–23. 8. Nakajima J, West M. Bayesian Analysis of Latent Threshold Dynamic Models. Journal of Business & Economic Statistics 2013;31:151–64. https://doi.org/10.1080/07350015.2012.747847. 9. Vaver J, Koehler J. Periodic Measurement of Advertising Effectiveness Using Multiple-Test-Period Geo Experiments. Google Inc.; 2012. References