Monte Carlo simulation function performance

New to Julia
1

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)
2

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

3

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:

  1. 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?
  2. Solve the PDE for the weak form if you only need a density solution.
  3. 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.

4

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.

5

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