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!

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