Simulate correlated Geometric Brownian Motion in the R programming language

Question

In response to this question: How to simulate correlated Geometric brownian motion for n assets?

One of the responses provides an implementation in MATLAB: http://www.goddardconsulting.ca/matlab-monte-carlo-assetpaths-corr.html

I attempted to port this code into the R programming language. The implementation is:

# Input Params:
S0 <- c(50, 48)
mu <- c(0.03, 0.03)
sig <- c(0.05, 0.05)
corr <- cbind(c(1,0.1), c(0.1,1))
dt <- 1/365
steps <- 10000
nsims <- 100

# Get the number of assets:
nAssets <- length(S0)

# Calculate the drift:
nu <- mu - sig * sig/2

# Do a Cholesky Factorization on the correlation matrix:
R <- chol(corr)

# pre allocate the output:

S <- array(1, dim=c(steps+1, nsims, nAssets))

# Generate correlated random sequences and paths:

for(idx in 1:nsims)
{
  # generate uncorrelated random sequence 
  x <- matrix(rnorm(steps * nAssets), ncol = nAssets, nrow = steps)

  # correlate the sequences 
  ep <- x %*% R

  #generate potential paths 

  S[,idx,] <- rbind(rep(1,nAssets), apply(exp(matrix(nu*dt,nrow=steps,ncol=2,byrow=TRUE) + (ep %*% diag(sig)*sqrt(dt))), 2, function(x) cumprod(x)) ) %*% diag(S0)
}

# TESTING: Compute Average Sample Correlation 
sum = 0
for(i in 1:nsims) 
{
  sum = sum + cor(S[,i,1], S[,i,2])
}
sampleCorrelation = sum / nsims 

To test if the implementation works as promised, I compute the average sample correlation between two assets across many simulations. By the law of large numbers, the sample correlation should be pretty close to the theoretical correlation matrix provided as an input parameter. In the example provided, the sample correlation between the two assets should be pretty close to 0.1 on average. However, this is not the case. So there must be an issue in the R code provided. I also tested the MATLAB code on a free online MATLAB simulator (http://octave-online.net/): it gave the correct result for the average sample correlation between the assets. Hence, there must be a porting issue, which I can't seem to determine. Can you identify the issue?

Answer 1

Have you solved it yet? For example in the drift parameter, the dt needs to be vector of time from 0 to 1 by dt.

My code is:

  GBM<-apply(BM,2,function(x) 100*exp((cumsum((r-0.5*sigma*sigma)*time)+sigma*x)))

where I'm using GBM on already cumsummed Brownian Motion (x).

Answer 2

You can use a multivariate standard normal distribution to achieve the desired results.

require(mvtnorm)
require(matrix)

Covariance.matrix <- # your covariance matrix here
Drift <- # your drift terms here
Vol <- # your vol terms here
n <- # desired number of samples
z <- rmvnorm(n, simga=nearPD(Covariance.matrix), method="chol")
GBM <- Drift + z * Vol

If you don't want to assume normality, you can use Cornish-Fisher or something similar to adjust for skew and kurtosis for each individual currency.

Answer 3

Make sure you do the RETURN correlation, not the price correlation. Also, if number is not converging, increase the amount of steps so that it is sufficient enough for the law of large numbers kick in.