基于ARMA-偏tGARCH和DCC-GARCH模型测算CoVaR——R语言实现

CoVaR是目前金融学界和管理实践中较为主流的测量一个机构（系统）对另一个机构（系统）风险溢出的指标，计算CoVaR的方法主要有分位数回归法、Coupla模型和DCC-GARCH型。本文主要介绍如何利用DCC-GARCH模型对CoVaR进行计算并利用R实现。代码见文末。

CoVaR

CoVaR这一概念由VaR衍生而来，其经济含义是当某一个机构发生风险时，另一机构所承担的总体在险价值。VaR则是指在一定的持有期内，在一定的置信水平下，机构所面临的最大潜在损失。VaR的一般表达式可以为：
$\mathrm{prob}(\Delta P>\mathrm{VaR})=1-T$

根据定义，CoVaR的一般表达式则为：

$\mathrm{Pr}\left({\mathrm{X}}^{\mathrm{i}}⩽{\mathrm{C}\mathrm{o}\mathrm{V}\mathrm{a}\mathrm{R}}_{\mathrm{q}}^{\mathrm{i}\mathrm{j}}\mid {\mathrm{X}}^{\mathrm{j}}={\mathrm{V}\mathrm{a}\mathrm{R}}_{\mathrm{q}}^{\mathrm{j}}\right)=\mathrm{q}$

机构i在机构j发生风险的情况下所承担的总体风险是${\mathrm{C}\mathrm{o}\mathrm{V}\mathrm{a}\mathrm{R}}_{\mathrm{q}}^{\mathrm{i}\mathrm{j}}$，如果要剥离机构i本身的自有风险，计算得到机构j给机构i带来的风险增量，则就采用$\Delta {\mathrm{C}\mathrm{o}\mathrm{V}\mathrm{a}\mathrm{R}}_{\mathrm{q}}^{\mathrm{i}\mathrm{j}}$指标：

$\Delta {\mathrm{C}\mathrm{o}\mathrm{V}\mathrm{a}\mathrm{R}}_{\mathrm{q}}^{\mathrm{i}\mathrm{j}}={\mathrm{C}\mathrm{o}\mathrm{V}\mathrm{a}\mathrm{R}}_{\mathrm{q}}^{\mathrm{i}\mathrm{j}}-{\mathrm{V}\mathrm{a}\mathrm{R}}_{\mathrm{q}}^{\mathrm{i}}$

DCC-GARCH的整体思想

DCC-GARCH模型主要用于描述不同变量间（本文主要讲双变量）的动态变化关系以及相互影响，其核心步骤在于：①先计算单变量的GARCH模型；②然后通过所得相关参数计算双变量间的动态相关系数$\rho$。因此，DCC-GARCH模型的出发点在于构建GARCH模型。

一般的GARCH模型表达式如下（以GARCH（1,1）为例）：
$r_t=u_t+{\epsilon }_t$

${\epsilon }_t={\sigma }_t{z}_t$

${\sigma }_t^2=\omega +\alpha {\epsilon }_{t-1}^2+\beta {\sigma }_{t-1}^2$

$r_t=u_t+{\phi }_1{r}_{t-1}+{\theta }_1{\epsilon }_{t-1}+{\epsilon }_t$

由上式可见，GARCH（1，1）模型主要由收益率均值模型和方差模型两部分组成，对于收益率均值模型而言，可以进一步考虑AR（p）、MA（q）以及ARMA（p，q）模型。在本文的代码中，对收益率采用自动定阶的方法取得ARMA（p，q）。同时，$z_t$这个参数的分布选取也是值得考究的，在金融市场当中，普遍认为金融风险具有“尖峰厚尾”的特征，因此对于$z_t$，可不必拘泥于完全理想的标准正态分布，可将其设定成服从偏t分布。如此一来，对于每一个单变量的收益率，我们就可以构建一个更符合现实情况的ARMA(p,q)-偏tGARCH模型。

如果涉及到多变量的DCC-GARCH，则涉及到矩阵和协方差的概念，如下面的式子：
$H_t=D_t{R}_t{D}_t$

