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

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

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

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

Симуляция независимых логарифмически нормальных случайных переменных тривиальна. Самый простой путь состоял бы в том, чтобы использовать lognrnd функция. Здесь, мы будем использовать mvnrnd функция, чтобы сгенерировать n пары независимых нормальных случайных переменных, и затем exponentiate их. Заметьте, что ковариационная матрица, используемая здесь, является диагональной, i.e., независимость между столбцами 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.'s также легко сгенерировать, с помощью ковариационной матрицы с ненулевыми недиагональными условиями.

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

Двумерное логарифмически нормальное распределение является простым решением в этом случае, и конечно легко делает вывод к более высоким размерностям и случаям, где предельные распределения являются различным lognormals. Другие многомерные распределения также существуют, например, многомерный t и распределения Дирихле используются, чтобы симулировать зависимого t и бета случайные переменные, соответственно. Но список простых многомерных распределений не длинен, и они только применяются в случаях, где marginals - все в том же семействе (или даже тех же самых распределениях). Это может быть действительным ограничением во многих ситуациях.

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

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

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

По определению применение нормального 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.'s с произвольными предельными распределениями. Поскольку преобразование работает над каждым компонентом отдельно, потребность двух получившихся r.v. даже не имеют те же предельные распределения. Преобразование задано как

  Z = [Z1 Z2] ~ N([0 0],[1 rho; rho 1])
  U = [PHI(Z1) PHI(Z2)]
  X = [G1(U1) G2(U2)]

где G1 и G2 является обратный CDFS двух возможно различных распределений. Например, мы можем сгенерировать случайные векторы от двумерного распределения с t (5) и Гамма (2,1) marginals.

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

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

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

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

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

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

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

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

Связки

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

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

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 распределение параметрируется с Ро, матрицей линейной корреляции, и ню, степенями свободы. Таким образом, например, мы можем говорить о 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, как Гауссова связка делает. Порядковая корреляция tau или rho_s между компонентами в t связке являются также той же функцией ро что касается Гауссова. Однако как эти графики демонстрируют, t (1), связка отличается вполне немного от Гауссовой связки, даже когда их компоненты имеют ту же порядковую корреляцию. Различие находится в их структуре зависимости. Не удивительно, когда ню параметра степеней свободы сделан больше, t (nu), связка приближается к соответствующей Гауссовой связке.

Как с Гауссовой связкой, любые предельные распределения могут быть наложены по t связке. Например, с помощью t связки с 1 степенью свободы, мы можем снова сгенерировать случайные векторы от двумерного распределения с Ножкой (2,1) и t (5) marginals:

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

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

Связки высшего порядка

Гауссовы и t связки известны как эллиптические связки. Легко обобщить эллиптические связки к более высокому количеству размерностей. Например, мы можем симулировать данные из trivariate распределения с Гаммой (2,1), Бета (2,2) и t (5) marginals использование Гауссовой связки можно следующим образом.

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

Заметьте, что отношение между ро параметра линейной корреляции и, например, 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

Связки и эмпирические предельные распределения

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

  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.