【時系列データ分析】東京都の平均気温を予測できるだろうか？

# はじめに
これまで機械学習に触れてきて、多くは画像系や自然言語系でしたが、ふと時系列データに触れてみたくなったので、時系列データ解析に入門しようと思いました。  
※記載しているコードは、全て1つのjupyter notebookでやってます。スクリプトの拡張子が`.py`になっているのは、単に`.ipynb`だと見えづらいからです。

## 時系列データとは  
言葉の通りで、時間の経過に応じて得られるデータを指します。  
例えば、有名なオープンデータセットの`iris`は、時間に関係なく、個体ごとのデータを用いて分析をすると思います。  
そうではなくて、タイトルにもありますが、ある地点の平均気温やあるお店の来場者数などは、時間(時期？)によってデータを得ることができます。

## どのようにして時系列データを分析するのか  
以下のサイトが参考になりました。数学的な説明もあって、非常にわかりやすかったです！
* [Pythonで時系列解析・超入門（その3） ARIMA系モデルで予測する方法](https://www.salesanalytics.co.jp/datascience/datascience087/)  
* [定常時系列の解析に使われるARMAモデル・SARIMAモデルとは？](https://ai-trend.jp/basic-study/time-series-analysis/sarima_model/)

## データセットの準備  
今回使用するデータセットは東京都の平均気温に関するデータセットです。  
[気象庁のHP](https://www.data.jma.go.jp/risk/obsdl/index.php)から、好きな都道府県の、どういった情報(雨や雪や日照時間など)かを選んで自由にダウンロードできます。ですが、一度にダウンロードできるデータには制限があるため、欲張ってあれこれと選ぶと一度にダウンロードできないです。  

今回使用するデータセットは以下の通りです。
* 地点は`東京`
* `2017/01/01〜2022/12/31`の6年分
* `平均気温`、`降水量の合計`、`日照時間`、`降雪量合計`、`平均雲量`の5つのカラムを選択
    * `平均気温`が目的変数
    * それ以外が説明変数

## 分析
### 使用したライブラリ/モジュール  
```sample.py
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_error, r2_score
```
### データの準備
```sample.py
csv_data_path = "data.csv" # ダウンロードしたデータセットのパス名
df = pd.read_csv(csv_data_path, header=3, encoding="shift_jis")
df
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/4993b804-12d7-dd5c-49bb-e0d17cbce52d.png)
オープンデータセットあるあるだと思いますが、案の定データセットの形がぐちゃぐちゃです(これでも、`header`で削っただけマシ)。カラム名に`Unnamed`と含まれているデータが、本来欲していたデータです。  
まずは、これをいい感じに整形するところからやりました。
```sample.py
df = pd.concat([df[["Unnamed: 0"]], df[["Unnamed: 1"]], df[["Unnamed: 4"]], 
                df[["Unnamed: 8"]], df[["Unnamed: 12"]], df[["Unnamed: 16"]]], axis=1)

columns = ["date", "mean_temp", "sum_rain", "sum_sunshine", "sum_snow", "mean_cloud"]
df = df.set_axis(columns, axis="columns")
df["date"] = pd.to_datetime(df["date"])
df = df.set_index("date")
df
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/06180e71-ab3e-c1fb-9f7e-941f93e826c4.png)
カラム名に関しては、僕の持ち前の英語センスで命名しました()。また、インデックスに日時を`datetime`型で付与しました。  
`datetime`にする一番の理由は、グラフを描画する際にx軸が`datetime`型になっていると、文字が潰れないように自動で軸の間隔を調整してくれるからです。

---
続いて、欠損値の確認をしてみました。
```sample.py
df.isnull().all().count() == df.shape[1]
>>> True
```
問題なさそうです。

---
次は、基本統計量を見ます。
```sample.py
df.describe().T
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/ac4512f6-a465-29b4-a21d-47a1d35d6051.png)
意外だと思ったことが、2017〜2022年の6年間で気温が0℃を下回った日がないんですね。当然、東京の中での観測地点にもよると思いますが。  

---
最後に、データを一度見てみることにします。
```sample.py
plt.plot(df.index, df["mean_temp"], color="b")
plt.xlabel("日付")
plt.ylabel("平均気温")
plt.tight_layout()
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/d3fc83ee-0861-7aa2-ea56-8a3b23aa6e44.png)
平均気温ということもあり、1年周期の規則性がありそうなことが目視できました。

### ARIMA
`ARIMA`モデルは、目的変数のみを用いて予測するモデルになります。そのため、訓練データとテストデータを`7:3`に分ける際には、目的変数のみ指定します。  
また、過去のデータを訓練データ、その先のデータをテストデータとするため、`shuffle=False`にしています。
```sample.py
arima_train, arima_test = train_test_split(df["mean_temp"], test_size=0.3, shuffle=False)
```
次に、最適なパラメータを探索してもらいました。
```sample.py
sm.tsa.arma_order_select_ic(arima_train, ic="aic")
>>> {'aic':               0            1            2
>>> 0  10722.418595  9139.041616  8290.307694
>>> 1   6711.856005  6699.668615  6582.189009
>>> 2   6705.837122  6715.227949  6550.301854
>>> 3   6647.056626  6548.559766  6580.593845
>>> 4   6623.662164  6550.391157  6551.075291,
>>> 'aic_min_order': (3, 1)}
```
`aic_min_order`の結果が`(3,1)`と探索してもらいました。このパラメータを用いて学習させます。 

---
```sample.py
arima_model = ARIMA(arima_train, order=(3,0,1)).fit()
```
`(3,0,1)`とパラメータを指定していますが、真ん中の0に関しては、今回手を加えていません。あくまで探索によって得られたパラメータは、`order`の第1, 3引数に対応します。

学習したモデルの概要を見ます。
```sample.py
arima_model.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/31477aa0-302b-9d3f-3ad6-6864322bafda.png)
画像の右上に`AIC`や`BIC`とありますが、このあたりが評価指数になります。ここらの値が小さい方が良い精度で学習できていることになります。

---
最後に学習した結果を用いて、テストデータとの予測を比較します。
```sample.py
pred_arima = arima_model.forecast(steps=arima_test.shape[0])

plt.plot(arima_test.index, arima_test, color="b", label="テストデータ")
plt.plot(arima_test.index, pred_arima, color="r", label="予測データ")
plt.xlabel("日付")
plt.ylabel("平均気温")
plt.title(f"MAE: {mean_absolute_error(arima_test,pred_arima):.2f}, R2: {r2_score(arima_test,pred_arima):.2f}")
plt.legend()
plt.tight_layout()
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/89837afe-31c3-3b8e-e0a0-e6cab8855dc8.png)
。。
むむっ。きびしっ。  
頑張って学習したのだろうと思いつつ、これが描画されたときには夜に一人で爆笑しました。

