沪深股市收益率的二元Copula模型
- 1. 案例描述
- 2.实现流程
- 2.1 确定随机变量的边缘分布
- 2.1.1 参数法计算流程
- 2.1.2 非参数法
- 2.2 选取适当的Copula函数
- 2.3 参数估计
- 3. 完整代码
- 参考资料
1. 案例描述
现有上海和深圳股市同时期日开盘价、最高价、最低价、收盘价、收益率等数据,跨度为2000年1月至2007年4月,各1696组数据。完整数据保存在文件hushi.xls和shenshi.xIs中,其中部分数据如表6-3和6-4所列。
其中,收益率=(收盘价-开盘价)/开盘价。根据收集到的1696组数据研究沪、深两市日收益率之间的关系,构建二元Copula模型,描述沪、深两市日收益率的相关结构。
2.实现流程
根据Copula理论,可以按照以下步骤构建Copula模型:
➢确定随机变量的边缘分布;
➢选取适当的,能够描述随机变量间相关结构的Copula函数;.
➢估计Copula模型中的未知参数。
2.1 确定随机变量的边缘分布
令X,Y分别表示沪、深两市的日收益率。首先，确定随机变量X,Y的分布。
确定随机变最分布的方法有两种：参数法和非参数法。
参数法：假定随机变量服从某种含有参数的分布,例如正态分布、t分布等常见分布，然后根据样本观测值估计分布中的参数,最后作出检验。
非参数法：基于经验分布和核光滑方法(核密度估计)把样本的经验分布函数作为总体随机变最的分布的近似，也可以根据样本观测数据,利用核密度估计的方法确定总体的分布。
2.1.1 参数法计算流程
为了确定随机变量X,Y的分布类型,首先作出它们的频率直方图。
%--------------------------------------------------------------------------
clear
clc
%******************************读取数据*************************************
% 从文件hushi.xls中读取数据
hushi = xlsread('hushi.xls');
% 提取矩阵hushi的第5列数据，即沪市的日收益率数据
X = hushi(:,5);
% 从文件shenshi.xls中读取数据
shenshi = xlsread('shenshi.xls');
% 提取矩阵shenshi的第5列数据，即深市的日收益率数据
Y = shenshi(:,5);
%****************************绘制频率直方图*********************************
% 调用ecdf函数和ecdfhist函数绘制沪、深两市日收益率的频率直方图
[fx, xc] = ecdf(X);
figure;
ecdfhist(fx, xc, 30);
xlabel('沪市日收益率'); % 为X轴加标签
ylabel('f(x)'); % 为Y轴加标签
[fy, yc] = ecdf(Y);
figure;
ecdfhist(fy, yc, 30);
xlabel('深市日收益率'); % 为X轴加标签
ylabel('f(y)'); % 为Y轴加标签
图1 频率直方图
%****************************计算偏度和峰度*********************************
% 计算X和Y的偏度
xs = skewness(X)
ys = skewness(Y)
% 计算X和Y的峰度
kx = kurtosis(X)
ky = kurtosis(Y)
结果如下：
xs = -0.025299871279193
ys = -0.003554260594235
kx = 6.377426684827699
ky = 6.633941978467045
图1以及X,Y的偏度、峰度反映出:随机变量X,Y的分布均是比较对称的，并且呈现出尖峰厚尾(或重尾)的特点。正态分布是轻尾(或薄尾)分布,可以初步断定X,Y不服从正态分布。
下面调用jbtest、kstest和llietest函数分别对X和Y进行正态性检验。
%******************************正态性检验***********************************
% 分别调用jbtest、kstest和lillietest函数对X进行正态性检验
[h_X_jb,p_X_jb] = jbtest(X) % Jarque-Bera检验
[h_X_ks,p_X_ks] = kstest(X,[X,normcdf(X,mean(X),std(X))]) % Kolmogorov-Smirnov检验
[h_X_lillie, p_X_lillie] = lillietest(X) % Lilliefors检验
% 分别调用jbtest、kstest和lillietest函数对Y进行正态性检验
[h_Y_jb,p_Y_jb] = jbtest(Y) % Jarque-Bera检验
[h_Y_ks,p_Y_ks] = kstest(Y,[Y,normcdf(Y,mean(Y),std(Y))]) % Kolmogorov-Smirnov检验
[h_Y_lillie, p_Y_lillie] = lillietest(Y) % Lilliefors检验
计算结果如下：
h_X_jb = 1
p_X_jb = 1.000000000000000e-03
h_X_ks = logical 1
p_X_ks = 9.280192003436800e-07
h_X_lillie = 1
p_X_lillie = 1.000000000000000e-03
h_Y_jb = 1
p_Y_jb = 1.000000000000000e-03
h_Y_ks = logical 1
p_Y_ks = 4.846695681424169e-06
h_Y_lillie = 1
p_Y_lillie = 1.000000000000000e-03
以上三种检验的h值均为1，p值均小于0.01,说明X和Y都不服从正态分布,而是服从某种对称的尖峰厚尾分布。遗憾的是,常见分布中难以找到这种类型的分布。下面利用非参数法确定X,Y的分布。
2.1.2 非参数法
当总体的分布不好确定时,可以调用ecdf函数求样本经验分布函数,作为总体分布函数的近似,或调用ksdensity函数,用核光滑方法估计总体的分布。
- 调用ecdf函数求样本经验分布函数
ecdf函数返回的向量Xsort和Ysort是各自经过排序后的样本数据，fx和fy是分别与向量Xsort和Ysort对应的经验分布函数值向量。为了求原始样本(未经排序的样本)观测值所对应的经验分布函数值，上面用到了样条插值法。当然也可以不用样条插值法，利用反排序也行,命令如下：
%****************************求经验分布函数值*******************************
% 调用ecdf函数求X和Y的经验分布函数
[fx, Xsort] = ecdf(X);
[fy, Ysort] = ecdf(Y);
% 调用spline函数，利用样条插值法求原始样本点处的经验分布函数值
U1 = spline(Xsort(2:end),fx(2:end),X);
V1 = spline(Ysort(2:end),fy(2:end),Y);
% 调用ecdf函数求X和Y的经验分布函数
[fx, Xsort] = ecdf(X);
[fy, Ysort] = ecdf(Y);
% 提取fx和fy的第2个至最后一个元素，即排序后样本点处的经验分布函数值
fx = fx(2:end);
fy = fy(2:end);
% 通过排序和反排序恢复原始样本点处的经验分布函数值U1和V1
[Xsort,id] = sort(X);
[idsort,id] = sort(id);
U1 = fx(id);
[Ysort,id] = sort(Y);
[idsort,id] = sort(id);
V1 = fy(id);
两次得到的U1完全一样,V1也完全一样。
- 调用ksdensity函数进行总体分布的估计
%*******************************核分布估计**********************************
% 调用ksdensity函数分别计算原始样本X和Y处的核分布估计值
U2 = ksdensity(X,X,'function','cdf');
V2 = ksdensity(Y,Y,'function','cdf');
调用ecdf函数得到的U1和调用ksdensity函数得到的U2不完全相同,V1和V2也不完全相同,但是U1和U2、V1和V2的差别都非常微小，如图2所示,经验分布函数图和核分布估计图几乎重合。
% **********************绘制经验分布函数图和核分布估计图**********************
[Xsort,id] = sort(X); % 为了作图的需要，对X进行排序
figure(2); % 新建一个图形窗口
subplot(1,2,1)
plot(Xsort,U1(id),'c','LineWidth',5); % 绘制沪市日收益率的经验分布函数图
hold on
plot(Xsort,U2(id),'k-.','LineWidth',2); % 绘制沪市日收益率的核分布估计图
legend('经验分布函数','核分布估计', 'Location','NorthWest'); % 加标注框
xlabel('沪市日收益率'); % 为X轴加标签
ylabel('F(x)'); % 为Y轴加标签
[Ysort,id] = sort(Y); % 为了作图的需要，对Y进行排序
subplot(1,2,2)
plot(Ysort,V1(id),'c','LineWidth',5); % 绘制深市日收益率的经验分布函数图
hold on
plot(Ysort,V2(id),'k-.','LineWidth',2); % 绘制深市日收益率的核分布估计图
legend('经验分布函数','核分布估计', 'Location','NorthWest'); % 加标注框
xlabel('深市日收益率'); % 为X轴加标签
ylabel('F(x)'); % 为Y轴加标签
图2 经验分布函数图
2.2 选取适当的Copula函数
在确定X的边缘分布U= F(x)和Y的边缘分布V=G(x)之后,就可以根据(U_i,V_i,)(i=1,2…n)的二元直方图的形状选取适当的Copula函数。绘制二元频数直方图的命令如下:
%****************************绘制二元频数直方图*****************************
% 调用ksdensity函数分别计算原始样本X和Y处的核分布估计值
U = ksdensity(X,X,'function','cdf');
V = ksdensity(Y,Y,'function','cdf');
figure; % 新建一个图形窗口
% 绘制边缘分布的二元频数直方图，
hist3([U(:) V(:)],[30,30])
xlabel('U（沪市）'); % 为X轴加标签
ylabel('V（深市）'); % 为Y轴加标签
zlabel('频数'); % 为z轴加标签
图3 二元频数直方图
以上命令作出的频数直方图如图3所示,在频数直方图的基础上还可以绘制频率直方图,并且频率直方图可以作为(U,V)的联合密度函数(即Copula密度函数)的估计。由频数直方图绘制频率直方图的命令如下:
%****************************绘制二元频率直方图*****************************
figure
hist3([U(:) V(:)],[30,30])
h = get(gca, 'Children');
cuv = get(h, 'ZData');
set(h, 'ZData', cuv * 30 * 30 / length(X));
xlabel('U（沪市）'); % 为X轴加标签
ylabel('V（深市）'); % 为Y轴加标签
zlabel('频数'); % 为z轴加标签
图4 二元频率直方图
作出的二元频率直方图如图4所示,可以看出,频率直方图具有基本对称的尾部，也就是说(U,V)的联合密度函数(即Copula密度函数)具有对称的尾部,因此可以选取二元正态Copula函数或二元t-Copula函数来描述原始数据的相关结构。
2.3 参数估计
考虑一般情况，边缘分布中可能含有未知参数，并且选取的Copula函数中也含有未知参数，因此需要进行参数估计。常用的参数估计方法有最大似然估计、分步估计和半参数估计。
对于选取的二元正态Copula和二元t-Copula,用核分布估计求随机变量X,Y的边缘分布,然后调用copulafit 函数估计Copula中的参数,命令如下:
%***********************求Copula中参数的估计值******************************
% 调用copulafit函数估计二元正态Copula中的线性相关参数
rho_norm = copulafit('Gaussian',[U(:), V(:)])
% 调用copulafit函数估计二元t-Copula中的线性相关参数和自由度
[rho_t,nuhat,nuci] = copulafit('t',[U(:), V(:)])
计算结果如下：
rho_norm =
1.000000000000000 0.926423396181286
0.926423396181286 1.000000000000000
rho_t =
1.000000000000000 0.932538838072873
0.932538838072873 1.000000000000000
nuhat =
4.008923194828739
nuci =
2.983921553182909 5.033924836474569
估计出Copula中的参数之后，可以调用copulapdf函数和copulacdf函数分别计算Copula密度函数和分布函数值，然后绘制Copula密度函数和分布函数图，相应的MATIAB命令如下：
%********************绘制Copula的密度函数和分布函数图************************
[Udata,Vdata] = meshgrid(linspace(0,1,31)); % 为绘图需要，产生新的网格数据
% 调用copulapdf函数计算网格点上的二元正态Copula密度函数值
Cpdf_norm = copulapdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);
% 调用copulacdf函数计算网格点上的二元正态Copula分布函数值
Ccdf_norm = copulacdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);
% 调用copulapdf函数计算网格点上的二元t-Copula密度函数值
Cpdf_t = copulapdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);
% 调用copulacdf函数计算网格点上的二元t-Copula分布函数值
Ccdf_t = copulacdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);
% 绘制二元正态Copula的密度函数和分布函数图
figure(5); % 新建图形窗口
subplot(1,2,1)
surf(Udata,Vdata,reshape(Cpdf_norm,size(Udata))); % 绘制二元正态Copula密度函数图
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('c(u,v)'); % 为z轴加标签
title('二元正态Copula的密度函数')
subplot(1,2,2)
surf(Udata,Vdata,reshape(Ccdf_norm,size(Udata))); % 绘制二元正态Copula分布函数图
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('C(u,v)'); % 为z轴加标签
title('二元正态Copula的分布函数')
% 绘制二元t-Copula的密度函数和分布函数图
figure(6); % 新建图形窗口
subplot(1,2,1)
surf(Udata,Vdata,reshape(Cpdf_t,size(Udata))); % 绘制二元t-Copula密度函数图
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('c(u,v)'); % 为z轴加标签
title('二元t-Copula的密度函数')
subplot(1,2,2)
surf(Udata,Vdata,reshape(Ccdf_t,size(Udata))); % 绘制二元t-Copula分布函数图
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('C(u,v)'); % 为z轴加标签
title('二元t-Copula的分布函数')
图5 二元正态Copula密度函数和分布函数
图6 二元t-Copula密度函数和分布函数
以上命令绘制的图形如图5和图6所示,可以看出，与线性相关参数为p=0.9264的二元正态Copula相比,线性相关参数为p=0.9325,自由度为k=4的二元t-Copula的密度函数具有更厚的尾部,更能反映变量之间的尾部相关性。
从图4可以看出沪、深两市日收益率之间有较强的尾部相关性，再将图4图3和图6密度函数图加以对比，可知线性相关参数为ρ=0.9325，自由度为k=4的二元t-Copula较好地反映了沪、深两市日收益率之间的尾部相关性，计算出的尾部相关系数为：0.6934，计算过程如下：
k = 4
pho = 0.932538838072873
x = sqrt(k+1)*sqrt(1-pho )/sqrt(1+pho )
lammda = 2 - tcdf(x, k+1)*2
估计出Copula中的参数之后,还可以调用copulastat 函数求Kendall 秩相关系数、Spearman秩相关系数的估计。
%**************求Kendall秩相关系数和Spearman秩相关系数***********************
% 调用copulastat函数求二元正态Copula对应的Kendall秩相关系数
Kendall_norm = copulastat('Gaussian',rho_norm)
% 调用copulastat函数求二元正态Copula对应的Spearman秩相关系数
Spearman_norm = copulastat('Gaussian',rho_norm,'type','Spearman')
% 调用copulastat函数求二元t-Copula对应的Kendall秩相关系数
Kendall_t = copulastat('t',rho_t)
% 调用copulastat函数求二元t-Copula对应的Spearman秩相关系数
% Spearman_t = copulastat('t',rho_t,nuhat,'type','Spearman')
Spearman_t = copulastat('t',rho_t,nuhat,'type','Spearman') % Yue
% 直接根据沪、深两市日收益率的原始观测数据，调用corr函数求Kendall秩相关系数
Kendall = corr([X,Y],'type','Kendall')
% 直接根据沪、深两市日收益率的原始观测数据，调用corr函数求Spearman秩相关系数
Spearman = corr([X,Y],'type','Spearman')
计算结果如下：
Kendall_norm =
1.000000000000000 0.754266435119928
0.754266435119928 1.000000000000000
Spearman_norm =
1.000000000000000 0.919818249600397
0.919818249600397 1.000000000000000
Kendall_t =
1.000000000000000 0.764823295419639
0.764823295419639 1.000000000000000
Spearman_t =
1.000000000000000 0.917923638116616
0.917923638116616 1.000000000000000
Kendall =
1.000000000000000 0.757242444481550
0.757242444481550 1.000000000000000
Spearman =
1.000000000000000 0.912625848233055
0.912625848233055 1.000000000000000
将以上求出的Kendall秩相关系数Kendall_ norm. Kendall t和Kendall加以对比，将求出的Spearman秩相关系数Spearman_ norm 、Spearman_ t和Spearman加以对比。可以看出，Kendall norm更接近Kendall, Spearman_norm更接近Spearman，说明了线性相关参数为p=0.9264的二元正态Copula较好地反映了沪.深两市日收益率之间的秩相关性。
对于沪深两市日收益率的观测数据，我们构建了二元正态Copula模型和二元t-Copula模型，为了评价两个模型的优劣，下面引入经验Copula的概念。
%******************************模型评价*************************************
% 调用ecdf函数求X和Y的经验分布函数
[fx, Xsort] = ecdf(X);
[fy, Ysort] = ecdf(Y);
% 调用spline函数，利用样条插值法求原始样本点处的经验分布函数值
U = spline(Xsort(2:end),fx(2:end),X);
V = spline(Ysort(2:end),fy(2:end),Y);
% 定义经验Copula函数C(u,v)
C = @(u,v)mean((U <= u).*(V <= v));
% 为作图的需要，产生新的网格数据
[Udata,Vdata] = meshgrid(linspace(0,1,31));
% 通过循环计算经验Copula函数在新产生的网格点处的函数值
for i=1:numel(Udata)
CopulaEmpirical(i) = C(Udata(i),Vdata(i));
end
figure(7); % 新建图形窗口
% 绘制经验Copula分布函数图像
surf(Udata,Vdata,reshape(CopulaEmpirical,size(Udata)))
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('Empirical Copula C(u,v)'); % 为z轴加标签
% 通过循环计算经验Copula函数在原始样本点处的函数值
CUV = zeros(size(U(:)));
for i=1:numel(U)
CUV(i) = C(U(i),V(i));
end
% 计算线性相关参数为0.9264的二元正态Copula函数在原始样本点处的函数值
rho_norm = 0.9264;
Cgau = copulacdf('Gaussian',[U(:), V(:)],rho_norm);
% 计算线性相关参数为0.9325，自由度为4的二元t-Copula函数在原始样本点处的函数值
rho_t = 0.9325;
k = 4.0089;
Ct = copulacdf('t',[U(:), V(:)],rho_t,k);
% 计算平方欧氏距离
dgau2 = (CUV-Cgau)'*(CUV-Cgau)
dt2 = (CUV-Ct)'*(CUV-Ct)
计算结果如下：
dgau2 = 0.018623588951791
dt2 = 0.014494967967151
图7 经验Copula分布函数图
以上命令绘制出的经验Copula分布函数图如图7所示。由计算出的平方欧氏距离可知，线性相关参数为0.9264的二元正态Copula与经验Copula的平方欧氏距离d=0. 0186；线性相关参数为0. 9325，自由度为4的二元t- Copula与经验Copula的平方欧氏距离d =0.014 5。因此在平方欧氏距离标准下，可以认为二元t - Copula模型能更好地拟合沪、深两市的日收益率观测数据。
3. 完整代码
%--------------------------------------------------------------------------
% Copula理论及应用实例
%--------------------------------------------------------------------------
clear
clc
%******************************读取数据*************************************
% 从文件hushi.xls中读取数据
hushi = xlsread('hushi.xls');
% 提取矩阵hushi的第5列数据，即沪市的日收益率数据
X = hushi(:,5);
% 从文件shenshi.xls中读取数据
shenshi = xlsread('shenshi.xls');
% 提取矩阵shenshi的第5列数据，即深市的日收益率数据
Y = shenshi(:,5);
%****************************绘制频率直方图*********************************
% 调用ecdf函数和ecdfhist函数绘制沪、深两市日收益率的频率直方图
[fx, xc] = ecdf(X);
figure;
subplot(1,2,1)
ecdfhist(fx, xc, 30);
xlabel('沪市日收益率'); % 为X轴加标签
ylabel('f(x)'); % 为Y轴加标签
[fy, yc] = ecdf(Y);
subplot(1,2,2)
ecdfhist(fy, yc, 30);
xlabel('深市日收益率'); % 为X轴加标签
ylabel('f(y)'); % 为Y轴加标签
%****************************计算偏度和峰度*********************************
% 计算X和Y的偏度
xs = skewness(X)
ys = skewness(Y)
% 计算X和Y的峰度
kx = kurtosis(X)
ky = kurtosis(Y)
%******************************正态性检验***********************************
% 分别调用jbtest、kstest和lillietest函数对X进行正态性检验
[h_X_jb,p_X_jb] = jbtest(X) % Jarque-Bera检验
[h_X_ks,p_X_ks] = kstest(X,[X,normcdf(X,mean(X),std(X))]) % Kolmogorov-Smirnov检验
[h_X_lillie, p_X_lillie] = lillietest(X) % Lilliefors检验
% 分别调用jbtest、kstest和lillietest函数对Y进行正态性检验
[h_Y_jb,p_Y_jb] = jbtest(Y) % Jarque-Bera检验
[h_Y_ks,p_Y_ks] = kstest(Y,[Y,normcdf(Y,mean(Y),std(Y))]) % Kolmogorov-Smirnov检验
[h_Y_lillie, p_Y_lillie] = lillietest(Y) % Lilliefors检验
%****************************求经验分布函数值*******************************
% 调用ecdf函数求X和Y的经验分布函数
[fx, Xsort] = ecdf(X);
[fy, Ysort] = ecdf(Y);
% 调用spline函数，利用样条插值法求原始样本点处的经验分布函数值
U1 = spline(Xsort(2:end),fx(2:end),X);
V1 = spline(Ysort(2:end),fy(2:end),Y);
% 调用ecdf函数求X和Y的经验分布函数
[fx, Xsort] = ecdf(X);
[fy, Ysort] = ecdf(Y);
% 提取fx和fy的第2个至最后一个元素，即排序后样本点处的经验分布函数值
fx = fx(2:end);
fy = fy(2:end);
% 通过排序和反排序恢复原始样本点处的经验分布函数值U1和V1
[Xsort,id] = sort(X);
[idsort,id] = sort(id);
U1 = fx(id);
[Ysort,id] = sort(Y);
[idsort,id] = sort(id);
V1 = fy(id);
%*******************************核分布估计**********************************
% 调用ksdensity函数分别计算原始样本X和Y处的核分布估计值
U2 = ksdensity(X,X,'function','cdf');
V2 = ksdensity(Y,Y,'function','cdf');
% **********************绘制经验分布函数图和核分布估计图**********************
[Xsort,id] = sort(X); % 为了作图的需要，对X进行排序
figure(2); % 新建一个图形窗口
subplot(1,2,1)
plot(Xsort,U1(id),'c','LineWidth',5); % 绘制沪市日收益率的经验分布函数图
hold on
plot(Xsort,U2(id),'k-.','LineWidth',2); % 绘制沪市日收益率的核分布估计图
legend('经验分布函数','核分布估计', 'Location','NorthWest'); % 加标注框
xlabel('沪市日收益率'); % 为X轴加标签
ylabel('F(x)'); % 为Y轴加标签
[Ysort,id] = sort(Y); % 为了作图的需要，对Y进行排序
subplot(1,2,2)
plot(Ysort,V1(id),'c','LineWidth',5); % 绘制深市日收益率的经验分布函数图
hold on
plot(Ysort,V2(id),'k-.','LineWidth',2); % 绘制深市日收益率的核分布估计图
legend('经验分布函数','核分布估计', 'Location','NorthWest'); % 加标注框
xlabel('深市日收益率'); % 为X轴加标签
ylabel('F(x)'); % 为Y轴加标签
%****************************绘制二元频数直方图*****************************
% 调用ksdensity函数分别计算原始样本X和Y处的核分布估计值
U = ksdensity(X,X,'function','cdf');
V = ksdensity(Y,Y,'function','cdf');
figure(3); % 新建一个图形窗口
% 绘制边缘分布的二元频数直方图，
hist3([U(:) V(:)],[30,30])
xlabel('U（沪市）'); % 为X轴加标签
ylabel('V（深市）'); % 为Y轴加标签
zlabel('频数'); % 为z轴加标签
%****************************绘制二元频率直方图*****************************
figure(4); % 新建一个图形窗口
% 绘制边缘分布的二元频数直方图，
hist3([U(:) V(:)],[30,30])
h = get(gca, 'Children'); % 获取频数直方图的句柄值
cuv = get(h, 'ZData'); % 获取频数直方图的Z轴坐标
set(h,'ZData',cuv*30*30/length(X)); % 对频数直方图的Z轴坐标作变换
xlabel('U（沪市）'); % 为X轴加标签
ylabel('V（深市）'); % 为Y轴加标签
zlabel('c(u,v)'); % 为z轴加标签
%***********************求Copula中参数的估计值******************************
% 调用copulafit函数估计二元正态Copula中的线性相关参数
rho_norm = copulafit('Gaussian',[U(:), V(:)])
% 调用copulafit函数估计二元t-Copula中的线性相关参数和自由度
[rho_t,nuhat,nuci] = copulafit('t',[U(:), V(:)])
%********************绘制Copula的密度函数和分布函数图************************
[Udata,Vdata] = meshgrid(linspace(0,1,31)); % 为绘图需要，产生新的网格数据
% 调用copulapdf函数计算网格点上的二元正态Copula密度函数值
Cpdf_norm = copulapdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);
% 调用copulacdf函数计算网格点上的二元正态Copula分布函数值
Ccdf_norm = copulacdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);
% 调用copulapdf函数计算网格点上的二元t-Copula密度函数值
Cpdf_t = copulapdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);
% 调用copulacdf函数计算网格点上的二元t-Copula分布函数值
Ccdf_t = copulacdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);
% 绘制二元正态Copula的密度函数和分布函数图
figure(5); % 新建图形窗口
subplot(1,2,1)
surf(Udata,Vdata,reshape(Cpdf_norm,size(Udata))); % 绘制二元正态Copula密度函数图
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('c(u,v)'); % 为z轴加标签
title('二元正态Copula的密度函数')
subplot(1,2,2)
surf(Udata,Vdata,reshape(Ccdf_norm,size(Udata))); % 绘制二元正态Copula分布函数图
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('C(u,v)'); % 为z轴加标签
title('二元正态Copula的分布函数')
% 绘制二元t-Copula的密度函数和分布函数图
figure(6); % 新建图形窗口
subplot(1,2,1)
surf(Udata,Vdata,reshape(Cpdf_t,size(Udata))); % 绘制二元t-Copula密度函数图
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('c(u,v)'); % 为z轴加标签
title('二元t-Copula的密度函数')
subplot(1,2,2)
surf(Udata,Vdata,reshape(Ccdf_t,size(Udata))); % 绘制二元t-Copula分布函数图
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('C(u,v)'); % 为z轴加标签
title('二元t-Copula的分布函数')
%**************求Kendall秩相关系数和Spearman秩相关系数***********************
% 调用copulastat函数求二元正态Copula对应的Kendall秩相关系数
Kendall_norm = copulastat('Gaussian',rho_norm)
% 调用copulastat函数求二元正态Copula对应的Spearman秩相关系数
Spearman_norm = copulastat('Gaussian',rho_norm,'type','Spearman')
% 调用copulastat函数求二元t-Copula对应的Kendall秩相关系数
Kendall_t = copulastat('t',rho_t)
% 调用copulastat函数求二元t-Copula对应的Spearman秩相关系数
% Spearman_t = copulastat('t',rho_t,nuhat,'type','Spearman')
Spearman_t = copulastat('t',rho_t,nuhat,'type','Spearman') % Yue
% 直接根据沪、深两市日收益率的原始观测数据，调用corr函数求Kendall秩相关系数
Kendall = corr([X,Y],'type','Kendall')
% 直接根据沪、深两市日收益率的原始观测数据，调用corr函数求Spearman秩相关系数
Spearman = corr([X,Y],'type','Spearman')
%******************************模型评价*************************************
% 调用ecdf函数求X和Y的经验分布函数
[fx, Xsort] = ecdf(X);
[fy, Ysort] = ecdf(Y);
% 调用spline函数，利用样条插值法求原始样本点处的经验分布函数值
U = spline(Xsort(2:end),fx(2:end),X);
V = spline(Ysort(2:end),fy(2:end),Y);
% 定义经验Copula函数C(u,v)
C = @(u,v)mean((U <= u).*(V <= v));
% 为作图的需要，产生新的网格数据
[Udata,Vdata] = meshgrid(linspace(0,1,31));
% 通过循环计算经验Copula函数在新产生的网格点处的函数值
for i=1:numel(Udata)
CopulaEmpirical(i) = C(Udata(i),Vdata(i));
end
figure(7); % 新建图形窗口
% 绘制经验Copula分布函数图像
surf(Udata,Vdata,reshape(CopulaEmpirical,size(Udata)))
xlabel('U'); % 为X轴加标签
ylabel('V'); % 为Y轴加标签
zlabel('Empirical Copula C(u,v)'); % 为z轴加标签
% 通过循环计算经验Copula函数在原始样本点处的函数值
CUV = zeros(size(U(:)));
for i=1:numel(U)
CUV(i) = C(U(i),V(i));
end
% 计算线性相关参数为0.9264的二元正态Copula函数在原始样本点处的函数值
rho_norm = 0.9264;
Cgau = copulacdf('Gaussian',[U(:), V(:)],rho_norm);
% 计算线性相关参数为0.9325，自由度为4的二元t-Copula函数在原始样本点处的函数值
rho_t = 0.9325;
k = 4.0089;
Ct = copulacdf('t',[U(:), V(:)],rho_t,k);
% 计算平方欧氏距离
dgau2 = (CUV-Cgau)'*(CUV-Cgau)
dt2 = (CUV-Ct)'*(CUV-Ct)
参考资料
《MATLAB统计分析与应用：40个案例分析》
