Case study: valid correlation matrices and correlated simulation

What a correlation matrix must be

Case study 3.2.5 built the covariance matrix of three assets as $\Sigma = DRD$. Here $D$ is the diagonal matrix of the volatilities 0.12, 0.18 and 0.25, and $R$ is the correlation matrix. The data are invented for illustration.

A variance $\mathbf{w}^\top\Sigma\mathbf{w}$ is never negative, so $\Sigma$ must be positive semidefinite (lesson 3.7.1). And $\mathbf{w}^\top\Sigma\mathbf{w} = (D\mathbf{w})^\top R\,(D\mathbf{w})$. With positive volatilities, every vector is $D\mathbf{w}$ for some $\mathbf{w}$, so $R$ must be positive semidefinite too.

The three assets pass, with $\rho_{AB} = 0.3$, $\rho_{AC} = 0.1$ and $\rho_{BC} = 0.5$. The eigenvalues of $R$ are 0.4571, 0.9124 and 1.6305 (lesson 3.7.1), its leading minors 1, 0.91 and 0.68, and Cholesky succeeds (lesson 3.7.2). So $R$ is positive definite: no portfolio of A, B and C is free of risk.

Three correlations that cannot hold together

Now the hook's correlations:

It is symmetric with 1 on its diagonal, and every entry is between $-1$ and 1. It fails the test that matters:

• Leading minors: 1, 0.19 and $-0.984$. Case study 3.4.3 noted that $\det R$ lies between 0 and 1; this one is negative.
• Eigenvalues: $-0.3767$, 1.2 and 2.1767. One is negative, so $R_{\text{bad}}$ is indefinite.
• Cholesky: np.linalg.cholesky raises LinAlgError: Matrix is not positive definite.

The eigenvector of $-0.3767$ is about $(-0.5192, 0.6788, -0.5192)$: long B against A and C. The round vector $\mathbf{x} = (1, -1, 1)$ points nearly the opposite way, which gives the same sign: $\mathbf{x}^\top R_{\text{bad}}\mathbf{x} = 3 - 4(0.9) + 2(-0.2) = -1$. With the volatilities, the weights $(1, -1, 1)$ give $\mathbf{w}^\top DR_{\text{bad}}D\,\mathbf{w} = -0.02258$, the hook's variance.

Why: correlations are angles

Case study 3.1.4 found that a correlation is the cosine of the angle between two centred return vectors. A correlation of 0.9 is an angle of $25.84°$.

So A is within $25.84°$ of B, and B within $25.84°$ of C. Then A and C are at most $51.68°$ apart, and their correlation is at least $\cos 51.68° = 0.62$. A correlation of $-0.2$ would need an angle of $101.54°$.

Keep $\rho_{AB} = \rho_{BC} = 0.9$ and move $\rho_{AC}$:

The matrix is valid exactly when $0.62 \le \rho_{AC} \le 1$. At either end it is singular: some combination of the three has zero variance.

Case study 3.2.5's desk note met a second impossible set, $\rho_{AB} = \rho_{AC} = 0.9$ with $\rho_{BC} = -0.9$. Its eigenvalues are $-0.8$, 1.9 and 1.9.

Where a desk meets this

A correlation matrix computed by np.corrcoef from one table of returns, with the same dates for every asset, is always valid. It is a matrix of the form $M^\top M$ (lesson 3.7.1), scaled. Trouble starts when the entries come from different places:

• pairwise estimates, each pair on the dates both assets traded, so every entry uses different data;
• hand-set values for a stress scenario, as in the hook;
• blocks from different sources stitched into one matrix, such as equities from one model and rates from another.

A simple repair

A desk cannot simulate or optimise with an impossible matrix, so it repairs it. The simplest repair uses the spectral theorem of lesson 3.6.3, $R = Q\Lambda Q^\top$:

1. 1

   Find the eigenvalues and eigenvectors: $R_{\text{bad}} = Q\Lambda Q^\top$.

   Why?
2. 2

   Replace each negative eigenvalue by 0, giving $\Lambda_+$, and rebuild $C = Q\Lambda_+Q^\top$.

   Why?
3. 3

   Divide each entry $c_{ij}$ by $\sqrt{c_{ii}\,c_{jj}}$, to put 1 back on the diagonal.

   Why?

Step 2 changes the diagonal, here to 1.1016, 1.1736 and 1.1016; step 3 puts the 1s back, as case study 3.4.3 turned $\Sigma$ into $R$. The result:

Its eigenvalues are 0, 1.0894 and 1.9106: valid. The portfolio $(1, -1, 1)$ now has variance 0.014058, a volatility of 11.86%.

