Numerically it appears that the following function is increasing (or at least non-decreasing) in $σ$, $$ f(x_1,x_2;σ)=∫_{-∞}^{∞}∫_{-∞}^{x_2-σx_1+σθ_1}[1-Φ(x_1-θ_1)]φ(θ_1)φ(θ_2)dθ_2dθ_1 $$ where $\phi$ and $\Phi$ are the density and distribution functions of the standard normal distribution. I would like to formally prove this. All variables are real and $\sigma$ is positive: $x_1, x_2, σ \in \mathbb{R}$ and $σ > 0$.
I have tried taking the derivative, $$ \frac{\partial f}{\partial \sigma} = ∫_{-∞}^{∞}[1-Φ(x_1-θ_1)]φ(θ_1)[(θ_1-x_1)φ(x_2-σx_1+σθ_1)]dθ_1 $$
but am struggling to show that this is positive. I am also familiar with the table of transformations listed in Owen (1980) and have tried writing my equation in terms of his equations 10,010.8 and 20,010.3, and the differentiating, but haven't made much progress that way either.
Substitute $y=x_1-\theta_1$, then the last integral becomes \begin{equation}\label{*}\tag{*}-\int_{-\infty}^\infty (1-\Phi(y))\phi(x_1-y) y \phi(x_2-\sigma y)\,\mathrm dy.\end{equation}
If we take for instance $\sigma=1$ and $x_1=x_2=1$, then this is (numerically*)
$$\approx-0.0255037<0.$$
So sadly $f$ is not always increasing with $\sigma$.
* Intuitively this also makes sense: When $x_1=x_2=1$, then the $\phi$'s are concentrated around $y=1$, so the $y=1$ part is contributing more to the integral than the $y<0$ part. (Note that $1-\Phi(y)$ is a term acting in the other direction, which is also why \eqref{*} is positive for $x_1=x_2=0, \sigma=1$.)