使用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.也可以使用具有非零对角线术语的协方差矩阵来易于生成。

rho=.7；SigmaDep=sigma.^2.[1 rho；rho 1]

SigmaDep = 0.2500 0.1750 0.1750 0.2500

ZDep=mvnrnd（[0]，SigmaDep，n）；XDep=exp（ZDep）；

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

plot（xdep（：，1），xdep（:,2），'。'）;轴平等的；轴([0 5 0 5]);Xlabel（X1的）;ylabel（“X2”）;

很明显，在第二个数据集中，对于大的X1值和大的X2值有更大的关联趋势，对于小的值也有类似的趋势。这种依赖性由潜在的二元正态的相关参数决定。仿真得出的结论很可能取决于X1和X2的生成是否具有相关性。

在这种情况下，二元对数正态分布是一个简单的解决方案，当然也很容易推广到高维数和边际分布的情况不同的对数正态。其他的多元分布也存在，例如，多元t和Dirichlet分布分别用来模拟相关的t和beta随机变量。但是简单的多元分布的列表并不长，而且它们只适用于边缘都在同一家族(甚至完全相同的分布)的情况。在许多情况下，这可能是一个真正的限制。

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

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

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

根据定义，将正态CDF(这里用PHI表示)应用到一个标准正态随机变量会得到一个在区间[0,1]上一致的r.v.。看这个，如果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）；hist（z-3.75:.5:3.75）；xlim（[-4]）；title(“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）；hist（x，.25:.5:9.75）；title('1000模拟N(0,1)值转换为Gamma(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);情节(X (: 1) X (:, 2),'。'); 轴（[012-88]）；h1=gca；头衔(“1000模拟相关的t和Gamma值”）;Xlabel（'X1~伽马（2,1）'）;ylabel（“X2 ~ t(5)”）;次要情节(2、2、4);酒吧(ctr1 n1 1);轴([0 12 -max(n1)*1.1 0]);轴(“关”）;甘氨胆酸h2 =;次要情节(2 2 1);barh (ctr2 n2 1);轴([-max(n2)*1.1 0 -8 8]);轴(“关”); h3=gca；h1.位置=[0.35 0.35 0.55 0.55]；h2.位置=[.35.1.55.15]；h3.位置=[1.35.15.55]；彩色地图（[8.81]）；

等级相关系数

在这种构造中X1和X2之间的依赖性是由二元正态的相关参数决定的。然而,它是不确实，X1和X2的线性相关性为ρ。例如，在原始对数正态情况下，该相关性有一个闭合形式：

cor（X1，X2）=（exp（rho.*sigma.^2）-1）。/（exp（sigma.^2）-1）

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

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

粗略地说，这些等级相关性度量了一个rv的大小值与另一个rv的大小值的关联程度。然而，与线性相关系数不同的是，它们仅根据排名来衡量关联。因此，在任何单调变换下，秩相关都保持不变。特别地，刚才描述的变换方法保留了秩相关性。因此，知道二元正态Z的秩相关就确定了最终变换r.v的秩相关。虽然仍然需要参数化潜在的二元正态，但在描述r.v之间的相关性时，Kendall的tau或Spearman的rho更有用。，因为它们对边际分布的选择是不变的。

事实证明，对于二元正态，Kendall's或Spearman's和线性相关系数之间存在一个简单的1-1映射:

rho = (2/pi)*arcsin(rho) or rho = sin(Tau *pi/2) rho = (6/pi)*arcsin(rho) or rho = 2*sin(rho_s*pi/6)

子地块（1,1,1）；rho=-1:01:1；τ=2.*asin（rho）。/pi；rho_=6.*asin（rho./2）。/pi；地块（rho，tau，' - ', rho_sρ' - '，[-1 1]，[ -  1 1]，“k:”)；轴（[-11-11]）；xlabel(的ρ）;ylabel（等级相关系数的）;传奇(肯德尔”年代τ那斯皮尔曼”年代rho_s那“位置”那“西北”）;

因此，无论X1和X2的边际分布如何，只要为Z1和Z2之间的线性相关选择正确的参数值，就可以很容易地建立X1和X2之间所需的秩相关。

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

Copulas.

上述构造的第一步定义了所谓的copula，特别是高斯copula。二元copula只是两个随机变量的概率分布，每个随机变量的边际分布是一致的。这两个变量可能完全独立，在确定性上相关（例如，U2=U1）二元高斯连接函数族由Rho=[1 Rho；Rho 1]线性相关矩阵参数化。U1和U2在Rho接近+/-1时接近线性相关性，在Rho接近零时接近完全独立性。

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

n = 500;Z = mvnrnd([0 0]， [1.8;1。8),n);U = normcdf (Z, 0,1);次要情节(2 2 1);情节(U (: 1), U (:, 2),'。'）;标题(‘rho=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；-.11]，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),'。'）;标题('rho=-0.8'）;Xlabel（‘U1’）;ylabel（“U2”）;

U1和U2之间的相关性与X1=G（U1）和X2=G（U2）的边际分布完全不同。可以给出X1和X2任何边际分布，并且仍然具有相同的秩相关。这是copula的主要吸引力之一——它们允许独立的依赖性和边际分布。

t连系动词

通过从二元t分布开始，并使用相应的t CDF进行变换，可以构造一个不同的连接函数族。二元t分布用Rho（线性相关矩阵）和nu（自由度）参数化。因此，例如，我们可以说t（1）或t（5）copula，分别基于一个和五个自由度的多元t。

不同水平rho的一些模拟随机值的散点图说明了t（1）copula的不同可能性范围：

n=500；nu=1；T=mvtrnd（[1.8；.81]，nu，n）；U=tcdf（T，nu）；子图（2,2,1）；绘图（U（：，1），U（：，2），'。'）;标题(‘rho=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]， n);U = tcdf (T,ν);次要情节(2、2、3);情节(U (: 1), U (:, 2),'。'）;标题('rho = -0.1'）;Xlabel（‘U1’）;ylabel（“U2”)T=mvtrnd（[1-.8；-.81]，nu，n）；U=tcdf（T，nu）；子图（2,2,4）；绘图（U（：，1），U（：，2），'。'）;标题('rho=-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 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); subplot(2,2,2); plot(X(:,1),X(:,2),'。'）;轴([0 15 -10 10]);甘氨胆酸h1 =;标题(“1000模拟相关的t和Gamma值”）;Xlabel（'X1~伽马（2,1）'）;ylabel（“X2 ~ t(5)”）;次要情节(2、2、4);酒吧(ctr1 n1 1);轴([0 15 -max(n1)*1.1 0]);轴(“关”）;甘氨胆酸h2 =;次要情节(2 2 1);barh (ctr2 n2 1);轴([-max(n2)*1.1 0 -10 10]);轴(“关”); h3=gca；h1.位置=[0.35 0.35 0.55 0.55]；h2.位置=[.35.1.55.15]；h3.位置=[1.35.15.55]；彩色地图（[8.81]）；

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

高阶连系动词

高斯copula和t copula被称为椭圆copula。很容易将椭圆copula推广到更高的维数。例如，我们可以使用高斯copula使用Gamma（2,1）、Beta（2,2）和t（5）边缘模拟来自三变量分布的数据，如下所示。

子批次（1,1,1）；n=1000；Rho=[1.4.2；.41-.8；.2-.81]；Z=mvnrnd（[0]，ρ，n）；U=标准CDF（Z，0,1）；X=[gaminv（U（：，1），2,1）betainv（U（：，2），2,2）tinv（U（：，3），5）]；图3（X（：，1），X（：，2），X（：，3），'。'）;网格在；视图([15]-55年);Xlabel（‘U1’）;ylabel（“U2”）;zlabel (U3的）;

请注意，线性相关参数和(例如，Kendall’s)之间的关系适用于这里使用的相关矩阵中的每一项。我们可以验证数据的样本秩相关性近似等于理论值。

tauTheoretical = 2π* asin(ρ)。/

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

tauSample = corr (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)每个变量的边际分布

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

负载储存者脑袋=大小(股票,1);次要情节(2,1,1);嘘(股票(: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）；楼梯（（1:nobs）/nobs，invCDF1，'B')；持有在；楼梯(invCDF2(1:脑袋)/脑袋,'r')；持有关；传奇(X1的那“X2”）;Xlabel（'累积概率'）;ylabel（“X”）;

对于模拟，我们可能想要实验不同的copula和相关性。在这里，我们将使用带有相当大的负相关参数的二元t(5)联系符。

n = 1000;ρ=。8;ν= 5;T = mvtrnd([1;Rho 1]， nu, n);U = tcdf (T,ν);X = [invCDF1(装天花板(n1 * U (: 1))) invCDF2(装天花板(n2 * U (:, 2))));[n1, ctr1] =嘘(X (: 1), 10);[n2, ctr2] =嘘(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-最大值（n1）*1.1 0]）；轴(“关”）;甘氨胆酸h2 =;次要情节(2 2 1);barh (ctr2 n2 1);轴([-max(n2)*1.1 0 -3.5 3.5]);轴(“关”); h3=gca；h1.位置=[0.35 0.35 0.55 0.55]；h2.位置=[.35.1.55.15]；h3.位置=[1.35.15.55]；彩色地图（[8.81]）；

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