Copula函数实战(Matlab)

一、介绍

采用两个椭圆族Copula函数(t-copula,Gaussian-copula)和三个阿基米德族Copula函数(Clayton-copula，Frank-copula和Gumbel-copula)拟合了两个随机变量，两个随机变量分别为玉米和大豆的产量数据。

二、变量和数据

2012-2020年玉米和大豆的产量数据，数据如下，玉米(spot),大豆(product)。

年份	玉米	大豆
2012	2393.88	5382.88
2013	2422.44	5683.66
2014	2468.71	5857.56
2015	2493.87	5814.57
2016	2495.43	5993.5
2017	2495.01	6316.27
2018	2519.5	6303.05
2019	2560.48	6706.5
2020	2545.16	6653.58

三、程序

%******************************读取数据*************************************
% 从文件lnspot.xls中读取数据
X= xlsread('lnspot.xls');
% 从文件product.xls中读取数据
Y = xlsread('product.xls');
Y=Y(:,2);%****************************绘制频率直方图*********************************
% 调用ecdf函数和ecdfhist函数绘制频率直方图
[fx, xc] = ecdf(X);
figure;
ecdfhist(fx, xc, 30);
xlabel('lnspot');  % 为X轴加标签
ylabel('f(x)');  % 为Y轴加标签
[fy, yc] = ecdf(Y);
figure;
ecdfhist(fy, yc, 30);
xlabel('product');  % 为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);
hold on
plot(Xsort,U2(id),'k-.','LineWidth',2);
legend('经验分布函数','核分布估计', 'Location','NorthWest'); % 加标注框
xlabel('lnspot');  % 为X轴加标签
ylabel('F(x)');  % 为Y轴加标签[Ysort,id] = sort(Y);  % 为了作图的需要，对Y进行排序
figure;  % 新建一个图形窗口
plot(Ysort,V1(id),'c','LineWidth',5);
hold on
plot(Ysort,V2(id),'k-.','LineWidth',2);
legend('经验分布函数','核分布估计', 'Location','NorthWest'); % 加标注框
xlabel('product');  % 为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（lnspot）');  % 为X轴加标签
ylabel('V（product）');  % 为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（spot）');  % 为X轴加标签
ylabel('V（product）');  % 为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的二元正态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)

Clayton-copula，Frank-copula和Gumbel-copula的实现方法类似，改函数中的参数即可。

部分说明

copulafit Fit a parametric copula to data.
RHOHAT = copulafit('Gaussian', U) returns an estimate RHOHAT of the matrix
of linear correlation parameters for a Gaussian copula, given data in U. U
is an N-by-P matrix of values in (0,1), representing N points in the
P-dimensional unit hypercube.

[RHOHAT, NUHAT] = copulafit('t', U) returns an estimate RHOHAT of the matrix
of linear correlation parameters for a t copula, and an estimate NUHAT of
the degrees of freedom parameter, given data in U. U is an N-by-P matrix of
values in (0,1), representing N points in the P-dimensional unit hypercube.

[RHOHAT, NUHAT, NUCI] = copulafit('t', U) returns an approximate 95%
confidence interval for the degrees of freedom parameter for a t copula,
given data in U.

PARAMHAT = copulafit(FAMILY, U) returns an estimate PARAMHAT of the copula
parameter for an Archimedean copula specified by FAMILY, given data in U. U
is an N-by-2 matrix of values in (0,1), representing N points in the unit
square. FAMILY is 'Clayton', 'Frank', or 'Gumbel'.

[PARAMHAT, PARAMCI] = copulafit(FAMILY, U) returns an approximate 95%
confidence interval for the copula parameter from an Archimedean copula
specified by FAMILY, given data in U.

[...] = copulafit(..., 'Alpha', ALPHA) returns an approximate 100(1-ALPHA)%
confidence interval for the parameter estimate.

copulafit uses maximum likelihood to fit the copula to U. When U contains
data transformed to the unit hypercube by parametric estimates of their
marginal cumulative distribution functions, this is known as the Inference
Functions for Margins (IFM) method. When U contains data transformed by
the empirical CDF, this is known as Canonical Maximum Likelihood (CML).

