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

 2019/09/23    時系列分析    


この記事では、Pythonで時系列解析・分析を行っていくうえでの基礎知識として、伝統的な時系列モデルの一つの形ともいえるSARIMAモデルについて説明します。平均・分散、正規分布などの統計学の初歩の知識を前提とした初心者向け入門記事です。手元にある分析したい時系列データに当てはめて、理解を深めてみてください。

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

まずは定常な時系列に対するモデルとしてARモデル・MAモデル・ARMAモデルを紹介します。 ARMAモデルはARモデルとMAモデルを組み合わせたものであるため、まずはARモデルとMAモデルについて理解することが重要です。

ARモデルとMAモデル

ARモデル (AutoRegressive model) は日本語で自己回帰モデルといいます。 自己共分散・自己相関と同様に、その系列の過去の値で回帰を行うため、そのような名前がついています。 ARモデルの一般形は次式で表されます。 ここで、 \(c\) は定数、 \(βi\) は回帰係数、 \(εt\) は \(t\) 時点でのノイズ、 \(σε\) は \(εt\) の分散です。
$$
y_t = c + \sum_{i=1}^p \phi_{i}y_{t-i} + \varepsilon_t,~\varepsilon_t\sim N(0,\sigma_\varepsilon^2).
$$
\(t\) 時点での値 \(yt\) が、 \(yt−1,yt−2,…,yt−p\) によって回帰されています。 このように、 \(p\) 時点過去までの値を使うARモデルをAR( \(p\) )モデルと言います。 AR( \(p\) )モデルは定常な時系列に対して有効ですが、推定された係数次第では必ずしも定常になるわけではありません。 AR( \(p\) )モデルの定常性に関しては次節で解説します。

MAモデル (Moving Average) は日本語で移動平均モデルといい、その一般形は次式で表されます。 ここで、 \(θj\) は回帰係数です。
$$
y_t = c + \sum_{j=1}^q \theta_{j}\varepsilon_{t-j} + \varepsilon_t,~\varepsilon_t\sim N(0,\sigma_\varepsilon^2).
$$
AR( \(p\) )の式と似ていますが、過去の値 \(yt−i\) で回帰しているのではなく、過去のノイズ \(εt−j\) で回帰しています。 このように、 \(q\) 時点過去までのノイズを使うMAモデルをMA( \(q\) )モデルといいます。 過去のデータで回帰するARモデルに比べ、観測不可能なノイズで回帰するMAモデルは直感的に理解しにくいですが、たとえば「今日の売上が予想より多ければ明日の売上は少なくなる」というような現象をモデリングすることができます。 ARモデルと異なり、MAモデルは常に定常になります。

AR( \(p\) )モデルとMA( \(q\) )モデルを組み合わせたARMA( \(p,q\) )モデルの一般形は次式で表されます。
$$
y_t = c + \sum_{i=1}^p \phi_{i}y_{t-i} + \sum_{j=1}^q \theta_{j}\varepsilon_{t-j} + \varepsilon_t,
~\varepsilon_t\sim N(0,\sigma_\varepsilon^2).
$$
ARにMAを加えたことにより柔軟性が高まり、定常な時系列に対しては非常に強力な説明力・予測力を持ちます。

ARの定常性とMAの反転可能性

たとえば次のようなARモデルを考えてみましょう。 誤差が加算されているものの、前期の値が2倍になっていくため、あっという間に発散してしまいます。
$$
y_t = c + 2y_{t-1} + \varepsilon_t
$$
このように、ARモデルは必ずしも定常になるわけではないのです。 一般に、次の \(z\) に関する方程式の解の絶対値が1よりも大きいとき、AR( \(p\) )モデルは定常になることが知られています。
$$
1 - \phi_1z^1-\dots-\phi_pz^p=0
$$
MAモデルは必ず定常であるため、定常なARモデルはMA( \(∞\) )モデルで表すことができます。

一方で、MAモデルをAR( \(∞\) )モデルで表せるとき、そのMAモデルは反転可能性 (Invertibility) を持つといいます。 MAが反転可能性を持つとARモデルで表せる、つまり誤差 \(εt\) を \(yt−k (k=1,2,…) \)で表せるため、誤差がどの程度の大きさになるかを過去のデータから判断できることになります。

