はじめに

統計的時系列モデルを使うと下図のような予測が可能です。これが、どういう仕組みで行われるかを明らかにします。

具体的には、

PythonのstatsmodelsというライブラリのARIMAを使用して、上図のpredictを求める機能を実装
Pythonでライブラリを使わずに、同様の機能を実装

を行い、1.と2.の結果が一致するかを確かめます。

ただ予測をするだけなら1.で十分ですが、2.と合わせることで、どういう仕組みで予測しているのかを理解することが狙いです。

ただし、2.では学習は対象外とします（スクラッチで書くのは難しすぎるので）。

定常的な時系列データの予測
トレンド成分がある時系列データの予測
トレンド成分と周期成分がある時系列データの予測

の順に解説していきます。

定常的な時系列データの予測

下図の予測を実装してみます。

データ作成

まずは青色のtrueのデータを作ります。

import numpy as np
np.random.seed(0)
n = 200
x = [0]
while len(x) < n:
    x.append(-0.8 * x[-1] + np.random.normal(0, 1))
x = np.array(x)

ライブラリを使って学習＆予測

次に、k=150サンプルまでを使って、パラメータが「order=(2, 0, 0)」のモデルを学習させます。

import statsmodels.api
import statsmodels.tsa.arima as arima
k = 150
model = arima.model.ARIMA(x[:k], order=(2, 0, 0)).fit()

続いて、ライブラリを使って151～200サンプルまでを予測します。オレンジ色のpredictのデータを求めています。その後、151サンプル目の予測値をprintします。

x_pred_mean = model.get_prediction(start=k, end=n-1).predicted_mean
print(x_pred_mean[0])
# 2.249875625186254

ライブラリを使わずに予測

まず、modelの中身を見てみましょう。

print(model.summary())
#         coef  std err    z    P>|z|  [0.025  0.975]
# const	 0.0718  0.048   1.505  0.132  -0.022  0.165
# ar.L1	-0.7731  0.080  -9.644  0.000  -0.930 -0.616
# ar.L2	 0.0234  0.083   0.283  0.777  -0.139  0.186
# sigma2 1.0381  0.144   7.217  0.000   0.756  1.320

上記のmodelのパラメータを既知として、ライブラリを使わずに151サンプル目の予測をしてみましょう。

試行錯誤した結果、次の通り実行すればライブラリと同じ結果が得られることが分かりました。

const = model.params[0]
arL1 = model.params[1]
arL2 = model.params[2]
sigma2 = model.params[3]
i = k
print(arL1 * (x[i - 1] -  const) + arL2 * (x[i - 2] -  const) + const)
# 2.249875625186254

つまり、このモデルでは、$i-1$、$i-2$番目のサンプルを使って$i$番目のサンプルを予測しています。つまり、過去の「2」サンプルを使って、次のサンプルを予測しているのです。

次に、予測式を解説します。上のコードによると、予測式は次の通りです。

x_i = a_1 (x_{i-1} - c) + a_2 (x_{i - 2} - c) + c

この式がどうやって得られるかを示します。

まず、元の時系列から「$-c$」します。

元の時系列	元の時系列-c
$x_i$	①$x _ i - c$
$x_{i-1}$	②$x_{i-1} - c$
$x_{i-2}$	③$x_{i-2} - c$
$\cdots$	$\cdots$

「元の時系列-c」に対して、定数項がない2次のARモデルを適用します。実際に、上表の②③を使ってARモデルで①を予測する式を構成してみます。

① = a_1\times②+a_2\times③ \\
\Leftrightarrow (x_i - c) = a_1(x_{i-1} - c) + a_2(x_{i-2} - c)\\
\Leftrightarrow x_i = a_1 (x_{i-1} - c) + a_2 (x_{i - 2} - c) + c

これは、コードの予測式と一致します。

学習時のパラメータである「order=(2, 0, 0))」の「2」はARモデルの次数です。

トレンド成分がある時系列データの予測

データ作成

まずは上昇トレンドのある青色のtrueのデータを作ります。

np.random.seed(0)
n = 40
x = (np.arange(n) / 5) ** 2 + np.random.normal(0, 1, n)

