先日、こんな論文を見つけて、なんかおもしろそうなのでちょっと試してみることにしました。
ちなみに著者はパッケージ開発者ご本人のようです。
手順としては普通に多変量自己回帰モデルを推定した後にその残差を動的条件付相関係数型GARCH（DCC-GARCH）モデル化するという流れです。
論文に書かれた例題は、まぁそのままRに入力していけば計算できるのですが、他のデータでは？という時にどうも私にはよくわからないことが。
5ページ目、いざDCC-GARCHの推定というところで、
「dcc.estimation関数であるが，この関数を利用するためには，尤度関数最大化のために利用されるパラメータの初期値を設定しなければならない．」
とあり、例題ではいくつかのベクトルと行列に「初期値」が代入されているのですが、この値がどこから来たものなのかがさっぱり・・・。
これはもう当然私の不徳の致すところ・・・と思いつつ、開発者が日本の方なのをいいことに直接お尋ねしたところ、早速お返事を下さいました。
どうもありがとうございます m(_ _)m
ここで必要な初期値というのは推定する条件付分散式と条件付相関係数式の尤度関数を最大化する際の初期値で、
「関数最適化の際の初期値については、optim関数のヘルプをご覧いただくのが良いと思います。」
なるほど、たしかに内部の実装を辿っていくとoptim()で最適化するようですが、
しかしこれは・・・最適解のあたりをつけなくてはならないということでしょうか？
それとも必要な要素数と行列の形状さえ満たしておけば何でもいいからとりあえず実数を放り込んでみるとか？
ちょっとこの辺知識が追い付いてなくてよくわからないです・・・。
・・・と、そんな私にも朗報！
「単にDCC-GARCHモデルを推定したいだけであれば、ccgarchをアップデートしたccgarch2というパッケージを用意してあります。ccgarch2は、モデル推定の際に初期値を与える必要はありません。CRANには公開していませんので、以下のコマンドでインストールしてください。」
わーい！
install.packages("ccgarch2",
repos = c("http://R-Forge.R-project.org", "http://cran.rstudio.com/"),
dep = TRUE)
#Mac or Linuxならtype = "source"を追加
該当する関数の引数を見てみると、
> args(estimateDCC)
function (inia = NULL, iniA = NULL, iniB = NULL, ini.dcc = NULL,
data, model = "diagonal", ...)
となっているので、初期値を指定することもできるようですが、どうやら内部でRsolnpパッケージと連係して初期値が与えられていなければランダムに(？)割り振って最適化してくれるようです。
ここまでの分析の流れは論文を参照してもらうとして、とりあえず新関数の引数にデータだけ与えて実行してみると・・・
dcc.results2 <- estimateDCC(data=eps)
summary(dcc.results2)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
mu1 0.015314 0.765497 0.0200 0.9840
mu2 0.002755 0.022497 0.1225 0.9025
mu3 0.006707 0.041638 0.1611 0.8720
a1 0.096745 0.095446 1.0136 0.3108
a2 0.077181 0.014088 5.4786 <2e-16 ***
a3 0.034602 0.006505 5.3188 <2e-16 ***
A11 0.083675 0.052808 1.5845 0.1131
A22 0.098773 0.013082 7.5502 <2e-16 ***
A33 0.060483 0.009412 6.4259 <2e-16 ***
B11 0.898779 0.068897 13.0452 <2e-16 ***
B22 0.888066 0.011462 77.4759 <2e-16 ***
B33 0.931986 0.009930 93.8534 <2e-16 ***
alpha 0.024019 0.005432 4.4221 1e-05 ***
beta 0.959803 0.012868 74.5890 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Number of Obs.: 5909
Log-likelihood: -33623.12032
Information Criteria:
AIC: 67274.241
BIC: 67367.820
CAIC: 67381.820
Optimization method: SQP by Rsolnp package
Convergence (1st, 2nd): 0 0
Iterations (1st) : 2 682
Iterations (2nd) : 2 57
無事に係数等が出力されました！
ところで同じデータ（論文とほぼ同じであろう株価データを作りました）で論文どおりの初期値による前バージョンでは、
a1 a2 a3
estimates 0.07604341 0.096108553 0.032569004
std.err 0.01632648 0.008005329 0.008501938
A11 A22 A33
estimates 0.07341938 0.10833925 0.05732549
std.err 0.01807429 0.01194691 0.01247122
B11 B22 B33
estimates 0.912673817 0.874512431 0.935390056
std.err 0.009697383 0.006827305 0.007863965
dcc alpha dcc beta
estimates 0.02296328 0.961902952
std.err 0.00211921 0.004695311
となっており、結果が異なります・・・。
これは初期値が異なるため結果の収束先が変化しているのか、推定方法の違いなのか？
見比べるとmuの有無とか、a1とa2が入れ替わってるようなのとか気になりますね・・・。
それぞれ相関係数をプロットして見比べてみると、動きが大きく変わるようなことはないけどややオフセットする感じでした。
試しに新バージョンに同じ初期値を与えてみたところ、
[1] "Bad initial values. Try again without initial values."
と怒られてしまいました。
最適化方法が変わったことにより初期値のレンジが制限されているのかも？
また、これもランダムな初期値の影響と思われますが、新バージョンでは実行毎に出力結果が微妙に変化しますので、この辺の解釈をどうすべきかも私的学習課題です。
（あと普通にVARを推定したその残差にGARCHをあてはめるのって・・・？とか）
さしあたり動くことが確認できただけ、という非常に中途半端なところですが、ここまでの経緯の中で関数のパラメーターの最適化に思いを馳せたり、optim関数やRsolnpパッケージの存在を知ったことも収穫でした。
開発者、中谷朋和さんにあらためましてお礼申し上げます。
どうもありがとうございました m(_ _)m
※中谷さんより補足説明を頂きましたのでコメント欄に追記しました。