Skip to content
All library documents

Simulating Correlated Geometric Brownian Motion Across Assets

Article Quant Q&A · Author: Ryan J. Shrott

Summary

The document examines how to simulate multiple geometric Brownian motion price paths with a specified correlation structure. The attempted R approach draws independent normal shocks, applies a Cholesky factor of the input correlation matrix, and uses the resulting shocks in asset paths. The author checks the output by averaging sample correlations across simulated paths, expecting them to reflect the input correlation.

The responses point to several considerations: accumulate Brownian increments over time and scale drift using the time grid; assess correlation on returns rather than price levels; and allow enough observations for sample estimates to stabilize. Another suggestion is to draw multivariate normal samples directly from a covariance matrix, with a possible adjustment if normality is not desired. The discussion does not settle on one verified correction to the posted implementation, and the example's average price correlation is not itself a reliable validation of the target return correlation.

Key ideas

  • Correlated asset shocks can be generated by transforming independent normal draws with a Cholesky factor.
  • Geometric Brownian motion requires cumulative Brownian increments and time-scaled drift and diffusion terms.
  • Validate the simulated dependence using returns rather than correlations between price levels.
  • Multivariate normal sampling is an alternative way to generate correlated shocks.
  • Finite samples may not reproduce the target correlation closely, and non-normal shock models require further assumptions.

Tags

Full text
# Simulate correlated Geometric Brownian Motion in the R programming language


# Simulate correlated Geometric Brownian Motion in the R programming language












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 by Jan Sila (score 3)

https://quant.stackexchange.com/a/28174

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 by milkmotel (score 1)

https://quant.stackexchange.com/a/34042

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 by Laxus (score 0)

https://quant.stackexchange.com/a/34040

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.

Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)

This summary was written by Stratmill's research agent from the original; it is not a copy of the source.