【FinE】正态分布和t分布下的CVaR

CVaR and VaR

CVaR(Conditional Value-at-Risk)也被称为Expected Shortfall(ES) 或者 Expected Tail Loss(ETL)，可以解释超过给定VaR值的期望损失，在很多风险分析中，CVaR是更具有解释力，当设置了VaR阈值后，CVaR说明了在之后的$h$天交易中的期望损失风险（即资产价值低于预先设置的VaR值）. 显然VaR和资产日回报率的分布有关.

当资产回报率的尾部较厚，使用整态分布刻画是不合适的，考虑使用t分布对资产回报率分布进行描述.

Model

normal distribution

令$X$表示$h$日的回报率，则
$$Va{R}_{h,\alpha }=-x_{h,\alpha }$$
其中，$P(X<x_{h,\alpha })=\alpha$.
CVaR使用条件期望的形式表达
$$CVa{R}_{h,\alpha }(X)=-\mathbb{E}(X\mid X<x_{h,\alpha })=-{\alpha }^{-1}{\int }_{-\infty }^{x_{h,\alpha }}=xf(x)dx$$
对于任意连续的概率密度函数$f(x)$，需要计算$x\sim 100(1-\alpha )\%h$天VaR的$xf(x)$积分.
考虑在$X\sim N({\mu }_h,{\sigma }_h^2)$的情况下的CVaR
$$CVa{R}_{h,\alpha }(X)={\alpha }^{-1}\phi ({\Phi }^{-1}(\alpha )){\sigma }_h-{\mu }_h$$
其中$\phi (z)$表示标准正态分布的pdf，$\Phi (\alpha {)}^{-1}$为$\alpha$分位数.
案例：计算$h=5$的CVaR，股价年化回报率符合$\sigma =41\%,\mu =0$的正态分布，设置一年的交易日为252天计算出
$${\sigma }_h=\sigma \sqrt{h/252}=0.41\sqrt{5/252}=0.05798$$
python计算gaussian分布VaR和CVaR代码如下

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from math import sqrt

def demo_1():
    h, D=5, 252
    mu, sig = 0, 0.41
    muh, sigh=mu*(h/D), sig*sqrt(h/D)
    alpha=0.01
    lev=100*(1-alpha)

    CVaRh=alpha**(-1)*norm.pdf(norm.ppf(alpha))*sigh-muh
    VaRh=norm.ppf(1-alpha)*sigh-muh

    print('{}% {} day Normal VaR= {:.2f}%'.format(lev, h, VaRh*100))
    print('{}% {} day Normal CVaR= {:.2f}%'.format(lev, h, CVaRh*100))

demo_1()

student t distribution

