Probability of three correlated gaussian variables Given three standard gaussian random variable $X,Y$ ($X$ and $Y$ are independent) and $Z$ with
$$X,Y = \mathcal{N}(0,1)$$
$$Z = Y-X$$
and  $a,b,c \in \Bbb R$

Is it possible to calculate analytically this probability $$\Bbb
  P(X<a,Y<b, Z<c)$$

(by analytically, I mean we can write this probability by using basic mathematical operations and standard normal probability functions $\Phi_n(\mathbf{z};\mathbf{0}_n, \mathbf{I}_n)$)
I compute the covariance matrix $\mathbf{\Sigma}$ of $(X,Y,Z)$, which is
$$\mathbf{\Sigma} = \begin{pmatrix}
   1 & 0 & -1\\
   0 & 1 & 1 \\
   -1 & 1 & 2 \\
   \end{pmatrix}$$
but this covariance matrix is unluckily singular.
Perhaps there is a clever way to decompose the region $\{X<a,Y<b,Z<c \}$ into several regions that can be analytically calculated?
 A: I'll share my progress until now.
For short, we will define $\Phi(x)=\Phi(x;0,1)$, that is, the cummulative distribuition function of the Standard Normal and we use $\phi(x)=\Phi'(x)$ for the pdf. Using Bayes' Theorem, the independence of $X$ and $Y$ and the fact that $Z=Y-X$ we have
\begin{align*}
\mathbb{P}(X<a,Y<b,Z<c) &= \mathbb{P}(X<a)\mathbb{P}(Y<b,Z<c|X<a) \\
&= \mathbb{P}(X<a)\mathbb{P}(Y<b,Y<c+X|X<a) \\
&= \mathbb{P}(X<a)\mathbb{P}(Y<b)\mathbb{P}(Y<c+X|X<a,Y<b).
\end{align*}
Let's calculate
$$\mathbb{P}(Y<c+X|X<a,Y<b)=\mathbb{P}(c+X\geq b)\mathbb{P}(Y<c+X|X<a,Y<b,c+X\geq b) + \mathbb{P}(c+X< b)\mathbb{P}(Y<c+X|X<a,Y<b,c+X< b).$$
See that
$$\mathbb{P}(Y<c+X|X<a,Y<b,c+X\geq b)\geq\mathbb{P}(Y<b|X<a,Y<b,c+X\geq b)=1,$$
i.e.,
$$\mathbb{P}(Y<c+X|X<a,Y<b,c+X\geq b)=1.$$
Moreover, define
$$L:=\mathbb{P}(Y<c+X|X<a,Y<b,c+X< b)=\mathbb{P}(Y<c+X|X<\min\{a,b-c\},Y<b).$$
Let's define $d:=\min\{a,b-c\}$. Then
\begin{align*}
L&=\mathbb{P}(Y<c+X|X<d,Y<b)\\
&= \frac{1}{2\pi}\int_{-\infty}^d \int_{-\infty}^b 1_{\{y-x<c\}}e^{-(x^2+y^2)/2}dy dx \\
&= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^d e^{-x^2/2} \int_{-\infty}^b 1_{\{y-x<c\}}\frac{1}{\sqrt{2\pi}}e^{-y^2/2}dy dx \\
&= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^d e^{-x^2/2} \mathbb{P}(Y<x+c) dx \\
&= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^d e^{-x^2/2} \Phi(x+c) dx \\
&= \int_{-\infty}^d \phi(x) \Phi(x+c) dx
\end{align*}
Finally, we have
\begin{align*}
\mathbb{P}(X<a,Y<b,Z<c) &= \mathbb{P}(X<a)\mathbb{P}(Y<b)\mathbb{P}(Y<c+X|X<a,Y<b) \\
&=\Phi(a)\Phi(b)(\mathbb{P}(c+X\geq b) + \mathbb{P}(c+X< b)\mathbb{P}(Y<c+X|X<a,Y<b,c+X< b)) \\
&=\Phi(a)\Phi(b)((1-\Phi(b-c)) +\Phi(b-c)L) \\
&=\Phi(a)\Phi(b)(1-\Phi(b-c)(L-1)) \\
\end{align*}
I don't quite know how to deal with this $L$ yet. I'll try to think about that more later. Maybe we can use some idea of this post.
