Документация

Определение зависимости между входами моделирования

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

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

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

n = 1000;

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

SigmaInd = 2×2

    0.2500         0
         0    0.2500

rng('default');  % For reproducibility
ZInd = mvnrnd([0 0],SigmaInd,n);
XInd = exp(ZInd);

plot(XInd(:,1),XInd(:,2),'.')
axis([0 5 0 5])
axis equal
xlabel('X1')
ylabel('X2')

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

SigmaDep = 2×2

    0.2500    0.1750
    0.1750    0.2500

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

plot(XDep(:,1),XDep(:,2),'.')
axis([0 5 0 5])
axis equal
xlabel('X1')
ylabel('X2')

Построение зависимых двухмерных распределений

Хотя конструкция, описанная в предыдущем разделе, создает двумерную логнормальную структуру, которая является простой, она служит для иллюстрации способа, который более широко применим.

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

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

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

n = 1000; 
rng default % for reproducibility
z = normrnd(0,1,n,1); % generate standard normal data

histogram(z,-3.75:.5:3.75,'FaceColor',[.8 .8 1]) % plot the histogram of data
xlim([-4 4])
title('1000 Simulated N(0,1) Random Values')
xlabel('Z')
ylabel('Frequency')

u = normcdf(z);  % compute the cdf values of the sample data

figure
histogram(u,.05:.1:.95,'FaceColor',[.8 .8 1]) % plot the histogram of the cdf values
title('1000 Simulated N(0,1) Values Transformed to Unif(0,1)')
xlabel('U')
ylabel('Frequency')

Заимствование из теории генерации одномерных случайных чисел, применяя обратную cdf любого распределения, F, к случайной величине Unif (0,1) приводит к случайной величине, распределение которой точно равно F (см. Методы инверсии). Доказательство по существу противоположно предыдущему доказательству для форвардного дела. Другая гистограмма иллюстрирует преобразование в гамма-распределение:

x = gaminv(u,2,1); % transform to gamma values

figure
histogram(x,.25:.5:9.75,'FaceColor',[.8 .8 1]) % plot the histogram of gamma values
title('1000 Simulated N(0,1) Values Transformed to Gamma(2,1)')
xlabel('X')
ylabel('Frequency')

Это двухшаговое преобразование можно применить к каждой переменной стандартной двухмерной нормы, создавая зависимые случайные величины с произвольными предельными распределениями. Поскольку преобразование работает для каждого компонента отдельно, две результирующие случайные величины не обязательно должны иметь одинаковые предельные распределения. Преобразование определяется как:

где ${\mathit{}}_{\mathrm{G1}}$ и ${\mathit{}}_{\mathrm{G2}}$ - обратные cdfs двух, возможно, различных распределений. Например, следующее генерирует случайные векторы из двумерного распределения с ${\mathit{}}_{\mathrm{маргиналами\; t5}}$ и 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)];

% draw the scatter plot of data with histograms 
figure
scatterhist(X(:,1),X(:,2),'Direction','out')

Этот график имеет гистограммы рядом с графиком рассеяния, чтобы показать как предельные распределения, так и зависимость.

Использование коэффициентов ранговой корреляции

Параметр корреляции, ρ, основного двумерного нормального определяет зависимость между X1 и X2 в этой конструкции. Однако линейная корреляция X1 и X2 не является start. Например, в исходном логнормальном случае закрытой формой для этой корреляции является:

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

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

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

Для двумерного нормального есть простое непосредственное отображение между τ Кендалла или ρ Копьеносца и линейным коэффициентом корреляции ρ:

На следующем графике показана взаимосвязь.

rho = -1:.01:1;
tau = 2.*asin(rho)./pi;
rho_s = 6.*asin(rho./2)./pi;

plot(rho,tau,'b-','LineWidth',2)
hold on
plot(rho,rho_s,'g-','LineWidth',2)
plot([-1 1],[-1 1],'k:','LineWidth',2)
axis([-1 1 -1 1])
xlabel('rho')
ylabel('Rank correlation coefficient')
legend('Kendall''s {\it\tau}', ...
       'Spearman''s {\it\rho_s}', ... 
       'location','NW')

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

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

