Finite Difference Methods and the Black-Scholes-Merton Equation: A Numerical Approach to Option Pricing

Introduction: Numerical Methods in Financial Engineering

In financial engineering, the valuation of options plays a central role in decision-making and risk management. A well-known model for pricing options is the Black-Scholes-Merton (BSM) equation, a partial differential equation (PDE) that provides a theoretical estimate for the price of options based on the underlying asset’s dynamics. Although the BSM equation can be solved analytically in some cases (e.g., for European options), most real-world scenarios, particularly American options or scenarios with dividends, require numerical methods to approximate solutions.

One popular numerical approach is Finite Difference Methods (FDM). These methods convert continuous differential equations into systems of algebraic equations that can be solved iteratively. By discretizing time and asset prices, FDM enables us to numerically solve the BSM equation in cases where analytical solutions are not feasible.

This article will explore key FDM techniques—explicit and implicit schemes—and their application to the Black-Scholes-Merton equation. We’ll also discuss Von Neumann stability analysis, a crucial tool for understanding the stability of these numerical methods.

Finite Difference Methods: Key Concepts

Approximating Derivatives with Finite Differences

The basic idea of FDM is to approximate the derivatives in the BSM equation using algebraic expressions. This transforms the continuous PDE into a discrete system, which can be solved numerically.

1. Forward Difference Approximation

The forward difference approximation estimates the first derivative of a function using the values at the current and next points. Given a small step size $h$, the first derivative of $f(x)$ at point $x$ is approximated by:

