仿真输入之间的相关性

蒙特卡罗模拟的设计决策之一是为随机输入选择概率分布。为每个变量选择一种分布往往很简单，但确定输入之间应该存在什么样的相关性却可能不那么简单。理想情况下，仿真的输入数据应反映要建模的实际数量之间已知的相关性。然而，判断仿真中的任何相关性时可以依据的信息可能很少或根本没有，在这种情况下，最好的做法是尝试不同的可能性，以确定模型的敏感性。

但是，当输入数据的分布不是标准多元分布时，很难真正生成具有相关性的随机输入。而且，有些标准多元分布只能对非常有限的几种相关性进行建模。我们可以始终将输入视为各自独立，这很简便，但并不总是合理的，有可能导致错误的结论。

例如，在有关金融风险的蒙特卡罗模拟中，可能有一些表示各种保险损失来源的随机输入。这些输入可以建模为对数正态随机变量。我们会想了解两两输入间的相关性对仿真结果有何影响。确实，从实际数据中可能已经知道，相同的随机条件可对两个损失来源都产生影响，如果在仿真中忽略这一点，可能会导致错误的结论。

独立的对数正态随机变量的仿真非常繁琐。最简单的方法是使用 lognrnd 函数。此处，我们将使用 mvnrnd 函数生成 n 对独立的正态随机变量，然后计算它们的幂。请注意，这里使用的协方差矩阵是对角矩阵，即 Z 的各列之间是独立的。

n = 1000;

sigma = .5;

SigmaInd = sigma.^2 .* [1 0; 0 1]

SigmaInd =

0.2500 0

0 0.2500

ZInd = mvnrnd([0 0], SigmaInd, n);

XInd = exp(ZInd);

plot(XInd(:,1),XInd(:,2),'.');

axis equal;

axis([0 5 0 5]);

xlabel('X1');

ylabel('X2');

生成具有相关性的二元对数正态随机变量也很简单，可使用非零项不在对角线上的协方差矩阵。

rho = .7;

SigmaDep = sigma.^2 .* [1 rho; rho 1]

SigmaDep =

0.2500 0.1750

0.1750 0.2500

ZDep = mvnrnd([0 0], SigmaDep, n);

XDep = exp(ZDep);

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

plot(XDep(:,1),XDep(:,2),'.');

axis equal;

axis([0 5 0 5]);

xlabel('X1');

ylabel('X2');

很明显，第二个数据集中较大的 X1 值更倾向于与较大的 X2 值关联，小值之间也是如此。这种相关性由基础二元正态分布的相关参数 rho 决定。从仿真中得出的结论可能很大程度上取决于生成的 X1 和 X2 是否存在相关性。

本示例中的二元对数正态分布是一个简单的解决方法，它当然可以泛化应用于更高维度以及边缘分布为其他对数正态分布的情况。也有一些多元分布，例如，多元 t 分布和 Dirichlet 分布，它们分别用于仿真相关的 t 随机变量和 beta 随机变量。但是，简单的多元分布并不多，而且仅适用于边缘分布都为同一族分布(甚至是完全相同的分布)的情形。在很多情况下，这可能会成为限制。

构造相关二元分布的更通用方法

尽管上面创建二元对数正态分布的构造很简单，但它说明了一种更普遍适用的方法。首先，我们从二元正态分布中生成值对。两两变量间存在统计相关性，而且每个变量都具有正态边缘分布。接下来，分别对每个变量进行变换(指数函数)，将边缘分布变成对数正态分布。变换后的变量仍具有统计相关性。

如果能够找到一种合适的变换，就可以将此方法泛化，从而为其他边缘分布生成相关二元随机向量。实际上，确实存在构造这种变换的一般方法，尽管不像求幂那样简单。

根据定义，将正态 CDF(这里用 PHI 表示)应用于标准正态随机变量将生成在区间 [0,1] 内均匀分布的随机变量。为说明这一点，假设 Z 具有标准正态分布，则 U = PHI(Z) 时，CDF 为

Pr{U <= u0} = Pr{PHI(Z) <= u0} = Pr{Z <= PHI^(-1)(u0)} = u0,

而这就是 U(0,1) 随机变量的 CDF。通过为仿真的正态值和仿真的变换值绘制的直方图证实了这一点。

n = 1000;

z = normrnd(0,1,n,1);

hist(z,-3.75:.5:3.75);

xlim([-4 4]);