Использование двухмерных копул

Первый этап конструкции, описанной в предыдущем разделе, определяет то, что известно как двумерная гауссова копула. Копула - многомерное распределение вероятностей, где каждая случайная величина имеет равномерное предельное распределение на единичном интервале [0,1]. Эти переменные могут быть полностью независимыми, детерминированно связанными (например, U2 = U1), или что-либо между ними. Из-за возможности зависимости между переменными можно использовать копулу для построения нового многомерного распределения для зависимых переменных. Путем преобразования каждой из переменных в копуле по отдельности с использованием метода инверсии, возможно с использованием различных cdfs, результирующее распределение может иметь произвольные предельные распределения. Такие многомерные распределения часто полезны при моделировании, когда известно, что различные случайные входные данные не независимы друг от друга.

Вычисляются функции набора инструментов для статистического и машинного обучения:

Например, используйте copularnd функция для создания диаграмм рассеяния случайных значений из двухмерной гауссовой копулы для различных уровней, для иллюстрации диапазона различных структур зависимости. Семейство двумерных гауссовых копул параметризуется матрицей линейной корреляции:

U1 и U2 приближаться к линейной зависимости по мере того, как startприближается к ± 1, и приближаться к полной независимости по мере того, как startприближается к нулю:

n = 500;

rng('default') % for reproducibility
U = copularnd('Gaussian',[1 .8; .8 1],n);
subplot(2,2,1)
plot(U(:,1),U(:,2),'.')
title('{\it\rho} = 0.8')
xlabel('U1')
ylabel('U2')

U = copularnd('Gaussian',[1 .1; .1 1],n);
subplot(2,2,2)
plot(U(:,1),U(:,2),'.')
title('{\it\rho} = 0.1')
xlabel('U1')
ylabel('U2')

U = copularnd('Gaussian',[1 -.1; -.1 1],n);
subplot(2,2,3)
plot(U(:,1),U(:,2),'.')
title('{\it\rho} = -0.1')
xlabel('U1')
ylabel('U2')

U = copularnd('Gaussian',[1 -.8; -.8 1],n);
subplot(2,2,4)
plot(U(:,1),U(:,2),'.')
title('{\it\rho} = -0.8')
xlabel('U1')
ylabel('U2')

Зависимость между U1 и U2 полностью отделен от предельных распределений X1 = G(U1) и X2 = G(U2). X1 и X2 могут быть даны любые предельные распределения, и все еще имеют одинаковую ранговую корреляцию. Это одно из главных обращений копул - они допускают эту отдельную спецификацию зависимости и предельного распределения. Можно также вычислить pdf (copulapdf) и cdf (copulacdf) для копулы. Например, на этих графиках показаны файлы pdf и cdf для

u1 = linspace(1e-3,1-1e-3,50);
u2 = linspace(1e-3,1-1e-3,50);
[U1,U2] = meshgrid(u1,u2);
Rho = [1 .8; .8 1];
f = copulapdf('t',[U1(:) U2(:)],Rho,5);
f = reshape(f,size(U1));

figure
surf(u1,u2,log(f),'FaceColor','interp','EdgeColor','none')
view([-15,20])
xlabel('U1')
ylabel('U2')
zlabel('Probability Density')

u1 = linspace(1e-3,1-1e-3,50);
u2 = linspace(1e-3,1-1e-3,50);
[U1,U2] = meshgrid(u1,u2);
F = copulacdf('t',[U1(:) U2(:)],Rho,5);
F = reshape(F,size(U1));

figure()
surf(u1,u2,F,'FaceColor','interp','EdgeColor','none')
view([-15,20])
xlabel('U1')
ylabel('U2')
zlabel('Cumulative Probability')

