Speeding up solution to system of ODES

Hi all,

function [y] = H2HWCF(u,S0,T,tau,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr)

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

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

[tausol, ysol] = ode23(EAODE,[0 T],[0 0]);              

E = ysol(end,1);
A = ysol(end,2);

y = exp(A+B(u,tau)*log(S0)+C(u,tau)*r0+D(u,tau)*v0 + E*sqrt(v0));

end

In this function, i have closed form solutions for B,C,D. this is not possible for E and A. So these are solved as a system of ODEs. I am solving the ODES from [0,T] with initial conditions at t = 0. I only need the values for E and A and time T though.

The problem im having is that i have to take this function and integrate it. Im using trapz to do my integration. I have to obviously evulate this function many times in order to get an accurate approximation for the integral. This is taking too long because of the ODE's that need to be solved. Is there anyway for me to speed up the run time?

Thanks!

[madhan ravi]

please provide all the datas(upload equation in latex form) and how you call the function with inputs

[Riaz Patel]

%Parameters
%Heston Parameters
S0 = 100;
K = 100;
T = 1;
k = 2.5;
sigma = 0.5;
v0 = 0.06;
vb = 0.06;

%Hull-White parameters
lambda = 0.05;
r0 = 0.07;
theta = 0.07;
eta = 0.01;

%correlations
pxv =  - 0.3;
pxr =  0.2;
pvr = 0;

%MC parameters 
N = 200;                  
dt = T/N;
t = (0:dt:T);
n = 100000;

alpha = 0.75;
vmax = 250;
H2HWCFtemp = @(u) H2HWCF(u,S0,T,T,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr);
psi = @(v) (H2HWCFtemp(v-(alpha+1)*1i))./(alpha^2+alpha-v.^2+1i*(2*alpha+1)*v);

vd =  0:0.1:vmax;
y = zeros(size(vd));

for i = 1:length(vd)
    y(i) = psi(vd(i));
end

This is my function call. So I am evaluating psi(v) at each v in vd. The psi function calls the original function which solves the ODEs.

The original equation is:

where

[madhan ravi]

B(u,0)=iu initial condition?

[Riaz Patel]

Yes. why?

Answer 1

[Torsten]

function main
%Parameters

%Heston Parameters
S0 = 100;
K = 100;
T = 1;
k = 2.5;
sigma = 0.5;
v0 = 0.06;
vb = 0.06;

%Hull-White parameters
lambda = 0.05;
r0 = 0.07;
theta = 0.07;
eta = 0.01;

%correlations
pxv =  - 0.3;
pxr =  0.2;
pvr = 0;

%MC parameters 
N = 200;                  
dt = T/N;
t = (0:dt:T);
n = 100000;

alpha = 0.75;
vmax = 250;
H2HWCFtemp = @(u) H2HWCF(u,S0,T,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr);
psi = @(v,y) (H2HWCFtemp(v-(alpha+1)*1i))./(alpha^2+alpha-v.^2+1i*(2*alpha+1)*v);
% Integrate psi from v=0 to v=vmax
[V,PSI]=ode15s(psi,[0 vmax],0);
% Display integral_{v=0}^{v=vmax} psi(v) dv
V(end,1)
end

function y = H2HWCF(u,S0,T,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr)

[tausol, ysol] = ode15s(@(tau,y)EAODE(tau,y,u,S0,T,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr),[0 T],[0 0]);  
E = ysol(end,1);
A = ysol(end,2);

D1 = sqrt((sigma*pxv*1i*u-k).^2 - sigma^2*1i*u.*(1i*u-1));
g = (k-sigma*pxv*1i*u - D1)./(k-sigma*pxv*1i*u + D1);

B = 1i:u;
C = (1i*u-1)*(1/lambda)*(1-exp(-lambda*T));
D = ((1 -exp(-D1*T))./(sigma^2*(1-g.*exp(-D1*T)))).*(k-sigma*pxv*1i*u-D1);

y = exp(A+B*log(S0)+C*r0+D*v0 + E*sqrt(v0));  

