A simple Monte Carlo in Julia and C++

Yifan_Liu · March 13, 2018, 11:56pm

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)

yuyichao · March 14, 2018, 12:23am

A few things:

passing int as const int& will only make your code slower. (Edit: and the same for double)
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.
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.

Yifan_Liu · March 14, 2018, 12:43am

Ah Ah! I missed the warning message! Thanks!

Yifan_Liu · March 14, 2018, 3:30pm

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

Seif_Shebl · March 14, 2018, 4:48pm

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

Yifan_Liu · March 14, 2018, 9:59pm

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

tshort · March 14, 2018, 10:39pm

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

Yifan_Liu · March 15, 2018, 1:54am

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

tshort · March 15, 2018, 11:24am

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

Yifan_Liu · March 15, 2018, 1:47pm




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

dpsanders · March 15, 2018, 2:11pm

You need

const rng = Xoroshiro128Star()

Yifan_Liu · March 15, 2018, 4:20pm

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

kristoffer.carlsson · March 15, 2018, 4:33pm

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

Yifan_Liu · March 15, 2018, 5:50pm

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

ChrisRackauckas · March 15, 2018, 6:09pm

Use one rng. Don’t keep remaking it.

dpsanders · March 15, 2018, 7:55pm

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