Solving ODE's
3 次查看(过去 30 天)
显示 更早的评论
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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0 个评论
回答(1 个)
madhan ravi
2018-11-26
编辑:madhan ravi
2018-11-26
option 1;
tau(~isnan(E))
E(~isnan(E)) %removes nan values
option 2;
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
3 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!