Python

時系列解析に出てくるARIMAモデルとSARIMAモデルを徹底解説

はじめに

今回は、様々な時系列データの解析手法のうち、ARIMAモデルとSARIMAモデルを紹介します。
ARIMAモデルとは、autoregressive integrated moving averageの略で、自己回帰モデル(ARモデル)、移動平均モデル(MAモデル)、和分モデル(Iモデル)の3モデルを組み合わせたモデルです。ARIMAモデルは、自己回帰和分移動平均モデルとも呼ばれます。
SARIMAモデルとは、Seasonal AutoRegressive Integrated Moving Averageの略で、ARIMAモデルに「季節的な周期パターン」を加えたモデルです。
ARIMAモデルとSARIMAモデルは、非定常過程のデータに対して適用できるのが大きな特徴です。
ARIMAモデルやSARIMAモデルを理解するには、MAモデル、ARモデル、ARMAモデルの理解が必要となりますので、事前にそれらを簡単に見ておきます。

MAモデル

MAモデルとは移動平均(Moving Average)モデルのことです。1次MAモデルを例に取ってMAモデルを考えてみましょう。
1次MAモデルは
yt = c + εt + θ1εt−1
の式で表されるモデルです。ここでεtは時点tにおけるホワイトノイズです。ホワイトノイズとは、時点tによって発生するランダムな数と思ってください。
1 次MAモデルでは、1時点前までのホワイトノイズを用い、n次MAモデルであれば、n時点前までのホワイトノイズを用います。

1次MAモデルでは、yt-1=c + εt-1 + θ1εt−2が成り立ち、ytとyt-1でεt-1という共通項を持つことが分かります。これは、ytとyt-1の間に相関が生じることを意味しており、1次自己相関を持つと言います。
では、1次MAモデルの例をいくつか見てみましょう(図1)。

import numpy as np
import matplotlib.pyplot as plt

theta = [-0.9, -0.5, 0.5, 0.9]
sigma = [0.5, 2] 
c = 0
sample = 100


fig = plt.figure(figsize=(15, 5))
fig.subplots_adjust(hspace=0.5, wspace=0.2)

for k in range(len(sigma)):
    for j in range(len(theta)):

        noize = np.random.normal(loc=0, scale=sigma[k], size=sample+1)
        y0 = 0
        y = np.zeros(sample+1)
        y[0] = y0
        
        for i in range(sample):
            y[i+1] = c+noize[i+1]+theta[j] * noize[i]
        y= y
        
        ax = fig.add_subplot(2, 4, 1+j+4*k)
        ax.set_title("c = 0, θ = {}, ε = {}".format(theta[j], sigma[k]))
        ax.plot([int(i) for i in range(1, sample+1)],y)
plt.show()

上段がε=0.5、下段がε=2なので、下段の方が振幅が大きくなっています。そして、左から右に向かってθ1が大きくなっています。符号が負であるとギザギザしており、符号が正であると滑らかになっています。これは符号が正の場合、1次の正の自己相関が生じるので、ytとyt-1は同じ方向に動く傾向となります。そして符号が負の場合、1次の負の自己相関が生じるので、ytとyt-1は逆方向に動く傾向となるためです。

また、ここでは詳しく述べませんが、MAモデルは期待値と自己共分散が時点tに依存していないので、MAモデルはいかなるパラメータを設定しても、定常過程であるという特徴があります。

3. ARモデル

ARモデルとは自己回帰(AutoRegressive)モデルのことです。1次ARモデルを例に取ってARモデルを考えてみましょう。

1次ARモデルは

yt = c + Φ1yt−1 + εt

の式で表されるモデルで、1時点前のデータとホワイトノイズから成り立っています。ytの式にyt-1の項が含まれているので、ytとyt-1が相関を持つことが明らかに分かります。

では、1次ARモデルの例をいくつか見てみましょう(図2)。

import numpy as np
import matplotlib.pyplot as plt
 
phi = [-1.1, -0.8, -0.2, 0.2, 0.8, 1.1]
sigma = [0.5, 2]
center = 0
sample = 100
fig = plt.figure(figsize=(16, 6))
fig.subplots_adjust(hspace=0.5, wspace=0.4)

