首先对VaR的定义做一回顾，上一篇提到，如果我们假设资产标准化的收益率符合正态分布，那么VaR的理论表达式为
上式右边第一项为资产收益率的波动率，第二项为正态分布分布函数的逆函数在p处的值。事实上，即使资产标准化的收益率不符合正态分布，只要右边第二项表示资产标准化的收益率的分布，等式依然成立。
因此，对VaR的估计可以分为两部分，对于波动率的估计和对于资产分布函数的估计。之前三种方法都假设资产标准化的收益率服从标准正态分布，因此只对波动率进行建模，HS和WHS方法过于简单不说，RM方法实际上可以看成是简化版的Garch模型。准确的说，是基于正态分布假设的Garch(1,1)模型的简化版。
这一篇的三种方法主要思路也是围绕上式中波动率估计和资产分布函数的估计展开。
1. 波动率估计
因为文章所用的数据中已经给出了波动率的估计结果，所以实际建模中这部分的模型就略过了。而实际上所用数据波动率的估计结果是书中第四章习题的结果，波动率的估计是通过带VIX指数项的Garch模型，要深究的话可以运行下面图片中的代码，得到的波动率结果和数据给出的是一致的。
因此本文介绍的三种方法主要是对标准化的收益率z的分布进行估计。
2. CF模型
CF（Cornish-Fisher）方法认为将标准化收益率的分布假设为正态分布是不太合理的，省略了过多信息。
因此直接用正态分布函数代替标准化收益率的分布D（0,1）误差会很大， 所以采用了类似泰勒展开的方法：
这种方法可以看做是在正态分布附近的泰勒展开到三次，
，
分别是标准化资产收益率的偏度和峰度。如果将他们看做为0，等价于假设标准化资产收益率的分布是标准正态分布。
用这种方法估计VaR的话，只需要多估计峰度和偏度两个参数，直接用矩估计的方法估计即可。（代码统一附在最后）
3. Garch模型
garch模型不过多介绍，资料网上一大堆，或者直接调用python的arch包，即使不懂原理也可以直接用，Garch模型实际上是对收益率的波动率进行建模。
Garch模型的参数估计一般采用极大似然估计方法（MLE）或者似极大似然方法（QMLE），对VaR问题来说，二者差别不大。这里我们采用两种Garch模型进行建模，一种是常见的基于正态性假设的Garch模型，一种是基于t分布的Garch模型。
基于正态性假设的Garch模型
模型的似然函数可以表示为
极大似然估计方法即最大化似然函数或对数似然函数。
python中可以利用optimize函数计算函数的极值。
基于t分布的Garch模型
金融数据分布最普遍的一个性质是尖峰厚尾，这使得用正态分布去拟合数据误差很大，因此大佬们找出了各种各样神奇的分布去代替正态分布拟合数据，之前提到的CF模型是一种方法，这里讲的t分布是另一种方法。
t分布有一个自由度参数d，随着d的增大，t(d)越来越逼近于正态分布。这里我们用基于t分布的Garch模型进行建模，首先推导Garch模型的极大似然估计方法。
t分布的密度函数为
z是标准化的收益率
z的密度函数为
其中，
在t分布假设下
对数似然函数为
因此，只需要求解使上述函数取最大值的d，这也可以通过optimize函数完成。
读入数据估计参数
代码语言：javascript复制import os
os.chdir('F:\python_study\python_friday\Aug')
import pandas as pd
import numpy as np
from scipy import optimize
from scipy.special import gammaln
import statsmodels.api as sm
from scipy import stats
from scipy.stats import norm,t
import matplotlib.pyplot as plt
data5 = pd.read_excel('Chapter6_Data.xls',sheetname='Question 6.5')
data5['return'] = np.log(1 data5.Close.pct_change(1))
data5['sigma'] = data5['sigma2']**0.5
data5['z'] = data5['return']/data5['sigma']
data5 = data5.dropna()
data5 = data5.reset_index(drop = True)
data_train = data5.loc[0:2008]
data_test = data.loc[2009:]这里为了比较各种方法的优劣，将数据分为两段，data_train和data_test,训练集用来估计参数计算VaR，测试集用于之前相同的方法构建投资组合计算收益。
optimize函数只能求函数的最小值，因此根据之前的公式定义负对数似然函数，最小化负对数似然函数。
代码语言：javascript复制def getNegativeLoglikelihood(d,r):
LogLikeLihood = r.shape[0]*(gammaln((d 1)/2) - gammaln(d/2) -np.log(np.pi)/2 - np.log(d-2)/2) - 1/2*(1 d)*np.log(r['z']**2/(d-2) 1).sum()
return -LogLikeLihood d_best = optimize.fmin(getNegativeLoglikelihood,np.array([10]),
args=(data_train,),ftol = 0.000000001)
print('QMLE估计结果为：',d_best)给定d的初值，d的极大似然估计结果为11.81，我们通过qq图来看参数的效果。
QQ-plot
代码语言：javascript复制s = sm.ProbPlot(data5['z'], dist=stats.t, distargs=(d_best[0],),a = 0.5)
x = s.qqplot(line = '45')
plt.xlim(-8,8)
plt.ylim(-8,8)
plt.show()
QQ图结果来看，参数效果还可以，估计出d之后，就可以计算VaR
除了用t分布之后，还可以用渐近t分布进行估计，方法类似，不再赘述。
4.EVT模型
EVT模型也可以称为极值理论，是风险模型中一类比较重要的方法，这种方法的思想也很有意思，它认为，风险的发生与否，与绝大多数数据是没有关系的，只与那些小概率事件即出现大波动的情况有关系，也就是尖峰厚尾的尾部，因此，我们可以过滤掉那些不重要的数据，只用它的尾部去建模，这样刻画风险的效果可能会更好。极值理论认为尾部既不服从正态分布也不服从t分布，而是服从广义帕累托分布（Generalized Pareto Distribution）
这里的
是GPD中唯一需要估计的参数，也可以通过MLE进行估计。
在估计
之前，还需要做的一件事情是定义尾部，设定一个阈值u，超出这个阈值就认为这是尾部，不超出就认为它不是尾部，估计GPD分布中参数的值只从尾部中取，即只从大于等于u的值中取。
这里u的选取是一个典型的bias and variance trade-off，如果u过大，尾部的值非常少，这样估计得到的
很不稳定，如果u很小，尾部值非常多，不符合EVT模型的假设，这样得到的结果是有偏的。 这个教材里也没有提很高大上的估计方法，只是给出了一个经验估计方法：选取的u保证尾部的值有50个即可。
而GPD分布的参数在MLE下有很简单的解析解表达式如下图，推导过程中运用了Hill estimator，过程略，细节见教材。
从而，极值部分的分布函数可以表示为
其中，T是样本总量，Tu是极值个数，本文中就是50。
从而，VaR可以表示为
其中，
EVT-QQ plot
首先估计参数
代码语言：javascript复制T = data_train.shape[0]
Tu = 50
u = data_train.z.sort_values().values[Tu]
print('u = ',u)
xi = 1/50*np.log(data_train.z.sort_values().values[:Tu]/u).sum()
print('xi = ',xi)
c = Tu/T*abs(u)**(1/xi)
print('c = ',c)这里的xi就是
，得到结果为
代码语言：javascript复制u = -2.05402902944
xi = 0.204467350277
c = 0.841148514934作出极值理论建模下的QQ图，，因为懒得自己写QQ-plot的函数，而scipy中qqplot函数又必须是scipy中的分布对象，可是scipy中又没有直接可以用的GPD分布函数，所以这里偷了个懒，从scipy中继承了正态分布，再把分布函数改成GPD的分布函数，然后直接调用scipy中的qqplot看结果，过程其实还是很tricky的，对我来说，这里可能是全文最难的点了，其他地方都很快写完了，这里写了蛮久，开始给这个对象定义了好多方法，后来发现其实只需要定义好做qqplot需要用的逆函数就行了。
代码语言：javascript复制from scipy.stats._continuous_distns import norm_gen
class EVT(norm_gen):
def _ppf(self,q):
u = -2.225
xi = 0.215
return u*(q**(-xi))结果
看起来比之前的t分布结果要好很多。
5. VaR
用训练集得到的参数计算VaR，结果如图，其中RM方法为上一篇中介绍过的。
看上去很接近，没什么差别，用与上一篇相同的方法构造投资组合，曲线图如下。（不想去看怎么构造的可以这样想，收益越高，表明估计效果最好）
可以看出，CF的效果最差，EVT和RM的效果差不多是最好的，后期RM的效果超过了EVT，两种Garch方法的结果居中。
本文整篇的流程都是按照作者在原书中的方法进行的，但是在Garch方法这里，感觉逻辑上是有一些问题的，Garch模型只在了对标准化收益的分布函数中，但实际上波动率也应该用估计得到的波动率替换，这样才是配套的，但作者使用的一直是另外一种方法得到的波动率，而EVT和CF中没有单独对波动率的估计，所以用别的方法估计波动率没有什么问题。
总的来说，RM方法是最简洁易行的，而且效果最好，EVT方法用到的数据量最少，仅50个，效果也很好。
6. 代码
var计算
代码语言：javascript复制p = 0.01
phi = norm(0,1).ppf(p)
data5['VAR_EVT'] = -data5['sigma']*u*(p/(Tu/T))**(-xi)
data5['VAR_norm'] = - data5['sigma']*phi
data5['VAR_t'] = - data5['sigma']*t(d_best[0]).ppf(p)*((d_best[0] - 2)/d_best[0])**0.5
data5['VAR_CF'] = - data5['sigma']*(phi data5.z.skew()/6*(phi**2-1)
data5.z.kurt()/24*(phi**3-3*phi) - data5.z.skew()**2/36*(2*phi**3-5*phi))
data5['sigma2_RM'] = np.nan
data5.loc[0,'sigma2_RM'] = 0
for i in range(data5.shape[0] - 2):
data5.loc[i 1,'sigma2_RM'] = data5.loc[i,'sigma2_RM']*0.94 0.06* data5.loc[i,'return']**2
data5['VAR_RM'] = -data5['sigma2_RM']**0.5 * phivar作图
代码语言：javascript复制xticklabel = data5.loc[:,'Date'].apply(lambda x:str(x.year) '-' str(x.month) '-' str(x.day))
xticks = np.arange(0,data5.shape[0] 1,np.int((data5.shape[0] 1)/4))
X = np.arange(data5.shape[0])
plt.figure(figsize=(16,5))
SP = plt.axes()
plt.plot(X,data5.VAR_EVT,label = 'VAR_EVT')
plt.plot(X,data5.VAR_norm,label = 'VAR_EVT')
plt.plot(X,data5.VAR_t,label = 'VAR_t(d)')
plt.plot(X,data5.VAR_CF,label = 'VAR_CF')
plt.plot(X,data5.VAR_RM,label = 'VAR_RM')
SP.set_xticks(xticks)
SP.set_xticklabels(xticklabel[xticks])
plt.grid()
plt.legend()
plt.show()组合收益计算
代码语言：javascript复制data51 = data5.loc[2009:].copy()
data51 = data51.reset_index(drop = True)
data51['position_EVT'] = 100000/(data51['VAR_EVT']*np.sqrt(10))
data51['position_norm'] = 100000/(data51['VAR_norm']*np.sqrt(10))
data51['position_t'] = 100000/(data51['VAR_t']*np.sqrt(10))
data51['position_CF'] = 100000/(data51['VAR_CF']*np.sqrt(10))
data51['position_RM'] = 100000/(data51['VAR_RM']*np.sqrt(10))
data51['return_EVT'] = (np.exp(data51['return']) - 1)*data51['position_EVT']
data51['return_norm'] = (np.exp(data51['return']) - 1)*data51['position_norm']
data51['return_t'] = (np.exp(data51['return']) - 1)*data51['position_t']
data51['return_CF'] = (np.exp(data51['return']) - 1)*data51['position_CF']
data51['return_RM'] = (np.exp(data51['return']) - 1)*data51['position_RM']
data51['cum_EVT'] = data51['return_EVT'].cumsum()
data51['cum_norm'] = data51['return_norm'].cumsum()
data51['cum_t'] = data51['return_t'].cumsum()
data51['cum_CF'] = data51['return_CF'].cumsum()
data51['cum_RM'] = data51['return_RM'].cumsum()组合收益曲线
代码语言：javascript复制X = np.arange(data51.shape[0])
xticklabel = data51.loc[:,'Date'].apply(lambda x:str(x.year) '-' str(x.month) '-' str(x.day))
xticks = np.arange(0,data51.shape[0],np.int((data51.shape[0])/5))
plt.figure(figsize = [20,5])
SP = plt.axes()
SP.plot(X,data51['cum_EVT'],label = 'EVT')
SP.plot(X,data51['cum_norm'],label = 'Garch-norm')
SP.plot(X,data51['cum_t'],label = 'Garch-t')
SP.plot(X,data51['cum_CF'],label = 'CF')
SP.plot(X,data51['cum_RM'],label = 'RM')
SP.set_xticks(xticks)
SP.set_xticklabels(xticklabel[xticks],rotation=45,size =20)
plt.legend()
plt.grid()
plt.show()参考文献
《Elements of Financial Risk Management》