ライブラリを使って学習＆予測

次に、k=30サンプルまでを使って、パラメータが「order=(2, 2, 0)」のモデルを学習させます。

k = 30
model = arima.model.ARIMA(x[:k], order=(2, 2, 0)).fit()

続いて、31～40サンプルまでを予測します。オレンジ色のpredictのデータを求めています。その後、31サンプル目の予測値をprintします。

x_pred_mean = model.get_prediction(start=k, end=n-1).predicted_mean
print(x_pred_mean[0])
# 38.028852151892224

ライブラリを使わずに予測

まず、modelの中身を見てみましょう。

print(model.summary())
#  coef     std    err      z    P>|z|  [0.025   0.975]
# ar.L1  -1.0913  0.145  -7.551  0.000   -1.375 -0.808
# ar.L2  -0.5686  0.134  -4.246  0.000   -0.831 -0.306
# sigma2  2.5972  0.692	  3.754  0.000    1.241  3.953

上記のmodelのパラメータを既知として、ライブラリを使わずに31サンプル目の予測をしてみましょう。

試行錯誤の結果、次の通り実行すればライブラリと同じ結果が得られることが分かりました。

arL1 = model.params[0]
arL2 = model.params[1]
sigma2 = model.params[2]
i = k
print(
    x[i - 1] + (x[i - 1] - x[i - 2]) +
    arL1 * ((x[i - 1] - x[i - 2]) - (x[i - 2] -  x[i - 3])) +
    arL2 * ((x[i - 2] - x[i - 3]) - (x[i - 3] -  x[i - 4])),
)
# 38.028852151892345

つまり、$i-1$、$i-2$、$i-3$、$i-4$サンプル目のデータを使って$i$サンプル目を予測しているのです。

予測式を解説します。上のコードによると、予測式は次の通りです。

x_i = x_{i-1} + (x_{i-1} - x_{i-2}) + a_1\{(x_{i-1} - x_{i-2})-(x_{i-2} - x_{i-3})\} + a_2\{(x_{i-2} - x_{i-3})-(x_{i-3} - x_{i-4})\}

この式がどうやって得られるかを示します。

まず、元の時系列の2回階差をとります（例えば、④ = ① - ②、⑦ = ④ - ③です）。

元の時系列	1回階差	2回階差
①$x_i$	④$x _ i - x_{i-1}$	⑦$( x_i - x_{i-1})-(x_{i-1} - x_{i-2}) $
②$x_{i-1}$	⑤$x_{i-1} - x_{i-2}$	⑧ $(x_{i-1} - x_{i-2})-(x_{i-2} - x_{i-3})$
③$x_{i-2}$	⑥$x_{i-2} - x_{i-3}$	⑨ $(x_{i-2} - x_{i-3})-(x_{i-3} - x_{i-4})$
$\cdots$	$\cdots$	$\cdots$

「2回階差」に対して、定数項がない2次のARモデルを適用します。実際に、上表の⑧⑨を使ってARモデルで⑦を予測する式を構成してみます。

⑦ = a_1\times⑧+a_2\times⑨ \\
\Leftrightarrow (x_i - x_{i-1})-(x_{i-1} - x_{i-2}) = a_1\{(x_{i-1} - x_{i-2})-(x_{i-2} - x_{i-3})\} + a_2\{(x_{i-2} - x_{i-3})-(x_{i-3} - x_{i-4})\} \\
\Leftrightarrow x_i = x_{i-1} + (x_{i-1} - x_{i-2}) + a_1\{(x_{i-1} - x_{i-2})-(x_{i-2} - x_{i-3})\} + a_2\{(x_{i-2} - x_{i-3})-(x_{i-3} - x_{i-4})\}

これは、コードの予測式と一致します。ライブラリに与えたパラメータorder=(α, β, 0)のうち、αはARモデルの次数、βは何回階差をとるかを示します。

尚、階差をとったあとにARモデルを適用する理由は、階差をとることによって非定常な時系列を定常な時系列に変換できる場合があるからです。

元の時系列
1回階差をとった時系列
2回階差をとった時系列

実際に、2回階差をとることで定常時系列に変換できていることが分かります。

