Normal orthant probability in six variables

Question

Let $\mathbf{X} \sim N(\mathbf{0}, \mathbf{\Sigma})$ be a $6$-dimensional Gaussian vector with covariance matrix of the form $$\mathbf{\Sigma} = \begin{pmatrix} 1 & c \\ c & 1 \end{pmatrix} \otimes \begin{pmatrix} 1 & 1/2 & 1/2 \\ 1/2 & 1 & 1/2 \\ 1/2 & 1/2 & 1 \end{pmatrix} = \begin{pmatrix} 1 & 1/2 & 1/2 & c & c/2 & c/2 \\ 1/2 & 1 & 1/2 & c/2 & c & c/2 \\ 1/2 & 1/2 & 1 & c/2 & c/2 & c \\ c & c/2 & c/2 & 1 & 1/2 & 1/2 \\ c/2 & c & c/2 & 1/2 & 1 & 1/2 \\ c/2 & c/2 & c & 1/2 & 1/2 & 1 \end{pmatrix}$$ Here $c$ is some constant between $0$ and $1$.

I am interested in $\Pr(\mathbf{X} \geq \mathbf{0})$, i.e. the probability that all of the six coordinates are non-negative/positive. This is also known as the orthant probability for $\mathbf{X}$, and explicit formulas for orthant probabilities for arbitrary covariance matrices $\mathbf{\Sigma}$ are known in $2, 3, 4$ dimensions. The analysis seems very tedious and nasty though, and I am not sure if the same techniques generalize easily to $6$ dimensions.

For my example of $6$ dimensions, my search has not yet returned any literature, attempting to solve this problem even in special cases with additional structure like above. I also tried computing this probability with Mathematica, but it cannot solve it analytically, and even numerically it seems to have a hard time to return exact results.

My question is: is there any way to find the orthant probability for such a structured $6$-dimensional matrix? Is there any literature I am missing? Or can someone solve this analytically?

---

Edit: The paper https://ieeexplore.ieee.org/document/1054159 describes a derivation of the orthant probability for four variables, where \begin{align} \mathbf{\Sigma} = \begin{pmatrix} 1 & c \\ c & 1 \end{pmatrix} \otimes \begin{pmatrix} 1 & d \\ d & 1 \end{pmatrix}. \end{align} Perhaps it is possible to use similar techniques for this case as well (using a similar path of integration approach), although the reduction at the bottom of page 389 would not go from a $4$-dimensional orthant probability to a $2$-dimensional one, which can readily be evaluated, but from $6$ dimensions to $4$ dimensions, and adding another integral over these orthant probabilities.

Still, this reference might be useful -- it might inspire similar techniques/approaches to this problem.

Answer 1

This is not a complete answer, but maybe this will help you to get a formula in terms of the distribution function of the multivariate normal distribution.

Lets denote by $X = (X_1,\ldots,X_6) \sim N(0,\Sigma)$ the original random vector. Set $I = \{1,2,3,4,5,6\}$. For $J \subset I$ I will denote by $F_J$ the distribution function of the random variables $X_{J_1},\ldots,X_{J_{m}}$, where $m = |J|$. For example taking $J = \{1,3,5\}$ we have $m = 3$ and $F_J(x) = \Pr[X_1 \leq x_1, X_3 \leq x_2, X_5 \leq x_3]$ denotes the distribution function of $(X_1,X_3,X_5)$ at $x=(x_1,x_2,x_3)$ which is again jointly normal. The joint distribution function of $X$ is then $F_I$ and for $J=\{2,6\}$ we have $F_J(x) = \Pr[X_2 \leq x_1, X_6 \leq x_2]$ and so forth.

Your orthant probability can be thought of as computing the measure of a box under the joint probability measure and is linked to the inclusion / exclusion principle. To express this in terms of the joint distribution function I will follow Chapter 10 in Nelsen (2006): An introduction to copulas.

For two vectors $a < b$ (meaning $a_i < b_i$ for all $i$) the corresponding box is simply $B = \times_{i=1}^6 [a_i,b_i]$, and $c=(c_1,\ldots,c_6)$ will denote a vertex (or corner) of $B$, that is: $c_k$ is either equal to $a_k$ or $b_k$. By equation (2.10.1) in Nelsen (2006) the measure of $B$ under $V_{F_I}$ (the measure induced by $F_I$) is now given as \begin{align} V_{F_I}(B) = \sum \text{sgn}(c)F_I(c), \end{align} where the sum is taken over all vertices $c$ of $B$, and $\text{sgn}(c)$ is given as \begin{align} \text{sgn}(c) = 1, \text{if } c_k = a_k \text{ for an even number of $k$'s and $-1$ otherwise.} \end{align} Here the sum will have $2^6 = 64$ terms which seems unmanageable, but in your case two simplifications are possible.

First, your orthant is a box with $a = (0,\ldots,0)$ and $b=(\infty,\ldots,\infty)$. If a corner $c$ has an infinity as a component, the joint distribution $F_I(c)$ will actually reduce to a lower dimensional one. For example in case of $c=(0,\infty,0,\infty,0,\infty)$ we have $F_I(c) = F_J(0,0,0)$ where $J=\{1,3,5\}$.

Second, with the specific covariance structure that you have in mind these lower dimensional distributions might be the same for some index sets, meaning that it is possible that $F_{J_1} = F_{J_2}$ for some $J_1$ and $J_2$ with $|J_1|=|J_2|$.

What you would need to check now is which lower dimensional distributions are the same (this depends on $\Sigma$) and which sign they get. Then you can collect equal terms in the sum and rewrite the sums in terms of the dimension of the lower dimensional distributions going from one to six.

In your case you will have for each dimension multiple terms with different coefficients, but if your random vector had the same covariance between all components you would end up with six terms only (and the coefficients would be related to the binomial coefficients).

I hope this helps or gives maybe a different direction to tackle your problem. Good luck!