The exponential probability distribution is very well understood and characterized. While there are many ways to derive the moments of the exponential distribution, there is one I enjoy in particular, so I am sharing it below. It relies on two simple concepts: integration by parts and recursion.

Problem setup

The exponential distribution with parameter $$\lambda$$ is described by the following probability density function (PDF): $$f_X(x) = \lambda e^{-\lambda x}, \; x \geq 0$$ and by the following cumulative distribution function (CDF): $$F_X(x) = 1-e^{-\lambda x}$$. The $$n^{th}$$ moment of this distribution is given by $$M_n = \int_0^\infty x^n f_X(x) dx$$ (by definition). We will show here that $$M_n = {n! \over \lambda^n}$$.

Derivation

$$M_n = \int_0^\infty x^n f_X(x) dx = \int_0^\infty \lambda x^n e^{-\lambda x} dx$$. 

Recall the integration by parts rule: $$\int_a^b u(x)dv(x) = u(b)v(b)-u(a)v(a) - \int_a^b v(x)du(x)$$. To compute $$M_n$$, let's set $$u(x) = x^n$$ (implying $$du(x) = nx^{n-1} dx$$), $$v(x) = -e^{-\lambda x}$$ (implying $$dv(x) = \lambda e^{-\lambda x} dx$$), $$a=0$$, and $$b=\infty$$ (with a slight abuse of notation). Substituting, we get:

$$M_n = -\lim_{ x \to \infty} x^n e^{-\lambda x} - u(0)v(0) + \int_0^\infty nx^{n-1}e^{-\lambda x} dx$$

It can be shown by applying L'Hospital's rule $$n$$ times to $${x^n \over e^{\lambda x}}$$ that $$\lim_{ x \to \infty} x^n e^{-\lambda x}=0$$. Alternatively, the limit can be derived by using the Taylor series expansion of $$e^{\lambda x}$$. Also, since $$u(0)=0$$ and $$v(0)=-1$$, we have $$u(0)v(0)=0$$. This leaves us with:

$$M_n = \int_0^\infty nx^{n-1}e^{-\lambda x} dx = {n \over \lambda} \underbrace{\int_0^\infty \lambda x^{n-1}e^{-\lambda x} dx}_{M_{n-1}}$$, i.e. the recursion $$M_n = {n \over \lambda}M_{n-1}$$. Applying the recursion repeatedly we get:

$$M_n = {n \over \lambda}M_{n-1} = {n \over \lambda}{n-1 \over \lambda}M_{n-2} = {n \over \lambda}{n-1 \over \lambda}{n-2 \over \lambda}M_{n-3} \; =...= \; \frac{n!}{\lambda^n}M_0$$. 

By definition, $$M_0 = \int_0^\infty f_X(x) dx = 1$$, since $$f_X(x)$$ is a probability density function and must integrate to $$1$$. Substituting $$M_0=1$$, we have our desired result, viz.

$$M_n = \frac{n!}{\lambda^n}$$.

Fairly simple derivation, but a couple of tricks you may find handy when manipulating integrands involving polynomials and exponentials.