There’s small things that you can do. You don’t need `SPresent`

and `SNext`

, instead just index differently. You can also generate a vector of random numbers with `rand!`

at a time and that might be helpful. You can place `@inbounds`

and manually get rid of some operations, like cache `r*DeltaT`

. Then you can parallelize it. Multithread it or make a GPU version.

Small performance optimizations only get so far. Basically, the Euler-Maruyama algorithm is a really awful algorithm. There’s very few cases where it’s a good choice and it’s usually just used because it’s easy. Things you should be doing:

- Geometric brownian motion has a closed-form solution for its probability density. The Black-Scholes problem has a closed-form solution. Why not use that?
- Solve the PDE for the weak form if you only need a density solution.
- The transformation to an additive noise problem is well-known for this (it’s part of the proof of the closed-form solution). On the transformed equation you can apply the algorithm from here (with adaptivity if you need it, depends on the parameters) and increase the time steps. If the parameters are wonky enough then this can decrease the overall work to get the same (strong and weak) accuracy (BS is simple enough that higher order algorithms may not pay off if you only need a digit or two of accuracy). (For reference the method is mentioned here).

Most of handling an SDE correctly is about utilizing the right algorithms, not trying to make a loop faster. So while you can probably snag like a 2x by trying a few things, you’re looking at like 1000x+ speedups if you instead just sample from a lognormal distribution. And if your actual problem doesn’t have a closed-form solution, you’re likely going to want high order adaptive methods and not Euler-Maruyama.