# A simple Monte Carlo in Julia and C++

**URL:** <https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [March 13, 2018, 11:56pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702 "2018-03-13T23:56:03Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [March 13, 2018, 11:56pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/1 "2018-03-13T23:56:03Z")

</div>

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

```julia
#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();
    S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm);
    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;

}

```

```julia
function monte_carlo_call_price(num_sims, S, K, r, v, T) 
    
    
        S_adjust = S * exp(T*(r-0.5*v*v))
        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 += 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)

```

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [March 14, 2018, 12:23am UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/2 "2018-03-14T00:23:06Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [March 14, 2018, 12:43am UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/3 "2018-03-14T00:43:18Z")

</div>

Ah Ah! I missed the warning message! Thanks!

---

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [March 14, 2018, 3:30pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/4 "2018-03-14T15:30:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [March 14, 2018, 4:48pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/5 "2018-03-14T16:48:18Z")

</div>

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:

```julia
#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();
    S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm);
    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:

```julia
using Random 
function monte_carlo_call_price(num_sims, S, K, r, v, T)     
    
	S_adjust = S * exp(T*(r-0.5*v*v))
	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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [March 14, 2018, 9:59pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/6 "2018-03-14T21:59:45Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![tshort](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tshort/32/43_2.png) [@tshort](https://discourse.julialang.org/u/tshort)\
**Post date:** [March 14, 2018, 10:39pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/7 "2018-03-14T22:39:33Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [March 15, 2018, 1:54am UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/8 "2018-03-15T01:54:36Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![tshort](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tshort/32/43_2.png) [@tshort](https://discourse.julialang.org/u/tshort)\
**Post date:** [March 15, 2018, 11:24am UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/9 "2018-03-15T11:24:24Z")

</div>

“Much slower” sounds suspicious. Can you post yor code?

---

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [March 15, 2018, 1:47pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/10 "2018-03-15T13:47:59Z")

</div>

```julia

function monte_carlo_call_price(num_sims, S, K, r, v, T) 
    
    
    S_adjust = S * exp(T*(r-0.5*v*v))
    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 += 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_adjust = S * exp(T*(r-0.5*v*v))
    S_cur = 0.0
    payoff_sum = 0.0
  
    for i= 1:num_sims
      gauss_bm = randn(rng)
      S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm)
      payoff_sum += max(S_cur - K, 0.0)
    end
    
    res = (payoff_sum /(num_sims)) * exp(-r*T) 
  
return res

end 

```

```julia
@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

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 15, 2018, 2:11pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/11 "2018-03-15T14:11:35Z")

</div>

You need

`const rng = Xoroshiro128Star()`

---

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [March 15, 2018, 4:20pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/12 "2018-03-15T16:20:20Z")

</div>

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

```julia
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_adjust = S * exp(T*(r-0.5*v*v))
    S_cur = 0.0
    payoff_sum = 0.0
  
    for i= 1:num_sims
      gauss_bm = randn_Xoroshiro128Star()
      S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm)
      payoff_sum += max(S_cur - K, 0.0)
    end    
    res = (payoff_sum /(num_sims)) * exp(-r*T) 
  
return res

end

```

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [March 15, 2018, 4:33pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/13 "2018-03-15T16:33:31Z")

</div>

> [@Yifan\_Liu](#):
>
> Thanks. How can I make rng a static local variable? like this

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

---

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [March 15, 2018, 5:50pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/14 "2018-03-15T17:50:08Z")

</div>

I got an error message here.

```julia
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_adjust = S * exp(T*(r-0.5*v*v))
    S_cur = 0.0
    payoff_sum = 0.0
  
    for i= 1:num_sims
      gauss_bm = randn_Xoroshiro128Star()
      S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm)
      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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 15, 2018, 6:09pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/15 "2018-03-15T18:09:17Z")

</div>

Use one rng. Don’t keep remaking it.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 15, 2018, 7:55pm UTC](https://discourse.julialang.org/t/a-simple-monte-carlo-in-julia-and-c/9702/16 "2018-03-15T19:55:31Z")

</div>

> [@Seif\_Shebl](#):
>
> monte\_carlo\_call\_price(1000000000, 100.0, 100.0, 0.05, 0.2, 1)

Something like this:

```julia

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

           S_adjust = S * exp(T*(r-0.5*v*v))
           S_cur = 0.0
           payoff_sum = 0.0

           for i= 1:num_sims
             gauss_bm = randn(rng)
             S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm)
             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

```