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?