Solving ODE's

Question

Riaz Patel 2018-11-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/431965-solving-ode-s

评论： Riaz Patel 2018-11-26

在 MATLAB Online 中打开

Hi all,

I am trying to solve the following ODE's for a given value of u:

where:

and t = T-tau.

I solved the ODE for E using ODE45. The solution returns Nan for the last few values of tau. I also dont know how to take this solution for E(u,tau) and use this to find A(u,tau).

My code is below:

clear;
clc;
%Parameters
%Heston Parameters
S0 = 100;
T = 1;
k = 1.5768;
sigma = 0.0571;
v0 = 0.0175;
vb = 0.0398;
%Hull-White parameters
lambda = 0.05;
r0 = 0.07;
theta = 0.07;
eta = 0.005;
%correlations
pxv =  - 0.5711;
pxr =  0.2;
pvr = 0;
tau = T; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cf = @(t) (1/(4*k))*sigma^2*(1-exp(-k*t));
d = (4*k*vb)/(sigma^2);
lambdaf = @(t) (4*k*v0*exp(-k*t))./(sigma^2*(1-exp(-k*t)));
lambdaC = @(t) sqrt(cf(t).*(lambdaf(t)-1) + cf(t)*d + (cf(t)*d)./(2*(d+lambdaf(t))));
D1 = @(u) sqrt((sigma*pxv*1i*u-k).^2 - sigma^2*1i*u.*(1i*u-1));
g = @(u) (k-sigma*pxv*1i*u - D1(u))./(k-sigma*pxv*1i*u + D1(u));
B = @(u,tau) 1i*u;
C = @(u,tau) (1i*u-1)*(1/lambda)*(1-exp(-lambda*tau));
D = @(u,tau) ((1 -exp(-D1(u)*tau))./(sigma^2*(1-g(u).*exp(-D1(u)*tau)))).*(k-sigma*pxv*1i*u-D1(u));
%ODE's that are solved numerically
muxi = @(t) (1/(2*sqrt(2)))*(gamma(0.5*(1+d))/sqrt(cf(t)))*(hypergeom(-0.5,0.5*d,-0.5*lambdaf(t))*(1/gamma(0.5*d))*sigma^2*exp(-k*t)*0.5 + hypergeom(0.5,1+0.5*d,-0.5*lambdaf(t))*(1/gamma(1+0.5*d))*((v0*k)/(1-exp(k*t))));
phixi = @(t) sqrt(k*(vb-v0)*exp(-k*t) - 2*lambdaC(t)*muxi(t));
u = 10;
EODE = @(tau,y) pxr*eta*B(u,tau)*C(u,tau) + phixi(T-tau)*pxv*B(u,tau)*y + sigma*phixi(T-tau)*D(u,tau)*y;
AODE = @(tau,y) k*vb*D(u,tau) + lambda*theta*C(u,tau) + muxi(T-tau)*E() +eta^2*0.5*C(u,tau)^2 + (phixi(T-tau))^2*0.5*E()^2;
%what do i put in for E() in the line above?
[tau, E] = ode45(EODE,[0 T],0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

madhan ravi 2018-11-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/431965-solving-ode-s#answer_348914

编辑：madhan ravi 2018-11-26

在 MATLAB Online 中打开

option 1;

tau(~isnan(E)) 
E(~isnan(E)) %removes nan values

option 2;

Let even function play the role it stops the evaluation when solution become a NaN

clear;
clc;
%Parameters
%Heston Parameters
S0 = 100;
T = 1;
k = 1.5768;
sigma = 0.0571;
v0 = 0.0175;
vb = 0.0398;
%Hull-White parameters
lambda = 0.05;
r0 = 0.07;
theta = 0.07;
eta = 0.005;
%correlations
pxv =  - 0.5711;
pxr =  0.2;
pvr = 0;
tau = T; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cf = @(t) (1/(4*k))*sigma^2*(1-exp(-k*t));
d = (4*k*vb)/(sigma^2);
lambdaf = @(t) (4*k*v0*exp(-k*t))./(sigma^2*(1-exp(-k*t)));
lambdaC = @(t) sqrt(cf(t).*(lambdaf(t)-1) + cf(t)*d + (cf(t)*d)./(2*(d+lambdaf(t))));
D1 = @(u) sqrt((sigma*pxv*1i*u-k).^2 - sigma^2*1i*u.*(1i*u-1));
g = @(u) (k-sigma*pxv*1i*u - D1(u))./(k-sigma*pxv*1i*u + D1(u));
B = @(u,tau) 1i*u;
C = @(u,tau) (1i*u-1)*(1/lambda)*(1-exp(-lambda*tau));
D = @(u,tau) ((1 -exp(-D1(u)*tau))./(sigma^2*(1-g(u).*exp(-D1(u)*tau)))).*(k-sigma*pxv*1i*u-D1(u));
%ODE's that are solved numerically
muxi = @(t) (1/(2*sqrt(2)))*(gamma(0.5*(1+d))/sqrt(cf(t)))*(hypergeom(-0.5,0.5*d,-0.5*lambdaf(t))*(1/gamma(0.5*d))*sigma^2*exp(-k*t)*0.5 + hypergeom(0.5,1+0.5*d,-0.5*lambdaf(t))*(1/gamma(1+0.5*d))*((v0*k)/(1-exp(k*t))));
phixi = @(t) sqrt(k*(vb-v0)*exp(-k*t) - 2*lambdaC(t)*muxi(t));
u = 10;
EODE = @(tau,y) pxr*eta*B(u,tau)*C(u,tau) + phixi(T-tau)*pxv*B(u,tau)*y + sigma*phixi(T-tau)*D(u,tau)*y;
AODE = @(tau,y) k*vb*D(u,tau) + lambda*theta*C(u,tau) + muxi(T-tau)*E() +eta^2*0.5*C(u,tau)^2 + (phixi(T-tau))^2*0.5*E()^2;
%what do i put in for E() in the line above?
opts = odeset('Events',@stopfunc) %it will stop integration when distance becomes zero
[tau, E] = ode45(EODE,[0 T],0,opts); %function call
function [position,isterminal,direction] = stopfunc(t,x) %function definition
position = x(1); % The value that we want to be zero
isterminal = 1;  % Halt integration 
direction = ~isnumeric(x(1));   % The zero can be approached from either direction
end

3 个评论
显示 1更早的评论隐藏 1更早的评论

madhan ravi 2018-11-26

Perfect Torsten ! , you are the man!

Riaz Patel 2018-11-26

在 MATLAB Online 中打开

Thanks guys!

So if im understanding correctly, you are writing this as a system of ODE's? what does the y(1) mean?

Secondly, why does this piece of code below not work? And how do you save the output? Do you get an output for E,A and tau? Sorry, I am new to solving differential equations numerically.

ode45(EAODE,[0 T],[0 0]);

请先登录，再进行评论。