Speeding up solution to system of ODES
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Riaz Patel
am 17 Dez. 2018
Kommentiert: Torsten
am 19 Dez. 2018
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!
4 Kommentare
Akzeptierte Antwort
Torsten
am 18 Dez. 2018
Bearbeitet: Torsten
am 18 Dez. 2018
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
6 Kommentare
Torsten
am 19 Dez. 2018
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.
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differentiation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!