for k in range(len(sigma)):
    for j in range(len(phi)):

        noize = np.random.normal(loc=0, scale=sigma[k], size=sample+1)
        y0 = 0
        y = np.zeros(sample+1)
        y[0] = y0
        
        for i in range(sample):
            y[i+1] = center + phi[j] * y[i] + noize[i]
        y= y
        
        ax = fig.add_subplot(2, 6, 1+j+6*k)
        ax.set_title("c = 0,  Φ= {}, ε = {}".format(phi[j], sigma[k]))
        ax.plot([int(i) for i in range(1, sample+1)],y)
plt.show()

上段がε=0.5、下段がε=2なので、1次MAモデルと同じく下段の方が振幅が大きくなっています。そして、両端のグラフと2~5列目のグラフで明らかに違うことが分かります。2~5列目のグラフは、1次MAモデルのように、上下方向に万遍なくばらついており、定常性を有していることがうかがえます。それに対して、両端のグラフは明らかな傾向が見て取れ、定常性を有していない、つまり、非定常なモデルとなっています。

1次ARモデルの場合、|Φ1|<1のとき定常過程となり、|Φ1|≧1のとき非定常となるという特徴がありますが、グラフを見れば一目瞭然です。

ARMAモデル

次は、ARMAモデルを見ていきます。これは自己回帰移動平均(AutoRegressive Moving Average)モデルのことで、名前の通り、これまで見てきたARモデルとMAモデルを組み合わせたモデルです。

(p, q)次ARMAモデルを一般式で表すと、

yt = c + Φ1yt−1 + … + Φpyt−p + εt + θ1εt−1 + … +θqεt−q 

で定義できます。

ARモデルにMAモデルを加えたことでモデル式の柔軟性が高まり、定常過程の時系列データに対しては、非常に強力な説明力・予測力を持っているのが大きな特徴です。

MAモデルは常に定常性を有していますので、ARMAモデルが定常過程かどうかはARモデルが定常過程かどうかに依存します。

では、(1,1)次ARモデルの例をいくつか見てみましょう(図3)。なお、Φ以外のパラメータは、θ1=0.4、ε=2に固定しています。

import numpy as np
import matplotlib.pyplot as plt

phi = [-1.1, -0.8, -0.2, 0.2, 0.8, 1.1]
theta = 0.4
sigma = 2
center = 50
sample = 100
fig = plt.figure(figsize=(16, 3))
fig.subplots_adjust(hspace=0.5, wspace=0.4)

for j in range(len(phi)):
    noize = np.random.normal(loc=0, scale=sigma, size=sample+1)
    y0 = center
    y = np.zeros(sample+1)
    y[0] = y0
    
    for k in range(sample):
        y[k+1] = center + phi[j] * y[k] + noize[k+1] + theta*noize[k]
    y= y

    ax = fig.add_subplot(1, 6, 1+j)
    ax.set_title("Φ= {}".format(phi[j]))
    ax.plot([int(i) for i in range(1, sample+1)],y)
plt.show()

MAモデルは常に定常のため、ARモデルと同じく、|Φ1|<1のとき定常過程となり|Φ1|≧1のとき非定常となっていることが分かります。ARMAモデルはMAモデルとARモデルを組み合わせているため、定常過程のデータに対しては、かなり精度よく予測が可能です。そこで、図3の左から2番目のデータを使って、各パラメータを予測してみましょう。

#データの作成

import numpy as np
import matplotlib.pyplot as plt

phi = -0.8
theta = 0.4
sigma = 2
center = 50
sample = 100

noize = np.random.normal(loc=0, scale=sigma, size=sample+1)
y0 = center
y = np.zeros(sample+1)
y[0] = y0
    
for k in range(sample):
    y[k+1] = center + phi * y[k] + noize[k+1] + theta*noize[k]
y= y

#パラメータの予測
#次数の推定

import statsmodels.api as sm

print(sm.tsa.arma_order_select_ic(y, max_ar=2, max_ma=2, ic='aic'))
{'aic':             0           1           2
0  579.283912  524.044364  508.941711
1  472.735539  467.073436  468.476308
2  469.536316  468.710843  469.834314, 'aic_min_order': (1, 1)}

まずは、ARMAモデルの次数をAICで推定します。今回の解析は、うまく(1,1)次と推定できました。

