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

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

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

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

Согласно теореме Фейнмана - Каца, u удовлетворяет дифференциальному уравнению в частных производных, $\frac{}{}\frac{{}^{\mathrm{}}}{\mathrm{}}\frac{{}^{\mathrm{}}}{{}^{\mathrm{}}}\frac{}{}\partial u\partial t+\sigma 22\partial 2u\partial x2+\mu \partial u\partial x-ur=0\mathrm{}$с конечным условием в момент времени 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 = 
$$\frac{{\sigma \left(T-t,X\right)}^2\hspace{0.17em}\frac{{\partial }^2}{\partial X^2}\mathrm{ }v\left(t,X\right)}{2}+\mu \left(T-t,X\right)\hspace{0.17em}\frac{\partial }{\partial X}\mathrm{ }v\left(t,X\right)-\frac{\partial }{\partial t}\mathrm{ }v\left(t,X\right)-v\left(t,X\right)\hspace{0.17em}r\left(T-t,X\right)$$

Решающее устройство pdepe требует дифференциального уравнения в частных производных в следующей форме. Здесь коэффициенты c, f, и s являются функциями x, t, vи $$\mathrm{\partial v/\partial x}$$.

Возможность решения уравнения 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 = 
$$\left(\begin{array}{ccc}-1& \frac{{\sigma \left(T-t,X\right)}^2}{2}& \mu \left(T-t,X\right)\hspace{0.17em}\frac{\partial }{\partial X}\mathrm{ }v\left(t,X\right)-v\left(t,X\right)\hspace{0.17em}r\left(T-t,X\right)\end{array}\right)$$

terms = $$\left(\begin{array}{ccc}\mathrm{dvt}& \mathrm{dvXX}& 1\end{array}\right)$$

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

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

Поэтому писать 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) удовлетворяет

Мы предполагаем, что логарифм цены задаётся одномерным аддитивным броуновским движением, то есть 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 = 
$$\frac{{\mathrm{e}}^X\hspace{0.17em}{{\sigma }_0}^2}{2}+{\mu }_0\hspace{0.17em}{\mathrm{e}}^X$$

sigma1 = $${\sigma }_0\hspace{0.17em}{\mathrm{e}}^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)

В качестве окончательного условия возьмите сопоставление идентификаторов. То есть, погашение в момент $$\mathrm{времени\; T}$$ определяется самой ценой актива. Эту строку можно изменить для изучения производных инструментов.

FeynmanKacIC = @(x) x;

Численное решение PDE может применяться только к конечной области. Поэтому необходимо указать граничное условие. Предположим, что актив продается в момент, когда его цена поднимается выше или ниже определенного предела, и, таким образом, решение v должен удовлетворять x - v = 0 в граничных точках x. Можно выбрать другое граничное условие, например, использовать v = 0 для моделирования вариантов выбивания. Нули на втором и четвертом выходных данных указывают, что граничное условие не зависит от $\frac{}{\mathrm{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');

Состояние геометрического броуновского движения с дрейфом $$\mathrm{мк1}$$ в момент времени t - логнормально распределенная случайная величина с ожидаемым значением $$\exp {(}_{\mathrm{}}мк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');

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

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

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

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

exitTimeEquation(X) = 
$$\frac{{\sigma \left(t,X\right)}^2\hspace{0.17em}\frac{{\partial }^2}{\partial X^2}\mathrm{ }y\left(X\right)}{2}+\mu \left(t,X\right)\hspace{0.17em}\frac{\partial }{\partial X}\mathrm{ }y\left(X\right)=-1$$

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

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

exitTimeGBM(X) = 
$$\left(\frac{X\hspace{0.17em}{{\sigma }_0}^2}{2}+X\hspace{0.17em}{\mu }_0\right)\hspace{0.17em}\frac{\partial }{\partial X}\mathrm{ }y\left(X\right)+\frac{X^2\hspace{0.17em}{{\sigma }_0}^2\hspace{0.17em}\frac{{\partial }^2}{\partial X^2}\mathrm{ }y\left(X\right)}{2}=-1$$

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

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

exitTime = 
$$\begin{array}{l}\frac{2\hspace{0.17em}{\mu }_0\hspace{0.17em}{\sigma }_5\hspace{0.17em}\log \left(b\right)-2\hspace{0.17em}{\mu }_0\hspace{0.17em}{\sigma }_4\hspace{0.17em}\log \left(a\right)+a^{{\sigma }_1}\hspace{0.17em}{{\sigma }_0}^2\hspace{0.17em}{\sigma }_5\hspace{0.17em}{\sigma }_4-b^{{\sigma }_1}\hspace{0.17em}{{\sigma }_0}^2\hspace{0.17em}{\sigma }_5\hspace{0.17em}{\sigma }_4}{{\sigma }_3}-\frac{\log \left(X\right)}{{\mu }_0}+\frac{{\sigma }_2\hspace{0.17em}\left({\sigma }_7-{\sigma }_6-a^{{\sigma }_1}\hspace{0.17em}{{\sigma }_0}^2\hspace{0.17em}{\sigma }_5+b^{{\sigma }_1}\hspace{0.17em}{{\sigma }_0}^2\hspace{0.17em}{\sigma }_4\right)}{{\sigma }_3}+\frac{X^{{\sigma }_1}\hspace{0.17em}{{\sigma }_0}^2\hspace{0.17em}{\sigma }_2}{2\hspace{0.17em}{{\mu }_0}^2}\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\frac{2\hspace{0.17em}{\mu }_0}{{{\sigma }_0}^2}\\ \\ \mathrm{ }{\sigma }_2={\mathrm{e}}^{-\frac{2\hspace{0.17em}{\mu }_0\hspace{0.17em}\log \left(X\right)}{{{\sigma }_0}^2}}\\ \\ \mathrm{ }{\sigma }_3=2\hspace{0.17em}{{\mu }_0}^2\hspace{0.17em}\left({\sigma }_5-{\sigma }_4\right)\\ \\ \mathrm{ }{\sigma }_4={\mathrm{e}}^{-\frac{{\sigma }_6}{{{\sigma }_0}^2}}\\ \\ \mathrm{ }{\sigma }_5={\mathrm{e}}^{-\frac{{\sigma }_7}{{{\sigma }_0}^2}}\\ \\ \mathrm{ }{\sigma }_6=2\hspace{0.17em}{\mu }_0\hspace{0.17em}\log \left(b\right)\\ \\ \mathrm{ }{\sigma }_7=2\hspace{0.17em}{\mu }_0\hspace{0.17em}\log \left(a\right)\end{array}$$

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

L = limit(subs(exitTime, sigma0, sigmavalue), mu0, muvalue)

L = $$-\left(\log \left(X\right)-\log \left(a\right)\right)\hspace{0.17em}\left(\log \left(X\right)-\log \left(b\right)\right)$$

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

assume(0 < X & X < 1)

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

limit(subs(L, b, 1), a, 0, 'Right')

ans = $$\infty$$

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