According to the Wikipedia article on the Gamma distribution:
If $X\sim\mathrm{Gamma}(a,\theta)$ and $Y\sim\mathrm{Gamma}(b,\theta)$, where $X$ and $Y$ are independent random variables, then $X+Y\sim \mathrm{Gamma}(a+b, \theta)$.
But I don't see any proof. Can anyone point me to its proof please?
Edit: Thanks to Zen a lot, and also I found the answer as an example on the Wikipedia page about characteristic functions.
The proof is as follows: (1) Remember that the characteristic function of the sum of independent random variables is the product of their individual characteristic functions; (2) Get the characteristic function of a gamma random variable here; (3) Do the simple algebra.
To get some intuition beyond this algebraic argument, check whuber's comment.
Note: The OP asked how to compute the characteristic function of a gamma random variable. If $X\sim\mathrm{Exp}(\lambda)$, then (you can treat $i$ as an ordinary constant, in this case)
$$\psi_X(t)=\mathrm{E}\left[e^{itX}\right]=\int_0^\infty e^{itx} \lambda\,e^{-\lambda x}\,dx = \frac{1}{1-it/\lambda}\, .$$
Now use Huber's tip: If $Y\sim\mathrm{Gamma}(k,\theta)$, then $Y=X_1+\dots+X_k$, where the $X_i$'s are independent $\mathrm{Exp}(\lambda = 1/\theta)$. Therefore, using property (1), we have $$ \psi_Y(t) = \left( \frac{1}{1-it\theta}\right)^k \, . $$
Tip: you won't learn these things staring at the results and proofs: stay hungry, compute everything, try to find your own proofs. Even if you fail, your appreciation of somebody else's answer will be at a much higher level. And, yes, failing is OK: nobody is looking! The only way to learn mathematics is by fist fighting for each concept and result.
Here is an answer that does not need to use characteristic functions, but instead reinforces some ideas that have other uses in statistics. The density of the sum of independent random variables is the convolutions of the densities. So, taking $\theta = 1$ for ease of exposition, we have for $z > 0$, $$\begin{align} f_{X+Y}(z) &= \int_0^z f_X(x)f_Y(z-x)\,\mathrm dx\\ &=\int_0^z \frac{x^{a-1}e^{-x}}{\Gamma(a)}\frac{(z-x)^{b-1}e^{-(z-x)}}{\Gamma(b)}\,\mathrm dx\\ &= e^{-z}\int_0^z \frac{x^{a-1}(z-x)^{b-1}}{\Gamma(a)\Gamma(b)}\,\mathrm dx &\scriptstyle{\text{now substitute}}~ x = zt~ \text{and think}\\ &= e^{-z}z^{a+b-1}\int_0^1 \frac{t^{a-1}(1-t)^{b-1}}{\Gamma(a)\Gamma(b)}\,\mathrm dt & \scriptstyle{\text{of Beta}}(a,b)~\text{random variables}\\ &= \frac{e^{-z}z^{a+b-1}}{\Gamma(a+b)} \end{align}$$
On a more heuristic level: If $a$ and $b$ are integers, the Gamma distribution is an Erlang distribution, and so $X$ and $Y$ describe the waiting times for respectively $a$ and $b$ occurrences in a Poisson process with rate $\theta$. The two waiting times $X$ and $Y$ are
and the waiting time for $a+b$ occurrences is distributed Gamma($a+b,\theta$).
None of this is a mathematical proof, but it puts some flesh on the bones of the connection, and can be used if you want to flesh it out in a mathematical proof.
