# A simple Monte Carlo in Julia and C++

I moved a small Monte Carlo example from C++ to Julia. I noticed that when my simulation number is 1e9, Julia is 2 times faster than C++.

#include <algorithm>    // Needed for the "max" function
#include <cmath>
#include <iostream>
#include <iomanip> // Needed for the "setprecision" function
#include <random> // Needed to generate random number

using namespace std;

double randn()
{
static random_device rd;
static mt19937 gen(rd());
static normal_distribution<double> dist(0.0, 1.0);

double res = dist(gen);

return res;

}

// Pricing a European vanilla call option with a Monte Carlo method
double monte_carlo_call_price(long num_sims, double S, double K, double r, double v, double T)
{
double S_adjust = S * exp(T*(r-0.5*v*v));
double S_cur = 0.0;
double payoff_sum = 0.0;

for (long i=0; i<num_sims; i++)
{
double gauss_bm = randn();
payoff_sum += max(S_cur - K, 0.0);
}

return (payoff_sum / num_sims) * exp(-r*T);
}

int main()
{

cout << setprecision(12)<<monte_carlo_call_price(100000000000, 100.0, 100.0, 0.05, 0.2, 1) << endl;

}
function monte_carlo_call_price(num_sims, S, K, r, v, T)

S_cur = 0.0
payoff_sum = 0.0

for i= 1:num_sims
gauss_bm = randn()
payoff_sum += max(S_cur - K, 0.0)
end

res = (payoff_sum /(num_sims)) * exp(-r*T)

return  res

end

@time monte_carlo_call_price(1_000_000_000_0, 100.0, 100.0, 0.05, 0.2, 1)

A few things:

1. passing int as const int& will only make your code slower. (Edit: and the same for double)
2. You should probably compare the actual time for each language as a function of iteration count. That can easily show which one is having problem.
3. 1e10 overflows int. The compiler should warn it to you and you should pay attention to all of the warnings. If not, make sure you have all of them enabled.
3 Likes

Ah Ah! I missed the warning message! Thanks!

I tried several different fast random number generator engines in C++, including pcg, sfmt, Xorshift+, etc., but still is slower than Julia.

Note that you can even make Julia much faster by dropping the expensive max function and use just a simple ifelse expression as below. This yields Julia near 10X faster than C++. The difference comes mainly from two things; randn in Julia is extremely fast, and the latest Julia 0.7 contains native Julia implementation of exp which is also very fast.

C++ code:

#include <algorithm>    // Needed for the "max" function
#include <cmath>
#include <iostream>
#include <iomanip> // Needed for the "setprecision" function
#include <random> // Needed to generate random number

using namespace std;
double randn() {
static random_device rd;
static mt19937 gen(rd());
static normal_distribution<double> dist(0.0, 1.0);
double res = dist(gen);
return res;
}

// Pricing a European vanilla call option with a Monte Carlo method
double monte_carlo_call_price(long long num_sims, double S, double K, double r, double v, double T) {
double S_adjust = S * exp(T*(r-0.5*v*v));
double S_cur = 0.0;
double payoff_sum = 0.0;

for (long long i = 0; i < num_sims; i++) {
double gauss_bm = randn();
payoff_sum += (S_cur - K > 0.0) ? S_cur - K : 0.0;
}

return (payoff_sum / num_sims) * exp(-r*T);
}

int main() {
cout << setprecision(12)<<monte_carlo_call_price(1000000000, 100.0, 100.0, 0.05, 0.2, 1) << endl;
}

Julia code:

using Random
function monte_carlo_call_price(num_sims, S, K, r, v, T)

S_cur = 0.0
payoff_sum = 0.0

for i = 1:num_sims
gauss_bm = randn()
S_cur = S_adjust * exp(sqrt(v*v*T) * gauss_bm)
payoff_sum += ifelse(S_cur-K > 0.0, S_cur-K, 0.0)
end

res = payoff_sum/num_sims * exp(-r*T)
end

@time y = monte_carlo_call_price(1000000000, 100.0, 100.0, 0.05, 0.2, 1)
println(y)

The timings:

C++:
10.4500318431

Process returned 0 (0x0)   execution time : 100.472 s
Press any key to continue.

Julia:
10.692408 seconds (5 allocations: 176 bytes)
10.45219256672127
5 Likes