\[f'(x) \approx \frac{f(x + h) - f(x)}{h}\]

This approximation is first-order accurate, meaning the error scales linearly with $h$.

2. Backward Difference Approximation

Similarly, the backward difference approximation uses the current and previous points to estimate the first derivative:

\[f'(x) \approx \frac{f(x) - f(x - h)}{h}\]

Like the forward difference, this method also has first-order accuracy.

3. Central Difference Approximation

A more accurate estimate of the first derivative comes from the central difference approximation, which averages the forward and backward differences:

\[f'(x) \approx \frac{f(x + h) - f(x - h)}{2h}\]

This method is second-order accurate, providing a better balance between approximation error and computation.

4. Second Derivative Approximation

For the second derivative, we add the Taylor expansions of $f(x + h)$ and $f(x - h)$:

\[f''(x) \approx \frac{f(x + h) - 2f(x) + f(x - h)}{h^2}\]

This approximation is also second-order accurate, essential for solving the BSM equation, which involves second-order spatial derivatives.

The Black-Scholes-Merton Differential Equation

The Black-Scholes-Merton (BSM) equation models the price $V(S, t)$ of an option as a function of the asset price $S$ and time $t$. The equation, derived from stochastic processes and the assumption of no arbitrage, is given by:

\[\frac{\partial V}{\partial t} + \frac{1}{2} \sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS \frac{\partial V}{\partial S} - rV = 0\]

where:

$V(S, t)$ is the option price as a function of asset price $S$ and time $t$,
$\sigma$ is the volatility of the asset,
$r$ is the risk-free interest rate.

The BSM equation is typically solved with boundary conditions, depending on whether we are pricing a call option or a put option.

Assumptions of the BSM Model

The BSM model operates under several key assumptions:

Lognormal distribution of asset prices: Asset prices follow geometric Brownian motion.
No arbitrage: It assumes no arbitrage opportunities exist in the market.
Constant volatility and risk-free rate: Both volatility and the risk-free interest rate remain constant throughout the life of the option.
Frictionless markets: No transaction costs, taxes, or liquidity constraints.
European-style options: The model primarily applies to European options, which can only be exercised at expiration.

Discretizing the Black-Scholes-Merton Equation

To solve the BSM equation numerically, we discretize the time and asset price domains. Let $\Delta t$ be the time step, and $\Delta S$ be the asset price step. We define grid points as follows:

$S_i = i \Delta S$ for $i = 0, 1, \dots, N$,
$t_n = n \Delta t$ for $n = 0, 1, \dots, M$.

Boundary Conditions

For a call option, the boundary conditions are:

At expiry: $V(S, T) = \max(S - K, 0)$,
At the lower boundary: $V(0, t) = 0$,
At the upper boundary: $V(S_{\text{max}}, t) \approx S_{\text{max}} - K$.

For a put option, the boundary conditions are:

At expiry: $V(S, T) = \max(K - S, 0)$,
At the lower boundary: $V(0, t) = K$,
At the upper boundary: $V(S_{\text{max}}, t) = 0$.

Numerical Schemes for Solving the BSM Equation

Explicit Scheme

The explicit scheme is the simplest approach to solving the BSM equation. It uses known values at the current time step to compute values at the next time step.

Steps in the Explicit Scheme:

Approximate the time derivative using a forward difference.
Approximate the spatial derivatives using central differences.
Use the finite difference approximations to update the option price at each grid point.

The update formula for the option price $V_i^n$ at the next time step is:

\[V_i^{n+1} = V_i^n + \Delta t \left( \frac{1}{2} \sigma^2 S_i^2 \frac{V_{i+1}^n - 2V_i^n + V_{i-1}^n}{(\Delta S)^2} + rS_i \frac{V_{i+1}^n - V_{i-1}^n}{2\Delta S} - rV_i^n \right)\]

While straightforward, the explicit scheme is conditionally stable, meaning that the time step $\Delta t$ must be sufficiently small to ensure accuracy and avoid divergence.

Implicit Scheme

The implicit scheme is more stable than the explicit scheme, allowing for larger time steps. Unlike the explicit method, the implicit scheme computes option prices by solving a system of linear equations at each time step.

Steps in the Implicit Scheme:

Approximate the time derivative using a backward difference.
Approximate spatial derivatives at the future time step using central differences.
Solve the resulting system of linear equations (often tridiagonal) to update the option prices.

The implicit method is unconditionally stable, meaning it can handle larger time steps without the risk of instability, but it requires more computation due to the need to solve the linear system.

Von Neumann Stability Analysis

Stability is a critical concern in numerical methods. The Von Neumann stability analysis helps determine whether a numerical scheme will produce stable solutions. For the BSM equation, we analyze the growth or decay of errors introduced by discretization.

The key idea is to express the error as a Fourier series and analyze whether the amplification factor $\xi$ satisfies $$

\xi

\leq 1$for all Fourier modes. For the explicit scheme, stability depends on the size of$\Delta t$$, while the implicit scheme remains stable regardless of the time step.

Choosing the Right Scheme

Both explicit and implicit schemes are useful for numerically solving the Black-Scholes-Merton equation, but they each have their strengths and weaknesses. The explicit scheme is simple and computationally efficient for small time steps but becomes unstable for larger steps. The implicit scheme, while more computationally intensive, is unconditionally stable and better suited for larger time steps or more complex option pricing scenarios.

For practitioners in financial engineering, understanding these numerical methods and their stability properties is crucial for accurately pricing options and managing risk in volatile markets.

Appendix: Python Code for Finite Difference Methods in Option Pricing

This appendix provides Python implementations for solving the Black-Scholes-Merton (BSM) equation using Finite Difference Methods (FDM). The examples include both the explicit scheme and implicit scheme for pricing a European call option.

Required Libraries

Make sure you have installed the necessary libraries before running the code. You can install them using pip:

pip install numpy scipy matplotlib

1. Explicit Scheme for the Black-Scholes Equation

This implementation uses the explicit finite difference method to price a European call option.

import numpy as np
import matplotlib.pyplot as plt

# Define option parameters
S_max = 100  # Maximum stock price
T = 1.0      # Time to maturity (1 year)
K = 50       # Strike price
r = 0.05     # Risk-free interest rate
sigma = 0.2  # Volatility

# Discretization parameters
M = 100  # Number of time steps
N = 100  # Number of price steps

# Grid step sizes
dt = T / M
dS = S_max / N

# Create the grid
S = np.linspace(0, S_max, N + 1)
V = np.zeros((M + 1, N + 1))

# Boundary conditions
V[-1, :] = np.maximum(S - K, 0)  # Payoff at maturity
V[:, 0] = 0  # V(S=0, t) = 0
V[:, -1] = S_max - K * np.exp(-r * (T - np.linspace(0, T, M + 1)))  # V(S=S_max, t)

# Coefficients for the explicit scheme
alpha = 0.5 * dt * (sigma**2 * S**2 / dS**2 - r * S / dS)
beta = -dt * (sigma**2 * S**2 / dS**2 + r)
gamma = 0.5 * dt * (sigma**2 * S**2 / dS**2 + r * S / dS)

# Explicit finite difference method
for m in range(M, 0, -1):
    for n in range(1, N):
        V[m-1, n] = alpha[n] * V[m, n-1] + beta[n] * V[m, n] + gamma[n] * V[m, n+1]

# Plot the result
plt.plot(S, V[0, :])
plt.title('Option Price using Explicit Scheme (European Call)')
plt.xlabel('Stock Price (S)')
plt.ylabel('Option Price (V)')
plt.grid(True)
plt.show()

2. Implicit Scheme for the Black-Scholes Equation

Next, we implement the implicit scheme, which requires solving a tridiagonal system of equations at each time step.

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

# Define option parameters (same as above)
S_max = 100
T = 1.0
K = 50
r = 0.05
sigma = 0.2

# Discretization parameters
M = 100  # Number of time steps
N = 100  # Number of price steps

# Grid step sizes
dt = T / M
dS = S_max / N

# Create the grid
S = np.linspace(0, S_max, N + 1)
V = np.zeros((M + 1, N + 1))

# Boundary conditions
V[-1, :] = np.maximum(S - K, 0)  # Payoff at maturity
V[:, 0] = 0
V[:, -1] = S_max - K * np.exp(-r * (T - np.linspace(0, T, M + 1)))

# Coefficients for the implicit scheme
alpha = 0.5 * dt * (sigma**2 * S**2 / dS**2 - r * S / dS)
beta = -dt * (sigma**2 * S**2 / dS**2 + r)
gamma = 0.5 * dt * (sigma**2 * S**2 / dS**2 + r * S / dS)

# Build the tridiagonal matrix for the implicit scheme
A = np.zeros((N-1, N-1))

for i in range(1, N):
    if i == 1:
        A[i-1, i-1] = 1 - beta[i]
        A[i-1, i] = -gamma[i]
    elif i == N-1:
        A[i-1, i-2] = -alpha[i]
        A[i-1, i-1] = 1 - beta[i]
    else:
        A[i-1, i-2] = -alpha[i]
        A[i-1, i-1] = 1 - beta[i]
        A[i-1, i] = -gamma[i]

# Solve using backward time stepping
for m in range(M, 0, -1):
    B = V[m, 1:N]
    V[m-1, 1:N] = la.solve_banded((1, 1), A, B)

# Plot the result
plt.plot(S, V[0, :])
plt.title('Option Price using Implicit Scheme (European Call)')
plt.xlabel('Stock Price (S)')
plt.ylabel('Option Price (V)')
plt.grid(True)
plt.show()

3. Stability Analysis

To ensure the explicit scheme is stable, we can perform a quick analysis by checking the time step condition $\Delta t \leq \frac{\Delta S^2}{2 \sigma^2 S_{\text{max}}^2}$. Here’s a snippet for checking the stability condition:

def stability_check(sigma, S_max, dS, dt):
    stability_threshold = dS**2 / (2 * sigma**2 * S_max**2)
    if dt <= stability_threshold:
        print(f"Stable: Time step {dt:.5f} <= {stability_threshold:.5f}")
    else:
        print(f"Unstable: Time step {dt:.5f} > {stability_threshold:.5f}")

# Run the stability check
stability_check(sigma, S_max, dS, dt)

Summary of Code

Explicit Scheme: This method computes the next time step based on known values at the current time step. It is simple but conditionally stable, so care must be taken to ensure that the time step is small enough.
Implicit Scheme: More robust and unconditionally stable, the implicit scheme requires solving a system of linear equations but allows for larger time steps.

Both implementations are useful for pricing European options under the Black-Scholes-Merton framework and can be extended to more complex scenarios, such as American options, by incorporating early exercise conditions.

By running the code provided, you can numerically solve the Black-Scholes equation and visualize the option prices for different asset prices using both finite difference methods.

Check out Data Science Books on Amazon

Share on

Twitter Facebook LinkedIn