モモノキ＆ナノネと一緒にPythonで時系列分析を覚えよう（9）
時系列分析の実践練習（SARIMAモデルで未来予測に挑戦しよう）
今回はSARIMAモデル（季節自己回帰和分移動平均）を実際に使って時系列分析をやってみるよ。
第9回まできて、やっと実践だね。
調べては試すの繰り返で、なかなか先に進まない。難しいね、時系列分析。
難しいけどガンバルよ。
じゃあ始めよう！ナノネ、飛行機乗客数のデータをPandasに読み込み込むまでのコードを準備して。これまでと同じでOKだよ。
# 飛行機乗客数の時系列データ（前回使用コード抜粋 -> 次回SARIAモデル実践に利用）
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot') # グラフのデザイン変更（お好みで利用）
plt.xkcd() # グラフのデザイン変更（お好みで利用、日本語は未対応）
import statsmodels.api as sm # version 0.8.0以上
# CSVファイル読み込み（事前にダウンロードしたCSVファイルを利用）
df = pd.read_csv('AirPassengers.csv')
# Pandas.Seriesにデータを格納（データに乗客数、インデックスは日付）
passengers = pd.Series(df['#Passengers'], dtype='float') # ①
passengers.index = pd.to_datetime(df['Month']) # ②
passengers.plot()
これね、何回も呼び出してるから慣れた。
OK、早速だけどSARIMAモデルを使ってみよう。
SARIMAモデルは、一派的にSARIMA(p, d, q)(P, D, Q)[s]と表記されることが多いよ。
最初の3つのパラメータp,d,qは、ARIMAに適用するAR,差分,MAの次数。次のP,D,Qは同じ対応順だけど、季節調整に適用する次数。そして最後sは季節調整に適用する周期を指定するよ。
PythonだとStatsModelsの『SARIMAX』でSARIMAモデルが使えるよ。最後にXがついてるけど気にしないで。Xの意味は外部変数もモデルに適用できることみたいだけど、今回は使わないよ（実のところ使い方もわからないけど）。
今回StatsModelsは上でsmという名前でインポート済みだから、sm.tsa.SARIMAX(複数パラメータ)とすればOK。
モデルをつくるときのコードや引数のパラメータは、どんな風に書けばいいの？
モデルを作るときはこんな感じで指定してみて。
モデル変数名 = sm.tsa.SARIMAX(passengers, order=(p, d, q), seasonal_order=(P, D, Q, s), enforce_stationarity = False, enforce_invertibility = False).fit()
まず、第1引数に時系列データを指定してね。前回も説明したけどSARIMAはトレンドや季節性を含むデータにも適用できるよ。
orderにはARIMAに適用する3つの次数を指定するよ。p,d,qの順番でタプル形式で値を渡してね。繰り返しになるけど、pはAR（自己回帰）、qはMA（移動平均）、dが差分の次数にだよ。
seasonal_orderはパラメータ4つをタプルで渡してね。順番はP,D,Q,sで前半の3つorderと一緒、最後4つ目は周期。季節調整に適用する各々の次数を指定してね。
そのあとの名前付き引数enforce_stationarity（定常性強制）とenforce_invertibility（可逆性強制）は、どっちもFalseを指定しておいてね。勉強不足で説明できなけど、計算過程のエラーを回避するのに必要みたい。
最後に.fit()とつなげて、作成されたモデルを変数で受け取ればOK。
コードは書けそうだけど、パラメータはどうやって決めればいいの？季節調整のsは月次データだし、前回周期もチェックしているから12でよさそうだけど。p,d,qやP,D,Qはどうやって求めるの？
。。。
。。。？！
実はそこが難しいとこなんだよ。
ARの次数p、MAの次数qは自動推定みたいな機能があったりすけど、確実ではないみたい。データをよく観察した上である程度当たりをつけて、最後は力技で組み合わせパターンの総当たりチェック。その中からベストな次数を選ぶとういう手法がよく使われているよ。
時系列データの内容にもよるけど、目安の次数範囲はあるみたい。ARMA(p,q)は0～3程度、差分は0か1（多くても2？）、P,Qはあまり大きな次数となるケースが少ないなども参考になるよ。
なんか、パラメータは試行錯誤で手探りっぽい選び方だね。
うん。まずはデータをチェックしながらパラメータの当たりをつけて、一度モデルを作ってみよう。
了解。データは何からチェックすればいい？
ADF検定というデータの定常性チェックを試してみよう。ADF検定は拡張ディッキー-フラー検定（Augmented Dickey-Fuller test, ADF test）のことだよ。時系列のサンプルデータが、単位根過程（非定常過程）であるかどうかを調べたい時に利用するよ。
単位根過程とういうのは非定常だけど、差分をとると（弱）定常になるようなデータだよ。ADF検定は帰無仮説『単位根過程（非定常過程）である』を棄却するこで、そのデータが定常過程とみなせるよ。
# ADF検定（原型列で確認だけ）
res_ctt = sm.tsa.stattools.adfuller(passengers, regression="ctt") # トレンド項あり（２次）、定数項あり
res_ct = sm.tsa.stattools.adfuller(passengers, regression="ct") # トレンド項あり（１次）、定数項あり
res_c = sm.tsa.stattools.adfuller(passengers, regression="c") # トレンド項なし、定数項あり
res_nc = sm.tsa.stattools.adfuller(passengers, regression="nc") # トレンド項なし、定数項なし
print(res_ctt)
print(res_ct)
print(res_c)
print(res_nc)
なんか一杯データでたけど、何か分かったの？
どの条件でもp-valueが高く帰無仮説を棄却できないから、定常とは言えないね。このあと、原型列に何らかの変換をして定常になる条件を確認していく流れみたいだけど、まだ勉強不足なので今回はここまでして、また時間があるときに詳しく調べてみよう。
了解、使い方はよくわからいなけど試してみたのね。
うん。練習、練習。。。
次はARMA(p,q)の次数を自動推定関数をつかって求めてみよう。原型列データがトレンドを持っているので、差分をとってから利用することにしょう。
passengers_diff = passengers - passengers.shift() # 差分(1階差) Pandasのdiff()でpassengers.diff()としてもOK
passengers_diff = passengers_diff.dropna() # 1個できるNaNデータは捨てる
passengers_diff.plot()
差分データ準備した。次は？
StatasModelsにあるsm.tsa.arma_order_select_ic()でARMAのパラメータ推定関数が使えるよ。引数は第1引数に時系列データを指定して、今回その他はic="aic", trend="nc"としてみよう。aicはモデルを評価する手法のひとつでよく使わているよ。aicの値が小さいほど良いモデルとみなせて、結果はaicがベストなp,qの組み合わせを出力してくるよ。
import warnings
warnings.filterwarnings('ignore') # 計算警告を非表示
# 自動ARMAパラメータ推定関数
res_selection = sm.tsa.arma_order_select_ic(passengers_diff, ic='aic', trend='nc')
res_selection
右下に'aic_min_order': (3, 2)とでてきた。
ベストaic（最小値）の組み合わせは、p=3,q=2だったみたい。ちなみに表示される行列データは縦がp、横がqの次数だよ。
ナノネ、復習でコレログラムもプロットしてみよう。
# 差分系列のコレログラム
fig = plt.figure(figsize=(8, 8))
# 自己相関(ACF)のグラフ
ax1 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(passengers_diff, lags=40, ax=ax1) #ACF計算とグラフ自動作成
# 偏自己相関(PACF)のグラフ
ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(passengers_diff, lags=40, ax=ax2) #PACF計算とグラフ自動作成
plt.tight_layout() # グラフ間スキマ調整
コレログラム描けた。まだ見かたはイマイチよくわからないけどだけど、周期はハッキリ読み取れる。
OK、コレログラムからも次数のメボシをつけられるみたいだけど難しいね。詳しい見かたは今後の課題にしよう。
自動推定でp=3,q=2。季節調整の周期sは12として、d,P,D,Q,の次数はどうする？
とりあえずだけど、差分は1回でd,Dを1、P,Qもあまり大きな次数をとらない想定で1としておこう。次数の最適化はあとでもう一度総当たりチェックね。
パラメータはこれで全部そろったから、上の条件で一度モデルを作ってみる。
sm.tsa.SARIMAX()を使って、第1引数に時系列データ（今回は原型列を使用）を指定。パラメータ次数はorder=(p, d, q)、seasonal_oreder=(P, D, Q, s)の順番。enforce＊の二つはFalseで最後に.fit()。
こんな感じで合ってる？
# SRIMAモデル作成その１
sarimax = sm.tsa.SARIMAX(passengers,
order=(3, 1, 2),
seasonal_order=(1, 1, 1, 12),
enforce_stationarity = False,
enforce_invertibility = False
).fit()
OK、無事作れたね。できたモデルの残差をコレログラムで確認してみよう。残差(residual)はモデル変数名のあとに.residと取得できるよ。
sarimax_resid = sarimax.resid # モデルの残差成分
# モデル残差のコレログラム
fig = plt.figure(figsize=(8, 8))
# 自己相関(ACF)のグラフ
ax1 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(sarimax_resid, lags=40, ax=ax1) #ACF計算とグラフ自動作成
# 偏自己相関(PACF)のグラフ
ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(sarimax_resid, lags=40, ax=ax2) #PACF計算とグラフ自動作成
plt.tight_layout() # グラフ間スキマ調整
モデル残差のコレログラムを描けた。残差の自己相関は全体的に小さめだね。
残差に周期性とか自己相関がたくさん残ってるいると、予測モデルとしてはダメでモデルの再考が必要だよ。完璧ではないけど、残差の自己相関はこれでヨシとしよう。
試しに季節調整なしで、モデルを作ってコレログラムを描くとこんな感じ。seasonal_order=(0, 0, 0, 0)として、ARIMAモデルにしたイメージだよ。
# SRIMAモデル（試しに）季節調整なし
sarimax_noseasonal = sm.tsa.SARIMAX(passengers,
order=(3, 1, 2),
seasonal_order=(0, 0, 0, 0),
enforce_stationarity = False,
enforce_invertibility = False
).fit()
sarimax_noseasonal_resid = sarimax_noseasonal.resid # 残差成分
fig = plt.figure(figsize=(8, 8))
# 自己相関(ACF)のグラフ
ax1 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(sarimax_noseasonal_resid, lags=40, ax=ax1) #ACF計算とグラフ自動作成
# 偏自己相関(PACF)のグラフ
ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(sarimax_noseasonal_resid, lags=40, ax=ax2) #PACF計算とグラフ自動作成
plt.tight_layout() # グラフ間スキマ調整
やっぱり季節調整つけないと、残差に周期性が残っちゃうね。
aicの値はモデル変数名.aicで確認できるよ。
print(sarimax.aic) # 季節調整あり -> 914
print(sarimax_noseasonal.aic) # 季節調整なし -> 1352
他にもモデル変数名.summaryでもモデルの情報が確認できるよ。解説はできなけど、気になるところは調べてみてね。
sarimax.summary()
わかった、気になる項目は今度調べてみる。
作ったモデルで予測値を求めたり、グラフにプロットするにはどうすればいいの？
モデル変数名.predict()が使えば結果を取得できるよ。引数は開始と終了時点を順番に指定してね。ナノネ、まずは1951年1月から1960年12月の条件で結果を取得してグラフも描いてみて。
# モデルの値でプロット
sarimax_pred = sarimax.predict('1951-01', '1960-12')
plt.figure(figsize=(8, 4))
plt.plot(sarimax_pred, c="b")
なんかオリジナルと似たようなデータがプロットされた。オリジナルデータと重ね書きしてみる。
sarimax_pred = sarimax.predict('1951-01', '1960-12')
plt.figure(figsize=(8, 4))
plt.plot(passengers, label="original")
plt.plot(sarimax_pred, c="b", label="model-pred", alpha=0.7)
plt.legend(loc='best')
ズレはあるけど、ある程度いい感じにフィットしているね。
次はオリジナルデータにはない未来の時点を指定して予測をしてみよう。指定する期間はオリジナルとラップしている必要があるから注意ね。1949-01～1960-12までデータで最後が1960-12だから、予測期間を1960-12～1962-12として2年分試してみよう。
sarimax_pred = sarimax.predict('1960-12', '1962-12')
plt.figure(figsize=(8, 4))
plt.plot(passengers, label="original")
plt.plot(sarimax_pred, c="b", label="model-pred", alpha=0.7)
plt.legend(loc='best')
それっぽい予測をして来たけど、それなりに合ってるの？
正解データがないから分からないよ。モデル作るときに（訓練データとして）データを全部使っちゃたからね。次はテスト用に最後の2年位を残してモデルを作ろう。モデル作成後に残しておいたテストデータを使って予測、どれくらいあっているか検証してみよう。
その前に、パレメータ次数の最適化を一度やってみよう。
候補とする次数を総当たりチェックする手法を試してみよう。力技でループ処理を回して、AICの値からベストな組み合わせを探すだけ。処理は下みたいな簡単なコードでも実行できるけど、組み合わせパターンを増やすとすごく時間がかかるのが難点。
古い32ビットPCでスペック低いから、力技の処理は苦手だけど。。。次数のチェック範囲はあまり増やさないで一度やってみる。
# SARIMAパラメター最適化（総当たりチェック）
import warnings
warnings.filterwarnings('ignore') # 警告非表示（収束：ConvergenceWarning）
# パラメータ範囲
# order(p, d, q)
min_p = 1; max_p = 3 # min_pは1以上を指定しないとエラー
min_d = 0; max_d = 1
min_q = 0; max_q = 3
# seasonal_order(sp, sd, sq)
min_sp = 0; max_sp = 1
min_sd = 0; max_sd = 1
min_sq = 0; max_sq = 1
test_pattern = (max_p - min_p +1)*(max_q - min_q + 1)*(max_d - min_d + 1)*(max_sp - min_sp + 1)*(max_sq - min_sq + 1)*(max_sd - min_sd + 1)
print("pattern:", test_pattern)
sfq = 12 # seasonal_order周期パラメータ
ts = passengers # 時系列データ
test_results = pd.DataFrame(index=range(test_pattern), columns=["model_parameters", "aic"])
num = 0
for p in range(min_p, max_p + 1):
for d in range(min_d, max_d + 1):
for q in range(min_q, max_q + 1):
for sp in range(min_sp, max_sp + 1):
for sd in range(min_sd, max_sd + 1):
for sq in range(min_sq, max_sq + 1):
sarima = sm.tsa.SARIMAX(
ts, order=(p, d, q),
seasonal_order=(sp, sd, sq, sfq),
enforce_stationarity = False,
enforce_invertibility = False
).fit()
test_results.iloc[num]["model_parameters"] = "order=(" + str(p) + ","+ str(d) + ","+ str(q) + "), seasonal_order=("+ str(sp) + ","+ str(sd) + "," + str(sq) + ")"
test_results.iloc[num]["aic"] = sarima.aic
print(num,'/', test_pattern-1, test_results.iloc[num]["model_parameters"], test_results.iloc[num]["aic"] )
num = num + 1
# 結果（最小AiC）
print("best[aic] parameter ********")
print(test_results[test_results.aic == min(test_results.aic)])
時間かかったけど、pは1～3、qは0～3、d、P,D,Qは0～1の範囲でチェックしたみた。
test_results.sort_values(by='aic').head(10) # aic top10
チェックした範囲では、order=(3,1,3), seasonal_order=(0,1,1)の組み合わせがAIC最小でベストみたい。AICの値も最初に作ったモデルより少しだけ小さくなって914が898に。
# SRIMAモデル（ちょっとだけパラメータ最適化・総当たりaicベストを適用）
sarimax_optimization = sm.tsa.SARIMAX(passengers,
order=(3, 1, 3),
seasonal_order=(0, 1, 1, 12),
enforce_stationarity = False,
enforce_invertibility = False
).fit()
sarimax_optimization_resid = sarimax_optimization.resid # 残差成分
fig = plt.figure(figsize=(8, 8))
# 自己相関(ACF)のグラフ
ax1 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(sarimax_optimization_resid, lags=40, ax=ax1) #ACF計算とグラフ自動作成
# 偏自己相関(PACF)のグラフ
ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(sarimax_optimization_resid, lags=40, ax=ax2) #PACF計算とグラフ自動作成
plt.tight_layout() # グラフ間スキマ調整
さっきの総当たりチェックでベストだったパレメータで、もう一度モデルを作ってみた。
残差もそこそこ小さいから、最適化の試し版としてはヨシとししょう。PCのスペックや時間に余裕ができたら、違う組み合わせも試してみてね。ただし、AICの値だけでなく、しっかりデータを分析して妥当な次数か見極めることも大事だよ。次数は無駄にあげると説明変数が増えてモデルも複雑化するので推奨されないよ。
OK、モデルはできるだけシンプルを心がける。次はテスト用データを残してモデル作成、予測結果の検証だね。
モデル作成用（訓練）データとテストデータを準備しよう。オリジナルが1949-01から1960-12なので、モデル作成用を1949-01から1958-12として、テスト用は1959-01～1960-12までの2年分としよう。
# オリジナル passengers 1949-01 ～ 1960-12
passengers_train = passengers['1949-01':'1958-12'] # モデル作成用データ（訓練）
print(passengers_train.head())
print(passengers_train.tail())
passengers_test = passengers['1958-01':'1960-12'] # テスト用データ2年分
print(passengers_test.head())
print(passengers_test.tail())
期間で分けて2つの変数に入れた。
OKだよ。モデルのパラメータは、さっき最適化で選んだ値を使うことにしよう。
# SRIMAモデル（テストデータを2年分除いてモデル作成）
sarimax_train = sm.tsa.SARIMAX(passengers_train,
order=(3, 1, 3),
seasonal_order=(0, 1, 1, 12),
enforce_stationarity = False,
enforce_invertibility = False
).fit()
テストデータ(passengers_train)を除いたデータでモデルを作った。
作ったモデルでテスト期間の値を予測をしてみよう。既知のデータだけど、モデルには教えてないから未来予測になるよね。
sarimax_train_pred = sarimax_train.predict('1958-12', '1960-12') # テストデータ2年分予測
予測まで完了。予測と正解を重ね書きでプロットしてみる。
plt.figure(figsize=(8, 4))
plt.plot(passengers, label="original")
plt.plot(sarimax_train_pred, c="b", label="model-pred", alpha=0.7)
plt.legend(loc='best')
予測のほうが山が少し小さいけど、それらしい形はしてるみたい。テストデータを1年分にして予測期間を短くすれば、もちょうっとフィットする？
どうかな？ナノネ、コピペで試してみたら。
# オリジナル passengers 1949-01 ～ 1960-12
passengers_train2 = passengers['1949-01':'1959-12'] # モデル作成用データ（訓練）1年テスト用残し
passengers_test2 = passengers['1960-01':'1960-12'] # テスト用データ1年分
# SRIMAモデル（テストデータ1年を除いてモデル作成）
sarimax_train = sm.tsa.SARIMAX(passengers_train2,
order=(3, 1, 3),
seasonal_order=(0, 1, 1, 12),
enforce_stationarity = False,
enforce_invertibility = False
).fit()
sarimax_train2_pred = sarimax_train.predict('1959-12', '1960-12') # テストデータ1年分予測
plt.figure(figsize=(8, 4))
plt.plot(passengers, label="original")
plt.plot(sarimax_train2_pred, c="b", label="model-pred", alpha=0.7)
plt.legend(loc='best')
plt.figure(figsize=(8, 4))
plt.plot(passengers_test2, label="actual") # 正解
plt.plot(sarimax_train2_pred, c="b", label="predict", alpha=0.7) # 予測
plt.legend(loc='best')
テストデータを1年分にして予測期間短くしたら、さっきよりもフィット感上がったみたい。テスト期間だけ拡大プロットもしてみた。
精度は別として、やっと予測まで実行できたね。
なんとか無事できた。
今回はここまで。完了でもいいんだけど、次回もう少し続きをやってみよう。
まだ続きがあるのね、ラジャー。
0 件のコメント :
コメントを投稿