# America option early exercice boundary via Monte Carlo simulation

I am trying to calculate an american option price via the simulation of the early exercise boundary using the method presented in this document: Monte Carlo Method For pricing a put Option.

I have inmplemented the boundary calculation algorithm (p15) but I only find this kind of results:

Which differ greatly from what is expected.

Can you take a look at my code ? (I have used an OptPayoff mtrix so that I won't have to look at the optimal payoff at each step).

r =0.05
sigma= 0.15
T = 1
K = 100
S0 = 100
N = 1000
m = 100

Normal<-matrix(0,N,m+1)

#simulation of normal(0,1)

for (i in 2:(m+1)){
Normal[,i]= rnorm(N)
}

#Black-Scholes trajectories

BS <-matrix(100,N,m+1)
for (i in 2:(m+1)){
BS[,i]= BS[,i-1] * exp((r-sigma*sigma/2)* T/m + sigma * Normal[,i]*sqrt(T/m))
}

# Instantaneous payoff actualised to 0
Payoff <-matrix(K-S0,N,m+1)
for (i in 2:(m+1)){
Payoff[,i]= (K - BS[,i])*exp ( -i*r*T/m)
}

for (i in 1:N){
for (j in 2:(m+1)){
Payoff[i,j] = max(Payoff[i,j],0)
}
}

#Optimal payoff after j, actualised to 0
OptPayoff <-matrix(0,N,m)

for (i in 1:N){
for (j in 1:m){
variable = Payoff[i,(j+1):m]
OptPayoff[i,j] = max(variable)
}
}

theta = matrix(K,m+1)

#Boundary algorithm

for (j  in (m-1):0 ){
theta[j+1,1]= theta[j+2,1]
thetahat = theta[j+1,1]
Phat = 0
for (i  in 1:N){
if (BS[i,j+1] < theta [j+1,1]) { Phat = Phat + Payoff[i,j+1] /N}
else { Phat = Phat + OptPayoff[i,j+1] / N }
}
for (i  in 1:N ){
if(BS[i,j+1] < theta [j+1,1]) {theta[j+1,1] =BS[i,j+1]
P= 0
for (ii  in 1:N ){
if (BS[ii,j+1] < theta [j+1,1]) { P = P + Payoff[ii,j+1] /N}
else { P =  P + OptPayoff[ii,j+1] / N }
}
if (P>Phat){ Phat = P
thetahat= theta[j+1,1]}
}
}
theta[j+1,1] = thetahat
}

plot(theta,type='l')

-