Chessboard Distributions and Random Vectors ... - Semantic Scholar

Chessboard Distributions and Random Vectors with Speci ed Marginals and Covariance Matrix Soumyadip Ghosh and Shane G. Henderson Department of Industrial and Operations Engineering University of Michigan Ann Arbor, MI 48109-2117, U.S.A. August 29, 2000 Abstract There is a growing need for the ability to specify and generate correlated random variables as primitive inputs to stochastic models. Motivated by this need, several authors have explored the generation of random vectors with speci ed marginals, together with a speci ed covariance matrix, through the use of a transformation of a multivariate normal random vector. A covariance matrix is said to be feasible for a given set of marginal distributions if a random vector exists with these characteristics. We develop a computational approach for establishing whether a given covariance matrix is feasible for a given set of marginals. The approach is used to rigorously establish that there are sets of marginals with feasible covariance matrix that the normal transformation technique referred to above cannot match. An important feature of our analysis is that we show that for almost any covariance matrix (in a certain precise sense), our computational procedure either explicitly provides a construction of a random vector with the required properties, or establishes that no such random vector exists. We also provide two new methodologies that may be used to deal with the situation where one cannot exactly match the desired marginals and covariance matrix using the normal transformation technique. The new methodologies possess certain advantages over other approaches that have been suggested in the past.

1

Introduction There is a growing need for the ability to specify and generate random vectors consisting of correlated observations as primitive inputs to stochastic models. For example, in a manufacturing setting, the processing times of a single job at dierent stations may be correlated due to characteristics of the job such as size. In determining reservoir release rules, the in ows of water to dierent reservoirs are invariably correlated. In generating random test problems for a given algorithm, it is advantageous to ensure that some elements of the problem are correlated (Hill and Reilly 1994, 2000). Further applications have recently been reported in cost analysis (Lurie and Goldberg 1998), and in decision and risk analysis (Clemen and Reilly 1999). Perhaps the \ideal" approach is to specify the full joint distribution of the random vector. This approach is typically limited to situations where the marginal distributions are all from the same parametric family. For methods of this type see, for example, Devroye (1986) and Johnson (1987). But the case where the marginals are not all from the same parametric family aords far greater modeling generality, and is perhaps the case of more interest from a practical standpoint. The primary diculty in this case is that a tremendous amount of information is typically required to specify (and t) such a joint distribution. Furthermore, special methods must be devised to generate random vectors with the given joint distribution, and this can be a practically insurmountable problem for a model of even moderate complexity (Law and Kelton 2000, p. 479). A practical alternative is to only specify the marginal distributions of the random variables, together with the correlation matrix or covariance matrix. (Note that this information does not necessarily uniquely specify the distribution.) The covariance measure could be Spearman's rank covariance, Pearson's product-moment covariance, Kendall's  , or any other convenient covariance measure. In this paper, we will focus on Pearson's product-moment covariance, and Spearman's rank covariance, because of their wide use and acceptance in application settings. It is important to note that we are restricting attention to the generation of nite-dimensional random vectors. As such, we are not attempting to generate a time series with given correlation properties. For such studies, see for example Cario and Nelson (1996), Melamed et al. (1992), and Lewis et al. (1989). Hill and Reilly (1994) describe a method for generating random vectors with speci ed marginals and covariances through mixtures of extreme correlation distributions. The approach is very eective for random vectors of low dimension (d  3 say), but the computational requirements quickly become excessive for higher dimensional random vectors. There is another diculty with this approach. We

2

say that a covariance matrix is feasible for a given set of marginal distributions if a random vector exists with the prescribed marginals and covariance matrix. We show (Section 1) that there are sets of marginals with feasible covariance matrix that cannot be matched using the technique of Hill and Reilly (see Section 1 below). Cario and Nelson (1997) described the \NORmal To Anything" (NORTA) method for generating random vectors with prescribed covariance matrix. The NORTA method basically involves a componentwise transformation of a multivariate normal random vector, and capitalizes on the fact that multivariate normal random vectors are easily generated; see Law and Kelton 2000, p. 480. Cario and Nelson traced the roots of the method back to Mardia (1970) who looked at bivariate distributions, and to Li and Hammond (1975) who concentrated on the case where all of the marginals have densities (with respect to Lebesgue measure). Iman and Conover (1982) implemented the same transformation procedure to induce a given rank correlation in the output. Their method is only approximate, in that the output will have only very approximately the desired rank correlation. Clemen and Reilly (1999) described how to use the NORTA procedure to induce a desired rank correlation in the context of decision and risk analysis. Lurie and Goldberg (1998) implemented a variant of the NORTA method for generating samples of a predetermined size. A natural question to ask is whether the NORTA procedure can match any feasible covariance matrix for a given set of marginals. Both Li and Hammond (1975) and Lurie and Goldberg (1998) give examples where this does not appear to be the case. However, the random vectors that are proposed in these papers as counterexamples are not proved to exist, and so the question has not yet been completely settled. In this paper, we prove that there are feasible covariance matrices for a given set of marginals that the NORTA method cannot match. To establish this result we derive a computational procedure for establishing whether or not a given covariance matrix is feasible for a given set of marginals, and if so, explicitly providing a joint distribution with the required properties. We call the constructed distributions \chessboard" distributions because of their structure; see Section 2. It should be noted though, that in the case of 2-dimensional random vectors, the class of feasible covariance matrices has been completely characterized; see Whitt (1976). However, for dimensions 3 and greater, little is known.

Remark 1 In the case where all of the marginal distributions have densities with respect to Lebesgue measure, the chessboard distribution we construct has a joint density with respect to d-dimensional Lebesgue measure. In this case, we can, and do, refer to a chessboard density.

3

If NORTA cannot precisely match a feasible covariance matrix, it is still possible to use NORTA to obtain the desired marginals exactly, and the desired covariance matrix approximately. Lurie and Goldberg (1998) gave an alternative approach to this problem, but we believe that our solution has properties that make it more desirable; see Section 5. We also provide a new method for generating random vectors with speci ed marginals and covariance matrix based on chessboard distributions. The new method requires a nontrivial amount of computational eort during initialization, but will typically be quite competitive with the NORTA method while actually generating the required random vectors. It has the potential to generate virtually any (in a certain precise sense; see Section 6) feasible rank covariance for a given set of marginals. The ability to more closely match certain covariance matrices than the NORTA method should be seen as a clear advantage of this approach. However, it should be noted that there are covariance matrices that can be matched by NORTA, but cannot be matched exactly by a chessboard distribution, even though chessboard distributions can come arbitrarily close to the desired covariance matrix. The philosophy of specifying marginals and correlations to model dependent random variates is clearly an approximate one, since the joint distribution is not completely speci ed. Therefore, one should be willing to live with reasonable discrepancies in the covariance matrix from that desired. What is \reasonable" depends on the application at hand, and this is the reason that we describe the new method for random vector generation alluded to above. If one is not willing to live with reasonable (again, this is a relative term) discrepancies from the desired covariance matrix, then perhaps a more careful approach to specifying the dependence structure is warranted. So while we present an alternative method that can be used in cases where NORTA does not exactly match the required covariances, we still believe that NORTA should typically remain the method of choice for this kind of variate generation. We view the primary contributions of this paper as follows. 1. We provide a computational procedure for determining whether a given covariance matrix is feasible or not for a copula, i.e., for a random vector with uniform(0; 1] marginals. (In this case the rank covariance and product-moment covariance are identical.) If the covariance matrix is feasible, then an explicit construction of a joint density with these properties, that we call a chessboard density, is provided. The method works for almost all covariance matrices; see Section 2. To the best of our knowledge, this is the rst example of such a characterization of feasible covariance matrices. This case is important because it is central to the analysis and use of rank covariance for continuous marginals; see Section 6 for example.

4

2. We provide a computational procedure for determining whether a given Pearson product-moment covariance matrix is feasible or not for a given set of more general marginal distributions. If the covariance matrix is feasible for the given set of marginals, we provide an explicit construction of a chessboard distribution with the desired properties. Again, we believe that this is the rst example of such a characterization. 3. We rigorously establish that there are feasible covariance matrices that cannot be matched using the NORTA method. 4. We provide a simple modi cation to the initialization phase of the NORTA method that enables one to use the NORTA method to closely approximate the desired covariance matrix. The modi cation involves the solution of a semide nite program, and works in both the rank correlation and Pearson product-moment correlation cases without any specialization. Based on a small computational study, it appears that when one cannot exactly match a desired covariance, the discrepancy between the desired and realized covariance matrices is quite small. 5. We develop a new method for generating random vectors with speci ed covariance structure using chessboard distributions that can exactly match covariances that the NORTA method cannot. The new method could involve a nontrivial amount of initialization eort, but certainly be competitive (computationally speaking) with the NORTA method while actually generating random vectors. The remainder of this paper is organized as follows. In Section 1 we review the NORTA method, and describe how it may be used to match a given Pearson product-moment covariance matrix, or a given rank covariance matrix. In Section 2 we develop the theory of chessboard distributions in the case where all of the marginal distributions are uniform(0; 1]. The chessboard distribution concept is then extended to more general marginals in Section 3. Next, in Section 4 we present a small computational study that sheds light on when we might expect the NORTA method to be unable to match a feasible covariance matrix. We also provide several examples where this occurs. Then, in Section 5, we present our modi cation of the NORTA method that involves semide nite programming. Finally, in Section 6 we present an alternative method for generating a random vector with speci ed marginals and covariance matrix, based on chessboard distributions. We are developing a software package that implements the overall NORTA process, incorporating the semide nite programming modi cation suggested in Section 5. The software is similar in spirit to the ARTAFACTS/ARTAGEN package developed by Cario and Nelson (1997b). Indeed, we are adapting their software for the task. When the software is ready, it will be freely available from the second

