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?

• $\delta$ is skew normal (with infinite skew parameter). Hence $X$ is also skew normal, see stats.stackexchange.com/a/236734/77222, – Jarle Tufto Jul 29 at 22:01
• Why can't you apply the regular convolution formula for the density of a sum (or difference)? – Xi'an Jul 29 at 22:19
• $\epsilon$ and $\delta$ are independent? For truncated normal, $N^+$ means truncating at 0, so $0\le \delta \lt \infty$? – user158565 Jul 30 at 3:29
• @Xi'an I have tried doing that but haven't been able to arrive at the above form. If you can show the steps, please include as an answer. – tkmckenzie Jul 30 at 14:50
• @JarleTufto I agree with you, but am reluctant to close because there is nothing in the duplicate thread that explicitly shows this truncated Normal distribution is Skew normal. If there were, that thread would contain a full and clear answer. – whuber Jul 31 at 13:03

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$$.
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).$$