Зависимость между входами симуляции

Одним из проекта решений для симуляции Монте-Карло является выбор распределений вероятностей для случайных входных параметров. Выбор распределения для каждой отдельной переменной часто прост, но решение, какие зависимости должны существовать между входами, может не быть. В идеале входные данные для симуляции должны отражать то, что известно о зависимости среди моделируемых действительных величин. Однако может быть мало или нет информации, на которой можно основывать любую зависимость в симуляции, и в таких случаях рекомендуется экспериментировать с различными возможностями, порядком определить чувствительность модели.

Однако может быть трудно фактически сгенерировать случайные входы с зависимостью, когда они имеют распределения, которые не из стандартного многомерного распределения. Кроме того, некоторые стандартные многомерные распределения могут моделировать только очень ограниченные типы зависимости. Всегда можно сделать входы независимыми, и, хотя это простой выбор, это не всегда разумно и может привести к неправильным выводам.

Например, симуляция финансового риска в Монте-Карло может иметь случайные входы, которые представляют различные источники страховых потерь. Эти входы могут быть смоделированы как логнормальные случайные переменные. Разумный вопрос, который нужно задать, заключается в том, как зависимость между этими двумя входами влияет на результаты симуляции. Действительно, из реальных данных может быть известно, что одни и те же случайные условия влияют на оба источника, и игнорирование этого в симуляции может привести к неправильным выводам.

Симуляция независимых логнормальных случайных переменных тривиальна. Самым простым способом было бы использование 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');

Зависимые двухмерные lognormal r.v. также легко сгенерировать, используя матрицу ковариации с ненулевыми off-диагональными терминами.

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 и Дирихле распределения используются для моделирования зависимых t и бета случайных переменных, соответственно. Но список простых многомерных распределений невелик, и они применяются только в тех случаях, когда маргиналы все находятся в одном семействе (или даже точно таких же распределениях). Это может быть реальным ограничением во многих ситуациях.

Более общий метод построения зависимых двухмерных распределений

Хотя вышеописанная конструкция, которая создает двухмерное логнормальное, проста, она служит для иллюстрации способа, который более широко применим. Во-первых, мы генерируем пары значений из двухмерного нормального распределения. Между этими двумя переменными существует статистическая зависимость, и каждая из них имеет нормальное предельное распределение. Затем преобразование (экспоненциальная функция) применяется отдельно к каждой переменной, меняя маргинальные распределения на логнормальные. Преобразованные переменные все еще имеют статистическую зависимость.

Если можно найти подходящее преобразование, этот метод может быть обобщен, чтобы создать зависимые двухмерные случайные векторы с другими маргинальными распределениями. На самом деле, общий метод построения такого преобразования действительно существует, хотя и не так просто, как просто экспоненция.

По определению применение нормальной CDF (обозначенной здесь PHI) к стандартной нормальной случайной переменной приводит к r.v. это равномерно по интервалу [0, 1]. Чтобы увидеть это, если Z имеет стандартное нормальное распределение, то CDF U = PHI (Z) является

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

и это CDF U (0,1) r.v. Гистограммы некоторых моделируемых нормальных и преобразованных значений демонстрируют этот факт.

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

Теперь заимствование из теории одномерной генерации случайных чисел, применение обратного CDF любого распределения F к U (0,1) случайной переменной приводит к r.v. чье распределение в точности F. Это известно как Метод Инверсии. Доказательство, по существу, противоположно вышеописанному доказательству для переднего случая. Другая гистограмма иллюстрирует преобразование в гамма- распределение.

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

