# Getting sets of random correlated variables

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)?

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
• thanks, g g. Instead of implementing the Gauss copula method that you described manually, is there a Python library that already does this? – ps0604 Jun 19 at 11:36
• I am not aware of any Python package which performs this transformation. But the process itself is not too difficult. The main challenge is to understand how the probability integral transform works. And you do not need an explicit formula for the Gauss Copula. You can use multivariate normals. – g g Jun 19 at 11:46
• I read that Gaussian copulas don't do a good job capturing tail dependence, isn't that an issue in my macro variables example? – ps0604 Jun 19 at 13:35
• Gauss copula have zero tail dependence. I suggested Gauss copula because it is the simplest and directly parametrized by correlations. Do your variables have tail dependence and do you know how much? If so, you can replace the Gauss copula by any other copula. The basic structure of the Algo remains the same. – g g Jun 19 at 13:43

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.

• Thanks Enrico, I work with Python, is there a library in Python similar to NMOF? – ps0604 Jun 19 at 11:37
• Rank correlations are a good point! With this exception this is exactly the same idea as I described (somewhat abstract) below. – g g Jun 19 at 12:01
• @ps0604: no idea; sorry. (Some efforts are underway to translate at least parts of the NMOF code into Python, from the R and MATLAB code that was used in the book.) – Enrico Schumann Jun 20 at 7:09
• @g g 19: Indeed; almost. Note my step 4: I map the then-correlated uniforms back via the (empirical) inverse distribution functions of the marginals. In this way the marginal distributions will be preserved: no matter if the samples in step 1 had a (truncated) normal, triangular, uniform or whatever distribution -- the resampled features will follow the same distribution. The size of the new sample needs not be the same as the original feature samples (it is resampling with replacement). – Enrico Schumann Jun 20 at 7:32

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.

• FYI, svd is a better way to do this vs cholesky. – will Jun 16 at 12:48
• @will Can you please provide a link to an example? I see that SVD is used widely, but I cannot find information on generating random variables. – ps0604 Jun 17 at 13:40
• The same way that you use the cholesky decomposition to calculate $\sqrt{\Sigma}$, just use svd instead, and then use the result the same way you would the cholesky decomposition. – will Jun 17 at 21:03
• @will, if the inputs and outputs are used the same way in SVD and Cholesky, what is the benefit of using SVD? – ps0604 Jun 18 at 15:25
• If you have a matrix which is not positive semi definite, you can use SVD to calcualte a "clsoe enough" matrix by flooring the eigen values at zero. It is also more stable in the case where the matrix is nearly positive semi definite. Further, when you have many factors, Cholesky makes the higher dimensions noisier (i.e. any variance reduction techniques you may have used will be somewhat nullified on the higher dimensions. This is less so when you use SVD. – will Jun 18 at 20:49