(cache) ■単変量/多変量GARCHモデル

GARCHモデルの計算を行うRのパッケージはいくつかあるけど、その中でもRcpp実装のrugarchとrmgarchを自分はよく使っている。

これまでコードをきちんと管理してなくて使う度に再実装していたので、使いたい時にいつでも使えるようにメモっておく。

概要

単変量GARCHがrugarchで多変量がrmgarchという構成になっている。これらのパッケージはAlexios Ghalanosという人物によって

作成されており、使い方やモデル構成などにある程度の一貫性がある。

また、この人はR/Financeとかでも発表しているので動向をチェックしておくと良いかもしれない。

加えて、この人は非線形最適化パッケージRsolnpの作者でもあり、精力的にパッケージ開発を行っている。

インストール

options(repos="http://cran.md.tsukuba.ac.jp")
install.packages(c("rugarch", "rmgarch"))

コード自体は特に問題がないんだけど、CRANとか見てもわかるようにDependsとかImportsしてるライブラリが大量にあるため

それらがインストールされていない場合は一気にインストールされるので注意。

単変量GARCH: rugarch

このパッケージで計算できるモデルは以下のような形式。(ただし、正規分布以外のモデルも使用可能。)

$r_t = \mu + \sum_{j=1}^{p}\phi_j(r_{t-j} - \mu) + \epsilon_t + \sum_{j=1}^q\theta_j\epsilon_{t-j}$

$\epsilon_t = \sigma_t z_t,$ $z_t \sim N(0, 1)$

$\sigma_t^2 = \left(\omega + \sum_{j=1}^m \zeta_j v_{jt} \right ) + \sum_{j=1}^q\alpha_j\epsilon_{t-j}^2 + \sum_{j=1}^p\beta_j\sigma_{t-j}^2$

リターンの部分がARMA(p, q)モデルで条件付きボラティリティ部分が普通のGARCH(p, q)に外生項vを足したものとなっている。

ARMAの移動平均次数qは普通なら整数に限定されるが、設定を切り替えれば実数まで拡張したARFIMAでも計算可能。

パラメータを推定するだけならugarchspecでモデル設定してugarchfitで計算を実行するだけでOK。

(この流れは今は亡きRHmmパッケージを彷彿とさせる書き方になっている。)

結果に再現性がほしい時はfixed.pars=list()のリストの中身に固定パラメータ値を入れれば良い。

で、具体的には以下のように書けばいい。

library(rugarch)
data(sp500ret)  # 適当なリターンデータをパッケージから取得
# デフォルトの設定、通常はspec.default <- ugarchspec() でOK
spec.default <- ugarchspec(variance.model = list(model="sGARCH", garchOrder=c(1, 1), submodel=NULL, external.regressors=NULL), 
                           mean.model     = list(armaOrder=c(1, 1), include.mean=T, archm=F, archpow=1, arfima=F, external.regressors=NULL, archex=F),
                           distribution.model="norm", start.pars=list())
# リターンをARMA(0, 0)とした場合
spec.garch <- ugarchspec(variance.model = list(model="sGARCH", garchOrder=c(1, 1)),
                         mean.model     = list(armaOrder=c(0, 0), include.mean=T))
out.default <- ugarchfit(data=sp500ret, spec=spec.default)
out.garch <- ugarchfit(data=sp500ret, spec=spec.garch)

結果は何も指定しなければ以下のようにバーっと帰ってくるけど、

個別にいろいろ値を取りたいときはスロット使ってout.default@model$parsみたいに取り出せばいい。

> out.default

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(1,1)
Mean Model	: ARFIMA(1,0,1)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000523    0.000087   6.0067        0
ar1     0.870134    0.072210  12.0501        0
ma1    -0.897344    0.064640 -13.8821        0
omega   0.000001    0.000000   5.2346        0
alpha1  0.087715    0.007724  11.3563        0
beta1   0.904935    0.008483 106.6826        0

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000523    0.000099   5.3055 0.000000
ar1     0.870134    0.085993  10.1186 0.000000
ma1    -0.897344    0.077195 -11.6244 0.000000
omega   0.000001    0.000001   2.1344 0.032809
alpha1  0.087715    0.029913   2.9324 0.003364
beta1   0.904935    0.029361  30.8211 0.000000

