Solution: Cholesky Factorisation and Simulating Correlated Gaussians

Exercise: Cholesky Factorisation and Simulating Correlated Gaussians

Part 1 — Cholesky by hand

Write $L = \begin{pmatrix}l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33}\end{pmatrix}$ and expand $LL^\top$ :

$l_{11}^2 = 0.04 \Rightarrow l_{11} = 0.2$ .
$l_{11}l_{21} = 0.02 \Rightarrow l_{21} = 0.1$ .
$l_{11}l_{31} = 0.01 \Rightarrow l_{31} = 0.05$ .
$l_{21}^2 + l_{22}^2 = 0.09 \Rightarrow l_{22}^2 = 0.08 \Rightarrow l_{22} = \sqrt{0.08} \approx 0.2828$ .
$l_{21}l_{31} + l_{22}l_{32} = 0.03 \Rightarrow 0.005 + 0.2828\cdot l_{32} = 0.03 \Rightarrow l_{32} = 0.0884$ .
$l_{31}^2 + l_{32}^2 + l_{33}^2 = 0.16 \Rightarrow 0.0025 + 0.00781 + l_{33}^2 = 0.16 \Rightarrow l_{33}^2 = 0.14969 \Rightarrow l_{33} = 0.3869$ .

L \approx \begin{pmatrix}0.2 & 0 & 0 \\ 0.1 & 0.2828 & 0 \\ 0.05 & 0.0884 & 0.3869\end{pmatrix}.

Verification by matrix multiplication recovers $\Sigma$ .

Part 2

If $Z \sim \mathcal{N}(0, I)$ and $X = \mu + LZ$ , then $X$ is an affine transformation of a gaussian, hence gaussian. Mean: $\mathbb{E}[X] = \mu + L\cdot 0 = \mu$ . Covariance: $\text{Cov}(X) = L\cdot \text{Cov}(Z)\cdot L^\top = L\cdot I\cdot L^\top = LL^\top = \Sigma$ . So $X \sim \mathcal{N}(\mu, \Sigma)$ .

Part 3 — Simulation

import numpy as np

mu = np.array([0.05, 0.08, 0.10])
Sigma = np.array([[0.04, 0.02, 0.01],
                  [0.02, 0.09, 0.03],
                  [0.01, 0.03, 0.16]])

L = np.linalg.cholesky(Sigma)
rng = np.random.default_rng(0)
N = 10_000

Z = rng.standard_normal((N, 3))
X = mu + Z @ L.T

print("sample mean:", X.mean(axis=0).round(4))
# sample mean: [0.0499 0.0801 0.1006]

print("sample cov:\n", np.cov(X.T).round(4))
# sample cov:
#  [[0.0399 0.0203 0.0098]
#   [0.0203 0.0906 0.0308]
#   [0.0098 0.0308 0.1576]]

Within sampling noise, both match $\mu$ and $\Sigma$ .

Part 4 — Portfolio variance

w = np.array([1, 1, 1]) / 3

# theoretical
var_theoretical = w @ Sigma @ w
print(var_theoretical)
# 0.04888888888888889

# empirical
port_returns = X @ w
var_empirical = port_returns.var(ddof=1)
print(var_empirical)
# 0.04936...

print("relative diff:", abs(var_empirical/var_theoretical - 1))
# ~1%

Theoretical: $w^\top\Sigma w = (1/3)^2\cdot (0.04 + 0.09 + 0.16 + 2\cdot 0.02 + 2\cdot 0.01 + 2\cdot 0.03) = (1/9)\cdot 0.44 = 0.0489$ . Empirical agrees to within 1%, as expected from $N = 10{,}000$ samples.

Takeaways

Cholesky factor is the computational engine for gaussian sampling in multi-dimensional Monte Carlo. Cheaper than eigen-decomposition when $\Sigma$ is positive definite.
Affine transformations preserve gaussianity. $X = \mu + LZ$ is gaussian with the desired mean and covariance — the most useful identity in simulation.
Portfolio variance $w^\top\Sigma w$ is the scalar summary that matters for mean-variance and VaR.
When $\Sigma$ is only PSD (not PD), Cholesky fails. Cures: pseudo-Cholesky (LDL decomposition), eigen-decomposition $X = \mu + Q\Lambda^{1/2}Z$ , or regularisation $\Sigma \gets \Sigma + \epsilon I$ .