end

function dy = EAODE(tau,y,u,S0,T,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr)

cf = (1/(4*k))*sigma^2*(1-exp(-k*(T-tau)));
d = (4*k*vb)/(sigma^2);
lambdaf = (4*k*v0*exp(-k*(T-tau)))./(sigma^2*(1-exp(-k*(T-tau))));
lambdaC = sqrt(cf.*(lambdaf-1) + cf*d + (cf*d)./(2*(d+lambdaf)));

D1 = sqrt((sigma*pxv*1i*u-k).^2 - sigma^2*1i*u.*(1i*u-1));
g = (k-sigma*pxv*1i*u - D1)./(k-sigma*pxv*1i*u + D1);

B = 1i*u;
C = (1i*u-1)*(1/lambda)*(1-exp(-lambda*tau));
D = ((1 -exp(-D1*tau))./(sigma^2*(1-g.*exp(-D1*tau)))).*(k-sigma*pxv*1i*u-D1);

muxi = (1/(2*sqrt(2)))*(gamma(0.5*(1+d))/sqrt(cf))*...
    (hypergeom(-0.5,0.5*d,-0.5*lambdaf)*(1/gamma(0.5*d))*sigma^2*exp(-k*(T-tau))*0.5 + ...
    hypergeom(0.5,1+0.5*d,-0.5*lambdaf)*(1/gamma(1+0.5*d))*((v0*k)/(1-exp(k*(T-tau)))));
phixi = sqrt(k*(vb-v0)*exp(-k*(T-tau)) - 2*lambdaC*muxi);

dy = zeros(2,1);
dy(1) = pxr*eta*B*C + phixi*pxv*B*y(1) + sigma*phixi*D*y(1);
dy(2) = k*vb*D + lambda*theta*C + muxi*y(1)+eta^2*0.5*C^2 + phixi^2*0.5*y(1)^2;

end

[Riaz Patel]

Hi Torsten.

Thanks for the reply!

i ran the revised code for H2HWCF just to check that i was getting the same answers as before, i got an error though.

function y = H2HWCFtest(u,S0,T,tau,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr)

[tausol, ysol] = ode15s(@(tau,y)EAODE(tau,y,u,S0,T,tau,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr),[0 T],[0 0]);  
E = ysol(end,1);
A = ysol(end,2);

D1 = sqrt((sigma*pxv*1i*u-k).^2 - sigma^2*1i*u.*(1i*u-1));
g = (k-sigma*pxv*1i*u - D1)./(k-sigma*pxv*1i*u + D1);

B = 1i:u;
C = (1i*u-1)*(1/lambda)*(1-exp(-lambda*tau));
D = ((1 -exp(-D1*tau))./(sigma^2*(1-g.*exp(-D1*tau)))).*(k-sigma*pxv*1i*u-D1);

y = exp(A+B*log(S0)+C*r0+D*v0 + E*sqrt(v0));  
end

function dy = EAODE(tau,y,u,S0,T,tau,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr)

cf = (1/(4*k))*sigma^2*(1-exp(-k*(T-tau)));
d = (4*k*vb)/(sigma^2);
lambdaf = (4*k*v0*exp(-k*(T-tau)))./(sigma^2*(1-exp(-k*(T-tau))));
lambdaC= sqrt(cf.*(lambdaf-1) + cf*d + (cf*d)./(2*(d+lambdaf)));

D1 = sqrt((sigma*pxv*1i*u-k).^2 - sigma^2*1i*u.*(1i*u-1));
g = (k-sigma*pxv*1i*u - D1)./(k-sigma*pxv*1i*u + D1);

B = 1i*u;
C = (1i*u-1)*(1/lambda)*(1-exp(-lambda*tau));
D = ((1 -exp(-D1*tau))./(sigma^2*(1-g.*exp(-D1*tau)))).*(k-sigma*pxv*1i*u-D1);

