Given $(X,Y)$ follows a bivariate normal distribution $\mathcal{N}\left(\begin{pmatrix} 0\\ 0 \end{pmatrix};\begin{pmatrix} 1 & r\\ r &1 \end{pmatrix} \right)$ with $-1<r<1$
Using Mathematica, I found this interesting identity formula
$$\Bbb E(X\mathbb{I}_{\{X<a\}}\mathbb{I}_{\{Y<b\}}) = -\varphi(a)\cdot\Phi\left(\frac{b-ar}{\sqrt{1-r^2}}\right) -r\cdot \varphi(b)\cdot\Phi\left(\frac{a-br}{\sqrt{1-r^2}} \right) \tag{1}$$
with
- $a,b\in \Bbb R$,
- $\varphi(x)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x^2}{2}\right)$ the standard normal density function
- $\Phi(x)=\int_{-\infty}^x\varphi(t)dt$ the standard normal probability function
- $E(X\mathbb{I}_{\{X<a\}}\mathbb{I}_{\{Y<b\}})$ the expectation of $X\mathbb{I}_{\{X<a\}}\mathbb{I}_{\{Y<b\}}$ with $\Bbb I_{\{...\}}$ is the identity function
How can we prove this identity $(1)$?
I denote the correlation by $\rho$ and the standard normal pdf by $\phi(\cdot)$.
The conditional distribution of $Y$ given $X=x$ is known to be $N(\rho x,1-\rho^2)$.
Then,
\begin{align} \mathbb E\left[X \mathbf1_{X<a}\mathbf1_{Y<b}\right]&=\mathbb E \,\mathbb E\left[X \mathbf1_{X<a}\mathbf1_{Y<b}\mid X\right] \\&=\mathbb E \left[X \mathbf1_{X<a} \,\mathbb P(Y<b\mid X)\right] \\&=\mathbb E \left[X \mathbf1_{X<a} \,\Phi\left(\frac{b-\rho X}{\sqrt{1-\rho^2}}\right)\right] \\&=\int_{-\infty}^a x\,\phi(x)\,\Phi\left(\frac{b-\rho x}{\sqrt{1-\rho^2}}\right) \,\mathrm dx \end{align}
Using $\phi'(x)=-x \phi(x)$, integrate by parts to get
\begin{align} \int_{-\infty}^a x\,\Phi\left(\frac{b-\rho x}{\sqrt{1-\rho^2}}\right) \phi(x)\,\mathrm dx &=\lim_{A\to -\infty}\left(-\Phi\left(\frac{b-\rho x}{\sqrt{1-\rho^2}}\right)\phi(x)\Big|_{A}^a\right)-\rho I \\&=-\Phi\left(\frac{b-\rho a}{\sqrt{1-\rho^2}}\right)\phi(a) - \rho I\,, \end{align}
where
$$I=\int_{-\infty}^a \frac1{\sqrt{1-\rho^2}}\phi\left(\frac{b-\rho x}{\sqrt{1-\rho^2}}\right) \phi(x)\,\mathrm{d}x$$
To find $I$, write down the expressions of $\phi(\cdot)$ and complete squares to get a normal pdf within the integral:
\begin{align} I&=\frac1{2\pi\sqrt{1-\rho^2}}\int_{-\infty}^a \exp\left[-\frac12\left\{\left(\frac{b-\rho x}{\sqrt{1-\rho^2}}\right)^2+x^2\right\}\right]\,\mathrm{d}x \\&=\frac1{2\pi\sqrt{1-\rho^2}}\int_{-\infty}^a \exp\left[-\frac1{2(1-\rho^2)}\left(x^2+b^2-2b\rho x\right)\right]\,\mathrm{d}x \\&=\frac{e^{-b^2/2}}{\sqrt{2\pi}} \int_{-\infty}^a \frac1{\sqrt{2\pi}\sqrt{1-\rho^2}}\exp\left[-\frac1{2(1-\rho^2)}(x-b\rho)^2\right]\,\mathrm{d}x \\&=\phi(b)\, \Phi\left(\frac{a-b\rho}{\sqrt{1-\rho^2}}\right) \end{align}
Thus,
$$\mathbb E\left[X \mathbf1_{X<a}\mathbf1_{Y<b}\right]=-\phi(a)\Phi\left(\frac{b-\rho a}{\sqrt{1-\rho^2}}\right) -\rho\, \phi(b)\, \Phi\left(\frac{a-b\rho}{\sqrt{1-\rho^2}}\right)$$