Debugging an Euler Monte Carlo Simulation for a European Call
Summary
The document presents a debugging question about simulating a European call under geometric Brownian motion and comparing the Monte Carlo estimate with the Black–Scholes price. The response rewrites the R simulation loop, removing redundant counters and correcting how the stock path is indexed and updated. It then calculates the terminal payoff for each simulated path, averages those payoffs, and discounts the result to the present.
The example illustrates that path construction and array bounds can undermine a simulation even when the intended model and pricing formula are familiar. However, the answer offers only code cleanup and index corrections; it does not systematically diagnose every source of bias or error, analyze convergence, or show results from a corrected run. It should therefore be read as a narrow implementation fix rather than a complete validation of the simulation.
Key ideas
- A Monte Carlo call estimate averages terminal payoffs across simulated price paths.
- Discount the average payoff at the risk-free rate over the option’s maturity.
- Careful path indexing is necessary to use the intended terminal stock value.
- Removing redundant loop updates can make simulation code easier to inspect.
- Code corrections alone do not establish convergence or validate every modeling assumption.
Tags
Full text
# Different Results Monte Carlo and Black-Scholes - where is my mistake?
# 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.
Thank you in advance!
## Answer by Raskolnikov (score 1)
https://quant.stackexchange.com/a/37028
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.Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.