**Monte Carlo Methods**

Monte Carlo methods, named after the area in Monaco famous for it casinos, are a collection of mathematical methods that use repeatedly sampling of random numbers to solve mathematical problems. They are commonly used to simplify problems that would otherwise be impossible or too time consuming to solve if we were forced to use analytical methods.

**An example**

Suppose we are given the option of playing the following game in a casino:

The game costs £5 per play.

The game involves you rolling 10 dice. If the sum of the dice is greater than 45, you win £10, otherwise you don't win anything.

Should you play this game?

Off the top of my head, I've got no idea whether this is good value or not. Given we are in a casino, I'm going to guess we would lose money over the long run, but what if we wanted to know for sure? We could work out the answer analytically, we could say things like the probability of rolling a 45 using 7 dice is impossible, but using 8 dice it is.... then work recursively from there and it would take a long time and would be quite tedious. An alternative approach would be to set up a model that would simulate the playing of this game, run the model a few thousand times (which should take seconds) and then see what the long term position is when we average across all the games. Much more interesting.

I quickly set up a test of the game, and it turns out it is very bad value. On average we would expect to lose about £4.60 every time we play! The interesting point though is that it took me a couple of minutes to work this out, compared to the hours it would have taken me to calculate it analytically.

**Monte Carlo Methods in insurance**

Monte Carlo methods are also commonly used in insurance. For example, in reinsurance pricing, when attempting to find an expected value for a loss to an excess of loss layer we may have a good idea of the size and frequency of the ground up losses, say we know that the frequency of losses is distributed according to a Poisson distribution with a known parameter, and the severity of the losses is distributed according to a Log-normal distribution with known parameters. What can we say about the loss to the layer however? It may be theoretically possible to derive this analytically in some cases, but generally this will be too time consuming and may not even be possible.

A Monte Carlo approach, like the dice game earlier, involves repeatedly playing the game, we imagine that we write the contract and we then sample values from the frequency and severity distributions to generate sample loss amounts, we can then apply all the contract details to these generated loss amounts. This process of sampling is repeated tens of thousands of times and by averaging across all our simulations we can derive estimates of the expected loss to the layer or any other piece of information we may be interested in.

**Other uses of Monte Carlo Methods**

We can also use Monte Carlo methods in other situations where we wish to estimate a value where analytical methods may fall short.

hold on

x = -1:0.001:1;

y = (1 - x.^2).^0.5;

plot(x,y,'linewidth',2);

plot(x,-y,'linewidth',2);

plot([1,-1],[-1,-1],'linewidth',2);

plot([-1,-1],[-1,1],'linewidth',2);

plot([-1,1],[1,1],'linewidth',2);

plot([1,1],[1,-1],'linewidth',2);hold on

x = -1:0.001:1;

y = (1 - x.^2).^0.5;

plot(x,y,'linewidth',2);

plot(x,-y,'linewidth',2);

plot([1,-1],[-1,-1],'linewidth',2);

plot([-1,-1],[-1,1],'linewidth',2);

plot([-1,1],[1,1],'linewidth',2);

plot([1,1],[1,-1],'linewidth',2);

Given the ratio of the circle to the square, we would expect the proportion of randomly selected sample points which end up inside the circle to approach $\pi /4$. And we can therefore use this value to estimate $\pi$.

**hold on**

x = -1:0.001:1;

y = (1 - x.^2).^0.5;

plot(x,y,'linewidth',2);

plot(x,-y,'linewidth',2);

plot([1,-1],[-1,-1],'linewidth',2);

plot([-1,-1],[-1,1],'linewidth',2);

plot([-1,1],[1,1],'linewidth',2);

plot([1,1],[1,-1],'linewidth',2);

SampleSize = 1000;

Sample = unifrnd (-1, 1, 2, SampleSize);

Results = 1:SampleSize;

for i = 1:SampleSize

Results(i) = Sample(1,i).^2 + Sample(2,i).^2;

endfor

for i = 1:SampleSize

if Results(i) <= 1

plot(Sample(1,i),Sample(2,i),'k');

else

plot(Sample(1,i),Sample(2,i),'r');

endif

endforx = -1:0.001:1;

y = (1 - x.^2).^0.5;

plot(x,y,'linewidth',2);

plot(x,-y,'linewidth',2);

plot([1,-1],[-1,-1],'linewidth',2);

plot([-1,-1],[-1,1],'linewidth',2);

plot([-1,1],[1,1],'linewidth',2);

plot([1,1],[1,-1],'linewidth',2);

SampleSize = 1000;

Sample = unifrnd (-1, 1, 2, SampleSize);

Results = 1:SampleSize;

for i = 1:SampleSize

Results(i) = Sample(1,i).^2 + Sample(2,i).^2;

endfor

for i = 1:SampleSize

if Results(i) <= 1

plot(Sample(1,i),Sample(2,i),'k');

else

plot(Sample(1,i),Sample(2,i),'r');

endif

endfor

The final code is:

**SampleSize = 100000;**

Sample = unifrnd (-1, 1, 2, SampleSize);

Results = 1:SampleSize;

for i = 1:SampleSize

Results(i) = Sample(1,i).^2 + Sample(2,i).^2;

endfor

Pi = 0;

for i = 1:SampleSize;

if Results(i) <= 1

Pi = Pi + 1/SampleSize;

endif

endfor

4*PiSample = unifrnd (-1, 1, 2, SampleSize);

Results = 1:SampleSize;

for i = 1:SampleSize

Results(i) = Sample(1,i).^2 + Sample(2,i).^2;

endfor

Pi = 0;

for i = 1:SampleSize;

if Results(i) <= 1

Pi = Pi + 1/SampleSize;

endif

endfor

4*Pi

Why not give it a go yourself!