東京都の新型コロナウィルス時系列データを時系列分析してみる

１．背景

時系列分析の勉強をしている最中に新型コロナウィルスが発生したので、ホットな話題ということもあり、時系列分析をぶつけてみた次第です。

２．方法

・使用データは東洋経済オンラインのサイトから取得しました。
toyokeizai.net

・言語はPython。時系列分析はstatsmodelsを使用します。tsa.seasonal_decomposeによるいわゆる乗法モデルベースの基本成分分解ですね。

３．結果と考察

まずはパッケージとデータのインポートから。

import statsmodels.api as sm
import pandas as pd
import numpy as np
import requests
import io
from matplotlib import pylab as plt
%matplotlib inline
# グラフを横長にする
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6

import warnings
warnings.simplefilter('ignore')

df = pd.read_csv('prefectures.csv')

データの中身。とりあえず、目的である東京だけ見てみる。

tokyo = df[df['prefectureNameJ']=='東京都']
tokyo.tail()

	year	month	date	prefectureNameJ	prefectureNameE	testedPositive	peopleTested	discharged	deaths
4524	2020	6	15	東京都	Tokyo	5592	16351.0	4887.0	314.0
4571	2020	6	16	東京都	Tokyo	5619	16351.0	4936.0	316.0
4618	2020	6	17	東京都	Tokyo	5633	61112.0	4978.0	317.0
4665	2020	6	18	東京都	Tokyo	5674	63083.0	5017.0	317.0
4712	2020	6	19	東京都	Tokyo	5709	64844.0	5071.0	320.0

新規感染者数を分析したいのですが、累積の感染者数しかありません。
なので、手作業で新規感染者数を作りにいきます。

tokyo['NewPositive'] = 0.0
for i in range(len(tokyo)-1):
    tokyo.iloc[i+1,9] = tokyo.iloc[i+1,5] - tokyo.iloc[i,5]
plt.figure(figsize=(15, 8))
plt.plot(tokyo.NewPositive)

よくみるグラフが出てきました。ぱっと見、上昇トレンドかなと。
次に自己相関係数と偏自己相関係数を見てみましょう。

fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(tokyo.NewPositive, lags=30, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(tokyo.NewPositive, lags=30, ax=ax2)

15日前との偏自己相関が大きい。しかもマイナスの偏自己相関になりました。
ということは約2週間前の結果が小さくなると、当日の結果が大きくなる関係にあります。
人々はある日の新規感染者数を見て警戒度を変えるので上記のような関係になるということでしょうか。
すなわち、2週間前の新規感染者数が小さい→油断する→当日の新規感染者数が増える、ということなのかもしれないです。
ただし、トレンドがある場合に偏自己相関係数をそのまま解釈してよいのか、勉強不足故、よくわかりません。
トレンドがある時系列なので、周期性を観察するために１階差分をとって自己相関係数をみてみます。

diff = tokyo.NewPositive.diff()
diff = diff.dropna()
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(diff, lags=30, ax=ax1)

まあ、7日ですかね。
次に、成分分解をしてみます。

fig, (ax1,ax2,ax3) = plt.subplots(3,1, figsize=(15,8))
seasonal_decompose_res = sm.tsa.seasonal_decompose(tokyo.NewPositive, freq=7)
seasonal_decompose_res.trend.plot(ax=ax1)
seasonal_decompose_res.seasonal.plot(ax=ax2)
seasonal_decompose_res.resid.plot(ax=ax3)

上昇トレンドにみえますね。
特段何もしなければ、今後も上昇トレンドは続いていくのか、注視が必要でしょう。