ドル円ボラティリティをGARCHで推定 その2(正規分布、t分布)

概要

・GARCH(1,1)パラメータ推定続き 正規分布とt分布の場合
・GARCHパラメータ標準誤差の計算
・fGarchパッケージとの推定結果の比較。概ね良さそうに見える

モデル

前回(http://nekopuni.holy.jp/?p=1146)で行なった推定をもう少し一般的に。
$$
y_t = \mu + \epsilon \\
\epsilon_t \sim iid, E \left( \epsilon_t \right) = 0, V \left( \epsilon_t \right) = \sigma_t^2 \\
\sigma_t^2 = \omega + \sum_{i=0}^p \beta_i \sigma_{t-i}^2 + \sum_{j=1}^q \alpha_j \epsilon_{t-j}^2
$$

\( \epsilon_t \)が正規分布だったのを、iidという仮定にゆるめます。今回は正規分布と、GARCH界隈では一般的なt分布の2パターンでのパラメータ推定を行います。

Rコード

データはドル円。サンプル数が多いと推定に時間がかかるので、2013年以降を使います。

Japan / U.S. Foreign Exchange Rate – FRED – St. Louis Fed
http://research.stlouisfed.org/fred2/series/DEXJPUS

変数yにドル円を格納しているとして、

という感じ。l.normが正規分布のもとでの対数尤度関数を計算する関数、l.tがt分布のもとでの対数尤度関数。
ちなみに一般化t分布(という名前があるのかは知りませんが)の密度関数はdstd関数を使っておりまして、これはfGarchパッケージにて定義されているものです。
非負制約のあるパラメータにはexpをかまして、t分布の自由度には推定結果が2以上になるように、関数eをかましています。

で、推定された結果が以下

上から順に\( \mu, \omega, \alpha, \beta \)でt分布の場合は、自由度がpar[5]に入っています。

次に推定係数の標準誤差の計算。
この計算は参考文献[5]のP33あたりを参考ということで。
対数尤度関数を推定パラメータで数値的に微分したヘシアンを使うと計算できるようです。

ただし、非負制約回避のためにexpをかましている部分があるためにそこだけゴニョゴニョと変数変換。
すると以下の様な値が出てきます。

分布を変えることによるボラティリティの違い

として、適当にプロットしてみると以下のようになります。
直近は正規分布とt分布にズレが出てきてて、t分布の方が上に来ていますね。

140721_1

fGarchパッケージとの比較

tseriesパッケージのgarch関数は正規分布にしか対応していなさそうなので、今回はfGarchパッケージを使用して同様のことをしてみます。

とすると簡単に推定をすることができまして、

ということで、大体あっていそうという感じでしょうか。

メモ

・推定に用いる系列は100倍して単位を%にしたほうが良さそう。
100倍しない場合、桁落ちか何かのせいで最適化が正常に終わらないことが多々^^;;

・optimにてmaxitをそれなりに大きい数字(今回は1000)にしたほうが良い。Nelder-Meadではmaxitのデフォルト値は500なのですが、これだとうまく収束しきれないのか何か知らんがfGarchの推定結果と遠い感じに。あまり深くは立ち入らないことにします笑

参考文献

T. Bollerslev: A Conditionally Heteroskedastic Time Series Model for Speculative Prices and Rates of Return (1987)
http://www.hss.caltech.edu/~camerer/SS280/BollerslevRES87.pdf

F. Klaassen: Improving GARCH Volatility Forecasts with Regime-Switching GARCH (2001)
http://dare.uva.nl/document/1743

package: fGarch
http://cran.r-project.org/web/packages/fGarch/fGarch.pdf

Student’s t-distribution – Wikipedia, the free encyclopedia
http://en.wikipedia.org/wiki/Student’s_t-distribution

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

R5クラスの定義+オプションGreeks(ベガ周りを主に)の計算

概要

・バニラオプションのGreeks(特にベガ周り)を知っておくと良い気がしてきたので個人的に気になる点のまとめ
・ついでにRでオブジェクト指向チックに書く練習(R5クラスっていうんですかね)

結果

Rのコードは長いので後掲するとして個人的に気になる結果を貼り付ける。
特にベガ周りに触るのは初めてなので違ってたら指摘してくださると有難いです。
対象はプレーンなバニラコールロングです。

1. ボラティリティ水準とガンマの比較
→ ATM付近ではガンマとボラティリティは逆の関係。ボラティリティが上昇するとガンマが下がる
140710_1

2. ボラティリティとオプション・プレミアムの関係
→ これは簡単に1対1の関係。ボラティリティが上昇するとプレミアムも上昇
140710_2

3. ATMFにおけるボラティリティとベガの関係
→ 一定レベル以上でのボラティリティの変動はベガにあまり影響を与えない
→ 低ボラティリティ環境下においては1ベガあたりのPLへの影響が大きい(本当か?)
140710_3

4. OTMにおけるボラティリティとベガの関係
→ ボラティリティが上昇するとベガも上昇。DeepOTMになるほどその影響は小さくなる
140710_4

5. 異なる原資産水準でのATMFの行使価格とベガの関係
→ 原資産の水準が高いほどベガは大きくなる
140710_5

6. 原資産価格とボルガ(ボラティリティの2階微分)の関係
→ 3で見たように、ATM近辺ではボラティリティが上昇してもベガはあまり変化しない(ボルガ≒0)
→ 4で見たように、ATM以外の領域ではボルガが正 = ボラティリティが上昇するとベガも上昇する
140710_6

7. 異なるボラティリティ水準でのボルガの変化
→ ATMではボラティリティ水準によらずボルガはほぼゼロ
→ OTMは知らん。というかボルガってオプショントレーダーでもあるまいし普通見るんかね。よう知らんけど。
140710_7

8. 原資産とゾンマの関係
→ あんのかなーと思ったら案の定あった謎のGreeksゾンマ。ガンマをボラティリティで微分したもの。
→ 1で見たように、ATM付近ではゾンマが負なのでボラティリティとガンマが逆関係。ATMから外れるとゾンマがプラスになり、ボラティリティが上昇するとガンマも上昇する。
140710_8

グリークスの名前メモ
$$
{\rm vega} = \frac{{\partial }V}{{\partial } \sigma} \\
{\rm volga} = \frac{{\partial }^2 V}{{\partial } \sigma^2} \\
{\rm zomma} = \frac{{\partial } \Gamma}{{\partial } \sigma} = \frac{{\partial }^3 V}{{\partial } S^2 {\partial } \sigma}
$$

他にも面白いGreeksがあるみたいなんで興味がある方はWikipedia見ると面白いと思います。
正直3次以降のGreeksをどう解釈したら良いのかよく分かりませんが。

Rコード

感度分析を行うとオブジェクトのパラメータを上書きしてしまうのが致命的ですがこのクラスは練習用なのでそのままにしておきます。

参考文献

Greeks (finance) – Wikipedia, the free encyclopedia
http://en.wikipedia.org/wiki/Greeks_(finance)

Rのオブジェクト指向について(R) – script of bioinformatics
https://sites.google.com/site/scriptofbioinformatics/r-tong-ji-guan-xi/rnoobujekuto-zhi-xiangnitsuite-r

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