Это двухэтапное преобразование может быть применено к каждой переменной стандартной двухмерной нормы, создавая зависимые r.v с произвольными маргинальными распределениями. Поскольку преобразование работает с каждым компонентом отдельно, два получившихся r.v.'s не должны даже иметь одинаковые маргинальные распределения. Преобразование определяется как

  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) и Гамма (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. Например, в исходном случае lognormal существует закрытая форма для этой корреляции:

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

который строго меньше rho, если rho не является именно одним. В более общих случаях, таких как конструкция Гамма/т выше, линейная корреляция между X1 и X2 трудно или невозможно выразить с точки зрения rho, но можно использовать симуляции, чтобы показать, что тот же эффект происходит.

Это потому, что линейный коэффициент корреляции выражает линейную зависимость между r.v.'s, и когда нелинейные преобразования применяются к этим r.v.'s, линейная корреляция не сохраняется. Вместо этого более подходящим является ранговый коэффициент корреляции, такой как тау Кендалла или ро Спирмана.

Грубо говоря, эти ранговые корреляции измеряют степень, до которой большие или маленькие значения одного r.v. связать с большими или маленькими значениями другого. Однако, в отличие от коэффициента линейной корреляции, они измеряют ассоциацию только в терминах рангов. Как следствие, ранг корреляция сохраняется при любом монотонном преобразовании. В частности, только что описанный способ преобразования сохраняет ранг корреляцию. Поэтому знание ранговой корреляции двухмерной нормы Z в точности определяет ранговую корреляцию конечного преобразованного r.v.'s X. В то время как rho все еще необходимо, чтобы параметризовать базовую двухмерную нормаль, тау Кендалла или rho Спирмана более полезны в описании зависимости между r.v.'s,

Оказывается, что для двухмерной нормы существует простое отображение 1-1 между тау Кендалла или ро Спирмана и коэффициентом линейной корреляции 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');

Таким образом, легко создать желаемую ранговую корреляцию между X1 и X2, независимо от их маргинальных распределений, путем выбора правильного значения параметров rho для линейной корреляции между Z1 и Z2.

Заметьте, что для многомерного нормального распределения ранговая корреляция Спирмана почти идентична линейной корреляции. Однако это не так, если мы преобразуемся в конечные случайные переменные.

Связки

Первый этап описанной выше конструкции определяет, что известно как копула, в частности, гауссову копулу. Двухмерная копула является просто распределение вероятностей на двух случайных переменных, каждое из маргинальных распределений которых равномерно. Эти две переменные могут быть полностью независимыми, детерминированно связанными (например, U2 = U1) или чем-либо средним. Семейство двухмерных гауссовых копул параметризовано Rho = [1 rho; rho 1], матрица линейной корреляции. U1 и U2 подходят к линейной зависимости, когда rho приближается к +/- 1, и подходят к полной независимости, когда rho приближается к нулю.

Графики поля точек некоторых моделируемых случайных значений для различных уровней rho иллюстрируют область значений различных возможностей для гауссовых копул:

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 могут быть даны любые маргинальные распределения, и все еще имеют ту же ранговую корреляцию. Это один из основных призывов копул - они допускают эту отдельную спецификацию зависимости и маргинального распределения.

t Копулас

Другое семейство копул может быть построено путем начала с двухмерной t- распределения и преобразования с использованием соответствующего t CDF. Двухмерное распределение t параметризовано с помощью Rho, линейной корреляционной матрицы и nu, степеней свободы. Таким образом, например, можно говорить о t (1) или t (5) копуле, основанной на многомерной t с одной и пятью степенями свободы, соответственно.

Графики поля точек некоторых моделируемых случайных значений для различных уровней rho иллюстрируют область значений различных возможностей для t (1) копул:

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

t копула имеет равномерные маргинальные распределения для U1 и U2, так же как Гауссова копула. Ранговая корреляция tau или rho_s между компонентами в t копуле также является такой же функцией rho, как и для Гауссова. Однако, как показывают эти графики, t (1) копула довольно немного отличается от гауссовой копулы, даже когда их компоненты имеют одинаковую ранговую корреляцию. Это различие находится в их структуре зависимости. Неудивительно, что, когда параметр nu степеней свободы сделан большим, t (nu) копула приближается к соответствующей Гауссовой копуле.

Как и в случае Гауссовой копулы, любые маргинальные распределения могут быть наложены на t копулу. Для примера, используя t-копулу с 1 степенью свободы, мы можем снова сгенерировать случайные векторы из двухмерного распределения с маргиналами 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]);

По сравнению с двухмерным распределением Гамма/т, построенным ранее, которое было основано на гауссовой копуле, построенное здесь распределение, основанное на t (1) копуле, имеет те же маргинальные распределения и ту же ранговую корреляцию между переменными, но очень разную структуру зависимости. Это иллюстрирует тот факт, что многомерные распределения не уникально заданы их маргинальными распределениями или их корреляциями. Выбор конкретной копулы в заявке может основываться на фактических наблюдаемых данных, или различные копулы могут использоваться в качестве способа определения чувствительности результатов симуляции к распределению входа.

Копулы более высокого порядка

Гауссовы и t копулы известны как эллиптические копулы. Легко обобщить эллиптические копулы к большему количеству размерностей. Для примера мы можем симулировать данные из трехмерного распределения с Гаммой (2,1), Бета (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 и, например, 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, 'type','Kendall')

tauSample =

    1.0000    0.2655    0.1060
    0.2655    1.0000   -0.6076
    0.1060   -0.6076    1.0000

Копулы и эмпирические маргинальные распределения

Чтобы симулировать зависимые многомерные данные с помощью копулы, мы увидели, что нам нужно задать

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

Для симуляции, мы, возможно, захотим экспериментировать с различными копулами и корреляциями. Здесь мы будем использовать двухмерную t (5) копулу с довольно большим параметром отрицательной корреляции.

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.