Generating Correlation Matrices Based on the Boundaries of Their Coefficients

Correlation coefficients among multiple variables are commonly described in the form of matrices. Applications of such correlation matrices can be found in many fields, such as finance, engineering, statistics, and medicine. This article proposes an efficient way to sequentially obtain the theoretical bounds of correlation coefficients together with an algorithm to generate n n correlation matrices using any bounded random variables. Interestingly, the correlation matrices generated by this method using uniform random variables as an example produce more extreme relationships among the variables than other methods, which might be useful for modeling complex biological systems where rare cases are very important.


Introduction
Many important properties of financial models, engineering problems, and biological systems can be represented as correlation matrices, which describe the linear relationships among variables. It is not always the case that these correlation matrices are known; therefore, correlation matrices are an integral part of simulation techniques for solving or analyzing problems in, for example, signal processing [1], portfolio selection [2], factor analytic research [3], genetic modeling [4], and neuroscience [5].
To create a correlation matrix, it is important to ensure that it is valid, meaning that the matrix must be symmetric and positive semi-definite, with the unit diagonal and other elements in the closed interval [21,1]. On the contrary, an invalid correlation matrix is one in which assets or variables cannot be correlated according to the specified relationship. The simplest method for constructing a correlation matrix is to use the rejection sampling method, which generates correlation coefficients using uniform random variables in the closed interval [21,1]. Subsequently, we check whether the matrix is semi-definite and, if not, another correlation matrix is generated. This procedure is repeated until a valid matrix is obtained. Further details of rejection sampling will be described later in this article. For a low-dimensional matrix, it is relatively easy to use rejection sampling, but when the dimension is greater than or equal to four, the chance of finding a valid correlation matrix becomes very low. However, the number of variables in physical or economic systems is normally considerably greater than four, and so the rejection sampling method is considered inefficient for the large-scale construction of correlation matrices.
Instead, for large-dimensional problems, there are several techniques for generating a correlation matrix. These can be classified, based on the relevant objectives or constraints, as follows: 1. Generating of a correlation matrix with predetermined eigenvalues and spectrum [6,7,8]; 2. Generating of a correlation matrix with a given mean value [9]; 3. Generating of a correlation matrix based on a random Gram matrix [10]; and 4. Generating of a correlation matrix in which each correlation coefficient is distributed within its boundaries [11].
This article focuses on the fourth method presenting an efficient algorithm to calculate the theoretical boundaries of correlation coefficients without the use of optimization techniques. Instead, the theoretical boundaries of each correlation coefficient are calculated from the mathematical structure of the correlation matrix constructed by hypersphere decomposition [12]. Although the theoretical work conducted in [11] is similar to the methodology presented here, its primary technique is the optimization approach, whereas our work uses a non-optimization technique. In addition, the sequence for computing the boundaries of each correlation coefficient is heavily reliant on the concept of adjusting the correlation matrix [13] and its boundaries [14]. After finding the theoretical bounds, we present the techniques for generating a correlation matrix.

Valid correlation matrix
It is important to have a common understanding of the definition of a valid correlation matrix. Such a matrix conforms to the following properties: 1. All diagonal entries must be equal to one; 2. Non-diagonal elements consist entirely of real numbers in the closed interval [21, 1]; 3. The matrix is symmetric; and 4. The matrix is positive semi-definite.
The first three requirements are relatively easy to satisfy. However, the final property of being positive semi-definite requires all eigenvalues to be greater than or equal to zero.
Interestingly, a valid correlation matrix (C) can be constructed using a method proposed in [12] in terms of trigonometric functions. The correlation matrix then becomes a function of angles (h(i,j)), which finally gives an efficient way of computing the correlation matrix boundaries without using an optimization method. According to [12], the valid correlation matrix can be described as: Generally, B is a square matrix with n dimensions whose elements are represented by the b i,j in (2). As explained in [15], (2) can be simplified by setting h i,i to zero for all i. B then reduces to a lower triangular matrix, and: It is evident from (4) that matrix B depends solely on h i,j , which is called the correlative angle. The square matrix of correlative angles (h) is defined as: : ð5Þ Thus, a valid correlation matrix can be calculated if the correlative angle matrix (h) in (2.5)is known.
Example 1. Let us assume that the four-dimensional correlative angle matrix is: The matrix B can then be expressed as: b4,1 b4,2 b4,3 b4,4 where which can be written in terms of the correlative angles as