トレンド成分と周期成分がある時系列データの予測

いよいよ下図の予測を実装してみます。

データ作成

まずは青色のtrueのデータを作ります。

np.random.seed(0)
n = 200
x = np.arange(n) / 10 + np.sin(2 * np.pi * np.arange(n) / 20) + np.random.normal(0, 0.2, n)

ライブラリを使って学習＆予測

160サンプルまでを使ってパラメータが「order=(1,1,0),seasonal_order=(1,1,0,20))」のモデルを学習させます。

k = 160
model = arima.model.ARIMA(x[:k], order=(1, 1, 0), seasonal_order=(1, 1, 0, 20)).fit()

続いて、161～200サンプルまでを予測します。オレンジ色のデータを求めています。その後、161サンプル目の予測値をprintします。

x_pred_mean = model.get_prediction(start=k, end=n-1).predicted_mean
print(x_pred_mean[0])
# 15.666874644120792

ライブラリを使わずに予測

まず、modelの中身を見てみましょう。

print(model.summary())
#  coef      std    err     z   P>|z|  [0.025   0.975]
# ar.L1     0.5291 0.077 -6.871 0.000  -0.680  -0.378
# ar.S.L20 -0.4154 0.077 -5.424 0.000  -0.566  -0.265
# sigma2    0.0952 0.012  8.099 0.000   0.072   0.118

次に、このパラメータを使って、ライブラリを使わずに161サンプル目の予測をしてみましょう。

試行錯誤の結果、次の通り実行すればライブラリと同じ結果が得られることが分かりました。

arL1 = model.params[0]
arSL2 = model.params[1]
sigma2 = model.params[2]
i = k
print(
    x[i - 1] + arL1 * (x[i - 1] - x[i - 2]) + ((x[i - 20] - x[i - 21]) - arL1 * (x[i - 21] - x[i - 22]))
    +  arSL2 * (((x[i - 20] - x[i - 21]) - arL1 * (x[i - 21] - x[i - 22])) - ((x[i - 40] - x[i - 41]) - arL1 * (x[i - 41] - x[i - 42])))
)
# 15.666874644123475

つまり、$i-1$、$i-2$、$i-20$、$i-21$、$i-22$、$i-40$、$i-41$、$i-42$サンプル目のデータを使って$i$サンプル目を予測しているのです。

予測式を解説します。上のコードによると、予測式は次の通りです。だいぶ複雑になってしまいました。

x_i = x_{i-1} + \alpha(x_{i-1} - x_{i-2}) + \{(x_{i-20}  - x_{i-21})-\alpha(x_{i-21} - x_{i-22})\} + \beta[\{(x_{i-20} - x_{i-21})-\alpha(x_{i-21} - x_{i-22})\} - \{(x_{i-40}  - x_{i-41})-\alpha(x_{i-41} - x_{i-42})\}]

この式がどうやって得られるかを示します。

元の時系列(※1)	1回階差	1次のARモデルで予測(※2)	1回階差-1回階差のAR予測(※3)
$x_i$	①$x_i - x_{i-1}$	④$\alpha(x_{i-1} - x_{i-2})$	⑦$(x_i - x_{i-1})-\alpha(x_{i-1} - x_{i-2})$
$x_{i-1}$	②$x_{i-1} - x_{i-2}$	⑤$\alpha(x_{i-2} - x_{i-3})$	⑧$(x_{i-1} - x_{i-2})-\alpha(x_{i-2} - x_{i-3})$
$x_{i-2}$	③$x_{i-2} - x_{i-3}$	⑥$\alpha(x_{i-3} - x_{i-4})$	⑨$(x_{i-2} - x_{i-3})-\alpha(x_{i-3} - x_{i-4})$
$\vdots$	$\vdots$	$\vdots$	$\vdots$

上表の補足
- (※2)はorder=(1,1,0)のモデルによる予測（トレンドを予測するモデル）
  - ①を②から予測したものが④。時刻$i$のサンプルを時刻$i-1$のサンプルを使って予測している。
  - ②を③から予測したものが⑤。時刻$i-1$のサンプルを時刻$i-2$のサンプルを使って予測している。