Thanks for the reply. Is it possible for Julia to further improve the performance of random number generator?

You could try a different random-number generator. I know @ChrisRackauckas likes Xoroshiro128Star from the RandomNumbers package.

1 Like

I tried the Xoroshiro128Star engine, but it seemed much slower in this case.

â€śMuch slowerâ€ť sounds suspicious. Can you post yor code?

1 Like

function monte_carlo_call_price(num_sims, S, K, r, v, T)

S_cur = 0.0
payoff_sum = 0.0

for i= 1:num_sims
gauss_bm = randn()
payoff_sum += max(S_cur - K, 0.0)
end

res = (payoff_sum /(num_sims)) * exp(-r*T)

return  res

end

using RandomNumbers.Xorshifts
rng = Xoroshiro128Star()

function monte_carlo_call_price_Xoroshiro128Star(num_sims, S, K, r, v, T)

S_cur = 0.0
payoff_sum = 0.0

for i= 1:num_sims
gauss_bm = randn(rng)
payoff_sum += max(S_cur - K, 0.0)
end

res = (payoff_sum /(num_sims)) * exp(-r*T)

return  res

end

@time monte_carlo_call_price(1_000_000_000, 100.0, 100.0, 0.05, 0.2, 1)
28.402888 seconds (27.19 k allocations: 1.455 MiB)
10.45169813930703

@time monte_carlo_call_price_Xoroshiro128Star(1_000_000_000, 100.0, 100.0, 0.05, 0.2, 1)
184.963270 seconds (10.00 G allocations: 149.012 GiB, 2.29% gc time)
10.450295954616529

You need

const rng = Xoroshiro128Star()

7 Likes

Thanks. How can I make rng a static local variable? like this

using RandomNumbers.Xorshifts

function randn_Xoroshiro128Star()
static rng = Xoroshiro128Star()
res = randn(rng)
return res
end

function monte_carlo_call_price_Xoroshiro128Star(num_sims, S, K, r, v, T)

S_cur = 0.0
payoff_sum = 0.0

for i= 1:num_sims
gauss_bm = randn_Xoroshiro128Star()
payoff_sum += max(S_cur - K, 0.0)
end
res = (payoff_sum /(num_sims)) * exp(-r*T)

return  res

end

You donâ€™t. Either pass it in through an argument, or use a global variable.

2 Likes

I got an error message here.

using RandomNumbers.Xorshifts

function randn_Xoroshiro128Star(rng = Xoroshiro128Star())
res = randn(rng)
return res
end

randn_Xoroshiro128Star() # works

function monte_carlo_call_price_Xoroshiro128Star(num_sims, S, K, r, v, T)

S_cur = 0.0
payoff_sum = 0.0

for i= 1:num_sims
gauss_bm = randn_Xoroshiro128Star()
payoff_sum += max(S_cur - K, 0.0)
end

res = (payoff_sum /(num_sims)) * exp(-r*T)

return  res

end

@time monte_carlo_call_price_Xoroshiro128Star(1_000_000_000, 100.0, 100.0, 0.05, 0.2, 1) #get error message

ERROR: SystemError: opening file /dev/urandom: Too many open files

Use one rng. Donâ€™t keep remaking it.

1 Like

Something like this:

julia> function monte_carlo_call_price(num_sims, S, K, r, v, T, rng=Base.GLOBAL_RNG)

S_cur = 0.0
payoff_sum = 0.0

for i= 1:num_sims
gauss_bm = randn(rng)
payoff_sum += max(S_cur - K, 0.0)
end
return payoff_sum
end
monte_carlo_call_price (generic function with 2 methods)

julia> monte_carlo_call_price(10, 100.0, 100.0, 0.05, 0.2, 1)
101.4072899016111

julia> using RandomNumbers
WARNING: using RandomNumbers.AbstractRNG in module Main conflicts with an existing identifier.

julia> xor = Xorshifts.Xoroshiro128Star()
RandomNumbers.Xorshifts.Xoroshiro128Star(0xa23c0ef6e788a562, 0x6059ddf575c4ac8a)

julia> monte_carlo_call_price(10, 100.0, 100.0, 0.05, 0.2, 1, xor)
107.77703550880462
3 Likes