I'm attempting to implement aforementioned SVCJ model by Duffie et al in MATLAB. so far without success. It's supposed to price vanilla (european) calls . parameters provided, the expected price is: ~6.0411. The function throws the error:
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.1e+03. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Using quadv instead of integral didn't yield the expected result, but instead just returned NaNs.
function price = SVCJDuffie2000Test(S0, K, V0, r, mus, muv, lambda, kappav, ...
thetav, sigmav, sigmas, rho, rhoj, q, t, T)
% this function should in theory calculate the analytical solution for
% the SVCJ model by Duffie et al
%
%
%
%
% T=0.25
% K = 100;
% S0 = 85;
% kappav=4;
% lambda=4;
% rhoj=-0.5;
% rho=-0.5;
% thetav=0.04;
% r=0.05;
% sigmas=0.06;
% mu = 0.2380
% muv=0.02;
% muJ = -0.04;
% sigmav=0.1;
% V0 = 0.25;
%
%
% mus = log((1+mu)*(1-rhoj*muv))-0.5*(sigmas^2);
% SVCJDuffie2000Test( S0, K, V0, r, mus, muv, lambda, kappav, thetav, sigmav, sigmas, rho, rhoj, 0,0,T )
% S0, V0, r, mu, muv, lambda, kappav,
% thetav, sigmav, sigmas, rho, rhoj,
y = log(S0);
X0 = y;
% nu, initial vola, nubar long run mean vola
nu = V0;
c = K;
nubar = thetav;
sigmav = sigmav;
sigmacy = sigmas;
mucy = mus;
mucv = muv;
rhoj = rhoj;
rhobar = rho;
kappav = kappav;
zetabar = q;
lambdac = lambda;
% specific to SVCJ model
lambdabar = lambda;
r = r;
% not needed for SVJJ, ie SVCJ
sigmay = 0;
lambday = 0;
lambdav = 0;
muv = 0;
muy = 0;
mubar = thetafunc(1,0) - 1
function retval = alphabar(tau,u)
thetainter = lambdabar^-1*(lambdac*fc(u,tau));
retval = alpha0(tau, u) - lambdabar*tau*(1+mubar*u)+lambdabar*thetainter;
end
function retval = betabar(tau,u)
a = u*(1-u);
b = sigmav*rhobar*u-kappav;
gamma = sqrt(b^2+a*sigmav^2);
retval = a*(1-exp(-gamma*tau)) / (2*gamma-(gamma+b)*(1-exp(-gamma*tau)));
end
function retval = alpha0(tau, u)
a = u*(1-u);
b = sigmav*rhobar*u-kappav;
gamma = sqrt(b^2+a*sigmav^2);
retval = -r*tau+(r-zetabar)*u*tau-kappav*nubar*...
(((gamma+b)/sigmav^2)*tau + ...
(2/sigmav^2)*log(1-((gamma+b)/(2*gamma))*(1-exp(-gamma*tau))));
end
function retval = fy(u,tau)
retval = tau*exp(muy*u+1/2*sigmay^2*u^2);
end
function retval = fv(u,tau)
a = u*(1-u);
b = sigmav*rhobar*u-kappav;
gamma = sqrt(b^2+a*sigmav^2);
c_ = 1 - rhoj*mucv*u;
retval = (gamma-b)/((gamma-b)+muv*a)*tau-(2*muv*a)/((gamma*c_)^2-(b*c_-muv*a)^2)...
*log(1-((gamma+b)*c_-muv*a)/(2*gamma)*(1-exp(-gamma*tau)));
end
function retval = fc(u,tau)
a = u*(1-u);
b = sigmav*rhobar*u-kappav;
c_ = 1 - rhoj*mucv*u;
gamma = sqrt(b^2+a*sigmav^2);
d_ = (gamma-b)/((gamma-b)*c_+mucv*a)*tau-(2*mucv*a)/((gamma*c_)^2-(b*c_-mucv*a)^2)...
*log(1-((gamma+b)*c_-mucv*a)/(2*gamma*c_)*(1-exp(-gamma*tau)));
retval = exp(mucy*u+sigmacy^2*u^2/2)*d_;
end
function retval = thetafunc(c1, c2)
retval = lambdabar^-1*(lambdac*thetacfunc(c1,c2));
end
function retval = thetayfunc(c)
retval = exp(muy*c+1/2*sigmay^2*c^2);
end
function retval = thetavfunc(c)
retval = 1/(1-muv*c);
end
function retval = thetacfunc(c1,c2)
retval = exp(mucy*c1+1/2*sigmacy^2*c1^2) / ...
(1-mucv*c2-rhoj*mucv*c1);
end
function retval = psi(y_, u, t, T)
retval = exp(alphabar(T-t,u)+u*y_+ betabar(T-t,u)*nu);
end
% extended transform p.13
function retval = psichi(nu_, u, y, t, T)
retval = exp(alphabar(T-t,u)+u*y+ betabar(T-t,u)*nu_);
end
function retval = Gab(a,b, y, X0_, T)
integratefunc = @(nu_) imag(psi(complex(a,nu_*b),X0_, 0, T) * exp(complex(0,-nu_*y))) / nu_;
retval = psi(a,X0_,0,T)/2 - 1/pi*integral(integratefunc, 0, Inf, 'ArrayValued',true);
%retval = 1;
end
% depends on payoff function, see p. 6, 18
%aT = -r*T;
%bT = 0;
%d = 1;
% could also be:
aT = 0;
bT = 1;
d = -r*T;
GbTplusdminusd = Gab(bT+d,-d,-log(c),X0,T)
GbTminusd =Gab(bT,-d,-log(c),X0,T)
price = exp(aT)*GbTplusdminusd-c*exp(aT)*GbTminusd;
end
Functions in the paper:
thetainer: p. 23
alpha0: p.22
alphabar: p.22
betabar: p.22
fy,fv,fc: p. 23
thetayfunc, thetacfunc, thetavfunc: p. 22
psi: p.21
Gab: p. 13
not a function, but depending on the derivative type:
GbTplusdminusd,GbTminusd, price: p. 18
lambdabar = lambda, according to: p. 24