Quantitative Finance Applications in R - 5: an Introduction to Monte Carlo Simulation

by Daniel Hanson

Last time, we looked at the four-parameter Generalized Lambda Distribution, as a method of incorporating skew and kurtosis into an estimated distribution of market returns, and capturing the typical fat tails that the normal distribution cannot.  Having said that, however, the Normal distribution can be useful in constructing Monte Carlo simulations, and it is still commonly found in applications such as calculating the Value at Risk (VaR) of a portfolio, pricing options, and estimating the liabilities in variable annuity contracts. 

We will start here with a simple example using R, focusing on a single security.  Although perhaps seemingly trivial, this lays the foundation for in more complexities such as multiple correlated securities and stochastic interest rates.  Discussion of these topics is planned for articles to come, as well as topics in option pricing.

Single Security Example

Under the oft-used assumption of Brownian Motion dynamics, the return of a single security (eg, an equity) over a period of time Δt is approximately [See Pelsser for example.]

μΔt + σZ・sqrt(Δt)                                                                (*)

where μ is the mean annual return of the equity (also called the drift), and σ is its annualized volatility (i.e., standard deviation).  Z is a standard Normal random variable, which makes the second term in the expression stochastic.  The time t is measured in units of years, so for quarterly returns, for example, Δt = 0.25.

As μ, σ, and Δt are all known values, generating a simulated distribution of returns is a simple task.  As an example, suppose we are interested in constructing a distribution of quarterly returns, where μ = 10% and σ= 15%.  In order to get a reasonable approximation of the distribution, we will generate n = 10,000 returns. 

n <- 10000
# Fixing the seed gives us a consistent set of simulated returns
set.seed(106)
z <- rnorm(n)        # mean = 0 and sd = 1 are defaults
mu <- 0.10
sd <- 0.15
delta_t <- 0.25
# apply to expression (*) above
qtr_returns <- mu*delta_t + sd*z*sqrt(delta_t)  

Note that R is “smart enough” here by adding the scalar mu*delta_t to each element of the vector in the second term, thus giving us a set of 10,000 simulated returns. Finally, let’s check out results.  First, we plot a histogram:

hist(qtr_returns, breaks = 100, col = "green")

This gives us the following:

The symmetric bell shape of the histogram is consistent with the Normal assumption.  Checking the annualized mean and variance of the simulated returns,

stats <- c(mean(qtr_returns) * 4, sd(qtr_returns) * 2)   # sqrt(4)
names(stats) <- c("mean", "volatility")
stats

We get:

      mean    volatility
0.09901252    0.14975805

which is very close to our original parameter settings of μ = 10% and σ= 15%.

Again, this is rather simple example, but in future discussions, we will see how it extends to using Monte Carlo simulation for option pricing and risk management models.

Comments

[Rob]

Thanks very much. I know this is a separate issue, but I was looking at GLD post that is referenced in this example and I did not understand the rgl() function. I could not find it and it did not work for me. The comments to that post are closed so I thought that I would try here. I am enjoying the series.

Rob