Cholesky Decomposition & Correlated Random Variables

This time we’re tackling Cholesky Decomposition. To get there, I have to do a little math first. The front half is going to be pretty much pure math, and then I’ll push the financial-engineering side to the back — why we use Cholesky in financial engineering, how we use it, that whole thing.

OK so. A real matrix is a matrix that comes back as itself when you take its complex conjugate.

(These terms… translating them into Korean just gives you stiff Sino-Korean jargon anyway, so I’ll mostly stick with the English names.)

So, a real matrix:

— you could write it like this, right? A matrix with that property is called a real matrix. (Which just means it only has real numbers in it.)

Now there’s also a kind of matrix that comes back as itself even after you take both the complex conjugate and the transpose:

Doing transpose-and-complex-conjugate together is called taking the dagger (which is why the symbol is, well, a dagger). If you dagger it and get the same matrix back, we call that matrix Hermitian (the complex cousin of a symmetric matrix).

But honestly — we’re dealing with a covariance matrix of real stock returns here. So instead of the broad Hermitian definition, I’ll carry on with real symmetric matrices, one notch narrower than Hermitian.

Now, if this A is real symmetric and

— that’s what we call it.

We’re not going to study all of these in detail, though. We’ll just take a quick look at positive definite and move on. Why? Because we’ll use a positive-definite covariance matrix in the examples. Covariance matrices in general are only guaranteed to be positive semidefinite — a zero-variance combination is possible.