5

author's website Henderson (2000).

1 The NORTA method Cario and Nelson (1997) described the \NORmal To Anything" (NORTA) method for generating i.i.d. replicates of random vectors with speci ed marginals and covariance structure. In this method, one starts by generating a random vector Z with a multivariate normal distribution and transforms Z to obtain a random vector X = (X1 ; : : : ; Xd) with the desired marginals and covariance structure. Let Fi be the distribution function of Xi , for i = 1; : : : ; d. The NORTA method generates i.i.d. replicates of X by the following procedure. 1. Generate an IRd valued standard normal random vector Z = (Z1 ; : : : ; Zd ) with mean vector 0 and covariance matrix Z = (Z (i; j ) : 1  i; j  d), where Z (i; i) = 1 for i = 1; : : : ; d. 2. Compute the vector X = (X1 ; : : : ; Xd ) via

Xi = Fi?1 ((Zi ));

(1)

for i = 1; : : : ; d, where  is the distribution function of a standard normal random variable, and

Fi?1 (u) = inf fx : Fi (x)  ug:

(2)

A vector X generated by this procedure will have the prescribed marginal distributions. To see this, note that each Zi has a standard normal distribution, so that (Zi ) is uniformly distributed on (0; 1), and so Fi?1 ((Zi )) will have the required marginal distribution. The covariance matrix Z should be chosen so that it induces the required correlation structure on X . There are many measures of correlation between two random variables, but perhaps the two most popular are Pearson's product-moment correlation, and Spearman's rank correlation.

1.1 Pearson's Product-Moment Correlation In this section we specialize to the case where one wishes to generate X with prespeci ed Pearson product-moment covariance matrix X , where X (i; j ) = EXi Xj ? EXi EXj for 1  i; j  d. This is the case that Cario and Nelson (1997) examined. Note that this is equivalent to prespecifying the correlation matrix, since the marginal distributions are also prespeci ed. To ensure

6

that the required correlations are de ned, we make the assumption that EXi2 < 1 for i = 1; : : : ; d. It turns out that choosing Z to ensure the correct covariance matrix X is a nontrivial problem. Each term cov(Xi ; Xj ) is a function of only cov(Zi ; Zj ) (Cario and Nelson 1997), and so the problem reduces to d(d ? 1)=2 separate problems of selecting cov(Zi ; Zj ) to match cov(Xi ; Xj ) = X (i; j ). Unfortunately, there is no general analytical expression relating these two quantities, and hence we cannot determine the exact Z that is to be used analytically. De ne the function cij (cov(Zi ; Zj )) = cov(Xi ; Xj ), where Xi and Xj are de ned via (1). Cario and Nelson (1997) established that under very mild conditions, the function cij is a continuous nondecreasing function of Z (i; j ). This result allows us to perform an ecient numerical search for values Z (i; j ) that yield cij (Z (i; j )) = X (i; j ) for i < j: (3) We take Z (i; i) = 1 for i = 1; : : : ; d, and for i > j , set Z (i; j ) = Z (j; i) to ensure that Z is symmetric. Alternatives to the numerical search suggested by Cario and Nelson (1997) include the use a stochastic root- nding algorithm (Chen 2000), or polynomial expansions (van der Geest 1998). Unless otherwise stated, we henceforth assume that a solution to (3) exists. One might hope that if the matrix Z satis es (3), then Z could be used in the NORTA method to generate i.i.d. replicates of X . Unfortunately, the results of this paper prove that this is not always the case. In fact, there exists a feasible covariance matrix for a 3-dimensional random vector with uniform marginals on (0; 1] that cannot be generated with the NORTA procedure. The problem arises when the matrix Z as determined from (3) is not positive semide nite, in which case it is not a valid covariance matrix. Li and Hammond (1975) suggested the following example to illustrate this important fact. Let X1 ; X2 and X3 be 3 uniformly distributed random variables on (0; 1] with correlation matrix

0 1 1 ? 0 : 4 0 : 2 B C RX = B B@ ?0:4 1 0:8 CCA : 0:2

0:8

1

In the special case when X has uniform marginals, the equations (3) can be solved analytically. In particular, Kruskal (1958) showed that the (unique) solution to (3) is given by Z (i; j ) = 2 sin[2X (i; j )]:

(4)

For the Li and Hammond example, the (unique) matrix Z that can be determined using (4) (with

7

X = RX =12) is not positive semide nite. It is important to observe though, that this is a counterexample only if the postulated random vector exists, and Li and Hammond did not show this.

Remark 2 It is straightforward to show that the Li and Hammond example cannot be generated using the extremal distributions method of Hill and Reilly (1994). One simply attempts to solve the linear program suggested by Hill and Reilly (1994), which turns out to be infeasible. Therefore, if the Li and Hammond example exists, it shows that there are feasible covariance matrices that cannot be matched using the extremal distributions technique.

Lurie and Goldberg (1998) gave an example with nonuniform marginals and positive de nite covariance matrix for which the solution to (3) is also not positive semide nite. They did not establish that the postulated random vector exists. When all of the marginals have continuous distribution functions, a natural alternative to the numerical search procedure mentioned earlier is to \work in the Gaussian space". In other words given a set of data with known (or tted) marginals with continuous distribution functions, we rst transform the data set into normal random variates using the inverse of the transformation (1). We can then compute an empirical covariance matrix Z and use this covariance matrix in the NORTA procedure. (If the distribution function F of a random variable X is not continuous, then F (X ) does not have a uniform distribution on (0; 1), and so one will not obtain a normally distributed random variable using ?1 (F (X )). Therefore, the continuity of the marginal distribution functions is needed.) This approach is certainly simpler than a numerical search procedure, but it has two important drawbacks. First, it requires a set of input data, which may not be available in general. But second, and perhaps more importantly, this procedure does not necessarily ensure that the resulting X variates will have the required covariance structure. To see why, observe that the transformed normal random variables mentioned above are unlikely to have a joint normal distribution. Therefore, the correlations of the jointly normal random variables used in the NORTA method using Z will be unlikely to transform through the NORTA procedure to yield the desired covariance matrix for X , as one might otherwise expect. This is a subtle point, but one that is worth bearing in mind.

1.2 Spearman's Rank Correlation In this section we specialize to the case where one wishes to generate X with prespeci ed Spearman's rank covariance matrix X , where X (i; j ) = rcov(Xi ; Xj ) = EFi (Xi )Fj (Xj ) ? EFi (Xi )EFj (Xj )

8

for 1  i; j  d. This is the case treated by Clemen and Reilly (1999). In contrast to product-moment covariance, the rank covariance is always de ned. Note that when the Xi 's all have nondegenerate distributions, specifying the rank covariance matrix is equivalent to prespecifying the rank correlation matrix with (i; j )th element rcor(Xi ; Xj ) = cor(Fi (Xi ); Fj (Xj )) =

rcov(Xi ; Xj ) : [varFi (Xi )varFj (Xj )]1=2

It turns out that choosing Z to ensure the correct rank covariance matrix X is easy in the special case that all of the distribution functions Fi are continuous (i = 1; : : : ; d). To see why this is the case, observe that if Xi and Xj have continuous distribution functions Fi and Fj , then Fi (Fi?1 (u)) = u for all u 2 (0; 1), and similarly for Fj . The rank covariance rcov(Xi ; Xj ) between Xi and Xj as generated by a NORTA procedure is therefore given by rcov(Xi ; Xj ) = cov(Fi (Xi ); Fj (Xj )) = cov(Fi (Fi?1 ((Zi ))); Fj (Fj?1 ((Zj ))))

(5)

= cov((Zi ); (Zj )):

(6)

But (6) is precisely the quantity X (i; j ) in (4). Therefore, given a desired rank covariance matrix X , we simply compute Z = Z via (4) and use this within the NORTA procedure. Observe that if the random vector in the Li and Hammond example (given above) exists, then it is again an example showing that there are feasible rank covariance matrices for a given set of marginals that cannot be matched using a NORTA procedure. In the case where Fi (say) is not continuous, equation (6) no longer holds, because Fi (Fi?1 (x)) 6= x for some x 2 (0; 1). Therefore, the analytical expression (4) cannot be used. However, one could use a numerical search procedure as in Cario and Nelson (1997) to identify the covariance Z (i; j ) that yields the required rank covariance rcov(Xi ; Xj ). This follows since the rank covariance (5) between Xi and Xj is a nondecreasing continuous function of the covariance between Zi and Zj . The nondecreasing property follows immediately from the proof of Theorem 1 in Cario and Nelson (1997), and the fact that the function Fi (Fi?1 ((
))) is nondecreasing. Continuity follows from Theorem 2 of Cario and Nelson.

1.3 To summarize We have just seen that the NORTA procedure can be used to generate a random vector with the desired product-moment, or rank, correlation. The same is true when correlation is measured by Kendall's  ,

9

as Clemen and Reilly (1999) note. The dierence lies in how the covariance matrix Z of the initial normal random vector is determined. We have already noted that if the Li and Hammond example exists, then we may conclude that there are sets of marginals with feasible (product-moment or rank) correlation matrices that cannot be matched using the NORTA procedure. In the following section, we introduce a linear programming based method that may be used to investigate the existence or not of random vectors with uniform marginals and speci ed covariance matrix.

