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

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

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

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

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

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

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. не обязательно имеют одинаковые предельные распределения. Преобразование определяется как

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

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

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

Оказывается, что для двумерной нормы существует простое 1-1 отображение между тау Кендалла или ро Спирмена, и коэффициентом линейной корреляции ро:

  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 приближается к нулю.

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

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

Графики рассеяния некоторых моделируемых случайных значений для различных уровней ро иллюстрируют диапазон различных возможностей для 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, как и гауссова копула. Ранговая корреляция тау или rho_s между компонентами в t-копуле также является той же функцией ро, что и для гауссова. Однако, как показывают эти графики, 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 и, например, тау Кендалла, сохраняется для каждой записи в корреляционной матрице 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.