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

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

Задействованы следующие шаги:

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

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

  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)(sigma (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 требует дифференциального уравнения с частными производными в следующей форме. Здесь коэффициенты c, f, и s являются функциями x, t, v, и v/x.

cvt=fx+s

Чтобы иметь возможность решить уравнение eq2 с pdepe решатель, map 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), sigma (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. Применить правило Ито

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

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

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

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

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;

Численное решение PDE может применяться только к конечной области. Поэтому необходимо задать граничное условие. Предположим, что актив продается в момент, когда его цена поднимается выше или падает ниже определенного предела, и, следовательно, решение 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 является lognormally распределенной случайной переменной с ожидаемым значением 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(sigma (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) *журнал (b) - 2*mu0*exp (-(2*mu0*log (b))/sigma0^2) *журнал (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)- (log (X) - log (a)) * (log (X) - log (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)

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

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