Correcting Crank-Nicolson Integration in Kou Jump-Diffusion Pricing
Article Quant Q&A · Author: pierrot
Summary
This note describes a numerical approach to pricing European and American options under Kou’s double-exponential jump-diffusion model. The implementation uses a Crank-Nicolson time scheme for the partial integro-differential equation, a nonuniform spatial grid, and a recurrence intended to evaluate the jump integral efficiently. The author initially reports a mismatch between an at-the-money European option price and the cited paper, then experiments with a fully implicit scheme that reproduces the paper’s stated price.
The follow-up identifies the discrepancy as a missing explicit contribution from the integral term in the Crank-Nicolson update. The corrected update weights both the implicit and explicit evaluations of that term according to the time-scheme parameter, and the revised code also handles a boundary contribution. This is a specific implementation diagnosis, not a general validation of the full code. The note supplies no convergence study or independent benchmark beyond the reported comparison, and its lengthy code is tied to particular model and grid settings.
Key ideas
- Kou jump-diffusion option pricing can be formulated as a partial integro-differential equation.
- The described scheme combines Crank-Nicolson time stepping with a nonuniform spatial grid.
- The jump integral requires both implicit and explicit contributions in the Crank-Nicolson update.
- Omitting the explicit integral term can lead to a discrepancy from a cited benchmark price.
- The reported price match is not accompanied by a convergence study or broader validation.
Tags
Full text
# Kou model — solving PIDE for European and American options in Python
# 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.
## Answer by pierrot (score 1)
https://quant.stackexchange.com/a/77578
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()
```
```Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.