asian geometric option valuation-- unable to get monte carlo simulation to converge to analytic value

Question

I'm trying to price asian put options in which the averaging window begins immediately (T=0). currently, I'm trying to match up geometric averaging between my Monte Carlo simulations and my attempt at replicating Kemna Vorst's 1990 analytic solution for put options.

Kemna Vorst:

def kemna_vorst_put_val(s0, k, r, T, sig):

    d_star = 0.5*(r-sig**2/6)*T
    d1 = (log(s0/k) + 0.5*(r + sig**2/6)*T)/(sig*sqrt(T/3))
    d2 = d1 - sig*sqrt(T/3)    
    put_val = exp(-r*T)*k*norm.cdf(-d2) - s0*exp(d_star)*norm.cdf(-d1)
    
    return put_val

Monte Carlo:

def asian_option_mc(s0, k, r, dt, sig, m, n):

 
    sig =  sig/sqrt(3)
    delta_t = dt / m  # length of time interval
    p = []
    for i in range(0, n):
        s = [s0]
        for j in range(0, m):
            s.append(s[-1] * exp((r - 0.5 * sig ** 2) * delta_t + (sig * sqrt(delta_t) * random.gauss(0, 1))))
        
        avg = scipy.gmean(s)         
        p.append(max((k - avg), 0))


    put_value = np.mean(p) * exp(-r * dt)

    return put_value

running kemna_vorst_put_val(100, 95, 0.05, 1, 0.20) gives me 0.5089522108680562

while running

asian_option_mc(100, 95, 0.05, 1, 0.20, 252, 100000) gives me 0.4064389976474143

I've run this several times and get answers that are substantially different from one another (On the order of 20-25%). Any idea what I may be doing wrong? Is there an issue with my Kemna Vorst expression? Or is there a mis-specification in my Monte Carlo simulation? I would expect them to be pretty close if I did this correctly. Thank you!