American Option Pricing via Longstaff and Schwartz (2001)

import numpy as np
import matplotlib.pyplot as plt
import numpy.random as npr

def gen_sn(M, I,  anti_paths=True, mo_match=True):
if anti_paths is True:
sn = npr.standard_normal((M+1, I/2))
sn = np.concatenate((sn,-sn), axis=1)
else:
sn = npr.standard_normal((M+1, I))

if mo_match is True:
sn = (sn-sn.mean())/sn.std()

return sn

def gbm_mcs_amer(K, option='C'):
# set model paramters
S0 = 100    # initial stock level
T = 1.0     # time to maturity
r = 0.05    # risk free rate
vol = 0.20  # volatility
M = 1000    # time stpes

dt = T/M
df = np.exp(-r*dt)

# simulation of stock levels
S = np.zeros((M+1, I))
S[0] = S0
rn = gen_sn(M, I)

for t in range(1, M+1):
S[t] = S[t-1]*np.exp((r-0.5*vol**2)*dt+vol*np.sqrt(dt)*rn[t])

if option == 'call':
h = np.maximum(S-K, 0)
else:
h = np.maximum(K-S, 0)

# LSM method
V = np.copy(h)
for t in range(M-1, 0, -1):
# calculate holding value: present value of next period payoff
# regression to get expected holding value
reg = np.polyfit(S[t], V[t+1]*df, 3)
C = np.polyval(reg, S[t])
# compare with exercise value
V[t] = np.where(C>h[t], V[t+1]*df, h[t])

# option value at time 0
C0 = df*np.sum(V[1])/I

return C0

print gbm_mcs_amer(100)