title: Why is Matlab 11 times faster than C++
body: <p>I am comparing the speed of a Monte Carlo pricing algorithm for a Vanilla call option between Matlab and C++. This is not the same as <a href="https://stackoverflow.com/questions/6058139">Why is MATLAB so fast in matrix multiplication?</a> since the speed-up is not due to matrix multiplication (there is only a dot product which is done quickly) but seems to be due to its highly efficient Gaussian random number generator. </p>
  
  <p>In Matlab the code has been vectorised and the code looks as follows</p>
  
  <pre><code>function [ value ] = OptionMCValue( yearsToExpiry, spot, strike, riskFreeRate, dividendYield, volatility, numPaths  )
  
      sd = volatility*sqrt(yearsToExpiry);
      sAdjusted = spot * exp( (riskFreeRate - dividendYield - 0.5*volatility*volatility) * yearsToExpiry);
  
      g = randn(1,numPaths);
      sT = sAdjusted * exp( g * sd );
      values = max(sT-strike,0);`
      value = mean(values);
      value = value * exp(-riskFreeRate * yearsToExpiry);
  
  end
  </code></pre>
  
  <p>If I run this with 10 million paths as follows</p>
  
  <pre><code>strike = 100.0;
  yearsToExpiry = 2.16563;
  spot = 100.0;
  volatility = 0.20;
  dividendYield = 0.03;
  riskFreeRate = 0.05;
  oneMillion = 1000000;
  numPaths = 10*oneMillion;
  
  tic
  value = OptionMCValue( yearsToExpiry, spot, strike, riskFreeRate, dividendYield, volatility, numPaths  );
  toc
  </code></pre>
  
  <p>I get </p>
  
  <pre><code>Elapsed time is 0.359304 seconds.
     12.8311
  </code></pre>
  
  <p>Now I do the same thing in C++ in VS2013</p>
  
  <p>My code is in an OptionMC class and is as follows</p>
  
  <pre><code>double OptionMC::value(double yearsToExpiry, 
                     double spot,
                     double riskFreeRate,
                     double dividendYield,
                     double volatility, 
                     unsigned long numPaths )
  {
      double sd = volatility*sqrt(yearsToExpiry);
      double sAdjusted = spot * exp( (riskFreeRate - dividendYield - 0.5*volatility*volatility) * yearsToExpiry);
      double value = 0.0;
      double g, sT;
  
      for (unsigned long i = 0; i &lt; numPaths; i++)
      {
          g = GaussianRVByBoxMuller();
          sT = sAdjusted * exp(g * sd);
          value += Max(sT - m_strike, 0.0);
      }
  
      value = value * exp(-riskFreeRate * yearsToExpiry);
      value /= (double) numPaths;
      return value;
  }
  </code></pre>
  
  <p>The BM code is as follows</p>
  
  <pre><code>double GaussianRVByBoxMuller()
  {
  double result;
  double x; double y;;
  double w;
  
  do
  {
      x = 2.0*rand() / static_cast&lt;double&gt;(RAND_MAX)-1;
      y = 2.0*rand() / static_cast&lt;double&gt;(RAND_MAX)-1;
      w = x*x + y*y;
  } while (w &gt;= 1.0);
  
  w = sqrt(-2.0 * log(w) / w);
  result = x*w;
  
  return result;
  }
  </code></pre>
  
  <p>I have set the Optimization option to Optimize for speed in Visual Studio.</p>
  
  <p>For 10m paths it takes 4.124 seconds.</p>
  
  <p>This is 11 times slower than Matlab.</p>
  
  <p>Can anyone explain the difference ?</p>
  
  <p>EDIT: On further testing the slow down does seem to be the call to GaussianRVByBoxMuller. Matlab seems to have a very efficient implementation - the Ziggurat method. Note that BM is sub-optimal here as it generates 2 RVs and I only use 1. Fixing just this would give a factor of 2 speed-up.</p>