# How To Make Monte Carlo Simulation Run Faster

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

• With SkewNormalDistribution[Mean[data2], StandardDeviation[data2], Skewness[data2]], data2 doesn't change during the calculation. Save it as a distribution once dist = SkewNormalDistribution[Mean[data2], StandardDeviation[data2], Skewness[data2]] so that you do not recalculate the statistics of the distribution needlessly many times. Then use dist in the RandomVariate. Feb 18 at 20:14
• Yes I missed that. Thanks! Feb 18 at 20:27

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.

• I implemented your recommendations and my code has sped up a lot. I'm now using FindDistribution instead of NormalSkewed and it is a lot slower. Is there anyway to make FindDistribution[data] run faster by compiling it or using an internal function. Thanks! Feb 18 at 23:54
• I would be surprised if FindDistribution could be sped up significantly by directly substituting internal functions into it. If you have knowledge about what kinds of distributions you're looking at, you might be best trying to fit each of them individually instead of using the general function. Might be worth looking at LearnDistribution too, but I haven't used it myself so I don't know about its speed. Feb 18 at 23:58

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.