LoginSignup
3
4
@T_death(T です)

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

Last updated at Posted at 2023-09-01

はじめに

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

時系列データとは

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

どのようにして時系列データを分析するのか

以下のサイトが参考になりました。数学的な説明もあって、非常にわかりやすかったです！

データセットの準備

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

今回使用するデータセットは以下の通りです。

  • 地点は東京
  • 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
オープンデータセットあるあるだと思いますが、案の定データセットの形がぐちゃぐちゃです(これでも、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
カラム名に関しては、僕の持ち前の英語センスで命名しました()。また、インデックスに日時をdatetime型で付与しました。
datetimeにする一番の理由は、グラフを描画する際にx軸がdatetime型になっていると、文字が潰れないように自動で軸の間隔を調整してくれるからです。

続いて、欠損値の確認をしてみました。

sample.py
df.isnull().all().count() == df.shape[1]
>>> True

問題なさそうです。

次は、基本統計量を見ます。

sample.py
df.describe().T

image.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
平均気温ということもあり、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
画像の右上に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
。。
むむっ。きびしっ。
頑張って学習したのだろうと思いつつ、これが描画されたときには夜に一人で爆笑しました。

SARIMA

SARIMAは、Seasonal+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になったりします(この引数の数字が大きいほど、学習に時間がかかるらしいです)。
第4引数が365になったおかげで時間がかかるため、seasonal_orderの他の引数はテキトーに決め打ちしました。ここの引数は、グリッドサーチなどで探す方法が一般的なんでしょうか？

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

モデルの概要を見てみます。

sample.py
sarima_model.summary()

image.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
結構ギザギザしている線にしては、頑張って予測できている気がします。$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
改めて、これまでのモデルと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
先ほど試したSARIMAは、$R^2$スコアが0.74だったため、精度が向上したと言えそうです。

東京都の平均気温を予測できるだろうか？

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

おわりに

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

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

3
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
4