title('1000 Simulated N(0,1) Random Values');

xlabel('Z');

ylabel('Frequency');

u = normcdf(z);

hist(u,.05:.1:.95);

title('1000 Simulated N(0,1) Values Transformed to U(0,1)');

xlabel('U');

ylabel('Frequency');

现在，借用一元随机数生成理论，将任意分布 F 的逆 CDF 应用于 U(0,1) 随机变量，可得到分布正好为 F 的随机变量。这就是逆变换法。其证明过程与上述正向情形的证明过程基本相反。另一个直方图说明了变换为 gamma 分布的结果。

x = gaminv(u,2,1);

hist(x,.25:.5:9.75);

title('1000 Simulated N(0,1) Values Transformed to Gamma(2,1)');

xlabel('X');

ylabel('Frequency');

这个两步变换方法可以应用于标准二元正态分布的每个变量，从而生成具有任意边缘分布的相关随机变量。由于变换分别作用于每个成分，因此生成的两个随机变量甚至不需要具有相同的边缘分布。变换定义为

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);

subplot(2,2,2);

plot(X(:,1),X(:,2),'.');

axis([0 12 -8 8]);

h1 = gca;

title('1000 Simulated Dependent t and Gamma Values');

xlabel('X1 ~ Gamma(2,1)');

ylabel('X2 ~ t(5)');

subplot(2,2,4);

bar(ctr1,-n1,1);

axis([0 12 -max(n1)*1.1 0]);

axis('off');

h2 = gca;

subplot(2,2,1);

barh(ctr2,-n2,1);

axis([-max(n2)*1.1 0 -8 8]);

axis('off');

h3 = gca;

h1.Position = [0.35 0.35 0.55 0.55];

h2.Position = [.35 .1 .55 .15];

h3.Position = [.1 .35 .15 .55];

colormap([.8 .8 1]);

秩相关系数

此构造中 X1 和 X2 之间的相关性由基础二元正态分布的相关参数 rho 决定。但是，X1 和 X2 的线性相关系数并不完全等于 rho。例如，在原来的对数正态分布中，该相关系数的闭式解为：

cor(X1,X2) = (exp(rho.*sigma.^2) - 1) ./ (exp(sigma.^2) - 1)

它严格小于 rho，除非 rho 刚好为 1。但更普遍的情形是，例如在上面的 gamma/t 分布中，X1 和 X2 之间的线性相关性很难或无法用 rho 表示，但可以通过仿真来说明存在同样的效应。

这是因为线性相关系数表达的是随机变量之间的线性相关性，当对这些随机变量应用非线性变换后，将不会保留线性相关性。在这种情况下，秩相关系数(如 Kendall 的 tau 或 Spearman 的 rho)更合适。

大体上来讲，这些秩相关衡量的是一个随机变量的大值或小值与另一个随机变量的大值或小值之间的相关程度。然而，与线性相关系数不同的是，秩相关系数只从秩的角度衡量相关性。因此，秩相关在任何单调变换下都会保留。前面所述的变换方法也会保留秩相关。因此，知道二元正态分布 Z 的秩相关就可以准确地确定最后变换的随机变量 X 的秩相关。尽管仍然需要使用 rho 来参数化基础的二元正态分布，但 Kendall tau 或 Spearman rho 在描述随机变量之间的相关性时更有用，因为它们不会随着选择的边缘分布不同而改变。

事实证明，对于二元正态分布，Kendall tau 或 Spearman rho 与线性相关系数 rho 之间存在简单的一一对应关系：

tau = (2/pi)*arcsin(rho) or rho = sin(tau*pi/2)

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

subplot(1,1,1);

rho = -1:.01:1;

tau = 2.*asin(rho)./pi;

rho_s = 6.*asin(rho./2)./pi;

plot(rho,tau,'-', rho,rho_s,'-', [-1 1],[-1 1],'k:');

axis([-1 1 -1 1]);

xlabel('rho');

ylabel('Rank correlation coefficient');

legend('Kendall''s tau', 'Spearman''s rho_s', 'location','northwest');

因此，通过为 Z1 和 Z2 之间的线性相关选择正确的 rho 参数值，很容易获得所需的 X1 和 X2 之间的秩相关系数，而不管它们的边缘分布是什么。

请注意，对于多元正态分布，Spearman 秩相关系数几乎与线性相关系数完全相同。但是，一旦我们变换为最终的随机变量，上面的情形便不再成立了。