2 Copulas A copula can be characterized as the joint distribution function of a random vector with uniform marginals on (0; 1] (Nelsen 1999). In this section, we develop a method for constructing the joint density of a 3-dimensional copula with prescribed covariance matrix. (Note that Pearson's productmoment covariance and Spearman's rank covariance coincide in the copula case.) The approach is easily carried over to the general case of a d-dimensional copula with prescribed covariance matrix. We will let X = (X1 ; X2 ; X3 ) denote a random variable with such a density and let  = (ij : 1  i; j  3) be the desired covariance matrix. We rst construct the probability mass function (pmf) of a random vector Y = (Y1 ; Y2 ; Y3 ) whose marginals are discretized versions of the marginals of X . The pmf will be constructed to try to ensure that Y has covariance matrix  (except for the diagonal entries). From this pmf, we then construct the density of X in such a way that the o-diagonal entries in the covariance matrix are maintained. (The diagonal elements are determined by the marginals of X .) This then yields the required construction. Our notation will appear partly redundant at times, but this is done to ensure consistency with Section 3 where we will extend these ideas to more general marginal distributions. Let n  1 be an integral parameter that determines the level of discretization that will be performed. Let yi;k = nk ; k = 0; : : : ; n be the set of points that divide the range (0,1] of the ith variable into n equal length sub-intervals. For k = 1; : : : ; n and i = 1; 2 and 3, let

Yi;k = E [Xi j Xi 2 (yi;k?1 ; yi;k ] ] = 2k2?n 1 ; be the conditional mean of Xi given that it lies in the kth sub-interval. The support of the random vector Y is the mesh of points

f(Y ;i ; Y ;j ; Y ;k ) : 1  i; j; k  ng: 1

2

3

10

(7)

Let

q(i; j; k) = P (Y1 = Y1;i ; Y2 = Y2;j ; Y3 = Y3;k )

be the probability that Y equals the (i; j; k)th point in the support of Y , so that q represents the pmf of the random vector Y . (Note that it is not the pmf itself, since the function q is de ned on integers, while the domain of the pmf is contained in the unit cube.) Consistent with the notion that Y is a discretized version of X , we also have that

q(i; j; k) = P (X 2 C (i; j; k)); where the cell C (i; j; k) represents the cube of points surrounding the (i; j; k)th point in the support of Y . More precisely,

C (i; j; k) = f(x1 ; x2 ; x3 ) : y1;i?1 < x1  y1;i ; y2;j?1 < x2  y2;j ; y3;k?1 < x3  y3;k g: We then see that n X

j;k=1 n X i;k=1 n X i;j =1

q(i; j; k) = P (Y1 = Y1;i ) = P (X1 2 (y1;i?1 ; y1;i ]) = n1 ; 8i = 1; : : : ; n;

(8)

q(i; j; k) = P (Y2 = Y2;j ) = P (X2 2 (y2;j?1 ; y2;j ]) = n1 ; 8j = 1; : : : ; n;

(9)

q(i; j; k) = P (Y3 = Y3;k ) = P (X3 2 (y3;k?1 ; y3;k ]) = n1 ; 8k = 1; : : : ; n;

(10)

q(i; j; k)  0 8i; j; k = 1; : : : ; n:

(11)

With these constraints satis ed, we then have that EYi = 1=2 = EXi for i = 1; : : : ; 3. To see this, note that for Y1 , we have that

EY1 = = =

n X

i;j;k=1 n X i;j;k=1 n X i=1

Y1;i q(i; j; k) E [X1 j X1 2 (y1;i?1 ; y1;i ] ]P (X 2 C (i; j; k))

E [X1 j X1 2 (y1;i?1 ; y1;i ] ]P (X1 2 (y1;i?1 ; y1;i ])

= EX1 : Recall that our intermediate goal is to match the covariance matrix of Y to that of X (with the exception of the diagonal elements), and we propose to do this using linear programming methods. If Cij = cov(Yi ; Yj ), then we want to minimize

jC ?  j + jC ?  j + jC ?  j: 12

12

13

13

11

23

23

(12)

Now

C12 =

n X i;j;k=1

Y1;i Y2;j q(i; j; k) ? EY1 EY2 ;

which is a linear function of the q(i; j; k)'s, with similar linear expressions for C13 and C23 . Furthermore, the matrix  is simply a parameter, and so, for example, C12 ? 12 is a linear function of the q(i; j; k)'s. Using a standard trick in linear programming, we can also represent jC12 ? 12 j in a linear fashion, and similarly for the other terms in (12) as follows. De ne Zij+ to be the positive part of the dierence Cij ? ij , i.e.,

Zij+ = (Cij ? ij )+ = maxfCij ? ij ; 0g; and Zij? to be

(Cij ? ij )? = ? minfCij ? ij ; 0g:

We can now attempt to match the covariances of Y to those of X using the LP min

P P 2 =1

i

+ ? j i (Zij + Zij ) 3 = +1

subject to Cij ? ij = Zij+ ? Zij?; i = 1 to 2 and j = i + 1 to 3

Zij+  0; Zij?  0; together with constraints (8), (9), (10) and (11).

This LP is always feasible since a product copula where the Yi 's are independent can be easily constructed by setting all q(i; j; k) = n?3 . Also, the objective function of the LP is bounded below by 0. Hence an optimal solution exists. We also know from standard results in linear programming that in any optimal solution to the LP, Zij+
Zij? = 0, i.e., at most one of Zij+ and Zij? is > 0. If the optimal objective value for the above LP turns out to be 0, then we have constructed a discretized joint probability mass function for Y that has the desired covariance structure. So, for i 6= j , we have cov(Yi ; Yj ) = ij :

(13)

Now, the discretized random vector Y does not possess continuous uniform (0,1] marginals. However, we can construct a random vector X with continuous uniform marginals from Y in such a way that cov(Yi ; Yj ) = cov(Xi ; Xj ) for i 6= j , i.e, the covariances are preserved. Assuming the optimal objective value of the LP is 0, this then yields an explicit construction of a random vector with the desired marginals and convariance matrix.

12

By conditioning on the cell containing X , we see that the requirement that cov(Y1 ; Y2 ) = cov(X1 ; X2) is equivalent to n X

i;j;k=1

q(i; j; k)
Y1;i Y2;j ? EY1 EY2 =

n X

i;j;k=1

E [X1 X2 jX 2 C (i; j; k)]
P (X 2 C (i; j; k)) ? EX1 EX2 (14)

But, EY1 = EX1 and EY2 = EX2 , and so (14) can be re-expressed as n X

i;j;k=1

q(i; j; k)
(E [X1 jX 2 C (i; j; k)]
E [X2 jX 2 C (i; j; k)] ? E [X1 X2 jX 2 C (i; j; k)]) = 0:

(15)

Equation (15) could be satis ed in many ways, but perhaps the simplest is to note that (15) will hold if, conditional on X lying in C (i; j; k), X1 ; X2 and X3 are independent. In that case, each term in the sum (15) is 0. One can ensure that this conditional independence holds, while simultaneously ensuring that X has the correct marginal distributions, by setting the density of X within the cell C (i; j; k) to that of independent, uniformly distributed random variables, scaled so that the total mass in the cell is q(i; j; k). To be precise, if f is the density of X , then for any x 2 C (i; j; k), we set

f (x) = n3 q(i; j; k):

(16)

In a sense, we are \smearing" the mass q(i; j; k) uniformly over the cell C (i; j; k). Theorem 1 below proves that if the optimal objective value of the LP is 0, then the density f so constructed has the desired marginals and covariance matrix.

Theorem 1 If the optimal objective value of the LP is 0, then the density f de ned via (16) has uniform (0,1] marginals and covariance matrix .

Proof: Clearly, f is nonnegative and integrates to 1. Next, we need to show that the marginals of f are uniform. Let the marginal density function of X be denoted by g (
), which we have not yet proved is equal to f as desired. For any x 2 (y ;i? ; y ;i ), we have that 1

1

1

g1(x)dx = = =

n X

j;k=1 n X j;k=1 n X

1

1

P (X1 2 [x; x + dx)jX 2 C (i; j; k))P (X 2 C (i; j; k)) P (X1 2 [x; x + dx)jX1 2 (y1;i?1 ; y1;i ])q(i; j; k)

R y dx1 dy q(i; j; k) y ? j;k n X =1

= n

1

j;k=1

1;i 1;i 1

q(i; j; k)dx = 1dx

13

The rst equation follows by conditioning on the cell in which the random vector lies, and the second by the conditional independence of X1 ; X2 and X3 given that X lies in C (i; j; k). The third follows from the assumption of uniform \smearing" of q(i; j; k) on the cell C (i; j; k). A similar result holds for the marginals of X2 and X3 , and so the joint density f has the right marginals. Next we need to show that the obtained covariances are indeed the desired ones. Take the case of cov(X1 ; X2). Starting with its de niton, we have cov(X1 ; X2 ) = EX1 X2 ? EY1 EY2 = EY1 Y2 ? EY1 EY2 = 12 : The rst equality follows from the fact that EY1 = EX1 and EY2 = EX2 , and the second is just a restatement of (15). The nal equation follows from the fact that the optimal objective value is 0, as noted in (13). The same follows for cov(X2 ; X3 ) and cov(X1 ; X3 ). Hence, f has the covariances as desired and this completes the proof. 2

Remark 3 The name \chessboard" distribution is motivated by the form of (16) in a 2 dimensional

problem. In this case, the unit square is broken down in n2 squares, and the density f is constant on each square, with value n2 q(i; j ).

Remark 4 There is no need for the cells used in the above construction to be of equal size. Indeed, Theorem 1 remains true for more general discretizations; see Theorem 10 in Section 3.

In practice, this LP can become quite expensive to solve even for moderate values of n like 20, since we will then have 8000 q(i; j; k)s to contend with. However, the feasible region of the LP can be reduced by the inclusion of some more constraints on the Zij+ s and Zij? s. This not only improves the performance of the LP, but as we shall see, also provides us with a new feasibility criterion to test for the existence of a random vector with the given covariance matrix. The constraints are developed by assuming that a random vector X with uniform marginals and covariance matrix  exists, discretizing X to obtain a new random vector X~ say, and then bounding the change in the covariances resulting from the discretization. So suppose that we discretize X to obtain X~ . Let

