えんじにあのじゆうちょう

勉強したことを中心にアウトプットしていきます。

コードで理解するARモデルの最尤推定

はじめに

前回に引き続き、ARモデルについて書きます。
t.marufeuille.dev

今回は最尤推定でのモデルのfitingを試していきます。

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

やってみる

最尤推定法

最尤推定は以前以下に書きましたが、ヒトコトで雰囲気を伝えるならf(x, \theta)を確率密度関数とすると、得られているデータ x_iに対して一番もっともらしいパラメータ\thetaを求めるということです。
t.marufeuille.dev

xを得られた値,  \thetaをその確率密度関数のパラメータとした時、尤度関数

 L(\theta) = \prod f(x_i; \theta)

を最大化するような \thetaを求めれば良いですが、このままだと計算しにくいので対数をとった対数尤度関数を使います。

 l(\theta) = \sum \log f(x_i; \theta)

今回確率密度関数が何になるかですが、ARモデルなので y_{t-1}, y_{t-2}, \cdots, y_{t-p}を元に y_tが得られる確率です。

今回以下で扱うのは前回同様AR(1)モデルなので、 \epsilon \sim Normal(0, \sigma^2)のとき、 y_t | y_{t-1} \sim Normal(c + \phi y_{t-1}, \sigma^2)に従うらしいので、これを f(x_i; \theta)として、式を整理すると、結局以下を最大化すること同義になります。
 -(y_t - (c + \phi y_{t-1}))^2

なので、これを解いて最大値を見つけることで最も良いらしいパラメータを導出することができます。

データについて

扱うデータは前回の記事同様ですが、以下に従い生成しています。

 y_t = 0.5 + 0.9y_{t-1} + \epsilon

念の為コードを張っておくとこんな感じです。

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)
for i in range(1, N):
    observed[i] = y_next(observed[i-1])

まずは力技で解いてみる

これを最大化するところを見つけるのに、まずは力技で総当り的に解いてみましょう。
とはいえ無限大に探すのは無理なので、以下のようにnumpyでmeshgridを作って当たりをつけて解を探します。

c_range = np.arange(0, 1, 0.01)
phi_range = np.arange(0, 1, 0.01)
def target(c, phi, y):
    return np.sum(-(y[1:] - (c + phi * y[:-1])) ** 2, axis=1)

c_, phi_ = np.meshgrid(c_range, phi_range)
max_likelyhood = 0
c, phi = 0, 0
idx = np.argmax(target(c_.reshape(c_.size,1) , phi_.reshape(c_.size, 1), observed))
c, phi = c_.reshape(-1)[idx], phi_.reshape(-1)[idx]

print(c, phi)
(0.47000000000000003, 0.91)

今回は c, \phiともに0 ~ 1の範囲を0.01刻みで探索しています。
本来この範囲はわからないのでずるいですね(笑)

微分して解いてみる

最大化する式を再掲します。
 -(y_t - (c + \phi y_{t-1}))^2

これよく見るとそれぞれのパラメータについて2次関数の形式になっています。

ということで、E資格の勉強で散々やったように微分して傾きがなる点に近づけていく、勾配降下法で解けるはずなので、やってみます。
それぞれ微分すると以下のようになるので、これでパラメータを更新すればよいですね。

 \frac{\partial f(x; \theta)}{\partial c} = 2(y_t - (c + \phi y_{t-1}))
 \frac{\partial f(x; \theta)}{\partial \phi} = 2y_{t-1}(y_t - (c + \phi y_{t-1}))


ではPythonで書いてみます。

phi, c = 0.5, 0.5
eta = 0.00001
for i in range(50000):
    dc = (np.sum(2 * (-observed[1:] + c + phi * observed[:-1])))
    dp = (np.sum(2 * observed[:-1] * (-observed[1:] + c + phi * observed[:-1])))
    c = c - eta * dc
    phi = phi - eta * dp
print(c, phi)
(0.47833273408680893, 0.9087585027352589)

ということでしっかり求められてますね。なお、phi, cは適当に与えています。

まとめ

ということで、今回は最尤推定でARモデルを推定してみました。
話がややずれますが、E資格(の講座)で嫌というほどやった勾配降下法を覚えていたおかげで案外楽に書くことができました。

次回はVARモデルのパラメータ推定をやってみようと思います。乞うご期待(?)