Aplicación del modelo Kou en Python

Question

Hey intento implementar el modelo Kou en Python. Este es mi código:

def Hh(n,x):
    if n<-1: return 0
    elif n==-1:
        return np.exp(-x**2/2)
    elif n==0:
        return math.sqrt(2*np.pi)*scs.norm.cdf(-x)
    else:
        return (Hh(n-2,x)-x*Hh(n-1,x))/n

def P(n,k):
    if k<1 or n<1: return 0
    elif n==k:
        return p**n
    else:
        suma=0
        i=k
        while i<=n-1:
            suma=suma+sc.special.binom(n-k-1, i-k)*sc.special.binom(n, i)*(n_1/(n_1+n_2))**(i-k)*(n_2/(n_1+n_2))**(n-i)*p**i*q**(n-i)
            i+=1
        return suma

def Q(n,k):
if k<1 or n<1: return 0
elif n==k:
    return q**n
else:
    suma=0
    i=k
    while i<=n-1:
        suma=suma+sc.special.binom(n-k-1, i-k)*sc.special.binom(n, i)*(n_1/(n_1+n_2))**(n-i)*(n_2/(n_1+n_2))**(i-k)*p**(n-i)*q**i
        i+=1
    return suma

def Pi(n):
    return (np.exp(-lam*T)*(lam*T)**n)/math.factorial(n)

def I(n,c,a,b,d):
    if b>0 and a!=0:
        suma=0
        i=0
        while i<=n:
            suma=suma+(b/a)**(n-i)*Hh(i,b*c-d)
            i+=1
        return -(np.exp(a*c)/a)*suma+(b/a)**(n+1)*(np.sqrt(2*np.pi)/b)*np.exp((a*d/b)+(a**2/(2*b**2)))*scs.norm.cdf(-b*c+d+a/b)
    elif b<0 and a<0:
        suma=0
        i=0
        while i<=n:
            suma=suma+(b/a)**(n-i)*Hh(i,b*c-d)
            i+=1
        return -(np.exp(a*c)/a)*suma-(b/a)**(n+1)*(np.sqrt(2*np.pi)/b)*np.exp((a*d/b)+(a**2/(2*b**2)))*scs.norm.cdf(b*c-d-a/b)
    else: return 0

def Y(mu,sigma,lam,p, n_1, n_2, a ,T):
    n=1
    suma1=0
    suma2=0
    while n<=10:
        k=1
        suma_1=0
        while k<=n:
            suma_1=suma_1+P(n,k)*(sigma*np.sqrt(T)*n_1)**k*I(k-1,a-mu*T,-n_1, -1/(sigma*np.sqrt(T)), -sigma*n_1*np.sqrt(T))
            k+=1
        suma1=suma1+Pi(n)*suma_1
        n+=1
    n=1
    while n<=10:
        k=1
        suma_2=0
        while k<=n:
            suma_2=suma_2+Q(n,k)*(sigma*np.sqrt(T)*n_2)**k*I(k-1,a-mu*T,n_2, 1/(sigma*np.sqrt(T)), -sigma*n_2*np.sqrt(T))
            k+=1
        suma2=suma2+Pi(n)*suma_2
        n+=1
    return np.exp((sigma*n_1)**2*T/2)/(sigma*np.sqrt(2*np.pi*T))*suma1+np.exp((sigma*n_2)**2*T/2)/(sigma*np.sqrt(2*np.pi*T))*suma2+Pi(0)*scs.norm.cdf(-(a-mu*T)/(sigma*np.sqrt(T)))

def Kou(r,sigma,lam,p,n_1,n_2,S_0,K,T):
    zeta=p*n_1/(n_1-1)+(1-p)*n_2/(n_2+1)-1
    lam2=lam*(zeta+1)
    n_12=n_1-1
    n_22=n_2+1
    p2=p/(1+zeta)*n_1/(n_1-1)
    return S_0*Y(r+1/2*sigma**2-lam*zeta,sigma,lam2,p2,n_12,n_22,math.log(K/S_0),T)-K*np.exp(-r*T)*Y(r-1/2*sigma**2-lam*zeta,sigma,lam,p,n_1,n_2,math.log(K/S_0),T)

y utilizar los siguientes datos (de Kou 2002)

S_0=100
sigma=0.16
r=0.05
lam=1
n_1=10
n_2=5
T=0.5
K=98
p=0.4

Por desgracia, mi resultado es 6,25, cuando en Kou debería ser 9,14732. ¿Puede alguien comprobar si mi código está bien, o si alguien tiene el código del modelo Kou, sería capaz de dar varios valores para diferentes funciones para que pueda comprobar en qué función tengo error?

Answer 1

¿Ayuda esta función de matlab?

La expresión para $\Upsilon$ ( Y ) se divide en tres partes. La suma infinita se trunca después de 10 iteraciones ( bound ). Ya hemos comparado nuestros resultados para PFunction , QFunction et IFunction .

function Y = Upsilon(x,T,mu,sigma,lambda,etaplus,etaminus,p,q)

    bound = 10;

    pi0 = exp(-lambda*T);
    pin = exp(-lambda*T) .*(lambda*T).^(1:bound)./factorial(1:bound);

    sump1 = zeros(bound,1);
    sumq1 = zeros(bound,1);

    for n=1:bound
        sump2 = zeros(n,1);
        sumq2 = zeros(n,1);
        for k=1:n
            sump2(k) = PFunction(n,k,p,q,etaplus,etaminus) * (sigma*sqrt(T)*etaplus)^k * IFunction(-etaplus,-1/(sigma*sqrt(T)),x-mu*T,-sigma*etaplus*sqrt(T),k-1);
            sumq2(k) = QFunction(n,k,p,q,etaplus,etaminus) * (sigma*sqrt(T)*etaminus)^k * IFunction(etaminus,1/(sigma*sqrt(T)),x-mu*T,-sigma*etaminus*sqrt(T),k-1);
        end
       sump1(n) = pin(n)*sum(sump2);
       sumq1(n) = pin(n)*sum(sumq2);
    end

    Y(1) = exp((sigma*etaplus)^2*T/2)/(sigma*sqrt(2*pi*T)) * sum(sump1);
    Y(2) = exp((sigma*etaminus)^2*T/2)/(sigma*sqrt(2*pi*T)) * sum(sumq1);
    Y(3) = pi0*normcdf(-(x-mu*T)/(sigma*sqrt(T)));
    Y = sum(Y);

end