Applying Cholesky Decomposition to Simulate Correlated Brownian Motions

# Cholesky decomposition reduces volatility of simulated Wiener Process / Brownian Motions


# Cholesky decomposition reduces volatility of simulated Wiener Process / Brownian Motions












I am trying to simulate $n$ correlated geometric brownian motions (GBM) given a specified correlation matrix $\Sigma$ by following this procedure which uses Cholesky decomposition.

However, when I implement the code in Python one of the highly correlated ($\rho=0.99$) Wiener Processes have almost no volatility as can be seen in the plot below in the red line.

I would expect a correlation of $\rho=1$ to imply that the two processes would be identical and for $\rho=-1$, I would expect a "mirror" around 0.

The plot shows the realization of $W$ across $t$.

Note: In my attempt I try to avoid using any loops.

```
# PARAMETERS
n = 2                                              # number of assets
T = 5                                              # number of years
N = 252                                            # number of steps pr. year (T)
r = np.array([0.03, 0.03])                         # drift-rate
q = np.array([0.0, 0.0])                           # divivend-rate
S0 = np.array([1.0, 1.0])                          # initial value
sigma = np.array([0.2, 0.2])                       # diffusion-coefficients 
corr_mat = np.array([[1.0, 0.99], [0.99, 1.0]])    # correlation-matrix

# Calculate step size
dt = 1 / (T * N)

# Vector of time index
t = np.array(range(N*T+1)).reshape((N * T + 1, 1)) * dt

# == PERFORM SIMULATION ==
# Draw normal distributed random variables
eps = np.random.normal(loc=0, scale=np.sqrt(dt), size=(N * T, n))

# Use Cholesky decomposition to get the matrix R satisfying: corr_mat = L x L^*
L = np.linalg.cholesky(corr_mat)

# Transform into realizations of correlated Wiener processes
W = eps @ L

# Calculate realizations of each stock as the realization of a GBM
# S_t = S0 * exp{ ((r-q) - ½sigma)*dt + sigma * W}
S = np.exp(
    (((r - q) - np.power(sigma, 2) / 2) * dt).T + W * sigma.T
)
S = np.vstack([np.ones(S.shape[1]), S])
S = S0 * S.cumprod(axis=0)
```

## Answer by jaehyukchoi49 (score 3, accepted)

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

Fix the following part

```
# Transform into realizations of correlated Wiener processes
W = eps @ L.T
```

Originally, it should be `L @ z` where `z` is a multivariate normal vector. Since each row of `eps` is such a vector, `eps @ L.T` is correct.

You get the desired result as below: