使用Copulas模拟依赖随机变量

打开脚本

当变量之间存在复杂的关系，或者单个变量来自不同的分布时，如何使用copula从多元分布中生成数据。

MATLAB®是一个理想的工具运行模拟，包括随机输入或噪声。统计和机器学习工具箱™提供了根据许多常见的单变量分布创建随机数据序列的函数。工具箱还包括一些从多元分布生成随机数据的函数，如多元正态分布和多元t。然而，没有内置的方法来为所有边际分布生成多元分布，或者在单个变量来自不同分布的情况下生成多元分布。

最近，Copulas在仿真模型中变得流行。Copulas是描述变量之间依赖性的函数，并提供创建分布以模拟相关多变量数据的方法。使用Copula，数据分析师可以通过指定边际单变量分布来构造多变量分布，并选择特定的Copula以提供变量之间的相关结构。双变量分布以及更高尺寸的分布。在此示例中，我们讨论如何使用Copulas使用统计和机器学习工具箱在Matlab中生成依赖多元随机数据。

模拟输入之间的依赖性

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

然而，当它们具有不来自标准多变量分布的分布时，可能难以实现具有依赖性的随机输入。此外，一些标准多变量分布可以仅模拟非常有限的依赖性。它总是可以使输入独立，而这是一个简单的选择，它并不总是明智的，可以导致错误的结论。

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

独立对数正态随机变量的模拟是平凡的。最简单的方法是使用lognrnd.函数。这里，我们用mvnrnd.函数生成n对独立的正常随机变量，然后对它们进行指数。请注意，这里使用的协方差矩阵是对角线，即Z列之间的独立性。

n = 1000;Sigma = .5;sigmaind = sigma。^ 2。* [1 0;0 1]

SigmaInd = 0.2500 00 0.2500