q~(i; j; k) = P (X~ = (Y1;i ; Y2;j ; Y3;k ));

14

and observe that q~ provides a feasible solution to the above LP. We now wish to bound the change in the covariance resulting from this discretization. Observe that cov(X~1 ; X~2 ) ? 12 = E X~1 X~2 ? EX1 X2 n X = (Y1;i Y2;j ? E [X1 X2 jX 2 C (i; j; k)])~q (i; j; k):

(17)

y1;i?1 y2;j?1  E [X1 X2 jX 2 C (i; j; k)]  y1;i y2;j :

(18)

i;j;k=1

But Combining (17) with (18) we see that cov(X~1 ; X~2 ) ? 12  cov(X~1 ; X~2 ) ? 12 

n X i;j;k=1 n X i;j;k=1

q~(i; j; k)(Y1;i Y2;j ? y1;i?1 y2;j?1 ) and

(19)

q~(i; j; k)(Y1;i Y2;j ? y1;i y2;j ):

(20)

To summarize then, we assumed that X exists, and under this assumption have shown that there is a feasible solution q~ to the LP satisfying the bounds (19) and (20). Therefore, we can add these bounds + ?. to our LP. In particular, (19) gives an upper bound on Z12 , and (20) gives an upper bound on Z12 Similar bounds may be obtained for the other covariances. After substituting in the explicit expressions for yi;k and Yi;k , we see that these bounds simplify to

Zij+  21n ? 4n1 2

and Zij?  21n + 4n1 2 1  i < j  3:

(21)

Once the above LP is augmented with the bounds (21), it is no longer guaranteed to be feasible. In fact, Theorem 2 below establishes that if the augmented LP is infeasible for any value of n  1, then the covariance matrix  is not feasible for uniform marginals. The proof is basically a summary of the above discussion, and is given to help clarify these ideas.

Theorem 2 If the augmented LP is infeasible for some n  1, then there cannot exist a random vector X with uniform marginals and the desired covariance matrix .

Proof: Suppose there exists a random vector X with uniform marginals and covariance matrix . Then, as above, we can construct a solution q~ by discretizing X that satis es all of the constraints, including the bounds (21). Thus the augmented LP is feasible, which is a contradiction. 2 In fact, one can prove a converse to Theorem 2.

15

Theorem 3 If the covariance matrix  is not feasible for uniform (0,1] marginals, then there exists an n  1 such that the augmented LP is infeasible. Proof: On the contrary, suppose that the augmented LP is feasible for all n  1. Let qn denote an

optimal solution to the nth augmented LP, and let n denote the probability measure corresponding to the density resulting from the smearing operation (16) applied to qn . Then each n is the distribution of a random vector with support contained in (0; 1]3 with uniform(0; 1] marginals. Hence, the sequence (n : n  1) is tight, and by Theorem 29.3 on p. 392 of Billingsley (1986), it possesses a weakly convergent subsequence (n(k) : k  1), converging to  say. Now,  has uniform (0; 1] marginals. This follows from Theorem 29.2, p. 391 of Billingsley (1986) since each n(k) has uniform(0; 1] marginals, n(k) )  as k ! 1, and the projection map j : IR3 ! IR that returns the j th coordinate of a vector in IR3 is continuous. Now, if C n is the covariance matrix of the distribution qn , then n XX 2

i=1 j =i+1

jCijn ? ij j  23n + 4n3 ! 0 2

as n ! 1. This follows from the bounds (21), and the fact that in any optimal solution, it is not the case that both Zij+ and Zij? are strictly positive. Finally, if X n(k) has distribution n(k) , then (Xin(k) Xjn(k) : k  1) is a bounded sequence of random variables, and therefore uniformly integrable. It immediately follows that the covariance matrix  of  is given by  = klim C n(k) = ij : !1 Thus,  has the required marginals and covariance matrix, which is a contradiction, and the result is proved. 2

Combining Theorems 2 and Theorem 3, we see that a covariance matrix is infeasible for uniform marginals if, and only if, the augmented LP is infeasible for some n  1. Given this very sharp characterization of infeasible covariance matrices, it is natural to ask whether a similar result holds for feasible covariance matrices. We would then have the result that a covariance matrix is feasible for a given set of marginals if and only if there is some nite n such that the optimal objective value of the augmented LP is zero. Unfortunately, this conjecture is false, as shown by the following counterexample in 2 dimensions. Suppose that X1 = X2 and hence cov(X1 ; X2 ) = var(X1 ) = 1=12. For given n, the covariance between Y1 and Y2 is maximized by concentrating all mass on the cells (i; i), and so q(i; i) = n?1 for

16

1  i  n. In that case, we have that cov(Y1 ; Y2 ) =

n  2i ? 1  X i=1

2n

2

1 ? 11 = 1 ? 1 n 2 2 12 12n2

Therefore, cov(Y1 ; Y2 ) < 1=12 for all nite n, and so the conjecture is false. Notice that the covariance matrix in this example is singular. This example is a special case of the following result.

Theorem 4 All chessboard densities have nonsingular covariance matrices. Proof: On the contrary, suppose that f is a chessboard density with singular covariance matrix , and

let X have density f . Since  is singular, there exists a nonzero vector  such that  = 0. Hence, var(0 X ) = 0  = 0, and so 0 X = 0 EX a.s. Since  is nonzero, we may, by relabelling variables if necessary, write X1 as a linear function of X2 ; X3 , say X1 = 0 + 2 X2 + 3 X3 . This equality must also hold conditional on X 2 C (i; j; k). But the components of X are conditionally independent given that X 2 C (i; j; k) because f is a chessboard density, which is the required contradiction. 2 The importance of Theorem 4 is that if  is feasible for the given marginals and singular, then no matter how large n may be, the optimal objective value of the LP will always be > 0, i.e., we cannot exactly match the covariance matrix . However, we can come arbitrarily close, as the following result shows.

Theorem 5 Suppose that the covariance matrix  is feasible for uniform (0; 1] marginals. Then for all n  1, the augmented LP is feasible, and if z (n) is the optimal objective value of the nth LP, then z (n) ! 0 as n ! 1. Proof: Since  is feasible for uniform marginals, the augmented LP is feasible for all n  1. (This is

just the contrapositive of Theorem 2.) Let qn denote an optimal solution to the nth LP, and let f n be the corresponding smeared density. If C n is the covariance matrix corresponding to f n , then the bounds (21) imply that 2 n X X jCijn ? ij j  23n + 4n3 2 ! 0 z (n) = i=1 j =i+1 as n ! 1. 2 Therefore, chessboard densities can come arbitrarily close to any required  that is feasible for uniform marginals. In fact, one can prove that chessboard densities can exactly match a (very) slightly restricted class of feasible covariance matrices. To state this result we need some notation.

17

We can and do easily state and prove Proposition 6 for a general dimension d (i.e., not just d = 3) without any notational diculty. Any covariance matrix  of a d dimensional random vector with uniform[0; 1) marginals can be characterized by d(d ? 1)=2 covariances, since the diagonal entries are determined by the marginals, and the matrix is symmetric. Hence we can, with an abuse of notation, think of  as a d(d ? 1)=2 dimensional vector in some contexts, and as a d  d matrix in others. Let  [?1=12; 1=12]d(d?1)=2 denote the space of feasible covariance matrices, so that  2 implies that there exists a random vector with uniform(0; 1] marginals, and covariance matrix . We will show below that is nonempty and convex (this is well-known), but also closed and full-dimensional (this appears to be new). In particular then, any covariance matrix on the boundary of is feasible. We will also show that  is contained in the interior of if, and only if, there is some nite n for which the augmented LP has objective value 0. The collective implications of this and our previous results will be discussed after the statement and proof of these results.

Proposition 6 The set is nonempty, convex, closed and full-dimensional. Proof: If the components of X are independent, then the covariance matrix  is diagonal, and so

contains the zero vector, and is therefore nonempty. It is well-known that is convex. For if 1 ; 2 2 , then there exist random vectors X; Y with uniform(0; 1] marginals, and covariance matrices 1 and 2 respectively. For  2 (0; 1), let Z be given by X with probability , and Y with probability 1 ? . Then Z has covariance matrix 1 + (1 ? )2 . The proof that is closed is virtually identical to that of Theorem 3 and is omitted. We use the NORTA method to prove that is full-dimensional. We will show that each of the vectors ek =12 are contained in , where ek is the vector whose components are all 0 except for a 1 in the kth position, for k = 1; : : : ; d(d ? 1)=2. The convexity of then ensures that is full-dimensional. Let Z be a multivariate normal random vector with mean 0 and covariance matrix consisting of 1's on the diagonal, and also in the (i; j )th and (j; i)th position (i 6= j ), with the remaining components being 0. That is, Z consists of 2 perfectly correlated standard normal random variables Zi and Zj , and d ? 2 independent standard normal random variables. Now let U be the random vector with uniform (0; 1) marginals obtained by setting Um = (Zm ) for m = 1; : : : ; d. Then Ui and Uj are perfectly correlated, and independent of all of the remaining components of U . Thus, U has covariance matrix whose components are all 0 except for the diagonal elements, and the (i; j ), and (j; i)th elements, which are equal to 1=12. Thus, ek =12 lies in , where k corresponds to the position (i; j ). A similar argument with perfectly negatively correlated Zi and Zj shows that ?ek =12 2 . Since i 6= j were arbitrary, the

18

proof is complete. 2 In Theorem 4 we showed that all chessboard densities have nonsingular covariance matrices. This is almost sucient to establish that all boundary points of do not have chessboard densities. However, it is certainly conceivable that the boundary of contains nonsingular, as well as singular, covariance matrices. So we strengthen Theorem 4 with the following result.

