I have implemented a monte carlo simulation for a plain vanilla European Option and I am trying to compare it to the analytical result obtained from the BS formula. Assuming my monte carlo pricer is correctly implemented, am I supposed to get the very same result with both methods (Monte Carlo and BS analytical formula)?
Fundamentally this is no different from other simulation-based estimation---see this little experiment in R:
R> set.seed(42) R> rowMeans(replicate(200,sapply(1:6, +> FUN=function(x) mean(rnorm(10^x)), simplify=TRUE)))  -2.47827e-02 -9.46800e-03 2.38226e-03 -1.08650e-03 9.41395e-05 1.06759e-05 R>
We are calculating the mean of a $N(0,1)$ vector for sample sizes from $10^1$ to $10^6$. That is then repeated 200 times, and we are calculating the mean of the 200 draws at the different sample sizes.
We find that by and large, the mean gets closer to zero. But even at $10^6$, repeated 200 times, we are still pretty far from 'zero'. That is the way it goes with simulation, and it pays to get a feel for this.
So while you have perfect benchmark with your analytical Black-Scholes result, you will be hard-pressed to get the difference to vanish completely.