copula

上述构造的第一步是定义一个 copula，具体来说，就是高斯 copula。二元 copula 就是两个随机变量的概率分布，每个随机变量的边缘分布都是均匀的。这两个变量之间的关系可以是完全独立、确定性相关(例如，U2 = U1)或者介于二者之间的任何状态。二元高斯 copula 族通过线性相关矩阵 Rho = [1 rho; rho 1] 进行参数化。当 rho 接近 +/- 1 时，U1 和 U2 接近线性相关；当 rho 接近 0 时，U1 和 U2 接近完全独立。

不同 rho 水平的仿真随机值散点图说明了高斯 copula 的不同可能性的范围：

n = 500;

Z = mvnrnd([0 0], [1 .8; .8 1], n);

U = normcdf(Z,0,1);

subplot(2,2,1);

plot(U(:,1),U(:,2),'.');

title('rho = 0.8');

xlabel('U1');

ylabel('U2');

Z = mvnrnd([0 0], [1 .1; .1 1], n);

U = normcdf(Z,0,1);

subplot(2,2,2);

plot(U(:,1),U(:,2),'.');

title('rho = 0.1');

xlabel('U1');

ylabel('U2');

Z = mvnrnd([0 0], [1 -.1; -.1 1], n);

U = normcdf(Z,0,1);

subplot(2,2,3);

plot(U(:,1),U(:,2),'.');

title('rho = -0.1');

xlabel('U1');

ylabel('U2');

Z = mvnrnd([0 0], [1 -.8; -.8 1], n);

U = normcdf(Z,0,1);

subplot(2,2,4);

plot(U(:,1),U(:,2),'.');

title('rho = -0.8');

xlabel('U1');

ylabel('U2');

U1 和 U2 之间的相关性完全独立于 X1 = G(U1) 和 X2 = G(U2) 的边缘分布。可以为 X1 和 X2 指定任何边缘分布，而它们的秩相关都将保持不变。这是 copula 的主要优势之一，即它们允许分别指定相关性和边缘分布。

t copula

基于二元 t 分布，使用对应的 t CDF 进行变换，可以构造一个不同的 copula 族。二元 t 分布使用线性相关矩阵 Rho 和自由度 nu 进行参数化。因此，假设多元 t 分布分别具有 1 个和 5 个自由度，我们可以把它们分别叫做 t(1) copula 和 t(5) copula。

不同 rho 水平时的仿真随机值散点图说明了 t(1) copula 的不同可能性的范围：

n = 500;

nu = 1;

T = mvtrnd([1 .8; .8 1], nu, n);

U = tcdf(T,nu);

subplot(2,2,1);

plot(U(:,1),U(:,2),'.');

title('rho = 0.8');

xlabel('U1');

ylabel('U2');

T = mvtrnd([1 .1; .1 1], nu, n);

U = tcdf(T,nu);

subplot(2,2,2);

plot(U(:,1),U(:,2),'.');

title('rho = 0.1');

xlabel('U1');

ylabel('U2');

T = mvtrnd([1 -.1; -.1 1], nu, n);

U = tcdf(T,nu);

subplot(2,2,3);

plot(U(:,1),U(:,2),'.');

title('rho = -0.1');

xlabel('U1');

ylabel('U2');

T = mvtrnd([1 -.8; -.8 1], nu, n);

U = tcdf(T,nu);

subplot(2,2,4);

plot(U(:,1),U(:,2),'.');

title('rho = -0.8');

xlabel('U1');

ylabel('U2');

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

与高斯 copula 一样，可以对 t copula 应用任何边缘分布。例如，使用具有 1 个自由度的 t copula，我们可以再次从具有 Gam(2,1) 和 t(5) 边缘的二元分布中生成随机向量：

subplot(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),'.');

axis([0 15 -10 10]);

h1 = gca;

title('1000 Simulated Dependent t and Gamma Values');

xlabel('X1 ~ Gamma(2,1)');

ylabel('X2 ~ t(5)');

subplot(2,2,4);

bar(ctr1,-n1,1);

axis([0 15 -max(n1)*1.1 0]);

axis('off');

h2 = gca;

subplot(2,2,1);

barh(ctr2,-n2,1);

axis([-max(n2)*1.1 0 -10 10]);

axis('off');

h3 = gca;

h1.Position = [0.35 0.35 0.55 0.55];