ARIMAモデル (AR Integrated MA model)

ARMAモデルは定常な時系列のモデルであり、非定常な時系列は定常な時系列に変換した後にARMAモデルに当てはめる必要がありました。 \(d\) 階差分をとった系列に対してARMA( \(p,q\) )を考えるモデルをARIMA( \(p,d,q\) )モデルといいます。 時系列が ( \(d\) 次式のトレンド) + (定常部分) といえる場合に有効なモデルです。

SARIMAモデル (Seasonal ARIMA model)

AirPassengersデータのように、明らかな周期性を持つデータが世の中には多く存在します。 ARIMAモデルに周期成分を取り入れたモデルをSARIMAモデルといいます。 SARIMAモデルのアイデアは、時系列方向の説明にARIMA( \(p,d,q\) ) モデルを使うだけでなく、周期方向の説明にもARIMA( \(P,D,Q\) )モデルを使おう、というものです。

SARIMAモデルでは合計7個の次数があります。 時系列方向のARIMA( \(p,d,q\) )に加え季節差分方向のARIMA( \(P,D,Q\) )、さらには周期 \(s\) があるためです。 これをSARIMA( \(p,d,q\) )( \(P,D,Q\) )[ \(s\) ]と表記することがあります。 それぞれの次数に対し0か1を考えるとしても、 27=128 通りのモデルを考えなければなりません。 このような組み合わせ爆発の問題を回避するため、周期 \(s\) は作図や自己相関関数をもとに決め打ちし、季節差分の \(P,D,Q \)は低く抑えることがよく行われます。

モデル選択とAIC

では、どのようにして (\(p,d,q,P,D,Q,s\)) の組み合わせを決めればよいのでしょうか。 何か数値的な「モデルの良さ」の指標があれば、様々な組み合わせの次数でモデルを作り、その指標に従って決定することができます。 本節ではその一つの指標として、赤池情報量基準 (Akaike Information Criterion: AIC)を紹介します。

尤度とは、モデルの尤(もっと)もらしさのことです。 あるパラメータを持つモデルにデータが従うと仮定したとき、そのデータが得られる確率はそのパラメータの関数となり、これを尤度 (Likelihood) といいます。 尤度の対数を最大にするようにパラメータを推定する方法は最尤推定とよばれ、その時の尤度を最大対数尤度といいます。 用いるパラメータを増やせば増やすほど、つまりモデルが複雑であればあるほど、基本的に当てはまりは良くなります。 一方で、パラメータが多いほどデータに過剰に適合してしまい、未知のデータに対する精度が悪くなる過学習が起こってしまいます。 AICは次式で定義され、モデルを複雑にすること(パラメータを増やすこと)にペナルティがかかっています。
$${\rm AIC}=-2(最大対数尤度)+2(パラメータ数)$$
符号より、AICは小さければ小さいほど良い指標です。 よって、様々な (\(p,d,q,P,D,Q,s\)) でSARIMAモデルを最尤推定し、AICが最も小さくなるようなモデルを選択する方針が考えられます。ただし、AICは最大対数尤度に基づいているため、データセットが同じならば異なるモデルを比較することができますが、データが異なる場合は比較できないことに注意してください。

 

実践!PythonによるSARIMAモデルの推定

実際にPythonを使ったSARIMAモデルの推定を行います。 1960年末までのAirPassengersデータのうち、1957年末までのデータを訓練データとし、それ以降のデータをテストデータとすることにします。 訓練データをもとにモデル選択とパラメータ推定を行い、テストデータで予測精度を検証します。 まずは、対数変換した訓練データをtrainという名前で作ります。 データの大きさをshapeで確認すると、144個のデータのうち108個を訓練データとして確保したことがわかりました。

train = np.log(df[:"1957-12-31"].Passengers)
print(train.shape)
print(Passengers_log.shape)

出力結果です。

(108,)
(144,)