ZInd = mvnrnd([0 0]， SigmaInd, n); / /输入XInd = exp (ZInd);情节(XInd (: 1) XInd (:, 2),'。'）;轴平等的；轴（[0 5 0 5]）;Xlabel（'x1'）;ylabel（'x2'）;

依赖于双变量Lognormal R.v.也可以使用具有非零对角线术语的协方差矩阵来易于生成。

ρ= 7;SigmaDep =σ。^2 .* [1;ρ1]

SIGMADEP = 0.2500 0.1750 0.1750 0.2500

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

第二散点图示出了这两个双变量分布之间的差异。

plot（xdep（：，1），xdep（:,2），'。'）;轴平等的；轴（[0 5 0 5]）;Xlabel（'x1'）;ylabel（'x2'）;

很明显，对于X1的大值，在X1的大值中有更多的趋势，与X2的大值相关，并且类似于小值。这种依赖是由相关的参数，RON，底层相当于正常的相关参数。从模拟中得出的结论很好地取决于依赖是否产生X1和X2。

在这种情况下，双变量的Lognormal分布是一种简单的解决方案，当然当然很容易地推广到更高的尺寸和下限分布的情况不同的浪核网络。其他多变量分布也存在，例如，多变量T和Dirichlet分布用于分别模拟相关的T和β随机变量。但简单的多变量分布的列表不长，它们仅适用于边际的案件，其中包括在同一个家庭（甚至是完全相同的分布）。这可以是许多情况下的真正限制。

构造相依二元分布的一种更一般的方法

虽然上面创建二元对数正态函数的构造很简单，但它用于说明一种更普遍适用的方法。首先，我们从二元正态分布中生成值对。这两个变量之间有统计依赖关系，每个变量都有正态边际分布。接下来，对每个变量分别应用一个变换(指数函数)，将边际分布变为对数正态分布。转换后的变量仍然具有统计相关性。

如果可以找到合适的变换，可以推广该方法以创建与其他边缘分布的相关的双变量随机向量。事实上，确实存在构建这种转化的一般方法，尽管不像只是指数一样简单。

根据定义，将正常的CDF（由PHI表示）应用于标准的正常随机变量导致R.V.间隔是均匀的[0,1]。要看到这一点，如果z具有标准的正态分布，则u = phi（z）的CDF是

Pr {u <= u0} = pr {phi（z）<= u0} = pr {z <= phi ^（ -  1）（u0）} = u0，

这就是U(0,1) r.v的CDF。一些模拟正态值和转换值的直方图证明了这一事实。

n = 1000;z = normrnd (0, 1, n, 1);嘘(z, -3.75: .5:3.75);xlim ([4 4]);标题（'1000模拟N(0,1)随机值'）;Xlabel（“Z”）;ylabel（“频率”）;

u = normcdf (z);嘘(u . 05。1:.95);标题（'1000模拟N（0,1）值转换为U（0,1）'）;Xlabel（'U'）;ylabel（“频率”）;

现在，从单变量随机数的理论借用，将任何分发F的反向CDF应用于U（0,1）随机变量的R.V。其分发正好是F.这被称为反演方法。证明基本上与前案件的上述证明相反。另一直方图说明了对伽马分布的转换。

x = gaminv (u 2 1);嘘(x,二十五分:.5:9.75);标题（'1000模拟n（0,1）值转换为伽马（2,1）'）;Xlabel（“X”）;ylabel（“频率”）;

这两步变换可以应用于标准双变量正常的每个变量，创建依赖的R.v.具有任意边际分布。因为转换分别在每个组件上工作，所以两个结果为R.V.不需要具有相同的边际分布。转型定义为

z = [z1 z2]〜n（[0 0]，[1 rho; rho 1]）U = [phi（z1）phi（z2）] x = [g1（u1）g2（u2）]

其中G1和G2是两种可能不同分布的cdf的逆。例如，我们可以从带有t(5)和Gamma(2,1)边缘的二元分布中生成随机向量。

n = 1000;ρ= 7;Z = mvnrnd([0 0]， [1 rho;ρ1],n);U = normcdf (Z);X = [gaminv(U(:，1)，2,1) tinv(U(:，2)，5)];

该曲线与散点图具有直方图，以显示边缘分布和依赖。

[n1, ctr1] =嘘(X (: 1), 20);[n2, ctr2] =嘘(X (:, 2), 20);次要情节(2,2,2);plot（x（：，1），x（:,2），'。'）;轴([0 12 -8 8]);甘氨胆酸h1 =;标题（'1000模拟相关性T和伽马值'）;Xlabel（“X1 ~γ(2,1)”）;ylabel（'x2〜t（5）'）;次要情节(2、2、4);酒吧(ctr1 n1 1);轴（[0 12 -MAX（N1）* 1.1 0]）;轴（“关闭”）;H2 = GCA;次要情节(2 2 1);Barh（Ctr2，-N2,1）;轴（[ -  max（n2）* 1.1 0 -8 8]）;轴（“关闭”）;甘氨胆酸h3 =;h1。Position = [0.35 0.35 0.55 0.55];h2。位置=[。35.0.55.15];h3。位置=[。1.35 . 15.55]; colormap([.8 .8 1]);

等级相关系数

在该构造中X1和X2之间的依赖性由相关参数Rho，潜在的二偏见性正常的相关参数确定。但是，它是不是确实X1和X2的线性相关是。例如，在原始对数正态情况下，这种相关性有一个封闭形式:

软木(X1, X2) = (exp(ρ*σ^ 2)- 1)。/ (exp(σ^ 2)- 1)

除非rho正好一，否则这严格少于rho。然而，在更通用的情况下，例如上面的伽马/吨施工，X1和X2之间的线性相关性难以或不可能在ROO方面表达，但是仿真可以用来表明发生相同的效果。

这是因为线性相关系数表示线性r.v之间的依赖。，并对其进行非线性变换。，则线性相关不被保留。相反，等级相关系数，如肯德尔的tau或斯皮尔曼的rho更合适。

粗略地说，这些等级相关测量一个R.V的大或小值的程度。与另一个人的大或小值相关联。然而，与线性相关系数不同，它们仅在等级中测量关联。结果，在任何单调转化下都保留了等级相关性。特别地，刚刚描述的转换方法保留了等级相关性。Therefore, knowing the rank correlation of the bivariate normal Z exactly determines the rank correlation of the final transformed r.v.'s X. While rho is still needed to parameterize the underlying bivariate normal, Kendall's tau or Spearman's rho are more useful in describing the dependence between r.v.'s, because they are invariant to the choice of marginal distribution.

事实证明，对于Bivariate Normal，Kendall的Tau或Spearman的Rho之间存在简单的1-1映射，以及线性相关系数Rho：

tau =（2 / pi）* Arcsin（Rho）或rho = sin（tau * pi / 2）rho_s =（6 / pi）* Arcsin（rho / 2）或rho = 2 * sin（rho_s * pi / 6）

次要情节(1 1 1);ρ= 1:.01:1;τ= 2 * asin(ρ)。/π;rho_s = 6。* asin (rho. / 2)。/π;情节(ρ,τ,' - '，rho，rho_s，' - '，[-1 1]，[ -  1 1]，凯西:”）;轴([-1 1 -1 1]);Xlabel（'rho'）;ylabel（'等级相关系数'）;传奇('肯德尔'的tau'那'spearman'的rho_s'那'地点'那'西北'）;

因此，通过选择Z1和Z2之间的线性相关性的正确的RHO参数值，易于在X1和X2之间创建所需的秩相关性。

请注意，对于多变量正态分布，Spearman的秩相关几乎与线性相关性相同。但是，一旦我们转换为最终随机变量，就不是这样。

Copulas.

上面描述的构造的第一步定义了所谓的关联，具体地说，是高斯关联。二元关联函数是两个随机变量的概率分布，每个随机变量的边际分布是均匀的。这两个变量可能是完全独立的，也可能是确定相关的(例如U2 = U1)，或者介于两者之间。二元高斯copulas族的参数化是Rho = [1 Rho;1]，线性相关矩阵。当趋于+/- 1时U1和U2趋于线性相关，当趋于0时趋于完全独立。

各种级别的一些模拟随机值的散点图揭示了高斯共用的不同可能性范围：

n = 500;z = mvnrnd（[0 0]，[1 .8; .8 1]，n）;U = normcdf (Z, 0,1);次要情节(2 2 1);情节（U（：，1），U（:,2），'。'）;标题（ρ= 0.8的）;Xlabel（‘U1’）;ylabel（“U2”）;z = mvnrnd（[0 0]，[1 .1; .1 1]，n）;U = normcdf (Z, 0,1);次要情节(2,2,2);情节（U（：，1），U（:,2），'。'）;标题（ρ= 0.1的）;Xlabel（‘U1’）;ylabel（“U2”）;Z = mvnrnd([0 0]， [1 - 1;-．1, n);U = normcdf (Z, 0,1);次要情节(2、2、3);情节（U（：，1），U（:,2），'。'）;标题（'rho = -0.1'）;Xlabel（‘U1’）;ylabel（“U2”）;Z = mvnrnd([0 0]， [1 -.8;-．8 1], n);U = normcdf (Z, 0,1);次要情节(2、2、4);情节（U（：，1），U（:,2），'。'）;标题（ρ= -0.8的）;Xlabel（‘U1’）;ylabel（“U2”）;

U1和U2之间的依赖关系完全独立于X1 = G(U1)和X2 = G(U2)的边际分布。X1和X2是已知的任何边际分布，并且仍然具有相同的秩相关。这是copula的主要吸引力之一——它们允许独立的依赖性和边际分布。

T Copulas.

从二元t分布出发，利用相应的t CDF进行变换，可以构造出不同的copulas族。二元t分布是用线性相关矩阵和自由度参数化的。因此，例如，我们可以说一个t(1)或t(5) copula，基于多元t分别有一个和五个自由度。

对于不同级别的rho，一些模拟随机值的散点图说明了t(1) copula的不同可能性范围:

n = 500;ν= 1;T = mvtrnd([1.8;.8 1]， nu, n);U = tcdf (T,ν);次要情节(2 2 1);情节（U（：，1），U（:,2），'。'）;标题（ρ= 0.8的）;Xlabel（‘U1’）;ylabel（“U2”）;t = mvtrnd（[1 .1; .1 1]，nu，n）;U = tcdf (T,ν);次要情节(2,2,2);情节（U（：，1），U（:,2），'。'）;标题（ρ= 0.1的）;Xlabel（‘U1’）;ylabel（“U2”）;t = mvtrnd（[1 -.1; -1 1]，nu，n）;U = tcdf (T,ν);次要情节(2、2、3);情节（U（：，1），U（:,2），'。'）;标题（'rho = -0.1'）;Xlabel（‘U1’）;ylabel（“U2”）;T = mvtrnd([1 -.8;-．8 1]， n);U = tcdf (T,ν);次要情节(2、2、4);情节（U（：，1），U（:,2），'。'）;标题（ρ= -0.8的）;Xlabel（‘U1’）;ylabel（“U2”）;

t copula对于U1和U2有一致的边际分布，就像高斯copula一样。t copula中各分量之间的秩相关tau或rho_s也是与高斯函数相同的函数。然而，正如这些图所示，t(1) copula与高斯copula有很大的不同，即使它们的成分有相同的秩相关。区别在于它们的依赖结构。毫不奇怪，当自由度参数nu变大时，一个t(nu) copula接近相应的高斯copula。

与高斯谱一样，任何边缘分布都可以施加在T Copula上。例如，使用具有1度自由度的T copula，我们可以再次从双方与Gab（2,1）和T（5）边缘的双变量分布产生随机向量：

次要情节(1 1 1);n = 1000;ρ= 7;ν= 1;t = mvtrnd（[1 rho; rho 1]，nu，n）;U = tcdf (T,ν);X = [gaminv(U(:，1)，2,1) tinv(U(:，2)，5)];[n1, ctr1] =嘘(X (: 1), 20);[n2, ctr2] =嘘(X (:, 2), 20);次要情节(2,2,2); plot(X(:,1),X(:,2),'。'）;轴（[0 15 -10 10]）;甘氨胆酸h1 =;标题（'1000模拟相关性T和伽马值'）;Xlabel（“X1 ~γ(2,1)”）;ylabel（'x2〜t（5）'）;次要情节(2、2、4);酒吧(ctr1 n1 1);轴（[0 15 -MAX（N1）* 1.1 0]）;轴（“关闭”）;H2 = GCA;次要情节(2 2 1);Barh（Ctr2，-N2,1）;轴（[ -  max（n2）* 1.1 0 -10 10]）;轴（“关闭”）;甘氨胆酸h3 =;h1。Position = [0.35 0.35 0.55 0.55];h2。位置=[。35.0.55.15];h3。位置=[。1.35 . 15.55]; colormap([.8 .8 1]);

与之前构建的基于高斯copula的二元Gamma/t分布相比，本文构建的基于t(1) copula的分布具有相同的边际分布和变量之间相同的秩相关，但依赖结构有很大的不同。这说明多元分布不是由它们的边际分布或它们的相关性唯一定义的事实。在应用程序中，可以根据实际观测数据选择特定的连接子，也可以使用不同的连接子来确定模拟结果对输入分布的敏感性。

高阶金属

高斯函数和t函数称为椭圆函数。很容易将椭圆联结推广到更高维度。例如，我们可以用下面的高斯copula来模拟带有Gamma(2,1)， Beta(2,2)和t(5)边缘的三元分布的数据。

次要情节(1 1 1);n = 1000;Rho = [1 4 .2;1。4。8;2 -。8 1];Z = mvnrnd([0 0 0]， Rho, n);U = normcdf (Z, 0,1);X = [gaminv (U(: 1)、2、1)betainv (U (:, 2), 2, 2) tinv (U (:, 3), 5)];plot3 (X (: 1) X (:, 2), X (:, 3),'。'）;网格在；查看（[ -  55,15]）;Xlabel（‘U1’）;ylabel（“U2”）;Zlabel（'U3'）;

请注意，线性相关参数RON与例如KENDALL的TAU之间的关系，在此处使用的相关矩阵RHO中的每个条目。我们可以验证数据的样本等级相关性大致等于理论值。

tauTheoretical = 2π* asin(ρ)。/

TautheOrative = 1.0000 0.2620 0.1282 0.2620 1.0000 -0.5903 0.1282 -0.5903 1.0000

tausample = cor（x，“类型”那'肯德尔'）

Tausample = 1.0000 0.2655 0.1060 0.2655 1.0000 -0.6076 0.1060 -0.6076 1.0000

Copulas和经验边缘分布

要使用copula模拟依赖多变量数据，我们已经看到我们需要指定

1）Copula系列（和任何形状参数），2）变量之间的等级相关性，以及3）每个变量的边缘分布

假设我们有两套库存回报数据，我们希望使用蒙特卡罗模拟，其中输入与我们的数据相同的分布。

加载储存者nobs =尺寸（股票，1）;子图（2,1,1）;SIGR（股票（：，1），10）;Xlabel（'x1'）;ylabel（“频率”）;次要情节(2,1,2);嘘(股票(:,2),10);Xlabel（'x2'）;ylabel（“频率”）;

（这两个数据向量具有相同的长度，但这并不至关重要。）

我们可以将参数模型分别适用于每个数据集，并使用这些估计作为我们的边际分布。然而，参数模型可能不具有足够的灵活性。相反，我们可以使用对边缘分布的实证模型。我们只需要一种方法来计算反向CDF。

这些数据集的经验反向CDF只是一个阶段函数，有一个步骤1 / nobs，2 / nobs，...... 1.步高度只是排序数据。

invCDF1 =排序(股票(:1));n1 =长度(股票(:1));invCDF2 =排序(股票(:,2));n2 =长度(股票(:,2));次要情节(1 1 1);楼梯(invCDF1(1:脑袋)/脑袋,'B'）;持有在；楼梯(invCDF2(1:脑袋)/脑袋,'r'）;持有从；传奇('x1'那'x2'）;Xlabel（'累积概率'）;ylabel（“X”）;

对于模拟，我们可能希望尝试不同的Copulas和相关性。在这里，我们将使用具有相当大的负相关参数的双相变T（5）copula。

n = 1000;rho = -.8;nu = 5;t = mvtrnd（[1 rho; rho 1]，nu，n）;U = tcdf (T,ν);x = [INVCDF1（CEIL（N1 * U（：1）））INVCDF2（CEIL（N2 * U（：，2）））];[n1，ctr1] = hist（x（：，1），10）;[n2，ctr2] = hist（x（:,2），10）;次要情节(2,2,2);plot（x（：，1），x（:,2），'。'）;轴（[ -  3.5 3.5 -3.5 3.5]）;甘氨胆酸h1 =;标题（'1000模拟依赖价值'）;Xlabel（'x1'）;ylabel（'x2'）;次要情节(2、2、4);酒吧(ctr1 n1 1);轴([-3.5 3.5 -max(n1)*1.1 0]);轴（“关闭”）;H2 = GCA;次要情节(2 2 1);Barh（Ctr2，-N2,1）;轴（[ -  max（n2）* 1.1 0 -3.5 3.5]）;轴（“关闭”）;甘氨胆酸h3 =;h1。Position = [0.35 0.35 0.55 0.55];h2。位置=[。35.0.55.15];h3。位置=[。1.35 . 15.55]; colormap([.8 .8 1]);

模拟数据的边缘直方图与原始数据的边缘直方图紧密匹配，并且随着我们模拟更多的数值，将变得相同。请注意，值从原始数据中绘制，并且由于每个数据集中只有100个观察，因此模拟数据有点“离散”。克服这一点的一种方法是添加少量随机变化，可能通常分布到最终的模拟值。这相当于使用经验逆CDF的平滑版本。