Numerical experiments, Tips, Tricks and Gotchas

## Generating correlated random variables

### 1. Problem statement

There are $n$ independent (uncorrelated) random variables $x_{1},x_{2},\ldots,x_{n}$ with zero means $$\left\langle \vec{x}\right\rangle =\vec{0}$$ and unit variances $$\left\langle \vec{x}\cdot\vec{x}^{\mathrm{T}}\right\rangle =\hat{I}\label{eq:unitvar}$$ Here the variables are arranged in n-dimensional vectors and the angle brackets denote averaging. The problem is to find a linear combination $\vec{y}$ of the variables $\vec{x}$ with a given covariance matrix $\hat{C}$: $$\hat{C}=\left\langle \vec{y}\cdot\vec{y}^{\mathrm{T}}\right\rangle \label{eq:cov}$$ In other words, it is necessary to find a transformation matrix $\hat{M}$ $$\vec{y}=\hat{M}\cdot\vec{x}\label{eq:Mx}$$ so that the equality (\ref{eq:cov}) takes place.

### 2. Solution

Substituting (\ref{eq:Mx}) into (\ref{eq:cov}) and taking into account (\ref{eq:unitvar}) we get: $$\hat{C}=\left\langle \hat{M}\cdot\vec{x}\cdot\vec{x}^{\mathrm{T}}\cdot\hat{M}^{\mathrm{T}}\right\rangle =\hat{M}\cdot\left\langle \vec{x}\cdot\vec{x}^{\mathrm{T}}\right\rangle \cdot\hat{M}^{\mathrm{T}}=\hat{M}\cdot\hat{M}^{\mathrm{T}}\label{eq:MMT}$$ Any matrix $\hat{M}$ that satisfies (\ref{eq:MMT}) solves the problem.

#### 2.1. Cholesky decomposition

The form of Eq.(\ref{eq:MMT}) suggests that we can use the Cholesky decomposition. Indeed, a covariance matrix is supposed to be symmetric and positive-definite. It always can be decomposed into the product of a lower triangular matrix $\hat{L}$ and its transpose [1]: $$\hat{C}=\hat{L}\cdot\hat{L}^{\mathrm{T}}$$

#### 2.1.1. Example

For the correlation matrix of two variables $$\hat{C}=\left[\begin{array}{cc} 1 & \rho\\ \rho & 1 \end{array}\right]\label{eq:C2x2}$$ The decomposition is $$\hat{L}=\left[\begin{array}{cc} 1 & 0\\ \rho & \delta \end{array}\right]\label{eq:chol}$$ Here $\delta=\sqrt{1-\rho^{2}}$.

#### 2.2. Square root decomposition

A symmetric and positive-definite matrix $\hat{C}$ can be presented as $$\hat{C}=\hat{S}^{2}$$ Here $\hat{S}$ is a square root of $\hat{C}$, which can be also chosen as symmetric $\hat{S}=\hat{S}^{\mathrm{T}}$and positive-definite [2].

#### 2.2.1. Example

For the matrix (\ref{eq:C2x2}) the square root is $$\hat{S}=\left[\begin{array}{cc} c & s\\ s & c \end{array}\right]\label{eq:sqrt}$$ where \begin{eqnarray} c & = & \sqrt{\frac{1+\delta}{2}}\\ s & = & \sqrt{\frac{1-\delta}{2}} \end{eqnarray}

#### 2.3. More solutions

If $\hat{M}$ is a solution of (\ref{eq:MMT}) then its product with an arbitrary orthogonal matrix [3] $\hat{U}$ $$\hat{M^{\prime}}=\hat{M}\cdot\hat{U}$$ is also a solution of (\ref{eq:MMT}). $$\hat{M^{\prime}}\cdot\hat{M^{\prime}}^{\mathrm{T}}=\hat{M}\cdot\hat{U}\cdot\hat{U}^{\mathrm{T}}\cdot\hat{M}^{\mathrm{T}}=\hat{M}\cdot\hat{M}^{\mathrm{T}}$$

### 3. Do we need many solutions?

In the important case of the multivariate normal (Gaussian) distribution we do not need this. The normal distribution is completely defined by its mean and covariance [4]. In general, however, different transformations result in different shape of the distribution functions.

#### 3.1. 2D example

Suppose variables $x_{1}$and $x_{2}$ are uniformly distributed in the intervals $[-\frac{1}{2},\frac{1}{2}]$. The vector $\vec{x}$ is uniformly distributed in the unit square:

Fig. 1. 2D uniform distribution.

#### 3.1.1. Cholesky decomposition

The transformation (\ref{eq:chol}) with $\rho=\frac{1}{2}$ results in the following distribution for $\vec{y}$ :

Fig. 2. Cholesky decomposition.

#### 3.1.2. Square root decomposition

The transformation (\ref{eq:sqrt}) with the same $\rho=\frac{1}{2}$ gives more symmetric distribution for $\vec{y}$ :

Fig. 3. Square root decomposition.

### 4. Implementation

#### 4.1. Cholesky decomposition

Cholesky decomposition [1] is a standard routine in many linear algebra packages.

#### 4.2. Square root decomposition

There are several iterative algorithms [2]. I use diagonalization $$\hat{C}=\hat{U}\cdot\hat{c}\cdot\hat{U}^{\mathrm{T}}$$ $$\hat{C}^{\frac{1}{2}}=\hat{U}\cdot\hat{c}^{\frac{1}{2}}\cdot\hat{U}^{\mathrm{T}}$$ The marix $\hat{c}$ is diagonal $\left(\hat{c}\right)_{i,j}=c_{i,i}\delta_{i,j}$ and $\left(\hat{c}^{\frac{1}{2}}\right)_{i,j}=\delta_{i,j}\,\sqrt{c_{i,i}}$ .

#### 4.3. Steps

It is convenient to split the implementation into the following steps.

0. Random numbers

Needless to say that one should start with good (pseudo-) random numbers $\vec{x}$.

1. Transformation to correlation

First transform $\vec{x}$ to get a correlation matrix $$\hat{R}=\hat{M}\cdot\hat{M}^{\mathrm{T}}$$ where $$R_{i,j}=\begin{cases} 1 & i=j\\ \rho_{i,j} & i\neq j \end{cases}$$ This is the most time-consuming step and can be performed once. Usually a correlation matrix is needed anyway.

2. Transform to covariance $$\hat{C}=\hat{\sigma}\cdot\hat{R}\cdot\hat{\sigma}$$ Here $\hat{\sigma}$ is a diagonal matrix of standard deviations: $\sigma_{i,j}=\sigma_{i}\delta_{i,j}$ The new transformation $$\hat{M^{\prime}}=\hat{\hat{\sigma}\cdot M}$$ gives $$C_{i,j}=\begin{cases} \sigma_{i}^{2} & i=j\\ \rho_{i,j}\sigma_{i}\sigma_{j} & i\neq j \end{cases}$$ Problem solved.

### References

1. Wikipedia. Cholesky decomposition.
2. Wikipedia. Square root of a matrix.
3. Wikipedia. Orthogonal matrix.
4. Wikipedia. Multivariate normal distribution