はじめに
以前k-meansをnumpyで実装するということをやりましたが、今回は時系列モデルの基礎としてARモデルを実装してみることにします。
以前一度やっているのですが、中途半端極まりなかったのでやり直し、というところです。
ちなみに今回は経済・ファイナンスデータの軽量時系列分析を勉強しながら書いたものなので、数式などの詳細はそちらを読むこと。
やってみる
データの生成
とりあえずAR(1)に従うデータを適当に作ってみます。
この時系列を表す式はとすると、平均
となるので、初期値は一旦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)
モデルの推定
では続いて、上で生成したデータを観測値としてモデルの推定、つまりの
を求めていきます。最小二乗法や最尤法を使うのですが、ここでは簡単な最小二乗法を利用した正規方程式から求める方法を使います。
mean = observed.mean() phi = np.sum((observed - mean)[1:] * (observed - mean)[:-1]) / np.sum((observed - mean)[:-1] ** 2) c = mean - phi * observed[:-1].mean()
単純ですね。今回の場合はと推定できました。
推定
このモデルを使って推定を行うことができます。
今回は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()
青い点が観測値、赤い点が予測、赤い線は予測と観測値の誤差です。
予測がピッタリ当たるなんてことはめったに無いですが、その誤差をどう表現するかと言うことは割と重要な問題です。
ピッタリ当たらないのであれば、だいたいどのくらいの範囲で当たるのか、ということを考えるのが普通なので。
信頼区間の計算
こんなときは信頼区間を計算します。
まず、先程見た観測値がある区間について。
ARモデルは正規分布に従うので、95%信頼区間は次のように表せます。
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()
これはつまり、ある時刻t+1のデータをモデルから1つ値を取り出したとき、それは点推定ですが本来はそこに正規分布に従う誤差(ノイズ)が乗ります。正規分布なので論理的にはから
の範囲を取りえますが、その確率を踏まえて、100回その施行を繰り返したとして95回くらいはこの区間内に収まるだろう、ということを表しています。
未来の予測と信頼区間
未来の予測はちょっとだけ手がかかります。
というのも、AR(1)の場合は1つ先の点を予測するにはその前のデータが必要なので、今回で言えばの最後の点を使って
を求め、さらに
というのはもちろん可能なのですが、単純に実行すると、以下のように平均に回帰してしまいます。
これはこういうものなので仕方ないですが、これでは将来の予測が立たなくて困るので、実際にはモデルから1つデータを取り出して、それに誤差を載せることで将来を予測します。
なんとなく平均(=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()
ここでも重要なのはこの信頼区間は真の平均値を95%の確率で含むと言ったことを言っているのではなく、モデルから同様の未来予測(施行)したときにまぁ100回中95回位はこの範囲に収まるだろうという範囲です。
先程の点での予測と違い、将来についてもある程度の傾向が見通せるようになったかなと思います。
終わりに
今回はARモデルを自前で書き起こして、中身を眺めてみました。
(完全にあっている自信は微妙ですが)とはいえなんとなくで理解していたよりは95%信頼区間の表しているものなどの意味もよく捉えられたと思います。
もしなにか気づいた点あればガッツリしてきくださると嬉しいです。