# Different Results Monte Carlo and Black-Scholes - where is my mistake?

as an exercise, I am trying to simulate the BS model via Monte Carlo Simulation in R to price a normal European-style call option. However, the code will give me results that are way higher than the BS results even in case of 100,000 simulations. Here is my Code:

nSim=1000
T=5
N=T*360
K=7589.42
r=0.0219
Y0=7846.33
dt=T/N
drift=0.0
sigma=0.197
a=1
b=1

callOptionPrice=simulateCall()
callOptionPrice

simulateCall=function(){
C<-vector(mode="double",nSim)
for(a in 1:nSim){
V=newPath()
C[a]=max(V[N]-K,0)
a=a+1
}
call=1/nSim*sum(C)*exp(-r*T)
return(call)
}

newPath=function() {
Y<-vector(mode="double", length=N)
i=0
t=1
Y[1]=Y0
for(i in 0:N){
dW=sqrt(dt)*rnorm(1)
Y[t+1] = Y[t] + r*Y[t]*dt + sigma*Y[t]*dW
t=t+1
}
return(Y)
}


Could anybody help me to find the mistake here? The actual result according to BS should be 1864.1388, however my code always returns numbers >2400.

• I tried your code several times, it gave results between 1800 and 1900 as you seem to expect. As @LocalVolatility said just below, you can directly compute the final value of $S_T$ without time discretization if you don't need to price path-dependant option thereafter. You should also consider vectorization in your code, it will be way faster! – AlexM Dec 22 '17 at 20:25
• I'm voting to close this question as off-topic because this is a programming question and should go to stackoverflow. – rbm Feb 20 '18 at 12:36

I've cleaned up your code a bit:

nSim=1000
T=5
N=T*360
K=7589.42
r=0.0219
Y0=7846.33
dt=T/N
drift=0.0
sigma=0.197

simulateCall=function(){
C<-vector(mode="double",nSim)
for(a in 1:nSim){
V=newPath()
C[a]=max(V[N]-K,0)
}
call=1/nSim*sum(C)*exp(-r*T)
return(call)
}

newPath=function() {
Y<-vector(mode="double", length=N)
Y[1]=Y0
for(i in 1:(N-1)){
dW=sqrt(dt)*rnorm(1)
Y[i+1] = Y[i] + r*Y[i]*dt + sigma*Y[i]*dW
}
return(Y)
}

callOptionPrice=simulateCall()
callOptionPrice


Check where I made changes, some stuff you put in there was useless like specifying the variable you loop over and adding code within the loop to increment that variable. Or defining superfluous variables.

• If you are only interested in pricing European options (as it is the case here), there is no need to make the grid finer if you directly simulate the solution of the GBM instead of an Euler discretization of the spot dynamics. Using more than one time step does not improve convergence in this case. – LocalVolatility Nov 22 '17 at 9:02
• Not sure I understand your comment, he's using an Euler discretization, right? – Raskolnikov Nov 22 '17 at 9:04
• Yes, he is. You originally suggested to use even more time steps to improve convergence. My point is that instead of doing this it would be better to just directly simulate the solution for $S_T$ instead and only use one time step. – LocalVolatility Nov 22 '17 at 9:05
• @LocalVolatility. OK, I understand. Maybe he wants to adapt the algorithm to other options though, in which case he might need to simulate the entire path. – Raskolnikov Nov 22 '17 at 9:06
• Now it shows the proper solution. The mistake seems to be that I used N instead of N-1 iterations. I am simulating the whole paths because I want to be able to visualize some sample paths for understanding of the dynamics. Also, I will add compound Poisson jumps with time-inhomogeneous frequency to the model and want to be able to see the difference to pure GBM diffusion dynamics. Finally, I in the next step, I might consider valuation of American options. If I understood it correctly, I will need to simulate full sample paths due to path-dependence of American options. – MrPefister Nov 22 '17 at 10:05