The effect. The repair moved every correlation, not only the doubtful one. The two estimates of 0.9 fell to 0.6748, and $-0.2$ rose to $-0.0894$.

The repaired matrix is singular. Its smallest eigenvalue is 0, so one combination of the assets has zero variance. In money, that combination is the weights $(1.7232, -1.5504, 0.8272)$, which add up to 1. A portfolio optimiser looking for low risk would pile into this "riskless" portfolio.

On a computer, that 0 comes out as a tiny number of either sign, so np.linalg.cholesky may accept the matrix on one computer and reject it on another.

Use a floor instead of 0. Raise each eigenvalue below 0.001 to 0.001, then rescale. The correlations become 0.6742 and $-0.0891$, and the eigenvalues 0.0009, 1.0891 and 1.91. Cholesky now works, with last diagonal entry 0.058, and the "riskless" portfolio has a volatility of 1.2%: small, but no longer 0.

If only $\rho_{AC}$ is doubtful, a better repair moves it alone, into its valid range from 0.62 to 1. Iterative methods find the closest valid matrix; this repair is simple, not the closest.

Simulating correlated returns

With a valid matrix, a desk can simulate: draw many possible years of returns for A, B and C, and see what a portfolio does across them.

Independent draws

rng.standard_normal returns numbers whose histogram follows the bell curve $\varphi$ of lesson 2.6.3. About 68.3% of them land between $-1$ and 1. Their mean is close to 0, and their standard deviation, the square root of their variance, close to 1.

Draw a $10{,}000 \times 3$ table $Z$ of them with the seed 37, one row per draw. Each number is drawn separately, so the columns are independent: one tells you nothing about the others (module 5.5 makes this precise).

In this table the column means are $-0.004$, $-0.015$ and 0.000, the standard deviations 1.001, 0.993 and 1.017, and 67.93% of the numbers lie between $-1$ and 1. The sample correlations between columns are $-0.003$, 0.016 and 0.006. So the sample covariance matrix of $Z$ is close to the identity $I$.

Mixing them with L

Take $L$, the Cholesky factor of $\Sigma$ (lesson 3.7.2), and set $\mathbf{x} = L\mathbf{z}$ for each draw $\mathbf{z}$:

Read it row by row. A's return is $0.12z_1$. B's is $0.054z_1 + 0.1717z_2$: part of A's draw, plus a draw of its own.

C's return uses all three draws. $L$ is $D$ times the Cholesky factor of $R$: each row of $R$'s factor scaled by that asset's volatility.

Why the covariance comes out right. In the table, row $t$ of $X$ is $\mathbf{x}_t^\top = \mathbf{z}_t^\top L^\top$, so $X = ZL^\top$. Centring each column, as case study 3.2.5 did, keeps this: $X_c = Z_cL^\top$.

Then $X_c^\top X_c = LZ_c^\top Z_cL^\top$ by the rule $(AB)^\top = B^\top A^\top$ of lesson 3.2.3. Dividing by $n - 1$ gives, exactly, for any sample:

Since $S_Z \approx I$, $S_X \approx LL^\top = \Sigma$. Lesson 5.5.3 gives the same rule for random variables.

From the 10,000 draws, the sample correlations for AB, AC and BC are about 0.30, 0.11 and 0.50, against 0.3, 0.1 and 0.5. The sample volatilities are 12.02%, 17.87% and 25.38%, and the portfolio $(0.5, 0.3, 0.2)$ has a sample volatility of 11.99%, against the model's 11.94%.

How case study 3.2.5 made its sample

Case study 3.2.5 drew 120 monthly returns with rng.multivariate_normal(mu / 12, Sigma / 12, size=120, method="cholesky") and the seed 7. That is this method: monthly means $\boldsymbol{\mu}/12$, and the Cholesky factor of $\Sigma/12$, which is $L/\sqrt{12}$.

Rebuilt by hand, with $Z$ drawn as a $120 \times 3$ table, the portfolio's first three monthly returns are 0.0095, $-0.0369$ and 0.0311: case study 3.2.5's numbers. Without method="cholesky", NumPy takes another square root of $\Sigma$, from the singular value decomposition (module 3.8): same covariance, different numbers.

import numpy as np

MU = np.array([0.06, 0.08, 0.11])
VOLS = 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]])
BAD = np.array([[1.0, 0.9, -0.2], [0.9, 1.0, 0.9], [-0.2, 0.9, 1.0]])


def _expect_value_error(f, why):
    try:
        f()
    except ValueError:
        return
    except np.linalg.LinAlgError:
        raise AssertionError(f"expected ValueError, got LinAlgError: {why}")
    raise AssertionError(f"expected ValueError: {why}")


