General Correlation Structure

We begin with the standard definition of the multivariate normal distribution when the variance-covariance matrix $\Sigma$ is positive definite. Define

$$F({\bf a},{\bf b},\Sigma)= \Pr\biggl[ \bigcap_{j=1}^m \biggl(... ...^{- \frac{1}{2} {\bf x}' \Sigma^{-1} {\bf x}} dx_m\cdots dx_1.$$

A key step in the development of the numerical methods described by Genz (1992) was to transform F using ${\bf x}= C{\bf y}$, where C is the Cholesky factor of $\Sigma$ ( $\Sigma = CC'$). After this transformation

$$F({\bf a},{\bf b},\Sigma)= (2\pi)^{-\frac{m}{2}} \int_{{\bf a... ...\bf y}\leq {\bf b}}e^{- \frac{1}{2} {\bf y}'{\bf y}} d{\bf y}.$$

We wish to generalise this approach for cases where $\Sigma$ is positive semi-definite, and show that a similar definition of F can also be used in these cases. Let $\Sigma = UDU'$ be a singular value decomposition for $\Sigma$, where U is an $m\times m$ orthogonal matrix and D is an $m\times m$ diagonal matrix with nonnegative diagonal entries $d_1, d_2, \ldots, d_m$. If k is the rank of $\Sigma$, and k<m, then $d_{k+1}=d_{k+2}=\ldots=d_m=0$. Now define $\Sigma(\epsilon) = \Sigma+\epsilon I$, for $\epsilon > 0$. Then

$$\Sigma(\epsilon)=UDU' +\epsilon I = U(D +\epsilon I)U'= UD(\epsilon)U',$$

where $D(\epsilon) = D+\epsilon I$, so $\Sigma(\epsilon)$ is positive definite, and $F({\bf a},{\bf b},\Sigma(\epsilon))$ is properly defined. Next let $V(\epsilon)=UD(\epsilon)^\frac{1}{2}$, so $\Sigma(\epsilon)= V(\epsilon)V'(\epsilon)$, and let ${\bf x}= V(\epsilon){\bf v}$. Then

$$F({\bf a},{\bf b},\Sigma(\epsilon))= (2\pi)^{-\frac{m}{2}} \i... ...\bf v}\leq {\bf b}}e^{- \frac{1}{2} {\bf v}'{\bf v}} d{\bf v}.$$

Taking the limit as $\epsilon$ approaches zero, we have

$$F({\bf a},{\bf b},\Sigma)= (2\pi)^{-\frac{m}{2}} \int_{{\bf a... ...\bf v}\leq {\bf b}}e^{- \frac{1}{2} {\bf v}'{\bf v}} d{\bf v},$$

where V = V(0). V can be written as $V = \hat{U}D^\frac{1}{2}$, where $\hat{U}$ is an $m\times m$ matrix with its first k columns as the first k columns of U and remaining m-k columns as columns of zeros, because $d_{k+1}=d_{k+2}=\ldots=d_m=0$. For a final step, we determine an $m\times m$ orthogonal matrix Qso that C=VQ' is lower triangular. The Q required for this step is an $m\times m$ identity matrix except for its principal $k \times k$ submatrix which is the $k \times k$ orthogonal matrix required to make the principal $k \times k$ submatrix of V lower triangular (see the book by Golub and Van Loan, 1996, for a description of how such a Q can be determined). Then, if ${\bf v}= Q'{\bf y}$, $F({\bf a},{\bf b},\Sigma)= (2\pi)^{-\frac{m}{2}} \int_{{\bf a}\leq C{\bf y}\leq {\bf b}}e^{- \frac{1}{2} {\bf y}'{\bf y}} d{\bf y}$, where C is lower triangular and c_ij = 0 if j>k. The integration region defined by ${\bf a}\leq C{\bf y}\leq {\bf b}$ places no constraint on the variables $y_{k+1}, y_{k+2}, \ldots, y_m$, and all possible values between $-\infty$ and $\infty$ for these variables are consistent with the definition of the integration region, so the innermost m-k integrals are all equal to one. Therefore

$$F({\bf a},{\bf b},\Sigma)= (2\pi)^{-\frac{k}{2}} \int_{{\bf a... ...\bf y}\leq {\bf b}}e^{- \frac{1}{2} {\bf y}'{\bf y}} d{\bf y},$$

where C is now the $m\times k$ matrix obtained by removing the original m-kzero columns from the original C and ${\bf y}= (y_1, y_2, \ldots, y_k)'$. This matrix C can be computed directly using the generalised Cholesky decomposition algorithm described by Healy (1968).

Following the procedure used by Genz (1992), the next step is to use the lower triangular structure of C to rewrite each of the inequalities ${\bf a}\leq C{\bf y}\leq {\bf b}$ and explicitly define the integration limits for F to give F in the form

$$F({\bf a},{\bf b},\Sigma) = (2\pi)^{-\frac{k}{2}} \int_{a'_1}... ...}^{b'_k(y_1, \ldots, y_{k-1})} e^{-\frac{y_k^2}{2}} d{\bf y},$$

where $a'_i(y_1,\ldots,y_{i-1})=(a_i-\sum_{j=1}^{i-1}c_{i,j}y_j)/c_{i,i}$and $b'_i(y_1, \dots, y_{i-1})=(b_i-\sum_{j=1}^{i-1}c_{i,j}y_j)/c_{i,i}$, for $i = 1, 2, \ldots, k$. But this definition of F is not complete if k<m, because we must take into account the m-k additional constraints $a_i \leq \sum_{j=1}^{i}c_{i,j}y_j \leq b_i$, for $i = k+1, k+2, \ldots m$, that the k integration variables must satisfy. There are various cases to consider. In order to introduce the general cases, we first consider the case where $c_{ik} \neq 0$ for all i>k. In this case we only need to place additional constraints on y_k, and modify the definitions of a'_k and b'_k. We first rewrite the last m-k+1 constraints in the form $a_i-\sum_{j=1}^{k-1}c_{i,j}y_j \leq c_{i,k}y_k \leq b_i-\sum_{j=1}^{k-1}c_{i,j}y_j$, for $i = k, k+1, \ldots m$. Then we divide these constraints into two groups, the first group consisting of the constraints where c_i,k > 0 and the second group consisting of the constraints where c_i,k < 0. For both groups, we divide by c_i,k to produce explicit constraints on y_k, but for the second group we must change the order of the inequalities. The revised limits for y_k can now be written as

$$\bar{a}'_k(y_1,\ldots,y_{k-1}) = \max\Big( \max_{c_{i,k}>0}\B... ...\Big(\frac{b_i-\sum_{j=1}^{k-1}c_{i,j}y_j}{c_{i,k}}\Big)\Big),$$

and

$$\begin{eqnarray*}\bar{b}'_k(y_1,\ldots,y_{k-1})= \max\Big(&& \bar{a}'_k(y_1,\ldo... ...rac{a_i-\sum_{j=1}^{k-1}c_{i,j}y_j)}{c_{i,k}}\Big)\Big) \ \Big). \end{eqnarray*}$$

In the cases where c_i,k=0 for some i's, with i>k, the associated constraints for these i's do not affect y_k, so the constraints for some of other variables need to be adjusted. For these general cases, we first reorder all of the constraints into groups of sizes $l_1, l_2, \ldots, l_k$, with $l_1+l_2+\cdots +l_k = m$, so that if we denote the reordered constraints by $\bar{{\bf a}} \leq \bar{C}{\bf y}\leq \bar{{\bf b}}$, then the first l₁ rows of $\bar{C}$ have the form $(*,0,\ldots,0)$, the next l₂ rows of $\bar{C}$ have the form $(?,*,0,\ldots,0)$, and so on with the last l_k rows of $\bar{C}$ in the form $(?,\ldots,?,*)$, where * denotes a nonzero and ? denotes zero or nonzero. Then, for each i, $i = 1, \ldots, k$, the set of l_i constraints are rewritten and merged to produce a single constraint on y_i using the procedure that we described for y_k. When this process is complete, the integral for F can be written as

$$F({\bf a},{\bf b},\Sigma) = (2\pi)^{-\frac{k}{2}} \int_{\bar{... ...r{b}'_k(y_1, \ldots, y_{k-1})} e^{-\frac{y_k^2}{2}} d{\bf y}.$$

In order to put F into a form that is easy to use with standard numerical integration methods, two more transformations are required. First, let $y_i = \Phi^{-1}(z_i)$, for i = 1,2,..kso that $\phi(y_i)dy_i = dz_i$. Then

$$F({\bf a},{\bf b},\Sigma) = \int_{d_1}^{e_1} \int_{d_2(z_1)}^... ...k(z_1, \ldots, z_{k-1})}^{e_k(z_1, \ldots, z_{k-1})} d{\bf z},$$

with

$$d_i(z_1,\ldots,z_{i-1}) = \Phi(\bar{a}'_i(\Phi^{-1}(z_1),\ldots,\Phi^{-1}(z_{i-1}))),$$

and

$$e_i(z_1,\ldots,z_{i-1}) = \Phi(\bar{b}'_i(\Phi^{-1}(z_1),\ldots,\Phi^{-1}(z_{i-1}))).$$

Finally, let z_i=d_i+(e_i-d_i)u_i, for $i = 1, \ldots, k$, so dz_i=(e_i-d_i)du_i. Then

$$F({\bf a},{\bf b},\Sigma) = (e_1-d_1) \int_0^1 (e_2(u_1)-d_2(... ...ldots, u_{k-1})-d_{k}(u_1, \ldots, u_{k-1}))\int_0^1 d{\bf u}.$$

The innermost integral is one, so the numerical problem involves the estimation of a (k-1) dimensional integral.