使用copula模拟相关随机变量

打开脚本

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

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

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

模拟输入之间的依赖关系

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

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

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

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

N = 1000;σ = .5;SigmaInd = sigma。^2 .* [10 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]);包含(X1的）;ylabel (“X2”）;

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

Rho = .7;SigmaDep = sigma。^2 .* [1 rho;ρ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值相关联，与较小的值相关联也是如此。这种依赖性是由相关参数决定的，rho，潜在的双变量正态。从模拟中得出的结论很可能取决于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)随机变量，得到一个分布恰好为F的r.v.，这被称为反演方法。对于正向情况，这个证明本质上与上面的证明相反。另一个直方图说明了gamma分布的转换。

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;rho 1) U = [PHI(Z1) PHI(Z2)] X = [G1(U1) G2(U2)]

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

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

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

[n1,ctr1] = hist(X(:，1)，20);[n2,ctr2] = hist(X(:，2)，20);次要情节(2,2,2);情节(X (: 1) X (:, 2),“。”）;轴([0 12 -8 8]);H1 = gca;标题(“1000个模拟相关t和Gamma值”）;包含('X1 ~ Gamma(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 = gca;h1。位置= [0.35 0.35 0.55 0.55];h2。位置=[。35 .1 .55 .15];h3。位置=[。1 .35 .15 .55]; colormap([.8 .8 1]);

等级相关系数

在这个结构中，X1和X2之间的相关性由相关参数rho决定，rho是基础二元正态的相关参数。然而，它是不X1和X2的线性相关性是正确的。例如，在原始的对数正态情况下，这种相关性有一个封闭形式:

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

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

这是因为线性相关系数表示线性r.v之间的依赖关系。，当非线性变换应用于这些r.v。，线性相关性不保留。相反，等级相关系数，如肯德尔的tau或斯皮尔曼的rho，是更合适的。

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

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

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

次要情节(1 1 1);Rho = -1:.01:1;Tau = 2.*asin(rho)./pi;Rho_s = 6.*asin(Rho_s ./2)./pi;情节(ρ,τ,“- - -”, rho_sρ“- - -”，[-1 1]，[-1 1]，凯西:”）;轴([-1 1 -1 1]);包含(的ρ）;ylabel (“等级相关系数”）;传奇(肯德尔”年代τ，斯皮尔曼”年代rho_s，“位置”，“西北”）;

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

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

连系动词

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

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

N = 500;Z = mvnrnd([0 0]， [1 .8;.8 1]， 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),“。”）;标题('rho = 0.1'）;包含(‘U1’）;ylabel (“U2”）;Z = mvnrnd([0 0]， [1 - 1;-．1 1]， n);U = normcdf(Z,0,1);次要情节(2、2、3);情节(U (: 1), U (:, 2),“。”）;标题('rho = -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进行转换，可以构造不同的copulas族。二元t分布用Rho(线性相关矩阵)和nu(自由度)参数化。因此，例如，我们可以谈论t(1)或t(5)联结，分别基于一个和五个自由度的多元t。

不同水平rho的一些模拟随机值的散点图说明了t(1) copulas的不同可能性范围:

N = 500;Nu = 1;T = mvtrnd([1 .8;.8 1]， 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,nu);次要情节(2,2,2);情节(U (: 1), U (:, 2),“。”）;标题('rho = 0.1'）;包含(‘U1’）;ylabel (“U2”）;T = mvtrnd([1 - 1;-．11 1]， nu, n);U = tcdf(T,nu);次要情节(2、2、3);情节(U (: 1), U (:, 2),“。”）;标题('rho = -0.1'）;包含(‘U1’）;ylabel (“U2”）;T = mvtrnd([1 -.8;-．8 1]， nu, n);U = tcdf(T,nu);次要情节(2、2、4);情节(U (: 1), U (:, 2),“。”）;标题('rho = -0.8'）;包含(‘U1’）;ylabel (“U2”）;

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

与高斯联结一样，任何边缘分布都可以施加在t联结上。例如，使用具有1个自由度的t联结，我们可以再次从具有Gam(2,1)和t(5)边缘的二元分布生成随机向量:

次要情节(1 1 1);N = 1000;Rho = .7;Nu = 1;T = mvtrnd([1 rho;Rho 1]， nu, n);U = tcdf(T,nu);X = [gaminv(U(:，1)，2,1) tinv(U(:，2)，5)];[n1,ctr1] = hist(X(:，1)，20);[n2,ctr2] = hist(X(:，2)，20); subplot(2,2,2); plot(X(:,1),X(:,2),“。”）;轴([0 15 -10 10]);H1 = gca;标题(“1000个模拟相关t和Gamma值”）;包含('X1 ~ Gamma(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 = gca;h1。位置= [0.35 0.35 0.55 0.55];h2。位置=[。35 .1 .55 .15];h3。位置=[。1 .35 .15 .55]; colormap([.8 .8 1]);

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

高阶连系动词

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

次要情节(1 1 1);N = 1000;Rho = [1 .4 .2;.4 1 -.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的）;

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

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

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

Copulas和经验边际分布

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

1) copula族(和任何形状参数)，2)变量之间的秩相关，3)每个变量的边际分布

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

负载stockreturnsNobs =规模(股票，1);次要情节(2,1,1);嘘(股票(:1)10);包含(X1的）;ylabel (“频率”）;次要情节(2,1,2);嘘(股票(:,2),10);包含(“X2”）;ylabel (“频率”）;

(这两个数据向量的长度相同，但这并不重要。)

我们可以为每个数据集分别拟合一个参数模型，并使用这些估计值作为我们的边际分布。然而，参数化模型可能不够灵活。相反，我们可以使用一个经验模型来计算边际分布。我们只需要一种计算CDF逆的方法。

这些数据集的经验逆CDF只是一个阶梯函数，阶梯值为1/nobs, 2/nobs，…1.阶高就是排序后的数据。

invCDF1 = sort(股票(:，1));N1 =长度(stocks(:，1));invCDF2 = sort(stocks(:，2));N2 = length(stocks(:，2));次要情节(1 1 1);楼梯(invCDF1(1:脑袋)/脑袋,“b”）;持有在；楼梯(invCDF2(1:脑袋)/脑袋,“r”）;持有从；传奇(X1的，“X2”）;包含(“累积概率”）;ylabel (“X”）;

对于模拟，我们可能需要使用不同的copulas和相关性进行实验。在这里，我们将使用具有相当大的负相关参数的二元t(5) copula。

N = 1000;Rho = -.8;Nu = 5;T = mvtrnd([1 rho;Rho 1]， nu, n);U = tcdf(T,nu);X = [invCDF1(装天花板(n1 * U (: 1))) invCDF2(装天花板(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 = gca;标题(“1000个模拟依赖值”）;包含(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 = gca;h1。位置= [0.35 0.35 0.55 0.55];h2。位置=[。35 .1 .55 .15];h3。位置=[。1 .35 .15 .55]; colormap([.8 .8 1]);

模拟数据的边缘直方图与原始数据的边缘直方图非常接近，并且随着我们模拟更多的值对而变得相同。请注意，这些值是从原始数据中提取的，由于每个数据集中只有100个观测值，因此模拟数据有些“离散”。克服这一问题的一种方法是在最终的模拟值中添加少量的随机变化，可能是正态分布的。这相当于使用经验逆CDF的平滑版本。