ブログ

You are currently viewing ARIMAモデルで時系列データの予測をします〜Python

ARIMAモデルで時系列データの予測をします〜Python

以前線形回帰を使って時系列データ予測の記事を書きました。但し、線形回帰は時間系列の関連性を考慮しない、つまり各データポイントは独立であるという仮説の基で提案されたアルゴリズムです。時系列データを扱うにはARIMA等の専用のアルゴリズムがあります。ARIMAモデルには線形回帰が導入されていますが、時系列データを処理できるように設計されています。

今回はARIMAモデルについて説明して行きます。

  • ARIMAモデルの概要とのARIMAモデルパラメータの決め方
  • ARIMAモデルを使って時系列データの解析方法
  • ARIMAモデルを使って先の時刻の予測方法

ARIMAモデルって何？

ARIMAモデルはAutoRegressive Integrated Moving Averageの略称です。時系列データ予測アルゴリズムにおいて非常に有名です。名前のとおり回帰（Regression）アルゴリズムです。但し、線形回帰が過去を考慮しないに対して、ARIMAモデルは過去の値を見て時系列データを解釈します。

ARIMAモデルに重要なのは名前にあるAR（AutoRegressive ）とMA（Moving Average）です。

  • AR＝自動回帰：観測値Y(t)（t時刻の値）と遅れた観測値Y(t-1)…Y(t-p)（t-p時刻の値）の依存関係のモデル
  • MA＝移動平均：観測値Y(t)（t時刻の値）と遅れた観測値の誤差E(t-1)…E(t-p)（t-p時刻の値）の依存関係のモデル

従って、ARIMAモデルは Y(p,d,q)＝定数＋ARモデル＋MAモデル の方程式で表すことができます。数式の詳細は割愛します。

p,d,qはARIMAモデルのパラメータですので、良いARIMAモデルを作るのに最適なp,d,qを探すのが重要な仕事になります。

  • p：pはARモデルのパラメータで、遅れた観測値の数を示します。
  • d：観測値の段差です。パラメータdはデータから周期性を取り除き、定常化します。
  • q：qはMAモデルのパラメータで、移動平均のウィンドサイズを表します。

おそらくここまで理論の話してもDeveloperの方には混乱しますので、実際のデータを見ながらp,d,qの値を決めていこうと思います。

東京電力電力量データ

今回扱うデータは東京電力の公開電力量データです。１時間間隔で電力の使用量が記録されています。データには時刻と電力量のみで、シンプルなデータセットです。2020/1/1から現在までのデータを使います。早速データを読み込みます。

import pandas as pd
# Read Energy data from tepco
energy_data_url="http://www.tepco.co.jp/forecast/html/images/juyo-2020.csv"
energy_data_df = pd.read_csv(energy_data_url, skiprows=[0, 1], encoding="shift-jis")
energy_data_df.columns = ["DATE", "TIME", "energy"]
energy_data_df["time"] = energy_data_df["DATE"] + " " + energy_data_df["TIME"] + ":00"

from datetime import datetime
energy_data_series = energy_data_df["energy"]
energy_data_series.index = [datetime.strptime(d, '%Y/%m/%d %H:%M:%S') for d in energy_data_df["time"].values]

from matplotlib import pyplot
energy_data_series.plot()
pyplot.show()

データには明らかのトレンドがあります。当然ですが、深夜と朝は寝る時間ですから電気をあんまり使わないし、May（五月）になると冷房が要らなりますから、トレンドが発生します。統計的にこういうデータを定常ではないデータと呼びます。時系列処理の前処理として通常、データを定常化する必要があり、定常化の基準はデータの段差を見ます。これをやるのがARIMAモデルのパラメータdです。

時系列データを定常化するためによく使われるのが自己相関（autocorrelation）図です。データを確認し、段差を決めます。自己相関図はpandasの関数を使えばすぐにできますので、やってみましょう。以下段差が1/2/3の場合のそれぞれの自己相関図をプロットしてあります。左が段差を取ったデータのプロットで、右が自己相関図になります。

from pandas.plotting import autocorrelation_plot
import numpy as np

# Original Series
fig, axes = pyplot.subplots(3, 2, figsize=(15,15))

# 1st Differencing
axes[0, 0].plot(np.diff(energy_data_series.to_numpy())); axes[0, 0].set_title('1st Order Differencing')
autocorrelation_plot(np.diff(energy_data_series.to_numpy()), ax=axes[0, 1])

# 2nd Differencing
axes[1, 0].plot(np.diff(np.diff(energy_data_series.to_numpy()))); axes[1, 0].set_title('2nd Order Differencing')
autocorrelation_plot(np.diff(np.diff(energy_data_series.to_numpy())), ax=axes[1, 1])

