時系列解析

時系列解析の手法

勉強中の時系列解析モデルです

・DID(差分の差分法)
・CausalImpact
・グレンジャー因果性検定
・ピアソン相関係数　例：https://pubmed.ncbi.nlm.nih.gov/35746459/
・ARIMAモデル
・SARIMAモデル
・GARCHモデル

DID（差分の差分法）

DIDはある地域での介入（広告配信、条例の施行など）の前後の比較を行うために用いられます。
非介入地域を比較群として扱うことで、地域特性がバイアスとして発生することがあります。
DIDの欠点は、分析のためには複数の地域や時間からのデータ取得が必要であり、その選択は担当者の仮説に依存することです。
平行トレンド仮定を満たすようなデータ取得や共変量の調整が必要です。
DIDは2時点の効果を比較するため、統合データの場合にただ差を取るだけになってしまうため、パネルデータ(同時に多数の個人、地域、事業所などを観察したデータ)のような個票があるものが良いです。

Causal Impact

Causal　Impactは介入前の期間のデータから目的変数を予測するモデルを作成します。
検索回数やアクセス数などの変数を含め、介入後のデータから介入がなかった場合の予測値を出し、その差分で効果を測定します。
ただし、CausalImpactでも平行トレンド仮定が重要であり、どのデータを使用するかは分析者の判断に委ねられます。
共変量(対照群)の時系列データはあった方が分析結果が正当化されやすいため、確実に介入の影響を受けないデータを入れる

Causal impact　結果の解釈

上段のパネル（ original ）
実線：実際の観測データ
点線：推論データ
青塗り：推論の 95% 信頼区間

中段のパネル（ pointwise ）
観測データと反実予測の違い

下段のパネル（ cumulative ）
中段のパネルで点ごとの寄与を合計した、介入の累積効果

Causal impact　レポートのデータ解釈

Actual: 実際のデータの平均と累積の値。
Prediction (s.d.): 介入がなかった場合に予測される値の平均とその標準偏差。
95% CI: 予測値の95%信頼区間。
Absolute effect (s.d.): 実際の値と予測値の差（介入の絶対的な効果）とその標準偏差。
Relative effect (s.d.): 絶対的な効果を予測値で割ったもの（介入の相対的な効果）とその標準偏差。
Posterior tail-area probability p: 介入が有意な効果を持たない場合のデータが観測される事後確率。この値が0.05より小さい場合、介入が有意な効果を持つと結論されることが多い。
Posterior prob. of a causal effect: 介入が何らかの効果を持つと結論される事後確率。

モデルの評価

ACF, PACF：後述

備考

・Causal impactの結果は毎回違う
これはベイズ構造時系列モデリングを行うBSTS(Bayesian Structural Time Series)が使うMCMC(マルコフ連鎖モンテカルロ法)による推定のためです
set.seedによるランダム性の排除は無効です※1

・シーズナリティの設定は一通りならば可能
nseasons. 季節成分の期間。季節成分を含めるには、1より大きい整数を設定します。たとえば、データが日次の観測を表す場合、曜日成分には7を使用します。このインターフェースは現在、1つの季節成分のみをサポートしています。複数の季節成分を指定するには、bstsを使用してモデルを直接指定し、bsts.modelとしてフィットしたモデルを渡します。デフォルトは1で、これは季節成分が使用されないことを意味します。

season.duration. 各シーズンの期間、つまり、各シーズンが範囲とするデータポイントの数。たとえば、日次の粒度のデータに曜日成分を追加するには、引数 model.args = list(nseasons = 7, season.duration = 1) を加えます。また、時刻ごとの粒度のデータに曜日成分を追加するには、model.args = list(nseasons = 7, season.duration = 24) を使用します。デフォルトは1です。

例えば以下は、週次データに対して1年ごとの周期性を持たせたものです
週ごとのデータそのまま（＝１）に対して1年ごと(=52)の周期性を持たせます

impact <- CausalImpact(zoo_data, pre_period, post_cov_period,
model.args = list(nseasons = 52, season.duration = 1))

※1　https://stats.stackexchange.com/questions/284204/q-r-causalimpact-order-of-control-time-series

グレンジャー因果性検定

