5.5.3 The multivariate normal distribution

In this lesson

Jointly normal variables, built as $\mathbf{X} = \boldsymbol{\mu} + L\mathbf{Z}$ from independent standard normals
The rule $\operatorname{Cov}(A\mathbf{X} + \mathbf{b}) = A\Sigma A^\top$ , so $\boldsymbol{\mu} + L\mathbf{Z}$ has covariance matrix $LL^\top$
The density written with $\Sigma^{-1}$ and $\det \Sigma$ , and its level curves, which are ellipses
Every linear combination of jointly normal variables is normal
Uncorrelated jointly normal variables are independent, and a pair that is normal one at a time but not jointly
Simulation with multivariate_normal

You need covariance, the covariance matrix and the variance rules from lesson 5.5.2, and independence and joint densities from lesson 5.5.1. You need the normal distribution and sums of independent normals from lesson 5.4.3. You also need quadratic forms and their ellipses from lesson 3.7.1, and the Cholesky factor from lesson 3.7.2.

Built from independent standard normals

Start with independent standard normal variables $Z_1, \ldots, Z_n$ , collected in a random vector $\mathbf{Z}$ . Their joint density is $\varphi(z_1)\cdots\varphi(z_n)$ , the product of lesson 5.5.1. Now mix them with a matrix $L$ and shift them by a vector $\boldsymbol{\mu}$ .

Rule: Multivariate normal

\begin{gathered} \mathbf{X} = \boldsymbol{\mu} + L\mathbf{Z} \\[4pt] \mathbf{X} \sim \mathcal{N}(\boldsymbol{\mu}, \Sigma), \quad \Sigma = LL^\top \end{gathered}

$\mathbf{X}$ is jointly normal, or multivariate normal, with mean vector $\boldsymbol{\mu}$ and covariance matrix $\Sigma$ . Here $L$ is an $n \times n$ matrix with $\det L \neq 0$ , and $\mathbf{Z}$ holds independent standard normals.

Take $\boldsymbol{\mu} = (1, 3)$ and $L = \begin{pmatrix} 2 & 0 \\ 1 & 2 \end{pmatrix}$ , so that

\begin{gathered} X_1 = 1 + 2Z_1 \\ X_2 = 3 + Z_1 + 2Z_2 \end{gathered}

Each is normal on its own. $X_1$ is a linear change of $Z_1$ , so $X_1 \sim \mathcal{N}(1, 4)$ (lesson 5.4.3). $X_2$ is 3 plus a sum of independent normals, so $X_2 \sim \mathcal{N}(3, 1 + 4) = \mathcal{N}(3, 5)$ .

They share $Z_1$ , so they are linked. By the bilinearity of lesson 5.5.2, $\operatorname{Cov}(X_1, X_2) = \operatorname{Cov}(2Z_1,\ Z_1 + 2Z_2) = 2\operatorname{Cov}(Z_1, Z_1) = 2$ . So

\begin{gathered} \Sigma = \begin{pmatrix} 4 & 2 \\ 2 & 5 \end{pmatrix} \\[6pt] \rho = \frac{2}{\sqrt{4 \cdot 5}} = 0.4472 \end{gathered}

Case study 3.7.4 simulated the three assets' returns in exactly this way, with the Cholesky factor of their covariance matrix as $L$ .

The covariance rule for random vectors

Two rules describe any linear change $A\mathbf{X} + \mathbf{b}$ of a random vector, with $A$ a matrix and $\mathbf{b}$ a vector. The mean moves by linearity, entry by entry (lesson 5.3.2). The covariance matrix follows from bilinearity: entry $(i, j)$ of $A\Sigma A^\top$ is the covariance of entries $i$ and $j$ of $A\mathbf{X}$ .

Rule: Linear change of a random vector

\begin{gathered} \mathbb{E}[A\mathbf{X} + \mathbf{b}] = A\boldsymbol{\mu} + \mathbf{b} \\[4pt] \operatorname{Cov}(A\mathbf{X} + \mathbf{b}) = A\,\Sigma\,A^\top \end{gathered}

