Speeding up solution to system of ODES
조회 수: 1 (최근 30일)
이전 댓글 표시
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
채택된 답변
Torsten
2018년 12월 18일
편집: Torsten
2018년 12월 18일
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
Torsten
2018년 12월 19일
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.
추가 답변 (0개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!