Correcting the High-Stock-Price Boundary in an Explicit Black–Scholes Scheme

# Black Scholes PDE explicit Scheme


# Black Scholes PDE explicit Scheme












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')
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-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_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 :)

## Answer by mirmo (score 0)

https://quant.stackexchange.com/a/75068

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.