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

Question

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