# Kou model — solving PIDE for European and American options in Python

Toivanen proposed$$^\color{magenta}{\star}$$ a method to solve the partial integro-differential equation (PIDE) with a numerical scheme based on Crank-Nicolson. In particular, he proposed an algorithm to estimate quickly the integral, saving a significant amount of time.

I have tried to implement it in Python (and NumPy) but I found out that the price of an European option at-the-money deviates from the price in the original paper ($$3.7755$$ versus $$3.97347$$). Option features are detailed in my code below and the double exponential distribution as well.

from IPython.display import display
import sympy
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

def display_matrix(m):
display(sympy.Matrix(m))

#Option characteristics
S0 = 100.; r = 0.05; K = 100.; Sigma = 0.15; T = 0.25

#Diffusion parameters
Lambda = .1; Alpha1 = 3.0465; Alpha2 = 3.0775; p = 0.3445; q = 1. - p
Zeta = p * Alpha1 /(Alpha1 - 1.) + q * Alpha2 /(Alpha2 + 1.) - 1.

#Grid Parameters
M = 5000 ; N = 400; Beta = 10.3; a = 4.; Gamma = 0.45

#Vector integral term
AMinus = np.zeros(M+1)
APlus  = np.zeros(M+1)

# Spatial Grid
index_x = np.linspace(start=0, stop=M, num=M+1)
gridX = K * (1 + np.sinh(Beta *(index_x/M - Gamma))/np.sinh(Beta*Gamma))
gridDelta =  gridX[1:] - gridX[:-1]

#Time Grid
def fillTimeGrid(k):
if k < N:
return T * (a**(-k/(N+2.))-1.)/(1./a - 1.)              , 0.5
else:
return T * (a**(-(k + N )/(2.*N+4.)) - 1.)/(1./a - 1.)  , 1.

index_t = np.linspace(start=0, stop=N+4, num=N+5)
gridT, Theta = np.vectorize(fillTimeGrid)(index_t)

gridDeltaT = gridT[1:] - gridT[:-1]

def multJ(v: np.ndarray, kprime: int) -> np.ndarray:

temp1 = 0.
temp2 = 0.
# Initialization
AMinus[0] = q / (1. + Alpha2) * (v0(gridT[kprime]) + Alpha2 * v[1])

for i in range(0,M-2):
temp1 = q * gridX[i+1]**(-Alpha2) / ((1. + Alpha2) * gridDelta[i]) * (gridX[i+1]**(1. + Alpha2) - (gridX[i+1] + Alpha2 * gridDelta[i]) * gridX[i]  **(Alpha2)) * v[i-1]
temp2 = q * gridX[i+1]**(-Alpha2) / ((1. + Alpha2) * gridDelta[i]) * (gridX[i  ]**(1. + Alpha2) - (gridX[i]   - Alpha2 * gridDelta[i]) * gridX[i+1]**(Alpha2)) * v[i]

AMinus[i+1] = (gridX[i] / gridX[i+1])**(Alpha2) * AMinus[i] + temp1 + temp2

# Initialization
temp1 = 0.
temp2 = 0.

APlus[M]   = 0.
APlus[M-1] = p * gridX[M-1]**(Alpha1) / ((Alpha1 - 1.) * gridDelta[M-1]) * (gridX[M]**(1. - Alpha1) - (gridX[M] - Alpha1 * gridDelta[M-1]) * gridX[M-1]**(-Alpha1)) * v[M-2]

for i in range(M-1,1,-1):
temp1 = p * gridX[i-1]**(Alpha1)/((Alpha1 - 1.) * gridDelta[i-1]) * (gridX[i  ]**(1. - Alpha1) - (gridX[i]     - Alpha1 * gridDelta[i-1]) * gridX[i-1]**(-Alpha1)) * v[i-2]
temp2 = p * gridX[i-1]**(Alpha1)/((Alpha1 - 1.) * gridDelta[i-1]) * (gridX[i-1]**(1. - Alpha1) - (gridX[i-1]   + Alpha1 * gridDelta[i-1]) * gridX[i  ]**(-Alpha1)) * v[i-1]
APlus[i-1] = (gridX[i-1] / gridX[i])**(Alpha1) * APlus[i] + temp1 + temp2

return -Lambda * (AMinus[1:M] + APlus[1:M])

def v0(tau: float) -> float:
return K * np.exp(-r*(T-tau))

#initialization of the v - Final value
v = np.maximum(0., K - gridX[1:-1])
v_previous = v; g = v

