Pythonによる時系列分析の基礎¶
https://logics-of-blue.com/python-time-series-analysis/¶
In [2]:
# 基本のライブラリを読み込む
import numpy as np
import pandas as pd
from scipy import stats
# グラフ描画
from matplotlib import pylab as plt
import seaborn as sns
%matplotlib inline
sns.set()
# グラフを横長にする
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
# 統計モデル
import statsmodels.api as sm
# ワーニングを表示させない
import warnings
warnings.filterwarnings('ignore')
時系列データの読み込み¶
In [3]:
# 普通にデータを読み込む
# https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/AirPassengers.html
dataNormal = pd.read_csv('AirPassengers.csv')
dataNormal.head()
Out[3]:
In [4]:
# 日付形式で読み込む（dtype=floatで読み込まないと、あとでARIMAモデル推定時にエラーとなる）
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
data_old_method = pd.read_csv('AirPassengers.csv', index_col='Month', date_parser=dateparse, dtype='float')
data_old_method.head()
Out[4]:
In [5]:
# 2018年7月7日追記
# 以下のコードの方が読み込みは簡単
data = pd.read_csv('AirPassengers.csv',
index_col='Month',
parse_dates=True,
dtype='float')
data.head()
Out[5]:
In [6]:
# 日付形式にする
ts = data['#Passengers']
ts.head()
Out[6]:
時系列データの取り扱い¶
In [7]:
# プロット
plt.plot(ts)
Out[7]:
In [8]:
# データの取得方法その1
ts['1949-01-01']
Out[8]:
In [9]:
# データの取得方法その2
from datetime import datetime
ts[datetime(1949,1, 1)]
Out[9]:
In [10]:
# 1949年のデータをすべて取ってくる
ts['1949']
Out[10]:
In [11]:
# シフト
ts.shift().head()
Out[11]:
In [12]:
# 差分系列
diff = ts - ts.shift()
diff.head()
Out[12]:
In [13]:
# 2018年7月7日追記
# こちらの方が簡単
ts.diff().head()
Out[13]:
In [14]:
# 対数差分系列
logDiff = np.log(ts) - np.log(ts.shift())
# NaNを取り除てから表示
logDiff.dropna().head()
Out[14]:
自己相関係数の推定¶
In [15]:
# 自己相関を求める
ts_acf = sm.tsa.stattools.acf(ts, nlags=40)
ts_acf
Out[15]:
In [16]:
# 偏自己相関
ts_pacf = sm.tsa.stattools.pacf(ts, nlags=40, method='ols')
ts_pacf
Out[16]:
In [17]:
# 自己相関のグラフ
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(ts, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(ts, lags=40, ax=ax2)
ARIMAモデルの推定¶
In [18]:
# たぶん和分過程なので、差分をとる
diff = ts - ts.shift()
diff = diff.dropna()
diff.head()
Out[18]:
In [19]:
# 差分系列のグラフ
plt.plot(diff)
Out[19]:
In [20]:
# 差分系列への自動ARMA推定関数の実行
resDiff = sm.tsa.arma_order_select_ic(diff, ic='aic', trend='nc')
resDiff
Out[20]:
In [21]:
# P-3, q=2 が最善となったので、それをモデル化
from statsmodels.tsa.arima_model import ARIMA
ARIMA_3_1_2 = ARIMA(ts, order=(3, 1, 2)).fit(dist=False)
ARIMA_3_1_2.params
Out[21]:
In [22]:
# 残差のチェック
# SARIMAじゃないので、周期性が残ってしまっている。。。
resid = ARIMA_3_1_2.resid
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
SARIMAモデルの推定¶
In [23]:
# SARIMAモデルを「決め打ち」で推定する
import statsmodels.api as sm
SARIMA_3_1_2_111 = sm.tsa.SARIMAX(ts, order=(3,1,2), seasonal_order=(1,1,1,12)).fit(method='bfgs', maxiter=300)
print(SARIMA_3_1_2_111.summary())
In [24]:
# 残差のチェック
residSARIMA = SARIMA_3_1_2_111.resid
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(residSARIMA.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(residSARIMA, lags=40, ax=ax2)
In [25]:
# 予測
pred = SARIMA_3_1_2_111.predict('1960-01-01', '1961-12-01')
print(pred)
In [26]:
# 実データと予測結果の図示
plt.plot(ts)
plt.plot(pred, "r")
Out[26]:
総当たり法によるSARIMAモデル次数の決定¶
In [27]:
# 総当たりで、AICが最小となるSARIMAの次数を探す
max_p = 3
max_q = 3
max_d = 1
max_sp = 1
max_sq = 1
max_sd = 1
pattern = max_p*(max_q + 1)*(max_d + 1)*(max_sp + 1)*(max_sq + 1)*(max_sd + 1)
modelSelection = pd.DataFrame(index=range(pattern), columns=["model", "aic"])
pattern
Out[27]:
In [28]:
# 自動SARIMA選択
num = 0
for p in range(1, max_p + 1):
for d in range(0, max_d + 1):
for q in range(0, max_q + 1):
for sp in range(0, max_sp + 1):
for sd in range(0, max_sd + 1):
for sq in range(0, max_sq + 1):
sarima = sm.tsa.SARIMAX(
ts, order=(p,d,q),
seasonal_order=(sp,sd,sq,12),
enforce_stationarity = False,
enforce_invertibility = False
).fit(method='bfgs', maxiter=300, disp=False)
modelSelection.ix[num]["model"] = "order=(" + str(p) + ","+ str(d) + ","+ str(q) + "), season=("+ str(sp) + ","+ str(sd) + "," + str(sq) + ")"
modelSelection.ix[num]["aic"] = sarima.aic
num = num + 1
In [29]:
modelSelection
Out[29]:
In [30]:
# AIC最小モデル
modelSelection[modelSelection.aic == min(modelSelection.aic)]
Out[30]:
In [31]:
# 参考：次数がすべて０だとエラーになる
# sarima = sm.tsa.SARIMAX(ts, order=(0,0,0), seasonal_order=(0,0,0,12), enforce_stationarity = False).fit()
In [32]:
bestSARIMA = sm.tsa.SARIMAX(ts, order=(3,1,3), seasonal_order=(0,1,1,12), enforce_stationarity = False, enforce_invertibility = False).fit()
In [33]:
print(bestSARIMA.summary())
In [34]:
# 残差のチェック
residSARIMA = bestSARIMA.resid
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(residSARIMA, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(residSARIMA, lags=40, ax=ax2)
In [35]:
# 予測
bestPred = bestSARIMA.predict('1960-01-01', '1961-12-01')
# 実データと予測結果の図示
plt.plot(ts)
plt.plot(bestPred, "r")
Out[35]:
