Python3でARIMAを使用した時系列予測のガイド

序章

時系列は、将来の値を予測する機会を提供します。以前の値に基づいて、時系列を使用して、経済、天気、および容量計画の傾向を予測できます。時系列データの特定の特性は、通常、特殊な統計手法が必要であることを意味します。

このチュートリアルでは、時系列の信頼できる予測を作成することを目指します。まず、自己相関、定常性、季節性の概念を紹介して説明し、時系列予測で最も一般的に使用される方法の1つであるARIMAの適用に進みます。

時系列の将来のポイントをモデル化および予測するためにPythonで使用できる方法の1つは、 SARIMAX として知られています。これは、 Seasonal AutoRegressive Integrated Moving Averages with eXogenousregressorsの略です。ここでは、主にARIMAコンポーネントに焦点を当てます。これは、時系列データを適合させて、時系列の将来のポイントをよりよく理解および予測するために使用されます。

前提条件

このガイドでは、ローカルデスクトップまたはリモートサーバーのいずれかで時系列分析を行う方法について説明します。大規模なデータセットの操作はメモリを大量に消費する可能性があるため、いずれの場合も、このガイドの計算の一部を実行するには、コンピュータに少なくとも2GBのメモリが必要です。

このチュートリアルを最大限に活用するには、時系列と統計にある程度精通していると役立ちます。

このチュートリアルでは、 JupyterNotebookを使用してデータを操作します。まだお持ちでない場合は、チュートリアルに従って、Python3用のJupyterNotebookをインストールおよびセットアップする必要があります。

ステップ1—パッケージのインストール

時系列予測用の環境をセットアップするために、最初にローカルプログラミング環境またはサーバーベースのプログラミング環境に移動しましょう。

cd environments

. my_env/bin/activate

ここから、プロジェクトの新しいディレクトリを作成しましょう。これをARIMAと呼び、ディレクトリに移動します。プロジェクトを別の名前で呼ぶ場合は、ガイド全体で必ずARIMAの代わりに自分の名前を使用してください

mkdir ARIMA
cd ARIMA

このチュートリアルでは、warnings、itertools、pandas、numpy、matplotlib、およびstatsmodelsライブラリが必要です。 warningsおよびitertoolsライブラリは標準のPythonライブラリセットに含まれているため、インストールする必要はありません。

他のPythonパッケージと同様に、これらの要件はpipでインストールできます。これで、pandas、statsmodels、およびデータプロットパッケージmatplotlibをインストールできます。それらの依存関係もインストールされます。

pip install pandas numpy statsmodels matplotlib

この時点で、インストールされたパッケージの操作を開始する準備が整いました。

ステップ2—パッケージのインポートとデータのロード

データの操作を開始するには、JupyterNotebookを起動します。

jupyter notebook

新しいノートブックファイルを作成するには、右上のプルダウンメニューから New > Python3を選択します。

これにより、ノートブックが開きます。

ベストプラクティスとして、ノートブックの上部に必要なライブラリのインポートから始めます。

import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')

また、プロットのmatplotlibスタイルをファイブサーティエイトに定義しました。

1958年3月から2001年12月までのCO2サンプルを収集した、「米国マウナロア天文台の連続大気サンプルからの大気CO2」というデータセットを使用します。このデータは次のように取り込むことができます。

data = sm.datasets.co2.load_pandas()
y = data.data

先に進む前に、データを少し前処理してみましょう。週次データは時間が短いため扱いにくい場合があるため、代わりに月次平均を使用しましょう。 resample関数で変換します。簡単にするために、 fillna（）関数を使用して、時系列に欠落値がないことを確認することもできます。

# The 'MS' string groups the data in buckets by start of the month
y = y['co2'].resample('MS').mean()

# The term bfill means that we use the value before filling in missing values
y = y.fillna(y.bfill())

print(y)

Outputco2
1958-03-01  316.100000
1958-04-01  317.200000
1958-05-01  317.433333
...
2001-11-01  369.375000
2001-12-01  371.020000

この時系列eをデータの視覚化として調べてみましょう。

y.plot(figsize=(15, 6))
plt.show()

