# How to evaluate multivariate normal integral with conditional upper bounds

Suppose I have independent normally distributed random variables: $$x_i \sim N(0,1)$$. In my actual application, $$i=1,\ldots,30$$, but for my example here I'll use $$i=1,2,3$$.

I want to evaluate (either analytically with math or approximately with a simulation technique) the probability that they occur in order, after being "adjusted" by some real-valued, non-random scalars that I will call $$a, b, c$$.

For example, what is: $$\ \Pr(x_3 + c < x_2 + b < x_1 + a)$$?

Here's what I've got so far, where $$I(\cdot)$$ is an indicator function, and $$\phi$$ and $$\Phi$$ are the standard normal pdf and cdf:

\begin{align*} &\Pr(x_3 + c < x_2 + b < x_1 + a) \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} I(x_3 + c < x_2 + b) \times I(x_2 + b < x_1 + a) \phi(x_3)\phi(x_2)\phi(x_1) \, dx_3 \, dx_2 \, dx_1 \\ &= \int_{-\infty}^{\infty} \phi(x_1) \left[ \int_{-\infty}^{\infty} I(x_2 + b < x_1 + a) \phi(x_2) \left( \int_{-\infty}^{\infty} I(x_3 + c < x_2 + b) \phi(X_3) \, dx_3 \right) \, dx_2 \right] \, dx_1 \\ &= \int_{-\infty}^{\infty} \phi(x_1) \int_{-\infty}^{x_1 + a - b} \phi(x_2) \int_{-\infty}^{x_2 + b - c} \phi(x_3) \, dx_3 \, dx_2 \, dx_1 \\ &= \int_{-\infty}^{\infty} \phi(x_1) \int_{-\infty}^{x_1 + a - b} \Phi(x_2 + b - c) \phi(x_2) \, dx_2 \, dx_1 \\ &= \ \ldots ? \end{align*}

A solution to this, or even just a pointer to a simulation method that I can learn about (i.e., does GHK work here?) would be very helpful - thank you!

OK - maybe I have an answer in the spirit of GHK simulation.

It would still be helpful if someone could comment on the correctness of this answer...

\begin{align} &\Pr(x_3+c < x_2+b < x_1+a) \\ &= \int_{-\infty}^\infty \Pr(x_3+c < x_2+b < x_1+a \big| x_1) \phi(x_1)\,dx_1 \\ &= \int_{-\infty}^\infty \Pr(x_3+c < x_2+b \big| x_2+b < x_1+a)\Pr(x_2+b < x_1+a \big| x_1) \phi(x_1)\,dx_1 \\ \end{align}

Approximate this integral as folows:

1. repeat steps 2-5 R times:
2. draw $$x_1^r$$ from $$\phi(x_1)$$
3. calculate $$\Phi(x_1^r+a-b)$$
4. draw $$x_2^r$$ from $$\bar{\phi}(x_2)$$ where $$\bar{\phi}$$ is a standard normal truncated above at $$x_1^r+a-b$$
5. calculate $$\Phi(x_2^r+b-c)$$
6. calculate $$\hat{\Pr}(x_3+c < x_2+b < x_1+a) \approx \frac{1}{R}\sum_{r=1}^R \Phi(x_1^r+a-b)\Phi(x_2^r+b-c)$$