Implementing Longstaff–Schwartz Monte Carlo for American Puts
Summary
The document describes an implementation problem with Least Squares Monte Carlo (LSM) pricing of American vanilla puts. The questioner follows Longstaff and Schwartz’s setup, including simulated paths with antithetic variates and regressions using weighted Laguerre basis functions, but obtains prices below those reported in the paper. Simpler polynomial bases and ordinary Laguerre functions produce reasonable results in the same implementation, prompting questions about the weighted basis and why many open implementations use powers of the underlying instead.
The code simulates asset paths under a geometric Brownian motion assumption, works backward through exercise dates, regresses discounted future values on in-the-money paths, and compares continuation value with immediate exercise value. However, the document contains no replies or diagnosis, so it does not establish the cause of the price discrepancy or confirm a matching implementation. It is best read as a concrete LSM setup and unresolved debugging question; details such as basis scaling and implementation conventions would need investigation.
Key ideas
- LSM estimates an American option’s continuation value by regression on simulated in-the-money paths.
- The example compares immediate exercise payoff with a fitted continuation value while stepping backward through time.
- The implementation uses antithetic paths and compares polynomial, Laguerre, and weighted Laguerre bases.
- The reported price mismatch is unresolved, and the document provides no validated fix.
Tags
Full text
# Struggling with Implementation of Longstaff-Schwartz (2001) LSM Method
# Struggling with Implementation of Longstaff-Schwartz (2001) LSM Method
I'm having difficulties implementing the Least Squares Monte Carlo (LSM) method from Longstaff and Schwartz's 2001 paper for American Vanilla Put option pricing. According to the paper, they used 100,000 simulated paths (50,000 plus 50,000 antithetic) and the first 3 weighted Laguerre polynomials with a constant term in the regressions. The option values I'm getting show significantly smaller prices compared to the results reported in the original Longstaff-Schwartz (2001) paper.
I've tried following the paper's methodology, but my results seem inconsistent with their published values.
- Has anyone successfully implemented the method using weighted Laguerre polynomials with results matching the paper?
- Why do most open-source implementations prefer simple polynomials (X, X^2, X^3)?
Note: I believe there's no issue with the algorithm itself since I get reasonable results using simple polynomials or standard Laguerre polynomials. The implementation works well with these basis functions. I'm using Python for implementation, but I'm open to suggestions in other languages.
```
def Vanilla_option_lsmc(S0, K, T, r, sigma, n_steps, n_paths, M, model='polynomial', option_type='call'):
dt = T / n_steps
if option_type == 'call':
otype = 1
else:
otype = -1
# path - antithetic
random = np.random.standard_normal((int(n_paths/2), n_steps))
S_ = S0 * np.exp(np.cumsum((r - 0.5 * sigma**2) * dt +
sigma * np.sqrt(dt) * random, axis=1))
S__ = S0 * np.exp(np.cumsum((r - 0.5 * sigma**2) * dt +
sigma * np.sqrt(dt) * -1 * random, axis=1))
S = np.append(S_, S__, axis=0)
S = np.insert(S, 0, S0, axis=1)
# value
V = np.zeros((n_paths, n_steps + 1))
V[:, -1] = np.maximum(otype * (S[:, -1] - K), 0) # Payoff
for t in range(n_steps - 1, 0, -1):
conditions = np.where(otype * (S[:, t] - K) > 0)[0] # ITM Conditions
# continuation value
x = S[:, t][conditions]
y = V[:, t+1][conditions] * np.exp(-r * dt)
if model == 'polynomial':
X = polynomial_design_matrix(x, M)
elif model == 'laguerre':
X = laguerre_design_matrix(x, M)
elif model == 'weighted laguerre':
X = weighted_laguerre_design_matrix(x, M)
else:
X = weighted_laguerre_design_matrix(x, M)
# OLS
beta_hat = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
continuation_value = np.zeros(n_paths)
continuation_value[conditions] = np.dot(X, beta_hat)
# exercise value
exercise_value = np.maximum(otype * (S[:, t] - K), 0)
# decision
V[:, t] = np.where(exercise_value > continuation_value,
exercise_value,
V[:, t+1] * np.exp(-r * dt))
option_price = V[:, 1] * np.exp(-r * dt)
price_mean = np.mean(option_price)
price_var = np.std(option_price)/np.sqrt(n_paths)
return price_mean, price_var
def polynomial_design_matrix(x, M):
if len(x.shape) == 1:
x = x.reshape(1, -1)
num_vars, len_vars = x.shape
one = np.ones((len_vars, 1))
square_sum_terms = []
for j in range(num_vars):
for deg in range(1, M):
square_sum_terms.append((x[j, :] ** deg)[:, np.newaxis])
X = np.concatenate([one] + square_sum_terms, axis=1)
return X
def laguerre_design_matrix(x, M):
# Preallocate the basis array
basis = np.zeros((len(x), M))
# Initialize the first two Laguerre polynomials
basis[:, 0] = 1 # L_0(x) = 1
if M > 1:
basis[:, 1] = 1 - x # L_1(x) = 1 - x
# Compute higher-order Laguerre polynomials iteratively
for i in range(2, M):
basis[:, i] = ((2 * (i - 1) + 1 - x) * basis[:, i - 1] - (i - 1) * basis[:, i - 2]) / i
return basis
def weighted_laguerre_design_matrix(x, M):
# Preallocate the basis array
basis = np.zeros((len(x), M))
# Initialize the first two Laguerre polynomials
basis[:, 0] = 1
deg = M - 1
if deg > 1:
basis[:, 1] = np.exp(-x/2) # L_0(x) = 1
if deg > 2:
basis[:, 2] = np.exp(-x/2)* (1 - x) # L_1(x) = 1 - x
# Compute higher-order Laguerre polynomials iteratively
for i in range(2, deg):
basis[:, i+1] = ((2 * (i - 1) + 1 - x) * basis[:, i] - (i - 1) * basis[:, i - 1]) / i
return basis
```Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.