Copulas:生成相关样本- MATLAB和Simulink - MathWorks Bene万博1manbetxlux

Copulas:生成相关样本

连系动词是描述变量之间依赖关系的函数，并提供了一种方法来创建建模相关多元数据的分布。使用copula，您可以通过指定边际单变量分布来构造一个多元分布，然后选择一个copula来提供变量之间的相关结构。二元分布以及高维分布都是可能的。

确定模拟输入之间的相关性

蒙特卡罗模拟的设计决策之一是随机输入的概率分布的选择。为每个单独的变量选择一个分布通常很简单，但是决定输入之间应该存在哪些依赖关系可能就不那么简单了。理想情况下，模拟的输入数据应该反映您所知道的正在建模的真实数量之间的依赖性。然而，在模拟中可能只有很少或根本没有信息来建立任何依赖关系。在这种情况下，为了确定模型的灵敏度，对不同的可能性进行实验是有用的。

当它们的分布不是来自标准多元分布时，生成具有依赖性的随机输入可能会很困难。此外，一些标准多元分布只能对有限类型的依赖进行建模。使输入相互独立总是有可能的，虽然这是一个简单的选择，但它并不总是明智的，可能会导致错误的结论。

例如，金融风险的蒙特卡洛模拟可以有两个随机输入，代表不同的保险损失来源。你可以把这些输入建模成对数正态随机变量。一个合理的问题是，这两个输入之间的依赖性如何影响模拟结果。事实上，你可能从真实数据中知道，相同的随机条件会影响两个来源;在模拟中忽略这一点可能会得出错误的结论。

生成和指数正态随机变量

打开实时脚本

的lognrnd函数模拟独立的对数正态随机变量。在下面的示例中，使用mvnrnd函数生成n对独立的正态随机变量，然后取幂。注意这里使用的协方差矩阵是对角线的。

N = 1000;σ = .5;SigmaInd = sigma。^2 .* [10 0;0 1]

SigmaInd =2×20.2500 00 0.2500

