Monte Carlo Method in R (with worked examples)

Monte Carlo method is a handy tool for transforming problems of probabilistic nature into deterministic computations using the law of large numbers. Imagine that you want to asses the future value of your investments and see what is the worst-case scenario for a given level of probability. Or that you want to plan the production of your factory given past daily performance of individual workers to ensure that you will meet a tough delivery plan with high enough probability. For such and many more real-life tasks you can use the Monte Carlo method.

Monte Carlo approximation of Pi

The method generally works in the following steps:

  1. Define a domain of possible inputs; e.g. what might be the price of a stock after 20 days; or how many pieces the worker can make in a day
  2. Generate inputs randomly from a probability distribution over the domain; depends strongly on the problem at hand; typically uniform, normal, beta, Bernoulli or Poisson distributions are used
  3. Perform a deterministic computation on the inputs; e.g. sum the productivity of individual workers
  4. Aggregate the results; e.g. take the lower bound of 99% confidence interval

Worst-case scenario for future stock price

Let’s see on a simple example how easy is to perform Monte Carlo method in R. Assume we want to calculate the worst-case scenario of a future stock price. This problem called value at risk is heavily used in risk management. By “worst-case scenario” we mean the value that the stock price will exceed with 99% probability (i.e., there is only 1% probability that the stock price will be below). The current price of our stock is 100 $. We want to see the possible future prices after 20 trading days.

As stated above, we need to make some assumptions about the future stock price. So let’s say that the drift over the 20 trading days is 10% (i.e. in average the price will go up to 110) and the volatility is 20%.

Second step is to generate some random inputs, which in our case means to model the future price given the current price, drift and the volatility. For our purpose we do not need anything fancy so we will use the standard stock price model called Geometric Brownian Motion:

St+1 = St * (1 + μΔt + σε√Δt)

Where Sis price in time t, μ is our drift, Δt means one period (in our case 1 trading day, i.e. one twentieth of the period), σ is our volatility and ε is a random number between 0 and 1. The core idea of Monte Carlo method is to generate the future price (which is random) high number of times to simulate what are all the situations that can occur.

In the third step we calculate the stock price at the end of 20th day given the model for all the randomly generated daily returns. This step (and also the previous ones) is demonstrated by the following code.

#' Stock price calculation
#' Calculates stock price after n periods using standard stock price model
#' @param stock_price original stock price
#' @param n number of periods
#' @param stock_mu expected percentual stock drift over n periods
#' @param stock_sigma expecter percentual stock volatility
#' @return stock price after n periods
f_stock_return <- function(stock_price, n, stock_mu, stock_sigma){
  delta_t <- 1/n # one period
  for (i in seq(n)){
    epsilon <- runif(n=1, min=0, max=1) # random generated number
    # calculate stock price (using quantile function of normal distribution)
    stock_price <- stock_price * (1 + qnorm(epsilon, 
                                            stock_mu * delta_t, 
                                            stock_sigma* sqrt(delta_t)))

This needs to be run again and again, e.g. 10 000 times.

# parameters
simulations <- 10000 # number of MC simulations
n <- 20 # trading days
stock_price <- 100
stock_mu <- .1 # drift 10%
stock_sigma <- .2 # volatility 20%

# Monte Carlo simulations
set.seed(42) # for reproducibility
stock_prices <- c()
for (i in seq(simulations)){
  stock_prices <- c(stock_prices,

The last step is to aggregate the results and answer our question.

quantile(stock_prices, c(.01, .05))

We can see that in 9 900 (99%) scenarios the price is larger than ca. 67 $; 95% of them were higher than ca. 77 $. Which gives us needful assessment of how bad our investment could go.

Approximation of Pi

Monte Carlo method can be applied also to problems that can be reformulated to have probabilistic interpretation. Very popular example is the approximation of the number Pi.

This example is based on the fact that if you randomly generate points in a square, π/4 of them should lie within an inscribed circle.

Assuming the perimeter of the circle is r, area of the square is equal to 4r2 and area of the inscribed circle is πr2. Therefore the probability of a point lying inside the inscribed circle is π/4.

# parameters
simulations <- 10000 # number of simulations
perimeter <- 1

# randomly generate a point and check if it is in circle
f_point_in_circle <- function(perimeter=1){
  x <- runif(n=1, min=-perimeter, max=perimeter)
  y <- runif(n=1, min=-perimeter, max=perimeter)
              in_circle=x^2 + y^2 <= perimeter^2))

# Monte Carlo simulations
pi_df <- data.frame(x=rep(NA, simulations),
                    y=rep(NA, simulations),
                    in_circle=rep(NA, simulations))
for (i in seq(simulations)){
  my_simulation <- f_point_in_circle()
  pi_df$in_circle[i] <- my_simulation$in_circle
  pi_df$x[i] <- my_simulation$x
  pi_df$y[i] <- my_simulation$y

my_pi <- 4 * sum(pi_df$in_circle) / nrow(pi_df)

This code generates an approximation of π equal to 3.1416. The more simulations we run, the more accurate it gets.

Following script generated the picture above.

ggplot(pi_df, aes(x=x, y=y, color=as.factor(in_circle))) +
  geom_point() +

2 thoughts on “Monte Carlo Method in R (with worked examples)

  1. Teddy

    The stock price example confuses me. I dont understand why we would need to perform monte carlo simulation to find out that in 95% of scenarios the price is larger than x. If we have the mean and standard deviation of a normal distribution then can we not just figure the answer out analytically? ie: 100 * (1 + qnorm(.095, stock_mu * delta_t, stock_sigma * sqrt(delta_t)))



Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.