データをプロットすると、いくつかの識別可能なパターンが表示されます。時系列には明らかな季節性パターンがあり、全体的に増加傾向があります。

時系列の前処理の詳細については、「 Python 3を使用した時系列の視覚化ガイド」を参照してください。ここでは、上記の手順について詳しく説明しています。

データを変換して調査したので、ARIMAを使用した時系列予測に移りましょう。

ステップ3—ARIMA時系列モデル

時系列予測で使用される最も一般的な方法の1つは、ARIMAモデルとして知られています。これは、 A utoreg R essive I統合M[ X167X] oving A平均。 ARIMAは、時系列データを適合させて、時系列の将来のポイントをよりよく理解または予測できるモデルです。

ARIMAモデルのパラメーター化に使用される3つの異なる整数（p、d、q）があります。そのため、ARIMAモデルはARIMA(p, d, q)という表記で表されます。これらの3つのパラメーターを合わせて、データセットの季節性、傾向、およびノイズを説明します。

pは、モデルの自己回帰部分です。これにより、過去の値の影響をモデルに組み込むことができます。直感的には、これは、過去3日間暖かくなった場合、明日は暖かくなる可能性があると述べるのと似ています。
dは、モデルの統合部分です。これには、差異の量を組み込んだモデル内の用語が含まれます（つまり、時系列に適用するために、現在の値から差し引く過去の時点の数）。直感的には、過去3日間の気温差が非常に小さければ、明日は同じ気温になる可能性が高いと言うのと似ています。
qは、モデルの移動平均部分です。これにより、モデルのエラーを、過去の以前の時点で観測されたエラー値の線形結合として設定できます。

季節的な影響を扱うときは、季節性ARIMAを使用します。これはARIMA(p,d,q)(P,D,Q)sと表記されます。ここで、(p, d, q)は上記の非季節パラメータですが、(P, D, Q)は同じ定義に従いますが、時系列の季節コンポーネントに適用されます。 sという用語は、時系列の周期性です（4は四半期ごと、12は年ごとなど）。

季節的なARIMA法は、複数の調整パラメーターが関係しているため、気が遠くなるように見える場合があります。次のセクションでは、季節性ARIMA時系列モデルの最適なパラメーターセットを特定するプロセスを自動化する方法について説明します。

ステップ4—ARIMA時系列モデルのパラメーター選択

時系列データを季節ARIMAモデルに適合させる場合、最初の目標は、対象のメトリックを最適化するARIMA(p,d,q)(P,D,Q)sの値を見つけることです。この目標を達成するための多くのガイドラインとベストプラクティスがありますが、ARIMAモデルの正しいパラメーター化は、ドメインの専門知識と時間を必要とする骨の折れる手動プロセスになる可能性があります。 Rなどの他の統計プログラミング言語はこの問題を解決する自動化された方法を提供しますが、それらはまだPythonに移植されていません。このセクションでは、Pythonコードを記述して、ARIMA(p,d,q)(P,D,Q)s時系列モデルに最適なパラメーター値をプログラムで選択することでこの問題を解決します。

「グリッド検索」を使用して、パラメーターのさまざまな組み合わせを繰り返し探索します。パラメータの組み合わせごとに、statsmodelsモジュールのSARIMAX()関数を使用して新しい季節ARIMAモデルを適合させ、その全体的な品質を評価します。パラメータの全体像を調査した後、最適なパラメータセットは、対象の基準に対して最高のパフォーマンスをもたらすものになります。評価したいパラメータのさまざまな組み合わせを生成することから始めましょう。

# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)

# Generate all different combinations of p, q and q triplets
pdq = list(itertools.product(p, d, q))

# Generate all different combinations of seasonal p, q and q triplets
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]

print('Examples of parameter combinations for Seasonal ARIMA...')
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))

OutputExamples of parameter combinations for Seasonal ARIMA...
SARIMAX: (0, 0, 1) x (0, 0, 1, 12)
SARIMAX: (0, 0, 1) x (0, 1, 0, 12)
SARIMAX: (0, 1, 0) x (0, 1, 1, 12)
SARIMAX: (0, 1, 0) x (1, 0, 0, 12)

