Why is there a difference in American option prices when comparing pricing methods (Python)?

Question

I have written a Python script to price American options using Least Squares Monte Carlo and added a QuantLib implementation below (analytical/binomial/finite difference) to compare. The problem is that my MCLS approach seems to slightly overprice calls and underprice puts and I can't seem to find the error in the code. Any help with this/advice on the best way to normalise the underlying's price would be greatly appreciated, thanks in advance!

    """ AMERICAN OPTION PRICING BY LEAST SQUARES MONTE CARLO, FINITE DIFFERENCE, ANALYTICAL AND BINOMIAL METHODS """

from numpy import zeros, concatenate, sqrt, exp, maximum, polyfit, polyval, shape, where, sum, argsort, random, \
    RankWarning, put, nonzero
from zlib import compress
import matplotlib.pyplot as plt
import os
import sys
from QuantLib import *
from pylab import *
import warnings
warnings.simplefilter('ignore', RankWarning)

plt.style.use('seaborn')

# Define global parameters
S0 = 100                                                       # Underlying price
K = 90                                                         # Strike
valuation_date = Date(1, 1, 2018)                              # Valuation date
expiry_date = Date(1, 1, 2019)                                 # Expiry date
t = ActualActual().yearFraction(valuation_date, expiry_date)   # Year fraction
T = 100                                                        # Time grid
dt = t / T                                                     # Delta time
r = 0.01                                                       # Interest rate
sig = 0.4                                                      # Volatility
sim = 10 ** 5                                                  # Number of MC simulations
DiscountFactor = exp(-r * dt)                                  # Discount factor

""" Least Squares Monte Carlo """


