Proving increasing function defined as bivariate normal Suppose $c>0,\sigma>0$ and $\tau>0$ are fixed real constants. Then I'd like to prove that the function $g_c:(-1,1)\mapsto\mathbb{R}$
 defined by
 \begin{equation}
 g_c(\rho)=\int_{-\infty}^\infty\int_{-\infty}^\infty\frac{ \{(-c)\vee x \wedge c\} \{(-c)\vee y \wedge c\}}{\sqrt{2\pi(1-\rho^2)}\sigma\tau}
 e^{-\frac{\sigma^2x-2\rho\sigma\tau xy+\tau^2y^2}{2(1-\rho^2)\sigma^2\tau^2}}dxdy
 \end{equation}
 is strictly increasing.
I have tried to prove by differentiating $g_c$ w.r.t. $\rho$. But, this doesn't help because the result is very ugly. Can anyone help me? Thanks in advance
PS: The function $g_c$ can be rewritten as $g_c(\rho)=E( (-c)\vee [(X_{\sigma,\rho}Y_{\tau,\rho})\wedge c])$ for some random variable $(X_{\sigma,\rho},Y_{\tau,\rho})^T\sim N_2(0,\Sigma)$ where
 $$
 \Sigma=
\begin{pmatrix}
\sigma^2 & \rho\sigma\tau \\
\rho\sigma\tau & \tau^2
\end{pmatrix}.
$$
But, I also don't know how to see my problem using this fact.
 A: EDIT: Ups I did a change of variables wrong.
EDIT 2: Ups also forgot to scale the pdf correctly
First define $V = XY$ and find the density of that. It will be given by 
$$\int_{-\infty}^\infty f_{X,Y}(x,v/x)\frac{1}{|x|}\mathrm{d}x$$
You can find this using mathematica (probably in Abr. & Steg. also if you're a purist)
$$f_V(v) = 2 \frac{1}{\sqrt{2 \pi  \left(1-\rho ^2\right)} \sigma  \tau } e^{\frac{\rho  v}{\sigma  \tau -\rho ^2 \sigma  \tau }} K_0\left(\frac{\left| v\right| }{\sigma  \tau -\rho ^2 \sigma  \tau }\right)$$
where $K_0$ is a modified Bessel function of the second kind.
Then you're interested in $g(\rho) = \mathbb{E}[\max(-c,\min(c,V))]$.
$$\frac{\sqrt{ \pi  \left(1-\rho ^2\right)}\sigma \tau}{\sqrt{2}}g(\rho) = -c \int_{-\infty}^{-c} e^{\rho v/((1-\rho^2)\sigma\tau)}K_0(-v/(1-\rho^2))\mathrm{d}v\\
+ \int_{-c}^0 v e^{\rho v/((1-\rho^2)\sigma\tau)}K_0(-v/((1-\rho^2)\sigma\tau))\mathrm{d}v\\
+ \int_0^c v e^{\rho v/((1-\rho^2)\sigma\tau)}K_0(v/((1-\rho^2)\sigma\tau))\mathrm{d}v\\
+ c\int_c^\infty e^{\rho v/((1-\rho^2)\sigma\tau)}K_0(v/((1-\rho^2)\sigma\tau))\mathrm{d}v$$
change variables to $w = -v/((1-\rho^2)\sigma\tau)$ in the first two integrals and $w = v/((1-\rho^2)\sigma\tau)$ in the second two.
$$\frac{\sqrt{\pi  \left(1-\rho ^2\right)}\sigma \tau}{\sqrt{2}}g(\rho) = -c ((1-\rho^2)\sigma\tau) \int_{c/((1-\rho^2)\sigma\tau)}^\infty e^{-\rho w}K_0(w)\mathrm{d}w\\
- ((1-\rho^2)\sigma\tau)^2 \int_0^{c/((1-\rho^2)\sigma\tau)} w e^{-\rho w}K_0(w)\mathrm{d}w\\
+ ((1-\rho^2)\sigma\tau)^2 \int_0^{c/((1-\rho^2)\sigma\tau)} w e^{\rho w}K_0(w)\mathrm{d}w\\
+ c ((1-\rho^2)\sigma\tau) \int_{c/(1-\rho^2)}^\infty e^{\rho w}K_0(w)\mathrm{d}w\\$$
$$ \frac{\sqrt{ \pi  \left(1-\rho ^2\right)}\sigma \tau}{2 \sqrt{2}}g(\rho)= c((1-\rho^2)\sigma\tau)\int_{c/((1-\rho^2)\sigma\tau)}^\infty \sinh(\rho w)K_0(w)\mathrm{d}w\\
+ ((1-\rho^2)\sigma\tau)^2 \int_0^{c/((1-\rho^2)\sigma\tau)} w \sinh(\rho w)K_0(w)\mathrm{d}w $$
EDIT: I think applying the Liebniz rules to the above might still work out like it did when I had a mistake in the change of variables. It's a little more messy because there's more product rules, but I think it will still work.
