作者：量化小白H Python愛好者社群專欄作者
個人公衆号：量化小白上分記
前文傳送門：
之前幾篇總結的方法都是對于向前一日VaR的模組化，都以是以VaR=波動率乘以分布函數逆函數為基礎。但如果要計算向前k日的VaR，如果還使用上述公式，波動率和分布函數應該換成k日滾動視窗的，好像還沒見過這樣的Garch模型。
是以本文總結一種可以對向前k日VaR進行估計的方法，蒙特卡洛方法。需要說明的是，通過蒙特卡洛方法估計VaR可以從很多角度入手，本文隻是其中一個角度，并非唯一模組化方法。
01
模型推導
Monte-Garch
前面提到，之前的基礎公式不太适用K日VaR的估計，需要尋找其他方法。蒙特卡洛方法從VaR定義出發，要計算資産在給定機率下，未來K日的最大損失，我們可以模拟出資産未來k日的日收益率，計算得到資産的K日收益模拟序列，多次模拟後，類似之前HS和WHS方法，取整個模拟序列的p分位數即可，具體操作如下
我們依然假設資産的日收益率序列服從正态分布，進而标準化收益率服從标準正态分布。
同時收益率序列的波動率可以用Garch族模型進行模組化，假設波動率服從Garch(1,1)模型
這樣，我們可以用标準正态随機數對向前一天的标準化收益率z進行模拟，記模拟次數為MC。
進而可以得到模拟的收益率為
同理對向前兩天的标準化收益率z進行模拟MC次
波動率通過Garch模型進行更新
之後可以得到向前兩天的日收益率
如果要估計向前K天的VaR，我們将上述過程重複K次，過程可以表述如下
這樣向前K日的總收益率可以把k日收益率加總得到（如果是算單利的話，複利用乘法）
這樣我們就得到了對于未來K日總收益的MC個模拟 值，我們用這MC個模拟值估計未來K日的VaR，類似HS方法，直接取p分位數即可
Monte-FHS
上面的方法是将蒙特卡洛模型與Garch組模型相結合，此外，也可以将蒙特卡洛方法與曆史模拟方法相結合，即不對未來收益率進行模組化，而是直接用過去一段時間的收益率作為未來每日收益率的估計量，對上述方法做一點小的修改，具體如下
對于過去一段時間的日收益率，我們可以通過Garch模型對波動率進行估計，得到過去一段時間的标準化收益率
随後，我們每次都從這個收益率序列中進行有放回抽樣，生成對未來k日的收益率模拟序列，再重複波動率的估計和收益率的加總過程
最後求出VaR
這裡做一點分析，這兩種方法可以認為屬于兩類方法：一類可以認為是完全沒有模型的方法，包括WHS,HS，FHS，對于收益率的分布不做任何估計，遵循曆史會重演的原則，用曆史收益率作為未來收益率的估計量，一類是有嚴密假設前提和模型推導的參數方法，包括Garch，EVT等等。
兩種方法各有優劣，從之前文章的結果來看，沒有模型的方法基本上是不能看的，不過從本文之後的兩種方法的對比來看還可以。而有假設有模型的方法結果可能會比較好，但前提是能對收益率序列作出合理假設，從正态分布，t分布，漸近t分布等等，理論研究中嘗試了各種分布，但顯然沒有一定最有效的分布，一切都在于嘗試。
Monte-DCC-Garch
剛才兩種方法都是對單個資産的VaR估計，也可以把蒙特卡洛方法與前一篇文章中的DCC方法相結合，估計組合的向前k日VaR。用Monte-Garch或者Monte-FHS都可以，過程差不多，這裡以Monte-Garch為例。
組合VaR估計與單資産VaR估計的不同在于組合不僅需要估計資産的波動率，還需要估計資産之間的相關性，換句話說，需要估計資産的協方差陣。和之前類似，我們可以模拟組合的日标準化收益率序列，用Garch模型更新波動率，用DCC模型跟新相關系數，然後計算組合的日收益率，K日總收益，最後計算組合的向前K日VaR。
但如果我們直接用之前的方法分别模拟每個但資産的收益率，顯然各個資産之間是不相關的，是以我們必須得到符合給定相關系數矩陣的模拟收益率序列，以2個資産為例進行說明
根據之前的方法我們可以得到如下的兩個不相關的标準化收益率序列
我們的目标是得到符合給定相關系數rho的标準化收益率序列
顯然通過如下轉換即可完成
兩資産情況下，相關系數矩陣的平方根很容易算出來
多資産情況下，相關系數矩陣的平方根不是很好算，但也可以算出來，用到cholesky分解，不過如果是寫代碼的話，也就是一行代碼調個函數的事情。
綜上，計算向前K日VaR整個過程可以表述如下：
1. 模拟生成不相關的資産标準化收益率序列，均服從标準正态分布
2. 将不相關的資産标準化收益率序列轉換為相關系數矩陣為給定值的标準化收益率序列
3. 将标準化收益率轉化為收益率，按照各資産的權重進行累加，得到組合的收益率
4. 通過Garch模型更新波動率，通過DCC模型跟新相關系數矩陣
5. 重複1-4過程K次，将K次的收益率序列相加總收益率
6. 重複1-5過程MC次，得到MC個總收益率，取p分位數即可。
02
實證分析
軟體：python3.6
資料：S&P500、US 10yr T-Note Fixed Term（同上一篇）
區間：2001-2010
蒙特卡洛模拟次數：10000次
資料和代碼在背景回複“VaR4”擷取
單資産的VaR估計均在資料序列最後一天向前估計，估計K=500天
波動率初始值使用當天的波動率，波動率模組化用之前文章中的Ngarch（1，1）模型，兩資産組合的VaR估計類似
Monte-Garch（代碼見VaR_Monte_NGarch）
上一篇中提到，Garch模型具有均值回歸的特性，最終會回到長期平均水準，如果我們将波動率初始值分别改為高初始波動率：3倍長期平均波動率和低初始波動率：0.5倍長期平均波動，結果如下
高初始值的結果先升高再降低，低初始值的結果持續升高，這個結果看起來不是特别明顯，文獻中的文獻中的結果要比我的好很多
Monte-FHS（代碼VaR_Monte_FHS）
對比Monte-Garch和Monte-FHS的結果，整體趨勢差不多，不過Monte-FHS得到的值要相對小一些。不同初始波動率的結果對比如下
此外，代碼中同時還對比了用RM，Monte-Garch、Monte-FHS計算出的向前一日VaR，其中
RM：0.045864978982555343；
Monte-Garch：0.052665510343511503；
Monte-FHS：0.053369497694956955
組合VaR估計（代碼見getVaR_Monte）
這裡與上一篇中DCC模型估計的向前一日VaR結果進行對比
黑色為Monte-Garch-DCC的結果，紅色為DCC的結果，幾乎是完全重合的，黑色被覆寫。但效率上講，蒙特卡洛的運作速度較慢，遠不如DCC。
03
代碼
這裡隻給出主要函數代碼，全部代碼背景擷取，另外代碼實作的過程實際是參考文獻中Chapter8中的習題，注釋寫的不多，看不明白的可以翻一翻文獻，如果想把代碼用到别處，需要做一些修改，比如組合收益率中，因為這裡考慮的是等權重組合，是以并沒有設定權重參數，直接取的平均，如果要考慮非等權，函數輸入參數中需加上權重參數。Monte-Garch1def VaR_Monte_NGarch(num,p,sigma_first,ndays,omega,alpha,beta,theta):2 np.random.seed(4)3 MC = num4 sigma2_Garch = sigma_first5 data_MC = pd.DataFrame(index = range(MC))6 data_R = pd.DataFrame(index = range(MC))7 data_MC['z_day1'] = np.random.normal(size = (10000,1))8 data_R['R_day1'] = sigma2_Garch**0.5 * data_MC['z_day1']9 data_MC['sigma2_day1'] = omega + alpha*(data_R['R_day1'] - theta*sigma2_Garch**0.5)**2 + beta*sigma2_Garch1011 # 低2-10天公式一樣，循環12 for i in range(2,ndays +1):13 exec("data_MC['z_day" + str(i) + "'] = np.random.normal(size = (10000,1))")14 exec("data_R['R_day" + str(i) + "'] = data_MC['sigma2_day" + str(i-1) + "']**0.5 * data_MC['z_day" + str(i) + "']")15 exec("data_MC['sigma2_day" + str(i) + "'] = omega + alpha*(data_R['R_day" + str(i) + "'] - theta*data_MC['sigma2_day" + \ str(i -1) + "']**0.5)**2 + beta*data_MC['sigma2_day" + str(i-1) + "']")16 VaR = pd.DataFrame(index = range(ndays))17 VaR['ndays'] = np.arange(1,ndays+1)18 VaR['VaR'] = 019 for i in range(ndays):20 R_ndays = data_R.iloc[:,:i+1].sum(axis = 1)21 VaR.loc[i,'VaR'] = - np.percentile(R_ndays,p) 22 return VaRMonte-FHS1def VaR_Monte_FHS(num,p,sigma_first,ndays,omega,alpha,beta,theta):2 np.random.seed(4)3 MC = num4 sigma2_FHS = sigma_first5 data_FHS = pd.DataFrame(index = range(MC))6 data_FHS_R = pd.DataFrame(index = range(MC))7 data_FHS['z_day1'] = data.loc[np.random.randint(0,data.shape[0]-1,MC),'z_Garch'].values8 data_FHS_R['R_day1'] = sigma2_Garch**0.5 * data_FHS['z_day1']9 data_FHS['sigma2_day1'] = omega + alpha*(data_FHS_R['R_day1'] - theta*sigma2_FHS**0.5)**2 + beta*sigma2_FHS1011 # 低2-10天公式一樣，循環12 for i in range(2,ndays +1):13 exec("data_FHS['z_day" + str(i) + "'] = data.loc[np.random.randint(0,data.shape[0]-1,MC),'z_Garch'].values")14 exec("data_FHS_R['R_day" + str(i) + "'] = data_FHS['sigma2_day" + str(i-1) + "']**0.5 * data_FHS['z_day" + str(i) + "']")15 exec("data_FHS['sigma2_day" + str(i) + "'] = omega + alpha*(data_FHS_R['R_day" + str(i) + "'] - theta*data_FHS['sigma2_day" + \ str(i -1) + "']**0.5)**2 + beta*data_FHS['sigma2_day" + str(i-1) + "']")1617 VaR = pd.DataFrame(index = range(ndays))18 VaR['ndays'] = np.arange(1,ndays+1)19 VaR['VaR'] = 020 for i in range(ndays):21 R_ndays = data_FHS_R.iloc[:,:i+1].sum(axis = 1)22 VaR.loc[i,'VaR'] = - np.percentile(R_ndays,p) 23 return VaRMonte-DCC-Garch1def getVaR_Monte(n,alpha_DCC,beta_DCC,sigma_first,rho_first,alpha_garch,beta_garch,omega_garch,theta_garch,p):2 MC = 100003 sigma2_Garch = np.dot(np.diag(sigma_first),np.ones([2,MC]))45 np.random.seed(4) 67 VaR = pd.DataFrame(index = range(n))8 VaR['ndays'] = np.arange(1,n+1)9 VaR['VaR'] = 0 10 VaR['rho'] = 011 VaR.loc[0,'rho'] = rho_first121314 # 記錄組合的模拟收益15 data_R = pd.DataFrame(index = range(MC))1617 # 模拟組合中各資産的标準化收益率18 z = np.random.normal(0,1,size = (2,MC))19 gamma_half = np.array([[1,0],[VaR.loc[0,'rho'],(1-VaR.loc[0,'rho']**2)**0.5]])20 # 轉換為指定相關系數的序列21 z = np.dot(gamma_half,z)2223 # 各資産的真實收益率 = 标準化收益率*波動率24 r = np.multiply(sigma2_Garch**0.5,z)2526 # 兩資産等權重，總收益率27 data_R['R_day1'] = np.mean(r,axis = 0)2829 if n==1:3031 R_ndays = data_R.iloc[:,:2].sum(axis = 1)32 VaR.loc[0,'VaR'] = - np.percentile(R_ndays,p) 33 return VaR3435 # 若計算向前多天的VaR,通過DCC,GARCH模型更新參數3637 alpha_garch = np.dot(np.diag(alpha_garch),np.ones([2,MC]))38 beta_garch = np.dot(np.diag(beta_garch),np.ones([2,MC]))39 omega_garch = np.dot(np.diag(omega_garch),np.ones([2,MC]))40 theta_garch = np.dot(np.diag(theta_garch),np.ones([2,MC])) 41 # 波動率更新 - Garch模型42 sigma2_Garch = omega_garch + np.multiply(alpha_garch,(r - np.multiply(theta_garch,sigma2_Garch**0.5)))**2 + np.multiply(beta_garch,sigma2_Garch)4344 # 相關系數更新 - DCC模型45 q11 = pd.DataFrame(index = range(MC))46 q12 = pd.DataFrame(index = range(MC))47 q22 = pd.DataFrame(index = range(MC))48 q11['day1'] = 1 + alpha_DCC*((z[0,]**2)-1) 49 q12['day1'] = VaR.loc[0,'rho'] + alpha_DCC*(z[0,]*z[1,]- VaR.loc[0,'rho']) + beta_DCC*(VaR.loc[0,'rho'] - VaR.loc[0,'rho'])50 q22['day1']= 1 + alpha_DCC*((z[1,]**2)-1)5152 VaR.loc[1,'rho'] = np.mean(q12['day1']/((q11['day1']*q22['day1'])**0.5))53 # gamma更新54 gamma_half = np.array([[1,0],[VaR.loc[1,'rho'],(1-VaR.loc[1,'rho']**2)**0.5]])555657 for i in range(2,n):58 # 生成不相關的随機序列，并轉換為指定相關系數的随機數59 z = np.random.normal(0,1,size = (2,MC))60 z = np.dot(gamma_half,z)61 r = np.multiply(sigma2_Garch**0.5,z)62 data_R['R_day' + str(i)] = np.mean(r,axis = 0) 6364 sigma2_Garch = omega_garch + np.multiply(alpha_garch,(r - np.multiply(theta_garch,sigma2_Garch**0.5)))**2 + np.multiply(beta_garch,sigma2_Garch)656667 q11['day'+str(i)] = 1 + alpha_DCC*((z[0,]**2)-1) + beta_DCC*(q11['day'+str(i-1)]-1)68 q12['day'+str(i)] = VaR.loc[0,'rho'] + alpha_DCC*(z[0,]*z[1,]- VaR.loc[0,'rho'])+beta_DCC*(q12['day' + str(i-1)] - VaR.loc[0,'rho'])69 q22['day'+str(i)] = 1 + alpha_DCC*((z[1,]**2)-1) + beta_DCC*(q22['day'+str(i-1)]-1)7071 VaR.loc[i,'rho'] = np.mean(q12['day'+str(i)]/((q11['day'+str(i)]*q22['day'+str(i)])**0.5))72 gamma_half = np.array([[1,0],[VaR.loc[i,'rho'],(1-VaR.loc[i,'rho']**2)**0.5]])7374 for i in range(n):75 R_ndays = data_R.iloc[:,:i+1].sum(axis = 1)76 VaR.loc[i,'VaR'] = - np.percentile(R_ndays,p) 7778 return(VaR)
04
總結
寫的有點多，總共三個模型，Monte-Garch，Monte-FHS，Monte-DCC-Garch，前兩個模型的實證結果與文獻結果一緻，大可放心，組合的實證文獻中沒有給出實證結果，是自己算出來的，不保證正确性，不過僅看一日的是和DCC一緻的，是以大機率沒什麼問題，如果覺得有問題，請和我聯系，謝謝。
此外還想說明的一點是，計算未來K日的VaR，K如果取的過大，結果沒有多大實際意義，越往後越趨近于平穩，因為在模型疊代過程中，每一天的收益率都是模拟的誤差不會變化太多，但是波動率和相關系數的誤差會随着預測天數增加逐漸增大，是以還是關注比較近的幾天比較有意義，比如，可以看一下K=500時，DCC模型得到的相關系數變化情況，不同曲線是分别取不同初始值時的結果，越往後預測相關性越小。
最後，本文所有的模型依然是在正态性假設下，沒有考慮非正态情況，非正态的情況下需要用到copula模型，在下一篇文章中給出。
參考文獻
Christoffersen, Peter F. Elements of financial risk management. Academic Press, 2011.
Python愛好者社群曆史文章大合集：
福利：文末掃碼立刻關注公衆号，“Python愛好者社群”，開始學習Python課程：
關注後在公衆号内回複“課程”即可擷取：
小編的Python入門免費視訊課程！！！
【最新免費微課】小編的Python快速上手matplotlib可視化庫！！！
崔老師爬蟲實戰案例免費學習視訊。
陳老師資料分析報告制作免費學習視訊。
玩轉大資料分析！Spark2.X+Python 精華實戰課程免費學習視訊。