def GBM(underlying, time, simulations, rate, sigma, delta_t):  # Geometric Brownian Motion
    GBM = zeros((time + 1, simulations))
    GBM[0, :] = underlying
    for t in range(1, time + 1):
        brownian = standard_normal(simulations // 2)
        brownian = concatenate((brownian, -brownian))
        GBM[t, :] = (GBM[t - 1, :] * exp((rate - sigma ** 2 / 2.) * delta_t + sigma * brownian * sqrt(delta_t)))
    return GBM


def Payoff(strike, paths, simulations):  # Define option type and respective payoff
    if OptionType == 'call':
        po = maximum(paths - strike, zeros((T + 1, simulations)))
    elif OptionType == 'put':
        po = maximum(strike - paths, zeros((T + 1, simulations)))
    else:
        print('Incorrect input')
        os.execl(sys.executable, sys.executable, *sys.argv)
    return po


def loadingBar(count, total, size):  # MC progress bar
    percent = float(count) / float(total) * 100
    sys.stdout.write("\r" + str(int(count)).rjust(3, '0') + "/" + str(int(total)).rjust(3, '0') + ' [' + '=' * int(
        percent / 10) * size + ' ' * (10 - int(percent / 10)) * size + ']')


# Graph the regression fit and simulations
OptionType = str(input('Price call or put:'))
print('Plotting fitted regression at T...')
GBM = GBM(S0, T, sim, r, sig, dt)
payoff = Payoff(K, GBM, sim)
ValueMatrix = zeros_like(payoff)
ValueMatrix[T, :] = payoff[T, :]
prices = GBM[T, :]
value = ValueMatrix[T, :]
regression = polyfit(prices, value * DiscountFactor, 4)
ContinuationValue = polyval(regression, prices)
sorted_index = argsort(prices)
prices = prices[sorted_index]
ContinuationValue = ContinuationValue[sorted_index]

ValueMatrix[T, :] = where(payoff[T, :] > ContinuationValue, payoff[T, :], ValueMatrix[T, :] * DiscountFactor)
ValueVector = ValueMatrix[T, :] * DiscountFactor
ValueVector = ValueVector[sorted_index]

plt.figure()
f, axes = plt.subplots(2, 1)
axes[0].set_title('American Option')
axes[0].plot(prices, ContinuationValue, label='Fitted Polynomial')
axes[0].plot(prices, ValueVector, label='Inner Value')
axes[0].set_ylabel('Payoff')
axes[0].set_xlabel('Asset Price')
axes[0].legend()
axes[1].set_title('Geometric Brownian Motion')
axes[1].plot(GBM, lw=0.5)
axes[1].set_ylabel('Asset Price')
axes[1].set_xlabel('Time')
f.tight_layout()
plt.show()

# MC results
print('Pricing option...')
for i in range(0, 100):
    loadingBar(i, 100, 2)
    for t in range(T - 1, 0, -1):
        ITM = payoff[t, :] > 0
        ITMS = compress(ITM, GBM[t, :])
        ITMP = compress(ITM, payoff[t + 1, :] * DiscountFactor)
        regression = polyval(polyfit(ITMS, ITMP, 4), ITMS)
        continuation = zeros(sim)
        put(continuation, nonzero(ITM), regression)
        payoff[t, :] = where(payoff[t, :] > continuation, payoff[t, :], payoff[t + 1, :] * DiscountFactor)
        price = sum(payoff[1, :] * DiscountFactor) / sim
print('\nLeast Squares Monte Carlo Price:', price)


""" QuantLib Pricing """

S0 = SimpleQuote(S0)
if OptionType == 'call':
    OptionType = Option.Call
elif OptionType == 'put':
    OptionType = Option.Put
else:
    print('Incorrect input')
    os.execl(sys.executable, sys.executable, *sys.argv)


def Process(valuation_date, r, dividend_rate, sigma, underlying):
    calendar = UnitedStates()
    day_counter = ActualActual()
    Settings.instance().evaluation_date = valuation_date
    interest_curve = FlatForward(valuation_date, r, day_counter)
    dividend_curve = FlatForward(valuation_date, dividend_rate, day_counter)
    volatility_curve = BlackConstantVol(valuation_date, calendar, sigma, day_counter)
    u = QuoteHandle(underlying)
    d = YieldTermStructureHandle(dividend_curve)
    r = YieldTermStructureHandle(interest_curve)
    v = BlackVolTermStructureHandle(volatility_curve)
    return BlackScholesMertonProcess(u, d, r, v)


def FDAmericanOption(valuation_date, expiry_date, OptionType, K, process):  # Finite difference
    exercise = AmericanExercise(valuation_date, expiry_date)
    payoff = PlainVanillaPayoff(OptionType, K)
    option = VanillaOption(payoff, exercise)
    time_steps = 100
    grid_points = 100
    engine = FDAmericanEngine(process, time_steps, grid_points)
    option.setPricingEngine(engine)
    return option


def ANAmericanOption(valuation_date, expiry_date, OptionType, K, process):  # Analytical
    exercise = AmericanExercise(valuation_date, expiry_date)
    payoff = PlainVanillaPayoff(OptionType, K)
    option = VanillaOption(payoff, exercise)
    engine = BaroneAdesiWhaleyEngine(process)
    option.setPricingEngine(engine)
    return option


def BINAmericanOption(valuation_date, expiry_date, OptionType, K, process):  # Binomial
    exercise = AmericanExercise(valuation_date, expiry_date)
    payoff = PlainVanillaPayoff(OptionType, K)
    option = VanillaOption(payoff, exercise)
    timeSteps = 10 ** 3
    engine = BinomialVanillaEngine(process, 'crr', timeSteps)
    option.setPricingEngine(engine)
    return option


def FDAmericanResults(option):
    print('Finite Difference Price: ', option.NPV())
    print('Option Delta: ', option.delta())
    print('Option Gamma: ', option.gamma())


def ANAmericanResults(option):
    print('Barone-Adesi-Whaley Analytical Price: ', option.NPV())


def BINAmericanResults(option):
    print('Binomial CRR Price: ', option.NPV())


# Quantlib results
process = Process(valuation_date, r, 0, sig, S0)
ANoption = ANAmericanOption(valuation_date, expiry_date, OptionType, K, process)
ANAmericanResults(ANoption)
BINoption = BINAmericanOption(valuation_date, expiry_date, OptionType, K, process)
BINAmericanResults(BINoption)
FDoption = FDAmericanOption(valuation_date, expiry_date, OptionType, K, process)
FDAmericanResults(FDoption)

os.system('say "completo"')

Answer 1

I've never dealt with Python, so I am just trying to understand what's going on, visually/logically. Overall it looks fine, apart from the fact that LS suggest to only use the in-the-money paths for the regression, can't see that in your code, or am I missing it? A few side notes, your variable naming is a little non-standard, I would change your T -> NT, t -> T. Also what you call discount rate is really usually called a discount factor. Finally you are trying to price an American(Bermudan) call with no dividends, which should have the same price as the equivalent European call.

The LS algorithm, if done properly, should underprice your Americans anyway. That's because the continuation value approximation via the basis functions, is just that, an approximation. Which means that the algo (which bases its decision on the cont. value) will not always take the correct (optimal) decision to exercise, which then means the option value will be slightly less than if you had always exercised optimally.

Now, as Bob Jansen says, it may well be that you simply have too much noise. I mean, who uses 10.000 MC paths nowadays! I'm guessing Python is way too slow for MC. So, I tried your problem with this tool and with 8.000 (Sobol) paths (for both the main simulation and the regression) and 6 basis functions and I get C=21.32 (in half a sec). Then tried with 131.000 paths and I get 21.02 (in 12 secs). The correct price is 21.046. So yes, it may well be that you just have too much noise with so few paths and you cannot draw any conclusion as to whether you have some other problem or not.