[...] = copulafit('t', U, ..., 'Method', 'ApproximateML') fits a t copula by
maximizing an objective function, as suggested by Bouye et al., that
approximates the profile log-likelihood for the degrees of freedom parameter
nu, for large sample sizes. This method can be significantly faster than
using maximum likelihood, however, it should be used with caution because
the estimates and confidence limits may not be accurate for small or
moderate sample sizes. copulafit('t', U, ..., 'Method', 'ML') is equivalent
to the default maximum likelihood fit.

[...] = copulafit(..., 'Options', OPTIONS) specifies control parameters for
the iterative algorithm used to compute estimates. This argument can be
created by a call to STATSET. See STATSET('copulafit') for parameter names
and default values. This argument does not apply to the 'Gaussian' family.
结果如下：

偏度和峰度

Spot和product的偏度

xs =-0.4609

ys =0.0475

Spot和product的丰度

kx =2.2937

ky =1.8913

spot的正态性检验

h =0 p =0.5000 Jarque-Bera检验

h =0 p =0.7589 Kolmogorov-Smirnov检验

h =0 p =0.3117 Lilliefors检验

product的正态性检验

h =0 p =0.5000 Jarque-Bera检验

h =0 p =0.9875 Kolmogorov-Smirnov检验

h =0 p =0.5000 Lilliefors检验

Copula函数	参数α	Kendall秩相关系数	Spearman相关系数	平方欧氏距离D
Gaussian-copula	0.9562	0.8110	0.9521	0.0026
t-copula	0.9696	0.8425	0.9547	0.0015
Clayton-copula	5.9728	0.3235	0.9097	0.0059
Frank-copula	19.2954	0.8104	0.9550	0.0035
Gumbel-copula	6.0996	0.8361	0.9608	0.0024

根据平方欧氏距离，选择t-copula作为拟合

Copula函数实战(Matlab)相关推荐

帮我用matlab写一段copula函数预测
非常抱歉,我需要更具体的信息才能帮助您编写适当的Copula函数.请提供以下信息: 您要预测的是什么?例如,您要预测股票价格.天气模式或其他类型的数据? 您想使用哪种Copula模型?例如,Gauss ...
matlab做copula,Copula理论及MATLAB应用实例
Copula理论及MATLAB应用实例 %-------------------------------------------------------------------------- % Co ...
R语言ggplot2可视化分面图使用facet_wrap函数和facet_grid函数实战
R语言ggplot2可视化分面图使用facet_wrap函数和facet_grid函数实战目录 R语言ggplot2可视化分面图使用facet_wrap函数和facet_grid函数实战
R语言数据纵向合并rbind函数实战（以及rbind.fill函数合并两个数据列不同的dataframe）
R语言数据纵向合并rbind函数实战(以及rbind.fill函数合并两个数据列不同的dataframe) 目录
R语言is.na函数实战（删除、替换、统计、条件判断等）
R语言is.na函数实战(删除.替换.统计.条件判断等) 目录 R语言is.na函数实战(删除.替换.统计.条件判断等) #NA.NaN.Nu
R语言NaN函数实战（计数、替换、删除）
R语言NaN函数实战(计数.替换.删除) 目录 R语言NaN函数实战(计数.替换.删除) #NA.NaN.Null区别?
R语言dplyr包排序及序号函数实战（row_number、ntile、min_rank、dense_rank、percent_rank、cume_dist）
R语言dplyr包排序及序号函数实战(row_number.ntile.min_rank.dense_rank.percent_rank.cume_dist) 目录 R语言dplyr包排序及序号函数实 ...
R语言dplyr包cumall函数、cumany函数和cummean函数实战
R语言dplyr包cumall函数.cumany函数和cummean函数实战目录 R语言dplyr包cumall函数.cumany函数和cummean函数实战 #导入dplyr包 #cumall函数
R语言数据横向合并cbind函数实战
R语言数据横向合并cbind函数实战目录 R语言数据横向合并cbind函数实战 #基本语法 # cbind横向为dataframe添加新的列
R语言case_when函数和cases函数实战
R语言case_when函数和cases函数实战目录 R语言case_when函数和cases函数实战 #仿真数据 # 导入dplyr包 # case_when函数

Copula函数实战(Matlab)

Copula函数实战(Matlab)相关推荐

最新文章

热门文章