SARIMAモデルを推定するには、statsmodelsのsarimax()を用います。 statsmodelsは巨大なモジュールでsarimax()まで辿るのが大変なため、ここではあらかじめSARIMAXという名前で読み込んでおきます。 これまでの流れを理解したうえで、sarimax()を適切に使えるようになることが本節の目標です。

from statsmodels.tsa.statespace.sarimax import SARIMAX

SARIMA( \(p,d,q\) )( \(P,D,Q\) )[ \(s\) ]の次数のうち、 \(s\) に関しては明らかに12の周期を持つとわかっているため、 \(s=12\) で決め打ちをします。 それ以外の次数については、それぞれに最小値と最大値を定めてその全通りを計算し、適切と思われるものを選択します。 この作業をグリッドサーチといい、機械学習におけるパラメータチューニングなどで広く行われます。 それぞれの次数の組み合わせを列挙することは面倒に感じられますが、range()とitertoolsのproduct()で簡単に行うことができます。 また、sarimax()では \(p,d,q\) をorderとして、 \(P,D,Q,s\) をseasonal_orderとして指定する必要があるため、別々に作成をしています。 ここでは、 \(p,d,q\) を0から2、 \(P,D,Q\) を0か1としました。 より広範囲を探索したい場合は、計算資源の余裕に応じて大きくすると良いでしょう。

試しに \(p,d,q\) のはじめの6通りを表示すると、確かに全通り列挙できていそうなことがわかります。

import itertools

# 各パラメータの範囲を決める
p = d = q = range(0, 3)
sp = sd = sq = range(0, 2)

# p, d, q の組み合わせを列挙するリストを作成
pdq = list(itertools.product(p, d, q))

# P, D, Q の組み合わせを列挙するリストを作成すると同時に、後ろに s = 12 を決め打ちでつけている。
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(sp, sd, sq))]
pdq[0:6]

出力結果です。

[(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 1, 0), (0, 1, 1), (0, 1, 2)]

ここからがsarimax()を使ってモデル推定をする核の部分になります。 (\(p,d,q\)) のリストと (\(P,D,Q,s\)) のリストのfor文を入れ子にすることを考え、より小さいAICの値になればbest_resultを更新する、ということを行っています。 また、sarimax()の引数enforce_stationarityとenforce_invertibilityに Trueを指定しています。 これは、「推定されるAR部分が定常性を持つように強制するか否か」と「推定されるMA部分が反転可能性を持つように強制するか否か」の指定になります。 前で述べた通り、これらが満たされないと十分な解析ができないので、基本的にTrueを指定することになります。 条件にそぐわない場合はエラーになりますが、sarimax()をtry文の中で使用することで、これを解決しています。 そのため、実行結果のAICには時折nanが入っています。

import warnings

warnings.filterwarnings("ignore") # warnings を表示させないようにする


best_result = [0, 0, 10000000]
for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            mod = SARIMAX(train,
                          order = param,
                          seasonal_order=param_seasonal,
                          enforce_stationarity=True,
                          enforce_invertibility=True)

            results = mod.fit()
            
            print('order{}, s_order{} - AIC: {}'.format(param, param_seasonal, results.aic))

            if results.aic < best_result[2]:
                best_result = [param, param_seasonal, results.aic]
        except:
            continue
            
