MATLAB利用Copula函数进行频率分析

0 Copula函数定义

1 Copula函数

按照构造方法的不同Copula函数可分为三类：椭圆Copula函数、Archimedean Copula函数和二次型Copula函数。

1.1 椭圆Copula函数

Fang等首先提出了椭圆Copula函数的概念，并对其进行了详细的阐述。Gaussian Copula函数和Student t Copula函数是最常用的椭圆Copula函数。

1.1.1 Gaussian Copula函数

1.1.2 t Copula函数

1.2 阿基米德Copula函数

Genest和Mackay(1986)给出了阿基米德Copula（Archimedean Copula）分布函数的定义。Archimedean Copula函数主要分为对称型和非对称型两大类。
Archimedean Copula函数仅含一个参数，结构简单、形式多样，在水文方面的应用极为广泛。
下面列出几种常用的二维Archimedean Copula函数：

Name	Note
Clayton	下尾相关，对分布的下尾部的变化比较敏感，能够捕捉到下尾相关的变化
Gumbel-Hougaard（GH）	上尾相关，对分布的上尾部的变化比较敏感，能够捕捉到上尾相关的变化
Ali-Mikhail-Haq（AMH）	尾相关
Frank	无尾相关，对称尾部，无法捕捉到随机变量之间的非对称的尾部相关
Joe	上尾相关

1.2.1 Clayton Copula函数

式中，u、v均为边缘分布函数；θ为Copula函数参数。

1.2.2 Gumbel-Hougaard（GH）Copula函数

式中，u、v均为边缘分布函数；θ为Copula函数参数。

1.2.3 Ali-Mikhail-Haq（AMH）Copula函数

式中，u、v均为边缘分布函数；θ为Copula函数参数。

1.2.4 Frank Copula函数

式中，u、v均为边缘分布函数；θ为Copula函数参数。

1.2.5 Joe Copula函数

1.2.6 其它

1.3 二次型Copula函数

2 变量相关性度量指标

2.1 Pearson线性相关系数

Pearson相关系数（Pearson Correlation Coefficient）是用来衡量两个数据集合是否在一条线上面，它用来衡量定距变量间的线性关系。

不同Pearson线性相关系数的示意图如下：

2.2 Kendall秩相关系数

2.3 Spearman秩相关系数

spearman相关系数是衡量两个变量的依赖性的非参数指标。 它利用单调方程评价两个统计变量的相关性。 如果数据中没有重复值， 并且当两个变量完全单调相关时，斯皮尔曼相关系数则为+1或−1。

2.4 尾部相关系数

在金融风险分析中，人们更关心的是随机变量的尾部相关性（Tail Dependence Coefficient, TDC），也就是当严格随机变量取较大的值或者较小的值时，它对另一个随机变量的取值是否有影响。
尾部相关包括上尾相关和下尾相关。

2.5 MATLAB实现

MATLAB函数如下：

% rho ρ Pearson线性相关系数
% tau τ Kendall秩相关系数
% rhos ρs Spearman秩相关系数
% lambda λ 尾部相关系数

rho = corr(X,Y ,'type','pearson');
tau =  corr(X,Y , 'type' , 'kendall');
rhos = corr(X,Y , 'type' , 'Spearman');  

根据上述函数，可依据所需，绘制相关系数图如下：

3 联合分布模型建立步骤

1.确定特征变量的边缘分布

2.确定Copula函数的参数

3.选择合适的Copula函数，并建立联合分布

4.进行频率分析，计算条件概率和重现期等

4 Copula函数实现代码

【建议】：可参考【MATLAB 统计分析与应用：40个案例分析】一书。

4.1 确定特征变量的边缘分布

4.1.1 K-S检验

K-S法进行检验：

[H,P,KSSTAT,CV] = kstest(X,cdf,alpha) 

H=0表示接受假设；H=1表示拒绝原假设，即拒绝指定分布。
P为原假设成立的概率
KSSTAT为测试统计量的值
CV为是否接受假设的临界值。

4.2 确定Copula函数的参数

4.3 优选Copula函数

4.3.1 Deviation Information Criterion, DIC

4.4 频率分析：条件概率分布、重现期等

4.5 The most-likely scenario

【理论学习】

the most-likely combination(MLE) point (d*,s*)

【实例探讨】
根据我的算例，我的实现结果如下：

绘图代码如下：

figure(1)
subplot(2,2,1)
hold on; box on;
contour(Yseries,Xseries,M',[2 2],'w');
DSS = GetContourMatrix(c,2);
h(2)=line( [ min(DSS{1,1}(2,:)) min(DSS{1,1}(2,:)) ],[ 50 200] ,'linestyle','--', 'Color',[0.4118 0.4118 0.4118], 'LineWidth', 1); 
h(3)=line( [ 0.5 1.8] ,[ min(DSS{1,1}(1,:)) min(DSS{1,1}(1,:)) ] ,'linestyle','--', 'Color',[0.4118 0.4118 0.4118], 'LineWidth', 1); 
[~,h(1)] = contour(Yseries,Xseries,M',[2 2],'k');
h(1).LineWidth = 1.5;
h(4)=plot(MLC(1,2),MLC(1,1), 'r.','markersize',25);
axis([ 0.5 1.75 50 190 ]);
xlabel("Severity");
ylabel("Duration");
hl = legend(h([1 4 ]),"Isoline of T_o_r","Most likely scenario");
set(hl,'Box','off');

subplot(2,2,2)
hold on; box on;
plot(ff1{1,1},aaa{1,1},'k:','linewidth',1.5);
index = find(aaa{1,1}==MLC(1,1));
plot(ff1{1,1}(index),aaa{1,1}(index),'r.','markersize',25);
axis([ 0 0.01 50 190  ]);
xlabel("f_D(d)");
ylabel("Duration");

subplot(2,2,3)
hold on; box on;
plot(bbb{1,1},ff2{1,1},'k:','linewidth',1.5);
index = find(bbb{1,1}==MLC(1,2));
plot( bbb{1,1}(index),ff2{1,1}(index),'r.','markersize',25 );
axis([ 0.5 1.75 0 2.2]);
xlabel("Severity");
ylabel("f_S(s)");

subplot(2,2,4)
hold on; box on;
plot(1:length(UVYFrank{1,1}),UVYFrank{1,1},'k:','linewidth',1.5);
index = find(bbb{1,1}==MLC(1,2));
plot( index ,UVYFrank{1,1}(index),'r.','markersize',25 );
ylabel("f(d,s)");

4.6 Bivariate uncertainty envelopes

总结

若需要完整代码需要可私信博主
此代码包括以下几个方面的内容：
1）边缘分布拟合：基于正态、指数、P-Ⅲ型等8种分布函数完成单变量拟合，包括相关拟合检验代码及拟合分布图
2）联合分布拟合：包括5种Copula函数，包括相关拟合检验代码及拟合分布图
3）成果部分：包括联合分布重现期计算、最可能场景及条件概率分布等
如有其它疑惑，也可私信博主或评论本博客

参考

1.书籍-《MATLAB统计分析与应用：40个案例分析》
2.博客-Copula函数
3.参考文献-Responses of precipitation and runoff to climate warming and implications for future drought changes in China
4.参考文献-Drought hazard transferability from meteorological to hydrological propagation