Difference between two i.i.d Laplace distributions? What is the PDF of the difference of two i.i.d Laplace distributed random variables?
I know that the difference of two i.i.d Normal variables is still the Normal distribution. Since the properties of the Laplace distribution are similar to the Normal distribution, I am guessing that the difference is also the Laplace distribution.
I have tried to solve this using Mathematica. Neither PDF nor Integrate functions give a correct answer. Thank you!
 A: The density of the difference of independent random variables $X$ and $Y$ is
given by $$f_{X-Y}(z) = \int_{-\infty}^\infty f_X(x)f_Y(x-z)\,\mathrm dx.$$
For iid Laplacian random variables, with $z>0$, this becomes
$$\begin{align}4f_{X-Y}(z) &= \int_{-\infty}^\infty \exp(-|x|)\exp(-|x-z|)\,\mathrm dx\\
&=\int_{-\infty}^0 \exp(x)\exp(x-z)\,\mathrm dx
+ \int_0^z \exp(-x)\exp(x-z)\,\mathrm dx\\
&\qquad\qquad\qquad+\int_z^\infty \exp(-x)\exp(-(x-z))\,\mathrm dx\\
&= \frac 12\exp(-z) + z\exp(-z)+\frac 12 \exp(-2z)\exp(z)\\
&= \exp(-z)+z\exp(-z).
\end{align}$$
Since the symmetry about the origin of the pdfs of $X$ and $Y$ implies
that pdf of their difference is also symmetric about the origin, we get that
$$\begin{align}
f_{X-Y}(z) &= \left.\left.\frac 14\right[\exp(-|z|) + |z|\exp(-|z|)\right],
\quad -\infty < z < \infty \tag{1}\\
&= \frac 12\left[\frac 12\exp(-|z|) + \frac 12|z|\exp(-|z|)\right], 
\quad -\infty < z < \infty \tag{2}
\end{align}$$
where in $(2)$ we can recognize that on the positive real line,
the density is one-half of a mixture density with equal weights,
where the mixture
is that of an exponential random variable and a Gamma random 
variable with shape parameter $2$, and of course, this is all
reflected on the negative real line.

Alternatively, as pointed out in my comment on the question,
for $z > 0$, a different answer of mine
arrives at 
$$\begin{align}
f_{X-Y}(z) &= \int_{0}^\infty (y+z)\exp(-(y+z))\cdot y\exp(-y)\,\mathrm dy\\
&= \exp(-z)\int_0^\infty (y^2+zy)\exp(-2y)\,\mathrm dy\\
&= \frac 12\exp(-z)\int_0^\infty \left(\frac{x^2}{4}+z\frac{x}2\right)\exp(-x)
\,\mathrm dx\\
&= \frac 12\exp(-z) \left(\frac{\Gamma(3)}{4}+z\frac{\Gamma(2)}{2}\right)\\
&= \left.\left.\frac 14 \right[\exp(-z) + z\exp(-z)\right]
\end{align}$$
by recognizing (as pointed out by Michael Mayer) that a Laplacian density
is the density of the difference of two exponential random variables, and
so regrouping $$(X_1-X_2)-(X_3-X_4) = (X_1+X_4)-(X_2+X_3)$$ allows us
to compute the desired density as that of the difference of two Gamma
random variables of shape parameter $2$.
