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

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

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

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

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

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

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

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

Используя коэффициенты порядковой корреляции

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

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

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

Примерно говоря, эти порядковые корреляции измеряют степень, до которой большие или маленькие значения одной случайной переменной сопоставляют с большими или маленькими значениями другого. Однако различающийся коэффициент линейной корреляции, они измеряют ассоциацию только в терминах рангов. Как следствие порядковая корреляция сохраняется при любом монотонном преобразовании. В частности, метод преобразования только описал, сохраняет порядковую корреляцию. Поэтому знание порядковой корреляции двумерного нормального 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, независимо от их предельных распределений, путем выбора правильного ρ значения параметров для линейной корреляции между 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 для ρ =.8:

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

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

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

Более высокие связки размерности

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

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