これで、上記で定義した3つのパラメーターを使用して、さまざまな組み合わせでARIMAモデルをトレーニングおよび評価するプロセスを自動化できます。統計と機械学習では、このプロセスはモデル選択のためのグリッド検索（またはハイパーパラメーター最適化）として知られています。

さまざまなパラメーターに適合した統計モデルを評価および比較する場合、データにどの程度適合しているか、または将来のデータポイントを正確に予測する能力に基づいて、それぞれを相互にランク付けできます。 AIC（赤池情報量基準）の値を使用します。これは、statsmodelsを使用して適合されたARIMAモデルで便利に返されます。 AICは、モデルの全体的な複雑さを考慮しながら、モデルがデータにどの程度適合しているかを測定します。多くの機能を使用しながらデータに非常によく適合するモデルには、同じ適合度を達成するために使用する機能が少ないモデルよりも大きなAICスコアが割り当てられます。したがって、AICの値が最も低くなるモデルを見つけることに関心があります。

以下のコードチャンクは、パラメーターの組み合わせを繰り返し、statsmodelsのSARIMAX関数を使用して、対応する季節ARIMAモデルに適合させます。ここで、order引数は(p, d, q)パラメーターを指定し、seasonal_order引数は季節ARIMAモデルの(P, D, Q, S)季節コンポーネントを指定します。各SARIMAX()モデルをフィッティングした後、コードはそれぞれのAICスコアを出力します。

warnings.filterwarnings("ignore") # specify to ignore warning messages

for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            mod = sm.tsa.statespace.SARIMAX(y,
                                            order=param,
                                            seasonal_order=param_seasonal,
                                            enforce_stationarity=False,
                                            enforce_invertibility=False)

            results = mod.fit()

            print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
        except:
            continue

一部のパラメータの組み合わせは数値の仕様ミスにつながる可能性があるため、警告メッセージの過負荷を回避するために、警告メッセージを明示的に無効にしました。これらの仕様の誤りは、エラーを引き起こして例外をスローする可能性もあるため、これらの例外をキャッチし、これらの問題の原因となるパラメーターの組み合わせを無視するようにしてください。

上記のコードは次の結果をもたらすはずです。これには時間がかかる場合があります。

OutputSARIMAX(0, 0, 0)x(0, 0, 1, 12) - AIC:6787.3436240402125
SARIMAX(0, 0, 0)x(0, 1, 1, 12) - AIC:1596.711172764114
SARIMAX(0, 0, 0)x(1, 0, 0, 12) - AIC:1058.9388921320026
SARIMAX(0, 0, 0)x(1, 0, 1, 12) - AIC:1056.2878315690562
SARIMAX(0, 0, 0)x(1, 1, 0, 12) - AIC:1361.6578978064144
SARIMAX(0, 0, 0)x(1, 1, 1, 12) - AIC:1044.7647912940095
...
...
...
SARIMAX(1, 1, 1)x(1, 0, 0, 12) - AIC:576.8647112294245
SARIMAX(1, 1, 1)x(1, 0, 1, 12) - AIC:327.9049123596742
SARIMAX(1, 1, 1)x(1, 1, 0, 12) - AIC:444.12436865161305
SARIMAX(1, 1, 1)x(1, 1, 1, 12) - AIC:277.7801413828764

私たちのコードの出力は、SARIMAX(1, 1, 1)x(1, 1, 1, 12)が277.78の最低のAIC値を生成することを示唆しています。したがって、これまで検討したすべてのモデルの中で、これが最適なオプションであると見なす必要があります。

ステップ5—ARIMA時系列モデルの適合

グリッド検索を使用して、時系列データに最適なモデルを生成するパラメーターのセットを特定しました。この特定のモデルをさらに詳しく分析することができます。

まず、最適なパラメータ値を新しいSARIMAXモデルにプラグインします。

mod = sm.tsa.statespace.SARIMAX(y,
                                order=(1, 1, 1),
                                seasonal_order=(1, 1, 1, 12),
                                enforce_stationarity=False,
                                enforce_invertibility=False)

results = mod.fit()

