Get the raw moments via the factorial moments

As with many other discrete distributions with simple factorial moments, obtaining the high-order raw moments is simplest when done through the factorial moments. For the Poisson distribution, the factorial moments have an extremely simple form:

$$\begin{align} \mathbb{E}((X)_r) &= \sum_{x=0}^\infty (x)_r \cdot \text{Pois}(x|\lambda) \\[6pt] &= \sum_{x=r}^\infty (x)_r \cdot \text{Pois}(x|\lambda) \\[6pt] &= \sum_{x=r}^\infty (x)_r \cdot \frac{\lambda^x}{x!} \cdot \exp(-\lambda) \\[6pt] &= \sum_{x=r}^\infty \frac{\lambda^x}{(x-r)!} \cdot \exp(-\lambda) \\[6pt] &= \lambda^r \sum_{x=r}^\infty \frac{\lambda^{x-r}}{(x-r)!} \cdot \exp(-\lambda) \\[6pt] &= \lambda^r \sum_{x=0}^\infty \frac{\lambda^x}{x!} \cdot \exp(-\lambda) \\[6pt] &= \lambda^r \sum_{x=0}^\infty \text{Pois}(x|\lambda) \\[6pt] &= \lambda^r. \\[6pt] \end{align}$$$$\begin{align} \mathbb{E}((X)_r) &= \sum_{x=0}^\infty (x)_r \cdot \text{Pois}(x|\lambda) \\[6pt] &= \sum_{x=r}^\infty (x)_r \cdot \text{Pois}(x|\lambda) \\[6pt] &= \sum_{x=r}^\infty (x)_r \cdot \frac{\lambda^x}{x!} \cdot \exp(-\lambda) \\[6pt] &= \sum_{x=r}^\infty \frac{\lambda^x}{(x-r)!} \cdot \exp(-\lambda) \\[6pt] &= \lambda^r \sum_{x=r}^\infty \frac{\lambda^{x-r}}{(x-r)!} \cdot \exp(-\lambda) \\[6pt] &= \lambda^r \sum_{x=0}^\infty \frac{\lambda^x}{x!} \cdot \exp(-\lambda) \\[6pt] &= \lambda^r \sum_{x=0}^\infty \text{Pois}(x|\lambda) \\[12pt] &= \lambda^r. \\[6pt] \end{align}$$

Using the relationship between the raw moments and the factorial moments, we get the following general formula for the raw moments of the Poisson distribution:

$$\begin{align} \mathbb{E}(X^k) &= \sum_{r=0}^k S(k,r) \cdot \mathbb{E}((X)_r) \\[6pt] &= \sum_{r=0}^k S(k,r) \cdot \lambda^r. \\[6pt] \end{align}$$$$\begin{align} \mathbb{E}(X^k) &= \sum_{r=0}^k S(k,r) \cdot \mathbb{E}((X)_r) \\[6pt] &= \sum_{r=0}^k S(k,r) \cdot \lambda^r, \\[6pt] \end{align}$$

where the terms $S(k,r)$ denote the Stirling numbers of the second kind. The fourth raw moment for the Poisson distribution is therefore given by:

$$\mathbb{E}(X^4) = \lambda^4 + 6 \lambda^3 + 7 \lambda^2 + \lambda.$$