1. 資産組合VaR模組化方法回顧 文章
中總結了通過DCC模型估計組合向前一日VaR的方法，整體思路如下：
● 通過Garch族模型估計各資産的波動率
● 通過DCC模型估計各資産間的相關系數，結合1得到資産組合的協方差矩陣
● 在各資産正态性假設的前提下，可以知道資産組合也服從正态分布，并且均值與協方差陣已在1,2中計算得到
● 在已知組合中各但資産權重w的情況下，根據下式計算組合VaR
中總結了通過蒙特卡洛方法估計組合向前K日VaR的方法，也可以僅計算組合向前一日VaR（本文隻考慮向前1日的情況），文章中也對比了蒙特卡洛方法與DCC方法得到的結果，差異并不大。蒙特卡洛方法的思路如下：
● 根據Garch族模型估計資産的波動率
● 根據DCC模型估計組合的相關系數
● 在1,2的基礎上，在正态性假設前提下，得到組合的分布函數，對組合收益率進行模拟，在給定各資産權重w的情況下，可以得到組合的總收益
● 重複1-3若幹次，可以得到組合總收益的模拟序列，類似HS方法，取p分位數即可
可以看出不論是DCC模型還是蒙特卡洛方法，都是在正态性假設的前提下，得到組合的分布函數再進行求解。事實上，也可以類比多元正态的概念建構多元t分布和多元漸進t分布，假設組合服從這樣的分布，求出分布的參數後，再用蒙特卡洛方法進行模拟，這些理論依據已經很成熟，推導過程見文獻[1]，這裡不再贅述。
但需要說明的是，多元t分布和多元漸近t分布都沒有邊際分布和線性組合依然多元t或者多元漸近t的性質。回憶多元正态的情況下，為了生成多元正态随機數，實際上是先産生不相關的n組一進制正态随機數向量，然後通過cholesky分解轉換為符合給定相關系數矩陣的組合收益率模拟序列。如果組合的分布不具有類似多元正态的性質，要根據分布函數模拟組合收益就比較困難，必須直接通過多元分布函數産生随機數，不能分解成單個資産去做，雖然也有相關的方法可以生成給定分布函數下随機數，但都比較麻煩，這是之前方法的一個局限性。
此外，多元正态假設所有的單個資産都是正态分布，多元t分布和多元漸近t分布的邊際分布并非t分布或者漸近t分布，而不同的資産可能服從不同的分布，需要用不同方法去模組化，已有的多元分布都不能滿足這一條件，這是之前方法的另一局限性。
比較理想的狀态是，我們可以用不同的方法對不同的單資産進行模組化，最終n各資産具有不同的分布函數
這也正是本文總結的Copula模型的邏輯。
2.Copula模型
Sklar定理
Copula模型整體來說比較複雜，這裡隻對關鍵的部分加以說明，模型中最重要的定理是Sklar定理，也就是上面所說的理想情況，具體叙述如下
G稱為copula CDF，在sklar定義的假設下，如果我們已經通過一些單變量模型得到了單資産的分布函數，隻需要确定出copula函數G，就相當于知道了組合的分布函數，進而把估計組合分布函數的問題轉化為估計copula函數的問題。當然copula函數也不是靠猜，有一些常用的copula函數可以選擇，在确定了copula函數之後，可以通過MLE等方法估計參數。
參數估計（MLE）
其中，序号1稱為Gumbel Copula函數，序号2稱為Clayton Copula函數，序号3稱為Frank Copula函數，之是以說明這三個，是因為這三個實際應用中比較多，python的copulalib包中也隻提供這三種方法，不過本文并未嘗試這幾種方法，有興趣的可以自己嘗試下。
VaR估計思路
從之前的叙述中可以看出，通過copula函數得到的組合分布函數沒有非常好的解析表達式，是以直接通過定義計算VaR的方法行不通，一般采取與蒙特卡洛方法相結合的方式，生成給定copula函數下的随機數，模拟資産組合的收益序列，再根據組合權重得到組合總收益，重複若幹次，取p分位數。
随機數構造
使用蒙特卡洛方法的難點在于生成給定copula函數下的随機數，需要用到Nelsen定理，詳見參考文獻[2]
服從二進制正态，可以直接模拟，然後再用标準正态分布函數作用，就可以得到符合給定多元正态copula的随機數，多元t-copula分布類似。
在得到符合給定copula分布的随機數u後，根據單個資産的分布F，可以得到單資産對應的随機數z
随後可以根據權重計算組合收益進而估計VaR。
綜上，可以将Copula函數估計VaR的過程總結如下
選擇copula函數，估計參數
第一步：根據單變量模型對所有單資産進行模組化，估計分布函數F；
第二步：根據所有的分布函數F和給定copula函數，最大化對數似然函數估計參數；
蒙特卡洛模拟估計VaR
第一步：生成符合copula函數的随機數;
第二步：通過随機數得到各資産收益的模拟序列；
第三步：根據各資産權重得到組合收益序列，取p分位數作為VaR估計值
3.實證分析
資料：S&P500、US 10yr T-Note Fixed Term（同上一篇）
區間：2001-2010
蒙特卡洛模拟次數：10000次
資料和代碼在背景回複“VaR5”擷取
僅估計最後一天的VaR。代碼中未給出太多注釋，可以參見文獻[1]第九章習題。
前兩道題首先通過threshold correlation說明正态性假設并不符合實際，threshold correlation定義如下，r(p)表示r的p分位數
1def getThre_cor(data,column1,column2,p):
2 cor = pd.DataFrame(p,columns = ['p'])
3 cor['thre_cor'] = 0
4 for i in range(cor.shape[0]):
5 if p[i] <=0.5:
6 condition1 = (data[column1] <= np.percentile(data[column1],p[i]*100))
7 condition2 = (data[column2] <= np.percentile(data[column2],p[i]*100))
8 else:
9 condition1 = (data[column1] > np.percentile(data[column1],p[i]*100))
10 condition2 = (data[column2] > np.percentile(data[column2],p[i]*100))
11 datas = data.loc[condition1 & condition2,:]
12 cor.loc[i,'thre_cor'] = np.corrcoef(datas[column1],datas[column2])[0,1]
13 return cor
14
15def Thre_cor_norm(num,rou):
16 np.random.seed(52)
17 data = pd.DataFrame(index = range(num))
18 data['r1'] = np.random.normal(size=(num,1))
19 data['r2'] = np.random.normal(size=(num,1))
20
21 data['r1_c'] = data['r1']
22 data['r2_c'] = data.r1*rou + data.r2*(1 - rou**2)**0.5
23 return data
24
25p = np.arange(0.15,0.90,0.05)
26rou = np.corrcoef(data1.Log_Return_SP,data1.Log_Return_US)[0,1]
27data_norm = Thre_cor_norm(30000,rou)
28cor_norm = getThre_cor(data_norm,'r1_c','r2_c',p)
29
30p = np.arange(0.15,0.86,0.01)
31cor = getThre_cor(data1,'Log_Return_SP','Log_Return_US',p)
32ax = plt.figure(figsize=(10,5))
33plt.plot(cor.p,cor.thre_cor,linewidth = 2)
34plt.plot(cor_norm.p,cor_norm.thre_cor,linewidth = 2,color = 'red')
35plt.grid()
36plt.show()
第三道題為用t-garch分别對兩個單資産進行模組化，估計參數d，不再說明；
第四道題為用第三問的結果建立二進制正态copula模型，估計組合VaR,過程前面已經說明，代碼如下
估計copula函數的參數
1def getNegativeLoglikelihood_copula(rou,r):
2 LogLikeLihood = -r.shape[0]*np.log(1 - rou**2)/2 - ((r.norm1**2 +r.norm2**2 - 2*rou*r.norm1*r.norm2)/(2*(1 - rou**2)) -
3 0.5*(r.norm1**2 + r.norm2**2)).sum()
4 return -LogLikeLihood
5
6rou_best = optimize.fmin(getNegativeLoglikelihood_copula,rou0, \
7 args=(copula_data,),ftol = 0.000000001)
8print('估計結果為：',rou_best)
模拟
1data4['u1c'] = data4['u1']
2data4['u2c'] = data4.u1*rou + data4.u2*(1 - rou**2)**0.5
3data4['F1'] = norm(0,1).cdf(data4['u1c'])
4data4['F2'] = norm(0,1).cdf(data4['u2c'])
5data4['z1'] = t(d_SP).ppf(data4.F1)*((d_SP-2)/d_SP)**0.5
6data4['z2'] = t(d_US).ppf(data4.F2)*((d_US-2)/d_US)**0.5
7data4['R1'] = data4.z1*sigma_SP**0.5
8data4['R2'] = data4.z2*sigma_US**0.5
9data4['R'] = 0.5*data4.R1 +0.5* data4.R2
10VaR = -np.percentile(data4.R,1)
最終估計結果為VaR = 0.0101，可以與上篇文章最後一日的結果相對比，基本上是一緻的。
原文釋出時間為：2018-10-1
本文作者：量化小白H
本文來自雲栖社群合作夥伴“
Python愛好者社群”，了解相關資訊可以關注“
”。