COVID-19のデータから感染者数のグレンジャー因果性を調べる

はじめに

世界を危機に陥れてるCOVID-19ですが、
昨今勉強を進めている時系列分析で何か知見が得られないか調べてみることにしました。
しかしながら、初学者であるため間違った解釈をしている可能性があります。
どしどし、ご指摘お願いいたしますm(_ _)m

やりたいこと

VARモデルを使って中国の感染者数から韓国、日本への感染者数に対してグレンジャー因果性があるかを検定する。検定には、Pythonのstatsmodelsを使いました。
www.statsmodels.org

VARモデル

ARモデル（自己回帰モデル）を多変量に拡張したものであり、以下の式で表されます。

$\begin{align}\begin{aligned}Y_t = \nu + A_1 Y_{t-1} + \ldots + A_p Y_{t-p} + u_t\\u_t \sim {\sf Normal}(0, \Sigma)\end{aligned}\end{align}$ 　　・・・1

$Y_t$ ：t時点の目的変数ベクトル

$Y_{t-i}$ ：t-i時点の目的変数ベクトル

$\nu$ ：n×1定数ベクトル

$A_p$ ：n×n定数ベクトル

$u_t$ ：ベクトルホワイトノイズ

Varモデルの優れている点は、自身の過去のデータだけでなく関係ありそうな他変数の過去のデータも回帰に使えることです。例えば、仮に中国、韓国、日本の感染者数にそれぞれ相関があると仮定します。
その場合、1式は以下のように書き下せます。

3変量のVAR(p)モデル

t時点の中国の感染者数=

$\nu_A + A_{11}×(t-1時点の中国の感染者数)+A_{12}×(t-1時点の韓国の感染者$
$数)+A_{13}×(t-1時点の日本の感染者数)+$
$A_{21}×(t-2時点の中国の感染者数)+A_{22}×(t-2時点の韓国の感染者+$
$数)+A_{23}×(t-2時点の日本の感染者数)+・・・+$
$A_{p1}×(t-p時点の中国の感染者数)+A_{p2}×(t-p時点の韓国の感染者数)+$
$A_{p3}×(t-p時点の日本の感染者数)+$
$u_{ta}$
・・・2-A

t時点の韓国の感染者数=

$\nu_B + B_{11}×(t-1時点の中国の感染者数)+B_{12}×(t-1時点の韓国の感染者$
$数)+B_{13}×(t-1時点の日本の感染者数)+$
$B_{21}×(t-2時点の中国の感染者数)+B_{22}×(t-2時点の韓国の感染者+$
$数)+B_{23}×(t-2時点の日本の感染者数)+・・・+$
$B_{p1}×(t-p時点の中国の感染者数)+B_{p2}×(t-p時点の韓国の感染者数)+$
$B_{p3}×(t-p時点の日本の感染者数)+$
$u_{tb}$
・・・2-B

t時点の日本の感染者数=

$\nu_C +C_{11}×(t-1時点の中国の感染者数)+C_{12}×(t-1時点の韓国の感染者$
$数)+C_{13}×(t-1時点の日本の感染者数)+$
$C_{21}×(t-2時点の中国の感染者数)+C_{22}×(t-2時点の韓国の感染者+$
$数)+C_{23}×(t-2時点の日本の感染者数)+・・・+$
$C_{p1}×(t-p時点の中国の感染者数)+C_{p2}×(t-p時点の韓国の感染者数)+$
$C_{p3}×(t-p時点の日本の感染者数)+$
$u_{tc}$
・・・2-C

沖本本pp.74によると

「モデルに含める変数は分析者の興味によって決めるものであり、目的に応じて様々です。例えば、国際株式市場の関係を分析したい場合は、日本、アメリカ、イギリスの株式収益率などからなるベクトルを考えます。（中略）重要なことは、分析者が、動学的関係に興味のある変数を全てモデルに含めることであるが、（中略）」

と書かれており、どんな変数をベクトルに加えるかは分析者の主観によるところもあるようです。

