Этот пример получает дифференциальное уравнение с частными производными, которое описывает ожидаемую окончательную цену актива, цена которого является стохастическим процессом, данным стохастическим дифференциальным уравнением.
Включенные шаги:
Задайте параметры модели Используя стохастические дифференциальные уравнения
Примените правило ITO
Решите уравнение Фейнмана-Каца
Вычислите ожидаемое время, чтобы продать актив
Модель за цену актива X(t)
, заданный во временном интервале [0,T]
, является стохастическим процессом, заданным стохастическим дифференциальным уравнением формы , где B(t)
является Винеровским процессом с модульным параметром отклонения.
Создайте символьную переменную T
и следующие символьные функции:
r(t)
является непрерывной функцией, представляющей точечную процентную ставку. Этот уровень определяет коэффициент дисконтирования для итоговой выплаты в то время T
.
u(t,x)
является ожидаемым значением обесцененной будущей цены, вычисленной как при условии X(t) = x
.
и дрейф и диффузия стохастического процесса X(t)
.
syms mu(t, X) sigma(t, X) r(t, X) u(t, X) T
Согласно теореме Фейнмана-Каца, u
удовлетворяет дифференциальное уравнение с частными производными с итоговым условием во время 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 =
Решатель pdepe
требует дифференциального уравнения с частными производными в следующей форме. Здесь коэффициенты c
, f
и s
являются функциями x
, t
, v
, и .
Чтобы смочь решить уравнение 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 =
terms =
Теперь, можно переписать 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);
Цены активов следуют за мультипликативным процессом. Таким образом, логарифм цены может быть описан с точки зрения SDE, но ожидаемое значение самой цены представляет интерес, потому что это описывает прибыль, и таким образом нам нужен SDE для последнего.
В целом, если стохастический процесс, который X
дан с точки зрения SDE, то в правиле ITO говорится, что преобразованный процесс G(t, X)
удовлетворяет
Мы принимаем, что логарифм цены дан одномерным аддитивным Броуновским движением, то есть, 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 =
sigma1 =
Замените 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});
Прежде чем можно будет преобразовать символьные выражения в указатели функции 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,(Y.^2.*dvdx)./2.0,(Y.*dvdx)./2.0)
Как итоговое условие, возьмите единичное отображение. Таким образом, выплата во время дан самой ценой активов. Можно изменить эту строку в порядке исследовать производные инструменты.
FeynmanKacIC = @(x) x;
Числовое решение УЧП может только быть применено к конечной области. Поэтому необходимо задать граничное условие. Примите, что актив продается в тот момент времени, когда его повышения цен выше или падают ниже определенного предела, и таким образом решения, v
должен удовлетворить x - v = 0
в граничных точках x
. Можно выбрать другое граничное условие, например, можно использовать v = 0
для образцовых опций нокаута. Обнуление во втором и четвертом выводе указывает, что граничное условие не зависит от .
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');
Состояние геометрического броуновского движения с дрейфом во время t
является логарифмически нормально распределенной случайной переменной с ожидаемым значением времена его начальное значение. Это описывает ожидаемую продажную цену актива, который никогда не продается из-за достижения предела.
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
Это - результат учебника, что ожидаемое выходное время, когда предел достигнут и актив, продается, дан следующим уравнением:
syms y(X)
exitTimeEquation(X) = subs(eq, {r, u}, {0, y(X)}) == -1
exitTimeEquation(X) =
Кроме того, y
должен быть нулем на контурах. Для mu
и sigma
, вставьте фактический стохастический процесс на рассмотрении:
exitTimeGBM = subs(subs(exitTimeEquation, {mu, sigma}, {mu1, sigma1}), Y, X)
exitTimeGBM(X) =
Решите эту проблему для произвольных границ интервала a
и b
.
syms a b exitTime = dsolve(exitTimeGBM, y(a) == 0, y(b) == 0)
exitTime =
Поскольку вы не можете заменить mu0 = 0
непосредственно, вычислите предел в mu0 -> 0
вместо этого.
L = limit(subs(exitTime, sigma0, sigmavalue), mu0, muvalue)
L =
L
не определен в a = 0
. Установите предположение что 0 < X < 1
.
assume(0 < X & X < 1)
Используя значение b = 1
для правильной границы, вычислите предел.
limit(subs(L, b, 1), a, 0, 'Right')
ans =