Pricing American options by Longstaff and Schwartz method

[mrmoellerjacob]

I am using monte carlo process to price american options, based on black-scholes framework. i expect the american call option prices equal to european prices when there is no dividend and larger than european call prices otherwise.
The simulated price i got for the american option is lower than the european. Can i anybody help me see where is the problem?

My program looks like this, and I'm using the program R

n <- 3 #antal af exercise points
N <- 10000 #numbers of simulated paths%
K <- 30 ##strike price
r <- 0.05 #interest rate
T <- 1 #time
sigma <- 0.4 #volatility
S0 <-36 #asset price


randn <- matrix(rnorm(n*N),N, n)

dt <- T/n
disc <- exp(-r*dt)

dX <- sqrt(dt)*randn

#Simulated paths
S <- matrix(ncol=n, nrow=N)
for (i in 2:n){
for (j in 1:N){
S[j,1] <- S0*exp((r-1/2*sigma^2)*dt+sigma*dX[j,1])
S[j,i] <- S[j,i-1]*exp((r-1/2*sigma^2)*dt+sigma*dX[j,i])
}
}


#Cashflow-matrix#
f <- function(x){max(K-x,0)}
P <- matrix(data=NA, N, n)
for (j in n:1) {
for (i in 1:N) {
P[i,j] <- f(S[i,j])
}
}

#Basis functions#
f1 <- function(x){exp(-x/2)}
f2 <- function(x){exp(-x/2)*(1-x)}
f3 <- function(x){exp(-x/2)*(1-2*x+x^2/2)}

X <- matrix(data=NA, N, 1)
Y <- matrix(data=NA, N, 1)
for (i in (n-1):1){
for (j in 1:N){
if (P[j,i] > 0) {
X[j,1] <- S[j,i]
Y[j,1] <- disc*P[j,i+1]
}
else{
Y[j,1] <- 0
X[j,1] <- 0
}
}

B1 <- ifelse(f1(X)[,1]==1,0,f1(X)[,1])
B2 <- ifelse(f2(X)[,1]==1,0,f2(X)[,1])
B3 <- ifelse(f3(X)[,1]==1,0,f3(X)[,1])
R <- cbind(Y,B1,B2,B3)
index <- which(P[,i]==0)
reg <- data.frame(R[-index,])

Beta <- lm(V1 ~ B1 + B2+ B3, data=reg)$coefficients
C <- ifelse(X==0,0, rep(Beta[1],N)+Beta[2]*R[,2] + Beta[3]*R[,3] + Beta[4]*R[,4])

for (j in 1:N){
if (C[j] < P[j,i]){
P[j,(i+1):n] <- 0
}
else{
P[j,i] <- 0
}
}
}


Z <- matrix(data=NA, N, n)
g <- function(x){x*disc^j}
for (i in 1:N){
for (j in 1:n){
Z[i,j] <- g(P[i,j])
}
}

AMPrice <- sum(Z)/N
AMPrice

 

[Jasmeet Singh Gujral]

Hi,
I am coding longstaff schwartz in python and facing the same issue. Can anyone please help.

def LSM_model(self):
_s_stimulated = self.stimulation_method()
_intrinsic_val = self.option_payoff(self.stimulation_method())
_option_value = np.zeros_like(_intrinsic_val)
_option_value[:, -1] = _intrinsic_val[:, -1]

for t in range(self.no_of_steps - 1, 0, -1):
_option_value[:, t] = _option_value[:, t + 1] * self.step_disc_fact
if len(_s_stimulated[_intrinsic_val[:, t] > 0, t]) == 0:
continue
else:

_ITM_check = _intrinsic_val[:, t] > 0
_x_axis = _s_stimulated[_ITM_check, t]
_y_axis = _option_value[_ITM_check, t + 1] * self.step_disc_fact
_rg = np.polyfit(_x_axis, _y_axis, 5)
_pv_rg = np.polyval(_rg, _x_axis)

_option_value[_ITM_check, t] = np.where(_intrinsic_val[_ITM_check, t] > _pv_rg,
_intrinsic_val[_ITM_check, t], _y_axis)

_option_value[:, 0] = _option_value[:, 1] * self.step_disc_fact
return np.mean(_option_value[:, 0])

 

[Daniel Duffy]

Is there a bug in your code? or wrong model (MC for early exercise??)

 

[Jasmeet Singh Gujral]

Hi Daniel,

Thanks for the reply. I am stuck here. Trying to follow the logic of early exercise but the code is throwing American option values lower than European, which definitely is wrong.

I have the complete code on git in the form of python package.

Will really appreciate if you could guide me on the possible source of error. What I am trying to replicate is Monte Carlo code here <quantsbin/pyfin>

Will really appreciate help on this.

 

[Daniel Duffy]

I have the complete code on git in the form of python package.

Your code or someone else's?

 

[Jasmeet Singh Gujral]

I have the complete code on git in the form of python package.

Your code or someone else's?

Git code under Quantsbin is mine. To start off with I am trying to replicate what Pyfin has done. But for Longstaff Schwartz, I am missing something in the code which is not giving me desired result.

 

[Daniel Duffy]

Git code under Quantsbin is mine. To start off with I am trying to replicate what Pyfin has done. But for Longstaff Schwartz, I am missing something in the code which is not giving me desired result.

Ok, Then you need to find that missing something. That's life.

Just posting raw, undocumented code is a waste of time and energy.

// Having said that, where's the orthogonal polynomial stuff?