print('\AICが最も良いモデル:', best_result)
order(0, 0, 0), s_order(0, 1, 1, 12) - AIC: -128.25128590729037
order(0, 0, 0), s_order(1, 1, 0, 12) - AIC: -164.81181183270064
order(0, 0, 0), s_order(1, 1, 1, 12) - AIC: -206.70631175832395
order(0, 0, 1), s_order(0, 0, 0, 12) - AIC: 530.3127699561581
order(0, 0, 1), s_order(0, 1, 0, 12) - AIC: -173.16257895032624
order(0, 0, 1), s_order(0, 1, 1, 12) - AIC: -186.86915329000297
order(0, 0, 1), s_order(1, 1, 0, 12) - AIC: -203.41594267184533
order(0, 0, 1), s_order(1, 1, 1, 12) - AIC: -242.44764861973928
order(0, 0, 2), s_order(0, 1, 0, 12) - AIC: nan
order(0, 0, 2), s_order(0, 1, 1, 12) - AIC: nan
order(0, 0, 2), s_order(1, 1, 0, 12) - AIC: nan
order(0, 0, 2), s_order(1, 1, 1, 12) - AIC: nan
order(0, 1, 0), s_order(0, 0, 1, 12) - AIC: -247.33926853799974
order(0, 1, 0), s_order(0, 1, 1, 12) - AIC: -335.0218187613289
order(0, 1, 0), s_order(1, 0, 0, 12) - AIC: -331.51534614835134
order(0, 1, 0), s_order(1, 0, 1, 12) - AIC: -354.48547395113894
order(0, 1, 0), s_order(1, 1, 0, 12) - AIC: -325.20551291030637
order(0, 1, 0), s_order(1, 1, 1, 12) - AIC: -333.0388295470783
order(0, 1, 1), s_order(0, 0, 0, 12) - AIC: -180.9217751062901
order(0, 1, 1), s_order(0, 0, 1, 12) - AIC: -246.12191914018157
order(0, 1, 1), s_order(0, 1, 0, 12) - AIC: -317.94780495592227
order(0, 1, 1), s_order(0, 1, 1, 12) - AIC: -344.5860243623633
order(0, 1, 1), s_order(1, 0, 0, 12) - AIC: -337.7402412377094
order(0, 1, 1), s_order(1, 0, 1, 12) - AIC: -361.1372283626219
order(0, 1, 1), s_order(1, 1, 0, 12) - AIC: -338.3558769302148
order(0, 1, 1), s_order(1, 1, 1, 12) - AIC: -342.7418283184463
order(0, 1, 2), s_order(0, 0, 0, 12) - AIC: -183.79809672445728
order(0, 1, 2), s_order(0, 0, 1, 12) - AIC: -244.83678568964754
order(0, 1, 2), s_order(0, 1, 0, 12) - AIC: -316.1078959233711
order(0, 1, 2), s_order(0, 1, 1, 12) - AIC: -342.8500892741997
order(0, 1, 2), s_order(1, 0, 0, 12) - AIC: -336.1241091596002
order(0, 1, 2), s_order(1, 0, 1, 12) - AIC: -359.63149238716346
order(0, 1, 2), s_order(1, 1, 0, 12) - AIC: -336.45222818008716
order(0, 1, 2), s_order(1, 1, 1, 12) - AIC: -341.0098789510462
order(0, 2, 0), s_order(0, 0, 1, 12) - AIC: -179.97216783915235
order(0, 2, 0), s_order(0, 1, 1, 12) - AIC: -239.7672632709707
order(0, 2, 0), s_order(1, 0, 0, 12) - AIC: -232.70638761901105
order(0, 2, 0), s_order(1, 0, 1, 12) - AIC: -254.25602507110472
order(0, 2, 0), s_order(1, 1, 0, 12) - AIC: -227.0329178233908
order(0, 2, 0), s_order(1, 1, 1, 12) - AIC: -237.9810151822921
order(0, 2, 1), s_order(0, 0, 0, 12) - AIC: -169.48899943697091
order(0, 2, 1), s_order(0, 0, 1, 12) - AIC: -236.16479154714102
order(0, 2, 1), s_order(0, 1, 0, 12) - AIC: -296.5772354982875
order(0, 2, 1), s_order(0, 1, 1, 12) - AIC: -322.2859269949689
order(0, 2, 1), s_order(1, 0, 0, 12) - AIC: -324.0342210823416
order(0, 2, 1), s_order(1, 0, 1, 12) - AIC: -347.28496813217384
order(0, 2, 1), s_order(1, 1, 0, 12) - AIC: -313.4625390157218
order(0, 2, 1), s_order(1, 1, 1, 12) - AIC: -320.34758072914957
order(0, 2, 2), s_order(0, 0, 0, 12) - AIC: -172.71141153656507
order(0, 2, 2), s_order(0, 0, 1, 12) - AIC: -237.72801620103127
order(0, 2, 2), s_order(0, 1, 0, 12) - AIC: nan
order(0, 2, 2), s_order(0, 1, 1, 12) - AIC: nan
order(0, 2, 2), s_order(1, 0, 0, 12) - AIC: -322.0261053978313
order(0, 2, 2), s_order(1, 0, 1, 12) - AIC: -345.25281641670404
order(0, 2, 2), s_order(1, 1, 0, 12) - AIC: nan
order(0, 2, 2), s_order(1, 1, 1, 12) - AIC: nan
order(1, 0, 0), s_order(0, 1, 0, 12) - AIC: -309.7053673802187
order(1, 0, 0), s_order(0, 1, 1, 12) - AIC: -334.6714173300306
order(1, 0, 0), s_order(1, 1, 0, 12) - AIC: -325.9199293221464
order(1, 0, 0), s_order(1, 1, 1, 12) - AIC: -332.6735219011772
order(1, 0, 1), s_order(0, 1, 0, 12) - AIC: -318.576765729886
order(1, 0, 1), s_order(0, 1, 1, 12) - AIC: -344.11585153256067
order(1, 0, 1), s_order(1, 1, 0, 12) - AIC: -338.2464599690521
order(1, 0, 1), s_order(1, 1, 1, 12) - AIC: -342.27855521245306
order(1, 0, 2), s_order(0, 1, 0, 12) - AIC: -316.68249937651336
order(1, 0, 2), s_order(0, 1, 1, 12) - AIC: -342.34397485549505
order(1, 0, 2), s_order(1, 1, 0, 12) - AIC: -336.3178269211952
order(1, 0, 2), s_order(1, 1, 1, 12) - AIC: -340.51839177450944
order(1, 1, 0), s_order(0, 0, 0, 12) - AIC: -179.2833871186191
order(1, 1, 0), s_order(0, 0, 1, 12) - AIC: -246.00491335905156
order(1, 1, 0), s_order(0, 1, 0, 12) - AIC: -316.8473621071101
order(1, 1, 0), s_order(0, 1, 1, 12) - AIC: -343.51032403473073
order(1, 1, 0), s_order(1, 0, 0, 12) - AIC: -337.1268196002543
order(1, 1, 0), s_order(1, 0, 1, 12) - AIC: -360.3781733242543
order(1, 1, 0), s_order(1, 1, 0, 12) - AIC: -336.57960574446133
order(1, 1, 0), s_order(1, 1, 1, 12) - AIC: -341.5629809069127
order(1, 1, 1), s_order(0, 0, 0, 12) - AIC: -183.59048483311363
order(1, 1, 1), s_order(0, 0, 1, 12) - AIC: -244.25575117254903
order(1, 1, 1), s_order(0, 1, 0, 12) - AIC: -316.2802618661529
order(1, 1, 1), s_order(0, 1, 1, 12) - AIC: -343.2987165793862
order(1, 1, 1), s_order(1, 0, 0, 12) - AIC: -335.28907784311053
order(1, 1, 1), s_order(1, 0, 1, 12) - AIC: -361.4645649737698
order(1, 1, 1), s_order(1, 1, 0, 12) - AIC: -336.518165442526
order(1, 1, 1), s_order(1, 1, 1, 12) - AIC: -341.41164350119175
order(1, 1, 2), s_order(0, 0, 0, 12) - AIC: -187.76815515299234
order(1, 1, 2), s_order(0, 0, 1, 12) - AIC: -250.1270053221665
order(1, 1, 2), s_order(0, 1, 0, 12) - AIC: -314.58577818182096
order(1, 1, 2), s_order(0, 1, 1, 12) - AIC: -340.58694683910556
order(1, 1, 2), s_order(1, 0, 0, 12) - AIC: -336.3229231783102
order(1, 1, 2), s_order(1, 0, 1, 12) - AIC: -360.1071325666829
order(1, 1, 2), s_order(1, 1, 0, 12) - AIC: -334.3817313525976
order(1, 1, 2), s_order(1, 1, 1, 12) - AIC: -338.7419400466974
order(1, 2, 0), s_order(0, 0, 0, 12) - AIC: -131.11853542980705
order(1, 2, 0), s_order(0, 0, 1, 12) - AIC: -197.09816034225756
order(1, 2, 0), s_order(0, 1, 0, 12) - AIC: -262.88892885849987
order(1, 2, 0), s_order(0, 1, 1, 12) - AIC: -288.76118602007386
order(1, 2, 0), s_order(1, 0, 0, 12) - AIC: -280.2963773241147
order(1, 2, 0), s_order(1, 0, 1, 12) - AIC: -301.87965172506296
order(1, 2, 0), s_order(1, 1, 0, 12) - AIC: -279.3205672751532
order(1, 2, 0), s_order(1, 1, 1, 12) - AIC: -287.023102534006
order(1, 2, 1), s_order(0, 0, 0, 12) - AIC: -171.04789919232468
order(1, 2, 1), s_order(0, 0, 1, 12) - AIC: -234.61367589904273
order(1, 2, 1), s_order(0, 1, 0, 12) - AIC: -305.33362313395315
order(1, 2, 1), s_order(0, 1, 1, 12) - AIC: -330.27227727470904
order(1, 2, 1), s_order(1, 0, 0, 12) - AIC: -329.44282701743106
order(1, 2, 1), s_order(1, 0, 1, 12) - AIC: -351.8621194543448
order(1, 2, 1), s_order(1, 1, 0, 12) - AIC: -324.07411002244356
order(1, 2, 1), s_order(1, 1, 1, 12) - AIC: -326.61540740897533
order(1, 2, 2), s_order(0, 0, 0, 12) - AIC: -175.33622852743667
order(1, 2, 2), s_order(0, 0, 1, 12) - AIC: -235.84861222217373
order(1, 2, 2), s_order(0, 1, 0, 12) - AIC: -303.3774845041204
order(1, 2, 2), s_order(0, 1, 1, 12) - AIC: -328.2363147792629
order(1, 2, 2), s_order(1, 0, 0, 12) - AIC: -327.73778538645985
order(1, 2, 2), s_order(1, 0, 1, 12) - AIC: -350.5331531053591
order(1, 2, 2), s_order(1, 1, 0, 12) - AIC: -321.98542591877975
order(1, 2, 2), s_order(1, 1, 1, 12) - AIC: -326.1660712404815
order(2, 0, 0), s_order(0, 1, 0, 12) - AIC: -317.94905838727584
order(2, 0, 0), s_order(0, 1, 1, 12) - AIC: -343.1007134437114
order(2, 0, 0), s_order(1, 1, 0, 12) - AIC: -336.7026026187913
order(2, 0, 0), s_order(1, 1, 1, 12) - AIC: -341.16381705131005
order(2, 0, 1), s_order(0, 1, 0, 12) - AIC: -316.8022705297305
order(2, 0, 1), s_order(0, 1, 1, 12) - AIC: -342.35270289334613
order(2, 0, 1), s_order(1, 1, 0, 12) - AIC: -336.36361548593186
order(2, 0, 1), s_order(1, 1, 1, 12) - AIC: -338.05167480345335
order(2, 0, 2), s_order(0, 1, 0, 12) - AIC: -314.6389704366976
order(2, 0, 2), s_order(0, 1, 1, 12) - AIC: -340.2469076138487
order(2, 0, 2), s_order(1, 1, 0, 12) - AIC: -334.5259618289924
order(2, 0, 2), s_order(1, 1, 1, 12) - AIC: -338.0155479390361
order(2, 1, 0), s_order(0, 0, 0, 12) - AIC: -180.98736936160856
order(2, 1, 0), s_order(0, 0, 1, 12) - AIC: -244.90919916330262
order(2, 1, 0), s_order(0, 1, 0, 12) - AIC: -314.8810489837429
order(2, 1, 0), s_order(0, 1, 1, 12) - AIC: -341.6019165053905
order(2, 1, 0), s_order(1, 0, 0, 12) - AIC: -335.12825318394755
order(2, 1, 0), s_order(1, 0, 1, 12) - AIC: -358.4231800371768
order(2, 1, 0), s_order(1, 1, 0, 12) - AIC: -334.86341580635593
order(2, 1, 0), s_order(1, 1, 1, 12) - AIC: -339.6776839950578
order(2, 1, 1), s_order(0, 0, 0, 12) - AIC: -188.1778735634631
order(2, 1, 1), s_order(0, 0, 1, 12) - AIC: -251.4907024592643
order(2, 1, 1), s_order(1, 0, 0, 12) - AIC: -337.7805156016141
order(2, 1, 1), s_order(1, 0, 1, 12) - AIC: -361.1752994327316
order(2, 1, 2), s_order(0, 0, 0, 12) - AIC: -186.42019891125844
order(2, 1, 2), s_order(0, 0, 1, 12) - AIC: -249.05298147090974
order(2, 1, 2), s_order(0, 1, 0, 12) - AIC: -314.74494025651165
order(2, 1, 2), s_order(0, 1, 1, 12) - AIC: -341.16829449523664
order(2, 1, 2), s_order(1, 0, 0, 12) - AIC: -337.4554235638826
order(2, 1, 2), s_order(1, 0, 1, 12) - AIC: -359.85795545879984
order(2, 1, 2), s_order(1, 1, 0, 12) - AIC: -333.95899368015245
order(2, 1, 2), s_order(1, 1, 1, 12) - AIC: -339.1811873368714
order(2, 2, 0), s_order(0, 0, 0, 12) - AIC: -141.54299624128384
order(2, 2, 0), s_order(0, 0, 1, 12) - AIC: -202.25054739245462
order(2, 2, 0), s_order(0, 1, 0, 12) - AIC: -267.7142653242864
order(2, 2, 0), s_order(0, 1, 1, 12) - AIC: -296.1904852998152
order(2, 2, 0), s_order(1, 0, 0, 12) - AIC: -285.21626711146513
order(2, 2, 0), s_order(1, 0, 1, 12) - AIC: -309.5595490358303
order(2, 2, 0), s_order(1, 1, 0, 12) - AIC: -287.2771580844477
order(2, 2, 0), s_order(1, 1, 1, 12) - AIC: -294.2084548600758
order(2, 2, 1), s_order(0, 0, 0, 12) - AIC: -172.68028009942793
order(2, 2, 1), s_order(0, 0, 1, 12) - AIC: -233.16167230406055
order(2, 2, 1), s_order(1, 0, 0, 12) - AIC: -327.4411237191209
order(2, 2, 1), s_order(1, 0, 1, 12) - AIC: -350.3683291437861
order(2, 2, 2), s_order(0, 0, 0, 12) - AIC: -173.33940670439299
order(2, 2, 2), s_order(0, 0, 1, 12) - AIC: -234.0270522978273
order(2, 2, 2), s_order(1, 0, 0, 12) - AIC: -325.94123287309867
order(2, 2, 2), s_order(1, 0, 1, 12) - AIC: -348.41075513083143
\AICが最も良いモデル: [(1, 1, 1), (1, 0, 1, 12), -361.4645649737698]

