black scholes - Monte Carlo option pricing with R

I am trying to implement a vanilla European option pricer with Monte Carlo using R. In the following there is my code for pricing an European plain vanilla call option on non dividend paying stock, under the assumption that the stock follows a GBM.

For teaching reasons, I have used both the analytic formula and the Euler-Maruyama approximation.

However, comparing the obtained results with those of B&S model, I found a quite big difference, therefore I would like to ask you if you can spot the mistake in my Monte Carlo code:

# Compute the Black-Scholes European option price on non-dividend paying stock
# Setting the  B&S parameters value
S <- 52 #stock price at time t
K <- 50 #strike price 
tau <- 0.25 #time to maturity T - t (in years) #0.25 = 3 months
r <- 0.05 #risk-free annual interest rate

sigma <- 0.3 #annual volatility of the stock price (standard deviation)

#call B&S fair value
d1 <- (log(S/K) + (r + 0.5*sigma^2)*tau)/(sigma*sqrt(tau))
d2 <- d1 - sigma*sqrt(tau)

V_BS_Call <- S*pnorm(d1) - K*exp(-r*(tau))*pnorm(d2) #fair value call


# Compute the Monte Carlo European option price on non-dividend paying stock 

# Assuming the non- dividend paying stock follows a Geometric Brownian Motion (GBM)

set.seed(2503) #set the seed
# Setting the Monte Carlo simulation and GBM  parameters
tau <- tau #time to expiry (we have already defined this variable)
N <- 250 #number of sub intervals
dt <- tau/N #length of each time sub interval
time <- seq(from=0, to=tau, by=dt) #time moments in which we simulate the process
length(time) #it should be N+1
nSim <- 10000 #number of simulations (paths) 


r <- r #GBM parameter 1
sigma <- sigma #GBM parameter 2
X0 <- S #initial condition (price of the underlying today)

#Monte Carlo with analytic formula
Z <-  matrix(rnorm(nSim*N, mean=0, sd=1),nrow = nSim, ncol = N) #standard normal sample of N elements
dW <- Z*sqrt(dt) #Brownian motion increments (N increments)x nSim simulations
W <- matrix(numeric(nSim*(N+1)), nrow = nSim, ncol = (N+1))
X_analytic <- numeric(nSim)

for(k in 1:nSim){
  W[k,] <- c(0, cumsum(dW[k,]))
  X_analytic[k] <- X0*exp((r - 0.5*sigma^2)*tau + sigma*W[k,ncol(W)]) #Analytic solution
}
payoff_expiry_call <-pmax(X_analytic-K,0) #pmax preserve the dimension of the matrix, so apply the max function to each element
expected_payoff_call <- sum(payoff_expiry_call)/length(payoff_expiry_call)
Monte_Carlo_call_price <- exp(-r*(tau))*expected_payoff_call

#Monte Carlo with Euler-Maruyama scheme
X_EM <- matrix(numeric(nSim*(N+1)), nrow = nSim, ncol = (N+1))

X_EM[,1] <- X0 #first element of X_EM is X0. with the for loop we find the other N elements

for(k in 1:nSim){
  for(i in 2:ncol(X_EM)){
    X_EM[k,i] <- X_EM[k,i-1] + r*X_EM[k,i-1]*dt + sigma*X_EM[k,i-1]*dW[k,i-1]
  }
}

payoff_expiry_call <-pmax(X_EM[,ncol(X_EM)]-K,0) #pmax preserve the dimension of the matrix, so apply the max function to each element
expected_payoff_call <- sum(payoff_expiry_call)/length(payoff_expiry_call)

Monte_Carlo_call_price <- exp(-r*(tau))*expected_payoff_call

So, using a 10,000 simulations:

• 
  

  the Monte Carlo price with analytical formula is approximately 4.535

  
  


• 
  

  the Monte Carlo price using Euler-Maruyama is approximately 4.536

  
  



• 
  

  the B&S price is 4.519

  
  

I think that the difference is too big, but I cannot spot the mistake.

Answer

Your code looks fine and it is encouraging that both MC simulations yield similar results. Please look at this simplified code for the analytical part of the Monte Carlo simulation. As you know, $$S_T=S_0\exp\left(\left(r-\frac{1}{2}\sigma^2\right)T+\sigma W_T\right).$$ A call is path-independent, so there is no need to simulate the entire path. I guess you want to teach your students to code as efficiently as possible. Since $W_T\sim N(0,T)$, you can directly simulate the final Brownian motion.

Z <- rnorm(nSim, mean=0, sd=1)
WT <- sqrt(tau) * Z
ST = X0*exp((r - 0.5*sigma^2)*tau + sigma*WT)

simulated_call_payoffs <- exp(-r*tau)*pmax(ST-K,0)
Call_price_MC_anal <- mean(simulated_call_payoffs)

If you play with this a bit around, you indeed obtain various prices which are not too close to the Black Scholes closed form solution. 10,000 sample values are just too few to accurately estimate the option price. Try one million simulations instead.

You could, in general, use this a motivation why variance reduction is so crucial for Monte Carlo simulations. The estimate may be consistent and unbiased but that does not help you if you have large standard errors. Recall that the confidence interval for the MC estimator is given by $$ \hat{C}_n \pm z_{\delta/2}\frac{s_C}{\sqrt{n}},$$ where $\hat{C}_n$ is the estimated call price with $n$ simulations and $s_c$ is the sample variance of the simulated call values. Obviously, the larger $n$, the smaller this interval. If nSim=1000000, I get an interval of $[4.51,4.53]$ (the BS price is $4.52$) but nSim=10000 only gives $[4.45, 4.69]$. The 95% confidence interval is computed via

lower_bound <- Call_price_MC_anal - 1.96*sd(simulated_call_payoffs)/sqrt(nSim)
upper_bound <- Call_price_MC_anal + 1.96*sd(simulated_call_payoffs)/sqrt(nSim)