VaR系列（三）：DCC模型估计组合VaR
1.模型推导
和单个资产类似，资产组合的VaR定义依然由下式给出
不同的地方在于，这里的波动率应换成组合的波动率，分布函数应换为组合的分布函数。
需要说明的一点是，如果我们假设所有的单个资产收益率都服从正态分布，资产组合的收益率是单个资产收益率的加权和，也服从正态分布，这种情况下，计算VaR只需要对组合的波动率给出估计。如果单个资产收益率不服从正态分布，比如t分布，渐进t分布等等，组合的分布并不是非常容易刻画，因此本篇的实证都假设单个资产收益率均服从正态分布，对于其他分布的建模需要用到copula等模型，在之后几篇中会总结。
假设组合中各资产的权重为w，收益率为r，则t+1时刻的组合收益率为
从而组合的波动率可以表示为
因此，对于组合波动率的估计可以转化为估计资产组合的协方差阵的估计，协方差阵的对角线——单个资产的波动率可以根据之前的模型进行估计，非对角线元素只需要估计资产间的相关系数。
对于相关系数的估计总体来讲有两类模型，一种是常相关系数模型，即认为资产间的相关系数不随时间变化，可以用一个常数来衡量，这种方法不过多赘述，比如可以用历史数据的相关系数作为估计量，另一种是时变相关系数模型，认为相关系数是随时间变化的，这方面的模型很多，这里只总结通过DCC模型对相关系数进行建模的方法。
题外话，如果单个资产数目过多，估计两两资产的相关系数耗时过多，这时可以考虑通过多因子模型，把对于组合协方差阵的估计，转化为对于因子协方差阵的估计，具体推导见文章，因子协方差阵的估计也可以用这里的方法。
Engle(2002)年提出了DCC模型，通过DCC模型对相关系数建的整体思路与之前对于单个资产波动率建模的思路基本一致，包括两种方法：
第一种类似RM方法， 通过指数平滑的方法对对相关系数建模，后文统称为DCC-RM模型，公式来自Engle的paper，需要文献在后台回复"VaR3"，公式中的
即为前文中的z，也就是标准化收益。
另一种方法类似Garch(1,1)模型，后文统称为DCC-Garch模型。
其中
是两两资产间的相关系数，也是q的标准化，通过对q标准化保证相关系数绝对值小于等于1.。
这两种方法最大的区别在于，RM方法不满足均值回归现象（右边两个系数和为1），Garch方法在（alpha+beta<1）的情况下满足均值回归现象，与实际相符。
矩阵形式可以表达为
2.参数估计
正态性假设下，上述两个模型均可以通过极大似然估计法估计参数，之后推导均以两资产情况为例，更高维的情况见文献。在仅有两个资产的情况下，对数似然函数可以表示为
其中，
基于上式，分别讨论两种模型的参数估计过程
- DCC-RM
两个资产下，模型可以表示为
带入对数似然函数，求使对数似然函数最大的lambda的值。
初始值：
关于初始值为什么这样设定，我确实也不太清楚，直接从文献中得到，这里可以讨论，代码见后文getNegLoglikelihood5。
- DCCC-Garch
两资产情况下
实际优化时，初始值设定同上，代码见后文getNegLoglikelihood6。
3.实证分析
软件：python3.6
数据：S&P500、US 10yr T-Note Fixed Term
区间：2001-2010
两种资产及资产组合的净值曲线，组合中两种资产等权重，其中，红色线为组合净值，黑色线为SP500，灰色线为US 10yr T-Note Fixed Term
通过Ngarch（1,1）模型分别对单个资产的波动率进行估计。模型代码见后文。
- S&P500波动率
- US 10yr T-Note Fixed Term波动率
黑色线为资产走势，红色线为波动率。
- 基于DCC-RM模型的时变相关系数
其中，黑色线与灰色线分别为两种资产的净值，红色线为时变相关系数。可以看出，资产同涨同增时，相关系数较大，资产净值走势背离时，条件相关系数较低，与实际比较符合。
- 基于DCC-RM模型的VaR
- 基于DCC-Garch模型的时变相关系数
其中，红色线为DCC-RM估计得到的相关系数，绿色线为DCC-Garch估计得到的相关系数，整体趋势一致。
- 基于DCC-Garch模型的VaR
其中，红色线为DCC-RM估计得到的VaR，绿色线为DCC-Garch估计得到的VaR，整体趋势一致。
4.代码
只包含主要代码，作图代码略，全部代码后台获取。
- 数据准备
1import os
2import pandas as pd
3import numpy as np
4from scipy import optimize
5import matplotlib.pyplot as plt
6from scipy.stats import norm
7
8data = pd.read_excel('VaR3.xlsx')
9data['SP_return'] = np.log(1+data.SP_Close.pct_change(1))
10data['US_return'] = np.log(1+data.US_Close.pct_change(1))
11data = data.dropna()
12data = data.reset_index(drop = True)
13data['nav'] = (data.SP_Close/data.SP_Close[0] + data.US_Close/data.US_Close[0])/2
14data.head()- 波动率Ngarch（1,1）建模，python中的optimze函数只能求最小值，因此对对数似然函数取负数再优化。
1def getNegativeLoglikelihood3(params,r):
2 omega,alpha,beta,theta = params
3 sigma2 = np.ones(len(r))
4 sigma2[0] = np.var(r)
5 for i in range(1,len(r)):
6 s = omega + alpha * (r[i - 1] - theta*sigma2[i-1]**0.5)**2 + beta * sigma2[i - 1]
7 sigma2[i] = (s>0)*s+(s<0)*100
8 LogLikeLihood = (np.log(2*np.pi) + np.log(sigma2) + r**2/sigma2).sum()/2
9 return LogLikeLihood
10
11# SP500
12params_MLE2 = optimize.fmin(getNegativeLoglikelihood3,np.array([0.0000015,0.05,0.8,1.25]), \
13 args=(data['SP_return'],), ftol = 0.00001)
14params_MLE2
15
16omega,alpha,beta,theta = params_MLE2
17data['SP_sigma2'] = np.var(data['SP_return'])
18for i in range(1,data.shape[0]):
19 data.loc[i,'SP_sigma2'] = omega + alpha * (data.loc[i - 1,'SP_return'] - theta*data.loc[i-1,'SP_sigma2']**0.5)**2 + beta * data.loc[i - 1,'SP_sigma2']
20
21# US
22params_MLE2_US = optimize.fmin(getNegativeLoglikelihood3,np.array([0.000005,0.03,0.97,0]), \
23 args=(data['US_return'],), ftol = 0.00001)
24params_MLE2_US
25
26omega,alpha,beta,theta = params_MLE2_US
27data['US_sigma2'] = np.var(data['US_return'])
28for i in range(1,data.shape[0]):
29 data.loc[i,'US_sigma2'] = omega + alpha * (data.loc[i - 1,'US_return'] - theta*data.loc[i-1,'US_sigma2']**0.5)**2 + beta * data.loc[i - 1,'SP_sigma2']
30
31data['SP_z'] = data['SP_return']/data['SP_sigma2']**0.5
32data['US_z'] = data['US_return']/data['US_sigma2']**0.5- DCC-RM
1def getNetLoglikelihood5(param,z1,z2):
2 lambd = param
3 if lambd >1:
4 return (10000)
5 else:
6 q11= np.ones(len(z1))
7
8 q12= np.ones(len(z1))
9 q12[0]= np.mean(z1*z2 )
10
11 q22= np.ones(len(z1))
12
13 rho= np.ones(len(z1))
14 rho[0]= q12[0]/np.sqrt(q11[0]*q22[0])
15
16 for i in range(1,len(z1)):
17 q11[i]= (1-lambd)*z1[i-1]**2 + lambd*q11[i-1]
18 q12[i]= (1-lambd)*z1[i-1]*z2[i-1] + lambd*q12[i-1]
19 q22[i]= (1-lambd)*z2[i-1]**2 + lambd*q22[i-1]
20
21 rho[i]= q12[i]/np.sqrt(q11[i]*q22[i])
22
23 LogLikeLihood = ( np.log(1-rho**2) + (z1**2+z2**2-2*rho*z1*z2)/(1-rho**2)).sum()*0.5
24 return LogLikeLihood
25
26# 参数估计
27params_MLE5 = optimize.fmin(getNetLoglikelihood5,0.94, args=(data['SP_z'],data['US_z'],), ftol = 0.000000001)
28params_MLE5
29
30data['q11'] = 1
31data['q12'] = np.mean(data.SP_z*data.US_z)
32data['q22'] = 1
33
34# 计算COR
35for i in range(1,data.shape[0]):
36 data.loc[i,'q11'] = (1-params_MLE5)*data.SP_z[i-1]**2 + params_MLE5*data.loc[i-1,'q11']
37 data.loc[i,'q12'] = (1-params_MLE5)*data.SP_z[i-1]*data.US_z[i-1] + params_MLE5*data.loc[i-1,'q12']
38 data.loc[i,'q22'] = (1-params_MLE5)*data.US_z[i-1]**2 + params_MLE5*data.loc[i-1,'q22']
39
40data['cor_dcc'] = data.q12/((data.q11*data.q22)**0.5)
41
42# 计算VaR
43data['VaR_DCC'] = -norm(0,1).ppf(0.01)*(data['SP_sigma2']*0.5**2 + data['US_sigma2']*0.5**2 + \
44 2*0.5*0.5*data['cor_dcc']*(data['SP_sigma2']**0.5)*(data['US_sigma2']**0.5))**0.5- DCC-Garch
1def getNegLoglikelihood6(params,z1,z2):
2
3 alpha,beta = params
4
5 if alpha + beta >1:
6 return (10000)
7 else:
8 q11= np.ones(len(z1))
9
10 q12= np.ones(len(z1))
11 q12[0]= np.mean(z1*z2 )
12
13 q22= np.ones(len(z1))
14
15 rho= np.ones(len(z1))
16 rho[0]= q12[0]/np.sqrt(q11[0]*q22[0])
17
18 rho_mean = np.mean(z1*z2)
19 for i in range(1,len(z1)):
20 q11[i]= 1 + alpha*((z1[i-1]**2)-1) + beta*(q11[i-1]-1)
21 q12[i]= rho_mean+ alpha*(z1[i-1]*z2[i-1]- rho_mean)+beta*(q12[i-1]-rho_mean)
22 q22[i]= 1+alpha*((z2[i-1]**2)-1)+beta*(q22[i-1]-1)
23
24 rho[i]= q12[i]/np.sqrt(q11[i]*q22[i])
25 LogLikeLihood = (np.log(1-rho**2) + (z1**2+z2**2-2*rho*z1*z2)/(1-rho**2)).sum()*0.5
26 return LogLikeLihood
27
28# 参数估计
29params_MLE6 = optimize.fmin(getNegLoglikelihood6,np.array([0.05,0.9]), args=(data['SP_z'],data['US_z']), ftol = 0.000000001)
30params_MLE6
31
32# cor
33data['q11_garch'] = 1
34data['q12_garch'] = np.mean(data.SP_z*data.US_z)
35data['q22_garch'] = 1
36rho_mean = np.mean(data.SP_z*data.US_z)
37
38for i in range(1,data.shape[0]):
39 data.loc[i,'q11_garch'] = 1 + alpha*((data.SP_z[i-1]**2)-1) + params_MLE6[1]*(data.loc[i-1,'q11_garch']-1)
40 data.loc[i,'q12_garch'] = rho_mean+ params_MLE6[0]*(data.SP_z[i-1]*data.US_z[i-1]- rho_mean)+ \
41 params_MLE6[1]*(data.loc[i-1,'q12_garch']-rho_mean)
42 data.loc[i,'q22_garch'] = 1+params_MLE6[0]*((data.US_z[i-1]**2)-1)+params_MLE6[1]*(data.loc[i-1,'q22_garch']-1)
43
44data['cor_dcc_garch'] = data.q12_garch/((data.q11_garch*data.q22_garch)**0.5)
45
46# var
47data['VaR_DCC_Garch'] = -norm(0,1).ppf(0.01)*(data['SP_sigma2']*0.5**2 + data['US_sigma2']*0.5**2 + \
48 2*0.5*0.5*data['cor_dcc_garch']*(data['SP_sigma2']**0.5)*(
data['US_sigma2']**0.5))**0.5参考文献
- Christoffersen, Peter F. Elements of financial risk management. Academic Press, 2011.
- Engle, R. F. "Dynamic conditional correlation-a simple class of multivariate GARCH models." Ssrn Electronic Journal 20.3(2000):339--350.