コードで理解するAutoRegressive モデル - えんじにあのじゆうちょう

はじめに

以前k-meansをnumpyで実装するということをやりましたが、今回は時系列モデルの基礎としてARモデルを実装してみることにします。

以前一度やっているのですが、中途半端極まりなかったのでやり直し、というところです。

ちなみに今回は経済・ファイナンスデータの軽量時系列分析を勉強しながら書いたものなので、数式などの詳細はそちらを読むこと。

やってみる

データの生成

とりあえずAR(1)に従うデータを適当に作ってみます。

この時系列を表す式は $y_{t+1} = 0.5 + 0.9 y_t + \epsilon$ とすると、平均 $\mu = \frac{c}{1-\phi} = \frac{0.5}{1-0.9} = 5$ となるので、初期値は一旦5を選択し、とりあえず100点分くらい作ってみます。。

np.random.seed(16)
def y_next(y_t):
    return 0.5+ 0.9 * y_t + np.random.normal(loc=0, scale=1, size=1)

N = 100
observed = np.zeros(N)
observed[0] = 0.5 / (1- 0.9)
#observed[1] = observed[0] + np.random.normal(loc=0, scale=1, size=1)
for i in range(1, N):
    observed[i] = y_next(observed[i-1])#, observed[i-2])
plt.scatter(range(N), observed)
plt.xlim(0, 110)
plt.ylim(0, 12)

f:id:marufeuillex:20200328124532p:plain

モデルの推定

では続いて、上で生成したデータを観測値としてモデルの推定、つまり $y_{t+1} = c + \phi y_t + \epsilon$ の $c, \phi$ を求めていきます。最小二乗法や最尤法を使うのですが、ここでは簡単な最小二乗法を利用した正規方程式から求める方法を使います。

$\phi = \frac{\sum_{t=1}^{T}(y_t - \bar{y_t})(y_{t-1}-\bar{y_{t-1}})}{\sum_{t=1}^{T}(y_{y-1}-\bar{y_{t-1}})^2}$
$c = \bar{y_t} - \phi \bar{y_{t-1}}$

mean = observed.mean()
phi = np.sum((observed - mean)[1:] * (observed - mean)[:-1]) / np.sum((observed - mean)[:-1] ** 2)
c = mean - phi * observed[:-1].mean()

単純ですね。今回の場合は $\phi = 0.9083593539321086, c = 0.4821382670070955$ と推定できました。

推定

このモデルを使って推定を行うことができます。
今回はAR(1)モデルなので、データがある区間に関しては、1時点前のデータを使って次時点の予測ができます。

いっぱい書くと見にくいので、とりあえず最初の10点位を予測してみます。

preds = pred(observed[:9], phi, c)
plt.scatter(np.arange(1, 10), preds, label="preds", color="red")
plt.scatter(np.arange(10), observed[:10], label="observed", color="blue")
for i in range(1, 10):
    plt.plot([i, i], [preds[i-1], observed[i]], alpha=0.2, color="red")
plt.legend()

f:id:marufeuillex:20200328130943p:plain

青い点が観測値、赤い点が予測、赤い線は予測と観測値の誤差です。
予測がピッタリ当たるなんてことはめったに無いですが、その誤差をどう表現するかと言うことは割と重要な問題です。
ピッタリ当たらないのであれば、だいたいどのくらいの範囲で当たるのか、ということを考えるのが普通なので。

信頼区間の計算

こんなときは信頼区間を計算します。
まず、先程見た観測値がある区間について。

ARモデルは正規分布に従うので、95％信頼区間は次のように表せます。

$\bar{X} \pm 1.96 \sqrt{MSE}$

MSE(Mean Squared Error)は観測値と実際の値の差の二乗をデータ数で割ったものです。
今回各時点には1つしかデータがないのでデータ数は1として計算して、先程のグラフに書き加えてみます。

preds = pred(observed[:9], phi, c)
plt.plot(np.arange(1, 10), preds, label="preds", color="red")
plt.scatter(np.arange(10), observed[:10], label="observed", color="blue")
mse = (preds - observed[1:10])**2
plt.fill_between(np.arange(1,10),  preds - np.sqrt(mse), preds + np.sqrt(mse),alpha=0.2, color="red", label="confidence")
plt.legend()

f:id:marufeuillex:20200328131853p:plain

これはつまり、ある時刻t+1のデータをモデルから1つ値を取り出したとき、それは点推定ですが本来はそこに正規分布に従う誤差(ノイズ)が乗ります。正規分布なので論理的には $-\infty$ から $+infty$ の範囲を取りえますが、その確率を踏まえて、100回その施行を繰り返したとして95回くらいはこの区間内に収まるだろう、ということを表しています。

未来の予測と信頼区間

未来の予測はちょっとだけ手がかかります。
というのも、AR(1)の場合は1つ先の点を予測するにはその前のデータが必要なので、今回で言えば $y_{99}$ の最後の点を使って $y_{100}$ を求め、さらに $y_{101}, y_{102}, \cdots$ というのはもちろん可能なのですが、単純に実行すると、以下のように平均に回帰してしまいます。

f:id:marufeuillex:20200328132342p:plain

これはこういうものなので仕方ないですが、これでは将来の予測が立たなくて困るので、実際にはモデルから1つデータを取り出して、それに誤差を載せることで将来を予測します。

f:id:marufeuillex:20200328133912p:plain

なんとなく平均(=5)方向に回帰しているものの、ばらついているのがわかるかと思います。今度はこのばらつきを範囲で表そうとしてみます。

ただし今回は観測値がないので先程の計算はそのまま使えません。そこで、先程の予測のときのようにM個の乱数(Normal(0,1))を用意して、それを予測値に足し合わせ、それを観測値として扱うことでMSEを算出します。

N = 100
M = 5
predicted_samples = np.zeros(M*20).reshape(M, 20)
predicted_samples[:, 0] = observed[-1]
for i in range(1, 20):
    predicted_samples[:, i] = c + phi * predicted_samples[:, i-1] + np.random.normal(loc=0, scale=1, size=(M))
predicted_samples = predicted_samples[:, 1:]

predicted_mean = predicted_samples.mean(axis=0)
MSEs = np.sum((predicted_samples-predicted_mean) ** 2, axis=0) / M


lower = predicted_mean - 1.96 * np.sqrt(MSEs)
mean= predicted_mean 
upper = predicted_mean+ 1.96 * np.sqrt(MSEs)

plt.scatter(range(N), observed, color="blue", label="observed")
plt.fill_between(np.arange(N, N+19), lower, upper, color="red", label="95%", alpha=0.1)
plt.plot(np.arange(N, (N+19)), mean, color="red", alpha=0.3, label="predict_mean")

plt.xlim(90, 120)
plt.ylim(0, 12)
plt.legend()