$\operatorname{Cov}(\mathbf{X})$ means the covariance matrix of $\mathbf{X}$ , here $\Sigma$ . The shift $\mathbf{b}$ changes the mean only. The rule holds for any random vector, normal or not.

Independent standard normals have mean $\mathbf{0}$ and covariance matrix $I$ , the identity. So $\boldsymbol{\mu} + L\mathbf{Z}$ has mean $\boldsymbol{\mu}$ and covariance matrix $LIL^\top = LL^\top$ . Case study 3.7.4 found $S_X = LS_ZL^\top$ for samples; this is the same rule for random variables, with $S_Z \approx I$ replaced by $I$ exactly.

For the example, $LL^\top = \begin{pmatrix} 4 & 2 \\ 2 & 5 \end{pmatrix}$ , as found term by term above. For the three assets, case study 3.7.4's $L$ gives back case study 3.2.5's $\Sigma$ .

Only $\boldsymbol{\mu}$ and $\Sigma$ matter: two matrices $L$ with the same $LL^\top$ give the same distribution. Given $\Sigma$ , the usual choice is its Cholesky factor (lesson 3.7.2); the symmetric square root of lesson 3.7.2 is another. Case study 5.5.5 uses the rule with a single row, $A = \mathbf{w}^\top$ .

LL^\top

import numpy as np mu = np.array([1.0, 3.0]) L = np.array([[2.0, 0.0], [1.0, 2.0]]) Sigma = L @ L.T print(Sigma) Z = np.random.default_rng(555).standard_normal((100_000, 2)) X = mu + Z @ L.T print(np.round(X.mean(axis=0), 3), np.round(np.cov(X, rowvar=False), 3)) lam, Q = np.linalg.eigh(Sigma) root = Q @ np.diag(np.sqrt(lam)) @ Q.T # the symmetric square root print(np.round(root, 4), np.allclose(root @ root.T, Sigma))

⌘/Ctrl + Enter runs the cell

The density

Case study 4.7.5 found the density of $\boldsymbol{\mu} + L\mathbf{Z}$ for two variables: $\frac{\varphi(z_1)\,\varphi(z_2)}{\det L}$ , with $\mathbf{z} = L^{-1}(\mathbf{x} - \boldsymbol{\mu})$ . Two facts rewrite it with $\Sigma$ alone.

The exponent. $\varphi(z_1)\varphi(z_2) = \frac{1}{2\pi}e^{-\mathbf{z}^\top\mathbf{z}/2}$ . Since $\Sigma^{-1} = (LL^\top)^{-1} = (L^\top)^{-1}L^{-1}$ (lesson 3.2.3), $\mathbf{z}^\top\mathbf{z} = (\mathbf{x} - \boldsymbol{\mu})^\top\Sigma^{-1}(\mathbf{x} - \boldsymbol{\mu})$ .
The constant. $\det \Sigma = (\det L)^2$ (lesson 3.7.2), so $\det L = \sqrt{\det \Sigma}$ .

Rule: Multivariate normal density

\begin{gathered} f(\mathbf{x}) = \frac{e^{-Q/2}}{(2\pi)^{n/2}\sqrt{\det \Sigma}} \\[6pt] Q = (\mathbf{x} - \boldsymbol{\mu})^\top\Sigma^{-1}(\mathbf{x} - \boldsymbol{\mu}) \end{gathered}

$Q$ measures how far $\mathbf{x}$ is from $\boldsymbol{\mu}$ in units set by $\Sigma$ . With $n = 1$ and $\Sigma = \sigma^2$ , $Q = \frac{(x - \mu)^2}{\sigma^2}$ and this is lesson 5.4.3's normal density.

For two variables, write $u = \frac{x - \mu_X}{\sigma_X}$ and $v = \frac{y - \mu_Y}{\sigma_Y}$ for the standardised values. Then $\det \Sigma = \sigma_X^2\sigma_Y^2(1 - \rho^2)$ , and

\begin{gathered} f(x, y) = \frac{e^{-Q/2}}{2\pi\sigma_X\sigma_Y\sqrt{1 - \rho^2}} \\[6pt] Q = \frac{u^2 - 2\rho uv + v^2}{1 - \rho^2} \end{gathered}

