Below is a summary of the answer to a similar question I posted on Mathematics Stack Exchange: [https://math.stackexchange.com/a/4891023/1307231](https://math.stackexchange.com/a/4891023/1307231).

Basically, if the Laplace distribution can be construed as the distribution of a difference of independent exponential variates with equal scales $\lambda$, then it only makes sense that a sum of $N$ Laplace rvs is distributed as a difference of independent gamma random variables with shape parameters $N$ and scale $\lambda$.

This is my derivation of the distribution:

$$
\begin{array}{ccl}
f_{S_{N}}\!\left(s\right) & = & \int_{\max\left(s,0\right)}^{\infty}\tfrac{1}{\lambda^{N}\left(N-1\right)!}\left(x-s\right)^{N-1}\exp\!\left(-\tfrac{x-s}{\lambda}\right)\cdot\tfrac{1}{\lambda^{N}\left(N-1\right)!}x^{N-1}\exp\!\left(-\tfrac{x}{\lambda}\right)\mathrm{d}x\\
 & = & \tfrac{1}{\lambda^{2N}\left\{ \left(N-1\right)!\right\} ^{2}}\int_{0}^{\infty}\left\{ \left(x+\max\left(s,0\right)\right)\left(x+\max\left(s,0\right)-s\right)\right\} {}^{N-1}\exp\!\left\{ -\tfrac{2\left(x+\max\left(s,0\right)\right)-s}{\lambda}\right\} \mathrm{d}x\\
 & = & \tfrac{1}{\lambda^{2N}\left\{ \left(N-1\right)!\right\} ^{2}}\int_{0}^{\infty}\left\{ \left(x+\tfrac{1}{2}\left|s\right|\right)^{2}-\left(\tfrac{s}{2}\right)^{2}\right\} {}^{N-1}\exp\!\left\{ -\tfrac{1}{\lambda/2}\left(x+\tfrac{1}{2}\left|s\right|\right)\right\} \mathrm{d}x\\
 & = & \tfrac{1}{\lambda^{2N}\left\{ \left(N-1\right)!\right\} ^{2}}\int_{0}^{\infty}\left\{ x\left(x+\left|s\right|\right)\right\} {}^{N-1}\exp\!\left\{ -\tfrac{1}{\lambda/2}\left(x+\tfrac{1}{2}\left|s\right|\right)\right\} \mathrm{d}x\\
 & = & \tfrac{1}{\lambda^{2N}\left\{ \left(N-1\right)!\right\} ^{2}}\exp\!\left\{ -\tfrac{1}{\lambda}\left|s\right|\right\} \int_{0}^{\infty}x^{N-1}\sum_{i=0}^{N-1}\binom{N-1}{i}\left|s\right|^{i}x^{N-1-i}\exp\!\left(-\tfrac{1}{\lambda/2}x\right)\mathrm{d}x\\
 & = & \tfrac{1}{\lambda^{2N}\left\{ \left(N-1\right)!\right\} ^{2}}\exp\!\left\{ -\tfrac{1}{\lambda}\left|s\right|\right\} \sum_{i=0}^{N-1}\binom{N-1}{i}\left|s\right|^{i}\left(\tfrac{\lambda}{2}\right)^{2N-i-1}\left(2N-i-2\right)!\\
 & = & \tfrac{1}{2^{2N-1}\lambda}\exp\!\left\{ -\tfrac{1}{\lambda}\left|s\right|\right\} \sum_{i=0}^{N-1}\binom{2N-i-2}{N-1}\tfrac{\left(\tfrac{2}{\lambda}\left|s\right|\right)^{i}}{i!}
\end{array}
$$

Beyond the demonstration itself, this general formula matches results shown in Kotz *et al.* (2001) for $N = 1, 2, 3, 4$.

The distribution can then easily be shifted and rescaled to match the sample mean distribution, say.

Hopefully this clarifies things a little bit.