print(results.summary().tables[1])

Output==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.3182      0.092      3.443      0.001       0.137       0.499
ma.L1         -0.6255      0.077     -8.165      0.000      -0.776      -0.475
ar.S.L12       0.0010      0.001      1.732      0.083      -0.000       0.002
ma.S.L12      -0.8769      0.026    -33.811      0.000      -0.928      -0.826
sigma2         0.0972      0.004     22.634      0.000       0.089       0.106
==============================================================================

SARIMAXの出力から得られるsummary属性は、かなりの量の情報を返しますが、係数の表に注目します。 coef列には、重量が表示されます（つまり、重要性）各機能の重要性と、各機能が時系列にどのように影響するか。 P>|z|列は、各特徴の重みの重要性を通知します。ここで、各重みのp値は0.05より低いか近いため、モデルにすべてを保持するのが妥当です。

季節性ARIMAモデル（およびその他のモデル）を適合させる場合、モデル診断を実行して、モデルによって行われた仮定に違反していないことを確認することが重要です。 plot_diagnosticsオブジェクトを使用すると、モデル診断をすばやく生成し、異常な動作を調査できます。

results.plot_diagnostics(figsize=(15, 12))
plt.show()

私たちの主な関心事は、モデルの残差が無相関であり、通常はゼロ平均で分布していることを確認することです。季節ARIMAモデルがこれらの特性を満たさない場合は、さらに改善できることを示しています。

この場合、モデル診断は、モデルの残差が以下に基づいて正規分布していることを示しています。

右上のプロットでは、赤いKDE線がN(0,1)線（N(0,1)）に密接に続いていることがわかります。これは、平均0および1の標準偏差）。これは、残差が正規分布していることを示す良い指標です。
左下のqq-plotは、残差（青い点）の順序付けられた分布が、N(0, 1)の標準正規分布から取得されたサンプルの線形傾向に従うことを示しています。繰り返しますが、これは残差が正規分布していることを強く示しています。
時間の経過に伴う残差（左上のプロット）は、明らかな季節性を表示せず、ホワイトノイズのように見えます。これは自己相関によって確認されます（すなわちコレログラム）右下のプロット。これは、時系列の残差がそれ自体の遅延バージョンとの相関が低いことを示しています。

これらの観察結果から、私たちのモデルは、時系列データを理解し、将来の値を予測するのに役立つ十分な適合を生成すると結論付けることができます。

十分な適合性がありますが、季節性ARIMAモデルの一部のパラメーターを変更して、モデルの適合性を向上させることができます。たとえば、グリッド検索ではパラメータの組み合わせの制限されたセットのみが考慮されたため、グリッド検索を拡大すると、より適切なモデルが見つかる可能性があります。

ステップ6—予測の検証

時系列のモデルを取得しました。これを使用して予測を作成できます。まず、予測値を時系列の実際の値と比較します。これは、予測の正確さを理解するのに役立ちます。 get_prediction()およびconf_int()属性を使用すると、時系列の予測の値と関連する信頼区間を取得できます。

pred = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=False)
pred_ci = pred.conf_int()

上記のコードでは、予測を1998年1月に開始する必要があります。

dynamic=False引数は、一歩先の予測を生成することを保証します。つまり、各ポイントの予測は、そのポイントまでの完全な履歴を使用して生成されます。

CO2時系列の実際の値と予測値をプロットして、どれだけうまくいったかを評価できます。日付インデックスをスライスして、時系列の最後にズームインしたことに注目してください。

ax = y['1990':].plot(label='observed')
pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)

ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.2)

ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
plt.legend()

plt.show()

全体として、私たちの予測は真の値と非常によく一致しており、全体的な増加傾向を示しています。

予測の正確さを定量化することも役立ちます。予測の平均誤差を要約するMSE（平均二乗誤差）を使用します。予測値ごとに、真の値までの距離を計算し、結果を2乗します。全体の平均を計算するときに正/負の差が互いに打ち消し合わないように、結果を2乗する必要があります。

y_forecasted = pred.predicted_mean
y_truth = y['1998-01-01':]