The example. $\det \Sigma = 4 \cdot 5 - 2 \cdot 2 = 16 = (\det L)^2$ , and $\Sigma^{-1} = \frac{1}{16}\begin{pmatrix} 5 & -2 \\ -2 & 4 \end{pmatrix}$ . At $\mathbf{x} = (3, 3)$ , $\mathbf{x} - \boldsymbol{\mu} = (2, 0)$ and $Q = \frac{4 \cdot 5}{16} = 1.25$ . The other route gives $\mathbf{z} = (1, -0.5)$ , and $\mathbf{z}^\top\mathbf{z} = 1.25$ too.

So $f(3, 3) = \frac{e^{-0.625}}{8\pi} = 0.021297$ . The peak, at $\boldsymbol{\mu}$ , is $\frac{1}{8\pi} = 0.039789$ .

Case study 4.7.5's returns. For A and B, $\sqrt{\det \Sigma} = 0.020605$ , the $\det L$ of that case study, and the peak $\frac{1}{2\pi\sqrt{\det \Sigma}}$ is its 7.7241. At $(0, -0.1)$ the density is 4.583039, the value of case study 4.7.5's formula there.

\Sigma^{-1}(\mathbf{x} - \boldsymbol{\mu})

import math import numpy as np from scipy.stats import multivariate_normal, norm def density(x, mu, Sigma): d = np.asarray(x, float) - mu Q = d @ np.linalg.solve(Sigma, d) return math.exp(-Q / 2) / ((2 * math.pi) ** (len(d) / 2) * math.sqrt(np.linalg.det(Sigma))) mu = np.array([1.0, 3.0]) L = np.array([[2.0, 0.0], [1.0, 2.0]]) Sigma = L @ L.T z = np.linalg.solve(L, np.array([3.0, 3.0]) - mu) print(z, round(density([3, 3], mu, Sigma), 6), round(norm.pdf(z[0]) * norm.pdf(z[1]) / np.linalg.det(L), 6), round(multivariate_normal(mu, Sigma).pdf([3, 3]), 6)) print(round(density(mu, mu, Sigma), 6), round(1 / (8 * math.pi), 6)) mu_ab = np.array([0.06, 0.08]) S_ab = np.array([[0.0144, 0.3 * 0.12 * 0.18], [0.3 * 0.12 * 0.18, 0.0324]]) print(round(math.sqrt(np.linalg.det(S_ab)), 6), round(density(mu_ab, mu_ab, S_ab), 4), round(density([0.0, -0.1], mu_ab, S_ab), 6))

⌘/Ctrl + Enter runs the cell

Level curves are ellipses

The density depends on $\mathbf{x}$ only through $Q$ , so it is constant where $Q = c$ . $\Sigma^{-1}$ is positive definite, so these level curves are ellipses centred at $\boldsymbol{\mu}$ (lesson 3.7.1).

Their axes lie along the eigenvectors of $\Sigma$ , which are also those of $\Sigma^{-1}$ . An eigenvalue $\lambda_i$ of $\Sigma$ is $\frac{1}{\lambda_i}$ for $\Sigma^{-1}$ , so lesson 3.7.1's semi-axis rule gives $\sqrt{c\,\lambda_i}$ .

Rule: Ellipses of the normal

\begin{gathered} Q = c: \text{ semi-axes } \sqrt{c\,\lambda_i} \\[4pt] n = 2:\ P(Q \le c) = 1 - e^{-c/2} \end{gathered}

$\lambda_i$ are the eigenvalues of $\Sigma$ , and the axes point along its eigenvectors. With two variables, the ellipse $Q = c$ holds the probability $1 - e^{-c/2}$ , whatever $\boldsymbol{\mu}$ and $\Sigma$ are.

The probability comes from case study 4.7.5. There $\mathbf{z}$ landed within distance $r$ of the origin with probability $1 - e^{-r^2/2}$ , and $Q = \mathbf{z}^\top\mathbf{z}$ is that distance squared.

