The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
60,505 characters
Moments by Integrating the Moment-Generating Function
\title{Moments by Integrating the Moment-Generating Function\thanks{Corresponding author: Chen Tong, Email: [email removed]. We are
grateful to Christian Berg, Markus Bibinger, Simon Broda, Raymond
Kan, and Michael Wolf for helpful and valuable comments. Chen Tong
acknowledges financial support from the Youth Fund of the National
Natural Science Foundation of China (72301227), the Ministry of Education
of China, Humanities and Social Sciences Youth Fund (22YJC790117),
and the Fujian Provincial Natural Science Foundation of China (2025J08008).}{\normalsize\emph{\medskip{}
}}}
\author{\textbf{Peter Reinhard Hansen}$^{\mathsection}$$\quad$\textbf{ }and$\quad$\textbf{Chen
Tong}$^{\ddagger}$\smallskip{}
\\
{\normalsize$^{\mathsection}$}{\normalsize\emph{Department of Economics,
University of North Carolina at Chapel Hill}}\\
{\normalsize$^{\ddagger}$}{\normalsize\emph{Department of Finance,
School of Economics \& Wang Yanan Institute}}\\
{\normalsize\emph{ for Studies in Economics (WISE), Xiamen University }}}
\date{{\normalsize\emph{\today}}}
\maketitle
\begin{abstract}
We introduce a novel method for obtaining a wide variety of moments
of any random variable with a well-defined moment-generating function
(MGF). We derive new expressions for fractional moments and fractional
absolute moments, both central and non-central moments. The expressions
are relatively simple integrals that involve the MGF, but do not require
its derivatives. We label the new method CMGF because it uses a complex
extension of the MGF and can be used to obtain complex moments. We
illustrate the new method with three applications where the MGF is
available in closed-form, while the corresponding densities and the
derivatives of the MGF are either unavailable or very difficult to
obtain.
\smallskip{}
\end{abstract}
{\small\textit{{\noindent}Keywords:}}{\small{} Moments, Fractional
Moments, Moment-Generating Function.}{\small\par}
\noindent{\small\textit{{\noindent}JEL Classification:}}{\small{}
C02, C40, C65 }{\small\par}
\clearpage{}
\section{Introduction}
Moments of random variables, including conditional moments, are central
to many areas of hard and social sciences. In this paper, we introduce
a new method for obtaining moments of random variables with a well-defined
moment-generating function (MGF). The proposed method is general,
computationally efficient, and can offer solutions to problems where
existing methods are inadequate.
It is well known that the $k$-th integer moment of a random variable
$X$, with MGF, $M_{X}(s)=\mathbb{E}[e^{sX}]$, is given by $\mathbb{E}[X^{k}]=M_{X}^{(k)}(0)$,
where $M_{X}^{(k)}(s)$ is the $k$-th derivative of $M_{X}(s)$.
Non-integer moments and fractional absolute moments, $\mathbb{E}|X|^{r}$,
$r\in\mathbb{R},$ have more complicated expressions, typically integrals
that involve derivatives of the MGF or derivatives of the characteristic
function (CF).\footnote{\citet{Kawata1972} provides an expression for $\mathbb{E}|X|^{r}$
that involves several derivatives of the CF and \citet{Laue1980}
provides an expression based on fractional derivatives, see also \citet{Wolfe:1975},
\citet{SamkoKilbasMarichev:1993}, \citet{MatsuiPawlas2016}, and
\citet{Tomovski2022}. For fractional moments of non-negative variables,
similar expressions were derived in \citet{CressieDavisFolksPolicello:1981},
\citet{CressieBorkent1986}, \citet{Jones1987}, \citet{SchurgerK:2002},
and \citet{Meng:2005}.} In this paper, we derive novel integral expressions for computing
fractional moments, including fractional absolute moments, and fractional
central moments. The new expressions stand out by being applicable
to all random variables with a MGF and can be used to compute moments
of $Y=g(X)$, whenever $g(x)$ can be expressed as a linear combination
of a finite number of fractional, integer, absolute, positive-part,
and possibly centralized powers of $x$. The new expressions are particularly
useful when derivatives are prohibitively difficult to obtain, as
is often the case in dynamic models, where the conditional MGF is
defined from recursive non-linear expressions. The new method is also
useful in situations where the density is unavailable, as can be the
case for compound distributions, exponential tilting, and some L�vy
processes.
The new method is labeled CMGF because it relies on a complex extension
of the MGF, and provides expressions for complex moments, $\mathbb{E}|X|^{r}$,
$r\in\mathbb{C}$ with $\operatorname{Re}(r)>-1$. Although complex
moments are not commonly used in econometrics, this generalization
can be included with minimal adaptation.\footnote{Complex moments are commonly used other fields, including quantum
physics, number theory, and statistical mechanics.}
The new CMGF expressions involve integrals, and while these cannot
be evaluated analytically in most cases, they facilitate new ways
to compute moments numerically, which are fast and accurate. We demonstrate
this with the normal-inverse Gaussian (NIG) distribution. Existing
expressions for absolute fractional moments of the NIG distribution
are limited to $\mu$-centered moments, and the expressions are infinite
series involving modified Bessel functions of the second kind. The
CMGF integral expressions are simpler and do not involve special functions.
We also show that the CMGF method offers fast and accurate moments
in dynamic models, where we are not aware of good alternative methods.
Specifically, we use the CMGF method to compute moments of cumulative
returns in the Heston-Nandi GARCH (HNG) model, to compute moments
of realized volatilities in the Heterogeneous Autoregressive Gamma
(HARG) model, and to compute conditional moments in the Autoregressive
Poisson Model.
A key step in our proof is the identity, $\int_{-\infty}^{+\infty}\frac{e^{-itx}}{(s+it)^{r+1}}\mathrm{d}t=0$,
that holds for $s>0$, $x>0$ and $\operatorname{Re}(r)>-1$. This
identity simplifies expressions involving an inverse Laplace transform,
such that the moment-generating function emerges. The identity, see
Lemma \ref{lem:IntegralZero}, appears to be new, which can explain
that these simple and general integral expressions were not discovered
earlier. We are not the first to derive integral expressions that
involve the MGF. \citet{Meng:2005} derived integral expressions of
$\mathbb{E}[X^{a}/Y^{b}]$, involving their joint MGF $M_{X,Y}(s_{1},s_{2})$
and its derivatives. The possibility of combining characteristic and
moment-generating functions, $\mathbb{E}[e^{sY+itX}]$, is mentioned
in the epilogue of \citet{Meng:2005}. For $X=Y$ this is an extension
of the MGF (or the CF) to the complex plane, $M_{X}(s+it)$, which
is at the core of our expression.
An auxiliary result that arises from applying the integral expression
to standard normal distribution, is a holomorphic integral representation
of Gamma function and its reciprocal, which do not rely on analytic
continuation, see \citet{HansenTong:2025ReciprocalGamma}.
The remainder of this paper is organized as follows. We present the
new moment expressions in Section 2. Theorem \ref{thm:factional-abs-moments}
has the expression for absolute fractional moments for general random
variables, Theorem \ref{thm:pos-part-moments} presents positive-part
moments, and Theorem \ref{thm:integer-moments} presents a new expression
for integer moments. Corollary \ref{cor:factional-moments-positive}
has expressions for fractional moments of nonnegative (positive) random
variables and Corollary \ref{cor:TailMoments} has expression for
(odd integers) negative tail moments. We apply CMGF to the normal-inverse
Gaussian distribution. Section 3 presents three applications of the
new expressions in dynamic models, and Section 4 concludes.
\section{A New Method for Computing Moments (CMGF)}
We use the following notation: $\mathbb{N}$ denotes the positive
integers and $\mathbb{N}_{0}=\{0\}\cup\mathbb{N}$. We use $z=s+it\in\mathbb{C}$,
where $s=\operatorname{Re}(z)>0$. We use $\mathbb{E}[X^{k}]$ to
denote integer moments, $k\in\mathbb{N}_{0}$, whereas $\mathbb{E}[X^{r}]$,
denotes general moments including fractional and complex moments,
$r\in\mathbb{C}$. The MGF is typically introduced as a real function
with a real argument, such that both $s$ and $M_{X}(s)$ take values
in $\mathbb{R}.$ In this paper, we establish our results by means
of complex arguments, $z\in\mathbb{C},$ such that $M_{X}(z)$ also
takes values in $\mathbb{C}$.
\textbf{Interpretation and scope:} For complex powers $z^{-(r+1)}$,
we fix the principal branch of the complex logarithm (branch cut on
$(-\infty,0]$) and interpret integrals as oscillatory (i.e., symmetric
Cauchy principal values) whenever $-1<\operatorname{Re}(r)\leq0$.
Assumptions on one-sided exponential moments of $X$ (i.e., $\mathbb{E}[e^{\pm sX}]<\infty$
for some $s>0$) are stated with each theorem. When $\operatorname{Re}(r)\in(-1,0)$,
we additionally assume $\mathbb{E},\lvert X-\xi\rvert^{r}<\infty$
to justify exchanging expectation and integration.
To gain some intuition for the complex extension of the MGF, we observe
that
\[
\psi(s,t)\equiv M_{X}(s+it)=\mathbb{E}[e^{(s+it)X}]\in\mathbb{C},
\]
nests both the standard moment-generating function, $\psi(s,0)=M_{X}(s)$,
and the characteristic function, $\psi(0,t)=\varphi_{X}(t)\equiv\mathbb{E}[e^{itX}]$,
as special cases. The complex-valued MGF, $M_{X}(z)$, is the bilateral
Laplace--Stieltjes transform (apart from the sign of $z$), and for
continuous random variables it is similar to the Laplace transform.\footnote{The Laplace transform of the density function, $f_{X}(x)$, is $\int_{0}^{\infty}e^{-zx}f_{X}(x)\mathrm{d}x$,
and it differs from $M_{X}(z)$ in terms of the domain of integration
and the sign of the exponent. Moreover, $M_{X}(z)$ is applicable
to discrete distributions and mixtures of continuous and discrete
distributions. }
The following result is key to several simplifications.
\begin{lem}
\label{lem:IntegralZero}
\begin{equation}
\int_{-\infty}^{+\infty}\frac{e^{-itx}}{(s+it)^{r+1}}\mathrm{d}t=0,\qquad\text{for}\quad x>0,\ s>0,\ \operatorname{Re}(r)>-1.\label{eq:Pos-s-identity}
\end{equation}
For $x=0$ the identity holds for $\operatorname{Re}(r)>0$.\footnote{Obviously, for $x<0$ the identity is $\int_{-\infty}^{+\infty}\frac{e^{itx}}{(s+it)^{r+1}}\mathrm{d}t=0$
for $s>0$ and $\operatorname{Re}(r)>-1$.}
\end{lem}
The identity in (\ref{eq:Pos-s-identity}) is obtained by a contour
argument that is valid for $\operatorname{Re}(r)>-1$. Note that $f\in L^{1}(\mathbb{R})$
if and only if $\operatorname{Re}(r)>0$, because
\[
\left|f(t)\right|=\left|\frac{e^{-itx}}{(s+it)^{r+1}}\right|=\frac{1}{|s+it|^{\operatorname{Re}(r)+1}}=\frac{1}{(s^{2}+t^{2})^{(\operatorname{Re}(r)+1)/2}}.
\]
For $-1<\operatorname{Re}(r)\leq0$ the integral is not absolutely
convergent; in this range we interpret the formula (\ref{eq:Pos-s-identity})
as the oscillatory improper integral,
\[
\lim_{T\to\infty}\int_{-T}^{T}\frac{e^{-itx}}{(s+it)^{r+1}}\mathrm{d}t,
\]
which exists and equals $0$ for $x>0$ by closing the contour in
the lower half-plane ($e^{-itx}$ decays there) and observing that
$s+it\neq0$ for $t\in\mathbb{R}$ when $s>0$.
For the special case $x=0$ and $r=0$, it is easy to establish the
well-known identity
\begin{equation}
\int_{-\infty}^{+\infty}\frac{1}{s+it}\mathrm{d}t=\pi\operatorname{sgn}(s),\qquad\text{for all }s\in\mathbb{R},\label{eq:x0r0sp}
\end{equation}
where $\operatorname{sgn}(s)\equiv\frac{s}{|s|}1_{\{s\neq0\}}\in\{-1,0,1\}$.
In this setting where $s>0$, we simply have $\int_{-\infty}^{+\infty}\frac{1}{s+it}\mathrm{d}t=\pi$,
which we will use in the first formula of following Lemma.
\begin{lem}
\label{lem:LemmaNew}Let $x\in\mathbb{R}$ and $z=s+it\in\mathbb{C}$,
where $s>0$ is an arbitrary constant. Then
\[
|x|^{r}=\frac{\Gamma(r+1)}{2\pi}\int_{-\infty}^{+\infty}\frac{e^{zx}+e^{-zx}}{z^{r+1}}\mathrm{d}t,\qquad\text{for}\quad\operatorname{Re}(r)>-1,\quad x\neq0,
\]
and this identity also holds for $x=0$ if $\operatorname{Re}(r)>0$
or $r=0$.\footnote{When $x^{r}$ is to be evaluated for $x=r=0$, we follow the standard
convention $0^{0}=1$.}
For $x>0$ we have
\[
x^{r}=\frac{\Gamma(r+1)}{2\pi}\int_{-\infty}^{+\infty}\frac{e^{zx}}{z^{r+1}}\mathrm{d}t,\qquad\text{for}\quad\operatorname{Re}(r)>-1,\quad x>0,
\]
and this identity also holds for $x=0$ if $\operatorname{Re}(r)>0$.
Similarly,
\begin{eqnarray*}
x^{k} & = & \frac{k!}{2\pi}\int_{-\infty}^{+\infty}\frac{e^{zx}+\left(-1\right)^{k}e^{-zx}}{z^{k+1}}\mathrm{d}t,\qquad k\in\mathbb{N}_{0},\quad x\in\mathbb{R}.
\end{eqnarray*}
\end{lem}
Here $x=0$, arises as a special case in the first two formulas of
Lemma \ref{lem:LemmaNew}, and this requiring results to include a
qualification for distributions with positive mass at zero.
Several results will now follows by taking expectation to both sides
of the identities in Lemma \ref{lem:LemmaNew}, and justify interchanging
expectation and integral. For $\operatorname{Re}(r)>0$ this is straight
forward by means of Fubini/Tonelli theorems. However, as detailed
above, the integrand is not absolute integrable for $\operatorname{Re}(r)\leq0$,
and interchanging expectation and integral requires a different line
of arguments, and the improper integrals in the form $\int_{-\infty}^{+\infty}[\cdot]\mathrm{d}t$
is to be interpreted as $\lim_{T\rightarrow\infty}\int_{-T}^{+T}[\cdot]\mathrm{d}t$
for $\operatorname{Re}(r)\in(-1,0]$.
Our first main result concerns fractional absolute moments.
\begin{thm}[Absolute moments]
\label{thm:factional-abs-moments}Suppose that $\mathbb{E}[e^{\pm sX}]<\infty$
for some $s>0$. For $\operatorname{Re}(r)\in(-1,0)$, assume additionally
that $\mathbb{E}|X-\xi|^{r}<\infty$. If $\Pr(X=\xi)=0$ then
\begin{equation}
\mathbb{E}|X-\xi|^{r}=\frac{\Gamma(r+1)}{\pi}\int_{0}^{\infty}\operatorname{Re}\left[\frac{e^{-\xi z}M_{X}(z)+e^{\xi z}M_{X}(-z)}{z^{r+1}}\right]\mathrm{d}t,\quad\text{for }r>-1,\label{eq:moment-r-abs-real}
\end{equation}
and for complex moments, $r\in\mathbb{C}$, we have
\begin{equation}
\mathbb{E}|X-\xi|^{r}=\frac{\Gamma(r+1)}{2\pi}\int_{-\infty}^{+\infty}\frac{e^{-\xi z}M_{X}(z)+e^{\xi z}M_{X}(-z)}{z^{r+1}}\mathrm{d}t,\quad\operatorname{Re}(r)>-1.\label{eq:moment-r-abs}
\end{equation}
Moreover, if $\Pr(X=\xi)>0$, the identities hold for $r\geq0$ and
$\ensuremath{\ensuremath{\operatorname{Re}(r)>0}}$ respectively.
\end{thm}
The finite moment requirement, $\mathbb{E}|X-\xi|^{r}<\infty$, is
implied by $\mathbb{E}[e^{\pm sX}]<\infty$ for $\operatorname{Re}(r)>0$,
but is not guaranteed for $\operatorname{Re}(r)<0$. The results in
Theorem \ref{thm:factional-abs-moments} include regular absolute
moments by setting $\xi=0$ and central moments by setting $\xi=\mathbb{E}X$.
Although analytical integration may often be impractical, (\ref{eq:moment-r-abs})
offers a novel method to evaluate moments numerically. For real moments,
$r\in\mathbb{R}$, it is our experience that it is computationally
advantageous to use (\ref{eq:moment-r-abs-real}) rather than (\ref{eq:moment-r-abs}).
\begin{rem*}
$\mathbb{E}[e^{\pm sX}]<\infty$ for some $s>0$ is equivalent to
assuming the MGF is well defined in a neighborhood of zero. In next
Theorem we only assume $\mathbb{E}[e^{sX}]<\infty$ for some $s>0$,
which does not guarantee the MGF is well-defined in a neighborhood
about zero.
\end{rem*}
Interestingly, Theorem \ref{thm:factional-abs-moments} can be used
to obtain new identities, by equating (\ref{eq:moment-r-abs}) with
an existing expression for absolute moments. For instance, the absolute
moment of a standard normal random variable, $X\sim N(0,1)$, is $\mathbb{E}|X|^{r}=\frac{\Gamma(\frac{r+1}{2})}{\sqrt{\pi}}2^{r/2}$,
for $r>-1$, and equating this with (\ref{eq:moment-r-abs}) yields
the following expression for the \emph{reciprocal gamma function}
\begin{equation}
\frac{1}{\Gamma(\frac{r}{2}+1)}=\frac{2^{r/2}}{\pi}\int_{-\infty}^{+\infty}\frac{e^{z^{2}/2}}{z^{r+1}}\mathrm{d}t,\qquad r>-1,\label{eq:ReciprocalGamma}
\end{equation}
for $z=s+it$, where $s>0$ is an arbitrary positive constant. From
Theorem \ref{thm:factional-abs-moments} we can only infer that (\ref{eq:ReciprocalGamma})
holds for real $r>-1$, but the identity actually holds for any $r\in\mathbb{C}$,
such that we obtain a new holomorphic integral representation of Gamma
function and its reciprocal that does not rely on analytic continuation,
see \citet{HansenTong:2025ReciprocalGamma}.
Below we will derive an expression for positive-part moments, where
we use the convention $x_{+}^{r}=x^{r}1_{\{x>0\}}$. This enables
us to accommodate negative powers, because $0_{+}^{r}=0$ for all
$r$.\footnote{Note that with this convention we have ($r=0$) we have $0=0_{+}^{0}$,
unlike $0^{0}=1$).} Importantly, for $\operatorname{Re}(r)\in(-1,0]$, the integral $\int_{-\infty}^{+\infty}\frac{e^{-\xi z}M_{X}(z)}{z^{r+1}}\mathrm{d}t$,
is to be interpreted as the symmetric Cauchy principal value, $\lim_{T\rightarrow\infty}\int_{-T}^{+T}\frac{e^{-\xi z}M_{X}(z)}{z^{r+1}}\mathrm{d}t$.
\begin{thm}[Positive-part moments]
\label{thm:pos-part-moments}Suppose that $\mathbb{E}[e^{sX}]<\infty$
for some $s>0$. For $\operatorname{Re}(r)\in(-1,0)$, assume additionally
that $\mathbb{E}|X-\xi|^{r}<\infty$. Then
\begin{equation}
\mathbb{E}[(X-\xi)_{+}^{r}]=\frac{\Gamma(r+1)}{\pi}\int_{0}^{\infty}{\rm Re}\left[\frac{e^{-\xi z}M_{X}(z)}{z^{r+1}}\right]\mathrm{d}t,\quad\text{for }r>0,\label{eq:pos-part-moment}
\end{equation}
and the the identity also holds for $r\in(-1,0]$ if $\Pr(X=\xi)=0$.
For complex moments, $r\in\mathbb{C}$, we have
\begin{equation}
\mathbb{E}[(X-\xi)_{+}^{r}]=\frac{\Gamma(r+1)}{2\pi}\int_{-\infty}^{+\infty}\frac{e^{-\xi z}M_{X}(z)}{z^{r+1}}\mathrm{d}t,\quad\text{for }\operatorname{Re}(r)>-1.\label{eq:pos-var-moment-real}
\end{equation}
\end{thm}
Note that $\mathbb{E}[e^{sX}]<\infty$ for some $s>0$ ensures that
$\mathbb{E}(X-\xi)_{+}^{r}<\infty$ for ${\rm Re}\left(r\right)>0$
but not for $\operatorname{Re}(r)\in\left(-1,0\right)$. Therefore,
for $\operatorname{Re}(r)\in\left(-1,0\right)$ where we rely on the
dominated convergence theorem, the condition $\mathbb{E}|X-\xi|^{r}<\infty$
is required for interchanging expectation and integral.
\begin{rem}
Here we assume $\mathbb{E}[e^{sX}]<\infty$ for some $s>0$, which
does not guarantee the existence of the MGF in an open neighborhood
of zero. So, $\mathbb{E}|X-\xi|^{r}<\infty$ is not guaranteed by
$\mathbb{E}[e^{sX}]<\infty$, as illustrated with the example below.
Note that $\mathbb{E}[e^{sX}]<\infty$ for some $s>0$ does ensure
$\mathbb{E}(X-\xi)_{+}^{r}<\infty$ for $\operatorname{Re}(r)>0$.
\end{rem}
\begin{example}[One-sided MGF and divergent moments.]
Let $Y\sim\mathrm{Exp}(1)$, $Z\sim\mathrm{Pareto}(\alpha)$ with
$\mathbb{P}(Z>z)=z^{-\alpha}$ for $z\geq1$, and choose $0<\alpha<1$.
Define
\[
X=\begin{cases}
+Y, & \text{with prob. }1/2,\\
-Z, & \text{with prob. }1/2,
\end{cases}
\]
with $Y$, $Z$ independent. Then, for $0<s<1$, $\mathbb{E}[e^{sX}]=\tfrac{1}{2}\mathbb{E}[e^{sY}]+\tfrac{1}{2}\mathbb{E}[e^{-sZ}]=\tfrac{1}{2}(1-s)^{-1}+\tfrac{1}{2}\mathbb{E}[e^{-sZ}]<\infty$.
However, for $\tilde{s}<0$ we have $\mathbb{E}[e^{\tilde{s}X}]\ge\tfrac{1}{2}\mathbb{E}[e^{|\tilde{s}|Z}]=\tfrac{1}{2}\int_{1}^{\infty}e^{|\tilde{s}|z}\alpha z^{-\alpha-1}\mathrm{d}z=\infty$.
Hence the MGF is not finite on any open neighborhood of 0. For $r\in(-1,\alpha)$
we have $\mathbb{E}|X|^{r}=\tfrac{1}{2}\Gamma(r+1)+\tfrac{1}{2}\frac{\alpha}{\alpha-r}$,
while for $r\geq\alpha$, $\mathbb{E}[|X|^{r}]\ge\tfrac{1}{2}\mathbb{E}[Z^{r}]=\infty$.
We can nevertheless apply Theorem \ref{thm:pos-part-moments} for
$r>-1$ because $\mathbb{E}(X)_{+}^{r}=\tfrac{1}{2}\Gamma(r+1)$.
\end{example}
The expression in Theorem \ref{thm:pos-part-moments} is identical
to one in \citet[theorem 1]{Pinelis:2011}, albeit we extend the range
of moments to include some negative moments, $r\in(-1,0)$, and complex
moments with $\operatorname{Re}(r)>-1$. The special case $r=1$ has
received considerable attention because $\mathbb{E}[(X-\xi)_{+}]$
plays a central role in option pricing, with $\xi$ representing the
strike price. For this special case, $r=1$, (\ref{eq:pos-var-moment-real})
coincides with \citet[lemma 7.1]{KimRachevBianchiFabozzi:2010} and
(\ref{eq:pos-part-moment}) coincides with \citet[eq. 3.1]{HuangOosterlee:2011}.\footnote{They derived the result under slightly stronger assumptions. In \citet{KimRachevBianchiFabozzi:2010}
it is assumed that the distribution is infinitely divisible and continuous,
while \citet[eq. 3.1]{HuangOosterlee:2011} assume continuity of the
distribution. Neither of these assumptions is required here.}
Note that the expression remains invariant to the value of $s>0$,
and numerical integration is, in our experience, insensitive to the
choice of $s$, so long as $\mathbb{E}[e^{sX}]<\infty$. For the positive-part
moments, we do not require $\mathbb{E}[e^{uX}]<\infty$ for any $u<0$.
For non-negative random variables the following is a simple implication
of Theorem \ref{thm:pos-part-moments} and Lemma \ref{lem:LemmaNew}.
\begin{cor}[Moments of nonnegative random variables]
\label{cor:factional-moments-positive}Suppose that $\mathbb{E}[e^{sX}]<\infty$
for some $s>0$ and for $\operatorname{Re}(r)\in(-1,0)$, assume additionally
that $\mathbb{E}|X-\xi|^{r}<\infty$. (Strictly positive). If $\Pr(X>\xi)=1$,
then
\begin{equation}
\mathbb{E}[(X-\xi)^{r}]=\frac{\Gamma(r+1)}{\pi}\int_{0}^{\infty}{\rm Re}\left[\frac{e^{-\xi z}M_{X}(z)}{z^{r+1}}\right]\mathrm{d}t\quad\text{for }r>-1.\label{eq:moment-pos-var-real}
\end{equation}
(Nonnegative). If $\Pr(X\geq\xi)=1$, then the identity holds for
$r>0$, and for complex moments, $r\in\mathbb{C}$, we have the identity
\begin{equation}
\mathbb{E}[(X-\xi)^{r}]=\frac{\Gamma(r+1)}{2\pi}\int_{-\infty}^{+\infty}\frac{e^{-\xi z}M_{X}(z)}{z^{r+1}}\mathrm{d}t,\quad\text{for }\operatorname{Re}(r)>-1.\label{eq:moment-pos-var}
\end{equation}
\end{cor}
The condition $\mathbb{E}[(X-\xi)^{r}]<\infty$ is redundant for $\operatorname{Re}(r)>0$,
because it is implied by $\mathbb{E}[e^{sX}]<\infty$, but for negative
moments this additional conditions is needed.
As an illustration, we apply (\ref{eq:moment-pos-var}) to an exponentially
distributed random variable. This example serves as an illustration
because integer moments, $k\in\mathbb{N}$, are easy to obtain from
the derivatives, $M_{X}^{(k)}(s)$, and other moments can be obtained
by evaluating $\int_{0}^{\infty}x^{r}\lambda e^{\lambda x}\mathrm{d}x$
directly.
\begin{example}
For an exponentially distributed random variable, with parameter $\lambda>0$,
$X\sim\operatorname{Exp}(\lambda)$, we have $M_{X}(z)=\frac{\lambda}{\lambda-z}$,
with region of convergence $\operatorname{Re}(z)<\lambda$. So, the
integral in (\ref{eq:moment-pos-var-real}) is $\int_{0}^{\infty}{\rm Re}\left[\frac{\lambda}{\lambda-z}\frac{1}{z^{r+1}}\right]\mathrm{d}t$.
For an integer moment, $k\in\mathbb{N}$, this integral equals
\begin{align*}
\int_{0}^{\infty}{\rm Re}\left[\lambda^{-k}\tfrac{1}{\lambda-z}+\lambda^{-k}\tfrac{1}{z}+\sum_{j=2}^{k+1}\lambda^{-k+j-1}\tfrac{1}{z^{j}}\right]\mathrm{d}t & =\lambda^{-k}\int_{0}^{\infty}{\rm Re}\left[\tfrac{1}{\lambda-z}+\tfrac{1}{z}\right]\mathrm{d}t\\
& =\lambda^{-k}\int_{0}^{\infty}{\rm Re}\left[\tfrac{\lambda-z^{\ast}}{(\lambda-z^{\ast})(\lambda-z)}+\tfrac{z^{\ast}}{z^{\ast}z}\right]\mathrm{d}t,
\end{align*}
as $\int_{0}^{\infty}{\rm Re}\left[z^{-j}\right]\mathrm{d}t=0$ for
$j>1$. Substituting $z=s+it$ in the last expression gives us
\begin{eqnarray*}
\lambda^{-k}\int_{0}^{\infty}\operatorname{Re}\left[\tfrac{\lambda-s+it}{(\lambda-s)^{2}-(it)^{2}}+\tfrac{s-it}{s^{2}-(it)^{2}}\right]\mathrm{d}t & = & \lambda^{-k}\int_{0}^{\infty}\left[\tfrac{\lambda-s}{(\lambda-s)^{2}+t^{2}}+\tfrac{s}{s^{2}+t^{2}}\right]\mathrm{d}t\\
& = & \lambda^{-k}\left[\left.\arctan(\tfrac{t}{\lambda-s})+\ensuremath{\arctan(\tfrac{t}{s})}\right|_{0}^{\infty}\right]=\pi\lambda^{-k}.
\end{eqnarray*}
Thus that $\mathbb{E}[X^{k}]=\lambda^{-k}k!$ as expected. Non-integer
moments, $r\notin\mathbb{N}$, can be derived using contour integrals
and the Cauchy Integral Theorem.
\end{example}
For random variables whose support may include negative numbers, we
have the following result for integer moments.
\begin{thm}
\label{thm:integer-moments}Suppose $\mathbb{E}[e^{\pm sX}]<\infty$
for some $s>0$. Then
\begin{eqnarray}
\mathbb{E}[(X-\xi)^{k}] & = & \frac{k!}{\pi}\int_{0}^{\infty}\operatorname{Re}\ensuremath{\left[\frac{e^{-\xi z}M_{X}(z)+(-1)^{k}e^{\xi z}M_{X}(-z)}{z^{k+1}}\right]\mathrm{d}t,\quad\text{for }k\in\mathbb{N}_{0}.}\label{eq:moment-k}
\end{eqnarray}
\end{thm}
Because $k$ is real, we have expressed the integral in the form,
$\frac{k!}{\pi}\int_{0}^{\infty}\operatorname{Re}[\cdot]\mathrm{d}t$,
which tends to be simpler to evaluate numerically than the equivalent
expression using the form, $\frac{k!}{2\pi}\int_{-\infty}^{+\infty}[\cdot]\mathrm{d}t$.
\begin{rem*}
In some situations, Theorems \ref{thm:factional-abs-moments}, \ref{thm:pos-part-moments},
and \ref{thm:integer-moments}, provide different expressions for
the same moments. For instance, (\ref{eq:moment-r-abs}) and (\ref{eq:moment-pos-var})
must agree for a non-negative random variable, because $\mathbb{E}|X|^{r}=\mathbb{E}[X^{r}]$,
and for integer moments, these expressions must also agree with (\ref{eq:moment-k}).
This is indeed the case, because the additional term in (\ref{eq:moment-r-abs})
and the additional term in (\ref{eq:moment-k}) are both zero for
non-negative variables. This is an implication of Lemma \ref{lem:IntegralZero},
which establishes that $\ensuremath{\int_{-\infty}^{\infty}\frac{e^{-\left(s+it\right)x}}{\left(s+it\right)^{r+1}}\mathrm{d}t=0}$
for all $x>0$ and all $s>0$.
\end{rem*}
Negative tail moments, which is part of the definition of \emph{expected
shortfall} (see below), have the following expression, where we can
combine Theorems \ref{thm:pos-part-moments} and \ref{thm:integer-moments}
and use the notation $x_{-}=\min(x,0)$.
\begin{cor}[Negative tail moments]
\label{cor:TailMoments}Suppose $\mathbb{E}[e^{-sX}]<\infty$ for
some $s>0$ and $\mathbb{E}|X-\xi|^{k}<\infty$. Then
\[
\mathbb{E}[(X-\xi)_{-}^{k}]=-\frac{k!}{\pi}\int_{0}^{\infty}\ensuremath{{\rm Re}\left[\frac{e^{\xi z}M_{X}(-z)}{z^{k+1}}\right]\mathrm{d}t,\quad\text{for }}k\text{ odd}.
\]
\end{cor}
We can also compute the cumulative distribution function (cdf), quantiles,
and expected shortfall using the new moment expressions. For a continuous
random variable, $X$, its cdf is explicitly given by, $F_{X}(x)=\ensuremath{\frac{1}{2}-\frac{1}{\pi}\int_{0}^{\infty}{\rm Re}\left[\frac{M_{X}(it)e^{-itx}}{it}\right]\mathrm{d}t},$
which is the classical Gil-Pelaez inversion formula that provides
a direct way to recover probability distributions from their transforms.
Here we have used that $M_{X}(it)$ is the characteristic function
for $X$. If unique, the $\alpha$-quantile is now given by $\xi_{\alpha}\equiv F_{X}^{-1}(\alpha)$,
and we can use Corollary \ref{cor:TailMoments} with $k=1$ to compute
expected shortfall,
\[
\operatorname{ES}_{\alpha}(X)=-\frac{1}{\alpha}\mathbb{E}[X1_{\{X<\xi_{\alpha}\}}]=-\frac{1}{\alpha}\mathbb{E}[(X-\xi_{\alpha})1_{\{X<\xi_{\alpha}\}}]-\xi_{\alpha}.
\]
\subsection{Example: Moments of NIG Distribution}
To illustrate the new method, we use Theorem \ref{thm:factional-abs-moments}
to compute fractional moments of the normal-inverse Gaussian (NIG)
distribution. The NIG distribution was introduced in \citet{Barndorff-Nielsen:1978}
and has four parameters, $\mu$ (location), $\delta$ (scale), $\alpha$
(tail heaviness), and $\beta$ (asymmetry). It has density
\[
f(x)=\tfrac{\alpha\delta}{\pi\sqrt{\delta^{2}+(x-\mu)^{2}}}K_{1}\left(\alpha\sqrt{\delta^{2}+(x-\mu)^{2}}\right)e^{\delta\gamma+\beta(x-\mu)},\qquad x\in\mathbb{R}
\]
where $\gamma=\sqrt{\alpha^{2}-\beta^{2}}$ and $K_{1}(\cdot)$ is
the modified Bessel function of the second kind, and the MGF is given
by $M_{X}(z)=\exp\left(\mu z+\delta\left[\gamma-\sqrt{\alpha^{2}-(\beta-z)^{2}}\right]\right),$
for $z=s+it$ with $\beta-\alpha<s<\beta+\alpha$.
Existing expressions for fractional moments of NIG distributions are
rather complicated. Let $X\sim\mathrm{NIG}(\alpha,\beta,\mu,\delta)$
and $Y\sim\mathrm{NIG}(\alpha,0,\mu,\delta)$, then \citet{Barndorff-NielsenStelzer:2005}
showed that
\[
\mathbb{E}|X|^{r}=e^{\delta(\gamma-\alpha)}\sum_{k=0}^{\infty}\frac{\beta^{k}}{k!}\mathbb{E}\left(|Y|^{r}(Y-\mu)^{k}\right),\qquad r>0,
\]
which expresses $\mathbb{E}|X|^{r}$ as an infinite sum of moments
involving the symmetric NIG, $Y$. For $\mu$-centered moments \citet[corollary 3]{Barndorff-NielsenStelzer:2005}
obtained the following explicit formula,
\begin{equation}
\mathbb{E}(|X-\mu|^{r})=\frac{\alpha}{\pi}\left(\frac{2\delta}{\alpha}\right)^{\tfrac{r+1}{2}}e^{\delta\gamma}\sum_{k=0}^{\infty}\frac{2^{k}\Gamma(k+\tfrac{r+1}{2})}{(2k)!}\left(\frac{\delta\beta^{2}}{\alpha}\right)^{k}K_{k+\tfrac{r-1}{2}}(\delta\alpha),\quad r>0.\label{eq:InfSum-NIG-expression}
\end{equation}
The corresponding CMGF integral is
\[
\mathbb{E}|X-\mu|^{r}=\frac{\Gamma(r+1)}{2\pi}e^{\delta\gamma}\int_{-\infty}^{+\infty}\frac{e^{-\delta\sqrt{\alpha^{2}-(\beta-z)^{2}}}+e^{-\delta\sqrt{\alpha^{2}-(\beta+z)^{2}}}}{z^{r+1}}\mathrm{d}t,\qquad r>-1,
\]
and Theorem \ref{thm:factional-abs-moments} yields similar expressions
for any other recentering ($\xi\neq\mu$), which is not possible with
(\ref{eq:InfSum-NIG-expression}). Evaluating the CMGF integral is
also substantially faster and more accurate than both simulation-based
methods and direct integration using the NIG density.
To illustrate this, we consider two standardized NIG distribution
(zero mean and unit variance), which can be characterized by the two
parameters, $(\xi,\chi)$, for $0\leq|\chi|<\xi<1$, where $\mu=-\tfrac{\chi}{\xi}\zeta$,
$\delta=\xi^{2}\sqrt{\xi^{2}-\chi^{2}}\zeta$, $\alpha=\xi\zeta$,
and $\beta=\chi\zeta$, with $\zeta=\sqrt{1-\xi^{2}}/(\xi^{2}-\chi^{2})$,
see \citet{BNBJS1985nig}. Specifically we consider $(\xi,\chi)=(1/2,-1/3)$
and $(\xi,\chi)=(1/8,-1/16)$, the corresponding densities are shown
in the left panel of Figure \ref{fig:NIG_dens_MGF}. The suitable
range for $s$ is $0<s<\alpha+\beta$, where the upper limit for these
distributions is $3\sqrt{3}/5\approx1.04$ and $2\sqrt{63}/3\approx5.29$,
respectively. So, we can set $s=1$ and the mapping, and $t\mapsto\operatorname{MGF}(1+it)$
is shown for the two distributions in the right panel -- solid lines
for the real part and dashed lines for the imaginary part.
\begin{figure}[tbh]
\centering{}\includegraphics[width=1\textwidth]{Figures/NIG_dens_cMGF}\caption{Densities of two standardized NIG{\small{} distributions (left panel)
and the $t\protect\mapsto$MGF($1+it)$ (right panel). The two distributions
are for $(\xi,\chi)=(1/2,-1/3)$ (blue lines) and }$(\xi,\chi)=(1/8,-1/16)${\small{}
(red lines).\label{fig:NIG_dens_MGF}}}
\end{figure}
We compute the absolute moments for $-0.85\leq r\leq4.2$ for the
two distributions using Theorem \ref{thm:factional-abs-moments} and
the moments are shown in the upper panel of Figure \ref{fig:MomentsNIG}.
The solid black dot is the second (absolute) moment, which is known
to be one for these standardized distributions, and the solid blue
and red dots represent the fourth (absolute) moment, which are known
to equal $3(1+4\chi^{2})/(1-\xi^{2})$. The values provided by the
CMGF method are very accurate. The first 12 digits of the fourth moments
are correct using a simple implementation in Julia.
\begin{figure}[tbh]
\begin{centering}
\includegraphics[width=1\textwidth]{Figures/NIGmoments}\caption{Absolute moments of the standardized NIG{\small{} distribution with
$(\xi,\chi)=(1/2,-1/3)$ (blue lines) and }$(\xi,\chi)=(1/8,-1/16)${\small{}
(red lines). These are shown in the left panels, where dots represent
the known moments for $r=2$ and $r=4$. The middle panel is a snippet
of the left panel where we have added x-crosses that represents 100
simulation-based estimates of $\mathbb{E}|X|^{r}$, for $r\in\{3.82,3.84,\ldots,4.18\}$.
Each estimate is based on $N=$1,000,000 random draws of the NIG distribution.
The right panel shows how many decimal places are accurate with a
one standard deviation simulation error.\label{fig:MomentsNIG}}}
\par\end{centering}
\end{figure}
Additionally, we estimate the absolute moments using $N=1,000,000$
independent draws from the two NIG distributions. The simulation-based
moment is given by $\hat{\mu}_{N}^{r}=\frac{1}{N}\sum_{j=1}^{N}|X_{i}|^{r}$
and we estimate the $\sigma^{2}(r)=\operatorname{var}(|X_{j}|^{r})$
with 100 million draws. The accuracy of a simulated moment is quantified
by its standard error. Assuming an error equal to \emph{one} standard
deviation, then the number of accurate decimal places is given by
$-\log_{10}(\sigma(r)/\sqrt{N})$.\footnote{For instance, a standard error equal to $0.00099$ and $N=10^{6}$
will map to $-\log_{10}(0.99\cdot10^{-3})\simeq3.004$.} We have plotted this number for $N=10^{6}$ in the lower left panel
of Figure \ref{fig:MomentsNIG}. The simulation-based moments based
on $N=1,000,000$ draws from the distribution are far less accurate.
At best, the first few decimal places are accurate, and the accuracy
deteriorates rapidly as $r$ increases. In the lower right panels
we have shown the moments based on the CMGF method for $3.8\leq r\leq4.2$
represented by the solid blue line. Each $\times$-cross is a simulation-based
moment estimate based on $N=1,000,000$ random draws.
\begin{table}[th]
\caption{CMGF Computation Time of $\mathbb{E}|X|^{r}$ for $X\sim\operatorname{NIG}$}
\begin{centering}
\vspace{0.2cm}
\begin{small}
\begin{tabularx}{\textwidth}{YYYYY}
\toprule
\midrule
$r$ & CMGF & $M_X^{(k)}$ & Integrate with density & Simulations $(N=10^6)$ \\
\midrule
\\[-0.2cm]
-0.5 & 28.3 & & 712.9 & 397,666 \\
0.5 & 29.2 & & 244.0 & 395,837 \\
1.0 & 20.3 & & 199.1 & 299,232 \\
1.5 & 29.3 & & 175.0 & 399,173 \\
2.0 & 20.4 & 55.7 & 108.3 & 298,635 \\
2.5 & 24.5 & & 149.3 & 399,182 \\
3.0 & 17.0 & & 122.5 & 291,924 \\
3.5 & 24.5 & & 114.7 & 402,552 \\
4.0 & 17.6 & 104.9 & 100.3 & 318,081 \\
\\[-0.2cm]
\midrule
\bottomrule
\end{tabularx}
\end{small}
\par\end{centering}
{\small Note: Computation time in microseconds (\textgreek{\textmu}s)
for evaluating $\mathbb{E}|X|^{r}$ using four methods: the CMGF method
of Theorem \ref{thm:factional-abs-moments}, the $k$-th derivative
of the MGF (for $r=2$ and $r=4$), by numerical integration of $\int_{-\infty}^{\infty}|x|^{r}f(x)\mathrm{d}x$,
and by simulations, where we generate $X_{1},\ldots,X_{N}$ independent
and identically NIG distributed, with $N=$1,000,000 and take the
average of $|X_{i}|^{r}$. Both CMGF and the density-based method
use numerical integration with the same tolerance threshold. Computation
times are evaluated with BenchmarkTools.jl for Julia, see \citet{BenchmarkTools.jl-2016}.
Computations were done with Julia v1.11.0, see \citet{Julia2017},
on a MacBook Pro M1 Max with 32 GM memory. \label{tab:NIGmoments}}{\small\par}
\end{table}
The CMGF method is not only far more accurate in this example, it
is also much faster than the simulation-based approach. This is shown
in Table \ref{tab:NIGmoments} where we report the computation time
obtained with a Julia implementation. In this case the density is
readily available, and we also compute the moment by evaluating $\int_{-\infty}^{\infty}|x|^{r}f(x)\mathrm{d}x$
using the same numerical integration method used to evaluate the integral
in Theorem \ref{thm:factional-abs-moments}. The CMGF method is between
5 and 25 times faster than the standard approach, and the CMGF method
is more than 10,000 times faster than averaging (and generating) $N=10^{6}$
independent draws from the NIG distribution. We have also compared
the methods using a standard normal distribution, for which the true
moment is known to be, $\mathbb{E}|X|^{r}=\frac{\Gamma(\frac{r+1}{2})}{\sqrt{\pi}}2^{r/2}$,
for all $r>-1$. For the Gaussian case the CMGF is also thousands
of times faster than simulation methods using 1 million draws and
far more accurate. For $|r|<1$, the CMGF is also faster and more
accurate than numerical integration based on the density (even if
we take advantage of the density being symmetric about zero). For
moments larger than one, integration with the density is about twice
as fast as CMGF, after taking advantage of the symmetric density.
A multivariate extension is beyond the scope of this paper. However,
some simple cross-moments, such as $\mathbb{E}[X_{1}X_{2}]$ and $\mathbb{E}[X_{1}X_{2}^{2}]$,
are easy to obtain from univariate results. Let $(X_{1},X_{2})$ have
MGF, $M_{X}(t_{1},t_{2})$, such that $M_{X_{1}+X_{2}}(t)=M_{X}(t,t)$,
$M_{X_{1}-X_{2}}(t)=M_{X}(t,-t)$, $M_{X_{1}}(t)=M_{X}(t,0)$, and
$M_{X_{2}}(t)=M_{X}(0,t)$, from which we obtain $\mathbb{E}[(X_{1}+X_{2})^{k}]$,
$\mathbb{E}[(X_{1}-X_{2})^{k}]$, etc. by Theorem 3, and the two cross
moments are given by $\mathbb{E}[X_{1}X_{2}]=\frac{1}{2}\left(\mathbb{E}[(X_{1}+X_{2})^{2}]-\mathbb{E}[X_{1}^{2}]\mathbb{E}[X_{2}^{2}]\right)$
and $\mathbb{E}[X_{1}X_{2}^{2}]=\frac{1}{6}\left(\mathbb{E}[(X_{1}-X_{2})^{3}]+\mathbb{E}[(X_{1}+X_{2})^{3}]-2\mathbb{E}[X_{1}^{3}]\right)$,
respectively.
In the next Section, we compute conditional moments in dynamic models
with the CMGF methods. The appropriate density is unknown in most
of these problems, and this makes the CMGF method particularly useful.
\section{Applications to Dynamic Models}
We consider three applications that demonstrate the usefulness of
the new CMGF expressions. In the first application, we use Theorems
\ref{thm:factional-abs-moments} and \ref{thm:integer-moments} to
obtain moments for cumulative returns from a Heston-Nandi GARCH model.
The second application uses Theorem \ref{cor:factional-moments-positive}
to compute fractional moments of realized volatility that follow an
Autoregressive Gamma process. The third application also applies Theorem
\ref{cor:factional-moments-positive} to obtain fractional moments,
but in this case for a discrete distribution that has probability
mass at zero.
A common feature of these three applications is that the relevant
MGF is available in closed form, while the corresponding density and
its derivatives are extremely difficult, if not impossible, to obtain.
Although simulation-based methods are always an option, they will
tend to be much slower than numerical integration of the new analytical
expressions. We will also implement simulation-based methods to compare
such with the new moment expressions, and compare their relative computational
burden.
\begin{comment}
We want to add a table that support this statement, what is the number
of simulations needed for accuracy (small standard error) and how
much time do then take relative to numerical CMGF integration.
\end{comment}
We use the following notation: $\mathbb{E}_{t}(\cdot)=\mathbb{E}(\cdot|\mathcal{F}_{t})$
is the conditional expectation where $\{\mathcal{F}_{t}\}$ is the
natural filtration. We will typically condition on $\mathcal{F}_{T}$,
and use $M_{X_{T+H}|T}(z)\equiv\mathbb{E}[\exp(X_{T+H}z)|\mathcal{F}_{T}]$
to denote the conditional MGF for the random variable $X_{T+H}\in\mathcal{F}_{T+H}$.
\subsection{Moments of Cumulative Returns in a GARCH Model}
The Heston-Nandi GARCH (HNG) model by \citet{HestonNandi2000} is
given by
\begin{equation}
\ensuremath{\begin{aligned}r_{t+1} & =r_{f}+\left(\lambda-\tfrac{1}{2}\right)h_{t+1}+\sqrt{h_{t+1}}z_{t+1},\qquad\text{with}\quad z_{t}\sim iidN(0,1),\\
h_{t+1} & =\omega+\beta h_{t}+\alpha\left(z_{t}-\gamma\sqrt{h_{t}}\right)^{2},
\end{aligned}
}\label{eq:HNG}
\end{equation}
where $r_{f}$ is the risk-free rate; $\lambda$ is the equity risk
premium; and $h_{t+1}={\rm var}\left(R_{t+1}|\mathcal{F}_{t}\right)$
is the conditional variance of daily log-returns $r_{t+1}=\log\left(S_{t+1}/S_{t}\right)$.
Among the many variations of GARCH models, the HNG model stands out
as having an analytical MGF for cumulative returns, $R_{T,H}=\sum_{t=T+1}^{T+H}r_{t}$,
where $H$ is the horizon in trading days. The underlying reason is
that its dynamic structure is carefully crafted for the purpose of
yielding closed-form option pricing formulae, and option pricing formulae
depend on the properties of cumulative returns, $R_{T,H}=\sum_{t=T+1}^{T+H}r_{t}$.
The conditional density function for $R_{T,H}$ does not have a known
analytical form, but the corresponding MGF is available in closed
form\footnote{We follow standard terminology and label an expression as ``closed
form'' if it can be expressed in terms of elementary mathematical
functions and operations, which is the case for our recursive expressions.} from simple recursive expressions, as stated below.
\begin{prop}
\label{prop:MGF_HNG}Let $r_{t}$, $t=1,2,\ldots$, be given by (\ref{eq:HNG})
and define cumulative returns, $R_{T,H}=r_{T+1}+\cdots+r_{t+H}$.
Then the conditional MGF has the affine form,
\begin{equation}
M_{R_{T,H}|T}(z)=\ensuremath{\exp\left(A(H,z)+B(H,z)h_{T+1}\right)},\label{eq:MGF_HNG}
\end{equation}
which is well-defined for $z\in\{\zeta\in\mathbb{C}:$$B(h,\operatorname{Re}(\zeta))<\frac{1}{2\alpha}$,
for all $h=1,\ldots,H\}$, where $A(H,z)$ and $B(H,z)$ are given
from the recursions,
\begin{align*}
A(h+1,z) & =A(h,z)+zr_{f}+B(h,z)\omega-\frac{1}{2}\log\left(1-2B(h,z)\alpha\right),\\
B(h+1,z) & =z(\lambda-\tfrac{1}{2})+B(h,z)\left(\beta+\alpha\gamma^{2}\right)+\frac{\left(z-2\alpha\gamma B(h,z)\right)^{2}}{2\left(1-2\alpha B(h,z)\right)},
\end{align*}
with initial values, $A(1,z)=zr_{f}$ and $B\left(1,z\right)=z(\lambda-\frac{1}{2})+\frac{z^{2}}{2}$.
\end{prop}
Obtaining the moments of cumulative return by way of the derivatives
of the MGF is nearly impossible, especially for higher moments and
for cumulative returns over many periods (large $H$). Instead, we
can compute the moments by Theorems \ref{thm:factional-abs-moments}
and \ref{thm:integer-moments}, where the former also enables us to
compute fractional absolute moments.
We will illustrate this with a simulation design that is based on
real data. We estimate the HNG model using daily log-returns for the
S\&P 500 index from January 1, 2000 to December 30, 2021. The data
were obtained from WRDS and the maximum likelihood estimates are presented
in Table \ref{tab:HNGstimation} along with their standard errors
(in parentheses).
\begin{table}[h]
\caption{Heston-Nandi GARCH Estimation Results for Daily Log-returns of SPX}
\begin{centering}
\vspace{0.2cm}
\begin{small}
\begin{tabularx}{\textwidth}{YYYYYYYYYYYY}
\toprule
\midrule
$\lambda$ & $\omega$ & $\beta$ & $\alpha$ & $\gamma$ & $\ell$ \\
\midrule
\\[-0.2cm]
1.9781 & 1.15 $\times 10^{-14}$ & 0.7593 & 5.67$\times 10^{-6}$ & 185.5 & \multicolumn{1}{c}{17,963} \\
(0.3166) & (0.85$\times 10^{-14}$) & (0.0190) & (1.18$\times 10^{-6}$) & (23.5) & \\
\\[-0.2cm]
\midrule
\bottomrule
\end{tabularx}
\end{small}
\par\end{centering}
{\small Note: Maximum likelihood estimates of the Heston-Nandi GARCH
model based on }daily S\&P 500 returns {\small with robust standard
errors in parentheses. Sample period: }January 1, 2000 to December
30, 2021{\small .\label{tab:HNGstimation}}{\small\par}
\end{table}
We can now illustrate Theorem \ref{thm:factional-abs-moments} by
computing the fractional absolute moments of cumulative returns, $R_{T,H}$
using the MGF we derived in Proposition \ref{prop:MGF_HNG}. Figure
\ref{fig:MomentsHNg-1} presents absolute moments, $\mathbb{E}|R_{T,H}|^{r}$,
for a range of $r\in(-1,4]$ and $H=21$ (1 month), $H=63$ (3 month),
and $H=126$ (6 months). The moments based on the new CMGF method
are shown with colored lines, and these agree with the dashed lines,
that are based on the average over 1,000,000 simulations. As for the
NIG distribution, the CMGF is more accurate and more than 100 times
faster than simulating 1,000,000 random and taking their average.
\begin{figure}
\begin{centering}
\includegraphics[width=0.8\textwidth]{Figures/HNGabs}
\par\end{centering}
\caption{{\small The conditional moments of absolute cumulative returns, $\mathbb{E}_{T}|R_{T,H}|^{r}$,
plotted against $r$, for horizons: one month ($H=21$), three months
($H=63$), and six months ($H=126$). The initial value of $h_{T+1}$,
is set as $\mathbb{E}(h_{t})$.\label{fig:MomentsHNg-1}}}
\end{figure}
We can also illustrate Theorem \ref{thm:integer-moments} in this
application, by computing the conditional integer moments, $\mathbb{E}_{T}[(R_{T,H})^{k}]$,
$k\in\{1,2,3,4\}$. From these moments, we compute the conditional
mean, $\ensuremath{\mu_{T,H}}$, standard deviation, $\sigma_{T,H}$,
skewness, and kurtosis, where $\ensuremath{\mu_{T,H}}=\mathbb{E}_{T}[R_{T,H}]$,
$\ensuremath{\sigma_{T,H}^{2}}=\mathbb{E}_{T}[R_{T,H}^{2}]-\mu_{T,H}^{2},$
\begin{eqnarray*}
\operatorname{Skew}_{T,H} & = & \frac{\mathbb{E}_{T}[R_{T+H}^{3}]-\mu_{T,H}^{3}}{\sigma_{T,H}^{3}}-3\frac{\mu_{t,H}}{\sigma_{T,H}},\\
\operatorname{Kurt}_{T,H} & = & \frac{\mathbb{E}_{T}[R_{T+H}^{4}]-\mu_{T,H}^{4}}{\sigma_{T,H}^{4}}-4\frac{\mu_{T,H}}{\sigma_{T,H}}\operatorname{Skew}_{T,H}-6\frac{\mu_{T,H}^{2}}{\sigma_{T,H}^{2}}.
\end{eqnarray*}
Figure \ref{fig:MomentsHNg} plots these quantities against $H$ (solid
lines) along with simulated quantities based on one million Monte
Carlo simulations (dashed lines), where the initial value of $h_{T}$
is set to be the unconditional mean, $\mathbb{E}(h_{t})$. We find
that the new expressions are in agreement with the simulated quantities.
\begin{figure}[th]
\begin{centering}
\includegraphics[width=0.8\textwidth]{Figures/HNG3}
\par\end{centering}
\caption{{\small Conditional moment quantities for cumulative returns, $R_{T,H}$
in the HNG model, plotted against $H$ that ranges from raging from
one day to six months. The initial value of $h_{T+1}$ is set to be
the unconditional mean, $\mathbb{E}(h_{t})$.\label{fig:MomentsHNg}}}
\end{figure}
\subsection{Moments in Autoregressive Gamma Model}
In this application, $X_{t}$ represents the daily realized variance,
which is computed from the intraday transaction data. We follow \citet{CorsiFusariVecchia2013}
and \citet{MajewskiBormettiCorsi2015} and adopt a Heterogeneous Autoregressive
Gamma (HARG) model for $X_{t}$. The HARG model is based on the Autoregressive
Gamma (ARG) model by \citet{GourierouxJasiak2006}, and both employ
a non-central Gamma distribution for the conditional density.\footnote{The transition density that is implied by the CIR model, see \citet{CoxIngersollRoss85},
is a non-central gamma density, and \citet{engle-gallo:06} employed
the standard Gamma distribution as the conditional distribution in
the MEM model.} What sets the two apart, is that the HARG employs a long-memory structure
for the location parameter. The conditional MGF conveniently has an
affine form, which is practical for evaluating the integral expressions
for the moments.
The HARG model has
\[
X_{t}|\mathcal{F}_{t-1}\sim f(x|\delta,\eta,\theta_{t})=\exp\left(-\frac{x}{\eta}-\theta_{t}\right)\left(\sum_{k=0}^{\infty}\frac{x^{\delta+k-1}}{\eta^{\delta+k}\Gamma(\delta+k)}\frac{\theta_{t}}{k!}\right),\quad x>0,
\]
which is the density for a non-central Gamma distribution with shape
parameter $\delta>0$, scale parameter $\eta>0$, and location parameter
$\theta_{t}>0$. The location parameter is time-varying, and given
by
\begin{equation}
\theta_{t}=\beta_{d}X_{t}^{(d)}+\beta_{w}X_{t}^{(w)}+\beta_{m}X_{t}^{(m)},\label{eq:HAR}
\end{equation}
where the variables, $X_{t}^{(d)}=X_{t-1}$, $X_{t}^{(w)}=\frac{1}{4}\sum_{i=2}^{5}X_{t-i}$,
and $X_{t}^{(m)}=\frac{1}{17}\sum_{i=6}^{22}X_{t-i}$ represent lagged
daily, ``weekly'', and ``monthly'' averages. Observe that $\theta_{t}$
is $\mathcal{F}_{t-1}$-measurable. The HARG model implies a restricted
AR(22) structure,
\begin{equation}
\mathbb{E}_{t-1}[X_{t}]=\eta\delta+\eta\theta_{t}=\mu+\sum_{j=1}^{22}\phi_{j}X_{t-j},\label{eq:RVmeanHAR}
\end{equation}
where $\mu=\eta\delta$ , $\phi_{1}=\eta\beta_{d}$, $\phi_{2}=\cdots=\phi_{5}=\tfrac{\eta}{4}\beta_{w}$,
and $\phi_{6}=\cdots=\phi_{22}=\tfrac{\eta}{17}\beta_{m}$, as in
the HAR model by \citet{Corsi:2009}.
The estimated parameters are reported in Table \ref{tab:HARGEstimation}
along with their standard errors, where $\phi_{d}=\eta\beta_{d}=\phi_{1}$,
$\phi_{w}=\eta\beta_{w}=\sum_{j=2}^{5}\phi_{j}$, and $\phi_{m}=\eta\beta_{m}=\sum_{j=6}^{22}\phi_{j}$,
such that $\phi_{d}+\phi_{w}+\phi_{m}=\sum_{j=1}^{22}\phi_{j}$ is
a measure of persistence.
The conditional MGF is conveniently given by
\begin{equation}
M_{X_{T+1}|T}(z)=\exp\left(\frac{\eta z}{1-\eta z}\theta_{T+1}-\delta\log(1-\eta z)\right),\qquad\text{for}\quad\operatorname{Re}(z)<\eta^{-1}.\label{eq:MGFone}
\end{equation}
From the estimated model we compute the term structure of conditional
moments, $\mathbb{E}_{T}[X_{T+H}^{r}]$. The first conditional moment,
$r=1$, is given directly from (\ref{eq:RVmeanHAR}), and we will
apply Theorem \ref{cor:factional-moments-positive} to obtain other
moments. The realized variance measures the second moment of returns,
such that $r=\frac{1}{2}$ corresponds to the conditional volatility
(standard deviation of returns). Similarly, $r=\frac{3}{2}$ and $r=2$
measures of the conditional skewness (of the absolute value) and the
conditional kurtosis, respectively. The inverse of volatility corresponds
to the negative moment, $r=-\frac{1}{2}$. This moment is interesting
because it is used for Sharpe ratio forecasting and asset allocation,
where the inverse of covariance matrix is employed. We compute these
four moments by Theorem \ref{cor:factional-moments-positive}. The
analytical form of the conditional MGF of $X_{T}$ is given in Proposition
\ref{prop:MGF_ARG}.
\begin{prop}
\label{prop:MGF_ARG}Suppose $X_{t}$ follows a HARG(p) process. Then
the conditional MGF for $X_{T+H}$, given $\mathcal{F}_{T}$, is given
by
\[
M_{X_{T+H}|T}(z)\equiv\mathbb{E}_{T}\left(e^{zX_{t+H}}\right)=\ensuremath{\exp\left(A(H,z)+\sum_{j=1}^{p}B_{j}(H,z)X_{T+1-j}\right)},
\]
where $A(H,z)$ and $B_{j}(H,z)$, are given from the initial values,
\[
A(1,z)=-\delta\log(1-\eta z)\quad\text{and}\quad B_{j}\left(1,z\right)=\frac{z}{1-\eta z}\phi_{i},\quad j=1,\ldots,p,
\]
and the recursions,
\begin{eqnarray*}
A(h+1,z) & = & A(h,z)-\delta\log[1-\eta B_{1}(h,z)],\\
B_{j}\left(h+1,z\right) & = & \tfrac{B_{1}(h,z)}{1-\eta B_{1}(h,z)}\phi_{j}+1_{\{j<p\}}B_{j+1}(h,z),\qquad j=1,\ldots,p.
\end{eqnarray*}
\end{prop}
Note that $M_{X_{T+H}|T}(z)$ is well defined for $z\in\{\zeta\in\mathbb{C}:B_{1}(h,\operatorname{Re}(\zeta))<\eta^{-1}\text{ for all }h\leq H\}$
and that $p=22$ in this model.
\begin{table}[th]
\caption{HARG Estimation Results for Daily Realized Variance}
\begin{centering}
\vspace{0.2cm}
\begin{small}
\begin{tabularx}{\textwidth}{YYYYYYYYYYYY}
\toprule
\midrule
$\tilde{\beta}_d$ & $\tilde{\beta}_w$ & $\tilde{\beta}_m$ & $\eta$ & $\delta$ & $\ell$ \\
\midrule
\\[-0.2cm]
0.4896 & 0.2789 & 0.0357 & 0.0053 & 0.9644 & {17,298} \\
(0.0279) & (0.0333) & (0.0120) & (0.0003) & (0.0226) & \\
\\[-0.2cm]
\midrule
\bottomrule
\end{tabularx}
\end{small}
\par\end{centering}
{\small Note: Maximum likelihood estimates for the HARG model with
robust standard errors in parentheses. We report $\tilde{\beta}_{j}=\eta\beta_{j}$
for $j=\{d,w,m\}$ because they (unlike $\beta_{j}$) can be interpreted
as the AR coefficients.\label{tab:HARGEstimation}}{\small\par}
\end{table}
We adopt the realized (Parzen) kernel estimator, $\operatorname{RK}_{t}$
by \citet{BNHLS-RK:2008}, as our realized measure of the variance.
The daily $\operatorname{RK}_{t}$ is estimated for the S\&P 500 index
over the sample period from January 1, 2000 to November 30, 2021 (5,490
trading days).\footnote{The data of realized kernel was obtained from the Realized Library
at the Oxford-Man Institute, which was discontinued in 2022.} The realized measures are annualized by the scaling, $X_{t}=252\times\operatorname{RK}_{t}$.
Table \ref{tab:HARGEstimation} presents the maximum-likelihood estimates
for the HARG models along with their corresponding standard errors
(in parentheses).
\begin{figure}[th]
\begin{centering}
\includegraphics[width=0.8\textwidth]{Figures/HARG2}
\par\end{centering}
\caption{{\small The conditional $r$-th moment, $\mathbb{E}_{T}[(X_{T+H})^{r}]$,
in HARG model for $r\in\{-\frac{1}{2},\frac{1}{2},\frac{3}{2},2\}$
and $H$ ranging from one day to six months. There is agreement between
the CMGF moments (solid lines) and the simulated moments based on
$N=10^{6}$ simulations (dashed line). The data generating process
is the HARG model with parameter values set to the estimates in Table
\ref{tab:HARGEstimation}. Conditional moments are compute for $X_{T}=X_{T-1}=\ldots=\frac{1}{10}\mathbb{E}[X_{t}]$,
where $\mathbb{E}\left(X_{t}\right)=\eta\delta/\left(1-\eta\beta_{d}-\eta\beta_{w}-\eta\beta_{m}\right)$.\label{fig:MomentsHARG}}}
\end{figure}
Under Proposition \ref{prop:MGF_ARG}, we compute the term structure
of moments $\mathbb{E}_{t}\left(X_{T+H}^{r}\right)$ for $H=1,\ldots,180$
and $r\in\{-\frac{1}{2},\frac{1}{2},\frac{3}{2},2\}$. Figure \ref{fig:MomentsHARG}
plots the term structure of these four moments when the initial value
of lagged $X_{t}$ components are all set as $\tfrac{1}{10}\mathbb{E}[X_{t}]$,
where $\mathbb{E}[X_{t}]=\eta\delta/\left(1-\eta\beta_{d}-\eta\beta_{w}-\eta\beta_{m}\right).$
We include both the simulated value (solid line) from one million
Monte-Carlo simulations and the numerical value (dashed line) from
Proposition \ref{prop:MGF_ARG}. We can find the numerical values
fit the simulated values very well.
\subsection{Moments in Autoregressive Poisson Model }
The autoregressive Poisson (ARP) model is given by
\begin{equation}
\Pr(Y_{t}=y|\mathcal{F}_{t-1})=\frac{\lambda_{t}^{y}}{y!}e^{-\lambda_{t}},\qquad y=0,1,\ldots\label{eq:Poisson}
\end{equation}
where the dynamic intensity parameter evolves according to
\begin{equation}
\lambda_{t+1}=\omega+\beta\lambda_{t}+\alpha Y_{t},\qquad t=1,2,\ldots,\label{eq:PoissonAR}
\end{equation}
such that $\lambda_{t+1}\in\mathcal{F}_{t}$. The average count over
the next $H$ periods is denoted,
\[
\bar{Y}_{T,H}=\frac{1}{H}\sum_{h=1}^{H}Y_{T+h},
\]
and we seek the conditional moments of $\bar{Y}_{T,H}$ given $\mathcal{F}_{T}$.
This could in principle be computed from the conditional distributions
of $(Y_{T+1},\ldots,Y_{T+H})$, which can be inferred from the ARP
model. However, there is substantial combinatorial complexity involved
with this, and the complexity increases rapidly with $H$. The new
method makes it simpler to compute moments, especially if $H$ is
large.
\begin{table}[th]
\caption{ARP Model Estimation Results for Count data of CBOE VIX Jumps}
\begin{centering}
\vspace{0.2cm}
\begin{small}
\begin{tabularx}{\textwidth}{YYYYYYYYYYYY}
\toprule
\midrule
$\omega$ & $\beta$ & $\alpha$ & $\mathbb{E}(\lambda)$ & $\ell$ \\
\midrule
\\[-0.2cm]
0.1548 & 0.7473 & 0.2043 & 3.2024 & -9722 \\
(0.0386) & (0.0349) & (0.0243) & & \\
\\[-0.2cm]
\midrule
\bottomrule
\end{tabularx}
\end{small}
\par\end{centering}
{\small Note: ML estimation results for the ARP model with robust standard
errors in parentheses. \label{tab:HARGEstimation-1}}{\small\par}
\end{table}
For $H=1$, the conditional MGF is
\begin{equation}
M_{Y_{T+1}|T}(z)=\exp\left(\lambda_{T+1}\left(e^{z}-1\right)\right).\label{eq:MGFone_ARP}
\end{equation}
More generally, the analytical form of the conditional MGF of $\bar{Y}_{T,H}$
is given in the following Proposition \ref{prop:MGF_ARP}.
\begin{prop}
\label{prop:MGF_ARP}Let $Y_{t}$ be given by (\ref{eq:Poisson})
and (\ref{eq:PoissonAR}). Then
\[
M_{\bar{Y}_{T,H}|T}(z)=\ensuremath{\exp\left(A(H,z)+B(H,z)\lambda_{T+1}\right)},
\]
where $A(H,z)$ and $B(H,z)$ are given from
\begin{align*}
A(h+1,z) & =A(h,z)+\omega B(h,z),\\
B(h+1,z) & =\beta B(h,z)+(e^{z/H+\alpha B(h,z)}-1),
\end{align*}
with initial value $A(1,z)=0$, and $B(1,z)=e^{z/H}-1$.
\end{prop}
We estimate an ARP model for daily volatility jumps. Let $Y_{t}$
be the number of daily volatility jumps, as defined by intraday jumps
in CBOE VIX index. We obtain high-frequency VIX data from Tick-data
for the sample period from July 1, 2003 to December 30, 2021. We use
the method by \citet{ABFN:2010} to identity the number of daily volatility
jumps.\footnote{A similar approach was used in \citet{AlitabBormettiCorsiMajewski:2020}
to determine the number of jumps in S\&P 500 index.} We estimated the ARP model with maximum likelihood. The estimate
parameters are presented in Table \ref{tab:HARGEstimation} along
with their standard errors (in parentheses). The average jump intensity
is about 3.2 jumps per day and $\hat{\pi}=\hat{\alpha}+\hat{\beta}=0.9517$
shows that the jump intensity is persistent.
\begin{figure}[th]
\begin{centering}
\includegraphics[width=0.8\textwidth]{Figures/ARP2}
\par\end{centering}
\caption{{\small Conditional moments, $\mathbb{E}_{t}[(\bar{Y}_{T,H})^{r}]$
in ARG model, where of $\bar{Y}_{T,H}=\frac{1}{H}\sum_{h=1}^{H}Y_{T+h}$.
We present the four moments, $r\in\{\frac{1}{2},1,\frac{3}{2},2\}$,
for $H$ ranging from 1 day to six months. We include both the simulated
value (solid line) from one million Monte-Carlo simulations and the
numerical value (dashed line) from Proposition \ref{prop:MGF_ARG}.
Model parameters are taken from Table \ref{tab:HARGEstimation-1}.
The initial value of $\lambda_{T+1}$ is set as $\frac{1}{10}\mathbb{E}\left(\lambda_{t}\right)$,
where $\mathbb{E}\left(\lambda_{t}\right)=\omega/\left(1-\beta-\alpha\right)$.
The horizontal axis indicates the calendar days.\label{fig:MomentsARP}}}
\end{figure}
Under Proposition \ref{prop:MGF_ARP}, we compute the term structure
of moments $\mathbb{E}_{T}[(\bar{Y}_{T,H})^{r}]$ for $r\in\{\frac{1}{2},1,\frac{3}{2},2\}$
and $H$ ranging from 1 day to six months. Figure \ref{fig:MomentsARP}
plots the term structure of these four moments when the initial value
of $\lambda_{t+1}$ set as $\frac{1}{10}\mathbb{E}\left(\lambda_{t}\right)$.
We include both the simulated value (solid line) from one million
Monte-Carlo simulations and the numerical value (dashed line) from
Proposition \ref{prop:MGF_ARP}. We can find the numerical values
fit the simulated values very well.
\section{Conclusion}
In this paper, we introduced a novel method for computing moments,
including fractional moments, of random variables using their moment-generating
function (MGF). A key advantage of our approach is that it avoids
the need for MGF derivatives, which can be computationally challenging
or unavailable in many models. We provided new integral expressions
for fractional moments, fractional absolute moments, and central moments
that extend the applicability of moment computation that is grounded
in the MGF. The CMGF method leverages a complex extension of the MGF
and is flexible enough to handle non-integer and complex moments.
The CMGF method may be valuable in structural models where moments
play an important role. Moments and conditional moments are also central
to inference methods. For instance, the generalized method of moments
(GMM) by \citet{Hansen1982} requires the mapping from parameter to
moments to be known. By offering solutions where other analytical
methods fall short, the CMGF method can broaden the applicability
of GMM to cover some problems that currently require simulated method
of moments (SMM), see \citet{McFadden:1989} and \citet{DufieSingleton:1993}.
We found the new method to be very fast and highly accurate for computing
moments of the normal-inverse Gaussian distribution. Moreover, the
CMGF method is especially useful in dynamic models where the MGF is
known but derivatives are difficult to obtain, as demonstrated by
the three applications in Section 3, where we computed moments of
cumulative returns in a Heston-Nandi GARCH model, moments of realized
volatilities in a Heterogeneous Autoregressive Gamma model, and moments
of number of volatility jumps in an Autoregressive Poisson model.
Future research could explore further extensions of this method to
other models with closed-form MGFs. An interesting extension is to
apply the CMGF method to multivariate distributions, allowing for
the computation of moments in multivariate distributions, including
cross-moments that capture dependencies, beyond the simplest cross-moments,
$\mathbb{E}[X_{1}X_{2}]$, discussed in Section 2.
\bibliographystyle{apalike}
\bibliography{prh}