Formula for the nth convolution of a Laplace random variable

Question

Let $x_1, x_2, ...$ be i.i.d. draws from a Laplace distribution with scale parameter $b$. Is there a relatively nice closed form for $x_1+x_2+...x_n$? I've seen a derivation floating around for when $b=1$, but I couldn't figure out a generalization.

Answer 1

The distribution of the $n$-th convolution of the Laplace distribution can be computed from the characteristic function (see on Wikipedia): $$\frac{\exp(i \mu t)}{1+b^2 t^2} \,.$$ The characteristic function of the $n$-th convolution becomes: $$\frac{\exp(i n \mu t)}{(1+b^2 t^2)^n} = \frac{\exp(i n \mu t)}{(1 - i b t)^n (1 + i b t)^n} \,.$$ The inverse Fourier transform can be computed using the residue theorem. The integration contour is closed from the upper or lower half plane according to the sign of $(x-n \mu)$.

Answer 2

One way to generate a Laplace random variable is to generate two IID (independent and identically distributed) exponential random variables and then subtract them: x_i = y_i - z_i with y_i and z_i ~ exponential(parameter=b), and of course everything independent. Then the sum of the x_i is simply (sum y_i) - (sum z_i); each of those two sums have Gamma distributions. To be more specific, since we are summing an integer number of terms, they have Erlang distributions. The difference of two Gammas is called "bilateral gamma", and there are a few papers out there on it. A quick search just found:

Bilateral gamma distributions and processes in financial mathematics Uwe Küchlera, Stefan Tappe

On the shapes of bilateral Gamma densities Uwe Küchlera, Stefan Tappe

It would be nice if someone would write a Wikipedia article about bilateral Gammas, I guess.

Answer 3

Below is a summary of the answer to a similar question I posted on Mathematics Stack Exchange: 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.