### SARIMA
`SARIMA`は、**S**easonal+ARIMAを意味していて、ARIMAによる予測時に周期性(=Seasonal)も考慮する手法になります。今回のデータにぴったりそうです。
```sample.py
%%time
sarima_model = sm.tsa.SARIMAX(endog=arima_train, order=(3,0,1), seasonal_order=(0,1,0,365)).fit()
```
`ARIMA`と同じく、目的変数のみを使って予測するため、引数には同じ訓練データを渡しています。また、`seasonal_order`というパラメータが新しく増えました。  
このパラメータの第4引数で、1周期がデータ何個分かを指定しています。今回の場合は、1年(=365日)周期なので`365`ですが、データセットが月別だったりする場合は`12`になったりします([この引数の数字が大きいほど、学習に時間がかかるらしいです](https://github.com/statsmodels/statsmodels/issues/6033))。  
第4引数が`365`になったおかげで時間がかかるため、`seasonal_order`の他の引数はテキトーに決め打ちしました。ここの引数は、グリッドサーチなどで探す方法が一般的なんでしょうか？  

`SARIMA`モデルなのに、使用しているモジュールが`SARIMAX`な点は問題ないです。`SARIMAX`として使用する場合は、第2引数にもう一つデータを指定する必要があり、それがない場合は、`SARIMA`モデルになるらしいです。

---
モデルの概要を見てみます。
```sample.py
sarima_model.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/815a50bd-6a51-9209-0fe0-dfe1086f5d1d.png)
`ARIMA`モデルと`SARIMA`モデルの評価指標である`AIC`を比較すると、それぞれ`6549`、`5730`となり、`SARIMA`モデルの方が良い結果になっていそうです。  
予測結果を見てみます。
```sample.py
pred_sarima = sarima_model.forecast(steps=arima_test.shape[0])

plt.plot(arima_test.index, arima_test, color="b", label="テストデータ")
plt.plot(arima_test.index, pred_sarima, color="r", label="予測データ")
plt.xlabel("日付")
plt.ylabel("平均気温")
plt.title(f"MAE: {mean_absolute_error(arima_test,pred_sarima):.2f}, R2: {r2_score(arima_test,pred_sarima):.2f}")
plt.legend()
plt.tight_layout()
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/4d66e214-9dfe-5217-b7bb-f06203f4af6e.png)
結構ギザギザしている線にしては、頑張って予測できている気がします。$R^2$スコアが`0.74`とそこまで高いわけでもないですが、パラメータをいじったりせずにこれなので、まだまだ改善はできそうです。

### SARIMAX
`SARIMAX`は、`説明変数を用いて目的変数を予測する`+`周期性を考慮する`モデルです。  
`ARIMA`や`SARIMA`は、目的変数のみで予測するモデルでしたが、末尾に`X`がつくことで、説明変数を用いて予測するモデルになります。今回は実施しませんでしたが、`ARIMAX`たるモデルもありました。

まずは、学習データセットとテストデータセットに分けるところから始めます。`ARIMA`や`SARIMA`のときと同様の割合で、過去データを訓練データ、未来データをテストデータにしました。
```sample.py
X = df.drop(columns=["mean_temp"])
y = df["mean_temp"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, shuffle=False)
```

---
次に学習します。
```sample.py
%%time
sarimax_model = sm.tsa.SARIMAX(endog=y_train, exog=X_train, order=(3,0,1), seasonal_order=(0,1,0,365)).fit()
```
`SARIMA`のときと異なり、第2引数に`X_train`を加えました。他のパラメータ等は`SARIMA`と全く同じです。  
余談ですが、`SARIMA`のときは学習に9分くらいかかったのですが、`SARIMAX`の学習には45分くらいかかりました。やはり前述した通り、1周期の数字が大きいと時間がかかりますね。  

モデルの概要を見てみます。
```sample.py
sarimax_model.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/bd830d51-22ad-140d-a7e8-0baaecf0a625.png)
改めて、これまでのモデルと`AIC`を比べてみると、
|モデル|AIC|
|:-:|:-:|
|ARIMA|6548|
|SARIMA|5729|
|SARIMAX|5482|

という結果になり、このモデルが一番評価指標が低いことが確認できました。  

結果を見てみます。
```sample.py
pred_sarimax = sarimax_model.forecast(steps=X_test.shape[0], exog=X_test)

plt.plot(arima_test.index, y_test, color="b", label="テストデータ")
plt.plot(arima_test.index, pred_sarimax, color="r", label="予測データ")
plt.xlabel("日付")
plt.ylabel("平均気温")
plt.title(f"MAE: {mean_absolute_error(y_test,pred_sarimax):.2f}, R2: {r2_score(y_test,pred_sarimax):.2f}")
plt.legend()
plt.tight_layout()
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1939085/ca72a285-fbad-f847-7ecd-b12207435a12.png)
先ほど試した`SARIMA`は、$R^2$スコアが`0.74`だったため、精度が向上したと言えそうです。

## 東京都の平均気温を予測できるだろうか？
できる気がしますが、今の僕の知識だと言い切れるわけではないので何とも言えないです..  
そもそも論、天気を予測してくれる気象予報士さんの存在のおかげで、過去のデータから予測ができなくても誰も困らないですね。毎日ありがとうございますmm

# おわりに
キャッチーなタイトルに憧れてタイトルを決めましたが、実際にはそれに及ぶほどの知識や知見がなかったため、理解を深めた上で、改めて挑戦してみたいですね。  

`ARIMA`、`SARIMA`、`SARIMAX`とモデルを紹介しましたが、これらの読み方ってそれぞれ`アリマ`、`サリマ`、`サリマックス`でよろしいのでしょうか？  
個人的には`サリマックス`がツボなんですが、~~ちょっとダサい気がしていて~~、カッコつけて`シーズナル・アリマ・エックス`と呼んでいます。有識者の方がいましたらご指摘して頂きたいです。