そこで、(1,1)次ARMAモデルを作成し、パラメータを推定します。

# ARMAモデルの作成と推定

model = sm.tsa.ARMA(y, order=(1, 1))
result = model.fit()

print(result.summary())
ARMA Model Results                              
========================================================================
Dep. Variable:                      y   No. Observations:            100
Model:                     ARMA(1, 1)   Log Likelihood          -229.537
Method:                       css-mle   S.D. of innovations        2.369
Date:                Fri, 03 Apr 2020   AIC                      467.073
Time:                        00:11:21   BIC                      477.494
Sample:                             0   HQIC                     471.291
                                                                              
========================================================================
           coef    std err         z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------
const    27.6869      0.162    170.739      0.000      27.369    28.005
ar.L1.y  -0.9846      0.021    -45.991      0.000      -1.027    -0.943
ma.L1.y  0.3547      0.109      3.250      0.002       0.141      0.569
                                    Roots                                    
========================================================================
             Real          Imaginary           Modulus         Frequency
------------------------------------------------------------------------
AR.1        -1.0157           +0.0000j            1.0157         0.5000
MA.1        -2.8190           +0.0000j            2.8190         0.5000
------------------------------------------------------------------------

データの実際のパラメータΦ1= -0.8, θ1=0.4, ε=2に対し、推定値はΦ1= -0.98, θ1=0.35, ε=2.369となり、当たらずといえども遠からずな結果が得られました。サンプル数が100だったので、もう少しデータを多くすれば、推定精度は上がるかもしてません。

以下のコードで、実測値と推定値のグラフを描けます。

result.plot_predict()

なお、図3の一番左または右の非定常のデータで、同じようにARMAモデルでパラメータの推定を行おうとすると、エラーとなります。ARMAモデルは、定常過程の時系列データしか解析できないことを覚えておいてください。

ARIMAモデルとは

定常性を有する時系列データに対しては、ARMAモデルがパラメータの推定に有効であることが分かりました。

しかし、定常過程でない、非定常なデータに対してはARMAモデルは無力です。そこで、非定常な時系列データに有効なモデルの一つである、ARIMAモデルを見ていきます。

世の中には、日経平均株価の推移など定常過程の性質を持たない、非定常なデータを得られることの方が多いので、それらのデータ解析に向いているモデルです。

これまで見てきた定常過程は、トレンドがないことと平均回帰性が主な性質ですが、定常過程の持つ性質を満たさないデータをモデル化する際に有効な考え方が、「単位根過程」です。単位根過程は非定常過程の代表であり、原系列(元データ)ytが非定常過程で、差分系列Δyt=yt-yt-1が定常過程であるとき、その過程を単位根過程と呼びます。

単位根過程は1次和分過程(1次Iモデル)とも呼ばれ、その差分系列が(p,q)次ARMAモデルで表されるとき、単位根過程は(p,1,q)次ARIMAモデル(自己回帰和分移動平均)と呼ばれます。一般に、d階差分を取った系列が(p,q)次ARMA過程に従う過程を(p,d,q)次ARIMAモデルと呼びます。

一言で述べると、データの差分を取ってからARMAを適用したモデルが、ARIMAモデルということになります。

それでは、単位根過程データを人工的に作り、そのデータをARIMAモデルで解析してみましょう。

まず、解析用のデータを、下記の式に従い(1, 1, 1)次ARIMAモデルで作成します。Φ1=1なので、このデータは非定常であることが分かります。

yt = 1+ yt−1 + εt - εt−1

そしてこの式を変形すると、

Δyt=yt-yt-1 = 1 + εt - εt−1

となるので、差分データΔytはc=1, θ1=-1の1次MAモデルであることが分かります。

#データの作成

import numpy as np
import matplotlib.pyplot as plt

phi = 1
theta = -1
sigma = 10
center = 1
sample = 200

noize = np.random.normal(loc=0, scale=sigma, size=sample+1)
y0 = center
y = np.zeros(sample+1)
delta_y = np.zeros(sample+1)
y[0] = y0
delta_y[0] = 0

for k in range(sample):
    y[k+1] = center + phi * y[k] + noize[k+1] + theta*noize[k]
    delta_y[k+1] = y[k+1] - y[k]