また、Varモデルで気を付けないといけないことは、非定常なデータには適用できないということです。例えば、Pythonで時系列分析を扱う際によく用いられるstatesmodelsのVarモデルのページには以下の注意書きがあります。非定常過程やトレンドを持つ過程は定常過程に変換してからVarモデルを適用しないと、不適切な結果を招くと書かれています。

f:id:banquet-kuma:20200405015818p:plain — Varモデルの注意書き

www.statsmodels.org

Varモデルの解説はこちらのサイトが分かり易かったです。
tjo.hatenablog.com

グレンジャー因果性

背後に明確な理論（経済・ファイナンス理論など）を仮定せずとも時系列データだけから因果性の有無を判断できる指標です。

沖本本pp.80には以下のように定義されています。

定義4.1(グレンジャー因果性)
現在と過去のxの値だけに基づいた将来のxの予測と、現在と過去のxとyの値に基づいた将来のxの予測をして、後者のMSE(Mean Squared Error ;平均2乗誤差)の方が小さくなる場合、 $y_t$ から $x_t$ へのグレンジャー因果性(Granger causality)が存在すると言われる。

定義4.2（一般的なグレンジャー因果性）
$x_t$ と $y_t$ をベクトル過程とする。また、xとyの現在と過去の値を含む、時点tにおいて利用可能な情報の集合を $Ω_t$ とし、 $Ω_t$ から現在と過去のyを取り除いたものを $Ω̃ _t$ とする。このとき、 $Ω̃ _t$ に基づいた将来のxの予測と、 $Ω_t$ に基づいた将来のxの予測を比較して、後者のMSEのほうが小さくなる場合、 $y_t$ から $x_t$ へのグレンジャー因果性が存在するといわれる。ここで、MSEの大小は行列の意味での大小であることに注意されたい。

例えば、日本の感染者数から中国の感染者数へのグレンジャー因果性が存在しないということは、2-A式に於いて $A_{13}=A_{23}=・・・=A_{p3}=0$ と同値になります。つまり、中国の感染者数を表す式に於いて日本の感染者数に関する係数が全て0になるということです。これは、以下の手順で判断できます。

沖本本pp.81

n変量VAR(p)モデルにおけるグレンジャー因果性検定の手順
(1)VARモデルにおける $Y_{kt}$ のモデルをOLSで推定し、その残渣平方和を $SSR_1$ とする。
(2)VARモデルにおける制約を課したモデルをOLSで推定し、その残渣平方和を $SSR_0$ とする。
(3)F統計量を $F=\frac{(SSR_0-SSR_1)/r}{SSR_1/(T-np-1)}$ で計算します。rはグレンジャー因果検定に必要な制約の数である。
(4) $rF$ を $χ^2(r)$ の95%点と比較し、 $rF$ に方が大きくければ、ある変数(群)から $Y_{kt}$ へのグレンジャー因果性は存在し、小さければグレンジャー因果性は存在しないと結論する。

グレンジャー因果性検定に於いて、
帰無仮説( $H_{0}$ )は「各変数群がその他の変数へ対してグレンジャー因果性がない」
になります。つまり、帰無仮説が棄却されればグレンジャー因果性があると判断されるわけです。

(2)の「制約を課した」という意味は、
2-A式に於いて $A_{12}=A_{22}=・・・=A_{p2}=0$ や $A_{13}=A_{23}=・・・=A_{p3}=0$ を課すという意味です。つまり、t時点の中国の感染者数をt-1〜t-p時点の中国の感染者数だけで表すと言うことです。日本、韓国の過去の感染者数を考慮したモデルでの残渣 $SSR_1$ に対して、制約を課したモデル（中国の過去の感染者数だけを考慮したモデル）での残渣 $SSR_0$ が十分に大きければ中国の感染者数と日本、韓国の感染者数に因果性があると判断するわけです。

グレンジャー因果性には長短・短所があり、まとめると以下のようになります。

長所
●定義が非常に明確である
●時系列データを用いて容易に検定できる

短所
●因果性が存在する必要条件であるが、十分条件でない
●真の因果性の方向とグレンジャー因果性の方向が一致するとは限らない
●定性的概念であり、関係の強さが測れない

