Симулируйте стохастический процесс формулой Фейнмана-Каца

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

Включенные шаги:

  1. Задайте параметры модели Используя стохастические дифференциальные уравнения

  2. Примените правило ITO

  3. Решите уравнение Фейнмана-Каца

  4. Вычислите ожидаемое время, чтобы продать актив

1. Задайте параметры модели Используя стохастические дифференциальные уравнения

Модель за цену актива X(t) заданный во временном интервале [0,T] стохастический процесс, заданный стохастическим дифференциальным уравнением формы dX=μ(t,X)dt+σ(t,X)dB(t), где B(t) Винеровский процесс модульным параметром отклонения.

Создайте символьную переменную T и следующие символьные функции:

  • r(t) непрерывная функция, представляющая точечную процентную ставку. Этот уровень определяет коэффициент дисконтирования для итоговой выплаты в то время T.

  • u(t,x) ожидаемое значение обесцененной будущей цены, вычисленной как X(T)exp(-tTr(t)dt) при условии X(t) = x.

  • μ(t,X) и σ(t,X) дрейф и диффузия стохастического процесса X(t).

syms mu(t, X) sigma(t, X) r(t, X) u(t, X) T

Согласно теореме Фейнмана-Каца, u удовлетворяет дифференциальному уравнению с частными производными ut+σ222ux2+μux-ur=0 с итоговым условием во время T.

eq = diff(u, t) + sigma^2*diff(u, X, X)/2 + mu*diff(u, X) - u*r;

Числовой решатель pdepe работает с начальными условиями. Чтобы преобразовать итоговое условие в начальное условие, примените реверсирование времени установкой v(t, X) = u(T - t, X).

syms v(t, X)
eq2 = subs(eq, {u, t}, {v(T - t, X), T - t});
eq2 = feval(symengine, 'rewrite', eq2, 'diff')
eq2 = 

σ(T-t,X)22X2 v(t,X)2+μ(T-t,X)X v(t,X)-t v(t,X)-v(t,X)r(T-t,X)(сигма (T - t, X) ^2*diff (v (t, X), X, 2))/2 + mu (T - t, X) *diff (v (t, X), X) - diff (v (t, X), t) - v (t, X) *r (T - t, X)

Решатель pdepe требует дифференциального уравнения с частными производными в следующей форме. Здесь коэффициенты cF, и s функции xTV, и v/x.

cvt=fx+s

Смочь решить уравнение eq2 с pdepe решатель, карта eq2 к этой форме. Во-первых, извлеките коэффициенты и условия eq2.

syms dvt dvXX
eq3 = subs(eq2, {diff(v, t), diff(v, X, X)}, {dvt, dvXX});
[k,terms] = coeffs(eq3, [dvt, dvXX])
k = 

(-1σ(T-t,X)22μ(T-t,X)X v(t,X)-v(t,X)r(T-t,X))[-sym (1), сигма (T - t, X) ^2/2, mu (T - t, X) *diff (v (t, X), X) - v (t, X) *r (T - t, X)]

terms = (dvtdvXX1)[dvt, dvXX, sym (1)]

Теперь можно переписать eq2 как k(1)*terms(1) + k(2)*terms(2) + k(3)*terms(3) = 0. Перемещая термин с производной времени к левой стороне уравнения, вы добираетесь

vt=k(2)2vX2+k(3)

Вручную интегрируйте правую сторону частями. Результат

X(k(2)vX)+k(3)-vXk(2)X

Поэтому записать eq2 в форме, подходящей для pdepe, используйте следующие параметры:

c = 1;
f = k(2) * diff(v, X);
s = k(3) - diff(v, X) * diff(k(2), X);

2. Примените правило ITO

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

В общем случае, если стохастический процесс X дан в терминах SDE, затем в правиле ITO говорится что преобразованный процесс G(t, X) удовлетворяет

dG=(μdGdX+σ22d2GdX2+dGdt)dt+dGdXσdB(t)

Мы принимаем, что логарифм цены дан одномерным аддитивным Броуновским движением, то есть, mu и sigma константы. Задайте функцию, которая применяет правило ITO, и используйте его, чтобы преобразовать аддитивное Броуновское движение в геометрическое броуновское движение.

