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

Question

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: \begin{equation} A_t = \frac{1}{t}\int_0^t S_u \, \mathrm{d}u \end{equation} 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.")

early_exercise_premium = american_asian_call_option_price - european_asian_call_option_price
print(f"The early exercise premium is {round(early_exercise_premium, 4)} dollars.")