rng (“默认”）;%用于再现性ZInd = mvnrnd([0 0]，SigmaInd,n);XInd = exp(ZInd);情节(XInd (: 1) XInd (:, 2),“。”)轴([0 5 0 5]平等的包含(X1的) ylabel (“X2”）

使用非零非对角线项的协方差矩阵也容易生成相关双变量对数正态随机变量。

Rho = .7;SigmaDep = sigma。^2 .* [1 rho;ρ1]

SigmaDep =2×20.2500 0.1750 0.1750 0.2500

ZDep = mvnrnd([0 0]，SigmaDep,n);XDep = exp(ZDep);

第二个散点图展示了这两个二元分布之间的差异。

情节(XDep (: 1) XDep (:, 2),“。”)轴([0 5 0 5]平等的包含(X1的) ylabel (“X2”）

很明显，在第二个数据集中有较大的值的趋势X1与…的大值相关联X2，对于较小的值也类似。相关参数 $ρ$ 基本的双变量正态决定了这种依赖性。从模拟中得出的结论很可能取决于您是否生成X1而且X2与依赖。在这种情况下，二元对数正态分布是一个简单的解;它很容易推广到更高维度的情况下，边际分布是不同的对数正态分布。

其他多元分布也存在。例如，多元变量t狄利克雷分布模拟依赖t和随机变量。但是简单多元分布的列表并不长，它们只适用于边缘值都属于同一族(甚至是完全相同的分布)的情况。在许多情况下，这可能是一个严重的限制。

构造相关二元分布

虽然上一节中讨论的构造创建了一个简单的二元对数正态函数，但它是为了说明一种更普遍适用的方法。

从二元正态分布生成值对。这两个变量之间有统计学上的依赖关系，每个变量都有正态的边际分布。
对每个变量分别应用变换(指数函数)，将边际分布变为对数正态分布。转换后的变量仍然具有统计依赖性。

如果能找到合适的变换，该方法可以推广到具有其他边际分布的依赖二元随机向量。事实上，构造这样一个变换的一般方法是存在的，尽管它不像单纯的求幂那么简单。

根据定义，将正态累积分布函数(cdf)(这里用Φ表示)应用于标准正态随机变量会得到在区间[0,1]上均匀的随机变量。看这个，如果Z有一个标准正态分布，那么CDFU=Φ(Z)是

$公关｛ U \leq u ｝＝公关｛ Φ （ Z ） \leq u ｝＝公关（ Z \leq Φ^{- 1} （ u ）｝＝ u$

这是Unif(0,1)随机变量的cdf。一些模拟正态值和转换值的直方图表明:

N = 1000;rng默认的再现率%Z = normrnd(0,1,n,1);%生成标准正常数据直方图(z, -3.75: .5:3.75,“FaceColor”,(。8 .8 1])绘制数据的直方图Xlim ([-4 4])1000个模拟N(0,1)随机值)包含(“Z”) ylabel (“频率”）

U = normcdf(z);%计算样本数据的CDF值图直方图(u . 05。1:.95,“FaceColor”,(。8 .8 1])绘制CDF值的直方图标题(1000个模拟N(0,1)值转换为Unif(0,1))包含(“U”) ylabel (“频率”）

借鉴单变量随机数生成理论，应用任意分布的逆cdf，F，赋给Unif(0,1)随机变量，得到的随机变量的分布恰好为F(见反演方法)．对于正向情况，这个证明本质上与前面的证明相反。另一个直方图说明了gamma分布的转换:

X = gaminv(u,2,1);%转换为gamma值图直方图(x)为:.5:9.75,“FaceColor”,(。8 .8 1])绘制伽马值的直方图标题(1000个模拟N(0,1)值转换为Gamma(2,1))包含(“X”) ylabel (“频率”）

您可以将此两步转换应用于标准双变量正态的每个变量，创建具有任意边际分布的相关随机变量。由于变换分别作用于每个组件，因此两个随机变量甚至不需要具有相同的边际分布。转换定义为:

$\begin{array}{l} Z ＝［ Z_{1} ， Z_{2} ］ \sim N （［ 0 ， 0 ］，［ \begin{matrix} 1 & ρ \\ ρ & 1 \end{matrix} ］） \\ U ＝［ Φ （ Z_{1} ）， Φ （ Z_{2} ）］ \\ X ＝［ G_{1} （ U_{1} ）， G_{2} （ U_{2} ）］ \end{array}$

在哪里G₁而且G₂是两个可能不同分布的逆CDFS。例如，下面从二元分布生成随机向量t₅Gamma(2,1)边缘值:

N = 1000;Rho = .7;Z = mvnrnd([0 0]，[1 rho;ρ1],n);U = normcdf(Z);X = [gaminv(U(:，1)，2,1) tinv(U(:，2)，5)];用直方图绘制数据的散点图图scatterhist (X (: 1), (:, 2),“方向”，“出”）

该图有直方图和散点图，以显示边际分布和相关性。

使用等级相关系数

相关参数，ρ，决定了之间的依赖关系X1而且X2在这个结构中。的线性相关X1而且X2不是ρ．例如，在原始的对数正态情况下，这种相关性的封闭形式是:

$c o r （ X 1 ， X 2 ）＝ \frac{e^{ρ σ^{2}} - 1}{e^{σ^{2}} - 1}$

哪个严格小于ρ,除非ρ正好是1。在更一般的情况下，比如/t构造，两者之间线性相关X1而且X2很难或不可能用…来表达ρ，但模拟显示同样的效果也会发生。

这是因为线性相关系数表示随机变量之间的线性依赖关系，当对这些随机变量进行非线性变换时，线性相关性就不存在了。相反，是一个等级相关系数，如肯德尔的τ或斯皮尔曼的ρ更合适。

粗略地说，这些等级相关性衡量的是一个随机变量的大值或小值与另一个随机变量的大值或小值相关联的程度。然而，与线性相关系数不同的是，它们只根据排名来衡量相关性。因此，在任何单调变换下，秩相关都保持不变。特别地，刚才描述的转换方法保留了秩相关。因此，知道二元正态的秩相关Z精确地确定了最终转换的随机变量的秩相关性，X．而线性相关系数，ρ，仍然需要参数化潜在的双变量正态，Kendall'sτ斯皮尔曼的ρ在描述随机变量之间的依赖关系时更有用，因为它们对边际分布的选择是不变的。

对于二元正态，Kendall之间有一个简单的一对一映射τ斯皮尔曼的ρ，为线性相关系数ρ:

$\begin{array}{l} τ ＝ \frac{2}{π} \arcsin （ ρ ）或 ρ =罪（ τ \frac{π}{2} ） \\ ρ_{年代} ＝ \frac{6}{π} \arcsin （ \frac{ρ}{2} ）或 ρ = 2罪（ ρ_{年代} \frac{π}{6} ） \end{array}$

下面的图表显示了这种关系。

Rho = -1:.01:1;Tau = 2.*asin(rho)./pi;Rho_s = 6.*asin(Rho_s ./2)./pi;情节(ρ,τ,“b -”，“线宽”, 2)在情节(ρ,rho_s“g -”，“线宽”，2) plot([-1 1]，[-1 1]，凯西:”，“线宽”，2)轴([-1 1 -1 1])的ρ) ylabel (“等级相关系数”)传说(肯德尔”年代{\ \τ}”，.．.斯皮尔曼”年代{\ \ rho_s}”，.．.“位置”，“西北”）

因此，很容易创建所需的排序相关性X1而且X2，不管它们的边际分布，通过选择正确的ρ参数值为之间的线性相关Z1而且Z2．

对于多元正态分布，斯皮尔曼秩相关几乎与线性相关相同。然而，一旦你转换到最终的随机变量，这就不成立了。

使用二元Copulas

上一节中描述的构造的第一步定义了所谓的二元高斯联结。copula是一个多元概率分布，其中每个随机变量在单位区间[0,1]上具有均匀的边际分布。这些变量可能是完全独立的，确定性相关的(例如，U2 = u1)，或介于两者之间。由于变量之间存在依赖关系的可能性，您可以使用copula为因变量构造一个新的多元分布。通过使用反演方法分别转换copula中的每个变量，可能使用不同的cdfs，得到的分布可以有任意的边际分布。当您知道不同的随机输入不是相互独立时，这种多元分布在模拟中通常是有用的。

统计和机器学习工具箱™功能计算:

概率密度函数(copulapdf)和累积分布函数(copulacdf)
由线性相关得出的秩相关(copulastat)，反之亦然(copulaparam）
随机向量(copularnd）
拟合数据的相关参数(copulafit）

例如，使用copularnd函数来创建随机值的散点图，从二元高斯联结的各个水平ρ，以说明不同依赖关系结构的范围。二元高斯关联函数族由线性相关矩阵参数化:

$Ρ ＝（ \begin{matrix} 1 & ρ \\ ρ & 1 \end{matrix} ）$

U1而且U2线性依赖的方法如下ρ趋于±1，ρ趋于0时趋于完全独立:

之间的依赖关系U1而且U2与的边际分布完全分离X1 = g (u1)而且X2 = g (u2)．X1而且X2可以给出任何边际分布，并且仍然具有相同的等级相关性。这是copula的主要吸引力之一——它们允许对依赖关系和边际分布进行单独的说明。你也可以计算pdf (copulapdf)及基金(copulacdf)表示一个连词。例如，这些图显示了的pdf和cdfρ=。8:

U1 = linspace(1e-3,1-1e-3,50);U2 = linspace(1e-3,1-1e-3,50);[U1,U2] = meshgrid(U1,U2);Rho = [1 .8;1。8);F = copulapdf(“t”, (U1 (:) U2(:)),ρ,5);f =重塑(f，大小(U1));图()冲浪(u1, u2,日志(f),“FaceColor”，的插值函数，“EdgeColor”，“没有”view([-15,20]) xlabel(‘U1’) ylabel (“U2”) zlabel (的概率密度）

U1 = linspace(1e-3,1-1e-3,50);U2 = linspace(1e-3,1-1e-3,50);[U1,U2] = meshgrid(U1,U2);F = copulacdf(“t”, (U1 (:) U2(:)),ρ,5);F =重塑(F，大小(U1));图()冲浪(u1, u2, F,“FaceColor”，的插值函数，“EdgeColor”，“没有”view([-15,20]) xlabel(‘U1’) ylabel (“U2”) zlabel (“累积概率”）

从一个二元变量出发，可以构造一个不同的联结函数族t分配与转换使用相对应tcdf。二元的t分布是参数化的P为线性相关矩阵，和ν，自由度。例如，你可以说at₁或者一个t₅基于多元变量的Copulat分别是1和5个自由度。

就像高斯联结，统计和机器学习工具箱函数t连系动词计算:

概率密度函数(copulapdf)和累积分布函数(copulacdf)
由线性相关得出的秩相关(copulastat)，反之亦然(copulaparam）
随机向量(copularnd）
拟合数据的相关参数(copulafit）

例如，使用copularnd函数从二元变量中创建随机值的散点图t₁各层次的Copulaρ，以说明不同依赖结构的范围:

N = 500;Nu = 1;rng (“默认”）再现率%U =共生(“t”,(1。8;8 1],ν,n);次要情节(2 2 1)情节(U (: 1), U (:, 2),“。”)标题('{\it\rho} = 0.8')包含(‘U1’) ylabel (“U2”) U = copularnd(“t”,(1。1;1。1],ν,n);次要情节(2,2,2)情节(U (: 1), U (:, 2),“。”)标题('{\it\rho} = 0.1')包含(‘U1’) ylabel (“U2”) U = copularnd(“t”约[1;-．11］，nu，n）;次要情节(2,2,3)情节(U (: 1), U (:, 2),“。”)标题('{\it\rho} = -0.1')包含(‘U1’) ylabel (“U2”) U = copularnd(“t”,(1。8;-．8 1]，nu, n);次要情节(2,2,4)情节(U (: 1), U (:, 2),“。”)标题('{\it\rho} = -0.8')包含(‘U1’) ylabel (“U2”）

一个tCopula具有均匀的边际分布U1而且U2，就像高斯联结一样。等级相关性τ或ρ_年代组件之间的tCopula也是相同的函数ρ对于高斯函数。然而，正如这些图所示，at₁连谱与高斯连谱有很大的不同，即使它们的分量具有相同的秩相关。区别在于它们的依赖结构。这并不奇怪，作为自由度参数ν做大了，at_νcopula接近相应的高斯copula。

与高斯关联函数一样，可以在a上施加任何边缘分布t连系动词。例如，使用at，您可以再次从具有Gamma(2,1)和的二元分布生成随机向量t₅的人使用copularnd:

N = 1000;Rho = .7;Nu = 1;rng (“默认”）再现率%U =共生(“t”,(1ρ;ρ1],ν,n);X = [gaminv(U(:，1)，2,1) tinv(U(:，2)，5)];图()scatterhist (X (: 1), (:, 2),“方向”，“出”）

与双变量Gamma/相比t之前构造的分布，是基于高斯联结的，这里构造的分布，是基于t₁Copula，具有相同的边际分布和变量之间相同的秩相关，但依赖结构非常不同。这说明了一个事实，即多元分布不是由它们的边际分布或相关性唯一定义的。在应用程序中选择特定的copula可以基于实际观测数据，或者可以使用不同的copula作为确定模拟结果对输入分布的灵敏度的方法。

高维Copulas

高斯和t接合点称为椭圆接合点。将椭圆联结推广到更高维度是很容易的。例如，模拟来自Gamma(2,1)、Beta(2,2)和的三变量分布的数据t₅边缘使用高斯联结和copularnd，详情如下:

N = 1000;Rho = [1 .4 .2;.4 1 -.8;2 -。8 1];rng (“默认”）再现率%U =共生(“高斯”ρ,n);X = [gaminv (U(: 1)、2、1)betainv (U (:, 2), 2, 2) tinv (U (:, 3), 5)];

绘制数据图。

次要情节(1 1 1)plot3 (X (: 1), (:, 2), X (:, 3),“。”网格)在View ([- 55,15]) xlabel(X1的) ylabel (“X2”) zlabel (“X3”）

注意，线性相关参数之间的关系ρ还有，比如，肯德尔的τ，对相关矩阵中的每一项都成立P这里使用。您可以验证数据的样本秩相关性近似等于理论值:

tauTheoretical = 2.*asin(Rho)./pi

tauTheoretical =3×31.0000 0.2620 0.1282 0.2620 1.0000 -0.5903 0.1282 -0.5903 1.0000

tauSample = corr(X，“类型”，“假象”）

tauSample =3×31.0000 0.2581 0.1414 0.2581 1.0000 -0.5790 0.1414 -0.5790 1.0000

阿基米德连系动词

统计和机器学习工具箱函数可用于三个二元阿基米德联结族:

克莱顿连系动词
弗兰克连系动词
甘力克连系动词

这些是直接根据其cdfs定义的单参数族，而不是使用标准的多元分布来构造性地定义。

比较这三个阿基米德联结函数与高斯和t二元联结，首先使用copulastat函数来寻找高斯或的秩相关tCopula与线性相关参数为0.8，然后使用copulaparam函数来找到对应于该等级相关性的克莱顿联结参数:

Tau = copulastat(“高斯”，.8日,“类型”，“假象”）

Tau = 0.5903

α = copulaparam(“克莱顿”τ,“类型”，“假象”）

Alpha = 2.8820

最后，绘制一个随机样本从克莱顿copula与copularnd．对Frank和Gumbel copulas重复相同的步骤:

N = 500;U =共生(“克莱顿”、α- n);次要情节(1,1)情节(U (: 1), U (:, 2),“。”）;标题(['Clayton Copula， {\it\alpha} = 'sprintf (' % 0.2 f 'α)])包含(‘U1’) ylabel (“U2”)阿尔法= copulaparam(“弗兰克”τ,“类型”，“假象”）;U =共生(“弗兰克”、α- n);次要情节(3、1、2)情节(U (: 1), U (:, 2),“。”)标题(“Frank Copula， {\it\alpha} =”sprintf (' % 0.2 f 'α)])包含(‘U1’) ylabel (“U2”)阿尔法= copulaparam(“甘力克”τ,“类型”，“假象”）;U =共生(“甘力克”、α- n);次要情节(3,1,3)情节(U (: 1), U (:, 2),“。”)标题(“Gumbel Copula， {\it\alpha} =”sprintf (' % 0.2 f 'α)])包含(‘U1’) ylabel (“U2”）

使用Copulas模拟相关的多元数据

要使用copula模拟依赖的多元数据，必须指定以下每一个:

copula族(以及任何形状参数)
变量之间的等级相关性
每个变量的边际分布

假设你有两只股票的回报数据，并且想要运行一个蒙特卡罗模拟，输入遵循与数据相同的分布:

负载stockreturnsNobs =规模(股票，1);次要情节(2,1,1)直方图(股票(:1),10日“FaceColor”,(。8 .8 1]) xlim([-3.5 3.5]) xlabel(X1的) ylabel (“频率”)子图(2,1,2)直方图(stocks(:，2)，10，“FaceColor”,(。8 .8 1]) xlim([-3.5 3.5]) xlabel(“X2”) ylabel (“频率”）

你可以为每个数据集分别拟合一个参数模型，并使用这些估计值作为边际分布。然而，参数化模型可能不够灵活。相反，您可以使用非参数模型来转换到边际分布。所需要的是一种计算非参数模型的逆cdf的方法。

最简单的非参数模型是经验cdf，由ecdf函数。对于离散的边际分布，这是合适的。然而，对于连续分布，应使用比步进函数计算更平滑的模型ecdf．一种方法是估计经验cdf，并用分段线性函数在步骤的中点之间进行插值。另一种方法是使用核平滑ksdensity．例如，将经验cdf与第一个变量的核平滑cdf估计进行比较:

[Fi,xi] = ecdf(stocks(:，1));图()楼梯(xi, Fi,“b”，“线宽”, 2)在Fi_sm = ksdensity(stocks(:，1)，xi，“函数”，“提供”，“宽度”,酒精含量);情节(xi, Fi_sm的r -，“线宽”(1.5)包含X1的) ylabel (“累积概率”)传说(“经验”，“平滑”，“位置”，“西北”网格)在

为了模拟，使用不同的copulas和相关性进行实验。这里，你将使用一个二元变量t具有相当小的自由度参数的Copula。对于相关性参数，可以计算数据的秩相关性。

Nu = 5;Tau = corr(stocks(:，1)，stocks(:，2)，“类型”，“假象”）

Tau = 0.5180

求对应的线性相关参数t连系动词使用copulaparam．

copulaparam(“t”，，“类型”，“假象”）

ρ = 0.7268

下一步,使用copularnd来生成随机值t使用非参数逆cdfs进行Copula和变换。的ksdensity函数允许你对分布进行核估计，并在copula点求逆CDF，这一切都在一步中完成:

N = 1000;U =共生(“t”,(1ρ;ρ1],ν,n);

X1 = ksdensity(stocks(:，1)，U(:，1)，.．.“函数”，“icdf”，“宽度”,酒精含量);X2 = ksdensity(stocks(:，2)，U(:，2)，.．.“函数”，“icdf”，“宽度”,酒精含量);

或者，当你有大量数据或需要模拟多个值集时，在区间(0,1)的值网格上计算逆cdf，并在关联点处使用插值来计算它，可能会更有效:

P = linspace(0.00001,0.99999,1000);G1 = ksdensity(stocks(:，1)，p，“函数”，“icdf”，“宽度”, 0.15);X1 = interp1(p,G1,U(:，1)，样条的）;G2 = ksdensity(stocks(:，2)，p，“函数”，“icdf”，“宽度”, 0.15);X2 = interp1(p,G2,U(:，2)，样条的）;scatterhist (X1, X2)“方向”，“出”）

模拟数据的边缘直方图是原始数据直方图的平滑版本。平滑的数量由输入的带宽控制ksdensity．

拟合Copulas到数据

打开实时脚本

这个例子展示了如何使用copulafit用数据校准接合器。生成数据Xsim矩阵中数据的分布“就像”(在边际分布和相关性方面)X的列拟合边际分布X，使用适当的CDF函数进行转换X来U，所以U值在0和1之间，使用copulafit使接合使与…接合U，生成新数据Usim从copula出发，并使用适当的CDF逆函数进行变换Usim来Xsim．

加载并绘制模拟股票收益数据。

负载stockreturnsX =股票(:，1);Y =股票(:，2);scatterhist (x, y,“方向”，“出”）

使用累积分布函数的核估计器将数据转换为关联尺度(单位平方)。

U = ksdensity(x,x，“函数”，“提供”）;V = ksdensity(y,y，“函数”，“提供”）;scatterhist (u, v,“方向”，“出”)包含(“u”) ylabel (“v”）

适合t连系动词。

[Rho,nu] = copulafit(“t”(u v),“方法”，“ApproximateML”）

ρ=2×21.0000 0.7220 0.7220 1.0000

Nu = 3.2017e+06

生成一个随机样本t连系动词。

R =共轭(“t”ρ,ν,1000);U1 = r(:，1);V1 = r(:，2);scatterhist (u1, v1,“方向”，“出”)包含(“u”) ylabel (“v”甘氨胆酸)组(get (,“孩子”)，“标记”，“。”）

将随机样本转换回数据的原始尺度。

X1 = ksdensity(x,u1，“函数”，“icdf”）;Y1 = ksdensity(y,v1，“函数”，“icdf”）;scatterhist (x1, y1,“方向”，“出”甘氨胆酸)组(get (,“孩子”)，“标记”，“。”）

如示例所示，copulas与其他分布拟合函数自然集成。