Double jumps stochastic volatility model (SVCJ, Duffie et al, 2000) - characteristic function for VIX

Currently I am working at my master's degree paper where I want to evaluate VIX options using stochastic volatility jump models.I got some MATLAB codes for the SVCJ model for the S&P, but as the VIX dynamic is different I can not use them. Does everybody know where I could find the MATLAB code for this characteristic function, or, could somebody help me implement it?

Thank you, Alex

This is what I've done, but the prices returned by the function are huge:

function price = CallDuffie(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; %mus = 0.02; %muJ = -0.04; %sigmav=0.1; % V0 = 0.25; % % % mus = log((1+mu)(1-rhojmuv))-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 = alpha(tau,u) thetainter = lambdabar^-1*(lambdayfy(u,tau)+lambdavfv(u,tau)+lambdac*fc(u,tau)); retval = -2*kappav*thetainter/sigmav^2*log(1-sigmav^2/2*kappav*(1-exp(-kappav*tau)*u)); end

function retval = beta(tau,u)
retval =2*kappav*u/(sigmav^2*u+(2*kappav-sigmav^2*u)*exp(kappav*tau));
end

function retval = gama(tau,u)
retval = 2*mucv/(2*kappav*mucv-sigmav^2)*log(1+((2*kappav*mucv-sigmav^2)*u)/2*kappav*(1-mucv*u)*(1-exp(-kappav*tau)));
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 = sigmavrhobaru-kappav; % gamma = sqrt(b.^2+asigmav^2); % retval = -r*tau+(r-zetabar)utau-kappavnubar... % (((gamma+b)/sigmav^2)*tau + ... % (2/sigmav^2)*log(1-((gamma+b)./(2*gamma)).(1-exp(-gammatau)))); %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


% equation 4.2, p. 21 function retval = psi(u, y, v, t, T) retval = exp(alpha(T-t,u)+beta(T-t,u)v+lambdabargama(T-t,u)); 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

Thank you, Alexenter preformatted text here

• The VIX is derived from observed call / put options on the index, no? If you want to write the Cf on the latent volatility, you might need the augmented AJD piece from the DPS2000 paper. May 8 '20 at 19:42
• Typo: I meant the „extended“ transform in their paper May 9 '20 at 8:22
• Have you tried implementing this yourself? Posting what you have done so far would be helpful. May 9 '20 at 13:49
• I am trying to adapt this matlab code quant.stackexchange.com/questions/14049/… with the closed form aproximation from this article gen.lib.rus.ec/scimag/… This is what I've done, but the prices returned by the function are negative:
– Alex
May 11 '20 at 16:44