Theorem 7 If f n is a chessboard density with covariance matrix , then f n is contained in the interior of .

Proof: Let X have density f n. We will show that we can both increase, and decrease, the covariance

between X1 and X2 . Symmetry then allows us to conclude that the same result holds for Xi and Xj with i 6= j . The convexity of then completes the proof. Let qn be the discretization of f n into its n3 cells, and let C (i; j; k) be a cell with qn (i; j; k) > 0. Now, divide the cell C (i; j; k) into 4 (equal size) subcells,

Cab (i; j; k) = f(x; y; z ) 2 C (i; j; k) : 2i ? 2(3n ? a) < x  2i ? 2(2n ? a) ; 2j ? (3 ? b) < y  2j ? (2 ? b) g; 2n 2n for 1  a; b  2. Now, generate a new density gn by the usual smearing (16) in all cells except C (i; j; k). Within the cell C (i; j; k), assign a mass of qn (i; j; k)=2 to each of the cells C11 (i; j; k), and C22 (i; j; k), and then uniformly smear within these cells. In other words, for (x; y; z ) contained in these two cells, set gn (x; y; z ) = 2n3qn (i; j; k) and set gn to be 0 in the cells Cab (i; j; k) for a 6= b. Then it is straightforward to show that gn has uniform marginals, that the (1; 2)th covariance is strictly increased, and that the other covariances remain unchanged. A similar argument placing the mass in the cells Cab (i; j; k) with a 6= b shows that the covariance can be strictly decreased, and so the proof is complete. 2 We have thus far shown that if a covariance matrix  is not in , then the augmented LP will be infeasible for some n  1, and if  is on the boundary of , then the LP approach will yield distributions with covariance matrices that arbitrarily closely approximate , but never actually achieve it. Our nal result shows that if  is contained in the interior of , then there is some n  1 for which the optimal objective value of the augmented LP is 0, and so one can exactly match  using a chessboard density. Before proving this result, we need the following lemma. This lemma basically states that given a xed

19

vector x, we can choose certain other vectors arbitrarily close to x, so that x is a convex combination of these \close" vectors, and if we perturb the close vectors slightly, then x is still a convex combination of the perturbed vectors. For x 2 IRm and  > 0, let B (x; ) denote the (open) set of vectors

fy 2 IRm : (x; y) < g; where  is the L1 distance

(x; y) =

m X i=1

jxi ? yi j:

Lemma 8 Let x 2 IRm, and let  > 0 be arbitrary. There exist m + 1 points x ; : : : ; xm 2 B(x; ), 1

and a  > 0 such that if

+1

(xi ; x0i ) <  8i = 1; : : : ; m + 1;

then x may be written as a convex combination of x01 ; : : : ; x0m+1 .

Proof: Basically, one chooses the xi 's to be the vertices of a simplex centered at x. To be precise, let

r > 0 be a parameter, and set x1 x2 x3 .. .

xm xm+1 where

( ?a1 ?a2 ( a1 ?a2 ( 0 2a2 .. .. .. . . . = ( 0 0 = ( 0 0

= = = .. .










.. .

?am? ?am? ?am?

1 1 1

?am )0 + x ?am )0 + x ?am )0 + x

.. .. . .


(m ? 1)am?1 ?am


0 mam

r m s ai = r

.. . )0 + x )0 + x;

1

m + 1 i(i + 1) :

Then, (Dantzig 1991), the xi 's de ne the vertices of an equilateral simplex whose center is x, and whose vertices are a (Euclidean) distance rm=(m + 1) from x. Choose r so that xi 2 B (x; ) for all i. Observe that the average of the xi 's is x. In fact, it is easy to show that the (m +1)  (m +1) matrix B consisting of the xi 's in columns, supplemented with a row of 1's is nonsingular, and so

y = B ?1 x = (m + 1)?1 (1; 1; : : : ; 1)0 : Now, observe that B ?1 is a continuous function of B , at least in a neighbourhood of B , and so y = B ?1 x is locally a continuous function of x1 ; : : : ; xm+1 . Hence, there is a  > 0 such that if (xi ; x0i ) <  for

20

all i = 1; : : : ; m + 1, and D consists of the vectors x0i in columns supplemented with a row of 1's, then y = D?1 x consists of all positive components, and the elements of y sum to 1. 2 We are now ready to state the nal result of this section. As in Proposition 6, there is no loss of clarity if we state this result for a general dimension d rather than just d = 3.

Theorem 9 If  is contained in the interior of , then there exists an n  1 such that the optimal objective value of the augmented LP is 0.

Proof: Let m = d(d ? 1)=2, and for now, consider  as an m-vector. Let  > 0 be such that B(; )  , and choose  ;  ; : : : ; m 2 B (; ) and  as in Lemma 8. Since i 2 , from Theorem 5 there exists an n(i) such that the augmented LP with target 1

2

+1

covariance matrix i has optimal objective value smaller than , for each i = 1; : : : ; m + 1. Let n = n(1)n(2)


n(m + 1), and let qi denote a solution to the augmented LP with target matrix i and discretization level n for i = 1; : : : ; m + 1. Then the optimal objective value corresponding to qi is also less than . (Note that if k; n  1 are integers, then the optimal objective values z (n) and z (kn) satisfy the relationship z (kn)  z (n), since the chessboard density obtained from the solution to the nth LP can also be obtained from the (kn)th LP.) Let 0i denote the covariance matrix corresponding to the chessboard density f i for the solution qi , for i = 1; : : : ; m + 1. Then, by Lemma 8, there exist nonnegative multipliers 1 ; 2 ; : : : ; m+1 summing to 1 such that mX +1 i 0i : (22) = i=1

If we set

f=

mX +1 i=1

i f i ;

then f is also a chessboard density with discretization level n, and from (22), its covariance matrix is exactly . 2 In summary then, we have shown that if  is infeasible for uniform marginals, then the augmented LP will be infeasible for some n  1. This includes the case where  is singular and infeasible for uniform marginals. Furthermore, we have shown that if  is contained in the interior of , then the augmented LP will have optimal objective value 0 for some n  1, and so one can construct a chessboard density from the solution to the augmented LP with the required marginals and covariance matrix. So if  is not contained in the boundary of , then we have an algorithm for determining, in nite time, whether

21

 is feasible for the given marginals or not. One simply solves the augmented LP for n = 1; 2; 3; : : : until the augmented LP is either infeasible, or has an optimal objective value of 0. In the latter case, we can deliver an explicit construction of the desired distribution. The case where  lies on the boundary of is more problematical. We have shown that in this case,  is feasible for uniform marginals, but that a chessboard density cannot be constructed with uniform marginals and covariance matrix . Therefore, for such matrices, the algorithm outlined above will not terminate in nite time. However, a chessboard distribution can come arbitrarily close to the required covariance matrix. These results shed a great deal of light on the class of feasible covariance matrices for random vectors with uniform marginals. It is natural to ask whether a similar approach may be applied to a more general class of marginal distributions, and this is the subject of the following section.

3 Joint Distributions with more general marginals The LP method used in Section 2 to evaluate the existence of a random vector with uniform marginals and given covariance matrix can be adapted to investigate the existence of random vectors having arbitrary marginal distributions and given Pearson product-moment covariance matrix. We will stick to the case of a 3-dimensional random vector for notational simplicity but it should be noted that the approach is easily extended to the general d-dimensional case. Let X = (X1 ; X2 ; X3) represent the random vector that is to be constructed, and let  be the desired covariance matrix. Let Fi (
) denote the distribution function of Xi , for i = 1; 2; 3. For ease of exposition we assume that each of the Fi 's has a density fi with respect to Lebesgue measure, although the approach applies more generally, and in particular, can be applied when some or all of the Xi 's have discrete distributions. In the spirit of the method developed in Section 2, we will rst construct the probability mass function of a discretized random vector Y = (Y1 ; Y2 ; Y3 ) with a covariance structure as close to the desired one as possible, and then derive a joint distribution function for X . Let n1 ; n2 and n3 represent the levels of discretization of the random variables X1 ; X2 and X3 respectively, and hence the number of points that form the support of Y1 ; Y2 and Y3 . Let the range of the variable Xi be divided into ni subintervals (which may, or may not, be equal in length) by the set of points fyi;0 ; yi;1 ; : : : ; yi;n g, with i

?1  yi; < yi; <


< yi;n  1: 0

1

22

i

Note that we explicitly allow yi;0 and yi;n to be in nite and the spacing between the yi;k s to be arbitrary. Let Yi;k denote the conditional mean of Xi , given that it lies in the subinterval (yi;k?1 ; yi;k ]. In other words, we set Z y fi(x) dx; x Yi;k = E [Xi jXi 2 (yi;k?1 ; yi;k ]] = y ?1 Pi (k ) where Pi (k) = Fi (yik ) ? Fi (yi;k?1 ) represents the probability that Xi lies in the kth subinterval. The support for the discretized random variable Y is then fY1;i ; Y2;j ; Y3;k : 1  i  n1 ; 1  j  n2 ; 1  k  n3 g. Let q(i; j; k) represent P (Y = (Y1;i ; Y2;j ; Y3;k )), i.e., the probability that X 2 C (i; j; k), where C (i; j; k) is de ned as in Section 2 to be the cell corresponding to q(i; j; k). We now give constraints on the q(i; j; k)s analogous to (8) through (11). Speci cally, we have that i

i;k

i;k

n X j;k=1 n X i;k=1 n X i;j =1

q(i; j; k) = P (Y1 = Y1;i ) = P (X1 2 (y1;i?1 ; y1;i ]) = P1 (i); 8i = 1; : : : ; n1 ;