What’s a covariance matrix??????? I’ll cover it in my corporate finance series. Haven’t written that one yet, but I’ll get to it in the next few days, and the second it’s up I’ll drop a link here. — Done, linked it :) Minimum Variance Portfolio (Corporate Finance I Studied #10)

OK, let me just touch on a couple of theorems about positive definite before moving on.

Theorem 1. A is real symmetric. Then the necessary and sufficient condition for A to be positive definite is that all eigenvalues of A are positive.

proof. Suppose A is positive definite, λ is an eigenvalue of A, and x is any eigenvector corresponding to λ. Then

Gotta check the arrow going the other way too.

Suppose A’s eigenvalues are all positive. For the set of A’s eigenvectors (orthonormal ones — a real symmetric A always has them) , take any nonzero vector and let

Then any vector can be written as a linear combination of the terms. That is,

Then is

OK — but isn’t requiring a matrix to be Hermitian or real symmetric kind of restrictive? Like, there aren’t even that many of those, so what’s the point…???

Nope, not at all….. In physics, Hermitian shows up everywhere. Especially in quantum mechanics….. ugh, I hated it…… ;_;

Anyway, a few more: if A is Hermitian and positive definite → it’s invertible; if A is Hermitian & positive definite, then det(A) > 0. That kind of stuff… I was going to write them up through property #6, but I think I can get to Cholesky decomposition without them, so I’m just going to skip them. heh heh. If you’re curious, crack open a linear algebra textbook.

On to Cholesky decomposition

We did LU decomposition before, right???? You can think of this as doing LU decomposition here too. Except! The assumption on the matrix A we’re going to hit with LU is different from before: the A we’re going to LU-decompose here is Hermitian (being a covariance matrix) and positive definite.

(From here on this ties into LU Decomposition (Financial Engineering Programming #19), so I’m leaving the link.)

Decompose that kind of A and it’ll look like this:

Since A is symmetric, U has to be L’s mirror image once L’s diagonal is pulled out: $u_{j,i} = l_{i,j}/l_{j,j}$. So I’ll write it like this:

— we can write it like that. In the end,

A can be decomposed like this. This decomposition is called Cholesky Decomposition.

The first reason Cholesky decomposition matters: when you want the solution vector x of the matrix equation $Ax = b$ and A is Hermitian and positive definite, Cholesky factorization uses about half the leading arithmetic work of general LU factorization. The actual running time still depends on the implementation. We want our code as lean as possible, and faster is better, so that’s why Cholesky decomposition matters.

Why we use it here: correlated random variables

But that’s not the meaning of Cholesky decomposition we’re going to put to use here.

Hmm… something I need to lay down first. A covariance matrix of real returns is symmetric, and every portfolio variance is nonnegative. That makes it positive semidefinite, not automatically positive definite!

Zero variance can happen — for example, the difference between two identical return series. For the Cholesky routine here, we additionally assume that no nonzero combination has zero variance. That gives us the positive-definite case.

So we’ll write this positive-definite covariance matrix as Σ, and put it to work. Cholesky decomposition is on the table. That is,

it can be decomposed like this. And the piece we’re actually going to use after the decomposition is L. The lower triangular matrix!!!!

Why are we using it!!!!

To inject correlation into stock price movements.

See, there’s a problem with everything we’ve learned up to this point. Stock price movement is random — yeah, sure — but when you look at two stocks at once, say Hyundai Motor and Kia, we want their random shocks to be correlated. In this model the shocks are Gaussian; the simulated stock prices themselves are lognormal. If you just sketch the two paths any old way like this:

— that probably doesn’t reflect reality very well. The two stock prices are each random on their own, true, but when you look at them together, it’s more accurate to see them as random in a way that’s correlated with each other…..

So now we want to draw two random variables, and we’d like these two random numbers to be random in a somewhat-correlated way. Random variables like that are called correlated random variables.

Back in the day we learned how to draw random numbers $\varepsilon \sim N(0,1)$. But here, if we just draw two of them

and use them as they are, it’s a little inappropriate. So we transform these two random variables appropriately

and want to draw $z_1, z_2$ — correlated random variables, random but with some correlation baked in!!!!

That’s exactly where we use the L matrix in !!!!!!!

We’re going to use it like this to make $z_1, z_2$!!!!!

But then why does doing

give you correlated random variables??!! Honestly — I myself can’t prove this cleanly with math……

BUT!!!!! If we draw $z_1, z_2$ a huge number of times, reconstruct $\mathrm{cov}(z_1, z_2)$, build the covariance matrix from those samples, and show that it comes back to Σ again — wouldn’t that be undeniable???

We can check it with this code!!!!!! We pin down one covariance matrix,

decompose it with a function procedure called Cholesky() (how to build that function procedure comes right after this), and then with MMult we multiply the Gaussian random numbers $\varepsilon_{1,2,3}$ by the L that came out of the Cholesky decomposition.

Then $z_{1,2,3}$ — called rv() here, right??? (I’m guessing it’s short for “random variable”) — we draw as many of them as the iteration count, have the code recompute their means and standard deviations plus the correlation for each pair (1,2), (1,3), (2,3), and the sample results come close to the targets in the covariance matrix, with some sampling noise.

Or they do once a slip in the Cholesky() function further down is fixed: its inner sum should read L(i, k) * L(j, k). As written, it printed the 0.477, 0.953 and 1.053 on the screen; fixed, those come out at about 0.500, 0.900 and 1.000.

What was the point of all that????

What we wanted to check was whether, if we really take the covariance matrix sigma, , Cholesky-decompose it like this, take the L we made there,

and use it like this, the random variables $Z_{1,2}$ that pop out really are correlated with each other???!?~~~

Result of the check: confirmed. heh heh heh

Oh yeah, on to the next part then lololol

The algorithm that spits out L

OK — now let’s study the thing I didn’t cover above, the algorithm that spits out L after the Cholesky decomposition, and wrap it up there. hehe hihi

Say the Cholesky decomposition goes like that. Then let’s look at $a_{11}$ first. Obvious at a glance: — like this?!?!?!

Let’s look at the $a_{21}$ part too. Multiplying row 2 of the front matrix by column 1 of the back matrix gives $a_{21}$, , and everything except that just dies.

Ack?!?!?! Come to think of it, — same deal with this one, — and this guy too, so — we can write it like this. Ahh, then

we can write it like this!!!!

Column 1 of L: done. Let’s look at column 2.

Multiplying row 2 of the front matrix by column 2 of the back matrix all the way through gives $a_{22}$, so

Next up would be $l_{32}$, but before that let me peek at $l_{33}$ and $l_{44}$ real quick.

Ah — now I can see the pattern.

We can write it like this. And looking at row 3 of the front matrix times column 2 of the back matrix,

somehow it feels like we can generalize column 2:

Let’s check column 3 too:

I think I can generalize column 3 too:

And looking at it all written out like this, it looks like I can write it for every entry below the diagonal of L!!!!!

— for $i>j$, we can just write it like this!!!! On the diagonal, use the square-root formula above. And since this is a lower triangular matrix, everything with $(j \gt i)$ is 0~~ — that goes without saying. heh

And this is something I noticed while deriving it just now: looks like we have to work out L in order — column 1, column 2, column 3…. one after another.

Alright, now that we roughly get the principle, let’s study the function procedure that spits out the matrix L after the Cholesky decomposition!!!

It’s short, and simple too, so… well;;;; there’s not really anything else to study here, I guess.

Generating correlated stock prices

So now we can generate the price movements of two or more stocks. For two or more stocks, work out the shock covariances, put a positive-definite covariance matrix into the code above, pull out L, and combine L with the Gaussian draws $\varepsilon$ to make Z. In the price-path code below we use a correlation matrix, then apply volatility separately, and generate the stock prices just like we did in Monte Carlo Method (Financial Engineering Programming #9).

Here’s one example you could try:

Putting in 0.9 as the correlation coefficient,

it came out like this. That was kind of cool, so I tried putting in -0.9.

Then the random shocks should tend to point in opposite directions, right? The price paths will not be exact mirror images — the shocks are not perfectly anticorrelated, and drift still matters.

All I did was put -0.9 in the rho slot, OK???

Whoa. Magic. Wild. LOL. I’m pumped lol

OK, so now instead of generating $z_{1,2}$ to make the stock prices, I just let the independent random numbers $\varepsilon_{1,2}$ drive them. The part of the code I touched is right there

in the for loop: I deleted the rv variables and just put the Gaussian variables in. That way $\varepsilon_{1,2}$ goes in instead of $z_{1,2}$, right?

Ah, this is fun lol. So cool lol

Anyway, I’m done now. Let’s see if I ever study financial engineering again. (-_-)

---

Heads-up. The figures here are the original notes, slips and all. A few things to watch for:

• A matrix is Hermitian when $A^\dagger = A$ (the first picture’s $\bar{A}^T = A^\dagger$ only defines the dagger). And in the definitions, $x^TAx \le 0$ for all $x$ is negative semidefinite, as the label (NSD) says; the word positive there is a slip.
• For a symmetric $A$, the LDU pictures need $l_{ij} = u_{ji}\,l_{jj}$ (for example $l_{21} = u_{12}\,l_{11}$), and in the second line the (2,1) entry of the unit lower factor is $u_{12}$. In the square-root form, the off-diagonal entries carry the root too ($u_{ji}\sqrt{l_{jj}}$), and the last diagonal entry is $\sqrt{l_{nn}}$, not $\sqrt{l_{22}}$.
• The shocks are $\varepsilon_1, \varepsilon_2 \sim N(0,1)$, not $N(1,0)$. In the element formulas, $l_{22} = \sqrt{a_{22} - a_{21}^2/a_{11}}$: the picture writes $a_{n1}^2$ for $a_{21}^2$.
• In the VBA function Cholesky(), the inner sum should read s = s + L(i, k) * L(j, k), as in the formula’s sum of l_ik l_jk; with L(k, j) it picks up the empty upper half, so for the 3 x 3 example the Immediate window shows 1.053, 0.477 and 0.953 instead of about 1.000, 0.500 and 0.900. The 2 x 2 runs and their Excel sheets are unaffected.
• The routine assumes a real symmetric positive-definite input. If a diagonal pivot is nonpositive, report failure rather than using the partly filled L returned after Exit For; a valid singular covariance matrix needs a method that handles semidefinite input.