Longstaff-Schwartz Algorithm: American Options in Python

In January 2001 Francis Longstaff and Eduardo Schwartz published Valuing American Options by Simulation: A Simple Least-Squares Approach. Their difficult example was not a vanilla put. It was a cancellable index amortising swap driven by a 20-factor string model of the term structure.

That choice made the point. Trees and finite-difference grids are excellent in one or two state variables, then suffer from dimensionality. Monte Carlo handles many factors but runs forwards, while exercise decisions must be solved backwards. Least-squares Monte Carlo, usually LSMC or LSM, joins the two.

Our Monte Carlo guide, option-pricing overview and Black-Scholes article cover the prerequisites.

Optimal stopping

For exercise dates $0=t_0<t_1<\cdots<t_M=T$ and payoff $h_i(X_i)$, the Bermudan value is

$$V_0=\sup_{\tau\in\{t_0,\ldots,t_M\}} \mathbb E^{\mathbb Q}\!\left[ e^{-\int_0^\tau r_sds}h_\tau(X_\tau) \right].$$

$X_i$ is the full Markov state needed to value the contract, not necessarily spot alone. For an Asian Bermudan it includes the running average; for a swaption it may include several curve factors.

At date $t_i$, the continuation value is

$$C_i(X_i)= \mathbb E^{\mathbb Q}\!\left[ D_{i,i+1}V_{i+1}(X_{i+1})\mid X_i \right],$$

where $D_{i,i+1}$ is the pathwise discount factor. Bellman's rule is simply

$$V_i(X_i)=\max\{h_i(X_i),C_i(X_i)\}.$$

The conditional expectation is unknown. Longstaff-Schwartz estimates it by regressing realised discounted future cashflows on basis functions of the state.

Backwards regression

Simulate risk-neutral paths through every exercise date and set the final cashflow to $h_M$. Then work backwards:

1. Discount each path's current future cashflow from $t_{i+1}$ to $t_i$.
2. On paths that are in the money at $t_i$, regress those discounted cashflows on basis functions $\psi_j(X_i)$.
3. Exercise where immediate value exceeds the fitted continuation value.
4. Replace the path's future cashflow with immediate exercise value where exercise occurs; otherwise retain the discounted future cashflow.

In symbols, ordinary least squares chooses

$$\widehat\beta_i =\arg\min_\beta \sum_{n\in\mathcal I_i} \left(Y_i^{(n)}-\sum_{j=0}^{J}\beta_j\psi_j(X_i^{(n)})\right)^2,$$

where $\mathcal I_i$ is the in-the-money set and $Y_i^{(n)}$ is the realised value from following the already-estimated policy at later dates.

Regressing on in-the-money paths follows the original algorithm and concentrates observations near the exercise decision. It is not a universal theorem that all-path regression is invalid. For some payoffs, alternative sampling, weighted regression or boundary-focused methods perform better. The test is out-of-sample policy value and boundary stability.

A correct Python implementation

Raw powers of a £100 spot create an ill-conditioned Vandermonde matrix. The implementation below scales spot by strike, fits on one path set and values the learned policy on an independent set. That separation turns the reported policy value into a valid lower-bound estimate, up to Monte Carlo error.

import numpy as np

def gbm_paths(rng, n, S0, T, r, q, vol, steps):
    dt = T / steps
    z = rng.standard_normal((n, steps))
    dx = (r - q - 0.5*vol**2)*dt + vol*np.sqrt(dt)*z
    future = np.exp(np.log(S0) + np.cumsum(dx, axis=1))
    return np.column_stack((np.full(n, S0), future))

def american_put_lsmc(S0, K, T, r, q, vol, steps=50,
                      n_train=250_000, n_price=500_000, degree=3, seed=42):
    rng, disc = np.random.default_rng(seed), np.exp(-r*T/steps)
    train = gbm_paths(rng, n_train, S0, T, r, q, vol, steps)
    value, policy = np.maximum(K - train[:, -1], 0.0), {}

    for i in range(steps - 1, 0, -1):
        y, x = disc*value, train[:, i]/K
        exercise = np.maximum(K - train[:, i], 0.0)
        itm = exercise > 0
        if itm.sum() <= degree:
            policy[i], value = None, y
            continue
        policy[i] = np.polynomial.polynomial.polyfit(x[itm], y[itm], degree)
        continuation = np.polynomial.polynomial.polyval(x, policy[i])
        value = np.where(itm & (exercise > continuation), exercise, y)

    paths = gbm_paths(rng, n_price, S0, T, r, q, vol, steps)
    value = np.maximum(K - paths[:, -1], 0.0)
    for i in range(steps - 1, 0, -1):
        y, x = disc*value, paths[:, i]/K
        exercise = np.maximum(K - paths[:, i], 0.0)
        if policy[i] is None:
            value = y
            continue
        continuation = np.polynomial.polynomial.polyval(x, policy[i])
        value = np.where((exercise > 0) & (exercise > continuation), exercise, y)

    pv = disc*value
    continuation_0 = pv.mean()
    price = max(K - S0, continuation_0)
    stderr = 0.0 if K - S0 >= continuation_0 else pv.std(ddof=1)/np.sqrt(n_price)
    return price, stderr

