自己回帰（AR）モデルを東京都のCOVID-19陽性者推移に当てはめてみた

#はじめに
最近本業で自己回帰モデルを勉強してるのですが、感染者数の時系列はちょうど題材として良いのではないかと思い、試してみました。

#使用するデータ
東京都が公開している、2021/4/30現在の新規陽性者数のデータ
https://stopcovid19.metro.tokyo.lg.jp/cards/number-of-confirmed-cases/
を用いて、報告日ごとに集計した上で、日ごとの新規陽性者データとしました。

- 図1：日ごとの新規陽性者データ
![cov19-tokyo.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/7bcd8b86-a249-acf6-8f3d-1d4a10e95984.png)


図1の横軸は、最初の国内の陽性確認から起算した日数、縦軸は、新規陽性者数です。

#使用する自己回帰モデル
自己回帰モデル（ARモデル）そのものについては、ほかにも紹介記事が多数あるのでここでは割愛しますが、ざっくりいうと、最新のデータが、過去の時系列データの線形結合によって、再帰的に表現されるモデルのことを言います。微分方程式になじみのある方なら、線形微分方程式を差分近似した漸化式、というとイメージがつかみやすいかもしれません。微分方程式との大きな違いは、漸化式の各係数（AR係数）が、実際の時系列データに基づいて決定されることです。AR係数の決め方にはいろいろありますが、本記事では最小二乗法による決定法を採用しました。

ARモデルでは、最新のデータを決めるために、何時点前までのデータの線形結合を考えるか（次数）と、AR係数を決定するために、いくつのデータを使用するかを決定する必要があります。前者（次数）については、情報量基準を用いて、最適な次数を自動的に決める方法もありますが、今回は、直近1週間のデータ（7日間）としました。また、後者については、ARモデルは定常過程を想定しているので、普通は得られた時系列データすべてを対象にすることが多いようです。しかし、新規陽性者の推移はあまり定常的とは言えないので、ここでは適当にウィンドウを区切って、ある種の移動平均的にAR係数を決めていきます。このウィンドウの大きさは、28日間（7日間の4倍、1か月弱）としました。

#自己回帰スペクトル
AR係数が計算できると、そこから周波数領域の情報であるスペクトルを計算することができます。本記事では、代表的な時点における自己回帰スペクトルを計算し、過去の感染拡大状況と自己回帰スペクトルの関係を見ていくことにします。ただし、生の自己回帰スペクトルは、AR係数決定の分散の影響を受けてしまうので、DC成分（周波数0の成分）で正規化したスペクトルで考察します。

#各時点での自己回帰スペクトル
## 50日経過時点
- 図2：50日経過時点での自己回帰スペクトル
![arspec-50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/fbc8850a-46ab-3ea3-ca9e-605bc557169c.png)

まだ、1日の新規陽性者が数人から数十人だったころです。あの頃が懐かしいですね。
## 60日経過時点
- 図3：60日経過時点での自己回帰スペクトル
![arspec-60.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/231ae5ae-85cd-ec66-c71a-099e23005c64.png)

2020年3月下旬頃、陽性者が徐々に増え始めた、第1波の入口くらいの時期です。この時期の特徴として、周波数 0.1～0.15[1/day]のピークが急激に強くなっています。ちなみに、周期1週間に相当する周波数が0.14[1/day]なので、この頃から、1週間単位の周期性が顕在化してきていたことが分かります。

## 75日経過時点
- 図4：75日経過時点での自己回帰スペクトル
![arspec-75.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/e616bb71-1dc7-191f-101b-493f691cfe8d.png)

東京都に1回目の緊急事態宣言が発出された日（4/7）時点でのスペクトルです。1週間周期に相当するピークが相変わらず強いですが、その2倍の周波数（3.5日周期）にも強いピークが見られます。

## 85日経過時点
- 図5：85日経過時点での自己回帰スペクトル
![arspec-85.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/873c1582-8980-0859-118d-ec0ac11bf066.png)

第1波で最も新規陽性者数が多かった、4/17時点でのスペクトルです。意外だったのが、この頃になると、1週間周期に相当するピークが消えていました。おそらく、DC成分で正規化しているので、新規陽性者の平均値が増えるほど、相対的に1週間周期の変動が見えづらくなるのでしょう。

## 95日経過時点
- 図6：95日経過時点での自己回帰スペクトル
![arspec-95.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/85e78d66-0d98-2611-da61-aaa910f59148.png)

昨年のGW直前です。感染状況としてはまだまだピークでしたが、やはり1週間周期に相当するピークは見えません。

## 105日経過時点
- 図7：105日経過時点での自己回帰スペクトル
![arspec-105.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/7138255a-8a8b-de4d-3c13-f16a76d1db7a.png)

第1波が収束した頃のスペクトルです。0.1～0.15[1/day]付近のピークが復活していますが、大きさとしてはごく小さいです。

## 150日経過時点
- 図8：150日経過時点での自己回帰スペクトル
![arspec-150.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/f495d964-023a-e52a-212c-bc0b9fe11f0b.png)

まだまだ第2波の兆候が見えなかった時期のスペクトルです。

## 160日経過時点
- 図9：160日経過時点での自己回帰スペクトル
![arspec-160.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/387c6ead-825d-b30e-deec-da3b4be89f2e.png)