$D_t=\mathrm{diag}\left\{\sqrt{h_t^j}\sqrt{h_t^j}\right\}$

$R_t={\left(Q_t^{\ast }\right)}^{-1}{Q}_t{\left(Q_t^{\ast }\right)}^{-1}$

$Q_t=(1-\alpha -\beta )\bar{Q}+\alpha \left({\epsilon }_{t-1}^i{\epsilon }_{t-1}^j,\right)+\beta Q_{t-1}$

其中，$H_t$是条件协方差矩阵，$D_t$是条件标准差$\sqrt{h_t}$组成的对角矩阵，条件标准差则来源于单变量GARCH拟合模型。$R_t$是动态相关系数矩阵，$Q_t$是协方差矩阵，$\bar{Q}$是经过残差标准化后的无条件协方差矩阵，$Q_t^{\ast }$是$Q_t$的对角矩阵。

DCC-GARCH（1,1）下两机构间的动态相关系数就可以表述为以下形式：
${\rho }_{ij,t}=\frac{(1-\alpha -\beta ){\bar{q}}_{ij}+\beta q_{ij,t-1}+\alpha {\delta }_{i,t-n}{\delta }_{j,t-n}}{{\left[(1-\alpha -\beta ){\bar{q}}_{ij}+\beta q_{ij,t-1}+\alpha {\delta }_{i,t-n}^2\right]}^{1/2}{\left[(1-\alpha -\beta ){\bar{q}}_{ij}+\beta q_{ij,t-1}+\alpha {\delta }_{j,t-1}^2\right]}^{1/2}}$

在DCC-GARCH的框架下，我们已经得到了动态相关系数，那么VaR和CoVaR的表达式就可以分别写成：
$Va{R}_{q,t}^i={\hat{u}}_t^i-Q(q){\hat{h}}_t^i$

${\mathrm{CoVaR}}_{q,t}^{ij}={\gamma }_t^{ij}Va{R}_{q,t}^i$

$\Delta {\mathrm{CoVaR}}_{q,t}^{ij}={\gamma }_t^{ij}\left(Va{R}_{q,t}^i-Va{R}_{50\%,t}^i\right)$

${\gamma }_t^{ij}={\rho }_{ij,t}{\sigma }_{i,t}^2/{\sigma }_{j,t}^2$

其中，Q(q)为置信水平为1-q时机构i收益率所符合的分布下的q分位数值，求解$Va{R}_{q,t}^i$的关键在于通过拟合单变量GARCH模型得到条件标准差${\hat{h}}_t^i$，至于Q(q)，可以利用现有样本取分位数，也可以通过蒙特卡洛的方法构建更大的样本取分位数，这个在本文中不讨论。而求解${\mathrm{CoVaR}}_{q,t}^{ij}$的关键则在于通过DCC函数找到动态相关系数${\rho }_{ij,t}$。

R语言代码实现（以测量我国银行对股市的风险溢出效应为例）

已有很多学术研究都在探讨我国银行系统对于股市的风险溢出效应，本文就以此例作为分析。从Wind上选取2010年1月1日到2020年12月31日申万银行指数以及上证指数共计11年的历史数据。

###导入数据###
##Bank_data:经对数化处理的银行收益率序列
##system_data:经对数化处理的上证指数收益率序列

###ARMA定阶部分###
auto.arima(Bank_data)    ##对银行收益率进行自动定阶处理，结果为ARMA(3,2)
auto.arima(system_data)  ##对上证指数收益率进行自动定阶处理,结果为ARMA(1,0)
arma_Bank=arima(Bank_data,order=c(3,0,2))    
arma_system=arima(system_data,order=c(1,0,0))
arma_Bank_resd=arma_Bank$residuals     #获得残差数据，下同
arma_system_resd=arma_system$residuals                  
resd_series=cbind(arma_Bank_resd,arma_system_resd)   #得到两个序列在ARMA模型下的残差序列
data_matrix=cbind(Bank_data,system_data)