Boundaries of the correlation coefficients
As shown in (6) to (10), a valid correlation matrix can be constructed from the matrix B, and the elements in B are determined by the correlative angles. Consequently, we can determine which elements of B are impacted by changes to the correlative angle in a four-dimensional correlation matrix, from which two important aspects can be inferred: 1. Correlation coefficients in the first column (c i,1 ) depend solely on h i,1 . 2. Other correlation coefficients (c i,j ) for j §2) can be calculated if h p,q are given, where pƒi and qƒjvi.
Because all h i,j are in the closed interval [0, p], the sine functions will produce non-negative values, whereas the cosine functions will output values in the range [21, 1]. Using the correlation coefficients in (10) as an example, it is straight forward to conclude that the boundaries of each correlation coefficient (c i,j ) can be calculated by setting cos h i,j to 21 or 1. Moreover, the boundaries require only h p,q where pƒi and qƒjvi, except for p~i and q~j (although not every h p,q is required), as shown in Table 1. As a result, if c i,j lies within its boundaries and the required h p,q are given, h i,j can be calculated by (11).
The same logic can easily be applied to higher-dimensional correlation matrices, albeit that longer formulas and computational procedures are obtained.

Algorithm for constructing a random correlation matrix
This section describes an algorithm to obtain a correlation matrix by sequentially computing the boundaries of each correlation coefficient, as described in earlier section, and generating uniform random variables (other bounded distributions can always be substituted) within these boundaries. Nevertheless, it is important to note that no optimization is needed to calculate the boundaries of each correlation coefficient. This non-optimization approach is the major difference between our work and that from presented in [11]. Let [0, 1] be the strictly lower triangular matrix of uniform random variables in the closed interval [0, 1], h be the strictly lower triangular matrix of correlative angles, and Y and Z be the strictly lower triangular matrix of lower and upper bounds of the correlation coefficients, respectively. The four-step algorithm for constructing an n|n correlation matrix is then: Step 1: Calculate correlation coefficients in the first column N For i = 1, . . ., n, set c i,1~{ 1z2|u i,1 , b i,1~ci,1 , and extract h i,1 .
Step 2: Calculate the remaining correlation coefficients from the third row to the last row and from the second column to the last column of each row. N The method for calculating these boundaries is explained in the earlier section. Please see Table 1 for an example of the upper and lower bounds using a four-dimensional correlation matrix.
N During our large numerical experiment, numerical instability occurs when the boundary gap (z i,j {y i,j ) becomes very small. As a result a threshold factor (K) is introduced. This reduces instability by forcing every correlation coefficient with a boundary gap of less than K to be centered within its boundaries. Larger value of K will produce a more stable system, but imply less randomization in the c i,j .
N Extract h i,j using similar formulas to those shown in (11).

End End
N Create a symmetric correlation matrix with unit diagonal elements based on all generated correlation coefficients.
Step 3: Randomly reorder the correlation matrix. The underlying concept of this step is to ensure that every correlation coefficient is equally distributed. Without this step, the cumulative distribution function (CDF) of correlation coefficients will not be the same (see Figure 1). After applying random reordering, the CDF of the same correlation coefficients will be almost identical, as displayed in Figure 2.
Step 4: Check the validity of the correlation matrix. Even though the above steps should theoretically generate a valid correlation matrix, in some cases numerical instability can still occur. We can detect two major causes of instability: Firstly, K is too low relative to the dimension of matrix; Secondly, generated correlation coefficients are very close to the boundaries. Based on our experiments, in which 1 million 100|100 correlation matrices were generated with K~0:01, there is only 0.0167% (or 167 matrices) probability that an invalid correlation matrix will occur. Although the probability of an invalid matrix is very small, it is non-zero. That is why this step is necessary, to ensure that invalid correlation matrices will be rejected. The two basic procedures of this step are: 1. Check the minimum eigenvalue. If it is negative, the correlation matrix is invalid. Otherwise, the correlation matrix is valid. 2. Reject the invalid correlation matrix, and regenerate the correlation matrix by returning to step 1 In addition, from (1) to (4), we can generate a valid correlation matrix directly from random sample of correlative angles. Unfortunately, based on our experiment, this direct method is not numerically stable. As a result, one may not be able to use the matrix generated from this method in some applications. Thus, we believe that our new algorithm is superior in terms of numerical stability.
Example 2. For a five-dimensional correlation matrix, let us assume that the uniform random matrix U described in step 1 of the algorithm is: : ð12Þ The lower-bound matrix Y , upper-bound matrix Z, and correlation matrix C (before being randomly reordered) can then be generated as follows:    : As the minimum eigenvalue of C in (15) is 0.00510, the correlation matrix is positive semi-definite. This confirms that C is a valid correlation matrix.

Results
All numerical tests in this study were conducted with MATLAB 7.8.0 (R2009a) on an Intel(R) Core TM 2 Duo CPU T6600 at 220 GHz with 3.50 GB of RAM. The computational performance and probability distribution function (PDF) of the proposed algorithm (NA) with K = 0.01 was evaluated and compared with the following two algorithms:

Rejection sampling method (RS)
The rejection sampling method uses uniform random variables in the closed interval [21, 1] to represent each correlation coefficient in the symmetric correlation matrix. The correlation matrix will be rejected if it is invalid.

Randcorr function of MATLAB (RC)
This algorithm is implemented as a MATLAB function, and is based on the work in [6] and [7].

Computational performance
The computational performance of each algorithm is primarily measured by the expected run time (T exp ), which can be calculated from the average run time (T avg ) divided by the probability of the generated correlation matrix being valid (P valid ). T avg includes the time taken to construct the correlation matrix and calculate the minimum eigenvalue. The performance summary of the three algorithms over 1 million simulations is illustrated in Table 2.
With a P valid score of 100% in all cases, both NA and RC algorithms are evidently stable. Moreover, the RC algorithm has the fastest expected run time when the dimension exceeds 35, although the RS algorithm is the fastest for dimensions of two and three. However, the RS method then becomes slower than the NA algorithm when n §4, and slower than RC for n §5. Even worse, the RS method cannot generate a valid correlation matrix for dimensions larger than seven, mainly due to the significant drop in Pvalid. Hence, the RS method is not very useful in practice. For dimensions from 4-35, the NA algorithm outperforms RS and RC in terms of expected run time.

Probability distribution function
To compare the PDF of the coefficients of correlation matrices, c 2,1 and c 5,4 are drawn from 100,000 valid 5|5 correlation matrices constructed by the above algorithms. Comparing  Figures 3 and 4, we can clearly see that the correlation coefficients generated by the RC algorithm have significant differences.This fact is verified by the kurtosis and standard deviation of the RC algorithm, which are given in Table 3. In general, correlation coefficients from the NA and RC algorithms are equally distributed, but the NA algorithm produces a higher standard deviation and lower kurtosis, which implies more extreme correlation coefficients than the other algorithms.

Discussion
In this paper, we have presented an efficient method to calculate the boundaries of correlation coefficients. We also demonstrated a   Table 3. Statistical summary of random correlation coefficients (c 2,1 and c 5,4 ). technique for generating correlation matrices using any bounded random variable distribution within the boundaries of each correlation coefficient. However, this method causes the correlation coefficients to be unevenly distributed. Thus, we incorporated a technique for random reordering to ensure the even distribution of all correlation coefficients. The performance of the proposed algorithm was compared to that of other algorithms. It was shown that the new algorithm could efficiently construct correlation matrices, particularly when the dimension of the matrix was in the range 4-35. In theory, our algorithm should always return valid correlation matrices. However, without setting a threshold factor and using rejection sampling logic, the algorithm exhibited some numerical instability when the dimension became large. It is possible to adjust invalid matrices to form valid ones; this method has been developed in many studies [16,17,18]. Therefore, we strongly believe that our new algorithm is useful in the many applications where extreme cases are very important. More importantly, the uniform distribution can be replaced with any bounded distribution.