Gamma Distribution out of sum of exponential random variables I have a sequence $T_1,T_2,\ldots$ of independent exponential random variables with paramter $\lambda$. I take the sum $S=\sum_{i=1}^n T_i$ and now I would like to calculate the probability density function.
Well, I know that $P(T_i>t)=e^{-\lambda t}$ and therefore $f_{T_i}(t)=\lambda e^{-\lambda t}$ so I need to find $P(T_1+\cdots+T_n>t)$ and take the derivative. But I cannot expand the probability term, you have any ideas?
 A: The usual way to do this is to consider the moment generating function, noting that if $S = \sum_{i=1}^n X_i$ is the sum of IID random variables $X_i$, each with MGF $M_X(t)$, then the MGF of $S$ is $M_S(t) = (M_X(t))^n$.  Applied to the exponential distribution, we can get the gamma distribution as a result.
If you don't go the MGF route, then you can prove it by induction, using the simple case of the sum of the sum of a gamma random variable and an exponential random variable with the same rate parameter.  Let's actually do this.  Suppose $Y \sim {\rm Gamma}(a,b)$ and $X \sim {\rm Exponential}(b)$ are independent, so that $$f_Y(y) = \frac{b^a y^{a-1} e^{-by}}{\Gamma(a)} \mathbb 1(y > 0), \quad f_X(x) = be^{-bx} \mathbb 1(x > 0), \quad a, b > 0.$$  Then, we notice that if $a = 1$, $Y$ would also be exponential (i.e., the exponential distribution is a special case of the Gamma with $a = 1$).  Now consider $Z = X+Y$.  The PDF is $$\begin{align*} f_Z(z) &= \int_{y=0}^z f_Y(y) f_X(z-y) \, dy \\ &= \int_{y=0}^z \frac{b^{a+1} y^{a-1} e^{-by} e^{-b(z-y)}}{\Gamma(a)} \, dy \\ &= \frac{b^{a+1} e^{-bz}}{\Gamma(a)} \int_{y=0}^z y^{a-1} \, dy \\ &= \frac{b^{a+1} e^{-bz}}{\Gamma(a)} \cdot \frac{z^a}{a} = \frac{b^{a+1} z^a e^{-bz}}{\Gamma(a+1)}. \end{align*}$$  But this is just a gamma PDF with new shape parameter $a^* = a+1$.  So, it is easy to see by induction that the sum of $n$ IID exponential variables with common rate parameter $\lambda$ is gamma with shape parameter $a = n$, and rate parameter $b = \lambda$.