(23)

q(i; j; k) = P (Y2 = Y2;j ) = P (X2 2 (y2;j?1 ; y2;j ]) = P2 (j ); 8j = 1; : : : ; n2 ;

(24)

q(i; j; k) = P (Y3 = Y3;k ) = P (X3 2 (y3;k?1 ; y3;k ]) = P3 (k); 8k = 1; : : : ; n3 ;

(25)

q(i; j; k)  0 8i = 1; : : : ; n1 ; j = 1; : : : ; n2 ; k = 1; : : : ; n3 :

(26)

When these constraints are satis ed, we have that EYi = EXi for each i, just as in the case of uniform marginals. We can now formulate an LP along the lines of that given in Section 2 to attempt to match the covariances of Y to those required of X . If Cij is de ned as in Section 2, then the formulation of the LP is min

P P 2 =1

i

+ ? j i (Zij + Zij ) 3 = +1

subject to Cij ? ij = Zij+ ? Zij? ; i = 1 to 2 and j = i + 1 to 3

Zij+  0; Zij?  0; together with constraints (23), (24), (25) and (26).

This LP is always feasible, since we can take the Yi 's to be independent, and the objective value is bounded below by 0. Hence, an optimal solution exists for every discretization. If the optimal objective value is 0, then we have been able to construct a probability mass function for Y with the required covariance structure. Constructing a joint distribution function for X from this pmf for Y is similar in avor to the method of uniform \smearing" used in Section 2.

23

Speci cally, the \smearing" process should be able to satisfy (15). Again, the easiest method for doing so is perhaps to ensure that the variables are conditionally independent, given that X lies within the cell. To ensure that this conditional independence holds while simultaneously ensuring that X has the right marginals, we set the density of X within the cell C (i; j; k) to be that of independent random variables with the right marginals, scaled so that the total mass in the cell is q(i; j; k). To be precise, if f is the density of X , then for any x = (x1 ; x2 ; x3 ) 2 C (i; j; k), we set

f (x) = fP1 (x(i1)) fP2 ((xj2)) fP3((xk3)) q(i; j; k): 1

2

(27)

3

We can now prove an analogous result to Theorem 1, which states that if the optimal objective value of the LP is 0, then the density f constructed using (27) has the desired marginals and covariance matrix. The proof is virtually identical to that of Theorem 1 and is omitted.

Theorem 10 If the optimal objective value of the LP is 0, then the density f de ned via (27) has the required marginals and covariance matrix .

We can also obtain analogues to the other results of Section 2. As in Section 2, let  IRd(d?1)=2 denote the set of feasible covariance matrices for the given marginals. Note that, as before, we think of a given covariance matrix as a vector in d(d ? 1)=2 dimensional space in some contexts, and as a d  d matrix in others. We have the following result on the structure of .

Proposition 11 The set is nonempty, convex and full-dimensional. The proof of this result is virtually identical to that of Proposition 6 and is omitted. We also have the following analogue of Theorems 4 and 7.

Theorem 12 Any chessboard density has a nonsingular covariance matrix. Furthermore, if f n is a chessboard density with covariance matrix , then  is contained in the interior of .

Again, the proofs are very similar to those in the previous section and omitted. To extend the other results of the previous section, we assume that yi;0 and yi;n are nite for all i, i.e., that all of the distribution functions Fi have bounded support. We further assume that all of the ni 's are equal to n say, and that all subintervals are of equal length. Thus, we will discretize on a regular grid containing n3 cells. In particular, suppose that a random vector X with the desired marginals and covariance matrix exists, and let X~ denote its discretization. Let q~(i; j; k) be the probability that X lies in the cell C (i; j; k). We can now bound the dierence in the covariance between X1 and X2 . But rst it is convenient to let i

24

ai = yi;0 and
i = yi;1 ? yi;0 , so that
i is the width of the cells in the ith coordinate direction, for i = 1; 2; 3. With this notation, we have that cov(X~1 ; X~2 ) ? 12 =

 = =

n X

i;j;k=1 n X i;j;k=1 n X i;j;k=1 n X i;j;k=1

q~(i; j; k)[Y1;i Y2;j ? E [X1 X2 jX 2 C (i; j; k)]] q~(i; j; k)[y1;i y2;j ? y1;i?1 y2;j?1 ] q~(i; j; k)[(a1 + i
1 )(a2 + j
2 ) ? (a1 + (i ? 1)
1 )(a2 + (j ? 1)
2 )] q~(i; j; k)[
1 (a2 + (j ? 1)
2 ) +
2 (a1 + (i ? 1)
1 ) +
1
2 ]


EX +
EX +

: 1

2

2

1

1

2

A similar lower bound can be derived, and so the LP can be augmented by the bounds

Zij+ ; Zij? 
i EXj +
j EXi +
i
j ; for 1  i < j  3. Observe that as n ! 1, these bounds converge to 0. This is the nal ingredient required to strengthen the other results of the previous section to the more general case of distributions with bounded support and densities. In particular, we now have the following results, which we state without proof because the proofs are similar to the case of uniform marginals. The implications of these results are given after the statement of the results.

Remark 5 The above bounds were derived assuming a regularly spaced grid of cells, so that the cells were all of identical size. However, similar bounds can be expected to hold when the cells are not of equal size. Indeed, one should be able to obtain bounds on Zij+ and Zij? which converge to 0 as long as the maximum sidelength of the cells converges to 0.

Proposition 13 Suppose that  is feasible for the given marginals. If all of the densities fi have bounded support, then as n ! 1, the optimal objective value of the LP converges to 0. Furthermore, the set is closed.

Theorem 14 Suppose that all of the densities fi have bounded support. Then 9n  1 such that the nth LP is infeasible if and only if the matrix  is infeasible for the given marginals.

25

Theorem 15 Suppose that all of the densities fi have bounded support. Then 9n  1 such that that the optimal objective value of the nth LP is 0 if and only if the matrix  is contained in the interior of .

So assuming that all of the densities fi have bounded support, and  does not lie on the boundary of

, then we have a nite algorithm for determining whether  is feasible or not, and if feasible, supplying an explicit joint density with the required properties. The algorithm is simply to solve a sequence of LPs for n = 1; 2; : : : until either the LP is infeasible, or has an optimal objective value of 0. If  lies on the boundary of , then we know that it is feasible, and we can approach it arbitrarily closely with chessboard distributions, but never exactly reach it. It is natural to ask how these results might be generalized to densities with in nite support. The primary diculty seems to lie in coming up with a satisfactory bound on the objective value of the augmented LP that converges to 0 in some sense, and we leave this as an open problem. In the case of distributions with in nite support then, the utility of this approach is limited to generating chessboard distributions for covariance matrices that lie in the interior of .

4 Unit Cube Exploration The characterizations of Section 3 are certainly potentially useful, especially given the wide use of product-moment correlation in practice. But it is a widely held contention (and one that we share) that rank covariance is a more meaningful measure of dependence than product-moment covariance when dealing with nonnormal distributions. Therefore, while we view the results of the previous section as useful and interesting, we prefer to emphasize the case where one is attempting to induce a given rank correlation matrix  over a set of marginal distributions. In the case where all of the marginal distribution functions are continuous, this problem reduces to that of inducing the required product-moment covariance matrix over uniform marginals; see (6), noting that (Zi ) has a uniform distribution on (0; 1). Therefore, it is important to gain some understanding of the covariance structure of uniform random variables, and in particular, to gain some understanding of the class of covariance matrices that can be achieved using the NORTA method. We present the results we have obtained by applying the LP method introduced in Section 2 to construct joint distributions for a random vector X = (X1 ; X2 ; X3 ) with uniform marginals and covariance

26

matrix . For notational convenience, we will use the correlation matrix

1 0 1   CC BB B R=B  1  C CA @ 12

12

13

23

13 23

1 instead of the covariance matrix  in the rest of this section. The matrix R is characterized by 12 ; 13 and 23 . The potential set  of all possible values of  = (12 ; 13 ; 23 ) that can constitute correlation matrices is just a rescaling of the set in Section 2, and is a proper subset of the cube [?1; 1]3 because a correlation matrix R is constrained to be positive semide nite. We examined all correlation matrices whose components 12 ; 13 and 23 lay on the grid f-1.0, 0.9, : : :, -0.1, 0, 0.1, : : :, 0.9, 1.0g3. Each matrix R was rst checked for positive semide niteness and discarded if not, thus leaving us 4897 matrices from the 9261 matrices that were tested. The remaining 4897 positive semide nite matrices were further tested to see whether they were NORTA feasible. We de ne the matrix R to be NORTA feasible if the covariance matrix  found via (4) is positive semide nite. In this case, a multivariate normal random vector with covariance matrix  will be transformed via the NORTA method to a multivariate uniform random vector with the required correlations . In this manner, a total of 160 sample matrices were identi ed to be NORTA defective. Note that since X1 ; X2 and X3 are identically distributed, many dierent 's form the same eective correlation matrix for X . For example,  =(0.5,-0.5,0.5), (-0.5,0.5,0.5) and (0.5,0.5,-0.5) constitute the same joint distribution for X up to a symmetry. If we eliminate such multiple occurences, the number of NORTA defective matrices reduces to 31 cases. The question that remains to be answered is whether these NORTA defective matrices are feasible for uniform marginals. This is where the LP method of Section 2 comes in handy. We apply the LP method to each NORTA defective case iteratively for increasing values of n, the level of discretization, to determine whether a chessboard density can be constructed in each case. For the sake of computational tractability, we limited the value of n to a maximum of 80. The results obtained have been tabulated in Table 1. Table 1 indicates that the LP method was able to construct a chess-board distribution for almost all of the cases for low to moderate values of n. In fact, cases 26 to 31 are the only cases where the method failed to tie down a chessboard density for X with the given correlation matrix R. This is expected from Theorem 4 since each of these R matrices are singular.

