Diagnosing Price Differences in Kou Jump-Diffusion Option Valuation

# Kou model - can't reproduce prices of European Option from Toivanen and Forsyth


# Kou model - can't reproduce prices of European Option from Toivanen and Forsyth












I have implemented the Kou option model for pricing vanilla option. I have checked that my program can replicate the price of the option in the original paper of 2002. However, when I use it to price other options, I observe deviation from the prices given by Forsyth (2004) "Robust Numerical Methods for Contingent Claims under Jump Diffusion Processes" or Toivanen "Numerical Valuation of European and American Options under Kou’s Jump-Diffusion Model".

For instance, they price a call at the money at 3.97 whereas with my code, I found 3.69. I attached my code below for your reference.

Just note that I have re-used multiple sources but the principal one was a previous discussion on this forum!

Thanks in advance.

I wonder whether I forgot something in my code.

```
import numpy as np, scipy as sp, math

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)*sp.stats.norm.cdf(-x)
    else:
        return (Hh(n-2,x)-x*Hh(n-1,x))/n

def P(n, k, p, q, eta_1, eta_2):
    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 + sp.special.binom(n-k-1, i-k) * sp.special.binom(n, i) * (eta_1/(eta_1 + eta_2))**(i-k)*(eta_2/(eta_1 + eta_2))**(n-i)*p**i*q**(n-i)
            i+=1
        return suma

def Q(n, k, p, q, eta_1, eta_2):
    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 + sp.special.binom(n-k-1, i-k) * sp.special.binom(n, i)*(eta_1/(eta_1+eta_2))**(n-i)*(eta_2/(eta_1+eta_2))**(i-k)*p**(n-i)*q**i
            i +=1
        return suma

def I(n,c,a,b,d):
    if n < 0:
        print(n)
    if b>0 and a!=0:
        suma=0
        i=0
        while i<=n:
            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 + 0.5 * a**2 / b**2) * sp.stats.norm.cdf(-b*c+d+a/b)
    elif b<0 and a<0:
        suma=0
        i=0
        while i<=n:
            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 + 0.5 * a**2 / b**2) * sp.stats.norm.cdf(b*c-d-a/b)
    else: 
        return 0

def PSI(mu, Sigma, Lambda, pp, qq, eta_1, eta_2, arg, Time):
    bound=15
    Aa = 0.
    Bb = 0.
    Cc = 0.
    
    pin=np.zeros(bound+1)

    sump1=0.
    sumq1=0.
    
    for i in range(1,bound+1):
        pin[i]=np.exp(-Lambda*Time)*(Lambda*T)**i/math.factorial(i)

    for n in range(1, bound):
        sump2=0.
        sumq2=0.
        
        for k in range(1, n+1):
            sump2 += P(n,k,pp,qq, eta_1, eta_2) * (Sigma * np.sqrt(Time) * eta_1)**k * I(k-1, arg-mu*Time, -eta_1, -1/(Sigma*np.sqrt(Time)), -Sigma * eta_1 * np.sqrt(Time))
            sumq2 += Q(n,k,pp,qq, eta_1, eta_2) * (Sigma * np.sqrt(Time) * eta_2)**k * I(k-1, arg-mu*Time, eta_2,  1/(Sigma*np.sqrt(Time)), -Sigma * eta_2 * np.sqrt(Time))
        
        sump1+=pin[n] * sump2
        sumq1+=pin[n] * sumq2

    Aa = np.exp((Sigma*eta_1)**2*Time/2)/(Sigma*np.sqrt(2*np.pi*Time)) * sump1
    Bb = np.exp((Sigma*eta_2)**2*Time/2)/(Sigma*np.sqrt(2*np.pi*Time)) * sumq1
    Cc = np.exp(-Lambda*Time) * sp.stats.norm.cdf(-(arg-mu*Time)/(Sigma*np.sqrt(Time)))
    
    return Aa + Bb + Cc

def Kou(r,sigma,Lambda,pp,eta_1,eta_2,S_0,K,Time):
    zeta = p * eta_1/(eta_1 - 1.) + (1. - p) * eta_2/(eta_2 + 1.) - 1.
    lam2 = Lambda * (zeta + 1.) 
    eta_12 = eta_1 - 1.
    eta_22 = eta_2 + 1.
    p2   = pp/(1. + zeta) * eta_1/(eta_1 - 1.)   

    return S_0 * PSI(r + 1/2 * sigma**2 - Lambda * zeta, sigma, lam2, p2, 1. - p2, eta_12, eta_22, math.log(K/S_0), Time) - K * np.exp(-r*T) * PSI(r - 1/2 * sigma**2 - Lambda * zeta, sigma, Lambda, pp, 1. - pp, eta_1, eta_2, math.log(K/S_0), Time)
    

S0 = 110.
r = 0.05
K = 100.
Sigma = 0.15
T = 0.25
Lambda_ = .1

Alpha1 = 3.0465
Alpha2 = 3.0775
                                             
p = 0.3445

print(Kou(r,Sigma,Lambda_,p,Alpha1,Alpha2,S0,K,T))

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
q =1.-p

Alpha2

print(Kou(r,sigma,lam,p,n_1,n_2,S_0,K,T))
```