I want to show that $$Y_1 + Y_2 \overset{d}= X_1 - X_2$$ where $Y_1, Y_2 \sim \text{Laplace}(1)$ and $X_1,X_2 \sim \Gamma(2,1)$, by checking their Moment Generating Functions.
For the left-hand side, since $Y_1, Y_2$ are iid, I use the fact that $S_{n=2}=Y_1+Y_2$, and that $$\psi_{S_{n=2}}(s)=(\psi_Y(t))^{n=2}=\left(\frac 1 {1+a^2t^2}\right)^{n=2}=\left[L(a=1)\right]=\left(\frac 1 {1+t^2}\right)^{2}$$ However I am not sure how to solve the right hand side, i.e. $X_1 - X_2$, nor how to find its MGF. Any help would be much appreciated!
A Laplacian distribution with parameter $\gamma$ can be written as $\gamma (E_1 - E_2)$ where $E_1, E_2$ are independent standard exponential random variables. This is easy to show by matching their MGF's, or by directly computation.
Thus, the left hand side is $E_1 - E_2 + E_3 -E_4$ where $E_1, E_2, E_3, E_4$ are i.i.d. Exp(1) random variables. Lump $E_1+ E_3$ together (this is $X_1$) and $E_2 + E_4$ together (this is $X_2$) to get the result.
Alternatively, you can note that the MGF of the right hand side (with $X_1, X_2$ independent) is given by $E[e^{t (X_1 - X_2)}] = E[e^{t X_1} e^{-t X_2}] = M_{X_1}(t) M_{X_2}(-t)$ and write this out by plugging in the MGF of a $\Gamma(2,1)$ distribution, and see that this matches the MGF of the LHS.