Finite Difference Methods for Option Pricing › Boundary Conditions

Michael BrenndoerferDecember 4, 20254 min read
Quantitative FinanceAll chapters
Part 6 of 15Numerical Methods for Derivative Pricing
Chapter 36Finite Difference Methods for Option Pricing › Boundary Conditions
Reading section 6 of 14 · Boundary Conditions · 4 minSection 6 of 14 · 4 min

Book contents

Section 6 of 14 · 4 min read

Finite Difference Methods for Option Pricing › Boundary ConditionsOne section · Sections in this chapter6 / 14 ↓
Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Proper boundary conditions are essential for finite difference methods. They determine the option value at the edges of the computational grid, where the standard update formulas cannot be applied directly because they would require values outside the grid.

Terminal Condition

At expiration (t=Tt = T), the option value equals its payoff. This is the fundamental condition that initializes our backward-in-time computation:

V(S,T)=max⁡(S−K,0)(call)V(S,T)=max⁡(K−S,0)(put)\begin{aligned} V(S, T) &= \max(S - K, 0) && \text{(call)} \\ V(S, T) &= \max(K - S, 0) && \text{(put)} \end{aligned}

where:

  • V(S,T)V(S,T): option value at expiration
  • SS: stock price at expiration
  • KK: strike price

This condition initializes the grid before we begin the backward time-stepping process. The non-smooth nature of the payoff function at S=KS = K can introduce complications, as the second derivative is undefined at that point. You might sometimes smooth the payoff near the strike to improve numerical behavior.

Boundary at S=0S = 0

When the underlying asset price is zero, it remains at zero forever (assuming geometric Brownian motion). This follows because geometric Brownian motion has the property that zero is an absorbing state: if dS=μSdt+σSdWdS = \mu S dt + \sigma S dW, then S=0S = 0 implies dS=0dS = 0. Therefore:

  • Call option:

    V(0,t)=0V(0, t) = 0

    where:

    • V(0,t)V(0, t): option value at asset price zero

    (a call is worthless if the stock is worthless)

  • Put option:

    V(0,t)=Ke−r(T−t)V(0, t) = Ke^{-r(T-t)}

    where:

    • V(0,t)V(0, t): option value at asset price zero
    • KK: strike price
    • rr: risk-free interest rate
    • T−tT-t: time to expiration

    (a put pays KK at expiration with certainty)

Boundary at S=Smax⁡S = S_{\max}

At very high stock prices, we need to specify the option value at the upper edge of the grid. The exact value depends on the type of option:

  • Call option:

    V(Smax⁡,t)≈Smax⁡−Ke−r(T−t)V(S_{\max}, t) \approx S_{\max} - Ke^{-r(T-t)}

    where:

    • V(Smax⁡,t)V(S_{\max}, t): option value at maximum grid price
    • Smax⁡S_{\max}: maximum asset price on the grid
    • KK: strike price
    • rr: risk-free interest rate
    • T−tT-t: time to expiration

    (deep in-the-money, behaves like a forward)

  • Put option:

    V(Smax⁡,t)≈0V(S_{\max}, t) \approx 0

    where:

    • V(Smax⁡,t)V(S_{\max}, t): option value at maximum grid price

    (deep out-of-the-money)

The boundary at Smax⁡S_{\max} requires some care. If Smax⁡S_{\max} is chosen too small, the boundary condition may distort the solution in the region of interest. A rule of thumb is to set Smax⁡≥3KS_{\max} \geq 3K to 5K5K for typical volatilities.

Implementing Put Option Boundaries

In[20]:
Code
def crank_nicolson_european_put(S_max, K, T, r, sigma, M, N):
    """
    Price a European put option using Crank-Nicolson.
    """
    dS = S_max / M
    dt = T / N

    S = np.linspace(0, S_max, M + 1)

    # Terminal condition: put payoff
    V = np.maximum(K - S, 0)

    # Coefficients
    j = np.arange(1, M)
    alpha = 0.25 * dt * (sigma**2 * j**2 - r * j)
    beta = 0.5 * dt * (sigma**2 * j**2 + r)
    gamma = 0.25 * dt * (sigma**2 * j**2 + r * j)

    a_left = -alpha[1:]
    b_left = 1 + beta
    c_left = -gamma[:-1]

    for n in range(N - 1, -1, -1):
        d = alpha * V[0 : M - 1] + (1 - beta) * V[1:M] + gamma * V[2 : M + 1]

        # Put boundary conditions
        boundary_0 = K * np.exp(-r * (N - n) * dt)  # PV of strike
        boundary_M = 0

        d[0] += alpha[0] * boundary_0
        d[-1] += gamma[-1] * boundary_M

        V[1:M] = thomas_algorithm(a_left, b_left, c_left, d)

        V[0] = boundary_0
        V[M] = boundary_M

    return S, V