考虑更适合肥尾的t分布进行建模，标准t分布的pdf函数为
$$f_{\nu }(x)=\underset{A}{\underbrace{((\nu -2)\pi {)}^{-1/2}\Gamma (\frac{v}{2}{)}^{-1}\Gamma (\frac{\nu +1}{2})}}(1+\underset{a}{\underbrace{(\nu -2{)}^{-1}}}{x}^2{)}^{\underset{b}{\underbrace{-(\nu +1)/2}}}$$
其中$\Gamma$表示gamma函数，$\nu$表示自由度，方程一般形式为
$$f_{\nu }(x)=A(1+a{x}^2{)}^b$$
代入到CVaR定义中得到
$$\begin{array}{rl}CVa{R}_{\alpha ,\nu }& =-{\alpha }^{-1}{\int }_{-\infty }^{x_{\alpha ,\nu }}x{f}_{\nu }(x)dx\\ & =-{\alpha }^{-1}{\int }_{-\infty }^{x_{\alpha ,\nu }}xA(1+a{x}^2{)}^{b}dx\\ & =-\frac{A}{\alpha }{\int }_{-\infty }^{x_{\alpha ,\nu }}x(\underset{y}{\underbrace{1+a{x}^2}}{)}^{b}dx\end{array}\text{\hspace{1em}(1)}$$
令$y=1+a{x}^2$进行变换，$dy=2axdx，B=1+(\nu -2{)}^{-1}{x}_{\alpha ,\nu }^2$，方程$(1)$变换得到
$$\begin{array}{rl}(1)& =-\frac{A}{\alpha }{\int }_{-\infty }^B\frac{x(1+a{x}^2{)}^b}{2ax}dy\\ & =-\frac{A}{2a\alpha }{\int }_{\infty }^B{y}^{b}dy\\ & =-\frac{A}{2a\alpha }\frac{B^{b+1}}{b+1}\end{array}$$
由于$b+1=(1-\nu )/2$代入得到$(1)$值为
$$-\frac{A}{(\nu -2{)}^{-1}\alpha }\frac{2B^{(1-\nu )/2}}{1-\nu }$$
根据$f_{\nu }(x)$的表达式可以知道
$$f_{\nu }(x)=A{B}^b\Rightarrow A=f_{\nu }(x)B^{-b}=f_{\nu }(x)B^{(1+\nu )/2}$$
代入方程$(1)$得到CVaR值为
$$\begin{array}{rl}(1)& =-{\alpha }^{-1}\frac{f_{\nu }(x)B^{(1+\nu )/2}}{2(\nu -2{)}^{-1}}\frac{2B^{(1-\nu )/2}}{1-\nu }\\ & =-{\alpha }^{-1}{f}_{\nu }(x)(\nu -2)(1-\nu {)}^{-1}B\\ & =-{\alpha }^{-1}(\nu -2)(1-\nu {)}^{-1}(1+(\nu -2{)}^{-1}{x}_{\alpha ,\nu }^2)f_{\nu }(x_{\alpha ,\nu })\\ & =-{\alpha }^{-1}(1-\nu {)}^{-1}[\nu -2+x_{\alpha ,\nu }^2]f_{\nu }(x_{\alpha ,\nu })\end{array}$$
所以，$h$天，t分布下CVaR值为
$$CVa{R}_{h,\alpha ,\nu }(X)=-{\alpha }^{-1}(1-\nu {)}^{-1}[\nu -2+x_{\alpha ,\nu }^2]f_{\nu }(x_{\alpha ,\nu }){\sigma }_h-{\mu }_h$$
案例计算自由度为$6$的t分布的VaR和CVaR
python计算t分布下VaR和CVaR代码如下

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, t
from math import sqrt

def demo_2():
    h, D = 10, 252
    mu, sig=0, 0.41
    muh, sigh=mu*(h/D), sig*sqrt(h/D)
    alpha=0.01
    lev=100*(1-alpha)
    nu=6 # 设置自由度
    xanu=t.ppf(alpha, nu)
    CVaRh=-1/alpha*(1-nu)**(-1)*(nu-2+xanu**2)*t.pdf(xanu, nu)*sigh-muh
    VaRh=sqrt(h/D*(nu-2)/nu)*t.ppf(1-alpha, nu)*sig-mu
    
    print('{}% {} day t VaR= {:.2f}%'.format(lev, h, VaRh*100))
    print('{}% {} day t CVaR= {:.2f}%'.format(lev, h, CVaRh*100))

demo_2()

可以发现，随着自由度的上升，t分布的VaR和CVaR逐渐收敛到gaussian分布的VaR和CVaR

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, t
from math import sqrt


def demo_3():
    h, D = 10, 252
    mu, sig=0, 0.41
    muh, sigh=mu*(h/D), sig*sqrt(h/D)
    alpha=0.01
    lev=100*(1-alpha)

    data=[]
    for nu in range(5, 101):
        xanu=t.ppf(alpha, nu)
        CVaRt=-1/alpha*(1-nu)**(-1)*(nu-2+xanu**2)*t.pdf(xanu, nu)*sigh-muh
        VaRt=sqrt(h/D*(nu-2)/nu)*t.ppf(1-alpha, nu)*sig-muh
        data.append([nu, VaRt, CVaRt])

    CVaRn=alpha**(-1)*norm.pdf(norm.ppf(alpha))*sigh-muh
    VaRn=norm.ppf(1-alpha)*sigh-muh

    data=np.array(data)

    fig, ax=plt.subplots(figsize=(8, 6))
    plt.plot(data[:, 0], data[:, 1]*100, 'b-', label='VaRt')
    plt.plot(np.arange(5, 100), VaRn*np.ones(95)*100, 'b:', label='VaRn')
    plt.plot(data[:, 0], data[:, 2]*100, 'r-', label='CVaRt')
    plt.plot(np.arange(5, 100), CVaRn*np.ones(95)*100, 'r:', label='CVaRn')

    plt.xlabel('student t. d.o.f')
    plt.ylabel('%')
    plt.legend()
    plt.show()

    
