用copula模拟相关随机变量

开放脚本

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

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

近年来，联结在仿真模型中已经很流行了。copula是描述变量之间依赖关系的函数，并提供了一种创建分布的方法来建模相关的多元数据。通过使用关联式，数据分析师可以通过指定边际单变量分布来构建多元分布，并选择一个特定的关联式来提供变量之间的相关结构。二元分布，以及高维分布是可能的。在这个例子中，我们讨论了如何在MATLAB中使用统计和机器学习工具箱使用copula来生成依赖的多元随机数据。

模拟输入之间的依赖性

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

然而，当随机输入具有非标准多元分布的分布时，很难实际生成具有相关性的随机输入。此外，一些标准多元分布只能对非常有限类型的相关性进行建模。始终可以使输入独立，而这是一个简单的过程如果选择，它并不总是明智的，可能会导致错误的结论。

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

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

n = 1000;σ= 5;SigmaInd =σ。^2 .* [1;0 1]

SigmaInd = 0.2500 00 0.2500

ZInd = mvnrnd([0 0]， SigmaInd, n); / /输入XInd = exp (ZInd);情节(XInd (: 1) XInd (:, 2),“。”）;轴平等的；轴（[0505]）；xlabel(X1的）;ylabel (“X2”）;

相依二元对数正态r。使用非零非对角项的协方差矩阵也很容易生成。

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）；

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

情节(XDep (: 1) XDep (:, 2),“。”）;轴平等的；轴（[0505]）；xlabel(X1的）;ylabel (“X2”）;

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

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

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

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

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

根据定义，将正态CDF（此处用φ表示）应用于标准正态随机变量会产生在区间[0，1]上均匀的r.v。要了解这一点，如果Z具有标准正态分布，则U=φ（Z）的CDF为

公关U < =情况}={公关{φ(Z) < =情况}=公关{Z <φ= ^(1)(情况)}=情况,

这就是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）随机值”）;包含(“Z”）;ylabel (“频率”）;

u=normcdf（z）；历史（u.05:.1:.95）；标题('1000模拟N(0,1)值转换为U(0,1)'）;包含(“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)'）;包含(“X”）;ylabel (“频率”）;

这个两步转换可以应用到一个标准二元正态的每个变量，创建相关的r.v。具有任意的边际分布。因为变换作用于每一个单独的成分，两个结果的r.v.。美国甚至不需要拥有相同的边际分配。转换定义为

Z = [Z1 Z2] ~ N([0 0]，[1 rho;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值”）;包含('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）

它严格小于除非正好是1。然而，在更一般的情况下，例如上面的Gamma/t构造，X1和X2之间的线性相关性很难或不可能用表示表示，但可以使用模拟来显示同样的效果。

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

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

事实证明，对于二元正态，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，rho_s，“- - -”，[-1 1]，[-1 1]，“k:”)；轴（[-11-11]）；xlabel(的ρ）;ylabel (等级相关系数的）;传奇(肯德尔”年代τ，斯皮尔曼”年代rho_s，“位置”，“西北”）;

因此，通过为Z1和Z2之间的线性相关性选择正确的rho参数值，很容易在X1和X2之间创建所需的秩相关性，而不管它们的边缘分布如何。

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

连系动词

上述构造的第一步定义了所谓的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’）;包含(‘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的）;包含(‘U1’）;ylabel (“U2”)Z=mvnrnd（[0,0]，[1-.1；-.11]，n）；U=normcdf（Z，0,1）；子图（2,2,3）；绘图（U（：，1），U（：，2），“。”); 头衔(ρ= -0.1的）;包含(‘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'）;包含(‘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’）;包含(‘U1’）;ylabel (“U2”）;T = mvtrnd([1 .1;1 1]， nu, n);U = tcdf (T,ν);次要情节(2,2,2);情节(U (: 1), U (:, 2),“。”); 头衔(ρ= 0.1的）;包含(‘U1’）;ylabel (“U2”); T=mvtrnd（[1-.1；-.11]，nu，n）；U=tcdf（T，nu）；子批次（2,2,3）；图（U（：，1），U（：，2），“。”); 头衔(ρ= -0.1的）;包含(‘U1’）;ylabel (“U2”)T=mvtrnd（[1-.8；-.81]，nu，n）；U=tcdf（T，nu）；子图（2,2,4）；绘图（U（：，1），U（：，2），“。”); 头衔('rho=-0.8'）;包含(‘U1’）;ylabel (“U2”）;

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

与高斯copula一样，任何边际分布都可以加在t copula上。例如，使用一个具有1个自由度的t联系符，我们可以再次从具有Gam(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值”）;包含('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年);包含(‘U1’）;ylabel (“U2”）;zlabel (“U3”）;

请注意，线性相关参数rho和（例如）Kendall的tau之间的关系适用于此处使用的相关矩阵rho中的每个条目。我们可以验证数据的样本秩相关性近似等于理论值。

tauTheoretical = 2π* asin(ρ)。/

tauTheoretical = 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

copula和经验边际分布

为了使用copula模拟相依的多元数据，我们已经看到需要指定

1） copula族（和任何形状参数），2）变量之间的秩相关性，以及3）每个变量的边际分布

假设我们有两组股票返回数据，并且我们希望运行一个蒙特卡罗模拟，其输入遵循与我们的数据相同的分布。

负载stockreturnsnobs=大小（库存，1）；子批次（2,1,1）；历史（库存（：，1），10）；xlabel(X1的）;ylabel (“频率”）;次要情节(2,1,2);嘘(股票(:,2),10);包含(“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”）;包含(“累积概率”）;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个模拟相关值”）;包含(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。