GARCHパラメータ推定

概要

・GARCHのパラメータを自力で最尤推定してみる
・tseriesのgarch関数とほぼ同じような推定を得る

GARCH尤度関数

前回のポスト(http://nekopuni.holy.jp/?p=1146)ではtseriesのgarch関数を使えばとりあえずモデルに当てはめてくれるということが分かりましたが、イマイチどのようにして推定されているのかがわからなかったため、ゼロから尤度関数を定義して最尤推定してみます。
GARCH(1,1)を想定した場合、尤度関数\( L \)は

$$
L = \prod_{t=1}^T \frac{1}{\sqrt{2 \pi \sigma_t^2}} {\rm exp}\left( -\frac{\epsilon_t^2}{2\sigma_t^2} \right) \\
\sigma_t^2 = w + \beta \sigma_{t-1}^2 + \alpha \epsilon_{t-1}^2 \\
\epsilon_0^2 = \sigma_0^2 = \frac{1}{T}\sum_{t=1}^T \epsilon_t^2
$$

という感じで定義されます。最終的には対数尤度を取ってそれを最大化するパラメータ\( w, \beta, \alpha \)を推定するという流れですね。

Rコード

とりあえず動けば良いやということで

対象データはSPを使っています。インデックス値が変数spに代入されています。

S&P 500© – FRED – St. Louis Fed
http://research.stlouisfed.org/fred2/series/SP500

とりあえず確認ということで、tseriesを使ったボラティリティを確認。
プロットを重ねて描画する方法は大仏様のを写経(参考文献2)
そしてプロットしたのが以下になります。

140705_1

よく言われるように、株価が下落するときにボラティリティは上昇する傾向にありますな。
で、これをゼロから構築すると以下

sp.garch.l関数で対数尤度を計算するようにして、それをoptim関数で最大化させています。
一つあるとすれば、推定されるパラメータは正でなければなりませんので、普通なら制約条件付きの最大化となりますが、
ここではexpをかませることで制約条件を回避しています。

推定されたパラメータを比較してみると以下

上がtseriesによる推定で、下が今回の自力での推定。
完全に一致というわけにはいきませんが、まあまあ同じような値。

140705_2

それぞれのボラティリティをプロットしてみるとほぼ一致しているのがわかります。初期値近くで若干ずれが発生していますが、これは恐らく初期値をどう置くかによる誤差ですかねえ。

140705_3

ちなみに差分を取ってみると、初期値近くで誤差が発生して以降はほぼゼロとなっているようです。
ということで、恐らくという条件付きですが、tseriesのgarch関数は渡された系列を素直に残差系列とみなして、GARCHに当てはめていそうということが分かりました。
したがって例えばARMA-GARCHモデルなどを適用したい場合には、まず原リターン系列に対してARMAを適用して残差系列(恐らくresidualsなどに格納される値)を算出し、それをgarchにぶっこむという流れになりますかな。

ちなみに

こんな回りくどいことをやらなくともtseriesのソース見に行けばよいのではと思われるかもしれないですが、実際見てみるとgarchのフィッティングを行っていると思われる関数はc言語で書かれており、それを見に行くとどうやらFortranに渡しているように見えて、そこで私は深入りするのをやめましたね(― ―)

参考文献

Amazon.co.jp: ボラティリティ変動モデル (シリーズ 現代金融工学): 渡部 敏明, 木島 正明: 本
http://www.amazon.co.jp/dp/4254275048

Rで2軸PLOT(y軸の左右に値を表記)する方法 – My Life as a Mock Quant
http://d.hatena.ne.jp/teramonagi/20120107/1325922512