107: Stat Introduction to Business and Financial Statistics
Class 4: Monte Carlo Simulation
1
Historical Note
The name Monte Carlo simulation comes
from the computer simulations performed
during the 1930s and 1940s to estimate the
probability that the chain reaction needed for
an atom bomb to detonate would work
successfully.
The physicists involved in this work were big
fans of gambling, so they gave the
simulations the code name Monte Carlo.
2
Uses of Monte Carlo Simulation
Ford, Proctor and Gamble, Pfizer, Bristol-Myers Squibb,
and Eli Lilly use simulation to estimate both the average
return and the risk factor of new products.
Proctor and Gamble uses simulation to model and
optimally hedge foreign exchange risk.
Sears uses simulation to determine how many units of
each product line should be ordered from suppliersfor
example, the number of pairs of Dockers trousers that
should be ordered this year.
Financial planners use Monte Carlo simulation to
determine optimal investment strategies for their clients
retirement.
3
The Monte Carlo Method
The method is usually used to understand how
a process or estimator will perform in practice.
The Monte Carlo method is often performed
when the computational resources of a
researcher are more abundant than the
researchers mental resources.
4
Wikipedia
fm3: Chapter 22
Monte Carlo
5
Warning
Monte Carlo is often time-consuming,
computationally wasteful
Our first example (computing ) is a good
illustration
Analytical methods (formulas ) are better if
available
But
Sometimes theres nothing else
Monte Carlo is fun!
6
Example: Estimating
Did you know that you can calculate
(3.14159..) by throwing darts at a dartboard ?
Unfortunately, you have to be pretty bad at
darts
darts for it to work, though.
7
Geometry Review
Suppose we have a circle of radius 1, with a
center at (0,0).
Then the unit circle is given by the equation
x2+y2=1.
8
The Set Up
Suppose our dartboard is a 2x2 square and
we inscribe a circle of radius 1 inside the
square, with a center at (0,0).
9
Computing :
Unit circle inscribed in unit square
1.0
The square has area 4
Unit circle has area
0.5
1.0
0.5
0.5
1.0
0.5
1.0
10
The Algorithm
We randomly throw darts at this square and
keep track of how many land inside the circle.
That would give us a relative estimation of what
percentage of the box is in the circle.
The area of the box is 4 and the area of the
circle is .
So the ratio of circle area to box area is /4.
Hence, the percentage of darts that land inside
the circle should be approximately /4.
11
The runif() command
To simulate throwing darts at the dartboard,
we need to learn about the R command
runif():
R Command: runif(n)
uniform random number generator that returns
n random numbers between 0 and 1.
12
Examples
> runif(5)
[1] 0.7000905 0.6717379 0.9857697 0.6081811 0.8611683
>
> runif(5)
[1] 0.1719788 0.5086680 0.6567634 0.1922734 0.8315596
>
> foo=runif(1000)
> hist(foo)
>
13
Our circle is as follows:
1
How do we simulate
darts thrown at this
board ?
The location of the
dart is given by its
(x,y) point.
(0,0)
y
-1
x
1
14
Code Fragment
Examine the following code fragment
x = 2*runif(1)-1
y = 2*runif(1)-1
if (x*x+y*y<=1) numincircle = numincircle+1
What are the possible values for x and y ?
15
The Complete Function
simulpi = function(n) {
numincircle = 0
Keep count of darts inside
circle
do n simulations
for(i in 1:n) {
x = 2*runif(1)-1
y = 2*runif(1)-1
if (x*x+y*y<1) numincircle <- numincircle+1
}
cat("Number of simulations = ",n,"\n")
cat("Estimate of PI = ",4*numincircle/n,"\n")
}
16
Output
> simulpi(10)
Number of simulations =
Estimate of PI = 3.2
> simulpi(100)
Number of simulations =
Estimate of PI = 3.04
> simulpi(1000)
Number of simulations =
Estimate of PI = 3.136
> simulpi(1000)
Number of simulations =
Estimate of PI = 3.184
> simulpi(100000)
Number of simulations =
Estimate of PI = 3.1472
10
100
1000
1000 iterations-why arent the
values the same each time ?
1000
1e+05
17
1.0
0.5
0.0
-0.5
-1.0
y1
-1.0
-0.5
0.0
x1
0.5
1.0
18
In Excel
A
B
C
D
E
COMPUTING PI USING MONTE CARLO METHODS
INITIAL EXPERIMENT
1
2 Number of data points
3 Inside circle
4 Pi?
5
Each cell in these columns contains the Excel
function =Rand()
6
7
8
9
10
30 <-- =COUNT(A:A)
23 <-- =COUNTIF(D:D,TRUE)
3.066666667 <-- =B3/B2*4
Experiment
1
2
3
Random1
0.23437
0.92322
0.95535
Random2
0.03080
0.61267
0.25644
In unit
circle?
TRUE
FALSE
TRUE
<-- =(B8^2+C8^2<=1)
Columns B and C have function Rand().
Column D uses Boolean function to determine if point is inside
unit circle
Cell B4 approximates p by
Computing ratio of inside/total points
Multiplying this ratio by 4
19
Different experiments give different
approximations
A
B
C
D
E
COMPUTING PI USING MONTE CARLO METHODS
INITIAL EXPERIMENT
1
2 Number of data points
3 Inside circle
4 Pi?
5
Each cell in these columns contains the Excel
function =Rand()
6
7
8
9
10
30 <-- =COUNT(A:A)
21 <-- =COUNTIF(D:D,TRUE)
2.8 <-- =B3/B2*4
Experiment
1
2
3
Random1
0.88255
0.37558
0.81044
In unit
circle?
TRUE
TRUE
FALSE
Random2
0.06206
0.41664
0.65635
A
B
<-- =(B8^2+C8^2<=1)
C
D
E
COMPUTING PI USING MONTE CARLO METHODS
INITIAL EXPERIMENT
1
2 Number of data points
3 Inside circle
4 Pi?
5
Each cell in these columns contains the Excel
function =Rand()
6
7
8
9
10
30 <-- =COUNT(A:A)
24 <-- =COUNTIF(D:D,TRUE)
3.2 <-- =B3/B2*4
Experiment
1
2
3
Random1
0.27544
0.81204
0.91731
Random2
0.04368
0.52793
0.90177
In unit
circle?
TRUE
TRUE
FALSE
<-- =(B8^2+C8^2<=1)
20
Stock Returns
Suppose you put $1000 into an S&P 500
index fund on January 1. How much will you
have at the end of the year?
We can easily simulate this using Monte
Carlo.
How do we simulate yearly stock returns?
There are several methods.
21
Simulating Stock Returns
Assuming normality of returns is not that
good an idea since stock returns demonstrate
tail behavior, but it is less common than what
the normal distribution would dictate.
We will produce the returns using a
procedure called resampling. We have a
dataset with the yearly S&P 500 returns from
1926 to 2010.
22
The Historical Returns
A snapshot
23
The Histogram of Historical Returns
24
Resampling
To produce a return in the future for the index,
we resample from our 83 historical observationsthat is, we just randomly pick one out of the 83
past yearly returns.
This scheme doesnt build in any ideas like good
years are more likely to be followed by good
years, etc but it is a technique used by several
software packages that do financial planning.
To do this in R, we use the sample command
25
The SAMPLE command
The command we use is
sample(x, size=<<see below>>, replace=F)
Collection of items to
sample from
How many times to
sample
Sample with or without
replacement
26
Example Output
> 0.2380297
> sample(rets,1)
[1] sample(rets,1)
[1] 0.09967052
> sample(rets,1)
[1] 0.06146142
> sample(rets,1)
[1] 0.3023484
> sample(rets,1)
[1] -0.01188566
> sample(rets,1)
[1] -0.01208205
27
Very basic idea
This is called a nonparametric routine-it does
not assume any structure on the returns.
This method assumes the historical data are
essential the population of possible returns.
The sampled data looks very much like the
historical data.
> summary(rets)
Min. 1st Qu.
Median
Mean 3rd Qu.
-0.43840 -0.02136 0.14220 0.11310 0.25490
> summary(sample(rets,100000,replace=TRUE))
Min. 1st Qu.
Median
Mean 3rd Qu.
-0.43840 -0.03064 0.14220 0.11220 0.25920
Max.
0.52560
Max.
0.52560
28
The Monte Carlo
How many
iterations
We are now ready to code
year.ret = function(rets,howmany) {
store = 1:howmany
Creates a place
to
to store our
results
for(i in 1:howmany)
store[i] = 1000*(1+sample(rets,1))
return(store)
}
29
Running the Monte Carlo
Cut and paste the code into R
> out=year.ret(rets,10000)
> hist(out)
The output is not normal.
Should it be?
30
Summarize the Results
> summary(out)
Min. 1st Qu.
561.6
969.4
>
Median
1142.0
Mean 3rd Qu.
1113.0 1259.0
Max.
1526.0
The probability of making a positive return is
> sum(out>1000)/10000
[1] 0.708
This makes sense since
> sum(rets>0)/83
[1] 0.7108434
71% of the historical returns are positive
so it makes sense that we made more
than $1000 71% of the time.
31
Convert to CAGRWhat is CAGR?
CAGR=compound annual growth rate
The compound annual growth rate is calculated
by taking the nth root of the total percentage
growth rate, where n is the number of years in
the period being considered.
32
Huh?
CAGR isn't the actual return in reality. It's an imaginary number that
describes the rate at which an investment would have grown if it
grew at a steady rate. You can think of CAGR as a way to smooth
out the returns.
Suppose you take a look at a stock price over ten years, and the
price
price went from P0 (10 years ago) to P1 (today).
Then the CAGR is defined so that
10
P1 = P0 (1 + CAGR)
33
Can Convert to CAGR if desired
Recall CAGR
> summary(out)
Min. 1st Qu.
561.6
969.4
>
CAGR= - 43.8%
Median
1142.0
Mean 3rd Qu.
1113.0 1259.0
Max.
1526.0
CAGR=52.6%
But this isnt an interesting example since it is
just what happens in one year.
34
Growth over 30 years
We can do the same idea for 30 years in a
row
money = 1000
for(i in 1:30) {
ret = sample(rets,size=1)
money = money*(1+ret)
}
cat("In 30 years you have ",round(money),"\n")
35
Lets make it into a function
mysim = function(numyears) {
money = 1000
for(i in 1:numyears) {
ret = sample(rets,size=1)
money = money*(1+ret)
}
cat("In ",numyears," years you have ",round(money),"\n")
}
36
Running this code in R
> mysim(30)
In 30 years
> mysim(30)
In 30 years
> mysim(30)
In 30 years
> mysim(30)
In 30 years
> mysim(30)
In 30 years
> mysim(30)
In 30 years
>
you have
66089
you have
8522
you have
20847
you have
38044
you have
22745
you have
148696
37
Keeping Track of the Results
simul=function(rets,howmany) {
values = 1:howmany
for(j in 1:howmany) {
money = 1000
for(i in 1:30) {
ret = sample(rets,size=1)
money = money*(1+ret)
This creates a vector to
store our results
This double looping
makes the problem
difficult to do in
Excel
Excel
}
values[j] = money
}
return(values)
}
38
Run the code in R
Cut and paste it into R
> summary(out)
Min. 1st Qu.
Median
249.1
7229.0 15370.0
>
CAGR= 9.5%
> (15370/1000)^(1/30) - 1
[1] 0.09535727
>
Mean
24810.0
3rd Qu.
Max.
30870.0 489900.0
CAGR=22.9%
> (489900/1000)^(1/30) - 1
[1] 0.229335
39
Run the code in R
hist(out,breaks=50)
40
Go Normal
What if we wanted to assume the returns
were normal?
Sample statistics
> mean(rets)
[1] 0.1131480
> sd(rets)
[1] 0.2021068
41
The rnorm command
simul=function(howmany) {
values = 1:howmany
for(j in 1:howmany) {
money = 1000
for(i in 1:30) {
ret=rnorm(1,mean=.1131,s=.2021)
money = money*(1+ret)
}
values[j] = money
}
return(values)
}
42
Different Results (of course)
> summary(out)
Min. 1st Qu.
119.2
7360.0
>
> summary(out)
Min. 1st Qu.
249.1
249.1
7229.0
>
Normal
Median
14490.0
Mean
24240.0
3rd Qu.
Max.
28670.0 454000.0
Resampled
Median
15370.0
Mean
24810.0
3rd Qu.
Max.
30870.0 489900.0
Which model for returns allows for slightly more
upside?
Which is more realistic, which would we show a
client?
What would logspline routine show?
43
Go Semiparametric
Density Estimation refers to determining a
probability density for a given set of data.
This can be as simple as assuming the data
is normal, and using the sample mean and
standard deviation as parameters.
Or it can involve some fancy math.
44
What do densities require
A probability density is a continuous function
f(x) so that
f(x) >=0 for all values of x
The function f(x) integrates to 1.
45
The logspline package
library(logspline)
Performs a semi-parametric density estimate,
then lets you sample from the estimate.
Pretty cool routine.
46
The historical stock returns
Just a few lines of code
fit=logspline(rets)
> hist(rets,probability=TRUE)
> x=seq(-.6,.6,.01)
> lines(x,dlogspline(x,fit))
47
The resulting density
48
Compare with the Normal
Longer tail
49
Can draw from this density
A powerful ability of the logspline routine is
that one can draw random values from the
fitted density function.
This is easily done using the rlogspline
routine.
50
Another Simulation
fit = logspline(rets)
simul=function(fit,howmany) {
values = 1:howmany
for(j in 1:howmany) {
money = 1000
for(i in 1:30) {
ret=rlogspline(1,fit)
money = money*(1+ret)
}
values[j] = money
}
return(values)
}
51
Results
More long tailed behavior
Min. 1st Qu.
-41290
6102
Median
14000
> plogspline(-.3,fit)
[1] 0.03610169
> pnorm(-.3,mean(rets),sd(rets))
[1] 0.02046658
Mean 3rd Qu.
24950
29360
Max.
485700
Almost double the
probability of a return being
less than -0.3!
52
Example: Retirement Planning (hw)
Joe Renk is a 35 year old freelance writer with $30,000
in an IRA. He figures he can put in $2000 (the maximum
allowed every year) for the next 30 years.
John wants to have a simple portfolio so he decides to
put all his money into an S&P500 index fund, which has
an annual management fee of 0.5%.
How much money will John have in 30 years ?
What is the probability he is a millionaire ?
53
A code fragment (modify for hw)
money = 30000
for(i in 1:30) {
ret = sample(rets,size=1)
money = (money+2000)*(1+ret)
fee = .005*money
money = money-fee
}
cat("In 30 years Joe has ",round(money),"\n")
Cut and pasting this into R, you get results such as
In 30 years Joe has
In 30 years Joe has
In 30 years Joe has
197919
364987
1040814
54
Run this code many times
values = 1:1000
for(j in 1:1000) {
money = 30000
for(i in 1:30) {
ret = sample(rets,size=1)
money = (money+2000)*(1+ret)
fee = .005*money
money = money-fee
}
values[j] = money
}
55
Results:
> summary(values)
Min. 1st Qu. Median
Mean 3rd Qu.
Max.
36930 207900 353700 501400 608200 9083000
> length(values)
[1] 1000
> sum(values >= 1000000)
[1] 108
> 108/1000
[1] 0.108
> sum(values<100000)
[1] 65
> sum(values>1400000)
[1] 46
11% chance john will have
a million or more.
6.5% chance john will have less
than 100000.
4.6% chance john will have more
than 1.4 million
56
300000
200000
100000
Value
400000
500000
600000
Middle 50% of Final Values
0
5
10
15
20
25
30
Time
57
Our numbers are close to financialengines.com, a website
founded by Nobel Prize winner Bill Sharpe:
58