Solving ODE's

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);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Answer 1

[madhan ravi]

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

[Torsten]

Use

EAODE = @(tau,y) [pxr*eta*B(u,tau)*C(u,tau) + phixi(T-tau)*pxv*B(u,tau)*y(1) + sigma*phixi(T-tau)*D(u,tau)*y(1);k*vb*D(u,tau) + lambda*theta*C(u,tau) + muxi(T-tau)*y(1)+eta^2*0.5*C(u,tau)^2 + (phixi(T-tau))^2*0.5*y(1)^2];

to solve simultaneously for E and A.

[madhan ravi]

Perfect Torsten ! , you are the man!

[Riaz Patel]

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]);