Return to Question

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

import numpy as np
import matplotlib.pyplot as plt
import os
import sys
from QuantLib import *

plt.style.use('seaborn')

# Define global parameters

 S0 = 100
K = 90
valuation_date = Date(1, 1, 2018)
expiry_date = Date(1, 1, 2019)
t = (expiry_date - valuation_date) / 365
T = 100
dt = t / T
r = 0.01
sig = 0.4
sim = 10 ** 4
discount_rate = np.exp(-r * dt)

""" Least Squares Monte Carlo """


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


def Payoff(strike, paths, simulations):
    if OptionType == 'call':
        po = np.maximum(paths - strike, np.zeros((T + 1, simulations), dtype=np.float64))
    elif OptionType == 'put':
        po = np.maximum(strike - paths, np.zeros((T + 1, simulations), dtype=np.float64))
    else:
        print('Incorrect input')
        os.execl(sys.executable, sys.executable, *sys.argv)
    return po


def ValueVector(payoff, time, GBM, discount):
    value_matrix = np.zeros_like(payoff)
    value_matrix[-1, :] = payoff[-1, :]
    for t in range(time - 1, 0, -1):
        regression = np.polyfit(GBM[t, :], value_matrix[t + 1, :] * discount, 6)
        continuation_value = np.polyval(regression, GBM[t, :])
        value_matrix[t, :] = np.where(payoff[t, :] > continuation_value, payoff[t, :],
                                      value_matrix[t + 1, :] * discount)
        ValueVector = value_matrix[1, :] * discount
    return ValueVector


def Price(ValueVector, simulations):
    return np.sum(ValueVector) / float(simulations)


OptionType = str(input('Call/put:'))
print('Pricing option...')
GBM = GBM(S0, T, sim, r, sig, dt)
payoff = Payoff(K, GBM, sim)
ValueVector = ValueVector(payoff, T, GBM, discount_rate)
price = Price(ValueVector, sim)
print('-'*4**3)
print('Least Squares Monte Carlo Price:', price)

# Graph the regression fit and simulations

# print('Fitting regression...')
# value_matrix = np.zeros_like(payoff)
# value_matrix[T, :] = payoff[T, :]
# prices = GBM[T, :]
# value = value_matrix[T, :]
# regression = np.polyfit(prices, value * discount_rate, 4)
# continuation_value = np.polyval(regression, prices)
# sorted_index = np.argsort(prices)
# prices = prices[sorted_index]
# continuation_value = continuation_value[sorted_index]

# value_matrix[T, :] = np.where(payoff[T, :] > continuation_value, payoff[T, :], value_matrix[T, :] * discount_rate)
# ValueVector = value_matrix[T, :] * discount_rate
# ValueVector = ValueVector[sorted_index]

# plt.figure()
# f, axes = plt.subplots(2, 1)
# axes[0].set_title('American Option')
# axes[0].plot(prices, continuation_value, label='Fitted Polynomial')
# axes[0].plot(prices, ValueVector, label='Value Vector')
# axes[0].set_ylabel('Payoff')
# axes[0].set_xlabel('Asset Price')
# axes[0].legend()
# axes[0].set_xlim(xmax=underlying * 1.06)
# 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()

""" 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):
    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):
    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):
    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())


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)

    """ 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"')

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

import numpy as np
import matplotlib.pyplot as plt
import os
import sys
from QuantLib import *

plt.style.use('seaborn')

# Define global parameters

S0 = 100
K = 90
valuation_date = Date(17, 4, 2017)
expiry_date = Date(17, 4, 2019)
t = (expiry_date - valuation_date) / 365
T = 100
dt = t / T
r = 0.015
sig = 0.4
sim = 10 ** 4
discount_rate = np.exp(-r * dt)

""" Least Squares Monte Carlo """


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


def Payoff(strike, paths, simulations):
    if OptionType == 'call':
        po = np.maximum(paths - strike, np.zeros((T + 1, simulations), dtype=np.float64))
    elif OptionType == 'put':
        po = np.maximum(strike - paths, np.zeros((T + 1, simulations), dtype=np.float64))
    else:
        print('Incorrect input')
        os.execl(sys.executable, sys.executable, *sys.argv)
    return po


def ValueVector(payoff, time, GBM, discount):
    value_matrix = np.zeros_like(payoff)
    value_matrix[-1, :] = payoff[-1, :]
    for t in range(time - 1, 0, -1):
        regression = np.polyfit(GBM[t, :], value_matrix[t + 1, :] * discount, 8)
        continuation_value = np.polyval(regression, GBM[t, :])
        value_matrix[t, :] = np.where(payoff[t, :] > continuation_value, payoff[t, :],
                                      value_matrix[t + 1, :] * discount)
        ValueVector = value_matrix[1, :] * discount
    return ValueVector


