Take the 2-minute tour ×
Quantitative Finance Stack Exchange is a question and answer site for finance professionals and academics. It's 100% free, no registration required.

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

share|improve this question
    
Whats the value of lambda for the example? –  Theja Jul 17 at 11:16
    
oh, forgot that one: lambda=4; –  jcfrei Jul 17 at 11:31
1  
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 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 at 10:10
1  
Thanks, unfortunately I won't be looking at it soon but this makes it more clear. –  Bob Jansen Jul 21 at 10:48

2 Answers 2

up vote 3 down vote accepted
+50
  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
share|improve this answer
    
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 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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.