# 3rd Differencing
axes[2, 0].plot(np.diff(np.diff(np.diff(energy_data_series.to_numpy())))); axes[2, 0].set_title('3rd Order Differencing')
autocorrelation_plot(np.diff(np.diff(np.diff(energy_data_series.to_numpy()))), ax=axes[2, 1])

pyplot.show()
  • まず段差のパラメータdを確定しましょう。ご覧のとおり、段差を2（2nd Order Differencing）にした場合に、データの変動が比較的に安定しているので、dを2にしましょう。
  • ARモデルのパラメータpを確定しましょう。段差を2にしましたので、段差が2の場合の自己相関図のみに注目しましょう。図を拡大して見てみましょう。set_xlimでLagの範囲を0~10に絞って自己相関図を再プロットします。Lagが1の時にAutocorrelationが有意水準（点線）を超えていましたが、Lagが2になると、有意水準以下になっていましたので、pを1にするのが妥当でしょう。
pyplot.figure(figsize=(8,8))
autocorrelation_plot(np.diff(np.diff(energy_data_series.to_numpy()))) \
    .set_xlim([0, 10])
pyplot.show()
  • MRモデルのパラメータqを確定しましょう。qは移動平均のサイズですが、株、仮想通貨等に詳しい方はすぐに理解できると思いますが、ここでは自己相関図がスムーズですので、とりあえずqを0にしましょう。

上記のとおり、p=1、d=2、q=0として初期設定を選びました。実際の案件ではモデルを最適化するためにp,d,qをチューニングしましょう。

では、p、d、qを確定しましたので、残りはstatsmodelsライブラリーのARIMAモデルを呼ぶだけです。ARIMA(p=1,d=2,q=0)で学習データにフィットします。

データを学習データとテストデータに分けますが、通常の6:2:2の交差検証を含めた分割方法ではなく、単純に後ほどモデルの予測結果を観測値と比較したいので、データの最後の数十件をテストデータセットとし、残り全部を学習データとします。

from statsmodels.tsa.arima_model import ARIMA
from matplotlib import pyplot
from pandas import DataFrame
from sklearn.model_selection import train_test_split

train, test = train_test_split(energy_data_series, test_size=0.005, shuffle=False)
# fit model
model = ARIMA(train, order=(1,2,0))
model_fit = model.fit(disp=0)
print(model_fit.summary())
# plot residual errors
residuals = DataFrame(model_fit.resid)
residuals.plot()
pyplot.show()
residuals.plot(kind='kde')
pyplot.show()
print(residuals.describe())
残差
残差の密度関数

ARIMA Model ResultsはARIMAモデルがデータにフィットした概要を示しています。残差（residual error）は0付近で変動しながらも幅がまだ大きい。残差の密度関数は正規分布の形をしていて、中心は0ではありません。

ARIMAモデルで予測します

ではARIMAモデルの学習は完了したところで、モデルを使って予測してみましょう。one-step forward予測法を使います。つまり、一回の予測では、一時刻先の予測しかできませんが、常に最新のデータで学習したモデルを使いますので、予測精度を向上できます。

print(len(test)) 
24

上記のようにテストデータセットに24件のデータがありますので、ARIMAモデルを使って24歩のone-step forward予測を行います。

history = [x for x in train]
predictions = list()
for t in range(len(test)):
    model = ARIMA(history, order=(1,2,0))
    model_fit = model.fit(disp=0)
    output = model_fit.forecast()
    yhat = output[0]
    predictions.append(yhat)
    obs = test[t]
    history.append(obs)
# plot
pyplot.plot(test.index, test.values)
pyplot.plot(test.index, predictions, color='red')
pyplot.show()

予測結果としては悪くないですね。チューニングは全然行っていませんので、ベストモデルではありません。興味がある方は時系列データの前処理、Cross-Validation、GridSearch等の手法を使ってチューニングしてみましょう。

もう一つ重要なチューニング基準は残差です。先ほど残差と残差の密度関数を出しましたね。残差が0付近で変化し分散が大きくないこと及び残差の密度関数の中心はできるだけ0に接近することを気にしながらチューニングを行えば良いです。

まとめ

この記事ではARIMAモデルを使って時系列の予測モデルを開発する方法を紹介しました：

  • ARIMAモデルの概要とのARIMAモデルパラメータの決め方
  • ARIMAモデルを使って時系列データの解析方法
  • ARIMAモデルを使って先の時刻の予測方法
タグ: ,