h2.Position = [.35 .1 .55 .15];

h3.Position = [.1 .35 .15 .55];

colormap([.8 .8 1]);

与之前基于高斯 copula 构造的二元 gamma/t 分布相比，这里基于 t(1) copula 构造的分布的变量之间具有相同的边缘分布和相同的秩相关，但相关性结构非常不同。这说明一个事实，即多元分布并不是仅由边缘分布或者仅由相关性定义的。在实际应用中，可以基于实际观测数据选择特定的 copula，也可以使用不同的 copula 来确定仿真结果对输入分布的敏感性。

高阶 copula

高斯 copula 和 t copula 被称为椭圆形 copula。将椭圆形 copula 泛化应用于更多的维数很容易。例如，我们可以按如下所述使用高斯 copula，来仿真具有 Gamma(2,1)、Beta(2,2) 和 t(5) 边缘分布的三元分布的数据。

subplot(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),'.');

grid on;

view([-55, 15]);

xlabel('U1');

ylabel('U2');

zlabel('U3');

请注意，对于此处使用的相关矩阵 Rho 中的每个条目，线性相关参数 rho 与 Kendall tau 等相关系数间的既定关系都成立。我们可以验证，数据的样本秩相关系数近似等于理论值。

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, 'type','Kendall')

tauSample =

1.0000 0.2655 0.1060

0.2655 1.0000 -0.6076

0.1060 -0.6076 1.0000

copula 和经验边缘分布

我们已看到，要使用 copula 来仿真相关多变量数据，我们需要指定

1) the copula family (and any shape parameters),

2) the rank correlations among variables, and

3) the marginal distributions for each variable

假设我们有两个股票收益数据集，我们想要使用分布与这两个数据集相同的输入来运行蒙特卡罗模拟。

load stockreturns

nobs = size(stocks,1);

subplot(2,1,1);

hist(stocks(:,1),10);

xlabel('X1');

ylabel('Frequency');

subplot(2,1,2);

hist(stocks(:,2),10);

xlabel('X2');

ylabel('Frequency');

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

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

这些数据集的经验逆 CDF 就是一个阶跃函数，阶跃点在 1/nobs、2/nobs 等值处。1.阶跃高度是简单的有序数据。

invCDF1 = sort(stocks(:,1));

n1 = length(stocks(:,1));

invCDF2 = sort(stocks(:,2));

n2 = length(stocks(:,2));

subplot(1,1,1);

stairs((1:nobs)/nobs, invCDF1,'b');

hold on;

stairs((1:nobs)/nobs, invCDF2,'r');

hold off;

legend('X1','X2');

xlabel('Cumulative Probability');

ylabel('X');

对于仿真，我们可能希望尝试不同的 copula 和相关性。这里，我们将使用一个二元 t(5) copula，并指定一个非常大的负相关参数。

n = 1000;

rho = -.8;

nu = 5;

T = mvtrnd([1 rho; rho 1], nu, n);

U = tcdf(T,nu);

X = [invCDF1(ceil(n1*U(:,1))) invCDF2(ceil(n2*U(:,2)))];

[n1,ctr1] = hist(X(:,1),10);

[n2,ctr2] = hist(X(:,2),10);

subplot(2,2,2);

plot(X(:,1),X(:,2),'.');

axis([-3.5 3.5 -3.5 3.5]);

h1 = gca;

title('1000 Simulated Dependent Values');

xlabel('X1');

ylabel('X2');

subplot(2,2,4);

bar(ctr1,-n1,1);

axis([-3.5 3.5 -max(n1)*1.1 0]);

axis('off');

h2 = gca;

subplot(2,2,1);

barh(ctr2,-n2,1);

axis([-max(n2)*1.1 0 -3.5 3.5]);

axis('off');

h3 = gca;

h1.Position = [0.35 0.35 0.55 0.55];

h2.Position = [.35 .1 .55 .15];

h3.Position = [.1 .35 .15 .55];

colormap([.8 .8 1]);

仿真数据的边缘直方图与原始数据的边缘直方图非常接近，而且随着我们仿真的值对组增多，二者会变得完全相同。请注意，这些值是从原始数据中抽取的，由于每个数据集中只有 100 个观测值，因此仿真数据看上去有些“离散”。要解决此问题，可在最后仿真的值中增加少量可能呈正态分布的随机变异。这就相当于使用经过平滑的经验逆 CDF。