Numerical Approximation of the Exercise Boundary of American Options

Problem Setup

An American put option gives the holder the right to sell one unit of an asset at strike price $E$ at any time $t \in [0, T]$. Unlike a European option, early exercise may be optimal when the asset price $S$ is sufficiently low. Let $\tau = T - t$ denote backward time (so $\tau=0$ at maturity). The option value $V(S,\tau)$ satisfies the Black–Scholes PDE in the region where exercise is not yet optimal:

Here $\sigma$ is volatility and $r$ is the risk-free interest rate. The payoff at maturity ($\tau = 0$) is $V^* = \max(E - S,\, 0)$, and the option must satisfy $V(S,\tau) \geq V^*$ at all times.

The Free Boundary

The domain splits into two regions separated by the unknown free boundary $S^*(\tau)$:

The free boundary $S^*(\tau)$ is the asset price below which immediate exercise is optimal. It is unknown and must be determined as part of the solution — this is what makes the problem a free boundary problem. Economically, $S^*(\tau)$ is monotone decreasing in $\tau$ (the boundary moves left as $\tau$ increases toward $T$).

The Penalty Method

Rather than tracking $S^*$ explicitly during the solve, we enforce the early-exercise constraint by adding a penalty term to the PDE [Forsyth & Vetzal, 2002]. The penalized PDE is:

where $\rho \gg 1$ is a large penalty parameter (we use $\rho = 10^7$).

The logic is elegant:

• If $V \geq E - S$: the constraint is satisfied, so $\max(E-S-V,0) = 0$ and the PDE is unchanged.
• If $V < E - S$: the term $\rho(E-S-V)$ becomes large and positive, driving $V$ back above the payoff floor.

We use Picard (penalty) iteration at each time step to handle the nonlinearity. The iteration converges when the residual is smaller than $\varepsilon_{\mathrm{tol}} \approx \rho^{-1} = 10^{-7}$.

Finite Difference Discretization

Spatial discretization on a nonuniform grid

We use a nonuniform grid in $S$, concentrated near the strike $E$ where the solution has a kink. At an interior node $S_i$ with left spacing $h_L = S_i - S_{i-1}$ and right spacing $h_R = S_{i+1} - S_i$, the central finite difference approximations are:

The grid has $N_x = 2^{k-1} \times 33$ spatial nodes at refinement level $k$, with nodes clustered near $S = E$ via a stretching parameter. Parameters: $E=100,\ r=0.05,\ \sigma=0.8,\ T=0.5,\ S \in [0,\,12E]$.

Crank–Nicolson time stepping

With $\theta = 1/2$, the time discretization is:

This is second-order accurate in both space and time.

Rannacher start-up and adaptive stepping

The non-smooth initial condition $V^* = \max(E-S,0)$ can degrade accuracy. Two techniques address this:

• Smoothing of initial conditions: The payoff is replaced by a smooth approximation (Phi1 or Phi2) to avoid order reduction [Kreiss et al., 1970].
• Rannacher start-up: The first 4 half-steps use fully implicit ($\theta=1$) stepping with $\Delta t/2$, damping oscillations from the kink before switching to Crank–Nicolson.
• Adaptive step selection: After step 4, the time step is chosen adaptively. At each node $j$, we compute a scale factor $$\lambda_j = \frac{D \cdot \max\!\bigl(R,\,|V_j^n|,\,|V_j^{n-1}|\bigr)}{|V_j^n - V_j^{n-1}|}$$ and set $\Delta t_{n+1} = \min_j(\lambda_j)\cdot\Delta t_n$, where $D = D_0/2^{k-1}$ is halved at each refinement level $k$.

Free Boundary Detection

Once the PDE is solved at a given time step, we extract $S^*(\tau)$ by post-processing the numerical solution. The workflow is:

Locating the interval

We find the leftmost index $i_0$ such that $V_{i_0} - \mathrm{payoff}_{i_0} \geq 2\,\varepsilon_{\mathrm{tol}}$. This is the first grid node in the continuation region. The stencil points for interpolation are then taken as $\{i_0, i_1, i_2\} = \{i_0, i_0+1, i_0+2\}$ (or shifted by one: $\{i_0+1, i_0+2, i_0+3\}$).

Free boundary conditions

At the true free boundary $S^*(\tau)$, two conditions hold simultaneously:

Value matching gives a double root when $V(S) - (E-S)$ is minimized, which makes fsolve unstable. Smooth pasting gives a single root and is much more stable numerically. All nine schemes below exploit the smooth pasting condition.

