Density of sum of truncated normal and normal distribution Suppose that $\varepsilon\sim N(0, \sigma_\varepsilon)$ and $\delta\sim N^+(0, \sigma_\delta)$. What is the density function for $X = \varepsilon - \delta$?
This proof apparently appeared in a Query by M.A. Weinstein in Technometrics 6 in 1964, which stated that the density of $X$ is given by
$$f_X(x) = \frac{2}{\sigma} \phi\left(\frac{x}{\sigma}\right) \left(1 - \Phi\left(\frac{x\lambda}{\sigma}\right)\right),$$
where $\sigma^2 = \sigma_\varepsilon^2 + \sigma_\delta^2$ and $\lambda = \sigma_\delta / \sigma_\varepsilon$ and $\phi$ and $\Phi$ are the standard normal density and distribution functions, respectively. However, that paper is very difficult to find online. What is the proof that the density of $X$ takes the above form?
 A: The zero-truncated, zero-mean variable $\delta\sim N^+(0,\sigma^2)$ can be seen as a limiting case of a skew-normal variable since its pdf
$$
f_\delta(\delta)
 =\lim_{\alpha_\delta\rightarrow\infty}\frac2{\sigma_\delta}\phi(\frac\delta{\sigma_\delta})\Phi(\alpha_\delta\sigma_\delta\delta)
$$
(except at $\delta=0$ but this do not matter when $\delta$ is continuously distributed).
Its mgf is given by the same limit of the skew-normal mgf,
\begin{align}
M_\delta(t) 
  &=\lim_{\alpha_\delta\rightarrow\infty}2 e^{\sigma_\delta^2 t^2/2}\Phi(\frac{\alpha_\delta}{\sqrt{1+\alpha_\delta^2}}\sigma_\delta t) 
\\&= e^{t^2/2}\Phi(\sigma_\delta t).
\end{align}
The rest of the proof is as in the post I link to in the comments.  The mgf of $\varepsilon$ is $M_\varepsilon(t)=e^{\sigma_\epsilon^2 t^2}$ and, because of independence, the mgf of $X = \varepsilon-\delta$ is
\begin{align}
M_{\varepsilon-\delta}(t)
  &=Ee^{(\varepsilon-\delta)t}
\\&=Ee^{\varepsilon t}Ee^{-t\delta}
\\&=M_\varepsilon(t)M_\delta(-t)
\\&=2e^{(\sigma_\delta^2+\sigma_\varepsilon^2) t^2/2}\Phi(-\sigma_\delta t).
\end{align}
This equals the mgf of a skew-normal,
$$
2e^{\mu t+\sigma^2t^2/2}\Phi(\sigma\frac{\alpha}{\sqrt{1+\alpha^2}}t),
$$
with location parameter $\mu=0$, scale parameter $\sigma=\sqrt{\sigma_\delta^2+\sigma_\varepsilon^2}$ and shape parameter $\alpha$ satisfying
$$
-\sigma_\delta=\sqrt{\sigma_\delta^2+\sigma_\varepsilon^2}\frac{\alpha}{\sqrt{1+\alpha^2}}
$$
that is, $\alpha=-\sigma_\delta/\sigma_\varepsilon$.  
A: Ultimately I needed to work through the algebra a bit more to arrive at the specified form. For posterity, the full proof is given below.
Proof
First consider the distribution function of $X$, which is given by
    $$F(x) = \Pr(X \leq x) = \Pr(\varepsilon - \delta \leq x)$$
$$= \int_{\varepsilon - \delta \leq x} f_\varepsilon(\varepsilon) f_\delta(\delta) d\delta d\varepsilon$$
$$= \int_{\delta\in\mathbb{R}^+} f_\delta(\delta) \int_{\varepsilon\in (-\infty, x + \delta]} f_\varepsilon(\varepsilon) d\varepsilon d\delta.$$
    Substituting in known density functions yields
    $$\int_0^\infty 2\phi(\delta | 0, \sigma_\delta) \int_{-\infty}^{x + \delta} \phi(\varepsilon | 0, \sigma_\varepsilon) d\varepsilon d\delta$$
