Debugging a Kou Double-Exponential Jump-Diffusion Option Pricer

# Kou model implementation Python


# Kou model implementation Python












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 by Kevin (score 3, accepted)

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

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