pair-copula.rar_pair copula程序_pair copula算法_pair-copula_paircopu

%******************************读取数据************************************* %从文件hushi.xls中读取数据 hushi = xlsread('五台山各月风电.xlsx');%风数据文件 % 提取矩阵hushi的第1列数据 X = hushi(:,2); % 从文件shenshi.xls中读取数据 shenshi = xlsread('汉川各月光伏.xlsx');%光数据文件 % 提取矩阵shenshi的第1列数据 Y = shenshi(:,2); %****************************求经验分布函数值******************************* % 调用ecdf函数求X和Y的经验分布函数 [fx, Xsort] = ecdf(X); [fy, Ysort] = ecdf(Y); % 调用spline函数，利用样条插值法求原始样本点处的经验分布函数值 U1 = spline(Xsort(2:end),fx(2:end),X);%Xsort,fx求函数，在X处求插值 V1 = spline(Ysort(2:end),fy(2:end),Y); % 调用ecdf函数求X和Y的经验分布函数 [fx, Xsort] = ecdf(X);%Xsort为按顺序由小到大排列的数，fx为其出现频率 [fy, Ysort] = ecdf(Y); % 提取fx和fy的第2个至最后一个元素，即排序后样本点处的经验分布函数值 fx = fx(2:end); fy = fy(2:end); %*******************************核分布估计********************************** % 调用ksdensity函数分别计算原始样本X和Y处的核分布估计值 U2 = ksdensity(X,X,'function','cdf','width',0.002); V2 = ksdensity(Y,Y,'function','cdf','width',0.01); % **********************绘制经验分布函数图和核分布估计图********************** [Xsort,id] = sort(X); % 为了作图的需要，对X进行排序 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轴加标签 %****************************绘制二元频数直方图***************************** %%利用核密度估计 figure; % 新建一个图形窗口 % 绘制边缘分布的二元频数直方图， hist3([U2(:) V2(:)],[30,30]) xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('频数'); % 为z轴加标签 %****************************绘制二元频率直方图***************************** figure; % 新建一个图形窗口 % 绘制边缘分布的二元频数直方图， hist3([U2(:) V2(:)],[30,30]) h = get(gca, 'Children'); % 获取频数直方图的句柄值 cuv = get(h, 'ZData'); % 获取频数直方图的Z轴坐标 set(h,'ZData',cuv*30*30/length(X)); % 对频数直方图的Z轴坐标作变换 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('频数100%'); % 为z轴加标签 %***********************求Copula中参数的估计值****************************** %%利用核密度估计 % 调用copulafit函数估计二元正态Copula中的线性相关参数 rho_norm = copulafit('Gaussian',[U2(:), V2(:)])%返回p估计值 % 调用copulafit函数估计二元t-Copula中的线性相关参数和自由度 [rho_t,nuhat,nuci] = copulafit('t',[U2(:), V2(:)]) %调用copulafit函数估计Glumbel中的线性相关参数 paramhatG=copulafit('Gumbel',[U2(:), V2(:)]) %调用copulafit函数估计Clayto中的线性相关参数 paramhatC=copulafit('Clayton',[U2(:), V2(:)]) %调用copulafit函数估计Frank中的线性相关参数 paramhatF=copulafit('Frank',[U2(:), V2(:)]) %********************绘制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); %调用copulapdf函数计算网格点上的Gumbel-Copula密度函数值 Cpdf_Gum = copulapdf('Gumbel',[Udata(:), Vdata(:)],paramhatG); % 调用copulacdf函数计算网格点上的Gumbel-Copula分布函数值 Ccdf_Gum = copulacdf('Gumbel',[Udata(:), Vdata(:)],paramhatG); %调用copulapdf函数计算网格点上的Clayton-Copula密度函数值 Cpdf_CL= copulapdf('Clayton',[Udata(:), Vdata(:)],paramhatC); % 调用copulacdf函数计算网格点上的Clayton-Copula分布函数值 Ccdf_CL = copulacdf('Clayton',[Udata(:), Vdata(:)],paramhatC); %调用copulapdf函数计算网格点上的Frank-Copula密度函数值 Cpdf_F = copulapdf('Frank',[Udata(:), Vdata(:)],paramhatF); % 调用copulacdf函数计算网格点上的Clayton-Copula分布函数值 Ccdf_F = copulacdf('Frank',[Udata(:), Vdata(:)],paramhatF); % 绘制二元正态Copula的密度函数和分布函数图 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Cpdf_norm,size(Udata))); % 绘制二元正态Copula密度函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('高斯Copula概率密度函数值'); % 为z轴加标签 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Ccdf_norm,size(Udata))); % 绘制二元正态Copula分布函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('高斯Copula累计分布值'); % 为z轴加标签 % 绘制二元t-Copula的密度函数和分布函数图 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Cpdf_t,size(Udata))); % 绘制二元t-Copula密度函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('t-Copula概率密度函数值'); % 为z轴加标签 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Ccdf_t,size(Udata))); % 绘制二元t-Copula分布函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('t-Copula累计分布值'); % 为z轴加标签 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Cpdf_Gum,size(Udata))); % 绘制二元t-Copula密度函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('Gumbel-Copula概率密度函数值'); % 为z轴加标签 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Ccdf_Gum,size(Udata))); % 绘制二元t-Copula分布函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('Gumbel-Copula累计分布值'); % 为z轴加标签 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Cpdf_CL,size(Udata))); % 绘制二元t-Copula密度函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('Clayton-Copula概率密度函数值'); % 为z轴加标签 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Ccdf_CL,size(Udata))); % 绘制二元t-Copula分布函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('Clayton-Copula累计分布值'); % 为z轴加标签 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Cpdf_F,size(Udata))); % 绘制二元t-Copula密度函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('Frank-Copula概率密度函数值'); % 为z轴加标签 figure; % 新建图形窗口 surf(Udata,Vdata,reshape(Ccdf_F,size(Udata))); % 绘制二元t-Copula分布函数图 xlabel('风电出力概率'); % 为X轴加标签 ylabel('光伏出力概率'); % 为Y轴加标签 zlabel('Frank-Copula累计分布值'); % 为z轴加标签 %**************求Kendall秩相关系数和Spearman秩相关系数*********************** % 调用copulastat函数求二元正态Copula对应的Kendall秩相关系数 Kendall_norm = copulastat('Gaussian',rho_norm) % 调用copulastat函数求二元正态Copula对应的Spearman秩相关系数 Spearman_norm = copulastat('Gaussian',rho_norm,'type','Spearman') % 调用copulastat函数求Gumbel正态Copula对应的Kendall秩相关系数 Kendall_G = copulastat('Gumbel',paramhatG) % 调用copulastat函数求Gumbel正态Copula对应的Spearman秩相关系数 Spearman_G= copulastat('Gumbel',paramhatG,'type','Spearman') % 调用copulastat函数求Gumbel正态Copula对应的Kendall秩相关系数 Kendall_CL = copulastat('Clayton',paramhatC) % 调用copulastat函数求Gumbel正态Copula对应的Spearman秩相关系数 Spearman_CL= copulastat('Clayton',paramhatC,'type','Spearman') % 调用copu