Longstaff–Schwartz Pricing Checks for an Asian Bermudan Call

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


# 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: \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.")
```