For the example, $\Sigma$ has eigenvalues $\frac{9 \mp \sqrt{17}}{2}$ : 2.4384 along $(0.7882, -0.6154)$ and 6.5616 along $(0.6154, 0.7882)$ . The long axis rises at 52.02° to the $x$ -axis. The ellipse $Q = 1$ has semi-axes 1.5616 and 2.5616, and $Q = 4$ twice those.

ellipse	$Q = 1$	$Q = 4$	$Q = 9$
probability inside, $1 - e^{-c/2}$	0.3935	0.8647	0.9889
share of 100,000 draws	0.3942	0.8646	0.9891

Draws of $\mathbf{X} = \boldsymbol{\mu} + L\mathbf{Z}$ and the ellipses $Q = 1$ , 4 and 9 of its density. The axes lie along the eigenvectors of $\Sigma$ , with semi-axes $\sqrt{c\lambda_i}$ ; correlation tilts them away from the coordinate axes.

With $\rho = 0$ , $\Sigma$ is diagonal and the axes are the coordinate axes. The ellipses of case study 4.7.5's figure tilted only once the correlation was 0.3.

\Sigma

import math import numpy as np mu = np.array([1.0, 3.0]) L = np.array([[2.0, 0.0], [1.0, 2.0]]) Sigma = L @ L.T lam, vecs = np.linalg.eigh(Sigma) vecs = vecs * np.sign(vecs[0]) # first entry positive print(np.round(lam, 4), np.round(vecs, 4), round(math.degrees(math.atan2(vecs[1, 1], vecs[0, 1])), 2)) X = mu + np.random.default_rng(555).standard_normal((100_000, 2)) @ L.T D = X - mu Q = np.sum(D * np.linalg.solve(Sigma, D.T).T, axis=1) for c in (1, 4, 9): print(c, np.round(np.sqrt(c * lam), 4), round(1 - math.exp(-c / 2), 4), round(np.mean(Q <= c), 4))

⌘/Ctrl + Enter runs the cell

Linear combinations are normal

Any combination $\mathbf{a}^\top\mathbf{X} = a_1X_1 + \cdots + a_nX_n$ of jointly normal variables is $\mathbf{a}^\top\boldsymbol{\mu} + (\mathbf{a}^\top L)\mathbf{Z}$ . That is a constant plus a combination of independent normals, so it is normal (lesson 5.4.3). Its mean and variance come from the rule above with the single row $A = \mathbf{a}^\top$ .

Rule: Combinations of jointly normal variables

\begin{gathered} \mathbf{a}^\top\mathbf{X} \sim \mathcal{N}\big(\mathbf{a}^\top\boldsymbol{\mu},\ \mathbf{a}^\top\Sigma\,\mathbf{a}\big) \\[4pt] A\mathbf{X} + \mathbf{b} \sim \mathcal{N}\big(A\boldsymbol{\mu} + \mathbf{b},\ A\Sigma A^\top\big) \end{gathered}

Sums of normal variables are normal when they are independent, or more generally jointly normal. Several combinations at once, $A\mathbf{X} + \mathbf{b}$ , are jointly normal too.

In the example, $X_1 + X_2 = 4 + 3Z_1 + 2Z_2$ , so $X_1 + X_2 \sim \mathcal{N}(4, 13)$ . Lesson 5.5.2's rule agrees: $4 + 5 + 2 \cdot 2 = 13$ . Likewise $X_1 - X_2 \sim \mathcal{N}(-2, 4 + 5 - 4) = \mathcal{N}(-2, 5)$ .

A probability then needs only lesson 5.4.3. $P(X_1 + X_2 > 10) = 1 - N\big(\frac{10 - 4}{\sqrt{13}}\big) = 1 - N(1.6641) = 0.04805$ , with $N$ the standard normal distribution function. The seeded draws give 0.04805 too.

Case study 4.7.5 used this rule before it was stated. Half in A and half in B is normal, with mean 7% and volatility 12.2229%. So it loses more than 10% with probability 0.082139, the same as that case study's double integral.

\mathbf{a}^\top\mathbf{X}

