In [1]:
import pandas as pd
import matplotlib.pyplot as plt
import io
import requests
In [2]:
# 飛行機乗客数データ
url = "https://www.analyticsvidhya.com/wp-content/uploads/2016/02/AirPassengers.csv"
stream = requests.get(url).content
content = pd.read_csv(io.StringIO(stream.decode('utf-8')), index_col='Month', parse_dates=True, dtype='float')
In [3]:
content.head()
Out[3]:
In [4]:
passengers = content['#Passengers'][:120]
passengers_plot = content['#Passengers']
plt.plot(passengers_plot)
plt.show()
ARモデル¶
STEP1：単位根検定(ADF)による定常性確認、自己相関・偏自己相関の確認¶
In [11]:
from statsmodels import api as s_api
In [8]:
result = s_api.tsa.stattools.adfuller(passengers)
result
Out[8]:
In [10]:
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('critical Values:')
for key, value in result[4].items():
print('臨界値\t%s: %.3f' % (key, value))
s_api.graphics.tsa.plot_acf(passengers, lags=35) # 自己相関表示
s_api.graphics.tsa.plot_pacf(passengers, lags=35) # 偏自己相関表示
plt.show()
- Autocorrelation図より、ラグ14まで自己相関係数が有意
- Partial Autocorrelation図より、ラグ1に大きな正の相関が、ラグ13に負の相関がある
- 検定結果も帰無仮説が棄却されないので、非定常性critical Values
STEP2：残差に対する自己相関分析、偏自己相関分析¶
In [12]:
from statsmodels import tsa
import warnings
warnings.filterwarnings('ignore')
In [29]:
# ARのorder選択
for i in range(10):
model = tsa.arima.model.ARIMA(passengers, order=(i+1, 0, 0)).fit()
print('order ',i,"'s AIC: ", model.aic)
- 上記より、AR(9)で実行
In [31]:
model = tsa.arima.model.ARIMA(passengers, order=(9, 0, 0)).fit()
In [33]:
# 残差の確認
resid = model.resid
resid[:5]
Out[33]:
残差の自己相関と偏自己相関を可視化¶
In [34]:
fig = plt.figure(figsize=(5,8))
ax1 = fig.add_subplot(211)
fig = s_api.graphics.tsa.plot_acf(resid, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = s_api.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
- 上記のように、非定常性のデータをそのままARモデルに適用すると、相関が残差に現れる。
- -> 逆に言えば、残差に相関が残っていることは、解析手法が元データに残っている非定常性を処理できていないことを示している。
STEP3：ARモデルの予測値を実測値をプロット¶
In [36]:
pred = model.predict('1955-01-01', '1958-12-01')
plt.figure(figsize=(6,5))
plt.plot(passengers[70:120],'--')
plt.plot(pred, "k")
plt.xticks(rotation=45)
pred = model.predict('1958-01-01', '1965-12-01')
plt.figure(figsize=(4,5))
plt.plot(passengers[40:],'--')
plt.plot(pred, "k")
plt.xticks(rotation=45)
plt.show()
- 点線：予測値、棒線：実測値の上グラフを見ると、予測値が右にずれていおり、見せかけの回帰が発生していることがわかる
- -> 予測できているように見えるが、単なる見せかけの回帰
- 下グラフより全く予測できていないことがわかる
MAモデル¶
STEP1：サンプルデータに対するADF検定と、自己相関分析・偏自己相関分析¶
In [ ]:
# ARモデルでの結果を参照
STEP2：残差に対する自己相関分析、偏自己相関分析¶
In [13]:
# ARのorder選択
for i in range(10):
model = tsa.arima.model.ARIMA(passengers, order=(0, 0, i+1)).fit()
print('order ',i,"'s AIC: ", model.aic, "'s BIC: ", model.bic)
- 上記より、MA(6)で実行
In [14]:
model = tsa.arima.model.ARIMA(passengers, order=(0, 0, 6)).fit()
In [15]:
resid = model.resid
resid
Out[15]:
In [16]:
fig = plt.figure(figsize=(5,8))
ax1 = fig.add_subplot(211)
fig = s_api.graphics.tsa.plot_acf(resid, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = s_api.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
- グラフより、残差のデータ同士はまだ相関が残っているのがわかる
- MA回帰条件は定常性を満たす必要があるため、非定常性のデータにそのままMAモデルを適用すると、相関を持った残差が表れる
STEP3：MAモデルによる予測値と実測値をプロット¶
In [19]:
pred = model.predict('1955-01-01', '1958-12-01')
plt.figure(figsize=(6,5))
plt.plot(passengers[70:120],'--', label='real value')
plt.plot(pred, "k", label='predicted value')
plt.xticks(rotation=45)
pred = model.predict('1958-01-01', '1965-12-01')
plt.figure(figsize=(6,5))
plt.plot(passengers[40:],'--', label='real value')
plt.plot(pred, "k", label='predicted value')
plt.xticks(rotation=45)
plt.legend()
plt.show()
- 予測にトレンドも周期性も全く反映されていない
ARMAモデル¶
STEP1：サンプルデータに対するADF検定と、自己相関分析・偏自己相関分析¶
In [ ]:
# ARモデルでの結果参照
STEP2：残差に対する自己相関分析、偏自己相関分析¶
In [23]:
# モデルの選定
sellected_order = tsa.stattools.arma_order_select_ic(passengers, max_ar=10, max_ma=10, ic='aic')
In [24]:
# 選定結果
sellected_order.aic_min_order
Out[24]:
In [25]:
model = tsa.arima.model.ARIMA(passengers, order=(8, 0, 9)).fit()
In [26]:
resid = model.resid
resid
Out[26]:
In [27]:
fig = plt.figure(figsize=(5,8))
ax1 = fig.add_subplot(211)
fig = s_api.graphics.tsa.plot_acf(resid, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = s_api.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
- 残差に12周期ごとの自己相関がのこっている
STEP3：ARMAモデルによる予測値と実測値をプロット¶
In [28]:
pred = model.predict('1955-01-01', '1958-12-01')
plt.figure(figsize=(6,5))
plt.plot(passengers[70:120],'--', label='real value')
plt.plot(pred, "k", label='predicted value')
plt.xticks(rotation=45)
pred = model.predict('1958-01-01', '1965-12-01')
plt.figure(figsize=(6,5))
plt.plot(passengers[40:],'--', label='real value')
plt.plot(pred, "k", label='predicted value')
plt.xticks(rotation=45)
plt.legend()
plt.show()
- 周期性は捉えられるようになってきている
- トレンドが捉えられていない
ARIMAモデル¶
STEP1省略¶
STEP2：残差に対する自己相関分析、偏自己相関分析¶
- ARIMAモデルのパラメータ探索は、まずデータを差分（微分）してから、arma_order_select_ic()で探索する
In [36]:
import numpy as np
In [37]:
# ARIMAのorder選択
# まずデータを差分して定常性にする
diffed_passengers = np.diff(passengers)
In [39]:
diffed_passengers
Out[39]:
In [40]:
# モデルの選定
sellected_order = tsa.stattools.arma_order_select_ic(diffed_passengers, max_ar=10, max_ma=10, ic='aic')
In [41]:
# 選定結果
sellected_order.aic_min_order
Out[41]:
In [42]:
# モデル
model = tsa.arima.model.ARIMA(passengers, order=(10, 1, 4)).fit()
In [43]:
resid = model.resid
resid
Out[43]:
In [44]:
fig = plt.figure(figsize=(5,8))
ax1 = fig.add_subplot(211)
fig = s_api.graphics.tsa.plot_acf(resid, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = s_api.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
STEP3：ARIMAモデルによる予測値と実測値をプロット¶
In [45]:
pred = model.predict('1955-01-01', '1958-12-01')
plt.figure(figsize=(6,5))
plt.plot(passengers[70:120],'--', label='real value')
plt.plot(pred, "k", label='predicted value')
plt.xticks(rotation=45)
pred = model.predict('1958-01-01', '1965-12-01')
plt.figure(figsize=(6,5))
plt.plot(passengers[40:],'--', label='real value')
plt.plot(pred, "k", label='predicted value')
plt.xticks(rotation=45)
plt.legend()
plt.show()
- 過去データはかなり一致してきているが、予測に関してはARMAとあまり変わらない？
In [47]:
# 残差のプロット
plt.plot(resid)
Out[47]:
SARIMAモデル¶
STEP1省略¶
STEP2：残差に対する自己相関分析、偏自己相関分析¶
- 今までの残差の自己相関から、12ごとの周期性があったので、seasonalに12を設定
- seasonalの(p,d,q)はとりあえず1で設定
In [59]:
# モデル
model = tsa.statespace.sarimax.SARIMAX(passengers, order=(10,1,4), seasonal_order=(1,1,1,12)).fit()
In [60]:
resid = model.resid
resid
Out[60]:
In [61]:
fig = plt.figure(figsize=(5,8))
ax1 = fig.add_subplot(211)
fig = s_api.graphics.tsa.plot_acf(resid, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = s_api.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
STEP3：SARIMAモデルによる予測値と実測値をプロット¶
In [62]:
pred = model.predict('1955-01-01', '1958-12-01')
plt.figure(figsize=(6,5))
plt.plot(passengers[70:120],'--', label='real value')
plt.plot(pred, "k", label='predicted value')
plt.xticks(rotation=45)
pred = model.predict('1958-01-01', '1965-12-01')
plt.figure(figsize=(6,5))
plt.plot(passengers[40:],'--', label='real value')
plt.plot(pred, "k", label='predicted value')
plt.xticks(rotation=45)
plt.legend()
plt.show()
- トレンドも考慮されるようになった
In [63]:
# 残差のプロット
plt.plot(resid)
Out[63]:
In [ ]: