# Sample many correlated random matrices all with the same pairwise correlation coefficient

I am looking to generate $$K$$ different correlated random matrices, of which the elements all have the same pairwise correlation coefficient. That is, let $$A_1, A_2, \ldots, A_K$$ be $$N \times N$$ random matrices with mean $$\mu$$ and variance $$\sigma^2$$. Then I want that the the pairwise correlation coefficient of elements abides: $$\text{cor}[(a_{ij})_k (a_{ij})_\ell] = \rho \quad \forall i, j$$ and for all $$k \neq \ell$$.

I know that for $$K=2$$, I can simply generate $$A_1 \sim \mathcal{N}(\mu,\sigma^2)$$ and then let $$A_2 = \rho A_1 + \sqrt{1 - \rho^2} B$$, where $$B \sim \mathcal{N}(\mu,\sigma^2)$$. However, I fail to see how to generalize this to $$K$$ matrices. My naive implementation simply uses the above formula for two matrices, but just repeat it, i.e. $$A_{k+1} = \rho A_k + \sqrt{1-\rho^2} B_k$$ where each $$B_k$$ is i.i.d. with mean $$\mu$$ and variance $$\sigma^2$$. But of course, this does not lead to the correct pairwise correlations between all matrices. I also tried Cholesky decomposition with some covariance matrix $$\Sigma$$, but I fail to see what dimensions $$\Sigma$$ should have when I want to generate $$N\times N$$ matrices.

How should I proceed?

• It might be helpful if you could add some code/data to show what you're talking about. Commented Feb 6 at 18:19

Most methods assume your $$\mathbf{A}_1, \mathbf{A}_2, \ldots, \mathbf{A}_K$$ matrices are instead $$n$$-length vectors $$\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_p$$, which form an $$n \times p$$ data matrix $$\mathbf{X}$$. Correlation between all possible $$(j,k)$$ pairs of vectors will result in a $$p \times p$$ measured correlation matrix $$\mathbf{R}$$.

Fundamentally, you really can't "say" or state what correlation coefficients should be via a specified (theoretical) multivariate correlation matrix $$\boldsymbol{\rho}$$ and expect them to appear (i.e., be the same) in a measured (empirical) correlation matrix $$\mathbf{R}$$. This is because:

1. your data vectors are not infinitely large
2. $$\boldsymbol{\rho}$$ may violate Cauchy-Schwarz inequality
3. $$\boldsymbol{\rho}$$ may contain pathologies
4. $$\boldsymbol{\rho}$$ may be positive semi-definite instead of
positive definite.

The measured correlation matrix $$\mathbf{R}$$ does not suffer from any of the above issues, since it's calculated from empirical data vectors $$\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_p$$.

Since the relationship you are using, i.e., $$z=\rho x + \sqrt{1-\rho^2} y$$ blows up with more than 2-3 vectors, if you want to continue to use what you have implemented thus far, then you will need to look at methodological procedures published in ref. 1 below (Rebonato et al.) in order to remove pathologies from a specified $$\boldsymbol{\rho}$$ matrix.

If you want to be able to simulate correlation matrices without pathologies, then use e.g. the Iman & Conover method2. But it still helps to remove pathologies in any simulated correlation matrix. If you want all the correlation matrices to be similar, then just run Iman & Conover's method multiple times using the same specified $$\boldsymbol{\rho}$$ matrix to generate sets of $$\mathbf{x}_j$$ which result in sets of measured correlation matrices $$\mathbf{R}_1, \mathbf{R}_2, \ldots, \mathbf{R}_K$$. Recall, however, none of the resulting $$\mathbf{R}_j$$ matrices will have the same correlation coefficients.

Summing everything up, when simulating a correlation matrix I either set $$\boldsymbol{\rho} \leftarrow \mathbf{R}$$ based on measured correlation and use Iman & Conover, or specify $$\boldsymbol{\rho}$$ and clean it up using methods described in Rebonato et al., then run Iman & Conover.

References:

1. R. Rebonato, P. J¨ackel. The most general methodology to create a valid correlation matrix for risk management and option pricing purposes. J. of Risk. 2:17–27; 2000.

2. R.L. Iman, W.J. Conover. A distribution-free approach to inducing rank correlation among input variables. Commun. Statist. Simulat. Computat. 11(3): 311–334; 1982.

• I do not understand the pathologies. Sure, when I have finite data points the empirical correlation will only be close to the theoretical one, but this should be possible no? What would be wrong about using the Cholesky decomposition to generate $n$ correlated vectors? E.g. with $\Sigma$ the correlation matrix, let $Z = YU$ with $Y$ random uncorrelated samples, and $U$ the Cholesky decomposition for which $U^TU = \Sigma$. See also here. Commented Feb 8 at 8:49
• Nothing's wrong with using the Cholesky decomp, but you do know that since the sample size of your $Y$ is not infinity, none of your correlation matrices will be the same? You're also not supposed to use the same $Y$ to generate another corr matrix - i.e., that's why it's called random. The pathologies will be apparent if you try setting all off-diagonal values of $\Sigma$ to e.g. 0.7 -- which won't work due to Cauchy-Schwarz violations. Commented Feb 8 at 15:53
• I don't understand why I should generate multiple correlation matrices. As I pre-define the correlation matrix, I just need one, right? Commented Feb 9 at 11:35
• Your first sentence states: "I am looking to generate $K$ different correlated random matrices..." If you want multiple $\mathbf{R}$ matrices (based on calculated correlation after simulation), then you only need one specified $\boldsymbol{\rho}$ matrix. But $\boldsymbol{\rho}$ needs to be "cleaned" with Rebonato's method before using Cholesky decomp (I always do this, without question). If you want one $\mathbf{R}$ matrix, you still need one $\boldsymbol{\rho}$ matrix. Commented Feb 9 at 15:40