American Option Pricing: Numerical Methods for Early Exercise Valuation
Introduction
American options, which permit exercise at any time before expiration, present one of the classic challenges in quantitative finance. Unlike European options with their elegant Black-Scholes closed-form solution, American options require numerical methods because the early exercise boundary is a free boundary that must be solved simultaneously with the option price. For the closed-form European case, our Options Calculator computes Black-Scholes pricing and Greeks directly for a given strike and expiration. This article covers three production-grade approaches: the Cox-Ross-Rubinstein binomial tree, the Longstaff-Schwartz least-squares Monte Carlo, and the finite difference method for the Black-Scholes PDE.The Early Exercise Problem
The value of an American option at any time $t$ is: $$ V_t = \max\left(h(S_t),\ \mathbb{E}^Q\left[e^{-r\Delta t} V_{t+\Delta t} \mid \mathcal{F}_t\right]\right) $$ where $h(S_t)$ is the intrinsic value (payoff from immediate exercise) and the expectation is the continuation value (holding the option). The option holder exercises when intrinsic value exceeds continuation value. For a put option, $h(S_t) = \max(K - S_t, 0)$. For a call on a non-dividend-paying stock, early exercise is never optimal (American call = European call). For calls on dividend-paying stocks and for all puts, the early exercise premium is material.Method 1: Binomial Tree (CRR)
The Cox-Ross-Rubinstein model discretizes the stock price process into a recombining binomial tree.Parameters
$$ u = e^{\sigma\sqrt{\Delta t}}, \quad d = e^{-\sigma\sqrt{\Delta t}} = \frac{1}{u}, \quad p = \frac{e^{(r-q)\Delta t} - d}{u - d} $$ where $\sigma$ is volatility, $r$ is the risk-free rate, $q$ is the dividend yield, and $\Delta t = T/N$.Implementation
``python
import numpy as np
def american_option_binomial(S: float, K: float, T: float, r: float,
sigma: float, q: float = 0.0,
N: int = 500, option_type: str = 'put') -> float:
"""
Price an American option using the CRR binomial tree.
Parameters
----------
S : float - Current stock price
K : float - Strike price
T : float - Time to expiration (years)
r : float - Risk-free rate
sigma : float - Volatility
q : float - Continuous dividend yield
N : int - Number of time steps
option_type : str - 'put' or 'call'
Returns
-------
float - Option price
"""
dt = T / N
u = np.exp(sigma * np.sqrt(dt))
d = 1 / u
p = (np.exp((r - q) * dt) - d) / (u - d)
disc = np.exp(-r * dt)
# Stock prices at maturity
stock_prices = S u np.arange(N, -1, -1) d np.arange(0, N + 1)
# Option values at maturity
if option_type == 'put':
option_values = np.maximum(K - stock_prices, 0)
else:
option_values = np.maximum(stock_prices - K, 0)
# Backward induction with early exercise check
for i in range(N - 1, -1, -1):
stock_prices = S u np.arange(i, -1, -1) d np.arange(0, i + 1)
# Continuation value
option_values = disc (p option_values[:-1] + (1 - p) * option_values[1:])
# Early exercise value
if option_type == 'put':
exercise = np.maximum(K - stock_prices, 0)
else:
exercise = np.maximum(stock_prices - K, 0)
# American: take max of continuation and exercise
option_values = np.maximum(option_values, exercise)
return option_values[0]
Example: American put
price = american_option_binomial(S=100, K=105, T=1.0, r=0.05, sigma=0.20, N=1000)
print(f"American Put Price: ${price:.4f}") # ~$8.72
`
With N=1000 steps, the binomial tree converges to within $0.01 of the true price. Richardson extrapolation (averaging results from N and N/2 steps) accelerates convergence.
Method 2: Longstaff-Schwartz Monte Carlo
The Longstaff-Schwartz (LSM) algorithm uses regression to estimate continuation values, enabling Monte Carlo simulation for American options.
Algorithm
- Simulate $M$ stock price paths over $N$ time steps
- At maturity, compute payoffs
- Working backward, at each time step:
- Identify paths that are in-the-money
- Regress discounted future cash flows on basis functions of the current stock price
- Compare the fitted continuation value with immediate exercise value
- Update cash flow matrix
`python
def american_option_lsm(S: float, K: float, T: float, r: float,
sigma: float, q: float = 0.0,
M: int = 100_000, N: int = 50,
option_type: str = 'put',
seed: int = 42) -> float:
"""
Longstaff-Schwartz least-squares Monte Carlo for American options.
"""
np.random.seed(seed)
dt = T / N
discount = np.exp(-r * dt)
# Simulate GBM paths
Z = np.random.standard_normal((M, N))
drift = (r - q - 0.5 sigma2) dt
diffusion = sigma * np.sqrt(dt)
log_returns = drift + diffusion * Z
log_prices = np.log(S) + np.cumsum(log_returns, axis=1)
stock_paths = np.exp(log_prices)
stock_paths = np.column_stack([np.full(M, S), stock_paths])
# Payoff function
if option_type == 'put':
payoff = lambda s: np.maximum(K - s, 0)
else:
payoff = lambda s: np.maximum(s - K, 0)
# Initialize cash flows at maturity
cash_flows = payoff(stock_paths[:, -1])
exercise_time = np.full(M, N)
# Backward induction
for t in range(N - 1, 0, -1):
S_t = stock_paths[:, t]
intrinsic = payoff(S_t)
itm = intrinsic > 0 # In-the-money paths
if itm.sum() < 10:
continue
# Discounted future cash flows for ITM paths
future_cf = np.zeros(M)
for i in range(M):
if itm[i]:
steps_ahead = exercise_time[i] - t
future_cf[i] = cash_flows[i] discount * steps_ahead
# Regression: continuation value ~ polynomial of stock price
X = S_t[itm]
Y = future_cf[itm]
# Laguerre polynomial basis (order 3)
basis = np.column_stack([
np.exp(-X / (2 * K)),
np.exp(-X / (2 K)) (1 - X / K),
np.exp(-X / (2 K)) (1 - 2 X / K + 0.5 (X / K)**2)
])
# OLS regression
coeffs = np.linalg.lstsq(basis, Y, rcond=None)[0]
continuation = basis @ coeffs
# Exercise if intrinsic > continuation
exercise_mask = np.zeros(M, dtype=bool)
exercise_mask[itm] = intrinsic[itm] > continuation
cash_flows[exercise_mask] = intrinsic[exercise_mask]
exercise_time[exercise_mask] = t
# Discount all cash flows to time 0
discounted = np.array([
cash_flows[i] np.exp(-r exercise_time[i] * dt)
for i in range(M)
])
return discounted.mean()
price = american_option_lsm(S=100, K=105, T=1.0, r=0.05, sigma=0.20)
print(f"LSM American Put: ${price:.4f}")
`
LSM is particularly valuable for multi-asset American options (e.g., rainbow options, basket options) where tree methods face exponential growth in dimensionality.
Method 3: Finite Difference Method
The finite difference approach directly solves the Black-Scholes PDE with the early exercise constraint.
The PDE
$$
\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + (r-q)S\frac{\partial V}{\partial S} - rV = 0
$$
subject to $V(S,t) \geq h(S)$ at all times (the American constraint).
Crank-Nicolson with PSOR
`python
def american_option_fd(S: float, K: float, T: float, r: float,
sigma: float, q: float = 0.0,
N_time: int = 1000, N_space: int = 200,
option_type: str = 'put') -> float:
"""
Finite difference (Crank-Nicolson + PSOR) for American options.
"""
S_max = 4 * K
dt = T / N_time
dS = S_max / N_space
S_grid = np.linspace(0, S_max, N_space + 1)
if option_type == 'put':
payoff = np.maximum(K - S_grid, 0)
else:
payoff = np.maximum(S_grid - K, 0)
V = payoff.copy()
# Coefficients for Crank-Nicolson
j = np.arange(1, N_space)
alpha = 0.25 dt (sigma2 j2 - (r - q) j)
beta = -0.5 dt (sigma2 * j2 + r)
gamma_coeff = 0.25 dt (sigma2 j2 + (r - q) j)
# Tridiagonal matrices
M1 = np.diag(1 - beta) + np.diag(-alpha[1:], -1) + np.diag(-gamma_coeff[:-1], 1)
M2 = np.diag(1 + beta) + np.diag(alpha[1:], -1) + np.diag(gamma_coeff[:-1], 1)
# Time stepping with Projected SOR
omega = 1.2 # Over-relaxation parameter
for n in range(N_time):
rhs = M2 @ V[1:-1]
# Boundary conditions
if option_type == 'put':
rhs[0] += alpha[0] (K np.exp(-r (N_time - n) dt))
else:
rhs[-1] += gamma_coeff[-1] (S_max - K np.exp(-r (N_time - n) dt))
# PSOR iteration
V_new = V[1:-1].copy()
for _ in range(100): # Max iterations
V_old = V_new.copy()
for i in range(len(V_new)):
residual = rhs[i] - M1[i] @ V_new
V_new[i] += omega * residual / M1[i, i]
V_new[i] = max(V_new[i], payoff[i + 1]) # American constraint
if np.max(np.abs(V_new - V_old)) < 1e-8:
break
V[1:-1] = V_new
# Interpolate to get price at S
return np.interp(S, S_grid, V)
``
Comparing the Methods
| Method | Complexity | Dimensions | Accuracy | Speed |
|---|---|---|---|---|
| Binomial Tree | O(N^2) | 1 | High (large N) | Fast |
| LSM Monte Carlo | O(M N) | Multi-asset | Moderate | Moderate |
| Finite Difference | O(N_t N_s) | 1-2 | Very high | Moderate |