y = y
delta_y = delta_y

fig = plt.figure(figsize=(10, 3))
ax1 = fig.add_subplot(1, 2, 1)
ax1.set_title("c=1, Φ=1, θ=-1, ε=10")
ax1.plot([int(i) for i in range(1, sample+1)],y)

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_title("delta_y data")
ax2.plot([int(i) for i in range(1, sample+1)],delta_y)

plt.show()

元データは右上がりのトレンドが見られますが、差分データは特にトレンドは見られず、0を中心に上下にばらついていることから、定常過程であることが分かります。

それでは、このデータを使ってARIMAモデルでパラメータを推定してみます。

import statsmodels.api as sm

for i in range(2): #0次と1次の和分過程で探索
    delta_y = np.diff(y, n=i)
    resdiff = sm.tsa.arma_order_select_ic(delta_y, ic="aic", trend="nc")
    print(resdiff)
0次和分過程
{'aic':              0            1            2
0          NaN  2237.244350  2073.501314
1  1615.698003  1530.025649  1530.031841
2  1552.847956  1529.986114  1533.685171
3  1553.196278  1552.800983  1532.931096
4  1544.370244  1533.666189  1535.477473, 'aic_min_order': (2, 1)}


1次和分過程
{'aic':              0            1            2
0          NaN  1515.417150  1515.424120
1  1538.361931  1515.378152  1517.307646
2  1525.871642  1517.371654  1519.240139
3  1518.737854  1519.060110  1520.164403
4  1520.168956  1519.165132  1521.732557, 'aic_min_order': (1, 1)}

AICで最適パラメータを計算すると、1次和分過程の(p,q)=(1,1)のときにAICが最小となることが分かりました。従って、与えられた時系列データは、(1, 1, 1)次ARIMAモデルで精度よく推定できると思われます。

(1, 1, 1)次ARIMAモデルのパラメータを推定します。

from statsmodels.tsa.arima_model import ARIMA

arima_1_1_1 = ARIMA(y, order=(1, 1, 1)).fit(dist=False)
arima_1_1_1.summary()
ARIMA Model Results
Dep. Variable:	D.y	No. Observations:	199
Model:	ARIMA(1, 1, 1)	Log Likelihood	-737.442
Method:	css-mle	S.D. of innovations	9.788
Date:	Fri, 03 Apr 2020	AIC	1482.885
Time:	08:13:42	BIC	1496.058
Sample:	1	HQIC	1488.216
coef	std err	z	P>|z|	[0.025	0.975]
const	0.9966	0.043	23.041	0.000	0.912	1.081
ar.L1.D.y	0.0103	0.085	0.121	0.904	-0.157	0.178
ma.L1.D.y	-0.9477	0.051	-18.404	0.000	-1.049	-0.847
Roots
Real	Imaginary	Modulus	Frequency
AR.1	96.9090	+0.0000j	96.9090	0.0000
MA.1	1.0551	+0.0000j	1.0551	0.0000

ARIMAモデルのパラメータは、c=1.0, θ1=-0.95, Φ1=0.01, ε=9.79と推定できました。実測値は、c=1, θ1=-0.1, Φ1=1, ε=10なので、Φ1がうまく予測できませんでした。データ数を増やすなどすれば、精度は上がる可能性があります。

SARIMAモデルとは

最後に、SARIMAモデルを見ます。単位根過程に対しては、ARIMAモデルでまずまず予測することができますが、季節変動がある非定常データはARIMAモデルではうまくフィットしません。そこで、季節変動があるデータに対して有効なのが、ARIMAモデルを拡張したSARIMAモデルです。具体的には、時系列方向にARIMAモデルを使い、かつ、周期方向にもARIMAモデルを使っているモデルです。

最終的に、SARIMAモデルでは7つの次数を決めます。それは、時系列方向の(p, d, q)次ARIMAのp, d, q、周期方向の(P, D, Q)次ARIMAのP, D, Q、そして周期sの7つです。周期sはグラフを見て判断できることが多いです。

それでは、Rの時系列データ解析の例でよく用いられるAirPassengersのデータを使って、SARIMAモデルで実際に解析してみます。

