Show that $Y_1+Y_2$ have distribution skew-normal 
Let $Y_1\sim SN(\mu_1,\sigma_1^2,\lambda)$ and $Y_2\sim
 N(\mu_2,\sigma_2^2)$ independents. Show that $Y_1+Y_2$ have
  a skew-normal distribution and find the parameters of this distribution.

Since the random variables are independent I tried to use convolution. Let $Z=Y_1+Y_2$     
$$f_Z(z)=\int_{-\infty}^{\infty}2\phi(y_1|\mu_1,\sigma_1)\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\phi(z-y_1|\mu_2,\sigma_2^2)\,\text{d}y_1$$
Here $\phi()$ and $\Phi()$ are the standard normal pdf and cdf, respectively.
$$f_Z(z)=\int_{-\infty}^{\infty}2\frac{1}{\sqrt{2\pi\sigma_1}}\frac{1}{\sqrt{2\pi\sigma_2}}exp\Big(-\frac{1}{2\sigma_1^2}(y_1-\mu)^2-\frac{1}{2\sigma_2^2}((z-y_1)^2-\mu)^2\Big)\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\,\text{d}y_1$$
For simplified notations, let $k=2\frac{1}{\sqrt{2\pi\sigma_1}}\frac{1}{\sqrt{2\pi\sigma_2}}$
\begin{align*}f_Z(z)&=k\int_{-\infty}^{\infty}\exp\Big(\frac{-1}{2\sigma_1^2\sigma_2^2}\Big(\sigma_1^2(y_1-\mu_1)^2+\sigma_2^2((z-y_1)-\mu_2)^2\Big)\Big)\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\,\text{d}y_1\\
&=k\int_{-\infty}^{\infty}\exp\Big(\frac{-1}{2\sigma_1^2\sigma_2^2}\Big(\sigma_2^2(y_1^2-2y_1\mu_1+\mu_1)+\sigma_1^2((z-y_1)^2-2(z-y_1)\mu_2+\mu_2^2)\Big)\Big)\\&\quad\times\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\,\text{d}y_1=k\int_{-\infty}^{\infty}
\exp\\&\Big(\frac{-1}{2\sigma_1^2\sigma_2^2}\Big(\sigma_2^2(y_1^2-2y_1\mu_1+\mu_1)+\sigma_1^2(z^2-2zy_1+y_1^2-2z\mu_2+2y_1\mu_2+\mu_2^2)\Big)\Big)\\&\quad\times\Phi\Big(\lambda(\frac{y_1-\mu_1}{\sigma_1})\Big)\,\text{d}y_1
\end{align*}
But I'm stuck at this point.
EDIT: Following the suggestions in the comments, taking $\mu_1=\mu_2=0$ and $\sigma_1^2=\sigma_2^2=1$
\begin{align*}
&\int_{-\infty}^\infty 2\frac{1}{\sqrt{2\pi}}\frac{1}{\sqrt{2\pi}}\exp\Big(-\frac{1}{2}[y_1^2+z^2-2zy_1+y_1^2]\Big)\Phi(\lambda y_1)dy_1
\\&\int_{-\infty}^\infty 2\frac{1}{\sqrt{2\pi}}\frac{1}{\sqrt{2\pi}}\exp\Big(-\frac{1}{2}y_1^2\Big)\Phi(\lambda y_1) \exp\Big(-\frac{1}{2}(z-y_1)^2\Big)dy_1\end{align*}
is skew-normal.
 A: Reparameterizing the skew in terms of $\delta=\lambda/\sqrt{1+\lambda^2}$ and using the mgf of the skew normal (see below), since $Y_1$ and $Y_2$ are independent, $Z=Y_1 + Y_2$ has mgf
\begin{align}
M_Z(t) &= M_{Y_1}(t)M_{Y_2}(t) \\
       &= 2e^{\mu_1 t +\sigma_1^2 t^2/2}\Phi(\sigma_1\delta t)e^{\mu_2t +\sigma_2^2 t^2/2} \\
&= 2e^{(\mu_1+\mu_2)t + (\sigma_1^2+\sigma_2^2)t^2/2}\Phi(\sigma_1 \delta t) \\
&= 2e^{\mu t + \sigma^2 t^2/2}\Phi(\sigma \delta' t),
\end{align}
that is, the mgf of a skew normal with parameters $\mu=\mu_1+\mu_2$, $\sigma^2=\sigma_1^2+\sigma_2^2$ and $\sigma\delta'=\sigma_1\delta$ where $\delta'$ is the new skew parameter.  Hence,
$$
\delta'=\delta\frac{\sigma_1}\sigma=\delta\frac{\sigma_1}{\sqrt{\sigma_1^2+\sigma_2^2}}.
$$ 
In the other parameterization, the new skew parameter $\lambda'$ can be written, after some algebra, e.g. as
$$
\lambda' = \frac{\delta'}{\sqrt{1-\delta'^2}}=\frac{\lambda}{\sqrt{1 + \frac{\sigma_2^2}{\sigma_1^2}(1+\lambda^2)}}.
$$
The mgf of a standard skew normal can be derived as follows:
\begin{align}
M_X(t)&=Ee^{tX} \\
&=\int_{-\infty}^\infty e^{xt}2\frac1{\sqrt{2\pi}}e^{-x^2/2}\Phi(\lambda x)dx \\
&=2\int_{-\infty}^\infty \frac1{\sqrt{2\pi}}e^{-\frac12(x^2-2tx)}\Phi(\lambda x)dx\\
&=2\int_{-\infty}^\infty \frac1{\sqrt{2\pi}}e^{-\frac12((x-t)^2-t^2)}\Phi(\lambda x)dx \\
&=2e^{t^2/2} \int_{-\infty}^\infty \frac1{\sqrt{2\pi}}e^{-\frac12(x-t)^2}P(Z\le \lambda x)dx, & \text{where }Z \sim N(0,1) \\
&=2e^{t^2/2} P(Z\le \lambda U), & \text{where }U \sim N(t,1)\\
&=2e^{t^2/2} P(Z - \lambda U \le 0) \\
&=2e^{t^2/2} P(\frac{Z - \lambda U +\lambda t}{\sqrt{1+\lambda^2}} \le \frac{\lambda t}{\sqrt{1+\lambda^2}}) \\
&=2e^{t^2/2}\Phi(\frac\lambda{\sqrt{1+\lambda^2}}t).
\end{align}
The mgf of a skew normal with location and scale parameters $\mu$ and $\sigma$ is then
$$
M_{\mu + \sigma X}(t)=Ee^{(\mu+\sigma X)t} = e^{\mu t}M_X(\sigma t) = 2e^{\mu t+\sigma^2 t^2/2}\Phi(\frac\lambda{\sqrt{1+\lambda^2}}\sigma t).
$$