# Compute the mean square error
mse = ((y_forecasted - y_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

OutputThe Mean Squared Error of our forecasts is 0.07

一歩先を行く予測のMSEは、0.07の値を生成します。これは、0に近いため、非常に低くなります。 MSEが0の場合、推定量はパラメーターの観測値を完全な精度で予測します。これは理想的なシナリオですが、通常は不可能です。

ただし、動的予測を使用すると、真の予測力をより適切に表すことができます。この場合、特定の時点までの時系列の情報のみを使用し、その後、以前の予測時点の値を使用して予測が生成されます。

以下のコードチャンクでは、1998年1月以降に動的予測と信頼区間の計算を開始するように指定しています。

pred_dynamic = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=True, full_results=True)
pred_dynamic_ci = pred_dynamic.conf_int()

時系列の観測値と予測値をプロットすると、動的予測を使用した場合でも、全体的な予測が正確であることがわかります。すべての予測値（赤い線）はグラウンドトゥルース（青い線）と非常によく一致しており、予測の信頼区間内に十分収まっています。

ax = y['1990':].plot(label='observed', figsize=(20, 15))
pred_dynamic.predicted_mean.plot(label='Dynamic Forecast', ax=ax)

ax.fill_between(pred_dynamic_ci.index,
                pred_dynamic_ci.iloc[:, 0],
                pred_dynamic_ci.iloc[:, 1], color='k', alpha=.25)

ax.fill_betweenx(ax.get_ylim(), pd.to_datetime('1998-01-01'), y.index[-1],
                 alpha=.1, zorder=-1)

ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')

plt.legend()
plt.show()

ここでも、MSEを計算することにより、予測の予測パフォーマンスを定量化します。

# Extract the predicted and true values of our time series
y_forecasted = pred_dynamic.predicted_mean
y_truth = y['1998-01-01':]

# Compute the mean square error
mse = ((y_forecasted - y_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

OutputThe Mean Squared Error of our forecasts is 1.01

動的予測から得られた予測値は、1.01のMSEをもたらします。これは、時系列からのより少ない履歴データに依存していることを考えると、予想される1ステップ先よりもわずかに高くなります。

一歩先を行く予測と動的な予測の両方で、この時系列モデルが有効であることが確認されます。ただし、時系列予測に関する関心の多くは、将来の値を前もって予測する機能です。

ステップ7—予測の作成と視覚化

このチュートリアルの最後のステップでは、季節性ARIMA時系列モデルを活用して将来の値を予測する方法について説明します。時系列オブジェクトのget_forecast()属性は、指定されたステップ数先の予測値を計算できます。

# Get forecast 500 steps ahead in future
pred_uc = results.get_forecast(steps=500)

# Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int()

このコードの出力を使用して、時系列とその将来の値の予測をプロットできます。

ax = y.plot(label='observed', figsize=(20, 15))
pred_uc.predicted_mean.plot(ax=ax, label='Forecast')
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.25)
ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')

plt.legend()
plt.show()

これで、生成した予測と関連する信頼区間の両方を使用して、時系列をさらに理解し、何を期待するかを予測できます。私たちの予測によると、時系列は着実なペースで増加し続けると予想されます。

私たちが将来をさらに予測するにつれて、私たちが自分たちの価値観に自信を失うのは当然です。これは、モデルによって生成された信頼区間に反映されています。信頼区間は、将来に向かってさらに進むにつれて大きくなります。

結論

このチュートリアルでは、Pythonで季節性ARIMAモデルを実装する方法について説明しました。 pandasおよびstatsmodelsライブラリを多用し、モデル診断の実行方法とCO2時系列の予測の作成方法を示しました。

あなたが試すことができる他のいくつかのことはここにあります：

動的予測の開始日を変更して、これが予測の全体的な品質にどのように影響するかを確認します。
パラメータの組み合わせを増やして、モデルの適合度を改善できるかどうかを確認してください。
別のメトリックを選択して、最適なモデルを選択してください。たとえば、AICメジャーを使用して最適なモデルを見つけましたが、代わりにサンプル外の平均二乗誤差を最適化することもできます。

さらに練習するために、別の時系列データセットをロードして、独自の予測を作成することもできます。