I am trying to compute a multivariate normal cdf where all but the last bounds of the integrals are symmetric:
$$F(a, \sigma, m ) = \frac1{\sqrt{(2\pi)^m|\Sigma|}}\int_{-a}^a{\int_{-a}^a...{\int_a^\infty}}e^{-\frac12X^t\Sigma^{-1}X}dX$$
Where X is a m-dimensional multivariate normal vector with mean $\mu_i=0$ and covariance matrix $\Sigma$ created by: a correlation matrix $\rho_{[m \times m]}$ : $\rho_{ij} = \sqrt{\frac ij}$ for i < j and vector of variances $v\in\Re^m$ : $v_i=\sigma i$, $\sigma \in \Re$.
Given the symmetry of the bounds and the special type of the covariance matrix, is there a way to simplify the expression for F() or to get a reasonably good numerical approximation in order to compute F() for values of $m \approx 300$? I have been also trying to find some recursive dependency on m: $F(a, \sigma, m)$ as a function of $F(a, \sigma, m-1)$ and/or a scaling on $\sigma$: $F(a, \sigma, m)$ as a function of $F(a, 1, m)$
Thank you for your help and ideas