用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),“。”）;轴平等的；轴([0 5 0 5]);包含(X1的）;ylabel (“X2”）;

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

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

SigmaDep = 0.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和Dirichlet分布分别用来模拟相关的t和beta随机变量。但是简单的多元分布的列表并不长，而且它们只适用于边缘都在同一家族(甚至完全相同的分布)的情况。在许多情况下，这可能是一个真正的限制。

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

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

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

根据定义，将正态CDF(这里用PHI表示)应用到一个标准正态随机变量会得到一个在区间[0,1]上一致的r.v.。看这个，如果Z是标准正态分布，那么U = PHI(Z)的CDF是

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

这就是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)随机值'）;包含(“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);嘘(x,二十五分:.5:9.75);标题('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),“。”）;轴([0 12 -8 8]);甘氨胆酸h1 =;标题(“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 =;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之间的依赖性是由二元正态的相关参数决定的。然而,它是不确实X1和X2的线性相关是。例如，在原始对数正态情况下，这种相关性有一个封闭形式:

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

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

这是因为线性相关系数表示线性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);ρ= 1:.01:1;τ= 2 * asin(ρ)。/π;rho_s = 6。* asin (rho. / 2)。/π;情节(ρ,τ,“- - -”, rho_sρ“- - -”，[-1 1]，[-1 1]，凯西:”）;轴([-1 1 -1 1]);包含(的ρ）;ylabel (等级相关系数的）;传奇(肯德尔”年代τ，斯皮尔曼”年代rho_s，“位置”，“西北”）;

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

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

连系动词

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

对于不同级别的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),“。”）;标题(ρ= 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;-．1, 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),“。”）;标题(ρ= -0.8的）;包含(‘U1’）;ylabel (“U2”）;

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

t连系动词

从二元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的）;包含(‘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;-．1 1]， n);U = tcdf (T,ν);次要情节(2、2、3);情节(U (: 1), U (:, 2),“。”）;标题(ρ= -0.1的）;包含(‘U1’）;ylabel (“U2”）;T = mvtrnd([1 -.8;-．8 1]， n);U = tcdf (T,ν);次要情节(2、2、4);情节(U (: 1), U (:, 2),“。”）;标题(ρ= -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 =;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),“。”）;网格在；视图([15]-55年);包含(‘U1’）;ylabel (“U2”）;zlabel (U3的）;

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

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

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

负载stockreturns脑袋=大小(股票,1);次要情节(2,1,1);嘘(股票(:1)10);包含(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);楼梯(invCDF1(1:脑袋)/脑袋,“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 -max(n1)*1.1 0]);轴(“关闭”）;甘氨胆酸h2 =;次要情节(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。