Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
60,509 characters · 8 sections · 39 citation commands
Moments by Integrating the Moment-Generating Function
{Keywords:}{ Moments, Fractional Moments, Moment-Generating Function.}
{JEL Classification:}{ C02, C40, C65 }
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{Kawata1972 provides an expression for $\mathbb{E}|X|^{r}$ that involves several derivatives of the CF and Laue1980 provides an expression based on fractional derivatives, see also Wolfe:1975, SamkoKilbasMarichev:1993, MatsuiPawlas2016, and Tomovski2022. For fractional moments of non-negative variables, similar expressions were derived in CressieDavisFolksPolicello:1981, CressieBorkent1986, Jones1987, SchurgerK:2002, and 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), 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. 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 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 HansenTong:2025ReciprocalGamma.
The remainder of this paper is organized as follows. We present the new moment expressions in Section 2. Theorem (ref) has the expression for absolute fractional moments for general random variables, Theorem (ref) presents positive-part moments, and Theorem (ref) presents a new expression for integer moments. Corollary (ref) has expressions for fractional moments of nonnegative (positive) random variables and Corollary (ref) 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.
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}$.
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.
The identity in ((ref)) 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)) 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
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.
Here $x=0$, arises as a special case in the first two formulas of Lemma (ref), 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), 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.
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) 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)) 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)) rather than ((ref)).
Interestingly, Theorem (ref) can be used to obtain new identities, by equating ((ref)) 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)) yields the following expression for the reciprocal gamma function
for $z=s+it$, where $s>0$ is an arbitrary positive constant. From Theorem (ref) we can only infer that ((ref)) 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 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$.
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.
The expression in Theorem (ref) is identical to one in 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)) coincides with KimRachevBianchiFabozzi:2010 and ((ref)) coincides with HuangOosterlee:2011.\footnote{They derived the result under slightly stronger assumptions. In KimRachevBianchiFabozzi:2010 it is assumed that the distribution is infinitely divisible and continuous, while 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) and Lemma (ref).
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)) 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.
For random variables whose support may include negative numbers, we have the following result for integer moments.
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$.
Negative tail moments, which is part of the definition of expected shortfall (see below), have the following expression, where we can combine Theorems (ref) and (ref) and use the notation $x_{-}=\min(x,0)$.
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) 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}. \]
To illustrate the new method, we use Theorem (ref) to compute fractional moments of the normal-inverse Gaussian (NIG) distribution. The NIG distribution was introduced in 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 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 Barndorff-NielsenStelzer:2005 obtained the following explicit formula,
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) yields similar expressions for any other recentering ($\xi\neq\mu$), which is not possible with ((ref)). 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 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). 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.
We compute the absolute moments for $-0.85\leq r\leq4.2$ for the two distributions using Theorem (ref) and the moments are shown in the upper panel of Figure (ref). 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.
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 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). 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.
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) 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). 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.
We consider three applications that demonstrate the usefulness of the new CMGF expressions. In the first application, we use Theorems (ref) and (ref) to obtain moments for cumulative returns from a Heston-Nandi GARCH model. The second application uses Theorem (ref) to compute fractional moments of realized volatility that follow an Autoregressive Gamma process. The third application also applies Theorem (ref) 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.
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}$.
The Heston-Nandi GARCH (HNG) model by HestonNandi2000 is given by
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.
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) and (ref), 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) along with their standard errors (in parentheses).
We can now illustrate Theorem (ref) by computing the fractional absolute moments of cumulative returns, $R_{T,H}$ using the MGF we derived in Proposition (ref). Figure (ref) 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.
We can also illustrate Theorem (ref) 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},$
Figure (ref) 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.
In this application, $X_{t}$ represents the daily realized variance, which is computed from the intraday transaction data. We follow CorsiFusariVecchia2013 and MajewskiBormettiCorsi2015 and adopt a Heterogeneous Autoregressive Gamma (HARG) model for $X_{t}$. The HARG model is based on the Autoregressive Gamma (ARG) model by 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 CoxIngersollRoss85, is a non-central gamma density, and 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
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,
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 Corsi:2009.
The estimated parameters are reported in Table (ref) 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
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)), and we will apply Theorem (ref) 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). The analytical form of the conditional MGF of $X_{T}$ is given in Proposition (ref).
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.
We adopt the realized (Parzen) kernel estimator, $\operatorname{RK}_{t}$ by 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) presents the maximum-likelihood estimates for the HARG models along with their corresponding standard errors (in parentheses).
Under Proposition (ref), 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) 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). We can find the numerical values fit the simulated values very well.
The autoregressive Poisson (ARP) model is given by
where the dynamic intensity parameter evolves according to
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.
For $H=1$, the conditional MGF is
More generally, the analytical form of the conditional MGF of $\bar{Y}_{T,H}$ is given in the following Proposition (ref).
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 ABFN:2010 to identity the number of daily volatility jumps.\footnote{A similar approach was used in 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) 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.
Under Proposition (ref), 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) 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). We can find the numerical values fit the simulated values very well.
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 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 McFadden:1989 and 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.