[proofplan]
We differentiate the moment generating function $M_X(\theta) = \mathbb{E}[e^{\theta X}]$ under the integral sign $r$ times and evaluate at $\theta = 0$. The interchange of differentiation and expectation is justified by dominated convergence, using the finiteness of $M_X$ in a neighbourhood of the origin. At $\theta = 0$, the factor $e^{\theta X}$ becomes $1$, leaving $\mathbb{E}[X^r]$.
[/proofplan]
[step:Differentiate $M_X(\theta)$ under the expectation sign]
Suppose $M_X(\theta) = \mathbb{E}[e^{\theta X}]$ is finite for all $\theta \in (-\delta, \delta)$ for some $\delta > 0$. For such $\theta$, if $X$ has density $f_X$,
\begin{align*}
M_X(\theta) = \int_{-\infty}^{\infty} e^{\theta x} f_X(x)\, dx.
\end{align*}
We differentiate with respect to $\theta$. The partial derivative of the integrand is $\frac{\partial}{\partial \theta} e^{\theta x} = x\, e^{\theta x}$. To justify differentiating under the integral sign, fix $\theta_0 \in (-\delta, \delta)$ and choose $\varepsilon > 0$ such that $[\theta_0 - \varepsilon, \theta_0 + \varepsilon] \subset (-\delta, \delta)$. For all $\theta \in [\theta_0 - \varepsilon, \theta_0 + \varepsilon]$,
\begin{align*}
|x\, e^{\theta x}| \le |x|\, e^{|\theta_0 + \varepsilon|\, |x|} \le |x|\, e^{(\delta - \varepsilon/2)|x|},
\end{align*}
and the right-hand side is integrable against $f_X$ because $M_X(\delta - \varepsilon/2) < \infty$ ensures $\mathbb{E}[|X| e^{(\delta-\varepsilon/2)|X|}] < \infty$. By the [dominated convergence theorem](/theorems/4) (applied to difference quotients), differentiation under the integral is valid, giving
\begin{align*}
M_X'(\theta) = \int_{-\infty}^{\infty} x\, e^{\theta x}\, f_X(x)\, dx = \mathbb{E}[X e^{\theta X}].
\end{align*}
[/step]
[step:Iterate to obtain the $r$-th derivative]
The same dominated convergence argument applies at each stage: the $r$-th partial derivative of $e^{\theta x}$ with respect to $\theta$ is $x^r e^{\theta x}$, and for $\theta$ in a compact subinterval of $(-\delta, \delta)$, the bound $|x^r e^{\theta x}| \le |x|^r e^{(\delta - \varepsilon)|x|}$ is integrable against $f_X$ (since all moments of $X$ are finite when the mgf exists in a neighbourhood of zero). By induction,
\begin{align*}
M_X^{(r)}(\theta) = \mathbb{E}[X^r e^{\theta X}] \quad \text{for all } \theta \in (-\delta, \delta).
\end{align*}
[/step]
[step:Evaluate at $\theta = 0$ to recover $\mathbb{E}[X^r]$]
Setting $\theta = 0$, the exponential factor becomes $e^{0 \cdot X} = 1$:
\begin{align*}
M_X^{(r)}(0) = \mathbb{E}[X^r \cdot 1] = \mathbb{E}[X^r].
\end{align*}
[/step]