ito = @(mu, sigma, G, t, X) ...
    deal( mu * diff(G, X) + sigma^2/2 * diff(G, X, X) + diff(G, t), ...
    diff(G, X) * sigma );

syms mu0 sigma0
[mu1, sigma1] = ito(mu0, sigma0, exp(X), t, X)
mu1 = 

eXσ022+μ0eX(exp (X) *sigma0^2)/2 + mu0*exp (X)

sigma1 = σ0eXsigma0*exp (X)

Замените exp(X) Y.

syms Y
mu1 = subs(mu1, X, log(Y));
sigma1 = subs(sigma1, X, log(Y));
f = f(t, log(Y));
s = s(t, log(Y));

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

r0 = 0;

c = subs(c, {mu, sigma}, {mu1, sigma1});
s = subs(s, {mu, sigma, r}, {mu1, sigma1, r0});
f = subs(f, {mu, sigma}, {mu1, sigma1});

3. Решите уравнение Фейнмана-Каца

Прежде чем можно будет преобразовать символьные выражения в указатели функции MATLAB, необходимо заменить вызовы функции, такие как diff(v(t, X), X) и v(t, X), с переменными. Можно использовать любые допустимые имена переменной MATLAB.

syms dvdx V;
dvX = diff(v, X);
c = subs(c, {dvX, v},  {dvdx, V});
f = subs(f, {dvX, v}, {dvdx, V});
s = subs(s, {dvX, v}, {dvdx, V});

Для плоской геометрии (симметрия перевода), устанавливает следующее значение:

m = 0;

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

muvalue = 0;
sigmavalue = 1;

c0 = subs(c, {mu0, sigma0}, {muvalue, sigmavalue});
f0 = subs(f, {mu0, sigma0}, {muvalue, sigmavalue});
s0 = subs(s, {mu0, sigma0}, {muvalue, sigmavalue});

Используйте matlabFunction создать указатель на функцию. Передайте коэффициенты c0, f0, и s0 в форме, требуемой pdepe, а именно, указатель на функцию с тремя выходными аргументами.

FeynmanKacPde = matlabFunction(c0, f0, s0, 'Vars', [t, Y, V, dvdx])
FeynmanKacPde = function_handle with value:
    @(t,Y,V,dvdx)deal(1.0,0.0,0.0)

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

FeynmanKacIC = @(x) x;

Числовое решение УЧП может только быть применено к конечной области. Поэтому необходимо задать граничное условие. Примите, что актив продается в тот момент времени, когда его повышения цен выше или падают ниже определенного предела, и таким образом решения v должен удовлетворить x - v = 0 за пределами точки x. Можно выбрать другое граничное условие, например, можно использовать v = 0 к опциям нокаута модели. Обнуление во втором и четвертом выходе указывает, что граничное условие не зависит от dvdx.

FeynmanKacBC = @(xl, vl, xr, vr, t) ...
    deal(xl - vl, 0, xr - vr, 0);

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

xgrid = linspace(0, 1, 101);

Выберите сетку времени. Из-за реверсирования времени, примененного в начале, это обозначает время, оставленное до итогового времени T.

tgrid = linspace(0, 1, 21);

Решите уравнение.

sol = pdepe(m,FeynmanKacPde,FeynmanKacIC,FeynmanKacBC,xgrid,tgrid);

Постройте решение. Ожидаемая продажная цена зависит почти линейно от цены во время t, и также слабо на t.

surf(xgrid, tgrid, sol)
title('Expected selling price of asset');
xlabel('price');
ylabel('time');
zlabel('expected final price');

Figure contains an axes. The axes with title Expected selling price of asset contains an object of type surface.

Состояние геометрического броуновского движения с дрейфом μ1 во время t логарифмически нормально распределенная случайная переменная с ожидаемым значением exp(μ1t) времена его начальное значение. Это описывает ожидаемую продажную цену актива, который никогда не продается из-за достижения предела.

Expected = transpose(exp(tgrid./2)) * xgrid;

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

sol2 = sol./Expected;
surf(xgrid, tgrid, sol2)
title('Ratio of expected final prices: with / without sales order at x=1')
xlabel('price');
ylabel('time');
zlabel('ratio of final prices');

Figure contains an axes. The axes with title Ratio of expected final prices: with / without sales order at x=1 contains an object of type surface.

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

