セーターの備忘録

大学編入やIT関連の内容を記していく予定です

ARIMA・SARIMAを用いてCOVID-19の感染者数を予測する

技術

はじめに

最近、新型コロナウイルスが世間を騒がせており、研究者が続々論文を発表したり、SIGNATEなどコンペサイトでも感染者数の予測コンペが行われるなど、様々な取り組みが行われています。
ここでは、時系列データ分析手法の一つであるARIMAとSARIMAを用いて日本の感染者数の予測を行いました。

ARIMAとSARIMAについて

ARIMAやSARIMAは複数の時系列データ分析手法を組み合わせたモデルになっています。
モデルの関係について、一般的に以下のような表で表されます。

f:id:resweater:20200424113436j:plain
（出典：未来を予測するビッグデータの解析手法と「SARIMAモデル」 - DeepAge）

ARIMA

ARモデル、MAモデル、ARMAモデルはいずれも定常過程を対象とした時系列モデルです。 定常過程に従う時系列データは、一定の周期で同程度の変動をしているデータのことです。 しかし、常に一定のデータは滅多になく、非定常過程のデータであることが多いため、非定常過程の時系列データに対応した手法を用いる必要があります。
ARIMAモデルは非定常過程に対応したモデルの一つです。 ARIMAモデルは前後のデータ間の差分dを定義しており、ARモデルに関する次数p、MAモデルに関連する次数qと合わせてARIMA(p, d, q)と表すことができます。

SARIMA

SARIMAはARIMAモデルに長期的な季節変動を取り入れたモデルになります。 実際の時系列データでは季節に合わせた周期性があるデータも多いため、それらを考慮したモデルとなっています。

予測で用いるデータ

COVID-19の感染者情報が載っているサイトはいくつかありますが、今回は以下のサイトからcsvファイルをインポートして用いました。

gis.jag-japan.com

読み込んだcsvファイルを加工し、1/15~4/20の日本全体の感染者数を抽出しました。

f:id:resweater:20200424115556p:plain
1/15 ~ 4/20 1日あたりの全国感染者数

コレログラムと季節成分分解

元のデータと元のデータから時間をずらしたデータとの相関のことを自己相関といい、これを表したグラフをコレログラムと言います。コレログラムを見ることで、データの周期性を調べることができます。
グラフを見た感じ、ほとんど周期性はなく、時間が経つにつれて相関性が小さくなっていることが分かります。

f:id:resweater:20200424120326p:plain
コレログラム

次に、季節成分分解をします。これにより、トレンド・季節性・残差に分解することができ、データにどのようなトレンド性があるのかなどについて知ることができます。
グラフを見ると、時間の後半にかけてトレンドの数値が大きくなっていることが分かります。

f:id:resweater:20200424121033p:plain
季節成分分解（上から順に元データ・トレンド・季節性・残差）

ARIMAでの予測

今回のデータは非定常過程のデータとみられることから、まずARIMA(p, d, q)を構築して予測してみることにしました（本当はデータが定常か非定常かについて検定をする必要がありますが、今回は省略しました、、、もっと勉強しなきゃ、、、）。
上記した通り、ARIMAにはp, d, qの3つのパラメータがあります。このうち、p, qはある程度自動推定してくれる機能があり、以下の通り実行します。

params = sm.tsa.arma_order_select_ic(diff, ic='aic', trend='nc')
print(params)

すると、以下のような結果が得られます。このうち、'aic_min_order': (4, 2)がp, qのパラメータになります。

{'aic':              0            1            2
0          NaN  1055.125479  1054.226782
1  1055.344019  1054.737810  1055.721728
2  1055.676846  1055.374114  1058.681375
3  1056.655811  1056.810369          NaN
4  1057.383730  1057.394580  1027.192160, 'aic_min_order': (4, 2)}

次に、差分dは自分で決めなければいけません。今回は差分回数は1で決め打ちします。差分データは以下の通り計算します。

diff = ts.diff()
diff = diff.dropna()

これでp, d, qの3つのパラメータが揃ったので、ARIMAモデルを構築します。tsは時系列データです。

from statsmodels.tsa.arima_model import ARIMA
arima_model = ARIMA(ts, order=(4,1,2)).fit(dist=False)

ここで、モデルから推測される値と実測値の差（残差）を表示してみます。後半になるにつれ残差が大きいことが分かります。

f:id:resweater:20200425180558p:plain
ARIMAモデルの残差

次に、構築したARIMAモデルを用いて予測をしてみます。青い線が1/15 ~ 4/20の実際の値、オレンジ色の線が4/20 ~ 5/10の予測値になります。かなり大きくぶれていて、あまり上手く予測できていなさそうなのが分かります。

f:id:resweater:20200425180905p:plain
4/20 ~ 5/10のARIMAモデルによる予測結果

SARIMAでの予測

SARIMA(p, d, q)(P, D, Q)[s]はARIMAのモデルのパラメータに加え、季節調整のsというパラメータとP, D, Qというパラメータが必要になります。
ここで、今回のデータはコレログラムでは周期性が見られませんでしたが、ドメイン知識として毎週月曜日は検査数が減るため感染者数も減るという事が分かっています。そこで、今回は7日ごと周期性があると仮定して、s=7でモデル構築を行いました。
P, D, Qは本来はグリッドサーチなどで最適な値を求める必要がありますが、今回はいずれも1で決め打ちしました。

sarima_model = sm.tsa.SARIMAX(ts, order=(4,1,2), seasonal_order=(1,1,1,7)).fit()

残差をプロットしてみると以下のようになりました。ARIMAモデルより振り幅が小さいことが分かります。

f:id:resweater:20200425181553p:plain
SARIMAモデルの残差

構築したSARIMAモデルを用いた予測結果は以下の通りです。ARIMAモデルと同じく、青い線が1/15 ~ 4/20の実際の値、オレンジ色の線が4/20 ~ 5/10の予測値になります。ARIMAモデルよりもそれっぽい予測が出来ていることが分かります。

f:id:resweater:20200425181825p:plain
4/20 ~ 5/10のSARIMAモデルによる予測結果

まとめ

時系列データに対してARIMAモデルとSARIMAモデルを用いて予測を行うことで、それっぽい結果を得ることができました。ただ、この結果はあくまで時系列データしか使っていないため、実際は様々な外的要因によって変化していくと思われます。
今回、初めて時系列データのモデル構築に挑戦したので、おかしいところがあればぜひご指摘下さい。 あと、グラフの時間軸が上手くプロットできなかったのが今後の課題です。。。

プログラム

分析で用いたプログラムはGitHubで公開しています。 github.com

参考