muxi = (1/(2*sqrt(2)))*(gamma(0.5*(1+d))/sqrt(cf))*...
    (hypergeom(-0.5,0.5*d,-0.5*lambdaf)*(1/gamma(0.5*d))*sigma^2*exp(-k*(T-tau))*0.5 + ...
    hypergeom(0.5,1+0.5*d,-0.5*lambdaf)*(1/gamma(1+0.5*d))*((v0*k)/(1-exp(k*(T-tau)))));
phixi = sqrt(k*(vb-v0)*exp(-k*(T-tau)) - 2*lambdaC*muxi);

dy = zeros(2,1);
dy(1) = pxr*eta*B*C + phixi*pxv*B*y(1) + sigma*phixi*D*y(1);
dy(2) = k*vb*D + lambda*theta*C + muxi*y(1)+eta^2*0.5*C^2 + phixi^2*0.5*y(1)^2;
end

i ran it:

u = 10;
tic
H2HWCFtest(u,S0,T,tau,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr)
toc

I got the following error

Error: File: H2HWCFtest.m

Line: 17 Column: 1

The variable tau was

mentioned more than once as

an input.

Also, I see that one of the changes you made was to ode15s, i changed that in my original formulation of the function H2HWCF and the run time increased dramatically. From around 15 seconds with ode23 to about 61 seconds.

[Torsten]

I think you didn't copy the code correctly. "tau" is removed from the list of input parameters to H2HWCF.

[Riaz Patel]

OH MY WORD!!!! run time is around 4 seconds per an evaluation of H2HWCF. Thank you so much!

Why did you choose ode15s instead of ode23. I see that ode23 is much faster than ode15s. Is ode15s more accurate?

Lastly, in my original question i did not fully show what i was doing. I actually have to integrate the following function:

fintegrand = @(v) exp(-1i*v*log(K)).*psi(v);

Now by using a for loop to evalation psi(v) and then saving the resulting vector. I can then evaluate exp(-1i*v*log(K)).*psi(v) for each value of K that i am interested in and use trapz to integrate the resulting vector. This means that for each K, i only need to evalute the psi(v) once.

My code looks something like this:

alpha = 0.75;
H2HWCFtemp = @(u) H2HWCF(u,S0,T,T,k,sigma,v0,vb,lambda,r0,theta,eta,pxv,pxr);
%this will now be your formulation of the function H2HWCF 
% which is much faster than my own 
psi = @(v) (H2HWCFtemp(v-(alpha+1)*1i))./(alpha^2+alpha-v.^2+1i*(2*alpha+1)*v);
 
vmax = 50;
vd = 0:0.1:vmax;
y = zeros(size(vd));
 
tic
 parfor i = 1:length(vd)
     y(i) = psi(vd(i)); %saved this vector as 'P1_6_0_50_01.mat'
 end
 toc
 
 K = [40 80 100 120 160];
 Res = zeros(length(K),2);
 
 for i = 1:length(K)
     load('P1_6_0_50_01')
     y = y.*exp(-1i*vd*log(K(i)));
     Res(i,1) = K(i);
     Res(i,2) = (1/pi)*exp(-alpha*log(K(i)))*real(trapz(vd,y));
 end

do you think this is a better approach than the first part of the function you gave me which uses the ode solver to integrate psi(v)?

[madhan ravi]

https://in.mathworks.com/help/matlab/math/choose-an-ode-solver.html

[Riaz Patel]

Just an update, i get these sorts of errors when i use ode15s.

Warning: Failure at t=1.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t.

> In ode15s (line 668)

I also get slightly different answers with the two solvers. how do i determine which one is more accurate for this specific system of odes?

[Torsten]

You will have to play with the solvers to see which one is the most appropriate for your problem. Usually, ODE15S is most accurate, but slower for non-stiff problems.

I suggested to circumvent application of "trapz" by using ODE15S also to integrate psi (or your "fintegrand" from above). This may save computation time because the ODE solver chooses its step in "v" automatically according to the tolerance you have chosen. At the moment, you "blindly" use 2500 function evaluation for psi. I leave it to you which method you prefer to implement.

Best wishes

Torsten.