plot(tgrid, sol2(:, xgrid == 1/2));
hold on
plot(tgrid, sol2(:, xgrid == 1/4));
legend('Limit at two times current price', 'Limit at four times current price');
xlabel('timespan the order is valid');
ylabel('ratio of final prices: with / without sales order');
hold off

Figure contains an axes. The axes contains 2 objects of type line. These objects represent Limit at two times current price, Limit at four times current price.

4. Вычислите ожидаемое время, чтобы продать актив

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

syms y(X)
exitTimeEquation(X) = subs(eq, {r, u}, {0, y(X)}) == -1
exitTimeEquation(X) = 

σ(t,X)22X2 y(X)2+μ(t,X)X y(X)=-1(сигма (t, X) ^2*diff (y (X), X, 2))/2 + mu (t, X) *diff (y (X), X) ==-1

Кроме того, y должен быть нуль на контурах. Для mu и sigma, вставьте фактический стохастический процесс на рассмотрении:

exitTimeGBM = subs(subs(exitTimeEquation, {mu, sigma}, {mu1, sigma1}), Y, X)
exitTimeGBM(X) = 

Xσ022+Xμ0X y(X)+X2σ022X2 y(X)2=-1((X*sigma0^2)/2 + X*mu0) *diff (y (X), X) + (X^2*sigma0^2*diff (y (X), X, 2))/2 ==-1

Решите эту задачу для произвольных границ интервала a и b.

syms a b
exitTime = dsolve(exitTimeGBM, y(a) == 0, y(b) == 0)
exitTime = 

2μ0σ5log(b)-2μ0σ4log(a)+aσ1σ02σ5σ4-bσ1σ02σ5σ4σ3-log(X)μ0+σ2σ7-σ6-aσ1σ02σ5+bσ1σ02σ4σ3+Xσ1σ02σ22μ02where  σ1=2μ0σ02  σ2=e-2μ0log(X)σ02  σ3=2μ02σ5-σ4  σ4=e-σ6σ02  σ5=e-σ7σ02  σ6=2μ0log(b)  σ7=2μ0log(a)(2*mu0*exp (-(2*mu0*log (a))/sigma0^2) *log (b) - 2*mu0*exp (-(2*mu0*log (b))/sigma0^2) *log (a) + a^ (2*mu0)/sigma0^2) *sigma0^2*exp (-(2*mu0*log (a))/sigma0^2) *exp (-(2*mu0*log (b))/sigma0^2) - b^((2*mu0)/sigma0^2)*sigma0^2*exp(-(2*mu0*log(a))/sigma0^2)*exp(-(2*mu0*log(b))/sigma0^2))/(2*mu0^2*(exp(-(2*mu0*log(a))/sigma0^2) - exp (-(2*mu0*log (b))/sigma0^2))) - журнал (X)/mu0 + (exp (-(2*mu0*log (X))/sigma0^2) * (2*mu0*log (a) - 2*mu0*log (b) - a^ (2*mu0)/sigma0^2) *sigma0^2*exp (-(2*mu0*log (a))/sigma0^2) + b^ (2*mu0)/sigma0^2) *sigma0^2*exp (-(2*mu0*log (b))/sigma0^2))) / (2*mu0^2* (exp (-(2*mu0*log (a))/sigma0^2) - exp (-(2*mu0*log (b))/sigma0^2))) + (X^ (2*mu0)/sigma0^2) *sigma0^2*exp (-(2*mu0*log (X))/sigma0^2)) / (2*mu0^2)

Поскольку вы не можете заменить mu0 = 0 непосредственно, вычислите предел в mu0 -> 0 вместо этого.

L = limit(subs(exitTime, sigma0, sigmavalue), mu0, muvalue)
L = -log(X)-log(a)log(X)-log(b)- (регистрируйте (X) - журнал (a)), * (журнал (X) - журнал (b))

L не определено в a = 0. Установите предположение что 0 < X < 1.

assume(0 < X & X < 1)

Используя значение b = 1 для правильной границы вычислите предел.

limit(subs(L, b, 1), a, 0, 'Right')
ans = sym(inf)

Ожидаемое выходное время бесконечно!

Для просмотра документации необходимо авторизоваться на сайте