• C++ Programming for Financial Engineering
    Highly recommended by thousands of MFE students. Covers essential C++ topics with applications to financial engineering. Learn more Join!
    Python for Finance with Intro to Data Science
    Gain practical understanding of Python to read, understand, and write professional Python code for your first day on the job. Learn more Join!
    An Intuition-Based Options Primer for FE
    Ideal for entry level positions interviews and graduate studies, specializing in options trading arbitrage and options valuation models. Learn more Join!

Pricing American options by Longstaff and Schwartz method

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

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]
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
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
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:

_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])
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

C++ author, trainer
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?
Last edited: