A simple Monte Carlo in Julia and C++

question

#1

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();
    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;


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

#2

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

Ah Ah! I missed the warning message! Thanks!


#4

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


#5

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();
    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:

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:

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

#6

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


#7

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


#8

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


#9

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


#10



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 



@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

#11

You need

const rng = Xoroshiro128Star()


#12

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_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

#13

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


#14

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_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


#15

Use one rng. Don’t keep remaking it.


#16

Something like this:


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