import math import numpy as np from scipy.stats import norm mu = np.array([1.0, 3.0]) L = np.array([[2.0, 0.0], [1.0, 2.0]]) Sigma = L @ L.T for a in (np.array([1, 1]), np.array([1, -1])): print(a, a @ L, a @ mu, a @ Sigma @ a) X = mu + np.random.default_rng(555).standard_normal((100_000, 2)) @ L.T s = X[:, 0] + X[:, 1] print(round(6 / math.sqrt(13), 4), round(1 - norm.cdf(6 / math.sqrt(13)), 5), round(np.mean(s > 10), 5), round(np.var(s, ddof=1), 3))

⌘/Ctrl + Enter runs the cell

Uncorrelated jointly normal variables are independent

When every correlation is 0, $\Sigma$ is diagonal, with the variances $\sigma_i^2$ on it. Then $Q = \sum_i \frac{(x_i - \mu_i)^2}{\sigma_i^2}$ and $\det \Sigma = \sigma_1^2 \cdots \sigma_n^2$ . The exponential of a sum is a product, so the density splits into $n$ one-variable normal densities: by lesson 5.5.1, the variables are independent.

For example, with $\boldsymbol{\mu} = (1, 0)$ and $\Sigma$ diagonal with entries 4 and 9, the density at $(2, -1)$ is 0.022144. That is the $\mathcal{N}(1, 4)$ density at 2 times the $\mathcal{N}(0, 9)$ density at $-1$ .

Normal one at a time is not enough. Let $Z$ be standard normal, and $S$ an independent fair coin flip, 1 or $-1$ . Set $Y = SZ$ .

$Y$ is standard normal, since flipping the sign of a symmetric variable at random changes nothing.
$\operatorname{Cov}(Z, Y) = \mathbb{E}[SZ^2] = \mathbb{E}[S]\,\mathbb{E}[Z^2] = 0$ : uncorrelated.
Yet $\lvert Y \rvert = \lvert Z \rvert$ always, so they are dependent.

The pair $(Z, Y)$ is not jointly normal. Whenever $S = -1$ , $Z + Y = 0$ , so $P(Z + Y = 0) = \frac{1}{2}$ , and no normal variable takes one value with probability $\frac{1}{2}$ .

Two pairs with standard normal coordinates and zero correlation. Only the first is jointly normal, and only there are the coordinates independent.

Rule: Zero correlation and independence

\begin{gathered} \text{jointly normal and } \Sigma \text{ diagonal} \\ \implies \text{independent} \end{gathered}

For jointly normal variables, uncorrelated and independent mean the same. Normal one at a time is not enough: $Z$ and $SZ$ are both normal and uncorrelated, but dependent.

Z

import numpy as np from scipy.stats import multivariate_normal, norm print(round(multivariate_normal([1, 0], np.diag([4.0, 9.0])).pdf([2, -1]), 6), round(norm.pdf(2, 1, 2) * norm.pdf(-1, 0, 3), 6)) rng = np.random.default_rng(556) z = rng.standard_normal(100_000) s = rng.choice([-1, 1], size=100_000) y = s * z print(round(y.mean(), 3) + 0.0, round(y.std(ddof=1), 3), round(np.mean(np.abs(y) < 1), 4)) print(round(np.corrcoef(z, y)[0, 1], 4), round(np.mean(z + y == 0), 4), round(np.corrcoef(np.abs(z), np.abs(y))[0, 1], 4))

⌘/Ctrl + Enter runs the cell

In Python

SciPy's multivariate_normal(mean, cov) gives the density with .pdf and draws with .rvs(size, random_state=rng). NumPy's rng.multivariate_normal(mean, cov, size, method="cholesky") draws too; case study 3.2.5 used it. Both take the covariance matrix, not the volatilities or the correlation matrix.

SciPy's .cdf finds probabilities by a randomised numerical method, so its last digits can change from run to run. For the probability of a region, use dblquad (case study 4.7.5) or a seeded simulation.

The covariance matrix must be valid: symmetric and positive semidefinite (case study 3.7.4). Repair a matrix that is not before drawing from it.

Try it The cell draws 100,000 years of the three assets of case study 3.2.5 with rvs and a fixed seed. In this run the sample means are 0.0597, 0.0788 and 0.1084, the volatilities 0.1203, 0.1804 and 0.2496, and the correlations 0.300, 0.100 and 0.499.

