Pricing SVCJ Options with Stable Characteristic Function Inversion

# SVCJ (SVJJ) Duffie et. al Model implementation in Matlab


# SVCJ (SVJJ) Duffie et. al Model implementation in Matlab












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

## Answer by pbr142 (score 4, accepted)

https://quant.stackexchange.com/a/14114

- 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.

- 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
```

## Answer by Kiwiakos (score 2)

https://quant.stackexchange.com/a/14115

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.