再び新規陽性者が増え始めた、第2波の入口です。1週間弱の周期のピークが出てきました。

## 170日経過時点
- 図10：170日経過時点での自己回帰スペクトル
![arspec-170.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/0d527381-f299-031d-ef4d-bb6368a589cf.png)

1週間と、その2倍の周波数のピークが顕在化し始めました。

## 180日経過時点
- 図11：180日経過時点での自己回帰スペクトル
![arspec-180.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/996f199e-770b-0be4-b1d2-de5925d96675.png)

1週間周期のピークがかなり明瞭になってきました。

## 190日経過時点
- 図11：190日経過時点での自己回帰スペクトル
![arspec-190.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/3a1160d7-b538-3533-8503-0741eab67cb6.png)

第2波の新規陽性者が最も多かったくらいの時点です。1週間相当のピークは見られますが、その大きさは弱まってきました。また、高周波のピークも見られません。

## 200日経過時点
- 図12：200日経過時点での自己回帰スペクトル
![arspec-200.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/860a7a8d-a7cb-9bba-13bb-2bd5c6a9e913.png)

1週間相当の弱いピークが見られます。

## 210日経過時点
- 図13：210日経過時点での自己回帰スペクトル
![arspec-210.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/09234ce6-e26b-5732-79fe-5655b7600324.png)

新規陽性者が減少傾向にあった時期です。1週間相当のピークはほとんど見えなくなりました。

## 230日経過時点
- 図14：230日経過時点での自己回帰スペクトル
![arspec-230.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/2d9d032b-d312-8379-963a-819b6d95bed3.png)

新規陽性者数が高止まりし始めた時期のスペクトルです。第1波収束時とは異なり、1週間ピークが比較的明瞭に見えます。

## 280日経過時点
- 図15：280日経過時点での自己回帰スペクトル
![arspec-280.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/b78c785a-e8c1-6bf2-c39c-ce57fb959d4b.png)

第3波に入る直前のスペクトルです。

## 300日経過時点
- 図16：300日経過時点での自己回帰スペクトル
![arspec-300.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/dbfcb11e-894f-6786-49b3-234fb0bf63af.png)

2020年11月下旬、新規陽性者が再び増加傾向にあった時期です。1週間相当のピークと、その整数倍周波数にピークが見られます。

## 320日経過時点
- 図17：320日経過時点での自己回帰スペクトル
![arspec-320.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/2570582d-bd51-cc73-8b71-73cc1d896362.png)

第3波の入口で一旦高止まりしていた時期のスペクトルです。ピークは小さくなっています。

## 330日経過時点
- 図18：330日経過時点での自己回帰スペクトル
![arspec-330.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/f9e84413-b6cb-a040-18f9-a816db67bc01.png)

第3波の拡大初期です。1週間ピークが顕在化しています。

## 340日経過時点
- 図19：340日経過時点での自己回帰スペクトル
![arspec-340.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/3f77f68c-9677-d1eb-c1c3-e9536ea3c860.png)

第3波拡大中の時期です。1週間と、その2倍周波数のピークが顕著に出ています。

## 350日経過時点
- 図20：350日経過時点での自己回帰スペクトル
![arspec-350.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/b38795d3-ac96-cf72-469b-232e50bd907a.png)

第3波のピーク時点です。1週間ピークは見られますが、拡大期と比べると弱まっています。

## 360日経過時点
- 図21：360日経過時点での自己回帰スペクトル
![arspec-360.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/89236ca1-361a-829c-5a3c-67e0dba25676.png)

第3波が収束し始めた時期です。1週間ピークはほとんど見られません。

#スペクトルの傾向と注意点
東京都の新規陽性者数の時系列データに対して、ARモデルを当てはめて計算したスペクトルには、概ね以下のような傾向がありました。

- 感染拡大初期から感染拡大期にかけては、1週間に相当する周波数、およびその整数倍の周波数のピークが強くなる。
- 感染の収束期には、これらの周波数帯のピークが弱まる。

ただし、これらの結果は、時系列データにARモデルを統計的に当てはめると、なんとなく以上のような傾向が見られたというだけであり、すべての場合において当てはまるとは限りません。また、本記事では、ある特定の時点におけるスペクトルのみを評価しましたが、ばらつきや特異値の影響を見るために、本当は前後数日間のデータを比較しながら傾向を分析したほうが良いでしょう。

ARモデルのひとつの適用事例として、時事的な話題を取り上げてみました。

#おわりに -気になる最新のスペクトル-

- 4/7時点
![arspec-440.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/78690680-76a1-df44-5f02-68cd7e0d81ef.png)

- 4/12時点
![arspec-445.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/e1c2e179-dd10-4ab8-b06d-da49eda40889.png)

- 4/17時点
![arspec-450.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/520a85c2-25d6-1d20-8548-0c305395b21f.png)

- 4/26時点
![arspec-459.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1064216/9388bc07-d858-ff4d-68de-a77c4d5d3fd5.png)

第4波においても、感染拡大初期から拡大期にかけて、ピーク強度が増えている傾向にあります。どのピークが卓越するかはデータを切り出す時点によって変わっていますが、4/26時点できれいに3つのピークが見えるので、これまでの傾向と一致すると仮定すると、もうしばらくは感染拡大傾向にありそうですね。。。