Explicit Finite Difference with Black-Scholes Formula (with code)

This is my first article from the series about options pricing and volatility modeling. I assume a foundational understanding in certain instances to move toward more practical examples of implementing theory into practice. I’ll try to mix the articles in terms of difficulty; I’ll cover volatility models, different options types (vanilla, exotics like multi-asset, Parisian, and other crazy ones), and some less popular topics (e.g. implied probabilities or pricing options under fat tails). Okay, enough with the introduction…

In this article, I’m diving into applying the Black-Scholes formula using the Finite Difference Method, taking cues from “Paul Wilmott on Quant Finance”. My goal is to walk you through this idea with code, emphasizing its importance and practical uses (and without fancy-schmancy terms). It is not supposed to be sophisticated, but rather the first step in the topic. I’ll expand on more advanced techniques in the future.

The method we’re employing here is referred to as “explicit” due to its direct and step-by-step approach to solving the involved equations. This directness makes it relatively easier to comprehend and implement in code compared to more complex implicit methods. It’s a great starting point to grasp the intuition behind the process.

However, while it’s easier to code, the explicit method lacks the robustness of the implicit method. It’s crucial to pay close attention to the outcomes, as this approach may lead to erratic results. The added complexity of the implicit method pays off, as it’s more dependable in convergence.

Black Scholes and explicit implementation (European options)

Binomial vs grid (our choice)

Finite Difference is an efficient and very powerful technique, especially on lower dimension problems (number of parameters). It can intuitively be grasped as a binomial tree, but fixed as a grid. Once implemented, it is easy to transform it to the other problems. It is a crucial concept to grasp if you want to implement many concepts in practice.

Intuitively, we are going to move in specific time steps (today, tomorrow, etc.) and inside these time steps in the range of possible asset values (columns). We are going to set up the equation in a convenient form, so we will have to perform a few operations.

Okay, so let’s see the equation. This is so-called Black-Scholes PDE:

Black Scholes PDE; PDE with introduced parameters

Here, we have three “Greeks” that represent the option’s sensitivity to parameters:
Delta: Measures the option price change concerning asset price, i.e., the rate at which the option price changes for a $1 movement in the underlying asset.
Gamma: Represents the rate at which Delta changes, i.e., the second derivative in terms of asset price.
Theta: Quantifies the option price’s time decay rate, showing how much the option’s price changes over time, all else constant (derivative in respect to time).

Above, we have our derivatives written in terms of discrete values. This is the clue of finite difference — changing continuous equations to discrete ones.

Here we have rearranged the equation to have the “k+1” on one side. It will be useful when implementing the equation. It might look scary, but it is easier than it looks (with time, of course).

Implementation

On the one-hand side, we’ll have theta, on the other gamma and delta. We’ll first compute gamma and delta to get theta. The value of the option at a specific asset level (row in the grid) is going to be:
Value today = Value tomorrow — time decay (theta)

First, we need to initialize our parameters inside

S = np.zeros(NAS + 1)  # Asset array
dS = 2 * strike / NAS  # 'Infinity' is twice the strike
dt = 0.9 / vol ** 2 / NAS ** 2  # For stability
NTS = int(expiration / dt) + 1  # Number of time steps
dt = expiration / NTS  # To ensure that expiration is an integer number of time steps away
V = np.zeros((NAS + 1, NTS + 1))  # Option value array
q = 1

The meaning of each parameter is as follows:
- S: This is an asset array. In other words, the number of possible asset values. We have to add 1, because if we have e.g. 1 asset step, then there are 2 states. The addition of 1 ensures the correctness
- dS: This is the upper boundary of the price. This is the “infinity”. In practice, you don’t have to compute the values up to “infinity”, the value 2–3x higher than the strike is satisfactory
- dt: This is our time step. The equation might look strange, but it is chosen to ensure the stability of the solution, which in simple words means: making sure that the result is reliable and realistic
- NTS: Number of time steps made as an integer. To have discrete time steps (i.e. +1, +1)
- dt: It adjusts ‘dt’ to ensure that the expiration time is an integer multiple of the calculated time step

Ok, now that we got to know all the variables and why they matter, let’s initialize our grid and set the initial

V = np.zeros((NAS + 1, NTS + 1))
q = 1
if p_type == "put":
    q = -1  # Test for call or put