因果関係の方向がはっきりしている場合を除いては、ある変数が他の変数の予測に有用かどうかと言う観点で解釈するのが良いみたいです。

データを確認する

各国の感染者数のデータは以下のサイトからcsvで入手しました。

data.humdata.org

1/22〜3/21に於ける、中国、韓国、日本、米国、イタリア、ドイツ、スペイン累積感染者数を可視化します。

# 必要なライブラリのインポート
import statsmodels as sm
import statsmodels.api as sap
from statsmodels.graphics import tsaplots
from statsmodels.tsa import stattools
from statsmodels.tsa.api import VAR, DynamicVAR
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
sns.set()

# データの読み込み
Confirmed_df = pd.read_csv("time_series-ncov-Confirmed0122-0321.csv",\
header=1,index_col="#date", parse_dates=True)

##affected+infected+value+num(感染者数)をint型にキャストしておく
Confirmed_df["#affected+infected+value+num"] = Confirmed_df["#affected+infected+value+num"].astype(int)

# 国ごとにデータを取得する
# 日付け（昇順）でソートしておく
China_pre_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "China"].sort_index()
US_pre_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "US"].sort_index()
Italy_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Italy"].sort_index()
Spain_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Spain"].sort_index()
Germany_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Germany"].sort_index()
Korea_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Korea, South"].sort_index()
Japan_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Japan"].sort_index()

# ChinaとUSは複数の州のデータがあるので国全体の合計にする
China_cumsum = China_pre_cumsum.groupby(['#date'])['#affected+infected+value+num'].sum()
US_cumsum = US_pre_cumsum.groupby(['#date'])['#affected+infected+value+num'].sum()

# 原系列の表示(感染者数の累積和)
plt.figure(figsize = (20,10))
ax = plt.subplot()

