SVCJ (SVJJ) Duffie et. al Modelo de implementación en Matlab

Question

Estoy tratando de implementar el mencionado SVCJ modelo por Duffie et al en MATLAB. hasta el momento sin éxito. Se supone que el precio de la vainilla (europea) de la llamada . los parámetros de siempre, el precio esperado es: ~6.0411. La función arroja el error:

Advertencia: Alcanzado el límite en el número máximo de intervalos en uso. Aproximado obligado en caso de error es 1.1 e+03. La integral puede no existir, o puede ser difícil para aproximar numéricamente a la solicitada la exactitud.

El uso de quadv en lugar de la integral no producir el resultado esperado, sino que simplemente devuelve Nan.

    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

Funciones en el papel:

thetainer: p. 23

alpha0: p.22

alphabar: p.22

betabar: p.22

el año fiscal,fv,fc: p. 23

thetayfunc, thetacfunc, thetavfunc: p. 22

psi: p.21

Gab: p. 13

no una función, pero dependiendo de los derivados de tipo:

GbTplusdminusd,GbTminusd, precio: p. 18

lambdabar = lambda, de acuerdo a: p. 24

Answer 1

1. Sobre el problema de integración: El integrando es muy oscilatorio, y la adaptación de cuadratura de Matlab no manejar tales integrands muy bien. En general, yo recomendaría Mathematica cuando Matlab estándar de procedimientos no funcionan bien. En este caso, una Levin-tipo de método funcionaría mucho mejor. La razón por la que quadv produce valores NaN es debido a la singularidad en 0. Usted tendría que comenzar la integración en un pequeño valor positivo.
2. En cualquier caso, hay un par de errores en su código. He intentado corregir el código. La integración funciona bien ahora, pero yo todavía no tienen el valor que estaban esperando (yo he usado t=0, q=0, mu=0 para los parámetros no especificados; mu es necesario en el mus de la fórmula). He añadido XXX como comentarios en las líneas que he cambiado. También, he hecho todas las operaciones de vector con valores para acelerar la integración.

   ---

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 2

Me gustaría sugerir que el uso de una forma más "moderna" método para recuperar la opción de precios de funciones características.

El enfoque de este (documentos para la práctica de los cálculos de la opción de precios) es un poco anticuado. La columna vertebral de los modelos afines (tales como SVJJ) es la característica de la función $\psi(u)$ de la log-precio de distribución, que es conocido en forma cerrada.

El Duffie/ Pan/ Singleton de papel, en los pasos de Heston y Bakshi/ Madan, expresa el precio de la opción como $$C = S_0\ \Pi_1 - K \exp(-rT)\ \Pi_2$$ Este es un formato intuitivo y que se asemeja a Black/ Scholes, con la diferencia de que $\Pi_{1,2}$ no son acumulables Gaussianas, sino acumulativa densidades de 'a medida' de las distribuciones. Que puede ser expresada como la integral de la función característica, que es lo que el papel que hace.

Expresar el precio de la opción a la Black/ Scholes, con una descomposición en términos de Delta ($\Pi_1$) y la probabilidad de acabar en-el-dinero ($\Pi_2$) es intutive y muy adecuado para la enseñanza o la investigación académica. Sin embargo, este enfoque no es numéricamente estable, como el integrando puede ser muy oscilatoria y decae muy lentamente. Esto se ilustra en Carr/ Madan, donde se proponen invertir la función característica y recuperar el precio de la opción directamente el uso de la FFT. Chourdakis (pdf aquí) propone un refinamiento adicional que mejora la robustez (fraccional FFT).

Estoy seguro de que no ha sido susequent de investigación que podría ofrecer mejores procedimientos para la inversión de la función característica (algunos nombres a Google incluyen Lewis, Lipton, Fang/ Oosterle, Levendorskii). Pero no tengo experiencia de primera mano de estos métodos. He implementado FFT/ FRFT en una configuración similar con mucho éxito, y puede responder por ellos.