LogLikelihood : 17902.41 

Information Criteria
------------------------------------
                    
Akaike       -6.4807
Bayes        -6.4735
Shibata      -6.4807
Hannan-Quinn -6.4782

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic   p-value
Lag[1]                      5.548 1.850e-02
Lag[2*(p+q)+(p+q)-1][5]     6.439 1.252e-05
Lag[4*(p+q)+(p+q)-1][9]     7.196 1.104e-01
d.o.f=2
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.104  0.2933
Lag[2*(p+q)+(p+q)-1][5]     1.499  0.7403
Lag[4*(p+q)+(p+q)-1][9]     1.957  0.9100
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]   0.01978 0.500 2.000  0.8882
ARCH Lag[5]   0.17495 1.440 1.667  0.9713
ARCH Lag[7]   0.53656 2.315 1.543  0.9750

Nyblom stability test
------------------------------------
Joint Statistic:  1.428
Individual Statistics:             
mu     0.2060
ar1    0.1477
ma1    0.1047
omega  0.1565
alpha1 0.1346
beta1  0.1134

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 1.49 1.68 2.12
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value      prob sig
Sign Bias           0.4296 6.675e-01    
Negative Sign Bias  2.9492 3.199e-03 ***
Positive Sign Bias  2.3921 1.678e-02  **
Joint Effect       28.9841 2.257e-06 ***


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     178.6    5.738e-28
2    30     188.1    3.136e-25
3    40     217.8    1.084e-26
4    50     227.9    7.064e-25


Elapsed time : 1.240071

ここではvariance.modelにsGARCHを指定したが、これがいわゆるvanillaなGARCHモデルを意味している。

他にも以下の拡張GARCHをvariance.modelに指定できる。

モデル名称	variance.modelに与える引数	paperへのリンク
GARCH	sGARCH	sGARCH
integrated GARCH	iGARCH
exponential GARCH	eGARCH	eGARCH
gjr-GARCH	gjrGARCH	gjrGARCH
asymmetric power ARCH	apARCH	apARCH
family GARCH	fGARCH	fGARCH
Component sGARCH	csGARCH
Multiplicative Component sGARCH	mcsGARCH
realized GARCH	realGARCH

さらに以下のようなplotコマンドでモデルの診断図やVaRを描画できる。

plot(out.default, which="all")

f:id:eizoo3010:20141103201530p:plain

加えて、このパッケージはバックテスティング機能が豊富に備わっている。

例えば、以下のRolling Estimationではwindowを1日ずつスライドさせながら

1期先の予測を行い種々の指標を計算してくれる。

また、parallelパッケージを用いた並列計算にも対応しているので処理を分割して計算ができる。

library(parallel)
cl <- makePSOCKcluster(4)
roll <- ugarchroll(spec.default, sp500ret, n.start=1000, refit.every=100,
                   refit.window="moving", solver="hybrid", calculate.VaR=T,
                   VaR.alpha=c(0.01, 0.05), cluster=cl, keep.coef=T)
stopCluster(cl)

# Backtest reportの出力
report(roll, type="VaR", VaR.alpha=0.05, conf.level=0.95)
plot(roll, which="all")

出力結果は以下。

> report(roll, type="VaR", VaR.alpha=0.05, conf.level=0.95)
VaR Backtest Report
===========================================
Model:			sGARCH-norm
Backtest Length:	4523
Data:				

==========================================
alpha:			5%
Expected Exceed:	226.2
Actual VaR Exceed:	251
Actual %:		5.5%

Unconditional Coverage (Kupiec)
Null-Hypothesis:	Correct Exceedances
LR.uc Statistic:	2.78
LR.uc Critical:		3.841
LR.uc p-value:		0.095
Reject Null:		NO

Conditional Coverage (Christoffersen)
Null-Hypothesis:	Correct Exceedances and
	Independence of Failures
LR.cc Statistic:	2.851
LR.cc Critical:		5.991
LR.cc p-value:		0.24
Reject Null:		NO

f:id:eizoo3010:20141103211916p:plain

ここでは過去データのテスト検証を見たが、推定したパラメータを用いた予測もできる。