Другое семейство копул может быть построено, начиная с двухмерного распределения t и трансформируя с использованием соответствующего t cdf. Двумерное распределение t параметризуется с помощью P, линейной корреляционной матрицы, и start, степеней свободы. Так, например, можно говорить о $${}_{\mathrm{}}$$$${}_{\mathrm{}}$$копуле t1 или t5, основанной на многомерной t с одной и пятью степенями свободы соответственно.

Так же, как и для гауссовых копул, для вычисления t копул используются функции набора инструментов статистики и машинного обучения:

Например, используйте copularnd функция для создания диаграмм рассеяния случайных значений из двухмерной $${}_{\mathrm{}}$$t1-копулы для различных уровней start, для иллюстрации диапазона различных структур зависимости:

n = 500;
nu = 1;

rng('default') % for reproducibility
U = copularnd('t',[1 .8; .8 1],nu,n);
subplot(2,2,1)
plot(U(:,1),U(:,2),'.')
title('{\it\rho} = 0.8')
xlabel('U1')
ylabel('U2')

U = copularnd('t',[1 .1; .1 1],nu,n);
subplot(2,2,2)
plot(U(:,1),U(:,2),'.')
title('{\it\rho} = 0.1')
xlabel('U1')
ylabel('U2')

U = copularnd('t',[1 -.1; -.1 1],nu,n);
subplot(2,2,3)
plot(U(:,1),U(:,2),'.')
title('{\it\rho} = -0.1')
xlabel('U1')
ylabel('U2')

U = copularnd('t',[1 -.8; -.8 1],nu, n);
subplot(2,2,4)
plot(U(:,1),U(:,2),'.')
title('{\it\rho} = -0.8')
xlabel('U1')
ylabel('U2')

T-копула имеет однородные предельные распределения для U1 и U2, как и гауссова копула. Корреляция разряда τ или $${}_{\mathrm{\rho s}}$$ между компонентами в t связке является также той же функцией ρ что касается Гауссовского. Однако, как показывают эти графики, $${}_{\mathrm{}}$$копула t1 немного отличается от гауссовой копулы, даже когда их компоненты имеют одинаковую ранговую корреляцию. Разница в их структуре зависимости. Не удивительно, как параметр степеней свободы $$\nu$$ сделан больше, $${}_{\mathrm{t\nu }}$$ связка приближается к соответствующей Гауссовской связке.

Как и в случае гауссовой копулы, любые предельные распределения могут быть наложены на t копулу. Например, используя t-копулу с 1 степенью свободы, можно снова создать случайные векторы из двухмерного распределения с Гамма (2,1) и $${}_{\mathrm{}}$$t5 маргиналами, используяcopularnd:

n = 1000;
rho = .7;
nu = 1;

rng('default') % for reproducibility
U = copularnd('t',[1 rho; rho 1],nu,n);
X = [gaminv(U(:,1),2,1) tinv(U(:,2),5)];

figure
scatterhist(X(:,1),X(:,2),'Direction','out')

По сравнению с двухмерным распределением Гамма/т, построенным ранее, которое было основано на гауссовой копуле, распределение, построенное здесь, основанное на $${}_{\mathrm{t1}}$$ копуле, имеет те же предельные распределения и ту же ранговую корреляцию между переменными, но очень разную структуру зависимости. Это иллюстрирует тот факт, что многомерные распределения не однозначно определяются их предельными распределениями или их корреляциями. Выбор конкретной копулы в заявке может быть основан на фактических наблюдаемых данных, или различные копулы могут быть использованы в качестве способа определения чувствительности результатов моделирования к входному распределению.

Копулы высших размеров

Гауссовы и t-копулы известны как эллиптические копулы. Легко обобщить эллиптические копулы до большего числа измерений. Например, смоделировать данные из трехмерного распределения с помощью Гамма (2,1), Бета (2,2) $${и}_{\mathrm{}}$$t5 маргиналов, используя гауссову копулу иcopularnd, следующим образом:

n = 1000;
Rho = [1 .4 .2; .4 1 -.8; .2 -.8 1];
rng('default') % for reproducibility
U = copularnd('Gaussian',Rho,n);
X = [gaminv(U(:,1),2,1) betainv(U(:,2),2,2) tinv(U(:,3),5)];