27

Correlations

No.

12

13

23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

-0.9 -0.9 -0.9 -0.9 -0.8 -0.8 -0.8 -0.8 -0.8 -0.8 -0.8 -0.8 -0.7 -0.7 -0.6 -0.5 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 -0.8 -0.8 -0.6 -0.5 -0.5 0.0

-0.6 -0.5 -0.2 -0.1 -0.8 -0.5 -0.4 -0.3 -0.3 0.1 0.2 0.3 -0.7 0.0 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 0.7 0.5 0.6 0.8 -0.6 0.0 0.0 -0.5 0.5 0.6

0.2 0.1 0.6 0.5 0.3 -0.1 -0.2 -0.3 0.8 0.5 0.4 0.3 0.0 0.7 0.9 0.9 0.8 0.8 0.8 0.8 0.8 0.7 0.9 0.9 0.8 0.0 0.6 0.8 -0.5 0.5 0.8

Final Optimal Level of Discretization Determinant Smallest Eigenvalue Objective Value n R of R j

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000159 0.0000159 0.0000159 0.0000159 0.0000159 0.0000159

18 11 18 11 12 11 10 10 12 11 10 10 12 12 18 11 11 10 10 10 11 12 11 18 12 80 80 80 80 80 80

j

0.0060 0.0200 0.0060 0.0200 0.0140 0.0200 0.0320 0.0360 0.0140 0.0200 0.0320 0.0360 0.0200 0.0200 0.0060 0.0200 0.0200 0.0320 0.0360 0.0320 0.0200 0.0200 0.0200 0.0060 0.0140 0 0 0 0 0 0

0.0034 0.0105 0.0034 0.0105 0.0087 0.0097 0.0151 0.0169 0.0087 0.0097 0.0151 0.0169 0.0101 0.0101 0.0034 0.0105 0.0097 0.0151 0.0169 0.0151 0.0097 0.0101 0.0105 0.0034 0.0087 0 0 0 0 0 0

Table 1: Results of the LP method when run for the NORTA defective R matrices

28

Hence, almost all of these NORTA defective matrices are feasible for uniform marginals. In particular, case 18 is the example given by Li and Hammond (1975). So we have rigorously established that there are feasible covariance matrices for a given set of marginals which cannot be generated via the NORTA method. It is interesting to question why the NORTA defective matrices are, in fact, NORTA defective. An inspection of the NORTA defective matrices in Table 1 shows that the determinants and the smallest eigenvalues of all these matrices are quite close to zero in magnitude. This means that the corresponding points  lie either on the boundary of  or in its close proximity, i.e., the NORTA defective matrices lie close to, or on, the boundary of the set of achievable correlations. Indeed, we can expect that this will be the case for more general distributions, and higher dimensional problems. We reason heuristically as follows. The set Z of feasible correlation matrices for a multivariate standard normal random vector is the set of all positive semide nite matrices with unit diagonal entries. The NORTA method transforms each element of this set into an element of . Therefore, the set Z is mapped into a subset of . Assuming the transformation is continuous, which indeed it is in great generality (see Cario and Nelson 1997), and not \too nonlinear", it is reasonable to expect that any elements of  that are not covered by the transformation of Z will be those that are close to the boundary of . Indeed, if one plots the set of NORTA defective vectors using a three dimensional plotting package, this is exactly what we see. The numbers quoted above seem to suggest that the occurence of NORTA defective R matrices is relatively rare. However, Lurie and Goldberg (1998) believe that singular and near singular correlation matrices actually represent a common situation in cost analysis for example. This is because correlations between cost elements are typically estimated from unbalanced data sets. This is likely to lead to inde nite target correlation matrices, so that any least adjustment to them is almost certainly going to result in an adjusted target matrix that is singular, or very nearly so. So then, we can either modify the NORTA method to generate random vectors with the required marginals and approximately (in NORTA defective cases) the right covariance matrix, or turn to some other method. In the following sections we describe such a modi cation to the initialization phase of the NORTA method, and outline an alternative approach to the NORTA method that relies on the generation of iterates with a chessboard distribution.

29

5 Semide nite Programming and NORTA As seen in the previous section, the Li and Hammond counterexample does exist, as do others, and so there are sets of marginals with feasible (product-moment or rank) correlation matrices that cannot be exactly matched using the NORTA method. Nevertheless, the NORTA method is very general, and quite easily implemented, and so we may wish to employ the method to generate random vectors with the required marginals, and at least approximately, the right covariance matrix. Lurie and Goldberg (1998) described a method for identifying a positive semide nite covariance matrix Z for use within the NORTA method, that yields approximately the desired product-moment covariance matrix X . Their approach involves a complicated nonlinear optimization, and must be specialized for approximating the rank correlation or product-moment correlation, depending on the case desired. In contrast, the method we present below does not require such specialization. Furthermore, although they report that their optimization procedure always converges in practice, they do not have a proof of this result. The technique below does not share such a problem because we are using standard semide nite programming techniques. Finally, their suggested technique appears to be limited to xed sample sizes, whereas our adjustment to the initialization phase of NORTA does not share this limitation. Let Z be the covariance matrix that we wish to use in the NORTA procedure. In what follows, we do not distinguish between the cases where Z is chosen to induce a given rank, product-moment, or other correlation in the output random vector X . Our procedure uses semide nite programming (SDP) to search for a positive semide nite matrix Z that is \close" (in a certain sense) to the matrix Z . We will then employ Z within the NORTA method to generate the required random vectors. Our formulation then falls under the broad class of matrix completion problems; see Alfakih and Wolkowicz (2000), or Johnson (1990). Why is this approach reasonable? In Theorem 2 of Cario and Nelson (1997), it is shown that under a certain moment condition, the target covariance matrix X is a continuous function of the input covariance matrix Z used in the NORTA procedure. This moment condition always holds when we are attempting to match rank covariances, and we can expect it to hold almost invariably when matching product-moment correlations. Therefore, it is eminently reasonable to try and minimize some measure of distance d(Z ; Z ) say, between Z and Z . Given Z as data, and assuming that we are operating in dimension d = 3, we wish to choose a

30

symmetric matrix Z to minimize jZ (1; 2) ? Z (1; 2)j + jZ (1; 3) ? Z (1; 3)j + jZ (2; 3) ? Z (2; 3)j subject to Z  0; Z (i; i) = 1; where the matrix inequality A  0 signi es a constraint that the matrix A be positive semide nite. We can use the same trick used in formulating the LPs in Section 2 and write each component of the objective function above as jZ (i; j ) ? Z (i; j )j = wij+ + wij? ; where wij+ and wij? are the positive and negative parts respectively of the dierence Z (i; j ) ? Z (i; j ), with wij+ ; wij?  0. In particular, we have that for 1  i < j  3, Z (i; j ) = Z (i; j ) + wij+ ? wij? . For a vector y, let diag(y) denote the diagonal matrix with diagonal entries equal to the elements of y. It is well known in SDP formulations that any set of linear inequalities of the form Ax + b  0 can be transformed to a matrix inequality of the form diag(Ax + b)  0. The nonnegativity constraints on the variables wij+ and wij? can be easily handled in this manner. The complete SDP formulation of the problem in the case of a 3 dimensional problem is then to min s/t

+ ? + w+ + w? + w+ + w? w12 + w12 13 13 23 23

2 4 Z 0

2 1 3 66 0 5 66 Z (1; 2) + w ? w? =6 64 Z (1; 3) + w ? w? D + 12

12

+ 13

13

+ Z (1; 2) + w12 ? w12? Z (1; 3) + w13+ ? w13? + 1 Z (2; 3) + w23 ? w23? + Z (2; 3) + w23 ? w23? 1

D

3 77 77  0; 77 5

+ ? ; w+ ; w? ; w+ ; w? g. The matrix we where D is the diagonal matrix with diagonal elements fw12 ; w12 13 13 23 23 desire is an optimal solution Z to this SDP. There is considerable exibility in formulating the problem using this SDP framework that allows us to include preferences on how the search for Z is performed. For example, if we believe that Z (i; j ) should not be smaller than Z (i; j ), we simply omit the variable wij? from the formulation. Alternatively, if we believe that the value Z (i; j ) should not be changed, we simply omit the variables wij+ and wij? . Alternatively, if we believe that Z (i; j ) should not be changed by more than  > 0 say, then we add the bounds wij+   and wij?   in the usual (diagonal) fashion. Ecient algorithms are available for solving semide nite problems; see Wolkowicz, Saigal and Vandenberghe (2000). All of the semide nite programs in this paper were solved using the CSDP solver, available at the NEOS Optimization Server (2000) website.

31