グリッドサーチにより、SARIMA(2, 1, 1)(1, 0, 1)[12]が良さそうだとわかりました。 では、実際に予測に使うモデルの学習をします。 モデルをsarimax()で決め、fit()で推定した結果に対しsummary()メソッドを適用すると、推定されたモデルの概要を知ることができます。 出力結果には様々な情報が含まれていますが、主なものにar.L1(AR部分のラグ1)などの各係数やAICなどの情報量基準、各種検定結果などがあります。

mod = SARIMAX(train,
              order=(2, 1, 1),
              #order = best_result[0],
              seasonal_order=(1, 0, 1, 12), 
              #seasonal_order=best_result[1],
              enforce_stationarity=True,
              enforce_invertibility=True)

results = mod.fit()

print(results.summary())

出力結果です。

                                 Statespace Model Results                                 
==========================================================================================
Dep. Variable:                         Passengers   No. Observations:                  108
Model:             SARIMAX(2, 1, 1)x(1, 0, 1, 12)   Log Likelihood                 186.588
Date:                            Sun, 14 Apr 2019   AIC                           -361.175
Time:                                    11:02:11   BIC                           -345.083
Sample:                                01-01-1949   HQIC                          -354.650
                                     - 12-01-1957                                         
Covariance Type:                              opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.5969      0.132      4.507      0.000       0.337       0.856
ar.L2          0.1624      0.105      1.542      0.123      -0.044       0.369
ma.L1         -0.9351      0.089    -10.565      0.000      -1.109      -0.762
ar.S.L12       0.9908      0.009    109.126      0.000       0.973       1.009
ma.S.L12      -0.6005      0.135     -4.440      0.000      -0.866      -0.335
sigma2         0.0013      0.000      6.652      0.000       0.001       0.002
===================================================================================
Ljung-Box (Q):                       34.95   Jarque-Bera (JB):                 4.62
Prob(Q):                              0.70   Prob(JB):                         0.10
Heteroskedasticity (H):               0.36   Skew:                             0.04
Prob(H) (two-sided):                  0.00   Kurtosis:                         4.02
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

