The Standard TVN Problem

Let the standard trivariate normal distribution function be defined by

$$\Phi({\bf b}, R) = \frac{1}{{(2\pi)}^{3/2}{\vert R\vert}^{1... ...ty}^{b_3} e^{-\frac{{\bf x}^TR^{-1}{\bf x}}{2}} dx_3dx_2dx_1,$$ (8)

where ${\bf b}= (b_1, b_2, b_3)$ and $R = (\rho_{ij})$ is a correlation matrix. More general TVN problems, where the covariance matrices do not have unit diagonal entries, can always be rescaled to be problems where the covariance matrix is a correlation matrix.

Owen (1956) briefly discussed the evaluation of the trivariate integral and gave the formula for the standard trivariate normal integral in terms of the bivariate normal integral as:

$$\Phi({\bf b}, R) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{b_1} e^{-x^2/2} F(x)dx,$$ (9)

where

$$F(x) = \Phi\bigg(\Big(\frac{b_2-\rho_{21}x}{(1-\rho^2_{21})^... ...(1-\rho^2_{21})^\frac{1}{2}(1-\rho^2_{31})^\frac{1}{2}}\bigg),$$

and $\Phi(h,k,\rho)$ is the standard bivariate normal distribution function.

Two methods were tested by the present author for calculating (9) using numerical integration. The first method that was tested transforms (9) using $x = \Phi^{-1}(t)$, so that

$$\Phi({\bf b}, R) = \int_{0}^{\Phi(b_1)} F(\Phi^{-1}(t)) dt,$$ (10)

and then a Gauss-Legendre numerical integration rule is used for the interval $[0, \Phi(b_1)]$.

The second method that was tested, which was also considered by Drezner (1992) as a method for general multivariate normal distribution computations, uses modified Gauss-Hermite rules (see Steen, Byrne and Gelhard, 1969). These rules are suitable for integrals in the form $\frac{1}{\sqrt{2\pi}}\int_0^\infty e^{-x^2/2} f(x) dx$, so (9) needs to be transformed to a $[0, \infty]$ integral. Let $y = x - b_1$, and then

$$\Phi({\bf b}, R) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^0 e^... ...^2/2}}{\sqrt{2\pi}}\int_{-\infty}^0 e^{-y^2/2-yb_1}F(y+b_1)dy.$$

Now let $z = -y$, so that

$$\Phi({\bf b}, R) = \frac{e^{-b_1^2/2}}{\sqrt{2\pi}} \int_0^{\infty}e^{-z^2/2+zb_1}F(b_1-z)dz.$$ (11)

It was found that both of these methods were often more accurate for a particular integration rule if the integration limits were arranged so that the shortest integration interval was associated with the outermost integration.