The problem is taken from the book Foundations and Methods of Stochastic Simulation by Barry Nelson where it is solved using VBA in section 4.5.

I will not repeat the whole discussion of the model, but focus on the numerical aspect. We are asked to calculate the value of the expression:

*X*is generated by a continuous state geometric Brownian motion. We will use Monte Carlo simulation to approximate the above expected value.

A single pass of the Monte Carlo simulation is approximated by a discrete sum:

The parameters we will use in the simulation are:

*T*=1,

*r*=0.05,

*K*=55, σ=0.3,

*m*=1000 and starting value of

*X*at time 0 is 50. We will run 100,000 replications of Monte Carlo simulation and calculate the average.

I want to compare five implementations of the above problem:

- Julia using loops
- Julia using vectorized code
- Python using loops
- Numba using loops
- NumPy using vectorized code

First go the codes in Julia (running time is given in the comment at the end):

And here are the codes in Python:

Julia implementation is a bit more fancy as it promotes the arguments to a common type on the run and in Python I pass all values as floats (which is simpler). Also vectorized code in Julia uses in-place updating. In general I tried to make the implementations a normal code in respective languages.

The key take-aways are the following:

The key take-aways are the following:

- standard Python is a no-go solution for this problem;
- loops in Julia are fastest;
- somewhat surprisingly vectorized Julia code is faster than Numba although the former has to allocate more memory;
- NumPy implementation is around three times slower than vectorized Julia;
- Vectorized Julia code hugely benefits from in-place operations (that avoid memory allocation); however, even without these optimizations it was faster than Numba.