We formulated SDPs as shown above for all of the 31 cases of NORTA defective correlation matrices that we identi ed in Section 4 and solved them using the CSDP solver. The optimal objective values d(Z ; Z ) obtained for each of these cases is shown in Table 2. The covariance matrix 0X for the vector X that can be achieved using Z can be calculated from (4). The values calculated for d(0X ; X ) are also shown in Table 2. Now, since the SDP essentially solves a convex optimization problem, the optimal Z obtained is guaranteed to be the nearest we can get to the desired matrix Z . In fact, Z is a projection of Z onto the set of positive semide nite matrices Z with diagonal elements equal to 1. It follows that Z lies on the boundary of Z and hence is singular, which indeed is true of all the 31 cases we solved. It does not immediately follow though, that the induced covariance matrix 0X is singular, since the NORTA transformation alters covariances in a nonlinear fashion. It should be pointed out, in fact, that 0X may not be the closest NORTA feasible covariance matrix to X , because the optimization was performed \in Gaussian space". This is in contrast to the Lurie and Goldberg (1998) procedure. But the values shown in Table 2 seem to suggest that the dierence d(0X ; X ) is usually very small. In summary then, we are simply proposing an additional step in the initialization phase of the NORTA method. The revised NORTA procedure is as follows. 1. Identify Z to match some aspect of the desired covariance structure of X in any fashion. 2. If Z is positive semide nite, then one can proceed directly with the NORTA procedure. 3. If not, solve an SDP of the form outlined above to identify a matrix Z that is \close" to Z . 4. Use the matrix Z in the NORTA procedure. The additional work involved in this modi cation only shows up in the initialization phase of the NORTA method, and so there is no additional computational overhead while the method is being used to generate replicates of X . Furthermore, public-domain algorithms are available for solving SDPs, and these algorithms can handle very large dimensional problems with relative ease. Finally, the method does not require tailoring to dierent correlation measures. The only step that depends on the correlation measure being used is the rst step when Z is identi ed.

32

Case No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 -

Optimal dierences Case d(Z ; Z ) d(X ; X ) No.

Optimal dierences d(Z ; Z ) d(X ; X )

0.009997 0.002898 0.009997 0.002899 0.006534 0.016755 0.011054 0.009247 0.006534 0.016755 0.011054 0.009247 0.019258 0.019258 0.009997 -

0.002899 0.016755 0.011054 0.009247 0.011054 0.016755 0.019258 0.002899 0.009997 0.006534 0.027322 0.027322 0.027322 0.052914 0.052914 0.027322

0

0.010699 0.003106 0.010699 0.003106 0.006827 0.017479 0.011539 0.009655 0.006827 0.017479 0.011539 0.009655 0.019679 0.019679 0.010699 -

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

0

0.003106 0.017479 0.011539 0.009655 0.011539 0.017479 0.019679 0.003106 0.010699 0.006827 0.028466 0.028466 0.028466 0.052249 0.052249 0.028466

Table 2: Results of the SDP when run for the NORTA defective R matrices

33

6 Generating From Chessboard Distributions Suppose, for convenience, that all of the marginals of X have continuous distribution functions, and that we wish to match a given (feasible) rank correlation matrix X . If the corresponding Z given by (4) is not positive semide nite, then the approximate method given in the previous section may be used to produce random vectors with approximately the desired rank covariance matrix. It may be the case, however, that one is not willing to live with such an approximation. An alternative method based on chessboard distributions described below could be used in such cases. Furthermore, the chessboard method could be competitive (from a numerical eciency standpoint) with the NORTA method in many cases, and so it is certainly worth exploring such an approach in this paper. Our assumption of continuous marginal distribution functions ensures that the given rank covariance matrix X can be matched by rst generating a random vector U of uniform(0; 1] random variables with Pearson product-moment covariance X , and then setting Xi = Fi?1 (Ui ) for i = 1; : : : ; d. Therefore, the main diculty lies in generating the random vector U . Since X is feasible for the given marginals, it is also feasible for U . Therefore, X lies in the set of feasible covariance matrices for uniform marginals. If X does not lie on the boundary of , then as in Section 2, we can identify a chessboard density that exactly matches X . If X lies on the boundary of , then we can get arbitrarily close to X , but never exactly match it. In any case, suppose that we have determined an acceptable chessboard density f n for U with uniform marginals, and covariance matrix X . How can we generate from this density? One approach is to rst generate the cell in which U will reside. Then, conditional on U lying in the selected cell, we generate the components of U independently. These components are simply uniform random variables on a range that depends on the cell. So the approach is as follows. 1. Formulate the problem of constructing a chessboard distribution as an augmented LP as indicated in Section 2. Solve this LP iteratively for increasing values of the level of discretization n until the optimal objective value is suciently close to 0. Record the nal optimal solution as the optimal chessboard distribution qn . 2. Generate i.i.d. replicates of X by rst generating the cell containing U , then generating the individual components of U within the cell independently, and nally transforming U into X as outlined above. The computational demands of this method of random vector generation may be competitive with that of the NORTA procedure, depending on the discretization level n involved. To see why, rst

34

note that perhaps the most taxing part of this method (apart from the initialization - see below) is to determine the cell containing U . This involves generating the cell from the discrete probability distribution speci ed by qn . An optimal extreme point solution to the augmented LP will yield a solution qn with at most nd + d(d ? 1)=2 strictly positive qn (i1 ; : : : ; id) terms, because this is the number of constraints in the LP. Notice that this bound will typically be far smaller than nd , the number of cells in the nth LP. Therefore, if n is not too large, we can use the Alias method (see Law and Kelton 2000, p. 474) for example, to generate the cell containing U in constant time. The storage requirements for this method are directly related to the number of \bins" needed, which is at most nd + d(d ? 1)=2. It is the initial setup for nding the optimal distribution qn that could potentially be computationally expensive, because this involves solving a succession of LPs in an attempt to match the covariance matrix X . Nevertheless, the approach is certainly feasible. It is worth noting that a similar method to that described above could also be used in conjunction with the LPs of Section 3 to generate random vectors with nonuniform marginals and speci ed Pearson product-moment correlation matrix. However, the same cautionary remarks apply to the computational requirements of the initialization phase.

Acknowledgments We would like to thank Marina Epelman for a discussion that was helpful in proving Lemma 8. This work was partially supported by NSF grant DMI-9984717.

References Alfakih, A. and H. Wolkowicz. 2000. Matrix completion problems. In Handbook of Semide nite Programming: Theory, Algorithms and Applications. H. Wolkowicz, R. Saigal, L. Vandenberghe, eds, Kluwer, Boston, 533{545. Billingsley, P. 1986. Probability and Measure. Wiley, New York. Cario, M. C. and B. L. Nelson. 1996. Autoregressive to anything: time-series input processes for simulation. Operations Research Letters 19:51{58. Cario, M. C. and B. L. Nelson. 1997. Modeling and generating random vectors with arbitrary marginal distributions and correlation matrix. Technical Report, Department of Industrial Engineering and

35

Management Sciences, Northwestern University, Evanston, Illinois. Cario, M. C. and B. L. Nelson. 1997b. Numerical methods for tting and simulating autoregressive-toanything processes. INFORMS Journal on Computing 10:72{81. Chen, H. 2000. Initialization for NORTA: generation of random vectors with speci ed marginals and correlations. Working paper, Department of Industrial Engineering, Chung Yuan Christian University, Chung Li, Taiwan. Clemen, R. T., and T. Reilly. 1999. Correlations and copulas for decision and risk analysis. Management Science. 45:208{224. NEOS Optimization Server. 2000. Available at http://www-neos.mcs.anl.gov/neos/ [accessed August 2, 2000]. Dantzig, G. B. 1991. Converting a converging algorithm into a polynomially bounded algorithm. Technical report 91-5, Systems Optimization Laboratory, Dept of Operations Research, Stanford University, Stanford, California. Devroye, L. 1986. Non-Uniform Random Variate Generation. Springer-Verlag, New York. Henderson, S. G. 2000. Personal Website. Department of Industrial and Operations Engineering, University of Michigan, Michigan. http://www-personal.umich.edu/ shaneioe [accessed August 2, 2000]. Hill, R. R., and C. H. Reilly. 1994. Composition for multivariate random vectors. In Proceedings of the 1994 Winter Simulation Conference, J. D. Tew, S. Manivannan, D. A. Sadowsky, A. F. Seila, eds. IEEE, Piscataway New Jersey, 332 { 339. Hill, R. R., and C. H. Reilly. 2000. The eects of coecient correlation structure in two-dimensional knapsack problems on solution procedure performance. Management Science, 46: 302{317. Iman, R. and W. Conover. 1982. A distribution-free approach to inducing rank correlation among input variables, Communications in Statistics: Simulation and Computation, 11: 311-334. Johnson, C. R. 1990. Matrix completion problems: a survey. Proceedings of Symposia in Applied Mathematics 40:171{198. Johnson, M. E. 1987. Multivariate Statistical Simulation. Wiley, New York. Kruskal, W. 1958. Ordinal measures of associaton. J. Amer. Statist. Assoc. 53:814{861. Law, A. M. and W. D. Kelton. 2000. Simulation Modeling and Analysis. McGraw Hill, Boston. Lewis, P. A. W., E. McKenzie, and D. K. Hugus. 1989. Gamma processes. Communications in Statistics:

36

Stochastic Models 5:1{30.

Li, S. T., and J. L. Hammond. 1975. Generation of pseudorandom numbers with speci ed univariate distributions and correlation coecients. IEEE Transactions on Systems, Man, and Cybernetics. 5:557{561. Lurie, P. M., and M. S. Goldberg. 1998. An approximate method for sampling correlated random variables from partially-speci ed distributions. Management Science. 44:203{218. Mardia, K. V. 1970. A translation family of bivariate distributions and Frechet's bounds. Sankhya. A32:119{122. Melamed, B., J. R. Hill, and D. Goldsman. 1992. The TES methodology: modeling empirical stationary time series. In Proceedings of the 1992 Winter Simulation Conference IEEE, Piscataway, New Jersey, 135{144. Nelsen, R. B. 1999. An Introduction to Copulas. Lecture Notes in Statistics, 139. Springer-Verlag, New York. van der Geest, P. A. G. 1998. An algorithm to generate samples of multi-variate distributions with correlated marginals Computational Statistics and Data Analysis 27: 271-289. Whitt,W. 1976. Bivariate distributions with given marginals. The Annals of Statistics. 4:1280{1289. Wolkowicz, H., R. Saigal, L. Vandenberghe, eds. 2000. Handbook of Semide nite Programming: Theory, Algorithms and Applications. Kluwer, Boston.

37