|
前面的推文介绍了各种关于资产组合收益率的模型，有助于主动投资的风险评估，但对于实际中我们新购买一支股票对于我们已有的资产组合有何影响？我们如何进行风险决策，我们还需要研究各独立资产的相关性。本部分内容共分三部分推出：一、构造关于动态协方差和相关系数的模型（假设服从正态分布）；二、使用模特卡罗等模拟工具进行多期风险评估；三、关于使用copula models 研究非正态分布和非线性相关性的风险评估【本文主要内容】。
关于资产间的相关性可借用蒙特卡罗方法进行模拟，代码如下： temp_1=pd.read_excel('DCC_Shocks_1.xlsx')
DCC_Shocks_1 = np.array(temp_1)
DCC_Shocks_1.shape
temp_2=pd.read_excel('DCC_Shocks_2.xlsx')
DCC_Shocks_2 = np.array(temp_2)
DCC_Shocks_2.shape
###### initial some paramters
def DCC_sim(alpha,beta,uncond_var,rho0):
day_p_list=[]
x1= [[float(1) for i in range(10)] for i in range(10000)]
q11=np.array(x1)
q12=np.array(x1)
q22=np.array(x1)
rho12=np.array(x1)
for i in range(10):
if i==0:
q11[:,i]=(1-alpha-beta)+alpha*(DCC_Shocks_1[:,i]**2)+beta*1
q12[:,i]=(1-alpha-beta)*uncond_var+alpha*DCC_Shocks_1[:,i]*DCC_Shocks_2[:,i]+beta*rho0
q22[:,i]=(1-alpha-beta)+alpha*(DCC_Shocks_2[:,i]**2)+beta*1
rho12[:,i]=q12[:,i]/(np.sqrt(q11[:,i]*q22[:,i]))
day_p_list.append(rho12[:,i].mean())
else:
q11[:,i]=(1-alpha-beta)+alpha*(DCC_Shocks_1[:,i]**2)+beta*q11[:,i-1]
q12[:,i]=(1-alpha-beta)*uncond_var+alpha*DCC_Shocks_1[:,i]*DCC_Shocks_2[:,i]+beta*q12[:,i-1]
q22[:,i]=(1-alpha-beta)+alpha*(DCC_Shocks_2[:,i]**2)+beta*q22[:,i-1]
rho12[:,i]=q12[:,i]/(np.sqrt(q11[:,i]*q22[:,i]))
day_p_list.append(rho12[:,i].mean())
return day_p_list
alpha,beta=[0.0676077250447349,0.895472849691596]
uncond_var=-0.306838867178967
#p_list=DCC_sim(alpha,beta,uncond_var,0.5)
p_dic={}
a=1
for rho0 in [-0.5,-0.25,0.25,0.5]:
p_list=DCC_sim(alpha,beta,uncond_var,rho0)
p_dic.update({a:p_list})
a=a+1
result_df= pd.DataFrame(None)
for key,value in p_dic.items():
result_df[key]=value
result_df.to_csv('correlatin.csv')
#### plot ####
data4=pd.read_csv('correlatin.csv')
HS = plt.subplot(111)
ax2 = HS
s1 = ax2.plot(data4['1'], 'red',label = '-0.5')
s2 = ax2.plot(data4['2'], 'yellow',label = '-0.25')
s3 = ax2.plot(data4['3'], 'blue',label = '0.25')
s4 = ax2.plot(data4['4'], 'black',label = '0.5')
s = s1+s2+s3+s4
lns = [l.get_label() for l in s]
plt.legend(s,lns)
plt.show()
结果如下： 
阈值相关：不同百分比下，资产序列间的相关性。对两个资产（SP500和 us_10）的收益率序列利用二元正态分布得出阀值相关性。 代码如下： #### question one
data1=pd.read_excel('question_1.xlsx')
data1['return_sp'] = np.log(data1['Close_sp']/data1['Close_sp'].shift(1))
data1['return_us'] = np.log(data1['Close_us']/data1['Close_us'].shift(1))
data1 = data1.dropna()
data1 = data1.reset_index(drop = True)
### 分位数点 百分位点 左边
return_sp = np.array(data1['return_sp'])
return_us = np.array(data1['return_us'])
sp_list=[]
for per in list(np.linspace(15,50,36)):
a=np.percentile(return_sp,per)
sp_list.append(a)
us_list=[]
for per in list(np.linspace(15,50,36)):
a=np.percentile(return_us,per)
us_list.append(a)
for i in range(len(sp_list)):
aa=data1.loc[(data1['return_sp'] <>'return_us'] <>
bb=data1.loc[(data1['return_sp'] <>'return_us'] <>
data1['sp_dot'+str(i)]=aa['return_sp']
data1['us_dot'+str(i)]=bb['return_us']
### 分位数点 百分位点 右边
sp_list=[]
for per in list(np.linspace(50,85,35)):
a=np.percentile(return_sp,per)
sp_list.append(a)
us_list=[]
for per in list(np.linspace(50,85,35)):
a=np.percentile(return_us,per)
us_list.append(a)
for i in range(len(sp_list)):
aa=data1.loc[(data1['return_sp'] > sp_list[i]) & (data1['return_us'] > us_list[i]),:]
bb=data1.loc[(data1['return_sp'] > sp_list[i]) & (data1['return_us'] > us_list[i]),:]
data1['sp_dot'+str(i+36)]=aa['return_sp']
data1['us_dot'+str(i+36)]=bb['return_us']
data1=data1.fillna(0)
#### 求相关系数
corr_list=[]
for i in range(71):
data_test=data1.loc[data1['sp_dot'+str(i)]!=0]
b= data_test['sp_dot'+str(i)].corr( data_test['us_dot'+str(i)])
corr_list.append(b)
结果如下： 
当我们对收益率用GRACH模型进行处理后，两资产都出现了较大损失之间的相关性比出现正的收益的相关性还大。结果如下： 
这暗示了我们使用二元正态分布来计算资产间的相关性是不合适的。因此我们需要重新建模。
本文主要介绍copula models模型 copula models模型关键在于可以连接我们之前学习的各种单变量模型，本文主要基于t分布的copula models模型。第一、利用极大似然估计t分布的自由度d值；第二，利用极大似然估计连接系数。代码如下： data3= pd.read_excel('question3.xlsx')
data3['re_log_sp5']=data3.close_sp500.pct_change(1)
data3['re_log_us']=data3.close_us_10.pct_change(1)
data3 = data3.dropna()
data3 = data3.reset_index(drop = True)
### sigma2 using RM
##### RiskMetrics ####
def RM_sigm2(data):
data['sigma2_RM_Sp'] = np.nan
return_need = list(data['re_log_sp5'])
var = np.var(return_need)
data.loc[0,'sigma2_RM_Sp'] = var
for i in range(1,(data.shape[0])):
data.loc[i,'sigma2_RM_Sp'] = data.loc[i-1,'sigma2_RM_Sp']*0.94 + 0.06* data.loc[i-1,'re_log_sp5']**2
data['st_return_RM']=data['re_log_sp5']/np.sqrt(data['sigma2_RM_Sp'])
#### SP500
RM_sigm2(data3)
#### US_10y
def RM_sigm2_1(data):
data['sigma2_RM_Us'] = np.nan
return_need = list(data['re_log_us'])
var = np.var(return_need)
data.loc[0,'sigma2_RM_Us'] = var
for i in range(1,(data.shape[0])):
data.loc[i,'sigma2_RM_Us'] = data.loc[i-1,'sigma2_RM_Us']*0.94 + 0.06* data.loc[i-1,'re_log_us']**2
data['st_return_RM_us']=data['re_log_us']/np.sqrt(data['sigma2_RM_Us'])
RM_sigm2_1(data3)
### estimate the d value
def likfunc_td(d,st_return):
'''
inputs: d - t(d) degrees of freedom
zs - standardarized returns
output: logliksum
'''
loglike = []
zs = st_return
for i in range(len(st_return)):
tmp=-scipy.special.gammaln((d+1)/2)+scipy.special.gammaln(d/2)+np.log(np.pi)/2+np.log(d-2)/2 + 0.5*(1+d)*np.log(1+zs[i]**2/(d-2))
loglike.append(tmp)
return np.array(loglike).sum()
##### 使用极大似然函数估计d值 #####
## sp500
arr_sp_return=np.array(data3['st_return_RM'])
params_MLE2 = optimize.fmin(likfunc_td,10, args=(arr_sp_return,), ftol = 0.00000001)
d_sp = params_MLE2[0]
## us_10
arr_us_return=np.array(data3['st_return_RM_us'])
params_MLE3 = optimize.fmin(likfunc_td,10, args=(arr_us_return,), ftol = 0.00000001)
d_us = params_MLE3[0]
### t分布的概率密度
###### Note that ui is simply the probability of observing a value below zi for asset i. ########
def getMu_Inverse_1(st_return):
mu_1_list=[]
inverse_u1_list=[]
for one in np.array(st_return):
print(one)
if one > 0:
mu_1=1-t.pdf(one*np.sqrt(d_sp/(d_sp-2)),d_sp)
mu_1_list.append(mu_1)
else:
mu_1=t.pdf((-1)*one*np.sqrt(d_sp/(d_sp-2)),d_sp)
mu_1_list.append(mu_1)
inverse_u1_list=[]
for one_1 in mu_1_list:
if one_1>0:
phi_inverse = norm.isf(1-one_1, loc=0, scale=1)
inverse_u1_list.append(phi_inverse)
else:
phi_inverse = norm.isf(one_1, loc=0, scale=1)
inverse_u1_list.append(phi_inverse)
data3['mu_1']=mu_1_list
data3['inverse_mu_1']=inverse_u1_list
def getMu_Inverse_2(st_return,d):
mu_list=[]
for one in np.array(st_return):
print(one)
if one > 0:
mu=1-t.pdf(one*np.sqrt(d/(d-2)),d)
mu_list.append(mu)
else:
mu=t.pdf((-1)*one*np.sqrt(d/(d-2)),d)
mu_list.append(mu)
inverse_u_list=[]
for one_1 in mu_list:
if one_1>0:
phi_inverse = norm.isf(1-one_1, loc=0, scale=1)
inverse_u_list.append(phi_inverse)
else:
phi_inverse = norm.isf(one_1, loc=0, scale=1)
inverse_u_list.append(phi_inverse)
data3['mu_2']=mu_list
data3['inverse_mu_2']=inverse_u_list
st_return = np.array(data3['st_return_RM'])
getMu_Inverse_1(st_return)
st_return = np.array(data3['st_return_RM_us'])
d= d_us
getMu_Inverse_2(st_return,d)
###########
b= data3['inverse_mu_1'].corr( data3['inverse_mu_2'])
### 估计相关系数
def likfunc_corr(b,inverse_mu_1,inverse_mu_2):
'''
inputs: inverse_mu_1
inverse_mu_2
b
output: logliksum
'''
loglike = []
for i in range(len(st_return)):
tempt=0.5*(1-b**2)+(inverse_mu_1**2+inverse_mu_2**2-2*inverse_mu_1*inverse_mu_2)/2*(1-b**2)-0.5*(inverse_mu_1**2+inverse_mu_2**2)
loglike.append(tempt)
return np.array(loglike).sum()
##### 使用极大似然函数估计相关系数b值 #####
## sp500
arr_mu_1=np.array(data3['inverse_mu_1'])
arr_mu_2=np.array(data3['inverse_mu_2'])
likfunc_corr = optimize.fmin(likfunc_corr,b, args=(arr_mu_1,arr_mu_2,), ftol = 0.00000001)
corr_2 = likfunc_corr[0]
利用copula models模型估计完的连接系数后，可以结合蒙特卡洛法预测未来的VaR和ES。
参考资料：Elements of Financial Risk Management 资料分享百度云盘（数据和教材）：（第八章后部分以及第九章） 链接: https://pan.baidu.com/s/1K444PvgzJhibDyNqL2EXFg 密码: bwje
【指导老师：首都经济贸易大学 余颖丰；作者：北京第二外国语学院 李燕】 |
