Reading my script from statistics, I have come across the following statement:
Consider the following random variables $$X_1,...,X_n$$ which are all independent and identically distributed exponential random variables with parameter $\theta$. Then, the distribution of $$\frac{1}{n} \sum^{n}_{i=1}X_{i}$$ is $\text{Gamma}(n,n\theta)$.
I may believe in that, but I have no clue how to calculate it and see for myself. I will be glad for help.
If $X_1,...,X_n\sim\mathrm{Exp}(\theta)$ and $S=X_1+...+X_n$ then: $$f_S(x)=f_{X_1}(x)*...*f_{X_1}(x)$$ Considering first just $X_1$ and $X_2$, the probability density function is given by the convolution: $$f_{X_1+X_2}(x)=\int_{\mathbb{R}}f_{X_1}(\xi)f_{X_2}(x-\xi)\mathrm d \xi=\int_0^xf_{X_1}(\xi)f_{X_2}(x-\xi)\mathrm d \xi=\int_0^x \theta e^{-\theta \xi} \theta e^{-\theta (x-\xi)}\mathrm d \xi=\theta^2xe^{-\theta x}$$ Now, considering also $X_3$, in $(X_1+X_2)+X_3$: $$f_{X_1+X_2+X_3}(x)=\int_0^xf_{X_1+X_2}(\xi)f_{X_3}(x-\xi)\mathrm d \xi=\int_0^x \theta^2\xi e^{-\theta \xi} \theta e^{-\theta (x-\xi)}\mathrm d \xi=\theta^3\frac{x^2}{2}e^{-\theta x}$$ Iterating this reasoning, you get: $$f_{X_1+...+X_n}(x)=\theta^n \frac{x^{n-1}}{(n-1)!}e^{-\theta x}=\frac{\theta e^{-\theta x} (\theta x)^{n-1}}{\Gamma(n)}$$ where $\Gamma(n)=(n-1)!$
This last expression is exactly the the p.d.f. of the Erlang distribution ( that is the Gamma distribution $\Gamma(\alpha, \theta)$ for $\alpha=n\in\mathbb{N}$).
Since you are looking for the average value, and it is simple to derive that, given $Y=kX$ then $f_Y(x)=\frac{1}{k}f_X(\frac{x}{k})$ therefore, in order to get the p.d.f. of $Y=\frac{1}{n}\sum\limits_{i=1}^n X_i$ you have just to apply this property where $k=\frac{1}{n}$:
$$f_Y(x) = n\theta^n \frac{(nx)^{n-1}}{(n-1)!}e^{-n\theta x}=(n\theta)^n \frac{x^{n-1}}{(n-1)!}e^{-n\theta x}=$$ which is exactly the result you are looking for, i.e. $Y \sim\Gamma(n,n\theta)$
Hint:
The CF (characteristic function) of sum of independent random variables is equal to the product of their individual CF.
Compare the CF of the exponential RV with that of the Gamma RV here. What can you conclude?
$\newcommand{\bbx}[1]{\,\bbox[8px,border:1px groove navy]{{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
The coveted PDF is given by the following multiple integration:
\begin{align} &\int_{0}^{\infty}\theta\expo{-\theta x_{1}}\ldots \int_{0}^{\infty}\theta\expo{-\theta x_{n}} \delta\pars{x - {1 \over n}\sum_{\ell = 1}^{n}x_{\ell}}\dd x_{1}\ldots\dd x_{n} \label{1}\tag{1} \end{align} Since integrations are performed over $\ds{\pars{0,\infty}^{\, n}}$, it's convenient to use a Dirac Delta Laplace Integral Representation. Namely, \begin{equation} \delta\pars{\xi} \equiv \int_{c -\infty\ic}^{c + \infty\ic}\expo{\xi s}\, {\dd s \over 2\pi\ic}\,,\qquad c > 0 \end{equation} Expression \eqref{1} becomes: \begin{align} &\int_{0}^{\infty}\theta\expo{-\theta x_{1}}\ldots \int_{0}^{\infty}\theta\expo{-\theta x_{n}} \int_{c -\infty\ic}^{c + \infty\ic} \exp\pars{s\bracks{x - {1 \over n}\sum_{\ell = 1}^{n}x_{\ell}}} \,{\dd s \over 2\pi\ic}\,\dd x_{1}\ldots\dd x_{n} \\[5mm] = &\ \int_{c -\infty\ic}^{c + \infty\ic}\expo{xs} \pars{\int_{0}^{\infty}\theta\expo{-\theta\xi}\expo{-s\xi/n}\,\dd\xi}^{n}\, {\dd s \over 2\pi\ic}\label{2}\tag{2} \end{align}
Note that $\ds{\int_{0}^{\infty}\theta\expo{-\theta\xi}\expo{-s\xi/n}\,\dd\xi = {n\theta \over s + n\theta}}$. Then, expression \eqref{2} becomes:
\begin{align} &\int_{c - \infty\ic}^{c + \infty\ic}\expo{xs} \pars{n\theta \over s + n\theta}^{n}\,{\dd s \over 2\pi\ic} = \pars{n\theta}^{n}\expo{-n\theta x}\ \overbrace{% \int_{c + n\theta - \infty\ic}^{c + n\theta + \infty\ic} {\expo{xs} \over k^{n}}\,{\dd s \over 2\pi\ic}} ^{\ds{x^{n - 1} \over \pars{n - 1}!}} = \bbox[#ffe,10px,border:1px dotted navy]{\ds{% {\pars{n\theta}^{n} \over \pars{n - 1}!}\,x^{n - 1} \expo{-n\theta x}}} \end{align}
