Getting sets of random correlated variables

Question

For the training of a machine learning model I need to add additional features (macro variables), and these features are correlated. I need to run the model N times, and for each time I have to add the random features.

Now, these features have ranges, for example:

Feature 1 (GDP growth) between 2.0 and 3.5
Feature 2 (Unemployment) between 4.4 and 5.1
Features 3 (Dow Jones index) between 22000 and 24000

With the Cholesky decomposition I can take one of the features and get the others using a correlation matrix. The input would be the fixed value of one of the features, and the result the values of the other correlated features.

But what I need is the system to calculate N sets of random features, within the ranges. So there's no one input to generate the features.

Is there a way to accomplish this (either with or without Cholesky)?

Answer 1

Cholesky (or SVD or any other approach based on matrix multiplication) only works for normal distributions, which your features cannot be, given that they have values within finite intervals.

To see why Cholesky does not work, assume two additional features, which are independent uniform $(U_1,U_2)$. Now you want to create features with correlation $\rho$ by multiplying $U$ with the Cholesky factor $C$ i.e. $Y = U C$. It is not difficult to see that the Cholesky factor is $$ C=\begin{pmatrix} 1 & \rho \\ 0 & \sqrt{1-\rho^2}\end{pmatrix}.$$

This means that $$ Y = (U_1, U_1 + \sqrt{1 - \rho^2}U_2).$$ While the first feature $Y_1=U_1$ is fine, your second variable $Y_2$ is no longer uniform and it will have a different range from $U_2$.

Now you might have the idea that you could start with $Y_2$ and then try to reverse engineer $U_2$ and $U_1$ from this to get the proper input for your simulations. The problem with this idea is that such a pair might not always exist. Starting with arbitrary margins (i.e. standalone independent features) a joint distribution with the prescribed correlation matrix might not exist. There are constraints as to which combinations of marginal distributions and correlation matrices are possible. Have a look at this paper by Embrechts, McNeil and Straumann.

So what CAN you do? If you are lucky and not too picky about precision, the following pragmatic approach might work:

1. Simulate from a Gauss copula see Wikipedia, take as parameters the correlation parameters you would like to obtain
2. Apply the probability integral transform to reorder your features according to 1.
3. Measure the correlations of the reordered features, if they are OK then you are done
4. If they are not OK, repeat 1. with slightly changed parameters of the Gauss copula

Answer 2

Here is one recipe, in case you can live with Spearman rank correlation. (Which you should: linear correlation is often not appropriate in the non-normal case. And in the normal case, there is almost no difference between the two correlation types.)

1. Generate samples of your $k$ features with all the desired attributes. These samples may be random or historical data. The samples need not have the same size.

2. Suppose you want to generate $n$ scenarios. Then generate $n$ $k$-dimensional normals with the desired correlation (i.e. a matrix of size $n$ times $k$). Here you may use Cholesky. (Strictly speaking, you would need to modify the correlations into $2\sin(\rho\pi/6)$. But this correction is so tiny that you may as well leave it out.)

3. Feed these correlated normal variates into the normal distribution function. The result will be $n$ vectors of uniforms which have the same rank correlation as the normals (because the distribution function is monotonically increasing).

4. Feed these uniforms into the inverses of the empirical distribution functions of your samples from step 1. (This is easy: for a uniform variate $u$, sort a given feature sample, and then pick element at position ceiling(u*length(feature_sample)).) Because the inverses are non-decreasing, the rank correlation will remain.

In case you use R: This whole procedure is implemented in the function resampleC in package NMOF (which I maintain). Here is an example.

library("NMOF")

gdp <- runif(10000, min = 2,     max = 3.5)
une <- runif(10000, min = 4.4,   max = 5.1)
dji <- runif(10000, min = 22000, max = 24000)

cor <- array(c(1.0, 0.5, 0.3,
               0.5, 1.0, 0.9,
               0.3, 0.9, 1.0), dim = c(3,3))

smp <- resampleC(gdp, une, dji, cormat = cor, size = 1000)
cor(smp)
##           var1      var2      var3
## var1 1.0000000 0.5140693 0.3041248
## var2 0.5140693 1.0000000 0.9012560
## var3 0.3041248 0.9012560 1.0000000

We can check the ranges.

apply(smp, 2, range)
##          var1     var2     var3
## [1,] 2.000432 4.401090 22005.30
## [2,] 3.499644 5.098867 23998.15

---

A completely-different approach is to create samples of your features (this time of the same length), put them into the columns of a matrix, and then rearrange the elements within the columns so that you come close to the desired correlation matrix. See this answer.

Answer 3

Since the correlation matrix $C=LL^{\top}$ is also $C=U^{\top}U$, where $U$ is the upper triangular matrix, rather than $L$ the lower triangular matrix, you can transform an uncorrelated features matrix $F$ containing features 1, 2 and 3 in its columns by multiplying this $F$ matrix with $U$, giving the correlated features matrix $F_c$:

$$F_c = F U $$

In other words, randomly generate the features in their respective ranges as column vectors without specifying the correlation, concatenate them together into a matrix, then transform them with Cholesky in the formula above into correlated features. I haven't tried, but hopefully the ranges are still retained in the transformed data.