R语言用多元ARMA,GARCH ,EWMA, ETS,随机波动率SV模型对金融时间序列数据建模（上）：https://developer.aliyun.com/article/1490539
我们最终可以比较静态误差和滚动误差：
barplot(rbind(error\_var\[, "out-of-sample"\], rolling\_error_var\[, "out-of-sample"\]) col = c("darkblue", "darkgoldenrod"), legend = c("静态预测", "滚动预测"),
我们可以看到，滚动预测在某些情况下是必须的。因此，实际上，我们需要定期进行滚动预测改进。
方差模型
ARCH和GARCH模型
对数收益率残差wt的ARCH（m）模型为
其中zt是具有零均值和恒定方差的白噪声序列，而条件方差σ2t建模为
其中，m为模型阶数，ω> 0，αi≥0为参数。
GARCH（m，s）模型使用σ2t上的递归项扩展了ARCH模型：
其中参数ω\> 0，αi≥0，βj≥0需要满足∑mi =1αi+ ∑sj = 1βj≤1的稳定性。
rugarch生成数据
首先，我们需要定义模型：
# 指定具有给定系数和参数的GARCH模型#> #> *---------------------------------*#> * GARCH Model Spec *#> *---------------------------------*#> #> Conditional Variance Dynamics #> ------------------------------------#> GARCH Model : sGARCH(1,1)#> Variance Targeting : FALSE #> #> Conditional Mean Dynamics#> ------------------------------------#> Mean Model : ARFIMA(1,0,0)#> Include Mean : TRUE #> GARCH-in-Mean : FALSE #> #> Conditional Distribution#> ------------------------------------#> Distribution : norm #> Includes Skew : FALSE #> Includes Shape : FALSE #> Includes Lambda : FALSE#> Level Fixed Include Estimate LB UB#> mu 0.005 1 1 0 NA NA#> ar1 -0.900 1 1 0 NA NA#> ma 0.000 0 0 0 NA NA#> arfima 0.000 0 0 0 NA NA#> archm 0.000 0 0 0 NA NA#> mxreg 0.000 0 0 0 NA NA#> omega 0.001 1 1 0 NA NA#> alpha1 0.300 1 1 0 NA NA#> beta1 0.650 1 1 0 NA NA#> gamma 0.000 0 0 0 NA NA#> eta1 0.000 0 0 0 NA NA#> eta2 0.000 0 0 0 NA NA#> delta 0.000 0 0 0 NA NA#> lambda 0.000 0 0 0 NA NA#> vxreg 0.000 0 0 0 NA NA#> skew 0.000 0 0 0 NA NA#> shape 0.000 0 0 0 NA NA#> ghlambda 0.000 0 0 0 NA NA#> xi 0.000 0 0 0 NA NA#> $mu#> \[1\] 0.005#> #> $ar1#> \[1\] -0.9#> #> $omega#> \[1\] 0.001#> #> $alpha1#> \[1\] 0.3#> #> $beta1#> \[1\] 0.65true_params#> mu ar1 omega alpha1 beta1 #> 0.005 -0.900 0.001 0.300 0.650
然后，我们可以生成收益率时间序列：
# 模拟一条路径hpath(garch\_spec, n.sim = T)#> num \[1:2000, 1\] 0.167 -0.217 # 绘图对数收益{ plot(synth\_log\_returns, main = "GARCH模型的对数收益", lwd = 1.5) lines(synth\_volatility
GARCH
现在，我们可以估计参数：
# 指定一个GARCH模型ugarchspec(mean.model = list(armaOrder = c(1,0)# 估计模型coef(garch_fit)#> mu ar1 omega alpha1 beta1 #> 0.0036510100 -0.8902333595 0.0008811434 0.2810460728 0.6717486402#> mu ar1 omega alpha1 beta1 #> 0.005 -0.900 0.001 0.300 0.650# 系数误差#> mu ar1 omega alpha1 beta1 #> 0.0013489900 0.0097666405 0.0001188566 0.0189539272 0.0217486402
我们还可以研究样本数量T对参数估计误差的影响：
# 循环for (T_ in T\_sweep) { garch\_fit error\_coeffs\_vs\_T <- rbind(error\_coeffs\_vs\_T, abs((coef(garch\_fit) - true\_params)/true\_params)) estim\_coeffs\_vs\_T <- rbind(estim\_coeffs\_vs\_T, coef(garch\_fit))# 绘图matplot(T\_sweep, 100*error\_coeffs\_vs\_T, main = "估计GARCH系数的相对误差", xlab = "T", ylab = "误差 (%)",
真实的ω几乎为零，因此误差非常不稳定。至于其他系数，就像在ARMA情况下一样，μ的估计确实很差（相对误差超过50％），而其他系数似乎在T = 800个样本后得到了很好的估计。
GARCH结果比较
作为健全性检查，我们现在将比较两个软件包 fGarch 和 rugarch的结果：
# 指定具有特定参数值的ARMA（0,0）-GARCH（1,1）作为数据生成过程garch\_spec #生成长度为1000的数据path(garch\_fixed\_spec, n.sim = 1000)@path$# 使用“ rugarch”包指定和拟合模型rugarch\_fit <- ugarchfit(spec = garch\_spec, data = x)# 使用包“ fGarch”拟合模型garchFit(formula = ~ garch(1, 1), data = x, trace = FALSE)# 比较模型系数#> mu omega alpha1 beta1 #> 0.09749904 0.01395109 0.13510445 0.73938595#> mu omega alpha1 beta1 #> 0.09750394 0.01392648 0.13527024 0.73971658# 比较拟合的标准偏差print(head(fGarch\_fi#> \[1\] 0.3513549 0.3254788 0.3037747 0.2869034 0.2735266 0.2708994print(head(rugar#> \[1\] 0.3538569 0.3275037 0.3053974 0.2881853 0.2745264 0.2716555
确实，这两个软件包给出了相同的结果。
使用rugarch包进行GARCH预测
一旦估计出GARCH模型的参数，就可以使用该模型预测未来的值。例如，基于过去的信息对条件方差的单步预测为
给定ω^ /（1-∑mi =1α^ i-∑sj =1β^ j）。软件包 rugarch 使对样本外数据的预测变得简单：
# 估计模型，不包括样本外garch\_fit coef(garch\_fit)#> mu ar1 omega alpha1 beta1 #> 0.0034964331 -0.8996287630 0.0006531088 0.3058756796 0.6815452241# 预测整个样本的对数收益garch\_fore@forecast$sigmaFor\[1, \]# 对数收益图plot(cbind("fitted" = fitted(garch\_fit), main = "合成对数收益预测", legend.loc = "topleft")
#波动率对数收益图plot(cbind("fitted volatility" = sigma(garch_fit), main = "预测合成对数收益率的波动性", legend.loc = "topleft")
不同方法
让我们首先加载S＆P500：
# 加载标准普尔500指数数据head(SP500\_index\_prices)#> SP500#> 2008-01-02 1447.16#> 2008-01-03 1447.16#> 2008-01-04 1411.63#> 2008-01-07 1416.18#> 2008-01-08 1390.19#> 2008-01-09 1409.13# 准备训练和测试数据x\_trn <- x\[1:T\_trn\]x\_tst <- x\[-c(1:T\_trn)\]# 绘图{ plot(x, main = "收益" addEventLines(xts("训练", in
常数
让我们从常数开始：
plot(cbind(sqrt(var\_constant), x\_trn) main = "常数")
移动平均值
现在，让我们使用平方收益的移动平均值：
plot(cbind(sqrt(var\_t), x\_trn), main = "基于简单滚动平方均值的包络线(时间段=20)
EWMA
指数加权移动平均线（EWMA）：
请注意，这也可以建模为ETS（A，N，N）状态空间模型：
plot(cbind(std\_t, x\_trn), main = "基于平方EWMA的包络")
乘法ETS
我们还可以尝试ETS模型的不同变体。例如，具有状态空间模型的乘性噪声版本ETS（M，N，N）：
plot(cbind(std\_t, x\_trn), col = c("red", "black") main = "基于平方的ETS（M，N，N）的包络"
ARCH
现在，我们可以使用更复杂的ARCH建模：
plot(cbind(std\_t, x\_trn), col = c("red", "black") main = "基于ARCH(5)的包络")
GARCH
我们可以将模型提升到GARCH：
plot(cbind(std\_t, x\_trn), col = c("red", "black") main = "基于GARCH（1,1）的包络")
SV随机波动率
最后，我们可以使用随机波动率建模：
或者，等效地，
plot(cbind(std\_t, x\_trn), col = c("red", "black"), main = "基于随机波动率的包络分析")
比较
现在，我们可以比较每种方法在样本外期间的方差估计中的误差：
#> MA EWMA ETS(M,N,N) ARCH(5) GARCH(1,1) SV #> 2.204965e-05 7.226188e-06 3.284057e-06 7.879039e-05 6.496545e-06 6.705059e-06barplot(error_all, main = "样本外方差估计中的误差"
滚动窗口比较
六种方法的滚动窗口比较：MA，EWMA，ETS（MNN），ARCH（5），GARCH（1,1）和SV。
#滚动窗口lookback <- 200len\_tst <- 40for (i in seq(lookback, T-len\_tst, by = len\_tst)) { # MA var\_t <- roll\_meanr(x\_trn^2, n = 20, fill = NA) var\_fore <- var(x\_trn/sqrt(var\_t), na.rm = TRUE) * tail(var\_t, 1) error\_ma <- c(error\_ma, abs(var\_fore - var\_tst)) # EWMA error\_ewma <- c(error\_ewma, abs(var\_fore - var\_tst)) # ETS(M,N,N) error\_ets\_mnn <- c(error\_ets\_mnn, abs(var\_fore - var\_tst)) # ARCH error\_arch <- c(error\_arch, abs(var\_fore - var\_tst)) # GARCH error\_garch <- c(error\_garch, abs(var\_fore - var\_tst)) # SV error\_sv <- c(error\_sv, abs(var\_fore - var\_tst))}barplot(error_all, main = "方差估计误差",
多元GARCH模型
出于说明目的，我们将仅考虑恒定条件相关（CCC）和动态条件相关（DCC）模型，因为它们是最受欢迎的模型。对数收益率残差wt建模为
其中zt是具有零均值和恒定协方差矩阵II的iid白噪声序列。条件协方差矩阵Σt建模为
其中Dt = Diag（σ1，t，...，σN，t）是标准化噪声向量C，协方差矩阵ηt=C-1wt（即，它包含等于1的对角线元素）。
基本上，使用此模型，对角矩阵Dt包含一组单变量GARCH模型，然后矩阵C包含序列之间的一些相关性。该模型的主要缺点是矩阵C是恒定的。为了克服这个问题，DCC被提议为
其中Ct包含等于1的对角元素。要强制等于1的对角元素，Engle将其建模为
Qt具有任意对角线元素并遵循模型
我们将生成数据，估计参数和预测。
从加载多元ETF数据开始：
- SPDR S＆P 500 ETF
- 20年以上国债ETF
- IEF：7-10年期国债ETF
# 下载数据prices <- xts()head(prices)#> SPY TLT IEF#> 2013-01-02 127.8779 99.85183 93.65224#> 2013-01-03 127.5890 98.49886 93.17085#> 2013-01-04 128.1493 98.88306 93.21463#> 2013-01-07 127.7991 98.92480 93.26714#> 2013-01-08 127.4314 99.57622 93.49468#> 2013-01-09 127.7553 99.48438 93.54719# 绘制三个对数价格序列plot(log(prices) main = "三个ETF的对数价格", legend.loc = "topleft")
首先，我们定义模型：
# 指定i.i.d.单变量时间序列模型ugarch_spec # 指定DCC模型spec( multispec(replicate(spec, n = 3))
接下来，我们拟合模型：
# 估计模型#> #> *---------------------------------*#> * DCC GARCH Fit *#> *---------------------------------*#> #> Distribution : mvnorm#> Model : DCC(1,1)#> No. Parameters : 44#> \[VAR GARCH DCC UncQ\] : \[30+9+2+3\]#> No. Series : 3#> No. Obs. : 1007#> Log-Likelihood : 12198.4#> Av.Log-Likelihood : 12.11 #> #> Optimal Parameters#> -----------------------------------#> Estimate Std. Error t value Pr(>|t|)#> \[SPY\].omega 0.000004 0.000000 11.71585 0.000000#> \[SPY\].alpha1 0.050124 0.005307 9.44472 0.000000#> \[SPY\].beta1 0.870051 0.011160 77.96041 0.000000#> \[TLT\].omega 0.000001 0.000001 0.93156 0.351563#> \[TLT\].alpha1 0.019716 0.010126 1.94707 0.051527#> \[TLT\].beta1 0.963760 0.006434 149.79210 0.000000#> \[IEF\].omega 0.000000 0.000001 0.46913 0.638979#> \[IEF\].alpha1 0.031741 0.023152 1.37097 0.170385#> \[IEF\].beta1 0.937777 0.016498 56.84336 0.000000#> \[Joint\]dcca1 0.033573 0.014918 2.25044 0.024421#> \[Joint\]dccb1 0.859787 0.079589 10.80278 0.000000#> #> Information Criteria#> ---------------------#> #> Akaike -24.140#> Bayes -23.925#> Shibata -24.143#> Hannan-Quinn -24.058#> #> #> Elapsed time : 0.8804049
我们可以绘制时变相关性：
# 提取时变协方差和相关矩阵dim(dcc\_cor)#> \[1\] 3 3 1007#绘图plot(corr\_t main = "时变相关", legend.loc = "left")
我们看到两个收益ETF之间的相关性非常高且相当稳定。与SPY的相关性较小，在小于0的区间波动。
