I have a function as below, I want to know how I can improve its performance

```
using Distributions
function OptionMC(S0, K, T, flag = "call")
sigma = 0.2
r = 0.05
NbTraj = 50_000
NbPas = 50_000
DeltaT = T/NbPas
SPresent = S0*ones(NbTraj)
SNext = zeros(NbTraj)
dist = Normal(0, 1)
dw = zeros(NbTraj)
for i = 1:NbPas
for j = 1:NbTraj
dw[j] = sqrt(DeltaT)*rand(dist)
SNext[j] = SPresent[j] + r*SPresent[j]*DeltaT + sigma*SPresent[j] * dw[j]
end
SPresent = SNext
end
if flag == "call"
return exp(-r*T)*mean(max.(0, SPresent - K))
else
return exp(-r*T)*mean(max.(0, K - SPresent))
end
end
@time OptionMC(100, 100, 0.5)
```

Start with the performance tips, especially profiling, then ask specific questions here.

(You may want to use less whitespace, it would make your code more readable.)

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.

I am using it as an example to show that how I can reduce the difference between monte carlo estimated price and analytic price with different techniques, and this is the starting point.

Yeah then the next thing to do is solve some PDEs .