#D Matrix
SigmaE = np.maximum(Sigma * Sigma, (r - Lambda * Zeta) * gridDelta[1:]/gridX[1:-1], -(r - Lambda * Zeta) * gridDelta[:-1]/gridX[1:-1])
DDMinus = gridX[1:-1] * (-SigmaE * gridX[1:-1] + (r - Lambda * Zeta) * gridDelta[1:]) / (gridDelta[:-1] * (gridDelta[1:] + gridDelta[:-1]))
DDPlus  = gridX[1:-1] * (-SigmaE * gridX[1:-1] - (r - Lambda * Zeta) * gridDelta[:-1]) / (gridDelta[1:] * (gridDelta[1:] + gridDelta[:-1]))
DDCenter = r + Lambda - DDMinus - DDPlus

D = np.diag(DDMinus[1:M], -1) + np.diag(DDCenter[0:M-1]) + np.diag(DDPlus[0:M-2], 1)
if M <= 100: display_matrix(D)

for k in range(N+3, -1,- 1):
# Equation (32), (33) and (34)
M1 = np.eye(M-1) + Theta[k] * gridDeltaT[k] * D
M2 = np.eye(M-1) - (1. - Theta[k]) * gridDeltaT[k] * D
M1_inv = np.linalg.inv(M1)

b = np.matmul(M2, v)
d = 1.
v_previous = v
v1 = v

#Jv and solve the equation (38) in the paper
normVectorb = np.linalg.norm(b)
while(d > 1.0e-10*normVectorb):
v1 = b - Theta[k] * gridDeltaT[k] * multJ(v_previous, k)
v1 = np.matmul(M1_inv, v1)
d  = np.linalg.norm(Theta[k] * gridDeltaT[k] * multJ(v1, k) -Theta[k] * gridDeltaT[k] * multJ(v_previous, k) )
v_previous = v1

#Update Price vector
v = v1

oPriceCall = np.interp(S0, gridX[0:M-1], v[0:M-1]) + S0 - K*np.exp(-r*T)
oPricePut  = np.interp(S0, gridX[1:-1], v)
print(f'The call Price with a spot value {S0} at the money is {oPriceCall} while the put price is {oPricePut}.')


The help of anyone familiar with this model and its implementation will be much appreciated because I have been looking into this problem for days (and nights). Any reference to any papers or books detailing the implementation will be helpful too.

EDIT Instead of a Crank-Nicolson scheme, I have tried a full implicit scheme. Then, I managed to retrieve the price of the price of the call and put option given in their paper. I am still not sure why. I will keep using this setting for developing the American option pricing in this framework.

Thanks to everyone!

$$\color{magenta}{\star}$$ Jari Toivanen, Numerical valuation of European and American options under Kou's jump-diffusion model.

• Please can we close this question since I answer myself to my question? Dec 2, 2023 at 6:30

The issue I described in my initial question is linked to the integral term. In the paper, this term is multiply by $$\theta \Delta \text{t}$$ but this is only the "implicit" part of the Crank-Nicolson. What is missing is the explicit part which should be multiplied by $$( 1- \theta) \Delta \text{t}$$. It appears that this term is missing in the paper and then in my implementation.

I attached a corrected code below:

from IPython.display import display
import sympy
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

def display_matrix(m):
display(sympy.Matrix(m))

#Option characteristics
S0 = 110.; r = 0.05; K = 100.; Sigma = 0.15; T = 0.25

#Diffusion parameters
Lambda = .1; Alpha1 = 3.0465; Alpha2 = 3.0775; p = 0.3445; q = 1. - p
Zeta = p * Alpha1 /(Alpha1 - 1.) + q * Alpha2 /(Alpha2 + 1.) - 1.

#Grid Parameters
M = 1000 ; N = 1000; Beta = 12.3; a = 4.; Gamma = 0.45

#Vector integral term
AMinus = np.zeros(M+1)
APlus  = np.zeros(M+1)

#Bounday limit
f = np.zeros(M-1)
f[0] = 1.

# Spatial Grid
index_x = np.linspace(start=0, stop=M, num=M+1)
gridX = K * (1 + np.sinh(Beta *(index_x/M - Gamma))/np.sinh(Beta*Gamma))
gridDelta =  gridX[1:] - gridX[:-1]

#Time Grid
def fillTimeGrid(k):
if k < N:
return T * (a**(-k/(N+2.))-1.)/(1./a - 1.)              , 0.5
else:
return T * (a**(-(k + N )/(2.*N+4.)) - 1.)/(1./a - 1.)  , 1.

index_t = np.linspace(start=0, stop=N+4, num=N+5)
gridT, Theta = np.vectorize(fillTimeGrid)(index_t)
gridDeltaT = gridT[1:] - gridT[:-1]

def multJ(v: np.ndarray, kprime: int) -> np.ndarray:

#Vector integral term
AMinus = np.zeros(M+1)
APlus  = np.zeros(M+1)

# Initialization
AMinus[0] = q / (1. + Alpha2) * (v0(kprime)+ Alpha2 * v[0])

