Counterfactualの推論
Counterfactual（もしそれが発生していなければ、どうなっていたか）の概念についてPythonのライブラリを使った検証例について熟考していたところ、同様のCounterfactualに関するのRパッケージがCRANに提供されているため、このライブラリを評価してみました。それは"CausalImpact"というベイズ構造化モデルを使った因果推論のRパッケージです。
このパッケージはグーグル社のKay H. BrodersenよりApacheライセンスで提供されています。内部でbstsライブラリ(Bayesian structual time series)、Boomライブラリ(Bayesian Object Oriented Model)を使っています。
時系列上の現象で作られた介入または干渉(intervention)[ 注 ]について因果の影響を予測するのは難しいことです。パッケージの目的は、interventionの後で、もしその介入または干渉（intervention）が発生していなければ、どのように反応が進んだかを予測するために、構造ベイズ時系列モデルを使うことで、この困難さにとり組むことです。
注：干渉または介入（intervention）
ここでいうinterventionとは、Judea Pearlのdo-operaterによる因果推論の intervention を指しています。
因果加法モデル(causal additive model)では、因果モデルの主要なツールとして、Pearlによってdo-calculus が導入されました。通常の条件確率 P[Y|X] はXの観測値 x とリンクしますが、do(.)条件確率 P[ Y| do(X=x) ]は強制してXに値 x を与えます。二つの条件確率はlookingとdoing の対立概念です。
参考文献
"Inferring causal impact using Bayesian structural time-series models"
Kay H. Brodersen Fabian Gallusser Jim Koehler Nicolas Remy Steven L. Scott
"Causality Models,Reasoning, and Inference" Judea Pearl
インストール
CausalImpactをCRANからインストールします。Boomを要求するのでBoomもインストールします。
> install.packages("Boom") --- このセッションで使うために、CRAN のミラーサイトを選んでください ---
URL 'https://mirrors.nics.utk.edu/cran/bin/macosx/big-sur-arm64/contrib/4.2/Boom_0.9.11.tgz' を試しています
Content type 'application/x-gzip' length 87420992 bytes (83.4 MB)
==================================================
downloaded 83.4 MB
ダウンロードされたパッケージは、以下にあります
/var/folders/z3/y2vb_4653kjcytsygg5bkv7w0000gn/T//Rtmpb5olu1/downloaded_packages > install.packages(CausalImpact) 要求されたパッケージ bsts をロード中です
要求されたパッケージ BoomSpikeSlab をロード中です
要求されたパッケージ Boom をロード中です
次のパッケージを付け加えます: ‘Boom’
以下のオブジェクトは ‘package:stats’ からマスクされています:
rWishart
次のパッケージを付け加えます: ‘BoomSpikeSlab’
以下のオブジェクトは ‘package:stats’ からマスクされています:
knots
要求されたパッケージ zoo をロード中です
次のパッケージを付け加えます: ‘zoo’
以下のオブジェクトは ‘package:base’ からマスクされています:
as.Date, as.Date.numeric
要求されたパッケージ xts をロード中です
次のパッケージを付け加えます: ‘bsts’
以下のオブジェクトは ‘package:BoomSpikeSlab’ からマスクされています:
SuggestBurnサンプルデータ
期間2021年5月12日〜2023年3月17日のドル円為替レートを例に評価してみます。データポイントは東欧の侵攻のタイムスパンに属しています。為替の時系列データはcsvファイルで読み込むことを想定しています。
> library(CausalImpact)
> usdfxr <- read.csv(file="./Downloads/xxx/USDJPY.csv",encoding="UTF-8", header=F, sep=",")
> tail(usdfxr)
V1 V2 V3 V4 V5 V6
432 2023-03-10 136.1675 136.99 134.119 134.86 0
433 2023-03-13 134.4365 135.037 132.298 133.211 0
434 2023-03-14 133.209 134.898 133.039 134.242 0
435 2023-03-15 134.242 135.112 132.2235 133.344 0
436 2023-03-16 133.344 133.826 131.729 133.483 0
437 2023-03-17 133.482 133.58 131.577 131.846 0ヘッダを定義します。
> names(usdfxr) <-c("date","open","high","low","close","volume")
> head(usdfxr)
date open high low close volume
1 date open high low close volume
2 2021-05-12 108.613 109.209 108.611 109.032 0
3 2021-05-13 109.669 109.785 109.482 109.5435 0
4 2021-05-14 109.4545 109.6535 109.21 109.348 0
5 2021-05-17 109.281 109.502 109.076 109.218 0
6 2021-05-18 109.21 109.285 108.839 108.946 0
> usd_close <- c(as.numeric(as.character(usdfxr[-1,]$close)),NA)
> head(usd_close)
[1] 109.0320 109.5435 109.3480 109.2180 108.9460 108.8980CausalImpactのトレーニングで用いるデータの終端の詳細を規定します。モデルのトレーニング用(Pre-intervention period)、Counterfactual予測(post-intervention period)
> pre.period <- c(1, 195)
> post.period <- c(196, 430)推論を実行するために、CausalImpactの分析を走らせます。
> impact <-CausalImpact(usd_close, pre.period,post.period)> plot(impact)
デフォルトでは三つのパネルを表示します。最初の表示が、事後の処置用のCounterfactual予測。2番目のパネルが観測データとCounterfactual予測の相違です。3番目のパネルは2番目のパネルからinterventionの蓄積の影響を図示しています。
分析結果
分析の数値結果を得るには
> summary(impact)Posterior inference {CausalImpact}
Average Cumulative
Actual 134 31594
Prediction (s.d.) 115 (2.7) 27103 (628.2)
95% CI [110, 120] [25893, 28285]
Absolute effect (s.d.) 19 (2.7) 4491 (628.2)
95% CI [14, 24] [3309, 5701]
Relative effect (s.d.) 17% (2.7%) 17% (2.7%)
95% CI [12%, 22%] [12%, 22%]
Posterior tail-area probability p: 0.00113
Posterior prob. of a causal effect: 99.88713%
For more details, type: summary(impact, "report")上記CausalImpactの事後推論の結果に表示されているActualは実際のデータ、Absolute effectは実際のデータと予測値の乖離、Relative effectはデータと予測値の乖離の％表示を意味しています。Averageは平均、Cumulativeは累積値です。
95% CIは 95%のConfidence Intervalの略です。確率分布の中心から離れた両端の極端な値（上位のテイルと下位のテイル部分）を外した部分を指しており、信頼度が95%の確率でこの範囲に入ることを意味していいます。
推論の結果は、予測値が平均115円で95%の確率で110円から120円の間にあることを示し、実際のデータ(134円)が平均で 17％乖離していることを示しています。
posterior tail-area probability p(事後テイル領域確率): は、この結果が偶然に発生する確率を意味し、その値は非常に小さくなっています(p=0.00113%)。 これは乖離の原因としては統計的に非常に高い因果性があることを意味します。
Posterior prob. of causal effect:がその確率を示しており、99.88713% の確率で因果性があることを示しています。
備考
CausalImpactではトレーニングに使うモデルは、ベイズ構造化時系列モデルを用いています。パラメータの設定は、データセットをトレーニング用と予測用に分離して設定するだけなので、簡易にCausal effectに関するベイズ推論を実行することができます。ただし、内部で使用するモデルはベイズ構造化時系列モデルに限定されます。上記の例では非常によく適合しています。
季節要因やその他固有の変数を追加し、詳細なモデルを特別に作り、Counterfactualの予測をするにはコードの量もMCMCの所要時間も要します。階層ベイズでは数日間MCMCを走らせることがあります。
状態空間モデルで適合できる時系列データでは、このライブラリを使うことで、所要時間もかからず簡単に因果性の検証ができるようです。
サンプルデータとして用いた為替レートの推移については、以下にも記述。
