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?
Best Answer
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.