SVCJ (SVJJ) Duffie et. al Model implementation in Matlab

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.

    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

-
Whats the value of lambda for the example? – Theja Jul 17 '14 at 11:16
oh, forgot that one: lambda=4; – jcfrei Jul 17 '14 at 11:31
I've looked at the code and the paper and it's hard for me to map the code to the formulas. Also, it's not clear to me how you derived the constants. Can you clarify the question? – Bob Jansen Jul 19 '14 at 16:05
the constants (lambda, thetav, etc.) are taken from another paper and the price mentioned should be correct (theses constants have to be calibrated usually). regarding mapping the formulas, I've added the page nr. for each formula. – jcfrei Jul 21 '14 at 10:10
Thanks, unfortunately I won't be looking at it soon but this makes it more clear. – Bob Jansen Jul 21 '14 at 10:48

1. About the integration problem: Your integrand is highly oscillatory, and the adaptive quadrature of Matlab doesn't handle such integrands very well. In general, I would recommend Mathematica when Matlab's standard procedures don't perform well. In this case, a Levin-type method would perform much better. The reason that quadv produces NaN values is because of the singularity at 0. You would have to start the integration at a small, positive value.
2. In any case, there are a few typos in your code. I have attempted to correct the code. The integration works well now, but I still do not get the value that you were expecting (I have used t=0, q=0, mu=0 for the unspecified parameters; mu is needed in the mus formula). I have added XXX as comments in the lines that I have changed. Also, I made all operations vector valued to speed up the integration.

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;
% 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; V0]; % XXX two-dimensional process
% nu, initial vola, nubar long run mean vola
%nu         = V0; % XXX
c          = K;
nubar      = thetav;
%sigmav     = sigmav;
sigmacy    = sigmas;
mucy       = mus;
mucv       = muv;
%rhoj       = rhoj;
rhobar     = rho;
%kappav     = kappav;
zetabar    = q;
lambdac    = lambda;

% not needed for SVJJ, ie SVCJ
sigmay    = 0;
lambday   = 0;
lambdav   = 0;
muv        = 0;
muy        = 0;

% specific to SVCJ model
lambdabar  = lambday + lambdav + lambdac;
r          = r;

mubar = thetafunc(1,0) - 1;

function retval = alphabar(tau,u)
thetainter = lambdabar^-1*(lambday*fy(u,tau)+lambdav*fv(u,tau)+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))); % XXX minus was missing
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);
retval =  (gamma-b)./((gamma-b)+muv*a)*tau-(2*muv*a)./(gamma.^2-(b-muv*a).^2)...
.*log(1-((gamma+b)-muv*a)./(2*gamma).*(1-exp(-gamma*tau))); % XXX equation had a number of mistakes
retval(a==0) = (gamma-b)./((gamma-b)+muv*a)*tau; % XXX take care of special case
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_;
retval(a==0) = (gamma-b)./((gamma-b).*c_+mucv*a)*tau; % XXX take care of special case
end

function retval = thetafunc(c1, c2)
retval = lambdabar^-1*(lambday*thetayfunc(c1)+lambday*thetavfunc(c2)+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

% equation 4.2, p. 21
function retval = psi(u, y, v, t, T)
retval =  exp(alphabar(T-t,u)+u*y+ betabar(T-t,u)*v);
end

function retval = psichi(u, X0, t, T)
y = X0(1);
v = X0(2);
retval = psi(u, y, v, t, T);
end

function retval = Gab(a,b, y, X0_, T)
integratefunc = @(nu_) imag(psichi(complex(a,nu_*b),X0_, 0, T) .* exp(complex(0,-nu_*y))) ./ nu_; % XXX
retval = psichi(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 = 0; % XXX values for European call option
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

-
I believe the formula is now correct as it stands. I forgot mu as well which should be mu = 0.2380. It seems to me that the conversion between mu and mus (mus = log((1+mu)*(1-rhoj*muv))-0.5*(sigmas^2);) isnt applicable here. Thanks a lot and I will upload the code in due time to my github page. – jcfrei Jul 21 '14 at 23:59

I would suggest that you use a more 'modern' method to recover option prices from characteristic functions.

The approach of this papers (for practical calculations of option prices) is somewhat outdated. The backbone of affine models (such as SVJJ) is the characteristic function $\psi(u)$ of the log-price distribution, which is known in closed form.

The Duffie/ Pan/ Singleton paper, in the footsteps of Heston and Bakshi/ Madan, expresses the option price as $$C = S_0\ \Pi_1 - K \exp(-rT)\ \Pi_2$$ This is an intuitive form that resembles Black/ Scholes, with the difference that $\Pi_{1,2}$ are not cumulative Gaussians, but rather cumulative densities of 'bespoke' distributions. They can be expressed as integrals of the characteristic function, which is what the paper does.

Expressing the option price a'la Black/ Scholes, with a decomposition in terms of Delta ($\Pi_1$) and the probability of ending up in-the-money ($\Pi_2$) is intutive and well suited for teaching or academic research. However, this approach is not numerically stable, as the integrand can be very oscillatory and decaying very slowly. This was illustrated in Carr/ Madan, where they propose to invert the characteristic function and recover the option price directly using FFT. Chourdakis (pdf here) proposes a further refinement that improves robustness (fractional FFT).

I am sure that there has been susequent research that might offer better procedures for the inversion of the characteristic function (some names to Google include Lewis, Lipton, Fang/ Oosterle, Levendorskii). But I don't have first hand experience of these methods. I have implemented FFT/ FRFT in similar settings very successfully, and can vouch for them.

-