Kou model implementation Python

Question

Hey I try to implement Kou model in Python. This is my code:

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)

and use following data (from 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

Unfortunately, my result is 6.25, when in Kou should be 9.14732. Can someone check if my code is OK, or if someone has code of Kou model, would he be able to give several values ​​for different functions so I can check in which function I have mistake?

Answer 1

Does this matlab function help?

The expression for $\Upsilon$ (Y) is split up in three parts. The infinite sum is truncated after 10 iterations (bound). We already compared our results for PFunction, QFunction and 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