・予測的因果関係を示すもの、ただしこの検定で示されれるのは名前の通りの因果ではなく相関関係という事に注意が必要
分かりやすい解説

実例

関連を見たい変数を選択(今回はVN, vac, covの3つ)
↓
VARモデル(ベクトル自己回帰モデル：p期前までの結果で未来を予測する)を作成(p＝13、今回は13週分のデータ)
↓
causarity関数でグレンジャー因果性テストを行う

実例コード

data_var <- df_weekly11[, c("n_VNpatients", "n_vacpatients", "n_covpatients")]
var_model <- VAR(data_var, p=13, type="both")
test_result_vac_to_VN <- causality(var_model, cause="n_vacpatients")
print(test_result_vac_to_VN)

※VARモデル（Vector Autoregression Model）の設定時に、type引数としてconstを指定すると、モデルに定数項が含まれることを意味します。VARモデルを推定する際のその他のtypeオプションには以下のようなものがあります：

この中でモデルに合ったものを選びます
0を起点に始まるものでなければ定数項を入れ、トレンドは大きな流れがある場合は入れておくべきです

const: 定数項のみを含む
trend: 定数項とトレンドを含む
both: 定数項とトレンドの両方を含む
none: 定数項もトレンドも含まない

※ここでグレンジャーの因果性の検定に使用するのは$Grander の部分のF-Test の値である($Instant の部分は無視してよい。これはグレンジャーの瞬時的因果性というものの検定でここでは関係ない)

※causality()関数によるグレンジャー因果性検定はcauseで指定した変数が他の残りの変数に対してグレンジャー因果性を持っているかどうかを同時に検定するので、例えば3変数VARのような場合ある1つの変数が他のもう1つの変数のみへグレンジャー因果性がないという帰無仮説の検定には使えない(causality()関数は他の2つの変数へグレンジャー因果性がないという帰無仮説を同時に検定している)
つまり2つの原因の影響を同時に推定する事が出来ない

ピアソン相関係数

・-1~+1の相関係数で表される

ARモデル

AR（自己回帰）モデルは、一定期間に遡る「過去の自分のデータ」を説明変数とします。そして、過去のデータに対し回帰を行うことで、現在の値を予測します。
そのため自己相関があるモデルにしか適応が出来ません。

p次のARモデルをAR(p)と表す。

AR(1)では1世代前のデータのみを用いて推測します。
AR(p)ではp世代前までのデータを用いて推測します。

MAモデル

MA（移動平均）モデルでは、ある時刻のデータを「過去の時刻での誤差項」を用いて表現します。別の表現を使うと、過去の誤差を考慮しながらデータの推移を予測します。
つまり、誤差の集積です。

q次のMAモデルをMA(q)と表す。

MA(1)では1世代前の誤差のみを用いて推測します。(日次データであれば前日、週次なら前週)
MA(q)では1-q世代の誤差を用いて推測します。

ARMAモデル

ARモデルとMAモデルの組み合わせです
前述のARのp，MAのqを用いてARMA（p，q）と表されます

ARモデルは非定常でも対応可能ですが、こちらは定常状態でのみ適応可能なモデルとなっています

ARIMAモデル

1次和分過程(1次Iモデル)とも呼ばれ、その差分系列が(p,q)次ARMAモデルで表されるとき、単位根過程は(p,1,q)次ARIMAモデル(自己回帰和分移動平均)と呼ばれます。一般に、d階差分を取った系列が(p,q)次ARMA過程に従う過程を(p,d,q)次ARIMAモデルと呼びます。
d階差分を取ることでトレンドを除去し、定常状態にすることでARMAモデルを適応可能としているイメージです

データの観察
↓
予測モデルの選択
↓
予測モデルの評価

このフローに沿ってモデル化します
詳しくはこちらに書いてあります
https://www.i-juse.co.jp/statistics/jirei/sympo/10/arima-model.html

SARIMAモデル

・季節変動があるデータに対して有効
・時系列方向にARIMAモデルを使い、かつ周期方向にもARIMAモデルを使っているモデル

GARCHモデル

・ARIMAモデルにボラティリティクラスタリングの要素を盛り込んだモデル

注意点