for i in range(NAS + 1):
    S[i] = i * dS
    V[i, 0] = max(q * (S[i] - strike), 0)
 
# Vecotorized solution, i.e. better solution
S = np.arange(0, NAS + 1) * dS
payoff_values = np.maximum(q * (S - strike), 0)
V[:, 0] = payoff_values

V: This is our actual grid. It is a matrix with asset values vertically, and time steps horizontally
q: This is important for determining our payoff function. It adjusts the payoff to either call or put.
S[i] = i * dS: It goes one by one of the asset steps and calculates the values. (E.g. 1*ds = 10, 2*ds = 20. It just determines the “start values” of our asset
V[i, 0]: This is, again, going through every asset price in the same column and setting the value of our option. This is just our payoff function of call/put)

for k in range(1, NTS + 1):  # Time loop
    for i in range(1, NAS):  # Asset loop
        delta = (V[i + 1, k - 1] - V[i - 1, k - 1]) / (2*dS)  # Central difference
        gamma = (V[i + 1, k - 1] - 2 * V[i, k - 1] + V[i - 1, k - 1]) / (dS*dS)  # Central difference
        theta = -0.5 * vol ** 2 * S[i] ** 2 * gamma - int_rate * S[i] * delta + int_rate * V[i, k - 1]
        V[i, k] = V[i, k - 1] - dt * theta
    
    V[0, k] = V[0, k - 1] * (1 - int_rate * dt)  # Boundary condition at S=0
    V[NAS, k] = 2 * V[NAS - 1, k] - V[NAS - 2, k]  # Boundary condition at S=infinity

Ok, now we have 2 loops. This is where the magic happens.
First, we determine the time point, which means, we go horizontally. Then in column (time step), we go through the asset values (second loop).

To determine the option value, we need 3 aforementioned Greeks from our diffusion process. They are computed from previous time steps (k-1).
Delta: This is the central difference, which means we take values i+1 and i-1.
Gamma: This is also a central difference, but the second derivative
Theta: This is computed as with our diffusion process. Theta is intuitive— how the option value decreases in the next time step with other parameters held constant.
V[i, k] = V[i, k — 1] — dt * theta: This is our option value at a specific grid point. The value of the option is computed as the option value from the previous time step (today's value given yesterday's value) minus (-) Time decay (theta).

V[0, k] = V[0, k - 1] * (1 - int_rate * dt)  # Boundary condition at S=0
V[NAS, k] = 2 * V[NAS - 1, k] - V[NAS - 2, k]  # Boundary condition at S=infinity

At the end of each time step, we have to apply the boundary conditions. At an asset price of 0 and our infinity (i.e. 2–3 x the Strike). These boundary conditions are crucial to ensure the numerical stability and accuracy of the finite difference method.
V[0, k]: At S = 0, we switch off the diffusion and drift terms. This means that on S0 the payoff is guaranteed, resulting in the condition

V[NAS, k]: As the option payoff for large asset values is linear (deep ITM, delta = 1, acts like underlying), we can use the extrapolation as the boundary condition, which mathematically is written below:

Now that we know how the finite difference works, let’s see the results:)

Results

The Finite Difference Method implementation of the Black-Scholes formula gives us a solid grip on how to value options with the method. It’s not just a technique, it’s the backbone of making sense of the financial derivatives.

sigma = 0.2
r = 0.05
K = 100
T = 1
NAS = 20 # number of asset steps

option_df = option_value_3d(sigma, r, "call", K, T, NAS)

The results:

Heatmap:

3D option surface:

American Options

The change in the case of the American option is straightforward. We just have to compare if the immediate exercise of the option is more profitable than holding it. So, if the option value is below the intrinsic. Don’t mistake it for the market value. If the option with our model is worth 10 and the market value is 1000, then you should not exercise. You just basically sell it to someone else.

The modification:

discounted_payoff = Payoff[i] * np.exp(-int_rate * k* dt)
V[i, k] = np.maximum(V[i, k], discounted_payoff)

I hope you enjoyed this short article. For any questions, feel free to comment, I’m happy to help.

For a code, refer to the Jupyter Notebook provided below:

https://github.com/antek0308/Volatility_notebooks/blob/main/Medium/Finite_diff_explicit.ipynb