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

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)