- (※3)はorder=(1,1,0)のモデルによる予測で表現しきれなかった残りの成分
  - ⑦=①-④
  - ⑧=②-⑤
  - ⑨=③-⑥

order=(1,1,0)のモデルで予測しきれなかった残りの成分は周期性を加味して予測します。周期が20サンプルだと分かっていたとすると、20サンプル毎のデータを使ってARモデルを予測させます。

20サンプル毎の(※1)	20サンプル毎の(※3)
$x_i$	➉$\{(x_i - x_{i-1})-\alpha(x_{i-1} - x_{i-2})\} - \{(x_{i-20} - x_{i-21})-\alpha(x_{i-21} - x_{i-22})\}$
$x_{i-20}$	⑪$\{(x_{i-20} - x_{i-21})-\alpha(x_{i-21} - x_{i-22})\} - \{(x_{i-40} - x_{i-41})-\alpha(x_{i-41} - x_{i-42})\}$
$x_{i-40}$	⑫$\{(x_{i-40} - x_{i-41})-\alpha(x_{i-41} - x_{i-42})\} - \{(x_{i-60} - x_{i-61})-\alpha(x_{i-61} - x_{i-62})\}$
$\cdots$

「20サンプル毎の(※3)」に対して、定数項がない1次のARモデルを適用します。実際に、上表の⑪を使って1次のARモデルで➉を予測する式を構成してみます。

➉ = \beta\times ⑪ \\
\\
\Leftrightarrow  \{(x_i - x_{i-1})-\alpha(x_{i-1} - x_{i-2})\} - \{(x_{i-20}  - x_{i-21})-\alpha(x_{i-21} - x_{i-22})\} =\beta[\{(x_{i-20} - x_{i-21})-\alpha(x_{i-21} - x_{i-22})\} - \{(x_{i-40}  - x_{i-41})-\alpha(x_{i-41} - x_{i-42})\}]
\\
\\
\Leftrightarrow  x_i = x_{i-1} + \alpha(x_{i-1} - x_{i-2}) + \{(x_{i-20}  - x_{i-21})-\alpha(x_{i-21} - x_{i-22})\} + \beta[\{(x_{i-20} - x_{i-21})-\alpha(x_{i-21} - x_{i-22})\} - \{(x_{i-40}  - x_{i-41})-\alpha(x_{i-41} - x_{i-42})\}]

これは、コードの予測式と一致します。

ライブラリに与えたパラメータ「order=(η,φ,0),seasonal_order=(γ,1,0,ξ))」のうち、ηは1番目のARモデルの次数、βは1番目のARモデル適用前に何回階差をとるかを示します。γは2番目のARモデルの次数、ξは2番目のARモデル適用前に何サンプル毎にデータを間引くかを表します。

補足

ちなみに、冒頭の図はこのコードで書くことができます。

import numpy as np
import statsmodels.api
import statsmodels.tsa.arima as arima
import matplotlib.pyplot as plt

np.random.seed(0)
n = 200
k = 160
x = np.arange(n) / 10 + np.sin(2 * np.pi * np.arange(n) / 20) + np.random.normal(0, 0.2, n)
model = arima.model.ARIMA(x[:k], order=(1, 1, 0), seasonal_order=(1, 1, 0, 20)).fit()

x_pred_mean = model.get_prediction(start=k, end=n-1).predicted_mean
x_pred_std = model.get_prediction(start=k, end=n-1).se_mean
x_pred_interval_upper = x_pred_mean + 2.58 * x_pred_std
x_pred_interval_lower = x_pred_mean - 2.58 * x_pred_std
y_max = max(np.concatenate([x_pred_interval_upper, x]))
y_min = min(np.concatenate([x_pred_interval_lower, x]))

plt.plot(range(len(x)), x, color='blue', label='true')
plt.plot(range(k, len(x)), x_pred_mean, color='orange', label='predict')
plt.fill_between(range(k, len(x)), x_pred_interval_lower, x_pred_interval_upper, color='silver')
plt.vlines(k - 0.5, y_min, y_max, color='black', linestyles='dashed')
plt.legend(loc='upper left')
plt.show()