Integrating a special skew normal -- the CDF of a convolution of a normal with a truncated normal I am having a little trouble trying to compute an integral.  In short, I wish to solve the following:
$$F(x) = \int_{-\infty}^x \phi(au-b)\,\Phi(au+b)\,du $$
My intuition is that this might be impossible.  So, failing that, I am trying to solve a special case:
$$ F(0) = \int_{-\infty}^0 \phi(au-b)\,\Phi(au+b)\,du $$
Edit: Where $$\phi(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} $$ 
And: $$ \Phi(x) = \int_{-\infty}^x \phi(t) \, dt $$
My solution method (so far) has been trying to use integration by parts (several times) but this is resulting in a rather large mess.  I was curious if anyone had a (perhaps elegant) way of solving this problem?
I have spent considerable time using the following resources:


*

*Wikipedia's list of integrals of Gaussian functions.  These are very helpful.

*This math overflow answer solves for the integral over the entire domain.  I am looking for either the indefinite integral or the solution given the limits $(-\infty,0)$.

*Here is a paper that includes many identities that I have found helpful.

*I think this is a special case of (some sort of) Skew Normal distribution: another way of framing the question is to look for the CDF of the convolution of a normal pdf with a truncated normal pdf.  This paper by Azzalini might be of passing interest.
Addition: I have noticed the following identity from Owen:
$$\int_{-\infty}^{0}\phi(ax)\Phi(bx)dx = \frac{1}{2\pi a} \arctan\left( \frac{a}{b} \right) $$
 A: For every positive $a$, the change of variable $v=au$ yields $F(0)=G(b)/a$ with

$$
G(b)=\int_{-\infty}^0\varphi(v-b)\Phi(v+b)\mathrm dv.
$$

First, some easy remarks: $G(0)=\frac12\left.\Phi(v)^2\right|_{-\infty}^0=\frac18$, and, since $\Phi\leqslant1$, $G(b)\leqslant\Phi(-b)$, in particular $G(+\infty)=0$. 
Now, the computation of $G$. Note that $\Phi'=\varphi$ and $\varphi'(u)=-u\varphi(u)$, hence $G'(b)=H(b)+K(b)$ with
$$
H(b)=\int_{-\infty}^0(v-b)\varphi(v-b)\Phi(v+b)\mathrm dv,
$$
and
$$
K(b)=\int_{-\infty}^0\varphi(v-b)\varphi(v+b)\mathrm dv.
$$
Integrating by parts $H(b)$ yields
$$
H(b)=\left.-\varphi(v-b)\Phi(v+b)\right|_{-\infty}^0+\int_{-\infty}^0\varphi(v-b)\varphi(v+b)\mathrm dv,
$$
that is,
$$
H(b)=-\varphi(b)\Phi(b)+K(b).
$$
Now, $\varphi(v-b)\varphi(v+b)=\varphi(b\sqrt2)\varphi(v\sqrt2)$, hence
$$
\sqrt2K(b)=\varphi(b\sqrt2)\left.\Phi(v\sqrt2)\right|_{-\infty}^0=\frac1{2}\varphi(b\sqrt2).
$$
Thus,
$$
G'(b)=-\varphi(b)\Phi(b)+2K(b)=-\varphi(b)\Phi(b)+\frac1{\sqrt2}\varphi(b\sqrt2).
$$
Together with the limit values computed above, this shows that, for every positive $a$ and for every $b$,

$$
\int_{-\infty}^0\varphi(au-b)\Phi(au+b)\mathrm du=\frac{\Phi(b\sqrt2)-\Phi(b)^2}{2a}.
$$

The same technique seems to yield, not explicit formulas for $G_x(b)$, but some funny-looking duality relations between $G_x(b)$ and $G_b(x)$, where
$$
G_x(b)=\int_{-\infty}^x\varphi(v-b)\Phi(v+b)\mathrm dv,
$$
namely, one might have (but this should be checked)

$$
G_x(b)+G_b(x)=\Phi(b\sqrt2)\Phi(x\sqrt2).
$$

A: So I think I have a "solution" to this problem: there is no nice solution.  
Because of $b$ in the equation, this integral is equivalent to a CDF from a bivariate normal.
From Owen, there is result 10,010.1, which states:
$$ \int_{-\infty}^Y \phi(x)\Phi(a+bx)\,dx = \text{BvN}\left[ \frac{a}{\sqrt{a+b^2}},Y; \rho = \frac{-b}{\sqrt{1+b^2}} \right] $$
For future reference (if one day in the future someone stumbles across this answer), the integral can be modified slightly to yield a very nice result:
In Owen, there is result 10,010.7:
$$ \int_{-\infty}^{0} \phi(x) \Phi(x+c)dx = \frac{1}{2}\Phi\left( \frac{c}{\sqrt{2}} \right)^2 $$
Sadly, a change of variable will not work in this case since it would mess up the limits.