それでは、AirPassengersのデータを読み込んで、生データを見てみましょう。

import pandas as pd
import matplotlib.pyplot as plt

data = pd.read_csv("AirPassengers.csv", 
                   index_col="Month", 
                   parse_dates=True, 
                   dtype="float")

fig = plt.figure()
ax1 = fig.add_subplot(1, 1, 1)
ax1.plot(data["#Passengers"])

ax1.set_title("Passengers trend")
ax1.set_xlabel("Year")
ax1.set_ylabel("Passengers")

plt.show()

乗客数は右上がりになっており、増加傾向にあることが分かります。また、周期的に乗客数が変化していることも分かります。非定常過程となる要因は2つあり、1つはトレンド、もう1つは季節変動ですが、このデータは両者を含んでいることになります。

SARIMAモデルを使えば、モデル中の次数を設定することで、これまで述べてきたMAモデル、ARモデル、ARMAモデル、ARIMAモデルも表すことができますので、汎用性の高いモデルと言えます。それでは最後は機械学習らしく、SARIMAモデルで訓練用データで学習して未知のデータを予測してみたいと思います。

data_2 = np.log(data) #対数変化したデータを使う
train = np.log(data[:"1957-12-31"]["#Passengers"])

# 最適の次数を決めるため、グリッドサーチを行うための準備

import itertools
 
# 各パラメータの範囲を決める
p = d = q = range(0, 3) #時系列方向のARIMAモデルの次数(p, d, q)検討範囲
s_p = s_d = s_q = range(0, 2) #周期方向のARIMAモデルの次数(P, D, Q)検討範囲
 
# p, d, q の組み合わせを表すリストを作成
pdq = list(itertools.product(p, d, q))
 
# P, D, Q の組み合わせを表すリストを作成。周期は、グラフからs = 12としている。
pdq_season = [(x[0], x, x, 12) for x in list(itertools.product(s_p, s_d, s_q))]

from statsmodels.tsa.statespace.sarimax import SARIMAX
import warnings
 
warnings.filterwarnings("ignore") # warnings を表示させないようにする

# グリッドサーチで、最適の次数をAICで探索
min_aic = [0, 0, 10000]
for param in pdq:
    for param_season in pdq_season:
        try:
            model_sarimax = SARIMAX(train,
                          order = param,
                          seasonal_order=param_season,
                          enforce_stationarity=True,
                          enforce_invertibility=True)
 
            results = model_sarimax.fit()
             
            print("order{}, s_order{} - AIC: {}".format(param, param_season, results.aic))
 
            if results.aic < min_aic:
                min_aic = [param, param_season, results.aic]
        except:
            continue
             
print("\AICが最も良いモデル:", min_aic)
order(0, 0, 0), s_order(0, 0, 0, 12) - AIC: 672.2469661806601
order(0, 0, 0), s_order(0, 0, 1, 12) - AIC: 563.9682273977863
order(0, 0, 0), s_order(0, 1, 0, 12) - AIC: -95.22009529373264
・・・
order(2, 2, 2), s_order(1, 1, 0, 12) - AIC: -320.52715587583884
order(2, 2, 2), s_order(1, 1, 1, 12) - AIC: -324.3508043456799
\AICが最も良いモデル: [(1, 1, 1), (1, 0, 1, 12), -361.46473041738244]


総当たりでの探索で、最適の次数は(p, d, q, P, D, Q, s)=(1, 1, 1, 1, 0, 1, 12)と決まりました。

SARIMAモデルの次数が決まったので、最適条件で予測モデルを作成します。

model_final = SARIMAX(train,
              order = min_aic[0],
              seasonal_order=min_aic, 
              enforce_stationarity=True,
              enforce_invertibility=True)
 
results = model_final.fit()
 
print(results.summary())
SARIMAX Results                                      
======================================================================
Dep. Variable:     #Passengers	No. Observations:              108
Model:  SARIMAX(1, 1, 1)x(1, 0, 1, 12)	Log Likelihood       185.732
Date:            Mon, 06 Apr 2020	AIC                       -361.465
Time:                    16:15:46	BIC                       -348.101
Sample:               01-01-1949   	HQIC                      -356.047
                     - 12-01-1957                                         
