Copula理论及MATLAB应用实例.doc

%-% copula理论及应用实例%-%*读取数据*% 从文件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和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,id = 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; % 新建一个图形窗口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 = 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秩相关系数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);endfigure; % 新建图形窗口% 绘制经验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的二元