ax.plot(China_cumsum,label = "Total Confirmed in China")
ax.plot(US_cumsum,label = "Total Confirmed in US")
ax.plot(Italy_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Italy")
ax.plot(Spain_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Spain")
ax.plot(Germany_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Germany")
ax.plot(Korea_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Korea")
ax.plot(Japan_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Japan")


plt.title('Coronavirus COVID-19',fontsize = 20)
plt.xlabel('date',fontsize = 20)
plt.ylabel('Total number of confirmed people',fontsize = 20)
plt.tick_params(labelsize=20)
plt.setp(ax.get_xticklabels(), rotation=90)
plt.legend(loc = "upper left",fontsize = 20)

f:id:banquet-kuma:20200405113852p:plain — 1/22 〜 3/21の各国の累積感染者数

今回使用したデータでは、中国と米国は州別にデータが入力されていたため、国全体の合計に変換しました。3/22時点では、中国の累積感染者数が最も多いですね。
日本、韓国だけを抜き出したのが以下のグラフです。
感染者数は増加傾向にあるものの、米国、イタリア、ドイツ、スペインに比べると明らかに傾きが緩いです。

f:id:banquet-kuma:20200405114325p:plain — 累積感染者数から日本、韓国のみ抜粋

次に、上記のデータに対して1階差分を取り日毎の感染者数を算出します。
1/22は前日のデータがないため、データが消失します。

# 日毎の感染者数を求める(1階差分を取る)
# 日毎の感染者数が負になっている日はデータの誤りとして除いている
China_diff = China_cumsum.diff().dropna()
China_diff = China_diff[China_diff>=0]
US_diff = US_cumsum.diff().dropna()
US_diff = US_diff[US_diff>=0]
Italy_diff = Italy_cumsum["#affected+infected+value+num"].diff().dropna()
Italy_diff = Italy_diff[Italy_diff>=0]
Spain_diff = Spain_cumsum["#affected+infected+value+num"].diff().dropna()
Spain_diff = Spain_diff[Spain_diff>=0]
Germany_diff = Germany_cumsum["#affected+infected+value+num"].diff().dropna()
Germany_diff = Germany_diff[Germany_diff>=0]
Korea_diff = Korea_cumsum["#affected+infected+value+num"].diff().dropna()
Korea_diff = Korea_diff[Korea_diff >=0]
Japan_diff = Japan_cumsum["#affected+infected+value+num"].diff().dropna()
Japan_diff  = Japan_diff[Japan_diff>=0]

#グラフの表示は累積感染者数と同じなので割愛

f:id:banquet-kuma:20200405122507p:plain — 1/22 〜 3/21の日毎の感染者数

中国のデータに於いて、2/13に巨大なピークがあります。これは湖北省が感染者の認定基準を変えたことによるものと思われます。こういう場合データをどう扱うか悩みますが（汗）、一旦このまま分析を進めることにします。

f:id:banquet-kuma:20200405124342p:plain — 中国湖北省に於ける認定基準の変更

ja.m.wikipedia.org

日毎の感染者数で韓国、日本だけを抜き出したのが以下のグラフ。
韓国はピークを過ぎているが、日本は増加トレンドにあるように見えます。

f:id:banquet-kuma:20200405181311p:plain — 日毎の感染者数から日本、韓国のみ抜粋

単位根検定を行う

前述したとおり、Varモデルは定常過程を前提としているため、過程が定常か調べておく必要があります。
これは単位根検定で調べることができます。単位根検定には主に以下の3つがあるようです。原理は理解できていないので割愛しますが、対立仮説と帰無仮説の対応は以下のとおりです。

名称	帰無仮説( $H_{0}$ )	対立仮説( $H_{1}$ )
Dickey-Fuller(DF)検定	単位根AR(1)過程である	定常AR(1)過程である
拡張DF(ADF)検定	単位根AR(p)過程である	定常AR(p)過程である
Phillips-Perron(PP)検定*	単位根AR(p)過程である	定常AR(p)過程である

（*ホワイトノイズの分散不均一性まで考慮した検定手法）

今回は日毎の中国、韓国、日本のデータに対してADF検定を行います。
結果は帰無仮説と対立仮説のモデルの置き方によって違ってくるため、3つのモデルを考えます。

-	帰無仮説( $H_{0}$ )	対立仮説( $H_{1}$ )	説明	regression
場合1	$𝑦_𝑡=𝑦_{𝑡−1}+𝑢_𝑡$	$𝑦_𝑡=𝜌𝑦_{𝑡−1}+𝑢_𝑡$ , \|𝜌\|<1	トレンド項なし、定数項なし	nc
場合2	$𝑦_𝑡=𝑦_{𝑡−1}+𝑢_𝑡$	$𝑦_𝑡=c+𝜌𝑦_{𝑡−1}+𝑢_𝑡$ , \|𝜌\|<1	トレンド項なし、定数項あり	c
場合3	$𝑦_𝑡=c+𝑦_{𝑡−1}+𝑢_𝑡$	$𝑦_𝑡=c+𝜌𝑦_{𝑡−1}+𝛿𝑡+𝑢_𝑡$ , \|𝜌\|<1	1次のトレンド項あり、定数項あり	ct

表中のregressionは以下で行うstattools.adfullerのオプションのregressionに対応します。
𝛿𝑡は1次のトレンドを表ます。

# 中国の日毎の感染者数に対してADF検定を行う
# 出力値は先頭からt値、p値、使用されたラグの数、計算に使用されたデータ数、各棄却域に対するt値、AICの最大値
# トレンド項あり（2次まで）、定数項あり
ctt = stattools.adfuller(China_diff, regression="ctt")
# トレンド項あり（1次）、定数項あり
ct = stattools.adfuller(China_diff, regression="ct")
# トレンド項なし、定数項あり
c = stattools.adfuller(China_diff, regression="c")
# トレンド項なし、定数項なし
nc = stattools.adfuller(China_diff, regression="nc")

# 韓国、日本も同様に行う

日毎の感染者数に対する単位根検定のp値をまとめると以下のようになりました。

国	場合1	場合2	場合3	結論
中国	0.051	1.0e-4	4.4e-05	場合2,3を仮定すると定常過程と言える
韓国	0.26	0.55	0.70	いずれの場合も単位根過程でないとは言えない
日本	0.74	0.80	2.7e-7	場合3を仮定すると定常過程と言える

韓国のデータは定常過程と言えないので除外し、中国の感染者数と日本の感染者数のグレンジャー因果性だけを検定することにします。

Varモデルの構築

中国と日本の日毎の感染者数のデータでからVarモデルを構築します。

# 中国と日本のデータで共通する日付けのデータだけを抜き出す
diff_CJ = pd.merge(China_diff, Japan_diff, on='#date', how='inner')
# 列名を変更する
diff_CJ.columns = ["China_diff","Japan_diff"]

#Varモデルを作成する
model_diff = VAR(diff_CJ)

フィッティングにあたってモデルのラグ（1式のp）を決めてやる必要があります。
現時点から何時点過去に遡って回帰式を構成するかを決める作業です。
これは、情報量基準(AIC、BIC）を指標に決めます。

# 最適なラグを探索する。ラグ10まで
model_diff.select_order(10).summary()

f:id:banquet-kuma:20200405210124p:plain:w300 — 各ラグに於ける情報量基準の計算

ラグ1まで考慮した場合に情報量基準が最小になりました。
実際は以下のようにmaxlags=10とic="aic"と指定してやることでラグ10までの範囲でAICが最小となるラグを自動的選んでフィッティングしてくれます。

# AIC基準で最適なラグを選択したモデルへのあてはめ
result_diff = model_diff.fit(maxlags=10, ic="aic")

#あてはめの結果を表示
result_diff.summary()

結果のサマリーを表示します。

f:id:banquet-kuma:20200406012134p:plain — あてはめの結果

サマリーの中に記載されている係数を用いて、
中国の感染者数と日本の感染者数は以下のように書けます。

t時点の中国の感染者数=

$const+ L1.China\_diff×(t-1時点の中国の感染者数)+$
$L1.Japan\_diff×(t-1時点の日本の感染者数)$

t時点の日本の感染者数=

$const+ L1.China\_diff×(t-1時点の中国の感染者数)+$
$L1.Japan\_diff×(t-1時点の日本の感染者数)$

それぞれの係数 $L1.China\_diff$ 、 $L1.Japan\_diff$ は、
中国の感染者数(China_diff)については、
「Results for equation China_diff」に記載されている値を
日本の感染者数(Japan_diff)については、
「Results for equation Japan_diff」に記載されている値を使います。

黒枠内に示した回帰係数のp値を見ると、
日本の感染者数に於いては、 $L1.China\_diff$ が有意になっておらず
（つまり、帰無仮説 $L1.China\_diff$ =0を棄却できない）、日本の感染者数に及ぼす
中国の感染者数の影響が有意でないことが示唆されました。

ここまでの結果で、中国の感染者数から日本の感染者数に対してグレンジャー因果性がないことが示唆されますが、一応グレンジャー因果性検定をやってみます。

result_diff.test_causality("Japan_diff", "China_diff", kind = "f").summary()

f:id:banquet-kuma:20200405215817p:plain:w300 — グレンジャー因果性検定の結果

帰無仮説：
「中国の感染者数から日本の感染者数に対してグレンジャー因果性が存在しない」を棄却できませんでした。したがって、昨日の中国の感染者数と今日の日本の感染者数にグレンジャー因果性がないと結論付けました。

この結果は以下を考えると妥当という気がします（分析手法が正しければですが、、）。
●1日という短期間で中国 → 日本間でウイルスの媒介が波及しそうにないこと
●感染が発覚した人のみのデータを使っているので相関が上手く拾えていない可能性があること

最後に

COVID-19の実データを使ってグレンジャー因果性検定をやってみましたが、かなり難しさを感じました。テキストに付随するデータだとすんなり結論を導けるのですが、実データとなると「そもそもこのモデルを使って良いデータなのか？前提条件を満たしているのか？」と言った疑念が浮かんできます（笑）。実務やコンペなどで経験を積んで正しく分析できるようになりたいです。