An interesting identity about the expectation $\Bbb E(X\mathbb{I}_{\{XGiven $(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)$?
 A: 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)$$
