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:
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.