Generating Correlated Normals From Independent Ones
You have a random number generator that produces independent standard normal draws $z_1, z_2, \ldots$
(a) How do you produce two standard normal random variables $x_1, x_2$ (each $N(0,1)$) with correlation $\rho$? Prove that your construction has the required means, variances and correlation.
(b) Generalize: given a mean vector $\mu \in \mathbb{R}^{n}$ and a covariance matrix $\Sigma \in \mathbb{R}^{n \times n}$ (symmetric positive definite), how do you generate $X \sim N(\mu, \Sigma)$ from a vector $Z$ of $n$ independent standard normals? Name the matrix decomposition involved.
Hints
- Linear combinations of independent normals are normal. Try $x_1 = z_1$ and $x_2 = a z_1 + b z_2$; choose $a, b$ so that $\operatorname{Var}(x_2) = 1$ and $\operatorname{cov}(x_1, x_2) = \rho$.
- $\operatorname{cov}(z_1, a z_1 + b z_2) = a$, so $a = \rho$; then $\operatorname{Var}(x_2) = a^2 + b^2 = 1$ forces $b = \sqrt{1 - \rho^2}$.
- In $n$ dimensions, if $X = \mu + LZ$ then $\operatorname{Cov}(X) = LL^{T}$. So you need any $L$ with $LL^{T} = \Sigma$; the Cholesky factorization gives a lower-triangular one, and the $2 \times 2$ case reproduces part (a).
Worked Solution
How to Think About It: Linear maps of Gaussian vectors are Gaussian, with covariance transformed as $\operatorname{Cov}(LZ) = L\,\operatorname{Cov}(Z)\,L^{T} = LL^{T}$. So generating $N(\mu, \Sigma)$ reduces to finding a square root $L$ of $\Sigma$; in two dimensions you can find it by hand, and in general Cholesky does it.
Approach: Build $x_2$ as a mix of $z_1$ and an independent $z_2$ with weights chosen to hit variance $1$ and covariance $\rho$; then observe that these weights are the second row of the Cholesky factor of the $2 \times 2$ correlation matrix and state the $n$-dimensional recipe.
Formal Solution:
Part (a): Two correlated standard normals
*Step 1 -- Construction.* Let $$x_1 = z_1, \qquad x_2 = \rho\,z_1 + \sqrt{1 - \rho^{2}}\,z_2 .$$
*Step 2 -- Distribution.* Each $x_i$ is a linear combination of independent normals, hence normal, and $(x_1, x_2)$ is jointly normal (any linear combination of them is a linear combination of $z_1, z_2$). Means: $E[x_1] = E[x_2] = 0$.
*Step 3 -- Variances.* $\operatorname{Var}(x_1) = 1$ and, by independence of $z_1, z_2$, $$\operatorname{Var}(x_2) = \rho^{2}\operatorname{Var}(z_1) + (1 - \rho^{2})\operatorname{Var}(z_2) = \rho^{2} + 1 - \rho^{2} = 1.$$
*Step 4 -- Correlation.* $$\operatorname{cov}(x_1, x_2) = \operatorname{cov}\!\left(z_1,\ \rho z_1 + \sqrt{1-\rho^2}\,z_2\right) = \rho\operatorname{Var}(z_1) + \sqrt{1-\rho^2}\operatorname{cov}(z_1, z_2) = \rho,$$ and since both variances are $1$, $\operatorname{corr}(x_1, x_2) = \rho$. The construction works for every $\rho \in [-1, 1]$.
Part (b): $n$ dimensions
*Step 5 -- The general recipe.* Let $Z = (z_1, \ldots, z_n)^{T}$ with $\operatorname{Cov}(Z) = I$. For any matrix $L$ with $LL^{T} = \Sigma$, set $$X = \mu + LZ.$$ Then $X$ is Gaussian with $E[X] = \mu$ and $\operatorname{Cov}(X) = L\,\operatorname{Cov}(Z)\,L^{T} = LL^{T} = \Sigma$.
*Step 6 -- Cholesky decomposition.* Since $\Sigma$ is symmetric positive definite it has a unique Cholesky factorization $\Sigma = LL^{T}$ with $L$ lower triangular and positive diagonal (equivalently $\Sigma = R^{T}R$ with $R = L^{T}$ upper triangular). Triangularity means $x_k$ depends only on $z_1, \ldots, z_k$: each new coordinate is its regression on the previous ones plus fresh noise. For the $2 \times 2$ correlation matrix $\begin{pmatrix} 1 & \rho \\ \rho & 1\end{pmatrix}$, $$L = \begin{pmatrix} 1 & 0 \\ \rho & \sqrt{1 - \rho^{2}} \end{pmatrix},$$ which is exactly the construction of part (a). (Any other square root, e.g. the symmetric $\Sigma^{1/2}$ from the eigendecomposition, also works; Cholesky is the cheapest at $O(n^3/3)$ flops.) If $\Sigma$ is only positive semidefinite, use a pivoted Cholesky or the eigendecomposition $\Sigma = V\Lambda V^{T}$ with $L = V\Lambda^{1/2}$.
```python import numpy as np
def correlated_pair(rho, size, rng=np.random.default_rng()): z1, z2 = rng.standard_normal(size), rng.standard_normal(size) return z1, rho * z1 + np.sqrt(1 - rho**2) * z2
def mvn(mu, Sigma, size, rng=np.random.default_rng()): L = np.linalg.cholesky(Sigma) # Sigma = L @ L.T, L lower triangular Z = rng.standard_normal((size, len(mu))) return mu + Z @ L.T # each row is mu + L z ```
Answer: (a) $x_1 = z_1$, $x_2 = \rho z_1 + \sqrt{1 - \rho^{2}}\,z_2$; both are $N(0,1)$ and $\operatorname{corr}(x_1, x_2) = \rho$. (b) Compute the Cholesky factor $L$ with $LL^{T} = \Sigma$ and set $X = \mu + LZ$; the two-dimensional case is $L = \begin{pmatrix} 1 & 0 \\ \rho & \sqrt{1-\rho^2}\end{pmatrix}$.
Intuition
A covariance matrix is a Gram matrix $\Sigma = LL^{T}$, and applying $L$ to white noise paints exactly that covariance onto it; Cholesky is simply the triangular way of choosing $L$, which is why the second variable is built as (a piece of the first) plus (fresh independent noise scaled to keep unit variance). This is the standard engine behind Monte Carlo pricing of basket and spread options, multi-asset risk simulation, and correlated Brownian motions $dW_2 = \rho\,dW_1 + \sqrt{1-\rho^2}\,dW_1^{\perp}$ in stochastic volatility models.