###构建单序列ARMA-GARCH模型部分###
module_Bank <- ugarchspec(variance.model = list(model = "sGARCH", 
                                           garchOrder = c(1, 1), 
                                           submodel = NULL, 
                                           external.regressors = NULL, 
                                           variance.targeting = FALSE), 
                     
                     mean.model     = list(armaOrder = c(3, 2), 
                                           external.regressors = NULL, 
                                           start.pars = list(), 
                                           fixed.pars = list()),
                     distribution.model = "sstd")   #选择偏t分布
fitgarch_Bank=ugarchfit(module_Bank,data=Bank_data,solver="nlminb")        #通过ARMA(3,2)-偏tGARCH（1,1）拟合银行业的收益率指数

module_system <- ugarchspec(variance.model = list(model = "sGARCH", 
                                                  garchOrder = c(1, 1), 
                                                  submodel = NULL, 
                                                  external.regressors = NULL, 
                                                  variance.targeting = FALSE), 
                            
                            mean.model     = list(armaOrder = c(1, 0), 
                                                  external.regressors = NULL, 
                                                  start.pars = list(), 
                                                  fixed.pars = list()),
                            distribution.model = "sstd")  #选择偏t分布
fitgarch_system=ugarchfit(module_system,data=system_data,solver="nlminb")          #通过ARMA(1,0)-偏tGARCH（1,1）拟合上证指数收益率指数

###构建双序列DccGARCH模型###
dcc.garch.spec = dccspec(uspec = multispec(c(module_Bank, 
                                             module_system)),
                         dccOrder = c(1,1),
                         distribution = "mvt")
dcc.fit = dccfit(dcc.garch.spec, data = data_matrix[,c(1,2)])   #构建双序列DccGARCH模型
dcc.fit    #显示模型拟合结果
cor<-data.frame(dcc.fit@mfit$R)
cor<-t(cor)
cor<-cor[,1]
del<-seq(1,length(cor),by=2)
cor<-cor[-del]
corr_DccGARCH=as.numeric(cor)        #得到传统DccGARCH模型下的动态相关系数
class(time_data)
corr_DccGARCH_plot<-zoo(corr_DccGARCH,as.Date(as.character(time_data)))  #生成传统DccGARCH模型下的动态相关系数图

###bayesDccGARCH部分###
out_data=bayesDccGarch(resd_series,nSim=5000,tail_ini=8,alpha_ini=rep(0.01,ncol(resd_series)),beta_ini=rep(0.97,ncol(resd_series)),errorDist=2)     #errorDist=2表示残差服从偏态t分布，该贝叶斯分布通过MH算法得到
bayes_var_x=out_data$H[,1]    
bayes_var_y=out_data$H[,4]     
bayes_cvar_xy=out_data$H[,3]
bayes_corr=sqrt((bayes_cvar_xy)/(bayes_var_x*bayes_var_y))          #贝叶斯DCCGARCH模型下的动态相关系数
bayes_corr_DccGARCH_plot<-zoo(bayes_corr,as.Date(as.character(time_data)))

###动态相关性作图对比###
plot(bayes_corr_DccGARCH_plot,ylab="corr",main="银行-金融市场")
lines(corr_DccGARCH_plot,col="blue")                     #动态相关性对比

###Dcc-GARCH方法下的CoVaR计算###
sigma_DccGARCH_data=dcc.fit@model[["sigma"]]
mu_DccGARCH_data=dcc.fit@model[["mu"]]
i_sigma=sigma_DccGARCH_data[,1]
system_sigma=sigma_DccGARCH_data[,2]
i_mu=mu_DccGARCH_data[,1]
VaR_i_0.05=i_mu+(quantile(ecdf(Bank_data),0.05))*(i_sigma)
VaR_i_0.5=i_mu+(quantile(ecdf(Bank_data),0.5))*(i_sigma)
gamma_DccGARCH=(corr_DccGARCH*((system_sigma)^2)/((i_sigma)^2))
CoVaR_system_DccGARCH=gamma_DccGARCH*VaR_i_0.05      #得到DCC-GARCH方法下的CoVaR值)