JuliaCon 2018: Asian option pricing using Monte-Carlo simulation using multi-threading in Julia

juliacon2018_presentation.md

Asian option pricing using Monte-Carlo simulation using multi-threading in Julia

JuliaCon 2018

Bogumił Kamiński

http://bogumilkaminski.pl/about

rand is not thread safe

# export JULIA_NUM_THREADS=4

function test0()
    x = zeros(10^6)
    for i in 1:10^6
        x[i] = rand()
    end
    length(unique(x))
end

function test1()
    x = zeros(10^6)
    Threads.@threads for i in 1:10^6
        x[i] = rand()
    end
    length(unique(x))
end

# julia> test0()
# 1000000
# 
# julia> test1()
# 422306

A correct way to use rand with threading

if Base.VERSION > v"0.6.4"
    using Random: MersenneTwister
    using Future: randjump
end

function threadlocal_rngs(rng)
    new_rng = similar(rng)
    Threads.@threads for i in 1:Threads.nthreads()
        new_rng[i] = deepcopy(rng[i])
    end
    new_rng
end

if Base.VERSION > v"0.6.4"
    const RNGS1 = accumulate(randjump, [big(10)^20 for i in 1:Threads.nthreads()],
                             init = MersenneTwister(1))
else
    const RNGS1 = randjump(MersenneTwister(1), Threads.nthreads()+1)[2:end]
end

const RNGS2 = threadlocal_rngs(RNGS1)

Asian option valuation benchmark code

using BenchmarkTools

function v_asian_sample(T, r, K, σ, X₀, m, rng)
    X = X₀
    x̂ = zero(X)
    Δ = T / m
    for i in 1:m
        X *= exp((r-σ^2/2)*Δ + σ*√Δ*randn(rng))
        x̂ += X
    end
    exp(-r*T)*max(x̂/m - K, 0)
end

function v_asian(T, r, K, σ, X₀, m, n, rs)
    res = Threads.Atomic{Float64}(0.0)
    total = Threads.Atomic{Int}(0)
    nt = Threads.nthreads()
    Threads.@threads for i in 1:nt
        while true
            v = v_asian_sample(T, r, K, σ, X₀, m, rs[Threads.threadid()])
            old_total = Threads.atomic_add!(total, 1)
            if old_total < n
                Threads.atomic_add!(res, v)
            else
                break
            end
        end
    end
    res[] / n
end

println("\n", Base.VERSION)
println("Threads: ", Threads.nthreads())

@btime v_asian(1.0, 0.05, 55.0, 0.3, 50.0, 100_000, 1_000, RNGS1)
@btime v_asian(1.0, 0.05, 55.0, 0.3, 50.0, 100_000, 1_000, RNGS2)

Sample output for 4 threads (AWS-EC2 c4.4xlarge, 16 logical cores)

0.6.4
Threads: 4
  779.282 ms (3 allocations: 144 bytes)
  379.049 ms (3 allocations: 144 bytes)

0.7.0-beta2.0
Threads: 4
  382.127 ms (2996 allocations: 46.89 KiB)
  382.385 ms (3000 allocations: 46.97 KiB)

More extensive benchmark results

version	threads	RNGS1	RNGS2
0.6.4	1	1532	1531
0.7.0-beta2.0	1	1524	1523
0.6.4	2	1947	755
0.7.0-beta2.0	2	762	763
0.6.4	3	1672	504
0.7.0-beta2.0	3	509	509
0.6.4	4	779	379
0.7.0-beta2.0	4	382	382
0.6.4	5	677	305
0.7.0-beta2.0	5	307	307
0.6.4	6	873	253
0.7.0-beta2.0	6	255	255
0.6.4	7	564	217
0.7.0-beta2.0	7	219	219
0.6.4	8	525	193
0.7.0-beta2.0	8	192	192
0.6.4	9	433	187
0.7.0-beta2.0	9	189	189
0.6.4	10	495	186
0.7.0-beta2.0	10	184	184
0.6.4	11	384	180
0.7.0-beta2.0	11	181	182
0.6.4	12	388	176
0.7.0-beta2.0	12	179	180
0.6.4	13	408	175
0.7.0-beta2.0	13	175	175
0.6.4	14	390	170
0.7.0-beta2.0	14	171	172
0.6.4	15	392	167
0.7.0-beta2.0	15	168	168
0.6.4	16	387	165
0.7.0-beta2.0	16	166	168

Julia 0.6.4

mutable struct MersenneTwister <: AbstractRNG
    seed::Vector{UInt32}
    state::DSFMT_state
    vals::Vector{Float64}
    idx::Int
end

# julia> sizeof(MersenneTwister)
# 32

# julia> diff(sort!(Int.(pointer_from_objref.(RNGS1))))
# 15-element Array{Int64,1}:
#  48
#  48
#  48
#  48
#  48
#  48
#  48
#  48
#  48
#  48
#  48
#  48
#  48
#  48
#  48
# 
# julia> diff(sort!(Int.(pointer_from_objref.(RNGS2))))
# 15-element Array{Int64,1}:
#  26060992
#    180080
#     16384
#     16384
#     16384
#     16384
#     16384
#     16384
#     16384
#     16384
#     16384
#     16384
#   1196032
#     81920
#     49152

Julia 0.7.0-beta2

mutable struct MersenneTwister <: AbstractRNG
    seed::Vector{UInt32}
    state::DSFMT_state
    vals::Vector{Float64}
    ints::Vector{UInt128}
    idxF::Int
    idxI::Int
end

# julia> sizeof(MersenneTwister)
# 48

# julia> diff(sort!(Int.(pointer_from_objref.(RNGS1))))
# 15-element Array{Int64,1}:
#  128
#   64
#   64
#   64
#   64
#   64
#   64
#   64
#   64
#  128
#   64
#   64
#   64
#   64
#   64
# 
# julia> diff(sort!(Int.(pointer_from_objref.(RNGS2))))
# 15-element Array{Int64,1}:
#  2260288
#    98304
#   114688
#   131072
#   114688
#   114688
#    98304
#    98304
#    81920
#    81920
#   131072
#   868352
#  2883584
#    81920
#  2818048

Appendix (how to run the example)

#!/bin/bash

for i in `seq 1 $(nproc)`
do
    export JULIA_NUM_THREADS=$i
    julia6 asianvaluation.jl
    julia7 asianvaluation.jl
done