推定されたモデルが妥当なものかどうか、説明しきれない残差(訓練期間におけるもの)の診断をしなければなりません。 もしモデルが妥当であれば、説明しきれない部分はホワイトノイズになっているはずだからです。 plot_diagnostics()メソッドを使用すると、残差に関する4つのプロットを見ることができます。 左上から順に、標準化した残差の時系列プロット、残差のヒストグラムと正規分布(とKDE分布)の密度関数、残差の正規QQプロット、残差の自己相関関数です。 残差の自己相関は低くまとまっており、ほぼ問題ないと言えるでしょう。

results.plot_diagnostics(lags=20);

いよいよ予測に入ります。 予測はget_prediction()メソッドで行うことができます。 主要な引数には予測する期間の始めstartと終わりendに加え、動的な予測を行うかどうかのdynamicがあります。 dynamic = Trueとすると、 \(yt\) の値を次時点 \(yt+1\) の予測に使用することができますが、今回は訓練データとテストデータを切り離しているため、Falseを指定しています。 また、予測の期待値はpredicted_meanで、区間予測はconf_int()で得ることができます。

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)

予測結果を可視化してみましょう。 以下のコードでは、まず正解データであるPassengers_logを時系列プロットしたのち、予測の期待値を赤線で、95%区間をブラケットで描画しています。

