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

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!

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

• This is a link to the earlier disscussion (2020) quant.stackexchange.com/questions/57498/… Nov 23, 2023 at 11:57
• Thanks, this is one of the discussion I used to write the code. After posting my code, I realized it is actually working. Nov 23, 2023 at 15:17
• Do you want to close the question ? Nov 23, 2023 at 15:20
• Yes please, I don't need any help on this. Nov 24, 2023 at 1:37