Specifically, if $p(S)$ is a local interpolant of $V$ (or of $V' = \partial V/\partial S$), we solve:

using MATLAB's fsolve with tolerance $10^{-9}$.

Nine Interpolation Schemes

We compared nine schemes organized into two families: those that interpolate $V$ directly (solving $p'(S^*)+1=0$), and those that interpolate $V' = \partial V/\partial S$ (solving $p(S^*)+1=0$). The derivative values $V'_i$ and $V''_i$ at grid nodes are obtained from the finite difference matrices applied to the PDE solution.

Scheme 6 in detail: Quadratic Lagrange on $V'$

Given three stencil nodes $S_1, S_2, S_3$ with derivative values $V'_1, V'_2, V'_3$ from the PDE solution, we fit the quadratic Lagrange interpolant:

Applying the smooth pasting condition $p(S^*) = -1$ gives:

This quadratic in $S^*$ is solved by fsolve with initial guess $S_1$ (index $i_0$).

Schemes 7 & 8 in detail: PCHIP

PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) is $C^1$, shape-preserving, and avoids Runge oscillations. Given nodes $S_{\mathrm{vec}} = (S_1,S_2,S_3)$:

• Scheme 7 (P3vC): Build $p(S) = \mathrm{pchip}(S_{\mathrm{vec}},\, V_{\mathrm{vec}})$, differentiate via fnder, and solve $p'(S^*)+1=0$.
• Scheme 8 (P3dvC): Build $p(S) = \mathrm{pchip}(S_{\mathrm{vec}},\, V'_{\mathrm{vec}})$ directly on the derivative values, and solve $p(S^*)+1=0$.

Scheme 9 in detail: Quintic Hermite–Birkhoff

Using data $V_j, V'_j, V''_j$ at both endpoints $S_1, S_2$ gives six interpolation conditions, determining a degree-5 polynomial. Define $L_j(S) = (S - S_{3-j})/(S_j - S_{3-j})$. The basis functions are:

The quintic interpolant is:

The free boundary satisfies $p'_5(S^*) + 1 = 0$.

Numerical Results

PDE solution convergence

The underlying PDE solver achieves clean second-order convergence in both option value (ErrV) and delta (ErrD). At refinement $k=6$ ($N_x=2112$) the option value error is below $5 \times 10^{-5}$.

Free boundary convergence: Schemes 6 and 8

Both derivative-based schemes (L3dvQ and P3dvC) achieve zero non-monotone timesteps across all refinement levels. The free boundary estimates at $k=5$ ($N_x = 528$) agree to five significant figures.

Comparison with literature (different parameters)

Using $K=100,\,r=0.1,\,\sigma=0.3,\,T=1$ year, our schemes match published reference values closely:

Summary across all nine schemes

Conclusions and Future Work

We implemented and compared nine interpolation-based schemes for locating the exercise boundary $S^*(\tau)$ in American put option pricing. All schemes use the smooth-pasting condition $\partial V/\partial S = -1$ (rather than value matching) for stability. The main takeaways are:

• The derivative family (schemes that interpolate $V'$) is more robust than the value family (schemes that interpolate $V$ and differentiate) because no numerical differentiation of the interpolant is needed.
• L3dvQ (quadratic Lagrange on $V'$) and P3dvC (PCHIP on $V'$) achieve zero non-monotone timesteps and closely match published reference values.
• Starting the stencil at $i_0$ (leftmost continuation node) gives smaller absolute error than $i_0+1$ at moderate refinement.
• Higher-order interpolants (quintic Hermite, degree 5) do not necessarily improve results — scheme 9 is the worst performer due to instability near the boundary.

Future directions:

• Analyze the residual $V_\tau - \mathcal{L}V$ near the free boundary and use its sign change to estimate $S^*$ more directly.
• High-order spatial discretizations (fourth-order compact differences).
• Extension to American call options and other payoff structures.

References

• Forsyth, P. A. and Vetzal, K. R. (2002). Quadratic convergence for valuing American options using a penalty method. SIAM Journal on Scientific Computing, 23(6):2095–2122.
• Kreiss, H.-O., Thomée, V., and Widlund, O. (1970). Smoothing of initial data and rates of convergence for parabolic difference equations. Communications on Pure and Applied Mathematics, 23(2):241–259.
• Christara, C. and Leung, N. C.-H. (2018). Analysis of quantization error in financial pricing via finite difference methods. SIAM Journal on Numerical Analysis, 56(3):1731–1757.
• Wang, D., Serkh, K., and Christara, C. (2023). A high-order deferred correction method for the solution of free boundary problems using penalty iteration, with an application to American option pricing. Journal of Computational and Applied Mathematics, 432:115272.
• Kim, B. J., Ahn, C., and Choe, H. J. (2013). Direct computation for American put option and free boundary using finite difference method. Japan journal of industrial and applied mathematics, 30(1):21–37.
• Wu, L. and Kwok, Y.-K. (1997). A front-fixing finite difference method for the valuation of American options. Journal of Financial Engineering, 6(4):83–97.
• Zhu, S.-P. (2006). An exact and explicit solution for the valuation of American put options. Quantitative Finance, 6(3):229–242.