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 = S0*2 / M
s = np.arange(S0*2, 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<= (ds**2/(2*(sigma*S0)**2))):
  print('Courant-Friedrichs-Lewy (CFL) condition is verified')
  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-K*np.exp(-r*T[0,:])

# Explicit scheme
for i in range(N-2,-1,-1): 
    V[1:-1,i] =  V[2:,i+1] * (0.5*dt*j*( r + sigma**2*j )) + V[1:-1,i+1]*(1- dt*( (sigma**2)*j**2 + r)) + V[0:-2,i+1]*(0.5*dt*j*( -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_zlabel('Option Value')
ax.set_title('Black-Scholes Option Value')

If someone is able to spot the mistake it would be a great help! Thanks a lot :)

  • $\begingroup$ when i tried running your code, it complained about t array being true or false. maybe try set T = 1.0 (float) instead of 1. Not sure if this is the cause. $\endgroup$
    – mirmo
    Mar 27, 2023 at 5:18
  • $\begingroup$ Solving the PDE in its original variables will inevitably cause problems for not even that large volatilities. It is better to apply the standard transformation first. For a detailed post showing the problems see here. $\endgroup$
    – Kurt G.
    Mar 27, 2023 at 10:55

1 Answer 1


Edit: based on comment.

In your question, you didn't specify the boundary conditions being considered.

V[0,:] = S0-K*np.exp(-r*T[0,:])

As this condition from the Put-Call parity is to approximate large S, it should be S here instead of the initial price, S0. In this case, it should be

V[0,:] = S0*2-K*np.exp(-r*T[0,:])

as per your maximum defined S.

  • $\begingroup$ Indeed, when S is large compared to the strike price, the call option price is : $$ C = S - K e^{-r\tau} $$ This can be understood using the call put parity. We see that if S is very large compared to K, the put price is zero, therefore: $$ C = C - P = S - K e^{-r \tau} $$ $\endgroup$ Apr 10, 2023 at 14:00
  • $\begingroup$ Yes, but you are using $S_0$ instead of $S$? $\endgroup$
    – mirmo
    Apr 12, 2023 at 2:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.