Covariance Type:                              opg                                         
======================================================================
                  coef      std err       z       P>|z|      [0.025      0.975]
------------------------------------------------------------------------
ar.L1      0.5238	     0.187     2.805     0.005      0.158      0.890
ma.L1     -0.8079      0.137    -5.918     0.000     -1.075     -0.540
ar.S.L12   0.9909      0.008    120.690    0.000      0.975      1.007
ma.S.L12  -0.6276      0.123    -5.115     0.000     -0.868     -0.387
sigma2     0.0014      0.000     7.204     0.000      0.001      0.002
======================================================================
Ljung-Box (Q):              37.92   Jarque-Bera (JB):           6.69
Prob(Q):                     0.56   Prob(JB):                   0.04
Heteroskedasticity (H):      0.38   Skew:                      -0.03
Prob(H) (two-sided):         0.00   Kurtosis:                   4.22
======================================================================

構築したモデルの妥当性を評価するために、残差を確認します。残差とは、予測値と実測値の差です。残差にクセがなくランダムであれば、作成したモデルはまずまずと言えます。

results.plot_diagnostics(figsize=(10,8), lags=20)

左上:標準化残差のプロット、右上:標準化残差のヒストグラム、左下:標準化残差のQ-Qプロット、右下:標準化残差の自己相関係数です。残差には特に特徴や傾向は見られないので、作成したモデルは特に問題ないと思われます。

最後に、構築したモデルで未知のデータを予測してみます。

pred = results.get_prediction(start=pd.to_datetime("1957-12-01"),
                              end = pd.to_datetime("1960-12-01"),
                              dynamic=False)

# 期待値と信頼区間を取り出す
pred_mean = pred.predicted_mean
pred_ci = pred.conf_int(alpha = .05)

# グラフの描画
data_2["1949-01-01":].plot(label="observed")
pred_mean.plot(label="forecast", alpha=.7, color = "r")
 
# 信頼区間の描画
plt.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='g', alpha=.2)
 
plt.xlabel("Date")
plt.ylabel("Passengers(log)")
plt.legend()

作成したSARIMAモデルで、乗客数をまずまずの精度で予測できていることが分かりました。

おわりに

非定常過程の説明に有効なARIMAモデルとSARIMAモデル、および、それらの基本となる定常過程に有効なMAモデル、ARモデル、ARMAモデルを見てきました。実際に得られるデータは非定常過程であることが多いので、ARIMAモデルとSARIMAモデルをマスターすることを目指してください。

参考サイト

Pythonによる時系列分析の基礎

https://logics-of-blue.com/python-time-series-analysis/

PythonでのARIMAモデルを使った時系列データの予測の基礎[前編]

https://blog.brains-tech.co.jp/entry/arima-tutorial-1

PythonでのARIMAモデルを使った時系列データの予測の基礎[後編]

https://blog.brains-tech.co.jp/entry/arima-tutorial-2

未来を予測するビッグデータの解析手法と「SARIMAモデル」

https://deepage.net/bigdata/2016/10/22/bigdata-analytics.html#sarima%E3%83%A2%E3%83%87%E3%83%AB

定常時系列の解析に使われるARMAモデル・SARIMAモデルとは?

https://to-kei.net/time-series-analysis/sarima_model/

時系列分析の話~時系列モデル1~

https://masamunetogetoge.com/mamodel

時系列分析の話~時系列モデル2~

https://masamunetogetoge.com/armodel

A comprehensive beginner’s guide to create a Time Series Forecast (with Codes in Python and R)

https://www.analyticsvidhya.com/blog/2016/02/time-series-forecasting-codes-python/


ABOUT ME
北爪 聖也
ダメ営業マンからデータサイエンティストへキャリアチェンジ。 技術とビジネスサイドの橋渡しが出来るため、ダメ営業マンの経験も役に立ちました。 広告代理店ADKにて3年勤務→データ分析受託の会社DATUM STUDIOにて1.2年勤務後、独立。
【事例集プレゼント】業務効率化したい医薬業界の方

株式会社piponでは医薬業界の企業様向けにDXの成功事例を集めた医薬DX事例集をe-bookとしてご提供しております。

ご興味ある方がいらっしゃいましたらこちらのフォームよりご連絡頂けると嬉しいです。