def Price(ValueVector, simulations):
    return np.sum(ValueVector) / float(simulations)


OptionType = str(input('Call/put:'))
print('Pricing option...')
GBM = GBM(S0, T, sim, r, sig, dt)
payoff = Payoff(K, GBM, sim)
ValueVector = ValueVector(payoff, T, GBM, discount_rate)
price = Price(ValueVector, sim)
print('Least Squares Monte Carlo Price:', price)


""" QuantLib Pricing """

S0 = SimpleQuote(S0)
if OptionType == 'call':
    put_or_call = Option.Call
elif OptionType == 'put':
    put_or_call = 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, put_or_call, K, process):
    exercise = AmericanExercise(valuation_date, expiry_date)
    payoff = PlainVanillaPayoff(put_or_call, 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, put_or_call, K, process):
    exercise = AmericanExercise(valuation_date, expiry_date)
    payoff = PlainVanillaPayoff(put_or_call, K)
    option = VanillaOption(payoff, exercise)
    engine = BaroneAdesiWhaleyEngine(process)
    option.setPricingEngine(engine)
    return option


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


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


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


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


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

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

import numpy as np
import matplotlib.pyplot as plt
import os
import sys
from QuantLib import *

plt.style.use('seaborn')

# Define global parameters

S0 = 100
K = 90
valuation_date = Date(1, 1, 2018)
expiry_date = Date(1, 1, 2019)
t = (expiry_date - valuation_date) / 365
T = 100
dt = t / T
r = 0.01
sig = 0.4
sim = 10 ** 4
discount_rate = np.exp(-r * dt)

""" Least Squares Monte Carlo """


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


def Payoff(strike, paths, simulations):
    if OptionType == 'call':
        po = np.maximum(paths - strike, np.zeros((T + 1, simulations), dtype=np.float64))
    elif OptionType == 'put':
        po = np.maximum(strike - paths, np.zeros((T + 1, simulations), dtype=np.float64))
    else:
        print('Incorrect input')
        os.execl(sys.executable, sys.executable, *sys.argv)
    return po


def ValueVector(payoff, time, GBM, discount):
    value_matrix = np.zeros_like(payoff)
    value_matrix[-1, :] = payoff[-1, :]
    for t in range(time - 1, 0, -1):
        regression = np.polyfit(GBM[t, :], value_matrix[t + 1, :] * discount, 6)
        continuation_value = np.polyval(regression, GBM[t, :])
        value_matrix[t, :] = np.where(payoff[t, :] > continuation_value, payoff[t, :],
                                      value_matrix[t + 1, :] * discount)
        ValueVector = value_matrix[1, :] * discount
    return ValueVector


def Price(ValueVector, simulations):
    return np.sum(ValueVector) / float(simulations)


OptionType = str(input('Call/put:'))
print('Pricing option...')
GBM = GBM(S0, T, sim, r, sig, dt)
payoff = Payoff(K, GBM, sim)
ValueVector = ValueVector(payoff, T, GBM, discount_rate)
price = Price(ValueVector, sim)
print('-'*4**3)
print('Least Squares Monte Carlo Price:', price)

# Graph the regression fit and simulations

# print('Fitting regression...')
# value_matrix = np.zeros_like(payoff)
# value_matrix[T, :] = payoff[T, :]
# prices = GBM[T, :]
# value = value_matrix[T, :]
# regression = np.polyfit(prices, value * discount_rate, 4)
# continuation_value = np.polyval(regression, prices)
# sorted_index = np.argsort(prices)
# prices = prices[sorted_index]
# continuation_value = continuation_value[sorted_index]

# value_matrix[T, :] = np.where(payoff[T, :] > continuation_value, payoff[T, :], value_matrix[T, :] * discount_rate)
# ValueVector = value_matrix[T, :] * discount_rate
# ValueVector = ValueVector[sorted_index]

# plt.figure()
# f, axes = plt.subplots(2, 1)
# axes[0].set_title('American Option')
# axes[0].plot(prices, continuation_value, label='Fitted Polynomial')
# axes[0].plot(prices, ValueVector, label='Value Vector')
# axes[0].set_ylabel('Payoff')
# axes[0].set_xlabel('Asset Price')
# axes[0].legend()
# axes[0].set_xlim(xmax=underlying * 1.06)
# 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()

""" 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):
    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):
    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):
    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())


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)

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 over price 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!

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!