DID、Causal impactの手法は、実験上の制約やABテストが不可能な場合に有用
しかし、複数の介入が同時に行われた場合（例：コロナとコロナワクチン）、どの介入が効果をもたらしたのか判断が難しい
共変量に価格などを含めることで一部対処可能だが、複数の介入が同時にある場合、効果推定は困難となる可能性
グレンジャー因果性検定はあくまで予測出来るかを検定するもので、因果関係を証明するものではない
ピアソン相関係数は抽出する期間によっては偶然の相関をとらえる可能性があり、また因果関係を証明するものではない
ARIMAモデルは季節変動がある非定常データでは予測性が落ちる
ARIMAモデルではボラティリティを不変として仮定しているため、ボラティリティクラスタリングがるような事象予測に弱い
標準化残差のプロット、ヒストグラム、Q-Qプロット、自己相関係数などで残差の周期性が無い事、正規分布している事を確認し、予測モデルを評価します
次のようなモデルで評価することは一般的です
ARMA (自己回帰移動平均) モデル: ARMAモデルは自己相関と移動平均を考慮に入れて時系列データをモデル化します。残差分析は、モデルのパラメーターが適切であるかどうかを確認するのに役立ちます。
ARIMA (自己回帰和分移動平均) モデル: ARIMAモデルは、非定常時系列データをモデル化するために使用されます。残差分析は、モデルがデータを適切に捉えているかどうかを評価するために重要です。
SARIMA (季節性自己回帰和分移動平均) モデル: SARIMAモデルは季節性コンポーネントを持つ時系列データをモデル化します。残差分析は、季節性パターンが適切にモデル化されているかどうかを評価するために必要です。
状態空間モデルやカルマンフィルター: 状態空間モデルやカルマンフィルターを使用する場合も、残差の非周期性や正規性をチェックすることが重要です。
その他の時系列モデル (例えば、VAR、GARCH、STLなど): これらのモデルでも、残差分析はモデルの適切さやパラメーターの適切さを評価するために必要です
Causal impact：残差分析の中でもACF、PACFについて分析し、自己相関がないことを確認スべきとされています

時系列解析を始める前に

始める前の事を一番最後に書くなよって話なんですが
時系列解析の前に元データを適切に処理する事が予測性に寄与します
以下のようなものです

季節性の除去: Rでの主な方法は、decompose()やstl()関数を使用することです。これらの関数は、時系列データを季節性、トレンド、そして残差に分解します。

トレンドの除去: トレンドを取り除くことで、短期的な変動や周期性をより鮮明に観察することができます。

変動の安定化: データの変動が時間とともに増加する場合（例：指数関数的成長）、ログ変換などの変換を使用して変動を安定化することができます。

ラグ効果の考慮: 特定のイベントや介入の影響が即座に現れるとは限らないため、遅延効果を考慮することが重要です。

外れ値の検出と処理: 時系列データには外れ値が含まれることがよくあります。外れ値は、モデルの適合度を低下させる可能性があるため、検出して適切に処理することが重要です。

シーズナリティの振幅: seasonal成分の最大値と最小値の差を計算することで、振幅の大きさを数値的に把握できます。
amplitude <- max(decomposed_data$seasonal) - min(decomposed_data$seasonal) print(amplitude)
月別のシーズナリティの影響: 月ごとの平均的なseasonal成分の値を計算することで、どの月が最も高い影響を受けているかや、その影響の大きさを確認できます。データが週次の場合、特定の週が季節性のピークや谷を持っているかを確認することもできます。
#データの開始日を入れる
start_date <- as.Date("2010-01-03")
#データの日数分の日付を生成
dates <- seq(start_date, by = "week", length.out =　length(decomposed_data$seasonal))
#月別の平均的なシーズナリティの影響を計算
monthly_effect <- tapply(decomposed_data$seasonal, months(dates), mean, na.rm = TRUE) print(monthly_effect)
プロット: seasonal成分をプロットして、繰り返しのパターンやピーク/谷のタイミングを視覚的に確認することができます。ピークや谷の時期を特定することで、シーズナリティが最も影響を及ぼす時期を明らかにすることができます。
plot(decomposed_data$seasonal, type = "l", main = "Seasonal Component")