[proofplan]
We prove that derivatives may be passed through the expectation locally in $t$. Around a fixed point $t \in (-a,a)$, the two-sided exponential moment assumption gives integrable exponential bounds at parameters slightly larger and slightly smaller than $t$. Polynomial factors $|X|^m$ are absorbed into those exponential bounds, so the difference quotients are dominated by an integrable [random variable](/page/Random%20Variable). Dominated convergence then gives the derivative formula inductively, and evaluating at $t=0$ gives the moment identity.
[/proofplan]
[step:Build a local integrable dominator for polynomially weighted exponentials]
Fix $t \in (-a,a)$. Choose $\varepsilon>0$ such that $t-\varepsilon \in (-a,a)$ and $t+\varepsilon \in (-a,a)$, and define the neighbourhood
\begin{align*}
I_t := (t-\varepsilon/2,t+\varepsilon/2).
\end{align*}
For each integer $m \geq 0$, define the finite constant
\begin{align*}
C_{m,\varepsilon}:=\sup_{u \geq 0} u^m e^{-\varepsilon u/2}.
\end{align*}
This supremum is finite because the [continuous function](/page/Continuous%20Function) $u \mapsto u^m e^{-\varepsilon u/2}$ on $[0,\infty)$ tends to $0$ as $u \to \infty$.
For $s \in I_t$ and $x \geq 0$, we have $s \leq t+\varepsilon/2$, and therefore
\begin{align*}
x^m e^{sx} \leq C_{m,\varepsilon} e^{(t+\varepsilon)x}.
\end{align*}
For $s \in I_t$ and $x<0$, set $u=-x>0$. Since $s \geq t-\varepsilon/2$,
\begin{align*}
|x|^m e^{sx}=u^m e^{-su} \leq C_{m,\varepsilon} e^{-(t-\varepsilon)u}=C_{m,\varepsilon}e^{(t-\varepsilon)x}.
\end{align*}
Thus, for all $s \in I_t$ and all $x \in \mathbb R$,
\begin{align*}
|x|^m e^{sx} \leq C_{m,\varepsilon}\left(e^{(t+\varepsilon)x}+e^{(t-\varepsilon)x}\right).
\end{align*}
Define the measurable map $G_m:\Omega \to [0,\infty)$ by
\begin{align*}
G_m(\omega):=C_{m,\varepsilon}\left(e^{(t+\varepsilon)X(\omega)}+e^{(t-\varepsilon)X(\omega)}\right).
\end{align*}
Since $t+\varepsilon$ and $t-\varepsilon$ belong to $(-a,a)$, the hypothesis gives
\begin{align*}
\int_\Omega G_m(\omega)\,d\mathbb P(\omega)=C_{m,\varepsilon}\left(M_X(t+\varepsilon)+M_X(t-\varepsilon)\right)<\infty.
\end{align*}
Hence $G_m \in L^1(\Omega,\mathcal F,\mathbb P)$, and for every $s \in I_t$,
\begin{align*}
|X|^m e^{sX} \leq G_m
\end{align*}
$\mathbb P$-a.s.
[guided]
The main local problem is that differentiating $e^{sX}$ produces powers of $X$, and these powers are not assumed integrable in advance. The two-sided exponential moment hypothesis supplies exactly the missing control: a polynomial factor can be absorbed into a slightly larger exponential on the positive half-line and into a slightly smaller exponential on the negative half-line.
Fix $t \in (-a,a)$. Because $(-a,a)$ is open, choose $\varepsilon>0$ such that both $t-\varepsilon$ and $t+\varepsilon$ still lie in $(-a,a)$. Define
\begin{align*}
I_t := (t-\varepsilon/2,t+\varepsilon/2).
\end{align*}
For each integer $m \geq 0$, define
\begin{align*}
C_{m,\varepsilon}:=\sup_{u \geq 0} u^m e^{-\varepsilon u/2}.
\end{align*}
The function $u \mapsto u^m e^{-\varepsilon u/2}$ is continuous on $[0,\infty)$ and tends to $0$ as $u \to \infty$, so the supremum is finite.
Now let $s \in I_t$ and $x \in \mathbb R$. If $x \geq 0$, then $s \leq t+\varepsilon/2$, so
\begin{align*}
x^m e^{sx}=x^m e^{-(\varepsilon/2)x}e^{(s+\varepsilon/2)x} \leq C_{m,\varepsilon}e^{(t+\varepsilon)x}.
\end{align*}
If $x<0$, write $u=-x>0$. Since $s \geq t-\varepsilon/2$,
\begin{align*}
|x|^m e^{sx}=u^m e^{-su}=u^m e^{-(\varepsilon/2)u}e^{-(s-\varepsilon/2)u} \leq C_{m,\varepsilon}e^{-(t-\varepsilon)u}=C_{m,\varepsilon}e^{(t-\varepsilon)x}.
\end{align*}
Combining the two cases gives the uniform bound
\begin{align*}
|x|^m e^{sx} \leq C_{m,\varepsilon}\left(e^{(t+\varepsilon)x}+e^{(t-\varepsilon)x}\right)
\end{align*}
for every $s \in I_t$ and every $x \in \mathbb R$.
Define the random variable $G_m:\Omega \to [0,\infty)$ by
\begin{align*}
G_m(\omega):=C_{m,\varepsilon}\left(e^{(t+\varepsilon)X(\omega)}+e^{(t-\varepsilon)X(\omega)}\right).
\end{align*}
This is integrable because the [moment generating function](/page/Moment%20Generating%20Function) is finite at $t+\varepsilon$ and at $t-\varepsilon$:
\begin{align*}
\int_\Omega G_m(\omega)\,d\mathbb P(\omega)=C_{m,\varepsilon}\left(M_X(t+\varepsilon)+M_X(t-\varepsilon)\right)<\infty.
\end{align*}
Thus $G_m$ is the desired local dominator for every random variable of the form $|X|^m e^{sX}$ with $s$ close to $t$.
[/guided]
[/step]
[step:Differentiate once under the expectation]
Let $m \geq 0$ be an integer, and define $F_m:I_t \to \mathbb R$ by
\begin{align*}
F_m(s):=\mathbb E[X^m e^{sX}].
\end{align*}
The previous step with exponent $m$ shows that $F_m(s)$ is finite for every $s \in I_t$.
Fix $s_0 \in I_t$. Choose $\delta>0$ such that $(s_0-\delta,s_0+\delta) \subset I_t$. For $h \in \mathbb R$ with $0<|h|<\delta$, define $Q_h:\Omega \to \mathbb R$ by
\begin{align*}
Q_h(\omega):=X(\omega)^m e^{s_0X(\omega)}\frac{e^{hX(\omega)}-1}{h}.
\end{align*}
For each $\omega \in \Omega$, the one-variable [mean value theorem](/theorems/186) applied to the map $r \mapsto e^{rX(\omega)}$ on the interval between $0$ and $h$ gives a number $\theta_h(\omega)$ between $0$ and $h$ such that
\begin{align*}
Q_h(\omega)=X(\omega)^{m+1}e^{(s_0+\theta_h(\omega))X(\omega)}.
\end{align*}
Since $s_0+\theta_h(\omega)\in I_t$, the previous step with exponent $m+1$ gives
\begin{align*}
|Q_h(\omega)|\leq G_{m+1}(\omega)
\end{align*}
for every $\omega \in \Omega$. Also,
\begin{align*}
\lim_{h\to 0}Q_h(\omega)=X(\omega)^{m+1}e^{s_0X(\omega)}
\end{align*}
for every $\omega \in \Omega$. By the [dominated convergence theorem](/theorems/4) applied on the [probability space](/page/Probability%20Space) $(\Omega,\mathcal F,\mathbb P)$ with dominating function $G_{m+1}$,
\begin{align*}
F_m'(s_0)=\lim_{h\to 0}\mathbb E\left[X^m e^{s_0X}\frac{e^{hX}-1}{h}\right]=\mathbb E[X^{m+1}e^{s_0X}].
\end{align*}
Thus $F_m$ is differentiable on $I_t$ and
\begin{align*}
F_m'(s)=\mathbb E[X^{m+1}e^{sX}]
\end{align*}
for every $s \in I_t$.
[/step]
[step:Iterate the differentiation identity to obtain all derivatives]
Apply the previous step first with $m=0$. Since $F_0(s)=M_X(s)$ for $s \in I_t$, we obtain
\begin{align*}
M_X'(s)=\mathbb E[Xe^{sX}]
\end{align*}
for every $s \in I_t$.
Now proceed by induction. Suppose that, for some integer $n \geq 1$,
\begin{align*}
M_X^{(n)}(s)=\mathbb E[X^n e^{sX}]
\end{align*}
for every $s \in I_t$. The right-hand side is $F_n(s)$, and the previous step with $m=n$ gives
\begin{align*}
\frac{d}{ds}F_n(s)=\mathbb E[X^{n+1}e^{sX}]
\end{align*}
for every $s \in I_t$. Therefore
\begin{align*}
M_X^{(n+1)}(s)=\mathbb E[X^{n+1}e^{sX}]
\end{align*}
for every $s \in I_t$.
By induction, for every integer $n\geq 1$ and every $s \in I_t$,
\begin{align*}
M_X^{(n)}(s)=\mathbb E[X^n e^{sX}].
\end{align*}
Since the original point $t \in (-a,a)$ was arbitrary, $M_X$ is infinitely differentiable on $(-a,a)$ and the displayed formula holds for every $t \in (-a,a)$.
[/step]
[step:Evaluate the derivative formula at the origin]
Because $0 \in (-a,a)$ and $e^{0X}=1$, the derivative formula with $t=0$ gives, for every integer $n\geq 1$,
\begin{align*}
M_X^{(n)}(0)=\mathbb E[X^n e^{0X}]=\mathbb E[X^n].
\end{align*}
This proves the stated identity for the derivatives of the moment [generating function](/page/Generating%20Function) and, in particular, the moment formula at the origin.
[/step]
