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