# Longstaff & Schwartz algorithm - Python: American option cheaper than European option

I have implemented the Longstaff & Schwartz algorithm for pricing American Option in Python, but I ran into an issue while doing some experiments: sometimes, for the same option, I get a higher price for the European option (compared to the American one), which make me think that I made a mistake, but I do not know where.

In the below reproducible example, I get a American Option price of $$21.517$$ dollars, and an European Option price of $$21.847$$ dollars. This looks absurd (as the early exercise premium is negative, $$=-0.33$$ dollars) but I could not find the error. I only noticed that the intermediate american option price is not increasing with the loop, which is strange.

Could you help me?

## Reproducible example

We take the example from the original paper, of an Asian American Bermuda call option. Let's consider a non-markovian option: the American-Bermuda-Asian Call Option. This option, of maturity $$2$$ years, pays $$(A_t - K)_+$$ at time $$0.25 \leq t \leq 2$$ (and $$0$$ otherwise), where: $$$$A_t = \frac{1}{t}\int_0^t S_u \, \mathrm{d}u$$$$ with $$S_u$$ the underlying value at time $$u$$.

import numpy as np
from numpy.polynomial.laguerre import Laguerre

# Step 1
maturity = 2 # maturity of the call option
strike = 100 # strike of the call option
sigma = 0.2 # volatility of the underlying asset
s0 = 110 # price of the asset at t = 0
N = 6 # number of simulated paths
times = np.array([0., 0.25, 0.5,0.75, 1., 1.25, 1.5, 1.75, 2.]).round(5) # times to record the simulation
american_exercising_times = np.array([0.25, 0.5,0.75, 1., 1.25, 1.5, 1.75, 2.]).round(5)
european_exercising_times = np.array([0.,2.]).round(5)
M = 3 # number of Laguerre polynomials to consider

def asian_call_payoff(paths, times, strike):
A = np.trapz(paths.T, times)/(times[-1] - times[0])

return np.maximum(A - strike, 0)

def spot(t1, t0, s0, r, sigma):
"""
Draws a value for S_t knowing S_0 = s0 at t0.
"""

W = np.random.normal(0, np.sqrt(t1-t0), s0.size)
st = s0 * np.exp((r-sigma**2/2)*(t1-t0) + sigma*W)

return st

def generatePaths(times, s0, N, r, sigma):
"""
Generates N sample paths and returns these paths under the form of a numpy array whose
entries are values of S at each time point times[i], starting from s = s0 at times[0].

Returns
-------
S: np.array of shape (times.size, N).
"""

sample_paths = np.zeros((times.size, N))
sample_paths[0] = s0
for i in range(1, times.size):
sample_paths[i] = spot(times[i], times[i-1], sample_paths[i-1], r, sigma)

return sample_paths

# Step 2
r = 0.06
sigma = 0.2

np.random.seed(10) # For seed = 10, early exercise premium is negative !!!
paths = generatePaths(times = times,
s0 = s0,
N = N,
r = r,
sigma = sigma)
print(f"Monte Carlo paths generated for s0 = {s0}, N = {N}, r = {r}, sigma = {sigma}.")


Pricing of the American option:

exercising_times = american_exercising_times.copy()
discounting_factors = np.exp(-r * (exercising_times - times[0])) # we are looking from time 0.

T = len(exercising_times)

exercising_indexes = np.where(np.in1d(times, exercising_times))[0]
cash_flows = np.zeros((T, N))
cash_flows[T-1] = asian_call_payoff(paths[:exercising_indexes[T-1]+1], times[:exercising_indexes[T-1]+1], strike)
print(f'Option price at iteration {T-1}: {(discounting_factors@cash_flows).mean()}\n')

for t in range(T-2, 0, -1):
ind_t = exercising_indexes[t]
ITM_paths = asian_call_payoff(paths[:ind_t+1], times[:ind_t+1], strike) > 0

X = paths[ind_t, ITM_paths]
y_not_discounted = cash_flows[t+1:, ITM_paths]
y = np.exp(- r * (exercising_times[t+1:] - exercising_times[t]))@y_not_discounted # discount factors over period T to t

cond_exp = Laguerre(0).fit(X, y, deg=M)
delayed_exercise = cond_exp(X)
early_exercise = asian_call_payoff(paths[:ind_t+1, ITM_paths], times[:ind_t+1], strike)

# Among the ITM paths, these are the paths were immediate exercise is optimal
early_is_optimal = early_exercise > delayed_exercise
ITM_early_paths = ITM_paths.copy()
ITM_early_paths[ITM_paths] = early_is_optimal

cash_flows[t, ITM_early_paths] = early_exercise[early_is_optimal]
cash_flows[t+1:, ITM_early_paths] = 0
print(f'Option price at iteration {t}: {(discounting_factors@cash_flows).mean()}\n')

american_cash_flows = cash_flows.copy()
american_asian_call_option_price = (discounting_factors@cash_flows).mean()
print(f"\nSettings: s0 = {s0}, sigma = {sigma}, maturity = {maturity}, risk-free rate = {r}")
print(f"The final American Asian Call Option price is: {american_asian_call_option_price} dollars.")


Pricing of the European option:

exercising_times = european_exercising_times.copy()
discounting_factors = np.exp(-r * (exercising_times - times[0])) # we are looking from time 0.
T = len(exercising_times)

exercising_indexes = np.where(np.in1d(times, exercising_times))[0]
cash_flows = np.zeros((T, N))
cash_flows[T-1] = asian_call_payoff(paths[:exercising_indexes[T-1]+1], times[:exercising_indexes[T-1]+1], strike)

european_cash_flows = cash_flows.copy()
european_asian_call_option_price = (discounting_factors@cash_flows).mean()
print(f"\nSettings: s0 = {s0}, sigma = {sigma}, maturity = {maturity}, risk-free rate = {r}")
print(f"The final European Asian Call Option price is: {european_asian_call_option_price} dollars.")