Faster geometric brownian motion

gbm.jl

function genS_jl(I)
    s0 = 600.0
    r = 0.02
    sigma = 2.0
    T = 1.0
    M = 100
    dt = T/M
    a = (r - 0.5*sigma^2)*dt
    b = sigma*sqrt(dt)

    paths = zeros(Float64, M, I)

    for i in 1:I
        paths[1, i] = st = s0
        for j in 2:M
            st *= exp(a + b*randn())
            paths[j, i] = st
        end
    end
    return paths
end

genS_jl(10) # Warm up JIT

@elapsed genS_jl(100000) # Outputs 0.538962298

Thank you! I'm still new to julia so this helps a lot.

Also note that if you just want to calculate the option price, you don't need to allocate any arrays at all because the option price only depends on the final value in each trajectory, and you can compute the average incrementally. You could write blazing fast straight scalar code that probably keeps everything in registers throughout.