3
$\begingroup$

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

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

2 Answers 2

4
+50
$\begingroup$
  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
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – jcfrei
    Jul 21, 2014 at 23:59
  • $\begingroup$ I have the impression that the code is not fully correct, with regards to the last statements to compute the price. I have written an implementation based on the above characteristic function in Julia at github.com/jherekhealy/CharFuncPricing.jl where I reproduce the price in published paper. I however use a different pricing formula. $\endgroup$
    – jherek
    Dec 5, 2023 at 14:48
2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Can You Send me the Paper by Carr and Madan. If you have some Code Can you share $\endgroup$ Mar 13, 2020 at 12:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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