Cholesky Decomposition

Examples by Prof. Jacob Escobar, EGADE Business School, Tec de Monterrey

Covariance Decomposition

Covariance Matrix \(\left(H\right)\):

\(\sigma_1^2\) \(\sigma_{1,2}\) \(\sigma_{1,3}\)
\(\sigma_{2,1}\) \(\sigma_2^2\) \(\sigma_{2,3}\)
\(\sigma_{3,1}\) \(\sigma_{3,2}\) \(\sigma_3^2\)


\[ \sigma_{i,j} = \rho_{i,j} \sigma_i \sigma_j \]

Traditional decomposition:

\[H = \sigma \rho \sigma\]


Where \(\sigma\) and \(\rho\) are matrices:

\(\sigma =\)

\(\sigma_1\) \(0\) \(0\)
\(0\) \(\sigma_2\) \(0\)
\(0\) \(0\) \(\sigma_3\)

\(\rho =\)

\(1\) \(\rho_{1,2}\) \(\rho_{1,3}\)
\(\rho_{2,1}\) \(1\) \(\rho_{2,3}\)
\(\rho_{3,1}\) \(\rho_{3,2}\) \(1\)

Cholesky decomposition:

\[H = A^TA\]

Where \(A\) is a half-matrix:

\(A =\)

\(a_{1,1}\) \(a_{1,2}\) \(a_{1,3}\)
\(0\) \(a_{2,2}\) \(a_{2,3}\)
\(0\) \(0\) \(a_{3,3}\)


The procedure to calculate \(a_{i,j}\) is sequential, going from top-left to bottom-right:

\[a_{1,1} = \sigma_1^2\]


\[a_{i,i} = \left( \sigma_i^2 - \sum_{k=1}^{i-1} a_{k,i}^2 \right)^{\frac{1}{2}}\]


\[a_{i,j} = \frac{1}{a_{i,i}} \left( \sigma_{i,j} - \sum_{k=1}^{i-1} a_{k,i} a_{k,j} \right)\]

\[j = i+1,i+2,...,n\]

Fortunately, we already have a function in R to calculate the Cholesky factor \(A\): chol( )

By the way, not all matrices have a Cholesky decomposition or factorization. For starters, they must be square matrices, but more formally speaking they must be symmetric positive semidefinite. In practical terms, if \(H\) isn’t “positive semi definite” the procedure will fail because at some point it will try to get a square root of a negative number. Fortunately too, most real-life covariances matrices that you find will be positive semidefinite.

H <- matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1), 3, 3)
H
##      [,1] [,2] [,3]
## [1,]  1.0  0.5  0.5
## [2,]  0.5  1.0  0.5
## [3,]  0.5  0.5  1.0
A <- chol(H)
A
##      [,1]      [,2]      [,3]
## [1,]    1 0.5000000 0.5000000
## [2,]    0 0.8660254 0.2886751
## [3,]    0 0.0000000 0.8164966
# Validation:
t(A) %*% A
##      [,1] [,2] [,3]
## [1,]  1.0  0.5  0.5
## [2,]  0.5  1.0  0.5
## [3,]  0.5  0.5  1.0

Correlated Random Numbers

We can use \(A\) (the Cholesky factor) to generate correlated random numbers. This is very important with any kind of simulation where we involve two or more stochastic elements, in order to produce realistic scenarios.

\[Y = EA\]

Where: * \(E\) is a matrix with \(n\) variables (columns), and each of them has \(m\) normally distributed random observations (rows). Since each series of random numbers is generated independently, \(E\) has the following distribution: ~\(N(0,I)\) - Remember that \(I\) is the identity matrix, so even if we believe we aren’t making any correlation assumption, if we use \(E\) then we are assuming \(0\) correlation among the variables:

  • Identity Matrix (\(I\)):
\(e_1\) \(e_2\) … \(e_n\)
1 0 … 0
0 1 … …
… … … …
0 0 … 1
  • \(Y\) is a matrix with \(n\) variables (columns), and each of them has \(m\) normally distributed random observations (rows). However, we reflected the covariance structure \(H\) in the random numbers \(E\) after multiplying it by the Cholesky factor \(A\). Hence, \(Y\) has the following distribution: ~\(N(0,H)\)

1. Generating Uncorrelated Random Numbers

Since we want to be able to “play” with the correlation matrix for didactic purposes, we will consider the case of only two variables (\(n = 2\)), to ensure the matrix is positive semidefinite.

We will generate 10 thousand normally distributed random observations for each variable (\(m = 10,000\)).

n <- 2
m <- 10000
E <- matrix(rnorm(n*m, mean = 0, sd = 1), nrow = m, ncol = n)
colnames(E) <- c("e1", "e2")
head(E)
##              e1         e2
## [1,] -2.3452683 -0.9509964
## [2,] -1.3614002  1.1780641
## [3,] -0.4398129 -0.8403616
## [4,] -0.6830152 -0.7822130
## [5,] -1.5728934 -1.2167661
## [6,] -1.5144486  1.1267274

Let’s look at the correlation of the variables generated:

cor(E)
##             e1          e2
## e1  1.00000000 -0.01166196
## e2 -0.01166196  1.00000000

In theory, the values we should see are 1’s in the main diagonal entries and 0’s in the off-diagonal entries. Even though we do see some values instead of 0’s, they are very small and statistically non-significant.

2. Correlations, Volatilities & Covariance Matrix

Let’s assume that the correlation below (rho) is the one we want to model. You can play with it by changing it value between 0 and 1.

rho <- 0.5 # You can change this value between 0 and 1
R <- matrix(c(1, rho, rho, 1), 2, 2)
R
##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  1.0

Let’s assume these are the standard deviations we want to model (v1 and v2), you can also play with them by changing their values.

v1 <- 1
v2 <- 1
V <- diag(c(v1,v2))
V
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

The final covariance matrix we have to reflect on \(E\) is:

H <- V %*% R %*% V
H
##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  1.0

3. Finding the Cholesky Factor

A <- chol(H)
A
##      [,1]      [,2]
## [1,]    1 0.5000000
## [2,]    0 0.8660254

Validation: \(H = A^TA\)

t(A) %*% A
##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  1.0

4. Incorporating the Covariance Matrix in the Random Numbers

Y <- E %*% A
colnames(Y) <- c("y1", "y2")
head(Y)
##              y1         y2
## [1,] -2.3452683 -1.9962212
## [2,] -1.3614002  0.3395333
## [3,] -0.4398129 -0.9476810
## [4,] -0.6830152 -1.0189239
## [5,] -1.5728934 -1.8401971
## [6,] -1.5144486  0.2185502

Validation: the simulated correlation and covariance matrices should be similar to the original ones.

cov(Y)
##           y1        y2
## y1 0.9757539 0.4778616
## y2 0.4778616 0.9897971
cor(Y)
##           y1        y2
## y1 1.0000000 0.4862492
## y2 0.4862492 1.0000000

Let’s look at the simulated correlation in a graph:

plot(x = Y[,1], y = Y[,2], col = "blue")

How is the visualization of the following correlations: * \(\rho = 0.50\) * \(\rho = 0.90\) * \(\rho = 0.9999\) * \(\rho = -0.9999\) * \(\rho = 0\)