How To Make Monte Carlo Simulation Run Faster

Question

I'm running a Monte Carlo simulation of stocks prices in which the data that is used to define the Skewed Normal Distribution updates after 1,000 simulations. Thus it is crucial that I get this code to run as fast as possible

listProduct[x_List] := Times @@ x
ParallelTable[listProduct[(1 + 
 Table[
  RandomVariate[
   SkewNormalDistribution[Mean[data2], StandardDeviation[data2], 
    Skewness[data2]]], {i, 1, 5}])] - 1, {i, 1, 1000}]

data2 is just a list of daily returns. I tried compiling it but that didn't alter the performance at all.

cf = Compile[{{x1, _Real}, {x2, _Real}, {x3, _Real}},RandomVariate[SkewNormalDistribution[x1, x2, x3]], CompilationTarget -> "C"]

AbsoluteTiming[
ParallelTable[
listProduct[(1 + 
   Table[
    cf[Mean[data2], StandardDeviation[data2], 
     Skewness[data2]], {i, 1, 5}])] - 1, {i, 1, 1000}];]

Does anyone know I make this run as fast as possible. Are there any internal functions I should look at.

Thanks

Answer 1

I think there are some clear places for performance improvements here:

I assumed a large data2 value using the following:

data2 = RandomReal[{5, 10}, 1000000];

Using the original expression (except for 2000 iterations instead of 1000, because that's what I did the calculation for after checking how the solution scales):

AbsoluteTiming[
 res = ParallelTable[
    listProduct[(1 + 
        Table[RandomVariate[
          SkewNormalDistribution[Mean[data2], 
           StandardDeviation[data2], Skewness[data2]]], {i, 1, 
          5}])] - 1, {i, 1, 2000}];]

{31.7116, Null}

Almost all of this time appears to be spent recalculating the statistics for the SkewNormalDistribution.

Let's pre-calculate that:

dist = SkewNormalDistribution[Mean[data2], StandardDeviation[data2], Skewness[data2]];

And run it again:

AbsoluteTiming[
 res = ParallelTable[
    listProduct[(1 + Table[RandomVariate[dist], {i, 1, 5}])] - 1, {i, 
     1, 2000}];]

{0.068085, Null}

This difference becomes even more notable when the number of iterations is large or when data2 is very large.

There is still a small gain which can be had with RandomVariate. You don't need to use Table to get multiple samples from it, you can just directly ask it to generate multiple.

AbsoluteTiming[
 res = ParallelTable[
    listProduct[(1 + RandomVariate[dist, 5])] - 1, {i, 1, 2000}];]

{0.034347, Null}

Note also that parallelization does not actually represent much of an improvement:

AbsoluteTiming[
 res = Table[
    listProduct[(1 + RandomVariate[dist, 5])] - 1, {i, 1, 2000}];]

{0.052233, Null}

This result is running on 1 kernel versus 16 parallel kernels, but the speedup factor is less than 2.

We can make things even faster by decreasing the number of times we call RandomVariate. One potential way to do this is:

listProduct2[l_] := Times @@ (1 + l) - 1;
AbsoluteTiming[res = listProduct2 /@ RandomVariate[dist, {2000, 5}];]

{0.007448, Null}

This creates the 2000 separate variates at once and then maps the effect of the original listProduct expression onto each of them at the same time. This appears to be about 5-10x faster again.

These all have somewhat different but very similarly looking results, and as such I am reasonably sure they're all the same process you were looking to simulate.

Answer 2

Assuming data2 is just a static list of numbers, which I'll call d, you want to avoid re-calculating the same thing multiple times, specifically the sample mean, standard deviation and skew, and the distribution itself.

For d I'll assume a range of {-100,100} and a length of 10000. It looks like you're using Table to call RandomVariate five times, then adding one and multiplying the resulting list together. You're probably better off letting MMA generate all of your numbers at once using the second argument of RandomVariate. You can generate a thousand lists of length five using RandomVariate[...,{1000,5}] and then do the rest of your operations on the resulting set of lists. You can use @@@ to map Times@@ to each of the lists generated without defining an additional function.

d = RandomReal[{-100,100},10000]; (*Dummy data needed for MWE*)
moments = {Mean@#,StandardDeviation@#,Skewness@#}&@d; (*define moments*)
sn = SkewNormalDistribution[Sequence@@moments];(*Define distribution*)
(Times @@@ (1 + RandomVariate[sn, {1000, 5}])) - 1

This evaluates in about a 1ms on my laptop, about ten times faster than your original method which seems to be dominated by the time it takes to launch the parallel kernels. Since I've only calculated the distribution once, there will be a constant overhead from that step and you shouldn't see any appreciable difference even if I've underestimated the size of the data.