Stability Tests done on Numerical Approaches under Option Pricing (Part 2)
Stability tests can also be performed for numerical methods used for option pricing in finance, such as the binomial model. The binomial model is a discrete-time model that approximates the continuous-time Black-Scholes model for option pricing. The binomial model assumes that the underlying asset can only move up or down in price by a fixed amount at each time step, and the option price is calculated recursively from the final time step back to the present time.
To perform a stability test for the binomial model, we can simulate the model for different time step sizes and compare the results with the true option price or other reference values. The test can be done using synthetic data, or real data can be used to test the model’s performance in real-world scenarios. Here is an example of how you can implement a stability test for the binomial model in Python using synthetic data:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Define the binomial option pricing function
def binomial_option(S, K, r, sigma, T, N, option_type):
dt = T / N
u = np.exp(sigma * np.sqrt(dt))
d = 1 / u
p = (np.exp(r * dt) - d) / (u - d)
q = 1 - p
S_values = np.zeros((N+1, N+1))
V_values = np.zeros((N+1, N+1))
S_values[0, 0] = S
for i in range(1, N+1):
S_values[i, 0] = S_values[i-1, 0] * d
for j in range(1, i+1):
S_values[i, j] = S_values[i-1, j-1] * u
if option_type == 'call':
for j in range(N+1):
V_values[N, j] = max(S_values[N, j] - K, 0)
elif option_type == 'put':
for j in range(N+1):
V_values[N, j] = max(K - S_values[N, j], 0)
for i in range(N-1, -1, -1):
for j in range(i+1):
V_values[i, j] = np.exp(-r * dt) * (p * V_values[i+1, j] + q * V_values[i+1, j+1])
return V_values[0, 0]
# Generate synthetic data
S0 = 100
K = 100
r = 0.05
sigma = 0.2
T = 1
N_list = [10, 20, 50, 100, 200]
dt_list = [T / N for N in N_list]
option_type = 'call'
true_price = 13.28
option_prices = {}
for N, dt in zip(N_list, dt_list):
price = binomial_option(S0, K, r, sigma, T, N, option_type)
option_prices[N] = price
# Plot the results
fig, ax = plt.subplots()
ax.plot(N_list, [true_price] * len(N_list), label='True Price')
ax.plot(N_list, list(option_prices.values()), label='Binomial Model')
ax.set_xlabel('Number of Time Steps (N)')
ax.set_ylabel('Option Price')
ax.legend()
plt.show()
The plot shows the option price calculated using the binomial model for different values of N, as well as the true option price. We can see that as the number of time steps N increases, the option price calculated using the binomial model converges to the true price, indicating that the model is stable.
We can also perform other stability tests, such as varying the other input parameters, such as S0, K, r, and sigma, or using different option types. By comparing the model's results with the true prices or other reference values, we can evaluate the model's stability and accuracy and make any necessary adjustments to improve its performance.
To compare the stability, accuracy, and speed of convergence of the explicit, implicit, and Crank-Nicolson methods in option pricing, we can perform a simulation study. In this simulation, we will generate synthetic data for different input parameters and use each method to calculate the option price. We will then compare the results in terms of stability, accuracy, and speed of convergence.
Let’s start by defining the functions to calculate the option price using the explicit, implicit, and Crank-Nicolson methods. We can use the numpy library to create arrays to store the option price and the coefficients of the differential equation.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
def explicit_fin_diff(S, K, T, sigma, r, q, N, Nj, CallPut):
'''
Explicit finite difference method for pricing European call and put options
Args:
S - intial price of underlying asset
K - strike price
T - time to maturity
sigma - volatility
r - risk-free rate
q - dividend rate
N - number of spacing points along the time partition (horizontal)
Nj - number of partition points from 0 to the upper/lower boundary (vertical)
CallPut - 'Call' or 'Put'
Returns the option price estimated by the finite difference grid
On the choice of N and Nj:
The following condition must be met to guarantee convergence of the explicit finite difference method:
dx >= sigma*sqrt(3*dt)
In fact, the best choice for dx is:
dx = sigma*sqrt(3*dt)
Therefore, finding the number N of time intervals is quite simple. For a given error epsilon, we set:
(sigma^2)*3*dt + dt = epsilon or dt = epsilon/(1 + 3*sigma^2)
Since dt = T/N, we can easily solve for N.
'''
dt = T/N
dx = sigma*np.sqrt(3*dt)
nu = r - q - 0.5*sigma**2
pu = 0.5*dt*((sigma/dx)**2 + nu/dx)
pm = 1.0 - dt*(sigma/dx)**2 - r*dt
pd = 0.5*dt*((sigma/dx)**2 - nu/dx)
grid = np.zeros((N+1, 2*Nj+1))
# Asset prices at maturity:
St = [S*np.exp(-Nj*dx)]
for j in range(1, 2*Nj+1):
St.append(St[j-1]*np.exp(dx))
# Option value at maturity:
for j in range(2*Nj+1):
if CallPut == 'Call':
grid[N, j] = max(0, St[j] - K)
elif CallPut == 'Put':
grid[N, j] = max(0, K - St[j])
# Backwards computing through grid:
for i in range(N-1, -1, -1):
for j in range(1, 2*Nj):
grid[i, j] = pu*grid[i+1, j+1] + pm * \
grid[i+1, j] + pd*grid[i+1, j-1]
# Boundary conditions
grid[i, 0] = grid[i, 1]
grid[i, 2*Nj] = grid[i, 2*Nj-1] + (St[2*Nj]-St[2*Nj-1])
return grid[0, Nj]
explicit_price = explicit_fin_diff(
S=100, K=100, T=1, sigma=0.2, r=0.05, q=0, N=10, Nj=20, CallPut='Put')
print(explicit_price)
def implicit_fin_diff(S, K, T, sigma, r, q, N, Nj, CallPut):
'''
Implicit finite difference method for pricing European call and put options
Args:
S - intial price of underlying asset
K - strike price
T - time to maturity
sigma - volatility
r - risk-free rate
q - dividend rate
N - number of time intervals (horizontal)
Nj - number of partition points from 0 to the upper/lower boundary (vertical)
CallPut - 'Call' or 'Put'
Returns the option price estimated by the finite difference grid
On the choice of N and Nj:
Similar to the case of the explicit method, we again choose dt and dx such that:
(dx)^2 + dt = epsilon
We find that setting each term to 0.5*epsilon and solving for N and Nj provides pretty good results.
'''
dt = T/N # dx = sigma*np.sqrt(3*dt)
dx = 1.0/(2*Nj+1)
nu = r - q - 0.5*sigma**2
pu = -0.5*dt*((sigma/dx)**2 + nu/dx)
pm = 1.0 + dt*(sigma/dx)**2 + r*dt
pd = -0.5*dt*((sigma/dx)**2 - nu/dx)
grid = np.zeros(2*Nj+1)
# Asset prices at maturity:
St = [S*np.exp(-Nj*dx)]
for j in range(1, 2*Nj+1):
St.append(St[j-1]*np.exp(dx))
# Option value at maturity:
for j in range(2*Nj+1):
if CallPut == 'Call':
grid[j] = max(0, St[j] - K)
elif CallPut == 'Put':
grid[j] = max(0, K - St[j])
# Boundary Conditions:
if CallPut == 'Call':
lambdaU = St[2*Nj] - St[2*Nj-1]
lambdaL = 0.0
elif CallPut == 'Put':
lambdaU = 0.0
lambdaL = -1.0*(St[1] - St[0])
# Backwards computing through grid
def tridiagonal(C, pU, pM, pD, lambda_L, lambda_U, nj):
'''
Helper function for solving the tridiagonal matrix system specified by the
implicit finite difference method
'''
C1 = np.zeros(2*nj+1)
pmp = [pM+pD]
pp = [C[1]+pD*lambda_L]
for j in range(2, 2*nj):
pmp.append(pM - pU*pD/pmp[j-2])
pp.append(C[j] - pp[j-2]*pD/pmp[j-2])
C1[2*nj] = (pp[len(pp)-1] + pmp[len(pmp)-1]
* lambda_U)/(pU + pmp[len(pmp)-1])
C1[2*nj-1] = C1[2*nj] - lambda_U
for j in range(2*nj-2, -1, -1):
C1[j] = (pp[j-1] - pU*C1[j+1])/pmp[j-1]
C1[0] = C1[1] - lambda_L
return C1
for i in range(N):
grid = tridiagonal(grid, pu, pm, pd, lambdaL, lambdaU, Nj)
return grid[Nj]
implicit_price = implicit_fin_diff(
S=100, K=100, T=1, sigma=0.2, r=0.05, q=0, N=10, Nj=20, CallPut='Put')
print(implicit_price)
def crank_nicolson(S, K, T, sigma, r, q, N, Nj, CallPut):
'''
Crank-Nicolson finite difference method for pricing European calls and puts
Args:
S - intial price of underlying asset
K - strike price
T - time to maturity
sigma - volatility
r - risk-free rate
q - dividend rate
N - number of spacing points along the time partition (horizontal)
Nj - number of partition points from 0 to the upper/lower boundary (vertical)
CallPut - 'Call' or 'Put'
Returns the option price estimated by the finite difference grid
On the choice of N and Nj:
dt and dx should be chosen specifically such that the following condition holds:
(dx)^2 + (0.5*dt)^2 = epsilon
Again, we find that setting each term to 0.5*epsilon and solving for N and Nj gives pretty good results
'''
dt = T/N
dx = 1.0/(2*Nj+1)
nu = r - q - 0.5*sigma**2
pu = -0.25*dt*((sigma/dx)**2 + nu/dx)
pm = 1.0 + 0.5*dt*((sigma/dx)**2) + 0.5*r*dt
pd = -0.25*dt*((sigma/dx)**2 - nu/dx)
grid = np.zeros(2*Nj+1)
# Asset prices at maturity:
St = [S*np.exp(-Nj*dx)]
for j in range(1, 2*Nj+1):
St.append(St[j-1]*np.exp(dx))
# Option value at maturity:
for j in range(2*Nj+1):
if CallPut == 'Call':
grid[j] = max(0, St[j] - K)
elif CallPut == 'Put':
grid[j] = max(0, K - St[j])
# Boundary Conditions:
if CallPut == 'Call':
lambdaU = St[2*Nj] - St[2*Nj-1]
lambdaL = 0.0
elif CallPut == 'Put':
lambdaU = 0.0
lambdaL = -1.0*(St[1] - St[0])
# Backwards computing through grid:
def tridiagonal(C, pU, pM, pD, lambda_L, lambda_U, nj):
'''
Helper function for solving the tridiagonal matrix system specified by the
Crank-Nicolson finite difference method
'''
C1 = np.zeros(2*nj+1)
pmp = [pM+pD]
pp = [-pU*C[2]-(pM-2)*C[1]-pD*C[0]+pD*lambda_L]
for j in range(2, 2*nj):
pmp.append(pM - pU*pD/pmp[j-2])
pp.append(-pU*C[j+1] - (pM-2)*C[j] - pD *
C[j-1] - pp[j-2]*pD/pmp[j-2])
# Boundary conditions:
C1[2*nj] = (pp[len(pp)-1] + pmp[len(pmp)-1]
* lambda_U)/(pU + pmp[len(pmp)-1])
C1[2*nj-1] = C1[2*nj] - lambda_U
# Back substitution
for j in range(2*nj-2, 0, -1):
C1[j] = (pp[j-1] - pU*C1[j+1])/pmp[j-1]
C1[0] = C[0]
return C1
for i in range(N):
grid = tridiagonal(grid, pu, pm, pd, lambdaL, lambdaU, Nj)
return grid[Nj]
cn_price = crank_nicolson(S=100, K=100, T=1, sigma=0.2,
r=0.05, q=0, N=10, Nj=20, CallPut='Put')
print(cn_price)
# Define parameters
S = 100
K = 100
r = 0.05
sigma = 0.2
T = 1
q = 0
Nj = 100
# Define time step sizes to test
N_values = [10, 20, 40, 80, 160, 320, 640]
# Define option type
CallPut = 'Call'
# Initialize arrays to store option prices
explicit_prices = np.zeros(len(N_values))
implicit_prices = np.zeros(len(N_values))
cn_prices = np.zeros(len(N_values))
# print(explicit_prices.shape, implicit_prices.shape, cn_prices.shape)
# Calculate option prices for different time step sizes using each method
for i, N in enumerate(N_values):
dt = T / N
explicit_prices[i] = explicit_fin_diff(
S, K, T, sigma, r, q, N, Nj, CallPut)
implicit_prices[i] = implicit_fin_diff(
S, K, T, sigma, r, q, N, Nj, CallPut)
cn_prices[i] = crank_nicolson(S, K, T, sigma, r, q, N, Nj, CallPut)
# Plot option prices as a function of time step size
plt.plot(N_values, explicit_prices, label='Explicit')
plt.plot(N_values, implicit_prices, label='Implicit')
plt.plot(N_values, cn_prices, label='Crank-Nicolson')
plt.xlabel('Time step size')
plt.ylabel('Option price')
plt.legend()
plt.show()
The resulting plot shows the option prices as a function of the time. Now that we have the implementations of the three methods, we can perform a comparison study of their stability, accuracy, and speed of convergence in option pricing. For this case, we will only consider ‘Call’ options. First, we will investigate the stability of the methods by varying the size of the time step and observing the behavior of the option prices. We will fix the other parameters as follows: S = 10, K = 10, r = 0.05, sigma = 0.2, T = 1, M = 100.
# Define parameters
S = 10
K = 10
r = 0.05
sigma = 0.2
Smax = 30
T = 1
M = 100
# Define time step sizes to test
N = [10, 20, 40, 80, 160, 320, 640]
# Define option type
CallPut = 'Call'
s = np.arange(10,30,0.5)
# Initialize arrays to store option prices
explicit_prices = np.zeros(len(N))
implicit_prices = np.zeros(len(N))
cn_prices = np.zeros(len(N))
# Calculate option prices for different time step sizes using each method
for i, N in enumerate(N):
dt = T / N
explicit_prices[i] = explicit_fin_diff(S, K, T, sigma, r, q, N, Nj, CallPut)
implicit_prices[i] = implicit_fin_diff(S, K, T, sigma, r, q, N, Nj, CallPut)
cn_prices[i] = crank_nicolson(S, K, T, sigma, r, q, N, Nj, CallPut)
print(explicit_prices.shape, implicit_prices.shape, cn_prices.shape)
print(explicit_prices)
# # Plot option prices as a function of time step size
plt.plot(explicit_prices, label='Explicit')
plt.plot(implicit_prices, label='Implicit')
plt.plot(cn_prices, label='Crank-Nicolson')
plt.xlabel('Time step size')
plt.ylabel('Option price')
plt.legend()
plt.show()
# Define option prices from Black-Scholes formula
d1 = (np.log(S0 / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
bs_price = S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
# Calculate option prices using each method
explicit_price = explicit_fin_diff(S, K, T, sigma, r, q, N, Nj, CallPut)
implicit_price = implicit_fin_diff(S, K, T, sigma, r, q, N, Nj, CallPut)
cn_price = crank_nicolson(S, K, T, sigma, r, q, N, Nj, CallPut)
# Print option prices
print('Black-Scholes price: {:.4f}'.format(bs_price))
print('Explicit price: {:.4f}'.format(explicit_price))
print('Implicit price: {:.4f}'.format(implicit_price))
print('Crank-Nicolson price: {:.4f}'.format(cn_price))
# Define grid sizes to test
M = [10, 20, 40, 80, 160, 320, 640]
# Initialize arrays to store option prices
explicit_prices = np.zeros(len(M))
implicit_prices = np.zeros(len(M))
cn_prices = np.zeros(len(M))
# Calculate option prices for different grid sizes using each method
for i, M in enumerate(M):
explicit_prices[i] = explicit_fin_diff(S, K, T, sigma, r, q, N, Nj, CallPut)
implicit_prices[i] = implicit_fin_diff(S, K, T, sigma, r, q, N, Nj, CallPut)
cn_prices[i] = crank_nicolson(S, K, T, sigma, r, q, N, Nj, CallPut)
# Calculate error relative to Black-Scholes price
explicit_error = np.abs(explicit_prices - bs_price) / bs_price
implicit_error = np.abs(implicit_prices - bs_price) / bs_price
cn_error = np.abs(cn_prices - bs_price) / bs_price
# Plot error as a function of spatial grid size
plt.loglog(explicit_error, label='Explicit')
plt.loglog(implicit_error, label='Implicit')
plt.loglog(cn_error, label='Crank-Nicolson')
plt.xlabel('M (grid size)')
plt.ylabel('Relative error')
plt.legend()
plt.show()
We can see that the explicit method becomes unstable as the time step size increases, while the implicit and Crank-Nicolson methods remain stable for all time step sizes. This is in line with our expectation, as the explicit method has a condition for stability that restricts the size of the time step. Therefore, for large values of N, the explicit method becomes unstable and produces nonsensical results. In contrast, the implicit and Crank-Nicolson methods do not have such a restriction on the time step size and are therefore more stable.
Next, we will investigate the accuracy and speed of convergence of the methods. We will do this by comparing the option prices obtained from the methods with the exact solution obtained from the Black-Scholes formula.
We can see that the explicit method produces the least accurate option price, while the implicit and Crank-Nicolson methods produce more accurate prices. However, the implicit method is slightly more accurate than the Crank-Nicolson method.
Finally, we will investigate the speed of convergence of the methods by varying the size of the spatial grid and observing the behavior of the option prices. We will fix the other parameters as follows: S = 100, K = 100, r = 0.05, sigma = 0.2, T = 1, N = 100.
The resulting plot shows the relative error of each method as a function of the spatial grid size. We can see that the error decreases for all methods as the grid size increases, indicating that the methods are converging to the exact solution. The Crank-Nicolson achieves the highest speed of convergence as seen by its line lower than that of explicit and implicit for grid size less than 4x10⁰. However, the Crank-Nicolson method diverges and we see an increment in the relative error for a grid size greater than 4x10⁰.
In conclusion, we have shown how to implement stability tests for numerical methods and how to compare the stability, accuracy, and speed of convergence of different numerical methods for option pricing. In our example, we found that the explicit method is unstable for large time step sizes, while the implicit and Crank-Nicolson methods are stable for all time step sizes. We also found that the Crank-Nicolson method is the most accurate and efficient method in terms of speed of convergence except for grid sizes greater that 4x10⁰. These results demonstrate the importance of selecting appropriate numerical methods for option pricing and highlight the trade-offs between stability, accuracy, and speed of convergence.
P.S. There’ll be a part 3 for the same series which will exclusively cover American Options. Thanks again!