The code assumes constant $r$, $q$ and volatility and allows exercise on the 50 grid dates. It deliberately excludes $t_0$ from regression because every path has the same state there; immediate exercise is compared with the estimated continuation value at the end.

For $S_0=K=100$, $T=1$, $r=5\%$, $q=0$ and $\sigma=20\%$, the stated seed gives approximately $6.06$ with Monte Carlo standard error $0.01$. A 2,000-step Cox-Ross-Rubinstein tree gives $6.0900$. The gap is plausible for a lower-bound policy estimate with a cubic basis; claiming automatic agreement to $0.1\%$ would be false precision.

Bias and validation

There are two effects that are often conflated.

Policy suboptimality gives a lower bound. Freeze the fitted policy, apply it to fresh paths and average the discounted payoff. Any feasible stopping rule is worth no more than the optimal one.

Reusing training paths creates look-ahead bias. The same observations influence both the fitted continuation value and their own exercise decision. This can bias an in-sample LSMC estimate upwards, while policy approximation pushes downwards. The net sign is not guaranteed. Independent valuation paths remove this reuse; leave-one-out LSMC is another approach.

An independent policy valuation is not an upper bound. Andersen and Broadie's primal-dual method constructs an upper bound from an approximate martingale, usually with nested simulation. A credible implementation reports the lower estimate, its sampling error and, for material books, a dual upper bound.

Validation should also include:

• convergence over path count, exercise-grid spacing and basis family;
• comparison with a tree or PDE in a low-dimensional special case;
• European lower bounds and American call identities where applicable;
• exercise-boundary plots, not only one headline price;
• independent random seeds and confidence intervals.

Basis and state design

Low-order polynomials, Laguerre functions and Hermite functions all work on textbook puts. Scaling and orthogonalisation matter more than the label. Add terms until the out-of-sample value stops improving, rather than until in-sample $R^2$ looks impressive.

The state must carry all information relevant to continuation. Regressing a path-dependent payoff on spot alone violates the Markov requirement and can produce a smooth but wrong boundary. For an interest-rate option, suitable curve state variables replace spot; the Vasicek model is a simple one-factor example.

LSMC avoids a spatial grid, but it does not defeat dimensionality. Polynomial tensor bases grow combinatorially, relevant states may be rare near the boundary and nested upper bounds are expensive. Sparse bases, local regression and neural networks can help, but each adds model-selection risk.

Where LSMC is the wrong tool

For a single-stock American put under Black-Scholes, a tree or finite-difference solver is faster, more accurate and easier to audit. LSMC earns its keep with several risk factors, path dependence or simulation-only dynamics.

It also prices only the exercise dates supplied. A Bermudan schedule is exact in that respect. Calling a fine Bermudan grid “American” introduces time-discretisation error, which must be measured rather than assumed negligible.

Frequently Asked Questions

What does LSMC stand for?

It stands for least-squares Monte Carlo, the regression-based simulation method introduced by Longstaff and Schwartz in 2001.

What is the regression target?

It is the realised discounted future cashflow from following the exercise policy already estimated at later dates, conditional on the current state.

Is an LSMC price always biased low?

An out-of-sample value from a fixed feasible policy is a lower bound apart from sampling error. An in-sample estimate can also contain upward look-ahead bias, so its net bias has no guaranteed sign.

Does a fresh valuation simulation give an upper bound?

No. It gives a cleaner lower-bound policy value. A dual martingale construction such as Andersen-Broadie is needed for a genuine upper bound.

Can LSMC handle many assets?

Yes, without a full spatial grid, but regression complexity and sparse coverage still worsen with dimension. The state and basis must represent cross-asset and path-dependent effects.

Which sources should I use?

Read the original Longstaff-Schwartz paper, Clément, Lamberton and Protter's convergence analysis and Andersen-Broadie's dual bounds. Simulation is the easy half. A defensible stopping policy is the product.