Monte Carlo Analysis

Monte Carlo analysis is a practical technique that has a long history and a ton of theory behind it. Fermi, Ulam and Von Neumann used statistical sampling ideas back in the 1930’s and 1940’s. The origins of statistical sampling date back to Laplace in the early 1800’s. The name Monte Carlo Analysis was suggested by Metropolis in 1946. Monte Carlo was used on the ENIAC computer to do neutron transport calculations in th mid 1940’s. Today, Monte Carlo analysis is utilized in all fields of research. The main assumption of this approach is that a randomly chosen sample tends to exhibit the same properties as the population from which it as drawn. Before we apply this technique to modeling stock prices, let’s take a look at a simple example.

Armed with our basic building block (Brownian Motion, we can go on to construct a plausible model for the behavior of actual stock prices. Before we proceed with constructing a model, let’s take a look at some of the stylized facts of actual stock prices.

plot(luck$Close, main = "The price of luck")

The first thing we notice is that this price series doesn’t appear to be stationary. In other words, there is no obvious mean price and it doesn’t make sense to talk about the standard deviation of the price. Working with such non-stationary timeseries is a hassle. If prices are not convenient to work with, then what should we use instead? Let’s take a look at the percentage returns of this stock.

returns_luck <- diff(log(as.numeric(as.character(luck$Close))))
plot(returns_luck, main = "LUCK % returns", col = "navyblue")

hist(returns_luck, breaks = 100, col="brown")

Apart from some clustering in the returns plot, it appears that the returns are distributed somewhat like a normal (Gaussian) distribution. This is an exciting fact since we already know how to work with normal distributions! How about independence? Are these returns independent of each other in time? Here’s a quick way to partially answer that question:

acf(returns_luck[-1], main = "Autocorrelation plot of returns")

Notice that there doesn’t seem to be any autocorrelation between consecutive returns. What are the mean and standard deviation of these returns?

mR  <- mean(returns_luck[-1])
sdR <- sd(returns_luck[-1])

mR
## [1] 0.001399149
sdR
## [1] 0.0251678

Leap of Faith

So, the typical argument goes as follows:

We want to deal with stationary timeseries since we have a ton of statistical tools available at our disposal that deal with such timeseries. Prices are surely non-stationary. Is there any other transformation of prices that, at least, looks like it might be stationary? It seems like percentage returns fit the bill. It also looks like percentage returns have a stable mean and standard deviation. So we can make the claim that percentage returns are normally distributed with mean ?? and standard deviation ??. Now, remember what our end goal is. We want a way to simulate stock prices. In order to do so, we need to come up with a model of how the prices behave (are distributed.) If returns are normally distributed, then how are prices distributed? The answer to this question is straightforward. A little math shows us the answer: Rt = log(Pt/Pt???1). The logarithm of the price is normally distributed. This means that price has a lognormal distribution. A straightforward method to simulate a stock price is to draw a random normal number with a certain mean and standard deviation value, and then exponentiate this number. Based on the formula from above: Pt = Pt???1*e^Rt. To summarize:

Draw a random number from N(??, ??). Exponentiate that number and mulitply it by Pt???1. Repeat for t = 1.N.

Portfolio simulation

We will use Lucky cement, habib bank and Fauji Fertilizer stock and simulate them.

psx_stocks <- c("OGDCL", "LUCK", "HBL", "FFC")

psx_data <- kse2 %>% select(c(1,2,6)) %>% filter(Symbols %in% psx_stocks) %>% 
            filter(date != "2008-09-22")

psx_data <- spread(psx_data,Symbols, Close)
library(zoo)
psx_data <- na.locf(psx_data)
library(lubridate)
psx_data$date <- dmy(psx_data$date)
psx_data <- psx_data %>% mutate_at(c(2,3,4), function(x) as.numeric(as.character(x)))


library(xts)
psx <- xts(psx_data[,-1], order.by = psx_data$date)

#compute returns matrix
rM <-  apply(psx,2,function(x) diff(log(x)))



#look at pairwise charts
pairs(coredata(rM))

#compute the covariance matrix
covR <- cov(rM, use = "pairwise.complete.obs")

library(MASS)
#use this covariance matrix to simulate normal random numbers
#that share a similar correlation structure with the actual data
meanV <- apply(rM, 2, function(x) mean(x, na.rm = TRUE))

rV    <- mvrnorm(n = nrow(rM), mu = meanV, Sigma = covR)

#simulate prices based on these correlated random variables

#calculate mean price
p0 <- sapply(psx, function(x) mean(x, na.rm = TRUE))
sPL <- list()
for(i in 1:ncol(rM)){
  sPL[[i]] <-round(p0[i]*exp(cumsum(rV[,i])),2)
}

#plot simulated prices
par(mfrow = c(1,3)) 
plot(sPL[[1]][1:500],main="FFC sim",type="l")
plot(sPL[[2]][1:500], main = "HBL sim",type = "l")
plot(sPL[[3]][1:500], main = "LUCK sim", type = "l")

In the prior example, we gather daily data for 3 pakistani stocks and we compute the covariance matrix of the returns, along with an average price for each security. Since the purpose of this exercise is to generate a realistic simulation of the portfolio, we use the function mvrnorm() to create a matrix of random normal variables that are correlated in a similar manner as the original data. The following graphs display the original stock prices, the pairwise plot of their returns and the simulated stock prices. One can also look at the correlation matrix of the actual returns and the simulated returns and verify that they are similar.