I am currently working on the implementation of classic schemes to solve the BS PDE and it seems that I make a mistake in my code because the result looks far from the result of the BS formula.
Here is the maths of the implicit scheme:
Substituting $\tau = T - t$ into the Black-Scholes equation and using the chain rule, we get:
$$-\frac{\partial V}{\partial \tau} + \frac{1}{2}\sigma^2S^2\frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0$$
To solve this equation using an explicit scheme, we first need to discretize it in both time and space. Let $\Delta t$ be the time step and $\Delta S$ be the space step. We define:
$$\tau_n = n\Delta \tau \qquad\text{and}\qquad S_j = j\Delta S$$
for integers $n$ and $j$. We also define the numerical solution at the grid points $(S_j, \tau_n)$ as $V_j^n \approx V(S_j, \tau_n)$.
We approximate the derivative with finite differences as follows:
$$\frac{\partial V}{\partial \tau}(S_j,\tau_n) \approx \frac{V_j^{n+1} - V_j^n}{\Delta \tau}$$
$$\frac{\partial V}{\partial S}(S_j,\tau_n) \approx \frac{V_{j+1}^n - V_{j-1}^n}{2\Delta S}$$
$$\frac{\partial^2 V}{\partial S^2}(S_j,\tau_n) \approx \frac{V_{j+1}^n - 2V_j^n + V_{j-1}^n}{\Delta S^2}$$
Substituting these approximations into the Black-Scholes equation, we get:
$$-\frac{V_j^{n+1} - V_j^n}{\Delta \tau} + \frac{1}{2}\sigma^2 S_j^2 \frac{V_{j+1}^n - 2V_j^n + V_{j-1}^n}{\Delta S^2} + rS_j\frac{V_{j+1}^n - V_{j-1}^n}{2\Delta S} - rV_j^n = 0$$
Rearranging this equation, we obtain an explicit scheme for $V_j^{n+1}$:
$$V_j^{n+1} = \frac{1}{2}\sigma^2 j^2\Delta S^2 \Delta \tau \frac{V_{j+1}^n - 2V_j^n + V_{j-1}^n}{\Delta S^2} + rj\Delta S\Delta \tau\frac{V_{j+1}^n - V_{j-1}^n}{2\Delta S} + (1-r\Delta \tau)V_j^n$$
So : $$V_j^{n+1} = V_{j+1}^n \left(\frac{\Delta \tau j}{2} \left(r+ \sigma^2 j\right)\right) + V_{j}^n \left(1 - \Delta \tau \left(\sigma^2 j^2 + r\right)\right) + V_{j-1}^n \left(\frac{\Delta \tau j}{2} \left(-r + \sigma^2 j\right)\right)$$
Here is the code :
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
Parameters
S0 = 100 # Initial stock price
K = 100 # Strike price
T = 1 # Time to maturity
r = 0.05 # Risk-free rate
sigma = 0.2 # Volatility
N = 2000 # Number of time steps
M = 100 # Number of stock price steps
Grid setup
dt = T / N
ds = S02 / M
s = np.arange(S02, 0, -ds)
t = np.arange(T, 0,-dt)
T, S = np.meshgrid(t, s)
j = np.arange(M-1,1,-1)
#Courant-Friedrichs-Lewy (CFL) condition,
if(dt<= (ds2/(2(sigmaS0)2))):
print('Courant-Friedrichs-Lewy (CFL) condition is verified')
else:
print('Courant-Friedrichs-Lewy (CFL) condition IS NOT verified')
Initial conditions
V = np.maximum(S-K, 0)
Boundary conditions
V[-1,:] = 0
V[0,:] = S0-Knp.exp(-rT[0,:])
Explicit scheme
for i in range(N-2,-1,-1):
V[1:-1,i] = V[2:,i+1] * (0.5dtj( r + sigma2j )) + V[1:-1,i+1](1- dt( (sigma2)*j2 + r)) + V[0:-2,i+1](0.5dtj( -r + (sigma*2)j ))
#Plot the result
fig = plt.figure(figsize=(20,12))
ax = plt.axes(projection='3d')
ax.plot_surface(S, T, V, cmap='coolwarm')
ax.view_init(10, 250)
ax.set_xlabel('Stock Price')
ax.set_ylabel('Time')
ax.set_zlabel('Option Value')
ax.set_title('Black-Scholes Option Value')
plt.show()
If someone is able to spot the mistake it would be a great help! Thanks a lot :)