Постройте график данных.

subplot(1,1,1)
plot3(X(:,1),X(:,2),X(:,3),'.')
grid on
view([-55, 15])
xlabel('X1')
ylabel('X2')
zlabel('X3')

Обратите внимание, что соотношение между линейным корреляционным параметром (, например, по Кендаллу) для каждой записи в корреляционной матрице (P), используемой здесь, сохраняется. Можно проверить, что ранговые корреляции выборки данных приблизительно равны теоретическим значениям:

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

tauTheoretical = 3×3

    1.0000    0.2620    0.1282
    0.2620    1.0000   -0.5903
    0.1282   -0.5903    1.0000

tauSample = corr(X,'type','Kendall')

tauSample = 3×3

    1.0000    0.2581    0.1414
    0.2581    1.0000   -0.5790
    0.1414   -0.5790    1.0000

Архимедовый копулас

Функции Statistics и Machine Learning Toolbox доступны для трех двухмерных семейств Archimedean copula:

Это однопараметрические семейства, которые определяются непосредственно в терминах cdfs, а не конструктивно с использованием стандартного многомерного распределения.

Чтобы сравнить эти три архимедовые копулы с гауссовой и t-бивариатной копулами, сначала используйте copulastat чтобы найти ранговую корреляцию для гауссовой или t-копулы с параметром линейной корреляции 0,8, а затем использовать copulaparam чтобы найти параметр копулы Клейтона, который соответствует этой ранговой корреляции:

tau = copulastat('Gaussian',.8 ,'type','kendall')

tau = 0.5903

alpha = copulaparam('Clayton',tau,'type','kendall')

alpha = 2.8820

Наконец, постройте график случайной выборки из пары Клейтона с copularnd. Повторите ту же процедуру для копул Фрэнка и Гумбеля:

n = 500;

U = copularnd('Clayton',alpha,n);
subplot(3,1,1)
plot(U(:,1),U(:,2),'.');
title(['Clayton Copula, {\it\alpha} = ',sprintf('%0.2f',alpha)])
xlabel('U1')
ylabel('U2')

alpha = copulaparam('Frank',tau,'type','kendall');
U = copularnd('Frank',alpha,n);
subplot(3,1,2)
plot(U(:,1),U(:,2),'.')
title(['Frank Copula, {\it\alpha} = ',sprintf('%0.2f',alpha)])
xlabel('U1')
ylabel('U2')

alpha = copulaparam('Gumbel',tau,'type','kendall');
U = copularnd('Gumbel',alpha,n);
subplot(3,1,3)
plot(U(:,1),U(:,2),'.')
title(['Gumbel Copula, {\it\alpha} = ',sprintf('%0.2f',alpha)])
xlabel('U1')
ylabel('U2')

Моделирование зависимых многомерных данных с помощью Copulas

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

Предположим, что у вас есть данные возврата для двух акций и вы хотите запустить моделирование Монте-Карло с входными данными, которые следуют тем же распределениям, что и данные:

load stockreturns
nobs = size(stocks,1);

subplot(2,1,1)
histogram(stocks(:,1),10,'FaceColor',[.8 .8 1])
xlim([-3.5 3.5])
xlabel('X1')
ylabel('Frequency')

subplot(2,1,2)
histogram(stocks(:,2),10,'FaceColor',[.8 .8 1])
xlim([-3.5 3.5])
xlabel('X2')
ylabel('Frequency')

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

Простейшей непараметрической моделью является эмпирическая cdf, вычисленная с помощью ecdf функция. Для дискретного маргинального распределения это уместно. Однако для непрерывного распределения используйте модель, которая является более гладкой, чем ступенчатая функция, вычисленная ecdf. Один из способов сделать это - оценить эмпирический cdf и интерполировать между серединами шагов с кусочно-линейной функцией. Другой способ - использовать сглаживание ядра с помощью ksdensity. Например, сравните эмпирический cdf со сглаженной ядром оценкой cdf для первой переменной:

