Copula理论MATLAB应用实例

.资料%--------------------------------------------------------------------------%Copula理论及其在matlab中的实现程序应用实例%--------------------------------------------------------------------------%******************************读取数据*************************************%从文件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轴加标签%****************************计算偏度和峰度*********************************%计算X和Y的偏度xs=skewness(X)ys=skewness(Y)%计算X和Y的峰度kx=kurtosis(X)ky=kurtosis(Y)%******************************正态性检验***********************************%分别调用jbtest、kstest和lillietest函数对X进行正态性检验[h,p]=jbtest(X)%Jarque-Bera检验[h,p]=kstest(X,[X,normcdf(X,mean(X),std(X))])%Kolmogorov-Smirnov检验.资料[h,p]=lillietest(X)%Lilliefors检验%分别调用jbtest、kstest和lillietest函数对Y进行正态性检验[h,p]=jbtest(Y)%Jarque-Bera检验[h,p]=kstest(Y,[Y,normcdf(Y,mean(Y),std(Y))])%Kolmogorov-Smirnov检验[h,p]=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;%新建一个图形窗口plot(Xsort,U1(id),'c','LineWidth',5);%绘制沪市日收益率的经验分布函数图holdonplot(Xsort,U2(id),'k-.','LineWidth',2);%绘制沪市日收益率的核分布估计图.资料legend('经验分布函数','核分布估计','Location','NorthWest');%加标注框xlabel('沪市日收益率');%为X轴加标签ylabel('F(x)');%为Y轴加标签[Ysort,id]=sort(Y);%为了作图的需要，对Y进行排序figure;%新建一个图形窗口plot(Ysort,V1(id),'c','LineWidth',5);%绘制深市日收益率的经验分布函数图holdonplot(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;%新建一个图形窗口%绘制边缘分布的二元频数直方图，hist3([U(:)V(:)],[30,30])xlabel('U（沪市）');%为X轴加标签ylabel('V（深市）');%为Y轴加标签zlabel('频数');%为z轴加标签%****************************绘制二元频率直方图*****************************figure;%新建一个图形窗口%绘制边缘分布的二元频数直方图，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;%新建图形窗口surf(Udata,Vdata,reshape(Cpdf_norm,size(Udata)));%绘制二元正态Copula密度函数图xlabel('U');%为X轴加标签ylabel('V');%为Y轴加标签zlabel('c(u,v)');%为z轴加标签figure;%新建图形窗口surf(Udata,Vdata,reshape(Ccdf_norm,size(Udata)));%绘制二元正态Copula分布函数图xlabel('U');%为X轴加标签ylabel('V');%为Y轴加标签zlabel('C(u,v)');%为z轴加标签%绘制二元t-Copula的密度函数和分布函数图figure;%新建图形窗口surf(Udata,Vdata,reshape(Cpdf_t,size(Udata)));%绘制二元t-Copula密度函数图xlabel('U');%为X轴加标签ylabel('V');%为Y轴加标签zlabel('c(u,v)');%为z轴加标签figure;%新建图形窗口surf(Udata,Vdata,reshape(Ccdf_t,size(Udata)));%绘制二元t-Copula分布函数图xlabel('U');%为X轴加标签ylabel('V');%为Y轴加标签zlabel('C(u,v)');%为z轴加标签%**************求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,'type','Spearman')%直接根据沪、深两市日收益率的原始观测数据，调用corr函数求Kendall秩相关系数Kendall=corr([X,Y],'type','Kendall')%直接根据沪、深两市日收益率的原始观测数据，调用corr函数求Spearman秩相关系数Sp