Copula理论及MATLAB应用实例.docx

1、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, x

2、c = 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和lil

3、lietest函数对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 = lillietes

4、t(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个至最后一个元素，即排序后样本点处的经验分布函数值

5、fx = fx(2:end);fy = fy(2:end);% 通过排序和反排序恢复原始样本点处的经验分布函数值U1和V1Xsort,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,i

6、d = sort(X); % 为了作图的需要，对X进行排序figure; % 新建一个图形窗口plot(Xsort,U1(id),c,LineWidth,5); % 绘制沪市日收益率的经验分布函数图hold onplot(Xsort,U2(id),k-.,LineWidth,2); % 绘制沪市日收益率的核分布估计图legend(经验分布函数,核分布估计, Location,NorthWest); % 加标注框xlabel(沪市日收益率); % 为X轴加标签ylabel(F(x); % 为Y轴加标签Ysort,id = sort(Y); % 为了作图的需要，对Y进行排序figure; % 新建

7、一个图形窗口plot(Ysort,V1(id),c,LineWidth,5); % 绘制深市日收益率的经验分布函数图hold onplot(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 = ksdensi

8、ty(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

9、*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的密度函数和分布函数图*Uda

10、ta,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 = copula

11、pdf(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轴加标签fig

12、ure; % 新建图形窗口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); %

13、 为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对应的Spearm

14、an秩相关系数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,Kend

15、all)% 直接根据沪、深两市日收益率的原始观测数据，调用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)

16、.*(V x=0:pi/10:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,x,y2)图5.1.1 函数plot 绘制的正弦曲线在绘制曲线图形时，常常采用多种颜色或线型来区分不同的数据组，MATLAB 软件专门提供了这方面的参数选项（见表5.1.1），我们只要在每个坐标后加上相关字符串，就可实现它们的功能。表5.1.1 绘图参数表色彩字符颜色线型字符线型格式标记符号数据点形式标记符号数据点形式y 黄- 实线. 点 大于号例如，在上例中输入 plot(x,y1,r+-,x,y2,k*:)则得图5.1.2图5.1.2 使用不同标记的plot 函数绘制的正弦曲线5.1.

17、2 图形修饰MATLAB 软件为用户提供了一些特殊的图形函数，用于修饰已经绘制好的图形。表5.1.2 图形修饰函数表函数 含义grid on (/off) 给当前图形标记添加（取消）网络xlable(string) 标记横坐标ylabel(string) 标记纵坐标title(string) 给图形添加标题text(x,y,string) 在图形的任意位置增加说明性文本信息gtext(string) 利用鼠标添加说明性文本信息axis(xmin xmax ymin ymax) 设置坐标轴的最小最大值例5.1.2 给例5.1.1 的图形中加入网络和标记。（见图5.1.3 和5.1.4） x=0:

18、pi/10:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,x,y2) grid on xlabel(independent variable X) ylabel(Dependent Variable Y1 & Y2) title(Sine and Cosine Curve) text(1.5,0.3,cos(x) gtext(sin(x) axis(0 2*pi -0.9 0.9)图5.1.3 使用了图形修饰的plot 函数绘制的正弦曲线5.1.3 图形的比较显示在一般默认的情况下，MATLAB 每次使用plot 函数进行图形绘制，将重新产生一个图形窗口。但有时希

19、望后续的图形能够和前面所绘制的图形进行比较。一般来说有两种方法：一是采用hold on(/off)命令，将新产生的图形曲线叠加到已有的图形上；二是采用subplot(m,n,k)函数，将图形窗口分隔成n m 个子图，并选择第k 个子图作为当前图形，然后在同一个视图窗口中画出多个小图形。例5.1.3 在同一窗口中绘制线段。 x=0:pi/10:2*pi; y1=sin(x); y2=cos(x); y3=x; y4=log(x); plot(x,y1,x,y2) hold on plot(x,y3) plot(x,y4) hold off例5.1.4 在多个窗口中绘制图形。（见图5.1.6） x

20、=0:pi/10:2*pi; y1=sin(x); y2=cos(x); y3=exp(x); y4=log(x); subplot(2,2,1); plot(x,y1); subplot(2,2,2); plot(x,y2); subplot(2,2,3); plot(x,y3); subplot(2,2,4); plot(x,y4);说明 （1）子窗口的序号按行由上往下，按列从左向右编号。（2）如果不用指令clf 清除，以后图形将被绘制在子图形窗口中。图5.1.6 图形的比较显示（图形窗口分割方法）5.2 三维立体图形5.2.1 三维曲线图与二维图形相对应，MATLAB 提供了plot3 函数，可以在三维空间中绘制三维曲线，它的格式类似于plot，不过多了z 方向的数据。plot3 的调用格式为：plot3(x1,y1,z1,x2,y2,z2,.)其中x1,y1,z1,x2,y2,z2,等分别为维数相同的向量，分别存储着曲线的三个坐标值，该函数的