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

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

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

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

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

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 =

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

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

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

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

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

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

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

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

Для двухмерной нормы существует простое взаимно-однозначное отображение между Kendall's

Следующий график показывает отношение.

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

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

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

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

Функции Statistics and Machine Learning Toolbox вычисляют:

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

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

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, линейной корреляционной матрицей, и Таким образом, например, можно говорить о $$t_1$$ или a $$t_5$$ копула, основанная на многомерной t с одной и пятью степенями свободы, соответственно.

Так же, как и для Гауссовых копул, функции Statistics and Machine Learning Toolbox для t копул вычисляют:

Для примера используйте copularnd функция для создания графиков поля точек случайных значений из двухмерной матрицы $$t_1$$ копула для различных уровней и для иллюстрации области значений различных зависимых структур:

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, так же, как это делает Гауссов копула. Ранговая корреляция $${\rho }_s$$ между компонентами в t копуле также является той же функцией, что и для Гауссова. Однако, как показывают эти графики, $$t_1$$ копула довольно немного отличается от гауссовой копулы, даже когда их компоненты имеют одинаковую ранговую корреляцию. Это различие находится в их структуре зависимости. Неудивительно, что как параметр степеней свободы $$\nu$$ сделано больше, а $$t_{\nu }$$ копула приближается к соответствующей Гауссовой копуле.

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

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

Копулы Высшей Размерности

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

Заметьте, что отношение между параметром линейной корреляции, и, например, Kendall's, удерживается для каждой записи в матрице корреляции 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 and Machine Learning Toolbox доступны для трех двухмерных семейств архимедовых копул:

Это однопараметрические семейства, которые заданы непосредственно в терминах их 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')

Симуляция зависимых многомерных данных с использованием копул

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

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

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

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

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-копулы и преобразовать с помощью непараметрического обратного cdfs. The 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')

Подгонка t копулы.

[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','.')

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