MATLAB: Speeding up solution to system of ODES

odetrapz

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!

Best Answer

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
Related Question