$$= 2\int_0^\infty \phi(\delta | 0, \sigma_\delta) \Phi(x + \delta | 0, \sigma_\varepsilon) d\delta.$$
    The density of $X$ is then given by
    $$f(x) = \frac{dF}{dx} = 2\int_0^\infty \phi(\delta | 0, \sigma_\delta) \phi(x + \delta | 0, \sigma_\varepsilon) d\delta.$$
    Using Sage to perform this integration, the result is given by
    $$f(x) = -\frac{{\left(\operatorname{erf}\left(\frac{\sigma_{\delta} x}{2 \, \sqrt{\frac{1}{2} \, \sigma_{\delta}^{2} + \frac{1}{2} \, \sigma_{\varepsilon}^{2}} \sigma_{\varepsilon}}\right) e^{\left(\frac{\sigma_{\delta}^{2} x^{2}}{2 \, {\left(\sigma_{\delta}^{2} \sigma_{\varepsilon}^{2} + \sigma_{\varepsilon}^{4}\right)}}\right)} - e^{\left(\frac{\sigma_{\delta}^{2} x^{2}}{2 \, {\left(\sigma_{\delta}^{2} \sigma_{\varepsilon}^{2} + \sigma_{\varepsilon}^{4}\right)}}\right)}\right)} e^{\left(-\frac{x^{2}}{2 \, \sigma_{\varepsilon}^{2}}\right)}}{2 \, \sqrt{\pi} \sqrt{\frac{1}{2} \, \sigma_{\delta}^{2} + \frac{1}{2} \, \sigma_{\varepsilon}^{2}}}.$$
    Defining $\lambda = \sigma_\delta / \sigma_\varepsilon$ and $\sigma^2 = \sigma_\varepsilon^2 + \sigma_\delta^2$, the following can be simplified:
    $$\frac{\sigma_{\delta} x}{2 \, \sqrt{\frac{1}{2} \, \sigma_{\delta}^{2} + \frac{1}{2} \, \sigma_{\varepsilon}^{2}} \sigma_{\varepsilon}} = \frac{\lambda x}{\sigma\sqrt{2}} = \frac{x}{(\sigma / \lambda) \sqrt{2}},$$
$$\frac{\sigma_{\delta}^{2} x^{2}}{2 \, {\left(\sigma_{\delta}^{2} \sigma_{\varepsilon}^{2} + \sigma_{\varepsilon}^{4}\right)}} = \frac{\lambda^2 x^2}{2\sigma^2} = \frac{x^2}{2(\sigma / \lambda)^2}.$$
    Thus,
    $$f(x) = -\frac{\exp\left(\frac{x^2}{2(\sigma / \lambda)^2}\right)\left(\operatorname{erf}\left(\frac{x}{(\sigma / \lambda) \sqrt{2}}\right) - 1\right)\exp\left(-\frac{x^2}{2\sigma_\varepsilon^2}\right)}{\sigma\sqrt{2\pi}}$$
$$= -\frac{\left(\operatorname{erf}\left(\frac{x}{(\sigma / \lambda) \sqrt{2}}\right) - 1\right) \exp\left(-x^2\left(\frac{1}{2\sigma_\varepsilon^2} - \frac{1}{2(\sigma / \lambda)^2}\right)\right)}{\sigma\sqrt{2\pi}}.$$
    Now,
    $$\operatorname{erf}\left(\frac{x}{(\sigma / \lambda) \sqrt{2}}\right) - 1 = \left(1 + \operatorname{erf}\left(\frac{x}{(\sigma / \lambda) \sqrt{2}}\right)\right) - 2 = 2\left(\frac{1}{2}\left(1 + \operatorname{erf}\left(\frac{x}{(\sigma / \lambda) \sqrt{2}}\right)\right) - 1\right)$$
$$= 2\left(\Phi\left(\frac{x\lambda}{\sigma}\right) - 1\right) = -2\left(1 - \Phi\left(\frac{x\lambda}{\sigma}\right)\right).$$
    Also,
    $$\frac{1}{2\sigma_\varepsilon^2} - \frac{1}{2(\sigma / \lambda)^2} = \frac{1}{2\sigma_\varepsilon^2} - \frac{\sigma_\delta^2}{2\sigma_\varepsilon^2(\sigma_\delta^2 + \sigma_\varepsilon^2)} = \frac{\sigma_\delta^2 + \sigma_\varepsilon^2 - \sigma_\delta^2}{2\sigma_\varepsilon^2(\sigma_\delta^2 + \sigma_\varepsilon^2)} = \frac{\sigma_\varepsilon^2}{2\sigma_\varepsilon^2(\sigma_\delta^2 + \sigma_\varepsilon^2)} = \frac{1}{2\sigma^2}.$$
    So,
    $$f(x) = 2\left(1 - \Phi\left(\frac{x\lambda}{\sigma}\right)\right)\frac{\exp\left(-\frac{x^2}{2\sigma^2}\right)}{\sigma\sqrt{2\pi}} = 2\left(1 - \Phi\left(\frac{x\lambda}{\sigma}\right)\right) \phi(x | 0, \sigma) = \frac{2}{\sigma}\phi\left(\frac{x}{\sigma}\right) \left(1 - \Phi\left(\frac{x\lambda}{\sigma}\right)\right).$$
