統計ソフトRのパッケージMASSに含まれているSP500の収益率のデータを使う
ret <- ts(SP500)
としてRコードで分析を進める。
ugarchでegarch(1,1)モデルを推定するEＧＡＲＣＨボラティリティモデ
ルをパッケージ rugarchを使って推定してみる。詳しい使い方については
Introduction to the rugarch package (Version 1.3-8)
Alexios Ghalanos に書かれている。cran.r-project.orgの
Vignettesnで読めるようになっている。
EGARCH モデルの概要
パッケージrugarchではeGARCH(1,1) を下記のような条件付分散の式で示される。
\(e_{t}\text{=}\sqrt{h_{t}}z_{t}\) 従って \(z_{t}\text{=}\frac{e_{t}}{\sqrt{h_{t}}}\)
\(z_{t}\ \ \ \ \ iid\ \ N(0,1)\) 期待値０、分散１の標準正規変数
\(ret_{t}=\mu+e_{t}\)
\(\ln(h_{t})=\omega+\alpha_{1}z_{t-1}+\gamma\left\{ \left(\left|z_{t-1}\right|-E\left(\left|z_{t-1}\right|\right)\right)\right\} +\beta_{1}\ln(h_{t-1})\)
ここで注目すべきは、この条件付分散の式に（標準正規変数の絶対値）の 期待値が組み込まれている点である。絶対値を利用して悪いニュースは良いニュースよりもボラティリティを増幅さすという非対称性をモデル化しようとしている。
egarchspec <- ugarchspec(variance.model=list(model="eGARCH",#各係数の推定値は以下のようになる。
## Estimate Std. Error t value Pr(>|t|)
## mu 3.278995e-02 0.011463333 2.860420538 4.230796e-03
## omega -1.182441e-05 0.002252076 -0.005250451 9.958108e-01
## alpha1 -8.274741e-02 0.010288179 -8.042960453 8.881784e-16
## beta1 9.822487e-01 0.002156778 455.424079226 0.000000e+00
## gamma1 1.260175e-01 0.012893593 9.773653326 0.000000e+00
# μ= 0.0327 ω= -0.00001 α=-0.0827 β=0.9822 γ=0.126
#と推定されている。
rhat <- egarchfit@fit$fitted.values
plot.ts(rhat,ylim=c(-0.2,0.2))
#収益率のグラフはretの平均値0.0327で一定となっている。
hhat <- ts(egarchfit@fit$sigma^2)
#
par(mfrow=c(2,1))
plot.ts(hhat)
plot.ts(ret)
# 収益率が大きなマイナスになった時には大きなプラスの収益率の時よりも
# 相対的にボラティリティが大きくなる傾向が見られる。
par(mfrow=c(1,1))
#
#
plot(egarchfit,which="all")
##
## please wait...calculating quantiles...
パッケージrugarchではeGARCH(1,1) を下記のような条件付分散の式で示される。
\(e_{t}\text{=}\sqrt{h_{t}}z_{t}\) 従って \(z_{t}\text{=}\frac{e_{t}}{\sqrt{h_{t}}}\)
\(z_{t}\ \ \ \ \ iid\ \ N(0,1)\) 期待値０、分散１の標準正規変数
\(ret_{t}=\mu+e_{t}\)
\(\ln(h_{t})=\omega+\alpha_{1}z_{t-1}+\gamma\left\{ \left(\left|z_{t-1}\right|-E\left(\left|z_{t-1}\right|\right)\right)\right\} +\beta_{1}\ln(h_{t-1})\)
\(\ln(h_{t})=-0.00001-0.082747z_{t-1}+0.12602\left\{ \left(\left|z_{t-1}\right|-E\left(\left|z_{t-1}\right|\right)\right)\right\} +0.98225\ln(h_{t-1})\)
ここで注意を要するのは、この条件付分散の式には標準正規変数の絶対値の
期待値\(\mathrm{\text{E}}\left(\left|\mathrm{z}_{\mathrm{t}-1}\right|\right)\)が組み込まれている点である。
\(\mathrm{\text{E}}\left(\left|\mathrm{z}_{\mathrm{t}-1}\right|\right)\)は\(\sqrt{\frac{2}{\pi}}\simeq\) 0.7978$
いう定数になるため、上式を書き換えると
\(\ln(h_{t})=-0.00001-0.082747z_{t-1}+0.12602\left\{ \left(\left|z_{t-1}\right|-0.79788\right)\right\} +0.98225{\ \ln(h}_{t-1})\)
\(=(-0.00001-0.12602\times0.79788)-0.082747z_{t-1}+0.12602\cdot\left|z_{t-1}\right|+0.98225\ln(h_{t-1})\)
\(=-0.1005-0.082747z_{t-1}+0.12602\cdot\left|z_{t-1}\right|+0.98225\ln(h_{t-1})\)
Rugarch以外のパッケージや他の統計ソフトで、もし条件付分散式を
\(\ln(h_{t})\text{=}\omega\text{+}\alpha_{1}z_{t-1}\text{+}\gamma\cdot\left|z_{t-1}\right|\text{+}\beta_{1}\ln(h_{t-1})\)として定数項を一本にまとめて推定して
Outputするような場合にはrugarchによるoutputと定数項の数値が0.79788だけ
差異が出てくる。もちろん式を整理すれば最終的には大体同じような結果を得られると思う。
やについても条件付分散式の算式が異なると係数の推定値も異なってくるので使用
しているソフトがどのような件付分散式を前提としているか確かめておく必要があるだろう。
\(\ln(h_{t})=-0.1005-0.082747z_{t-1}+0.12602\cdot\left|z_{t-1}\right|+0.98225\ln(h_{t-1})\)
\(z_{t-1}\geq0\) ならば\((\alpha+\gamma)z_{t-1}\) つまり（-0.082747+0.12602\(\fallingdotseq\)0.043）
\(z_{t-1}<0\) ならば\(\left(\alpha-\gamma\right)z_{t-1}\) つまり（-0.082747-0.12602\(\fallingdotseq\)-0.2087）
\(\ln(h_{t})=-0.1005+0.98225\ln(h_{t-1})+\begin{bmatrix}0.043z_{t-1}\\ -0.2087z_{t-1} \end{bmatrix}\ \ \ \ \ \ \ \ \ \ \begin{matrix}z_{t-1}\\ z_{t-1} \end{matrix}\ \begin{matrix}\geq0\\ <0 \end{matrix}\)
指数関数表示にすれば
\(h_{t}=\ h_{t-1}^{0.98225}\centerdot e^{-0.1005}\times\begin{bmatrix}e^{0.043z_{t-1}}\\ e^{-0.2087z_{t-1}} \end{bmatrix}\ \ \ \ \ \ \ \ \ \ \begin{matrix}z_{t-1}\\ z_{t-1} \end{matrix}\ \begin{matrix}\geq0\\ <0 \end{matrix}\)
\(E\left(\left|x\right|\right)=\sqrt{\frac{2}{\pi}}\)の求め方
xを標準正規変数とすると、その絶対値の期待値は、期待値の定義から
\(E(x)\text{=}\frac{1}{\sqrt{2\pi}}\int^{\infty}_{-\infty}\left|x\right|\cdot e^{-\frac{x^{2}}{2}}dx\text{=}\frac{2}{\sqrt{2\pi}}\cdot\int^{\infty}_{0}x\cdot e^{-\frac{x^{2}}{2}}dx\)
ここで\(\text{-}\frac{x^{2}}{2}\)=tと置き換えると
\(-xdx\text{=}dt\)を使い、\(\int^{\infty}_{0}x\cdot e^{-\frac{x^{2}}{2}}dx\text{=}\int^{\infty}_{0}(-e^{t})dt\)
\(\sqrt{\frac{2}{\pi}}\cdot\int-e^{-\frac{x^{2}}{2}}dx\)
\(\text{=}\sqrt{\frac{2}{\pi}}\)[0-(-1)]=\(\sqrt{\frac{2}{\pi}}\)
単位根検定の背景（ディッキー・フラー分布のシミュレーション with R）