demo_3()

Case Study

使用IBM数据进行实证分析， 使用gaussian分布和t分布拟合IBM日回报收益率得到拟合图像如下

可以发现，t分布更适合拟合这种存在高噪声的尖峰厚尾的数据，不同分布计算的CVaR和VaR值如下

99.0% 1 day gaussian VaR= 2.81%
99.0% 1 day gaussian CVaR= 3.22%
99.0% 1 day student t VaR= 2.22%
99.0% 1 day student t CVaR= 3.91%

案例python代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, t, skew, kurtosis, kurtosistest
from math import sqrt
import pandas_datareader.data as web
import pickle

# Fetching Yahoo Finance for IBM stock data
def data_fetch():
    data = web.DataReader('IBM', data_source='yahoo', start='2010-12-31', end='2015-12-31')['Adj Close']
    adj_close=np.array(data.values)
    ret=adj_close[1:]/adj_close[:-1]-1
    file=open('ibm_ret', 'wb')
    pickle.dump(ret, file)

with open('ibm_ret', 'rb') as file:
    ret=pickle.load(file)

dx=0.0001
x=np.arange(-0.1, 0.1, dx)

# gaussian fit
mu_norm, sig_norm=norm.fit(ret)
pdf_norm=norm.pdf(x, mu_norm, sig_norm)

# student t fit
nu, mu_t, sig_t=t.fit(ret)
nu=np.round(nu)
pdf_t=t.pdf(x, nu, mu_t, sig_t)

# VaR and CVaR
h=1
alpha=0.01
lev=100*(1-alpha)
xanu=t.ppf(alpha, nu)

CVaRn=alpha**(-1)*norm.pdf(norm.ppf(alpha))*sig_norm-mu_norm
VaRn=norm.ppf(1-alpha)*sig_norm-mu_norm
CVaRt=-1/alpha*(1-nu)**(-1)*(nu-2+xanu**2)*t.pdf(xanu, nu)*sig_t-h*mu_t
VaRt=sqrt((nu-2)/nu)*t.ppf(1-alpha, nu)*sig_t-h*mu_t

print('{}% {} day gaussian VaR= {:.2f}%'.format(lev, h, VaRn*100))
print('{}% {} day gaussian CVaR= {:.2f}%'.format(lev, h, CVaRn*100))
print('{}% {} day student t VaR= {:.2f}%'.format(lev, h, VaRt*100))
print('{}% {} day student t CVaR= {:.2f}%'.format(lev, h, CVaRt*100))

plt.figure(figsize=(12, 8))
grey=0.75, 0.75, 0.75
plt.hist(ret, bins=50, density=True, color=grey, edgecolor='none')
plt.axis('tight')
plt.plot(x, pdf_norm, 'b-.', label='Guassian Fit')
plt.plot(x, pdf_t, 'g-.', label='Student t Fit')
plt.xlim([-0.2, 0.1])
plt.ylim([0, 50])
plt.legend(loc='best')
plt.xlabel('IBM daily return')
plt.ylabel('ret. distr.')

# plt.savefig('ibm_fit')
# inset local
sub=plt.axes([0.22, 0.35, 0.3, 0.4])
plt.hist(ret, bins=50, density=True, color=grey, edgecolor='none')
plt.plot(x, pdf_norm, 'b')
plt.plot(x, pdf_t, 'g')
plt.plot([-CVaRt, -CVaRt], [0, 3], 'g:')
plt.plot([-CVaRn, -CVaRn], [0, 4], 'b:')
plt.text(-CVaRt-0.015, 3.1, 'Stu. t', color='g')
plt.text(-CVaRn-0.015, 4.1, 'Gaussian CVaR', color='b')
plt.xlim([-0.09, -0.02])
plt.ylim([0, 5])

# plt.savefig('ibmcvar')
plt.show()

从数据和图像可以看出，Gaussian拟合计算的CVaR值存在低估风险的倾向.

Reference

CVaR Quant at risk