さらに、1期先の予測データを使って更に先の予測を行うブートストラップ的なボラ計算や

シミュレーションも行うことができる。

欠点としては、入力データの長さが100に満たない場合はugarchfitが機能しないという点である。

パラメータ数的にも最低100サンプルの系列データがないと結果に信頼性が持てないということなのかもしれない。

多変量GARCH: rmgarch

rugarchでも一部のモデルは多変量計算ができるが、より本格的なモデルはrmgarchで実装されている。

扱えるモデルは以下。

DCC-GARCH

GARCH-Copula

GO-GARCH

とりあえずDCC-GARCHを計算してみると以下のようになる。

単変量GARCHのspecを複製してrugarch::multispecに渡してその結果をrmgarch::dccspecに入れている。

library(rmgarch)
data(dji30retw)

spec.multi <- multispec(replicate(5, ugarchspec(mean.model=list(armaOrder=c(1, 0)))))
spec.dcc <- dccspec(spec.multi, dccOrder=c(1, 1), distribution="mvnorm")
fit.dcc <- dccfit(spec.dcc, data=dji30retw[, 1:5], solver="solnp", fit.control=list(eval.se=TRUE))

で結果は以下のような形になって出てくる。

> fit.dcc

*---------------------------------*
*          DCC GARCH Fit          *
*---------------------------------*

Distribution         :  mvnorm
Model                :  DCC(1,1)
No. Parameters       :  37
[VAR GARCH DCC UncQ] : [0+25+2+10]
No. Series           :  5
No. Obs.             :  1141
Log-Likelihood       :  10645.03
Av.Log-Likelihood    :  9.33 

Optimal Parameters
-----------------------------------
              Estimate  Std. Error   t value Pr(>|t|)
[AA].mu       0.002679    0.001255  2.135026 0.032759
[AA].ar1     -0.025964    0.030103 -0.862493 0.388416
[AA].omega    0.000068    0.000054  1.253456 0.210040
[AA].alpha1   0.116360    0.055375  2.101317 0.035613
[AA].beta1    0.867593    0.061163 14.185032 0.000000
[AXP].mu      0.002733    0.001001  2.729992 0.006334
[AXP].ar1    -0.077552    0.035184 -2.204194 0.027511
[AXP].omega   0.000011    0.000010  1.081608 0.279427
[AXP].alpha1  0.064660    0.013358  4.840447 0.000001
[AXP].beta1   0.934340    0.016757 55.756992 0.000000
[BA].mu       0.003206    0.001026  3.126333 0.001770
[BA].ar1     -0.041330    0.030913 -1.336977 0.181230
[BA].omega    0.000010    0.000009  1.081883 0.279305
[BA].alpha1   0.054330    0.014704  3.694920 0.000220
[BA].beta1    0.943693    0.011622 81.198629 0.000000
[BAC].mu      0.002577    0.000972  2.650778 0.008031
[BAC].ar1    -0.029489    0.033589 -0.877936 0.379978
[BAC].omega   0.000010    0.000010  0.948148 0.343054
[BAC].alpha1  0.066366    0.015806  4.198818 0.000027
[BAC].beta1   0.932634    0.020835 44.762431 0.000000
[C].mu        0.002273    0.007420  0.306304 0.759373
[C].ar1      -0.077359    0.035212 -2.196977 0.028022
[C].omega     0.000032    0.000711  0.044824 0.964248
[C].alpha1    0.087233    0.576555  0.151300 0.879739
[C].beta1     0.911767    0.819853  1.112111 0.266091
[Joint]dcca1  0.007705    0.005888  1.308554 0.190685
[Joint]dccb1  0.981505    0.011801 83.173449 0.000000

Information Criteria
---------------------
                    
Akaike       -18.594
Bayes        -18.431
Shibata      -18.596
Hannan-Quinn -18.533


Elapsed time : 8.944512

単変量と同じようにplotでいろんな分析図を見ることができる。

その他、詳しいことは以下の参考文献/webを見てほしい。

参考文献/web

rugarch vignettes

rmgarch vignettes

eizoo3010の日記

社畜の備忘録

■単変量/多変量GARCHモデル

概要

インストール

単変量GARCH: rugarch

多変量GARCH: rmgarch

参考文献/web