2019/09/23
2020/04/14
定常時系列の解析に使われるARMAモデル・SARIMAモデルとは?
この記事では、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\)) にすると予測精度がどのように変わるのか試してみるとよいでしょう。
Recommended