python编程估计Copula并计算拟合优度
前⾔：Use Copula to define the correlation structure between two variables and compute their joint distribution.
⽬的：计算Gaussian Copula、Gumbel Copula、Clayton Copula
简单思路介绍和步骤：
1. Gaussian Copula
假设要模拟两个时间序列｛x｝、｛y｝之间的相关性结构，定义Gaussian Copula如下：
其中u1、u2分别为对应变量分布函数的值。
u1 = F1(x)
u2 = F2(y)
⽽其中的参数p为待估计变量，也就是待检验的x与y的分布函数之间的相关系数，可以看成x与y之间独⽴性的衡量指标。本⽂采⽤模拟估计的⽅法进⾏求解。‘简单来说，就是从[-1, 1]按⼀定步长取p值，计算出联合分布函数，检验其对应的联合分布函数与实际样本的经验联合分布函数之间的偏差，如果偏差在⾜够⼩的范围内，就认为其拟合程度达到要求。
2. Gumbel Copula
3. Clayton Copula
2、3的估计⽅法与1相同。下⾯，从2003年1⽉6⽇⾄2019年2⽉15⽇，取3202个有效数据，计算中国与美、德、法、英、⽇五个
国家的Copula系数值，⽐较相关性⼤⼩。数据来⾃国泰安CSMAR数据库。
Gaussian Copula相关系数p：代表被分析的两个分布函数之间的相关程度。
Gumbel theta：代表两个分布之间的独⽴性程度，theta>=1，并且越接近1，表明两个分布的独⽴性越强，反之则越弱。
Clayton theta：同样代表两个分布的独⽴程度，theta越接近于0，代表两个分布函数之间的独⽴性越强，反之则越弱。
代码如下：
import math
import numpy as np
import xlrd
from scipy.integrate import tplquad, dblquad, quad
#⽤于计算Gaussian copula的函数值. 输⼊x, y以及相关系数p, 得出copula值C(F(x), G(y))
def Gaussian_integrate(marginal_x, marginal_y, p):
Gaussian_copula = []
for i in range(len(marginal_x)):
val2,err2 =dblquad(lambda y,x:1/(2*np.pi*math.sqrt(1-p*p))*pow(math.e,-(x*x-2*p*x*y+y*y)/(2*(1-p*p))),float("-inf"),marginal_x[i],float("-inf"),marginal_y[ i])
#print(val2, '\t', marginal_x[i], '\t', marginal_y[i])
Gaussian_copula.append(val2)
return Gaussian_copula
#Student-t分布
def Student_t_integrate(marginal_x, marginal_y, p):
for i in range(len(marginal_x)):
val2,err2 =dblquad(lambda y,x:1/(2*np.pi*math.sqrt(1-p*p))*pow(math.e,1+(x*x-2*p*x*y+y*y)/(2*(1-p*p))),float("-inf"),marginal_x[i],float("-inf"),marginal_ y[i])
print(val2, '\t', marginal_x[i], '\t', marginal_y[i])
print(val2, '\t', marginal_x[i], '\t', marginal_y[i])
#计算copula的拟合精度
def goodness_of_fit(theory_copula, empirical_copula):
error = 0
for i in range(len(theory_copula)):
error += pow(theory_copula[i] - empirical_copula[i], 2)
return error
#计算联合分布函数
def cumulative_function(marginal_x, marginal_y):
empirical_copula = []
for i in range(len(marginal_x)):
count = 0
for j in range(len(marginal_y)):
if marginal_x[j] <= marginal_x[i] and marginal_y[j] <= marginal_y[i]:
count += 1
empirical_copula.append(count/len(marginal_x))
return empirical_copula
def Gumbel_copula(marginal_x, marginal_y, cumulative_x, cumulative_y, theta, probability):
probability_x = []
lambda编程probability_y = []
for i in range(len(marginal_x)):
for i2 in range(len(cumulative_x)):
if cumulative_x[i2] == marginal_x[i]:
probability_x.append(probability[i2])
for j in range(len(marginal_y)):
for j2 in range(len(cumulative_y)):
if cumulative_y[j2] == marginal_y[j]:
probability_y.append(probability[j2])
Gumbel_copula = []
for i in range(len(probability_x)):
G_copula_x = -math.log(probability_x[i], math.e)
G_copula_y = -math.log(probability_y[i], math.e)
#计算Gumbel_copula的对应值
Gumbel_copula.append(pow(math.e, -pow((pow(G_copula_x, theta)+pow(G_copula_y, theta)), 1/theta))) return Gumbel_copula
def Clayton_copula(marginal_x, marginal_y, cumulative_x, cumulative_y, theta, probability):
probability_x = []
probability_y = []
#将marginal_x与对应的分布函数概率统⼀起来, ⽅便查得到u1, u2(对应边际分布函数概率)
for i in range(len(marginal_x)):
for i2 in range(len(cumulative_x)):
if cumulative_x[i2] == marginal_x[i]:
probability_x.append(probability[i2])
for j in range(len(marginal_y)):
for j2 in range(len(cumulative_y)):
if cumulative_y[j2] == marginal_y[j]:
probability_y.append(probability[j2])
Clayton_copula = []
for i in range(len(probability_x)):
Clayton_copula.append(pow((pow(probability_x[i], -theta)+pow(probability_y[i], -theta)-1), -1/theta))
return Clayton_copula
if __name__=="__main__":
data = xlrd.open_workbook(".\\中国与亚洲国家风险相关性研究.xlsx")
#data = xlrd.open_workbook(".\\中国股市内部相关性分析.xlsx")
#data = xlrd.open_workbook(".\\中欧美⽇-风险相关性研究.xlsx")
#data = xlrd.open_workbook(".\\中欧美⽇-风险相关性研究-近3年.xlsx")
#data = xlrd.open_workbook(".\\中欧美⽇-风险相关性研究-前3年.xlsx")
table = data.sheets()[4]
table = data.sheets()[4]
#⽤于计算Gaussian_copula
print("/*****************************Gaussian_copula*****************************")
marginal_y = l_values(8)
"""
for k in range(1,8):
marginal_x = l_values(k)
empirical_copula = cumulative_function(marginal_x, marginal_y)
for i in range(-9, 10):
print(i/10, goodness_of_fit(Gaussian_integrate(marginal_x,marginal_y,i/10), empirical_copula))
for i in range(1, 28):
marginal_y = l_values(i)
for j in range(i+1, 28):
marginal_x = l_values(i+1)
empirical_copula = cumulative_function(marginal_x, marginal_y)
for i in range(-9, 10):
print(goodness_of_fit(Gaussian_integrate(marginal_x,marginal_y,i/10), empirical_copula), i, j)
"""
#marginal_x⽤于表⽰实际的数据
marginal_x = l_values(1)
#table_CDF⽤于查询分布函数概率
table_CDF = data.sheets()[6]
cumulative_x = l_values(1)
cumulative_y = l_values(8)
probability = l_values(9)
#⽤于计算Gumbel_copula
print("/*****************************Gumbel_copula*****************************/")
for j in range(2, 3):
marginal_x = l_values(j)
cumulative_x = l_values(j)
for i in range(100, 180):
print(goodness_of_fit(Gumbel_copula(marginal_x, marginal_y, cumulative_x, cumulative_y, i/100, probability), cumulative_function(marginal_x, ma rginal_y)), "\t", i/100)
"""
#⽤于计算Clayton_copula
print("/*****************************Clayton_copula*****************************/")
for j in range(2, 3):
marginal_x = l_values(j)
cumulative_x = l_values(j)
for i in range(1, 200):
print(goodness_of_fit(Clayton_copula(marginal_x, marginal_y, cumulative_x, cumulative_y, i/100, probability), cumulative_function(marginal_x, mar ginal_y)), "\t", i/100)
"""
结论：
Gumbel和Clayton的拟合效果相当，并且明显⾼于Gaussian Copula。
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论