def test_valid_matrix_unchanged():
    """The three assets' R, already valid, comes back unchanged"""
    expect_close(repair_correlation(R), R, what="repair_correlation(R)", rel=1e-9, abs_tol=1e-9)


def test_repair_values():
    """The impossible matrix becomes 0.6742 for AB and BC and −0.0891 for AC (floor 0.001)"""
    F = np.asarray(repair_correlation(BAD), dtype=float)
    expect_close([F[0, 1], F[1, 2], F[0, 2]], [0.674246, 0.674246, -0.089102], what="the repaired correlations AB, BC, AC", rel=1e-4)


def test_repair_is_valid():
    """The repaired matrix is symmetric, has 1s on its diagonal and passes Cholesky"""
    F = np.asarray(repair_correlation(BAD), dtype=float)
    expect_close(F, F.T, what="the repaired matrix against its transpose")
    expect_close(np.diag(F), np.ones(3), what="the diagonal of the repaired matrix")
    assert np.linalg.eigvalsh(F).min() > 5e-4, f"the smallest eigenvalue should be about 0.0009, got {np.linalg.eigvalsh(F).min()}"
    np.linalg.cholesky(F)


def test_repair_floor_zero():
    """With floor=0 the negative eigenvalue is clipped to 0: 0.6748 for AB"""
    F = np.asarray(repair_correlation(BAD, floor=0.0), dtype=float)
    expect_close(F[0, 1], 0.674764, what="the repaired AB correlation with floor=0", rel=1e-4)
    expect_close(np.linalg.eigvalsh(F).min(), 0.0, what="the smallest eigenvalue with floor=0", abs_tol=1e-9)


def test_repair_rejects_bad_input():
    """A matrix that is not symmetric, or has 2 on its diagonal, raises ValueError"""
    _expect_value_error(lambda: repair_correlation([[1, 0.5], [0.4, 1]]), "the matrix is not symmetric")
    _expect_value_error(lambda: repair_correlation([[2, 0.5], [0.5, 1]]), "a diagonal entry is 2, not 1")


def test_simulate_shape_and_seed():
    """simulate_returns gives one row per draw and one column per asset, and repeats with the same seed"""
    X = np.asarray(simulate_returns(MU, VOLS, R, 500, 1), dtype=float)
    expect_equal(X.shape, (500, 3), what="the shape of simulate_returns(..., n=500)")
    expect_close(np.asarray(simulate_returns(MU, VOLS, R, 500, 1)), X, what="a second run with the same seed")


def test_simulate_matches_case_study():
    """Case study 3.2.5's sample: seed 7, 120 months, gives portfolio returns 0.0095, −0.0369, 0.0311"""
    X = np.asarray(simulate_returns(MU / 12, VOLS / np.sqrt(12), R, 120, 7), dtype=float)
    r_p = X @ np.array([0.5, 0.3, 0.2])
    expect_close(np.round(r_p[:3], 4), [0.0095, -0.0369, 0.0311], what="the first three monthly portfolio returns")


def test_simulate_statistics():
    """20,000 draws reproduce the means, volatilities and correlations closely"""
    X = np.asarray(simulate_returns(MU, VOLS, R, 20_000, 11), dtype=float)
    expect_close(X.mean(axis=0), MU, what="the sample means", abs_tol=0.01, rel=0)
    expect_close(X.std(axis=0, ddof=1), VOLS, what="the sample volatilities", rel=0.03)
    C = np.corrcoef(X, rowvar=False)
    expect_close([C[0, 1], C[0, 2], C[1, 2]], [0.3, 0.1, 0.5], what="the sample correlations AB, AC, BC", abs_tol=0.03, rel=0)


def test_simulate_rejects_impossible_matrix():
    """The impossible correlations raise ValueError, and their repair simulates"""
    _expect_value_error(lambda: simulate_returns(MU, VOLS, BAD, 100, 1), "the correlation matrix has a negative eigenvalue")
    X = np.asarray(simulate_returns(MU, VOLS, repair_correlation(BAD), 100, 1), dtype=float)
    expect_equal(X.shape, (100, 3), what="the shape after repair")


def test_simulate_rejects_bad_input():
    """A negative volatility, or sizes that do not match, raise ValueError"""
    _expect_value_error(lambda: simulate_returns(MU, [0.12, -0.18, 0.25], R, 100, 1), "a volatility is negative")
    _expect_value_error(lambda: simulate_returns(MU, VOLS[:2], R, 100, 1), "two volatilities for three assets")