for i in range(0,M-1):
temp1 = (gridX[i+1]**(1. + Alpha2) - (gridX[i+1] + Alpha2 * gridDelta[i]) * gridX[i]  **(Alpha2)) * v[i-1]
temp2 = (gridX[i  ]**(1. + Alpha2) - (gridX[i]   - Alpha2 * gridDelta[i]) * gridX[i+1]**(Alpha2)) * v[i]
AMinus[i+1] = (gridX[i] / gridX[i+1])**(Alpha2) * AMinus[i] + q * gridX[i+1]**(-Alpha2) / ((1. + Alpha2) * gridDelta[i]) * (temp1 + temp2)

APlus[M]   = 0.
APlus[M-1] = p * gridX[M-1]**(Alpha1) / ((Alpha1 - 1.) * gridDelta[M-1]) * (gridX[M]**(1. - Alpha1) - (gridX[M] - Alpha1 * gridDelta[M-1]) * gridX[M-1]**(-Alpha1)) * v[M-2]

for i in range(M-1,1,-1):
temp1 = (gridX[i  ]**(1. - Alpha1) - (gridX[i]     - Alpha1 * gridDelta[i-1]) * gridX[i-1]**(-Alpha1)) * v[i-2]
temp2 = (gridX[i-1]**(1. - Alpha1) - (gridX[i-1]   + Alpha1 * gridDelta[i-1]) * gridX[i  ]**(-Alpha1)) * v[i-1]
APlus[i-1] = (gridX[i-1] / gridX[i])**(Alpha1) * APlus[i] + p * gridX[i-1]**(Alpha1)/((Alpha1 - 1.) * gridDelta[i-1]) * (temp1 + temp2)

return -Lambda * (AMinus[1:M] + APlus[1:M])

def v0(k) -> float:
return np.maximum(0., K * np.exp(-r*(T-gridT[k]))  )

#initialization of the v - Final value
v = np.maximum(0., K - gridX[1:-1])
v_previous = v; g = v
if M <= 100: display_matrix([gridX[1:-1].round(4), v.round(4)])

#D Matrix
SigmaE = np.maximum(Sigma * Sigma, (r - Lambda * Zeta) * gridDelta[1:]/gridX[1:-1], -(r - Lambda * Zeta) * gridDelta[:-1]/gridX[1:-1])
#SigmaE = Sigma * Sigma
DDMinus = gridX[1:-1] * (-SigmaE * gridX[1:-1] + (r - Lambda * Zeta) * gridDelta[1:]) / (gridDelta[:-1] * (gridDelta[1:] + gridDelta[:-1]))
DDPlus  = gridX[1:-1] * (-SigmaE * gridX[1:-1] - (r - Lambda * Zeta) * gridDelta[:-1]) / (gridDelta[1:] * (gridDelta[1:] + gridDelta[:-1]))
DDCenter = r + Lambda - DDMinus - DDPlus

D = np.diag(DDMinus[1:M], -1) + np.diag(DDCenter[0:M-1]) + np.diag(DDPlus[0:M-2], 1)
if abs(DDMinus[0]) > 0:print(f'DD is not null {DDMinus[0], 4}')
if M <= 100: display_matrix(D.round(4))

for k in range(N+3, -1,- 1):
# Equation (32), (33) and (34)
M1 = np.eye(M-1) + Theta[k] * gridDeltaT[k] * D
M2 = np.eye(M-1) - (1. - Theta[k]) * gridDeltaT[k] * D
M1_inv = np.linalg.inv(M1)

b = np.matmul(M2, v)

#Boundary limit
b = b - DDMinus[0] * gridDeltaT[k] * (Theta[k] * v0(k) + (1. - Theta[k]) * v0(k+1)) * f

d = 1.
v_previous = v
v1 = v

#Jv and solve the equation (38) in the paper
normVectorb = np.linalg.norm(b)
while(d > 1.0e-10*normVectorb):

v1 = b - gridDeltaT[k] * (Theta[k] * multJ(v_previous, k) + (1. - Theta[k]) * multJ(v, k+1))
v1 = np.matmul(M1_inv, v1)
d  = np.linalg.norm(v1 - v_previous)
v_previous = v1

#Update Price vector
v = v1

print(k," ", gridDeltaT[k], " ", gridT[k],",", np.interp(S0, gridX[0:M-1], v[0:M-1]) + S0 - K*np.exp(-r*(T-gridT[k]))," ", DDMinus[0])

oPriceCall = np.interp(S0, gridX[1:-1], v) + S0 - K * np.exp(-r*T)
oPricePut  = np.interp(S0, gridX[1:-1], v)
print(f'The call Price with a spot value {S0} at the money is {oPriceCall} while the put price is {oPricePut}.')

fig, ax = plt.subplots()
ax.scatter(gridX[0:M-2], v[0:M-2], s=0.5)
ax.scatter(gridX[0:M-2], g[0:M-2], s=0.5)
plt.show()

$$$$
`