[Fi,xi] = ecdf(stocks(:,1));

figure()
stairs(xi,Fi,'b','LineWidth',2)
hold on

Fi_sm = ksdensity(stocks(:,1),xi,'function','cdf','width',.15);

plot(xi,Fi_sm,'r-','LineWidth',1.5)
xlabel('X1')
ylabel('Cumulative Probability')
legend('Empirical','Smoothed','Location','NW')
grid on

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

nu = 5;
tau = corr(stocks(:,1),stocks(:,2),'type','kendall')

tau = 0.5180

Найти соответствующий параметр линейной корреляции для t-пары с помощью copulaparam.

rho = copulaparam('t', tau, nu, 'type','kendall')

rho = 0.7268

Далее используйте copularnd для генерации случайных значений из t copula и преобразования с использованием непараметрических обратных cdfs. ksdensity функция позволяет сделать оценку распределения ядра и оценить обратное cdf в точках копулы все за один шаг:

n = 1000;
U = copularnd('t',[1 rho; rho 1],nu,n);
X1 = ksdensity(stocks(:,1),U(:,1),...
               'function','icdf','width',.15);
X2 = ksdensity(stocks(:,2),U(:,2),...
               'function','icdf','width',.15);

В качестве альтернативы, если имеется большое количество данных или необходимо смоделировать несколько наборов значений, может быть более эффективным вычислить обратное cdf по сетке значений в интервале. (0,1) и использовать интерполяцию для его оценки в точках совокупностей:

p = linspace(0.00001,0.99999,1000);
G1 = ksdensity(stocks(:,1),p,'function','icdf','width',0.15);
X1 = interp1(p,G1,U(:,1),'spline');
G2 = ksdensity(stocks(:,2),p,'function','icdf','width',0.15);
X2 = interp1(p,G2,U(:,2),'spline');

scatterhist(X1,X2,'Direction','out')

Краевые гистограммы моделируемых данных являются сглаженной версией гистограмм для исходных данных. Величина сглаживания контролируется вводом полосы пропускания в ksdensity.

Подбор копул к данным

В этом примере показано, как использовать copulafit для калибровки копул с данными. Создание данных Xsim с распределением «точно как» (в терминах предельных распределений и корреляций) распределение данных в матрице X, нужно подогнать маргинальные распределения к столбцам X, используйте соответствующие функции cdf для преобразования X кому U, так что U имеет значения от 0 до 1, используйте copulafit чтобы вписать копулу в U, создать новые данные Usim из копулы и использовать соответствующие обратные cdf-функции для преобразования Usim кому Xsim .

Загрузите и запишите смоделированные данные возврата запаса.

load stockreturns
x = stocks(:,1);
y = stocks(:,2);

scatterhist(x,y,'Direction','out')

Преобразование данных в шкалу копулы (единичный квадрат) с использованием оценщика ядра кумулятивной функции распределения.

u = ksdensity(x,x,'function','cdf');
v = ksdensity(y,y,'function','cdf');

scatterhist(u,v,'Direction','out')
xlabel('u')
ylabel('v')

Подойди к копуле.

[Rho,nu] = copulafit('t',[u v],'Method','ApproximateML')

Rho = 2×2

    1.0000    0.7220
    0.7220    1.0000

nu = 3.2727e+06

Создайте случайную выборку из t-пары.

r = copularnd('t',Rho,nu,1000);
u1 = r(:,1);
v1 = r(:,2);

scatterhist(u1,v1,'Direction','out')
xlabel('u')
ylabel('v')
set(get(gca,'children'),'marker','.')

Преобразуйте случайную выборку обратно в исходный масштаб данных.

x1 = ksdensity(x,u1,'function','icdf'); 
y1 = ksdensity(y,v1,'function','icdf');

scatterhist(x1,y1,'Direction','out')
set(get(gca,'children'),'marker','.')

Как показано в примере, копулы естественным образом интегрируются с другими функциями фитинга распределения.