import numpy as np from scipy.stats import multivariate_normal mu = np.array([0.06, 0.08, 0.11]) sigma = np.array([0.12, 0.18, 0.25]) R = np.array([[1.0, 0.3, 0.1], [0.3, 1.0, 0.5], [0.1, 0.5, 1.0]]) Sigma = np.diag(sigma) @ R @ np.diag(sigma) draws = multivariate_normal(mu, Sigma).rvs(size=100_000, random_state=np.random.default_rng(557)) C = np.corrcoef(draws, rowvar=False) print(np.round(draws.mean(axis=0), 4), np.round(draws.std(axis=0, ddof=1), 4)) print(np.round([C[0, 1], C[0, 2], C[1, 2]], 3))

⌘/Ctrl + Enter runs the cell

Case study 5.5.5 uses this distribution to ask what one asset's return says about another's.

Common mistakes

Normal one at a time is not jointly normal. For $Z$ and $SZ$ , zero correlation does not give independence, and $Z + SZ$ is not normal.
The density uses $\Sigma^{-1}$ and $\sqrt{\det \Sigma}$ . In two variables the constant is $2\pi\sqrt{\det \Sigma}$ , here $8\pi$ , not $2\pi \det \Sigma$ .
Combinations need the covariances. $\operatorname{Var}(X_1 + X_2) = 13$ , not $4 + 5 = 9$ .
cov means the covariance matrix. Passing volatilities, standard deviations or the correlation matrix gives the wrong distribution, often without an error.
The ellipses tilt. Their axes follow the eigenvectors of $\Sigma$ , not the coordinate axes, unless $\rho = 0$ .
Check $\Sigma$ first. A matrix that is not positive semidefinite is not a covariance matrix (case study 3.7.4).

The multivariate normal distribution

0 of 12 correct

1.
$X_1 = 3Z_1$ and $X_2 = Z_1 + Z_2$ , with $Z_1, Z_2$ independent standard normals. Find $\operatorname{Var}(X_2)$
2.
For the same pair, find $\operatorname{Cov}(X_1, X_2)$
3.
$L = \begin{pmatrix} 3 & 0 \\ 1 & 1 \end{pmatrix}$ . Find the entry $\Sigma_{11}$ of $\Sigma = LL^\top$
4.
A bivariate normal has $\det L = 3$ . Find $\det \Sigma$
5.
A bivariate normal has $\det \Sigma = 16$ . Find its density at the mean (type pi for $\pi$ )
6.
$\boldsymbol{\mu} = \mathbf{0}$ and $\Sigma$ is diagonal with entries 4 and 9. Find $Q$ at $\mathbf{x} = (2, 3)$
7.
A bivariate normal: what is the probability that $Q \le 4$ ? Type an exact answer
8.
$\Sigma$ has eigenvalue 9. How long is the semi-axis of the ellipse $Q = 1$ along its eigenvector?
9.
$X_1$ and $X_2$ are jointly normal with variances 4 and 9 and covariance 3. Find $\operatorname{Var}(X_1 + X_2)$
10.
$X_1 + X_2 \sim \mathcal{N}(4, 16)$ . How many standard deviations above its mean is 10?
11.
$X$ and $Y$ are jointly normal, each with mean 0, and uncorrelated. Find $P(X > 0, Y > 0)$
12.
$Z$ is standard normal and $S$ an independent fair sign, $\pm 1$ . Find $P(Z + SZ = 0)$

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

\mathbf{X} = \boldsymbol{\mu} + L\mathbf{Z}

The multivariate normal distribution

Built from independent standard normals

Rule: Multivariate normal

The covariance rule for random vectors

Rule: Linear change of a random vector

The density

Rule: Multivariate normal density

Level curves are ellipses

Rule: Ellipses of the normal

Linear combinations are normal

Rule: Combinations of jointly normal variables

Uncorrelated jointly normal variables are independent

Rule: Zero correlation and independence

In Python

Common mistakes

The multivariate normal distribution

Summary card

Free further reading