# 折れ線の描画
Passengers_log['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='r', alpha=.2)

plt.xlabel('Date')
plt.ylabel('Passengers')
plt.legend();

予測の期待値は常に実測値の上にありやや過大に見積もっていますが、それなりの予測が行えていることがわかります。

SARIMAモデルをあてはめたのは対数系列に対してだったので、指数をとって原系列に戻してみましょう。 numpyのexp()を使います。 図より、原系列に戻してもほぼ問題なく予測ができていると言えそうです。

pred_ci2 = np.exp(pred.conf_int(alpha = .05))


np.exp(Passengers_log)['1949-01-01':].plot(label='observed')
np.exp(pred_mean).plot(label='forecast', alpha=.7, color = "r")

plt.fill_between(pred_ci2.index,
                pred_ci2.iloc[:, 0],
                pred_ci2.iloc[:, 1], color='r', alpha=.2)

plt.xlabel('Date')
plt.ylabel('Passengers')
plt.legend();

最後に予測精度を出してみましょう。 ここではRoot Mean Squared Error (RMSE) を使った精度を算出してみます。 予測値を \(y^t\) とすると予測誤差は \(yt−y^t\) なので、 \(t=1,…,T\) におけるRMSEは
$$
{\rm RMSE} = \sqrt{
\frac{1}{T} \sum_{t=1}^T (y_t - \hat y_t)^2
}
$$
で定義されます。

y_forecasted = np.exp(pred.predicted_mean)
y_truth = np.exp(Passengers_log['1958-01-01':])

# RMSEの計算
rmse = np.sqrt(((y_forecasted - y_truth) ** 2).mean())
print('RMSE = {}'.format(round(rmse, 4)))

出力結果です。

RMSE = 42.1929

予測精度は必ずテストデータを用いて算出してください。 訓練データへの当てはまりはいくらでも良くすることができるためです。

まとめ

本記事ではAICによるモデル選択を行いましたが、他の (\(p,d,q,P,D,Q,s\)) にすると予測精度がどのように変わるのか試してみるとよいでしょう。

  • 人気の投稿とページ

  • コメントを残す

    メールアドレスが公開されることはありません。 * が付いている欄は必須項目です