Fractional Moments by The Moment-Generating Function: Peter Reinhard Hansen Chen Tong
Fractional Moments by The Moment-Generating Function: Peter Reinhard Hansen Chen Tong
Fractional Moments by The Moment-Generating Function: Peter Reinhard Hansen Chen Tong
Function∗
arXiv:2410.23587v1 [econ.EM] 31 Oct 2024
§
Department of Economics, University of North Carolina at Chapel Hill
‡
Department of Finance, School of Economics & Wang Yanan Institute
for Studies in Economics (WISE), Xiamen University
November 1, 2024
Abstract
∗
Corresponding author: Chen Tong, Email: tongchen@xmu.edu.cn. Chen Tong acknowledges fi-
nancial support from the Youth Fund of the National Natural Science Foundation of China (72301227),
and the Ministry of Education of China, Humanities and Social Sciences Youth Fund (22YJC790117).
1
1 Introduction
Moments of random variables, including conditional moments, are central to much of
economics, finance, and econometrics. 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 provide
solutions in settings where standard approaches are insufficient.
It is well known that integer moments of a random variable X, with MGF MX (s),
(k) (k)
are given by E[X k ] = MX (0), where MX (s) denotes the k-th derivative of the MGF.
Non-integer moments and fractional absolute moments, E|X r |, r ∈ R, have more com-
plicated expressions, typically integrals that involve derivatives of the MGF or deriva-
tives of the characteristic function (CF).1 In this paper, we derive novel integral expres-
sions for computing fractional moments, including fractional absolute moments, and
fractional central moments. Unlike existing expressions, the new expressions do not
involve derivatives (of the MGF or the CF). Thus, the new expressions are particularly
useful when derivatives are prohibitively difficult to obtain. This is often the case in
dynamic models, where the conditional MGF is defined from recursive non-linear ex-
pressions. The new method is labelled CMGF because it relies on a complex extension
of the MGF, and provides expressions for complex moments, E|X|r , r ∈ C. Although
complex moments are not commonly used in econometrics, this generalization can be
included in our results without much adaptation. Complex moments are widely applied
in fields such as statistical mechanics, number theory, and quantum physics.2 The new
expressions can also be used for establishing new identities, such as the expression for
the reciprocal Gamma function, which we infer from the absolute moment of Gaussian
random variables.
The new CMGF expressions involve integrals, and while these cannot be evaluated
analytically in many situations, they facilitate new ways to compute moments numeri-
cally, which are fast and accurate. We demonstrate this with the normal-inverse Gaus-
sian distribution, an we also illustrate the CMGF method with three well-motivated
applications. First, for the Heston-Nandi GARCH (HNG) model, we obtain moments
of cumulative returns. Second, we compute the moments of realized volatilities that
are given by the Heterogeneous Autoregressive Gamma (HARG) model. Third, within
1
Kawata (1972) provides an expression for E|X|r that involves several derivatives of the CF and
Laue (1980) provides an expression based on fractional derivatives, see also Wolfe (1975), Samko et al.
(1993), Matsui and Pawlas (2016), and Tomovski et al. (2022). For fractional moments of non-negative
variables, similar expressions were derived in Cressie et al. (1981), Cressie and Borkent (1986), Jones
(1987), and Schürger (2002).
2
In part, because complex moments are related to Mellin transforms.
2
the Autoregressive Poisson Model (APM), we compute the conditional moment for the
averages number of volatility jumps in the S&P 500 index. In these applications, we
are not aware of alternative methods that can compete with the CMGF method.
The remainder of this paper is organized as follows. We present the new moment
expressions in Section 2. Theorem 1 has the expression for absolute fractional moments
for general random variables, Theorem 2 has the expression for fractional moments of
non-negative random variables, and Theorem 3 presents a new expression for integer
moments. As an example, 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.
The MGF is typically used as a real function with a real argument, such that both s
and MX (s) take values in R. In this paper, we establish our results by means of complex
arguments, z ∈ C, such that MX (z) also takes values in C.
To gain some intuition for the complex extension of the MGF, we observe that
nests both the standard moment-generating function, ψ(s, 0) = MX (s), and the charac-
3
teristic function, ψ(0, t) = φX (t) ≡ E[eitX ], as special cases. The complex valued MGF,
MX (z), is somewhat similar to the Laplace transform when X is a continuous random
variable.3
The results are based on the following identities.
for any x ̸= 0, and this identity also holds for x = 0 if Re(r) > 0. For x > 0 we have
ˆ +∞
r Γ(r + 1) ezx
x = dt, Re(r) > −1.
2π −∞ z r+1
Similarly,
ˆ
k k! +∞
ezx + (−1)k e−zx
x = dt, k ∈ N0 ,
2π −∞ z k+1
for Re(r) ≥ 0, and if Pr(X = ξ) = 0, then the identity holds for all Re(r) > −1.
If r is real, then (1) is equivalent to
ˆ ∞
e−ξz MX (z) + eξz MX (−z)
" #
Γ(r + 1)
r
E|X − ξ| = Re dt. (2)
π 0 z r+1
´∞
3
The Laplace transform of the density function, fX (x), is 0 e−zx fX (x)dx, and it differs from
MX (z) in terms of the domain of integration and the sign of the exponent. Moreover, MX (z) is
obviously also applicable to discrete distributions and distributions that are mixtures of continuous
and discrete distributions.
4
The result obviously covers regular and central absolute moments by setting ξ = 0
for regular moments and ξ = EX for central moments. Note that Theorem 1 provides
a continuum of expressions, one for each s ∈ N+ . It is our experience that it makes
little difference which s ∈ N+ is used to evaluate the moment. Although analytical
integration may be impractical, (1) provides a new and alternative method to eval-
uate moments numerically. For real moments, r ∈ R, it is our experience that it is
computationally advantageous to use (2) rather than (1).
Interestingly, Theorem 1 can be used to obtain new identities, by equating (1) with
an existing expression for absolute moments. For instance, the absolute moment for
Γ( r+1 )
standard normal random variable, X ∼ N (0, 1), is E|X|r = √2π 2r/2 , for r > −1, and
equating this with (1) yields the following expression for the reciprocal gamma function
ˆ ∞ 2
1 2r/2 ez /2
= dt, r > −1, (3)
Γ( 2r + 1) π −∞ z r+1
for z = s + it, where s > 0 is an arbitrary positive constant. While we can only infer
that (3) holds for real r > −1 , the identity actually holds for any r ∈ C, see Hansen
and Tong (2024).
For non-negative random variables we have the following simpler expression for
fractional moments.
5
λ
X ∼ Exp(λ), we have MX (z) = λ−z , with N = (−∞, λ). So, the integral in (5) is
´∞ h
λ 1
i
0
Re λ−z z r+1
dt. For an integer moment, k ∈ N, this integral equals
ˆ ∞
k+1
ˆ ∞ h i
Re λ−k λ−z λ−k z1 λ−k+j−1 z1j dt −k
X
1 1 1
+ + =λ Re λ−z
+ z
dt
0 j=2 0
ˆ ∞
λ−z ∗ z∗
h i
−k
=λ Re (λ−z ∗ )(λ−z)
+ z∗ z
dt,
0
´∞
because 0
Re [z −j ] dt = 0 for j > 1. Substituting z = s + it in the last expression
gives us
ˆ ∞ h i ˆ ∞ h i
−k λ−s+it s−it −k λ−s s
λ Re (λ−s)2 −(it)2
+ s2 −(it)2
dt = λ (λ−s)2 +t2
+ s2 +t2
dt
0 0
h ∞i
= λ−k arctan( λ−s
t
) + arctan( st ) = πλ−k .
0
For random variables whose support may include negative numbers, we have the
following result for integer moments.
where z = s + it with s ∈ N+ .
´∞
Because k is real, we have expressed the integral on the form, k!
π 0
Re[·]dt, which
tends to be simpler to evaluate numerically than the equivalent expression using the
k!
´ +∞
form, 2π −∞
[·]dt.
Remark. In some situations, Theorems 1, 2, and 3, provide different expressions for the
same moments. For instance, (4) and (1) must agree for a non-negative random variable,
because E|X|r = E[X r ] in this case, and for integer moments, these expressions must
also agree with (6). This is indeed the case, because the additional term in (1) and the
additional term in (6) are both zero for non-negative variables. This is an implication
´ ∞ e−(s+it)x
of Lemma A.2, which establishes that −∞ (s+it) r+1 dt = 0 for all x > 0 and all s > 0.
We can combine Theorems 1 and 3 to compute odd moments over tail events.
6
Corollary 1 (Tail Moments). Let X satisfy Assumption 1, ξ ∈ R, and let k be odd
number. Then
ˆ
k! ∞ e−ξz MX (z)
" #
k
E[(X − ξ) 1{X>ξ} ] = + Re dt,
π 0 z k+1
ˆ
k! ∞ e−ξz MX (z)
" #
k
E[(X − ξ) 1{X<ξ} ] = − Re dt,
π 0 z k+1
We can also compute the cdf, quantiles, and expected shortfall using the new mo-
ment expressions. We can obtain the cdf for a continuous random variable, X, from
the MGF by Gil-Pelaez inversion,
ˆ ∞
MX (it)e−itx
" #
1 1
FX (x) = − Re dt,
2 π 0 it
where we used that MX (it) is the characteristic function for X. Now the α-quantile is
given by ξα ≡ FX−1 (α), and we can use Corollary 1 (with k = 1) to compute,
1 1
ESα (X) = − E[X1{X<ξα } ] = − E[(X − ξα )1{X<ξα } ] − ξα ,
α α
√
where γ = α2 − β 2 and K1 (·) is the modified Bessel
q
function of the second kind, and
h i
the MGF is given by MX (z) = exp λz + δ γ − α − (β − z)2 , for all z ∈ C. We
2
consider two standardized NIG distribution (zero mean and unit variance), which can
be characterized by the two parameters, (ξ, χ), for 0 ≤ |χ| < ξ < 1, where λ = − χξ ζ, δ =
√ √
ξ 2 ξ 2 − χ2 ζ, α = ξζ, and β = χζ, with ζ = 1 − ξ 2 /(ξ 2 − χ2 ), see Barndorff-Nielsen
7
et al. (1985). Specifically we consider (ξ, χ) = (1/2, −1/3) and (ξ, χ) = (1/8, −1/16),
the corresponding densities are shown in the left panel of Figure 1. The corresponding
MGFs, t 7→ MGF(s + it) are shown in the right panel for s = 1 – solid lines for the real
part and dashed lines for the complex part.
Figure 1: Densities of two standardized NIG distributions (left panel) and the t 7→MGF(1+
it) (right panel). The two distributions are for (ξ, χ) = (1/2, −1/3) (blue lines) and (ξ, χ) =
(1/8, −1/16) (red lines).
We compute the absolute moments for −0.85 ≤ r ≤ 4.2 for both distributions using
Theorem 1 and the moments are shown in the upper panel of Figure 2. 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χ2 )/(1 − ξ 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 compute the absolute moments using N = 1, 000, 000 independent
draws from the two NIG distributions. The simulation based moment is given by
µ̂rN = N1 N r 2 r
j=1 |Xi | and we estimate the σ (r) = var(|Xj | ) with 100 million draws. The
P
accuracy of a simulated moment can be computed from its standard error. Assuming
an error equal to one standard deviation, then the number of accurate decimal places is
√
given by − log10 (σ(r)/ N ).4 We have plotted this number for N = 106 in the lower left
4
For instance, an error where the first three decimal places are zero, say 0.00099, will map to
− log10 (0.99 · 10−3 ) ≃ 3.004.
8
Figure 2: Absolute moments of the standardized NIG distribution with (ξ, χ) =
(1/2, −1/3) (blue lines) and (ξ, χ) = (1/8, −1/16) (red lines). These are shown in the
left panels, where dots are represent the known moments for r = 2 and r = 4. The middle
panel is a snippet of the left panel which can be compare to simulation-based moments for 15
values of r. For each r we computed E|X|r 100 times as represented by the x-crosses. Each
of these simulation-based moments are based on N =1,000,000 draws of NIG distribution.
The right panel shows how many decimal places are accurate with a one standard deviation
simulation error.
panel of Figure 2. 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 ≤ r ≤ 4.2 represented by the
solid blue line. The ×-crosses are simulation-based moment estimates, each 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 1 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 −∞ |x|r f (x)dx using
the same numerical integration method used to evaluate the integral in Theorem 1.
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 = 106 independent draws from the NIG distribution. We have also compared the
methods using a standard normal distribution, for which the true moment is known
Γ( r+1 )
to be, E|X|r = √2π 2r/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
9
Table 1: CMGF Computation Time of E|X|r for X ∼ NIG
accurate. For |r| < 1, the CMGF is also faster and more accurate than integration with
the density, even after exploiting the symmetry of the density about zero. With this
modification for higher moments, density integration is about twice as fast as CMGF
for moments larger than one.
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.
10
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: Et (·) = E(·|Ft ) is the conditional expectation where
{Ft } is the natural filtration. We will typically condition on FT , and use MXT +H |T (z) ≡
E[exp(XT +H z)|FT ] to denote the conditional MGF for the random variable XT +H ∈
FT +H .
where rf is the risk free rate; λ is the equity risk premium; and ht+1 = var (Rt+1 |Ft ) is
the conditional variance of daily log-returns rt+1 = log (St+1 /St ).
Among the many variations of GARCH models, the HNG model stands out as having
an analytical MGF for cumulative returns, RT,H = Tt=T
P +H
+1 rt . 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, RT,H = Tt=T
P +H
+1 rt . The conditional density function for RT,H does
not have a known analytical form, but the corresponding MGF is available in closed
form5 from simple recursive expressions, as stated below.
1
which is well-defined for z ∈ {ζ ∈ C :B(h, Re(ζ)) < 2α
, for all h = 1, . . . , H}, where
5
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.
11
A(H, z) and B(H, z) are given from the recursions,
1
A(h + 1, z) = A(h, z) + zrf + B(h, z)ω − log (1 − 2B(h, z)α) ,
2
1
2
(z − 2αγB(h, z))2
B(h + 1, z) = z(λ − 2 ) + B(h, z) β + αγ + ,
2 (1 − 2αB(h, z))
z2
with initial values, A(1, z) = zrf and B (1, z) = z(λ − 21 ) + 2
.
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 1 and 3,
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 es-
timate 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 maxi-
mum likelihood estimates are presented in Table 2 along with their standard errors (in
parentheses).
λ ω β α γ ℓ
Note: Maximum likelihood estimates of the Heston-Nandi GARCH model based on daily
S&P 500 returns with robust standard errors in parentheses. Sample period: January 1,
2000 to December 30, 2021.
12
Figure 3: The conditional moments of absolute cumulative returns, ET |RT,H |r , plotted
against r, for horizons: one month (H = 21), three months (H = 63), and six months
(H = 126). The initial value of hT +1 , is set as E(ht ).
the conditional mean, µT,H , standard deviation, σT,H , skewness, and kurtosis, where
2 2
µT,H = ET [RT,H ], σT,H = ET [RT,H ] − µ2T,H ,
Figure 4 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
hT is set to be the unconditional mean, E(ht ). We find that the new expressions are in
agreement with the simulated quantities.
13
Figure 4: Conditional moment quantities for cumulative returns, RT,H in the HNG model,
plotted against H that ranges from raging from one day to six months. The initial value of
hT +1 is set to be the unconditional mean, E(ht ).
model is based on the Autoregressive Gamma (ARG) model by Gourieroux and Jasiak
(2006), and both employ a non-central Gamma distribution for the conditional density.6
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δ+k−1
! !
x X θt
Xt |Ft−1 ∼ f (x|δ, η, θt ) = exp − − θt δ+k
, x > 0,
η k=0 η Γ(δ + k) k!
which is the density for a non-central Gamma distribution with shape parameter δ > 0,
scale parameter η > 0, and location parameter θt > 0. The location parameter is
6
The transition density that is implied by the CIR model, see Cox et al. (1985), is a non-central
gamma density, and Engle and Gallo (2006) employed the standard Gamma distribution as the con-
ditional distribution in the MEM model.
14
time-varying, and given by
22
X
Et−1 [Xt ] = ηδ + ηθt = µ + ϕj Xt−j , (10)
j=1
η
where µ = ηδ , ϕ1 = ηβd , ϕ2 = · · · = ϕ5 = η4 βw , and ϕ6 = · · · = ϕ22 = 17 βm , as in the
HAR model by Corsi (2009).
The estimated parameters are reported in Table 3 along with their standard errors,
where ϕd = ηβd = ϕ1 , ϕw = ηβw = 5j=2 ϕj , and ϕm = ηβm = 22
P P
j=6 ϕj , such that
P22
ϕd + ϕw + ϕm = j=1 ϕj is a measure of persistence.
The conditional MGF is conveniently given by
!
ηz
MXT +1 |T (z) = exp θT +1 − δ log(1 − ηz) , for Re(z) < η −1 . (11)
1 − ηz
From the estimated model we compute the term structure of conditional moments,
ET [XTr +H ]. The first conditional moment, r = 1, is given directly from (10), and we
will apply Theorem 2 to obtain other moments. The realized variance measures the
second moment of returns, such that r = 21 corresponds to the conditional volatility
(standard deviation of returns). Similarly, r = 32 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 = − 12 . This moment is interesting
because it is used for Shape ratio forecasting and asset allocation, where the inverse of
covariance matrix is employed. We compute these four moments by Theorem 2. The
analytical form of the conditional MGF of XT is given in Proposition 2.
15
where A(H, z) and Bj (H, z), are given from the initial values,
z
A(1, z) = −δ log(1 − ηz) and Bj (1, z) = ϕi , j = 1, . . . , p,
1 − ηz
Note that MXT +H |T (z) is well defined for z ∈ {ζ ∈ C : B1 (h, Re(ζ)) < η −1 for all h ≤
H} and that p = 22 in this model.
Note: Maximum likelihood estimates for the HARG model with robust standard errors in
parentheses. We report β̃j = ηβj for j = {d, w, m} because they (unlike βj ) can be interpreted
as the AR coefficients.
16
Figure 5: The conditional r-th moment, ET [(XT +H )r ], in HARG model for r ∈ {− 21 , 12 , 32 , 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 = 106 simulations (dashed line). The
data generating process is the HARG model with parameter values set to the estimates in
1
Table 3. Conditional moments are compute for XT = XT −1 = . . . = 10 E[Xt ], where E (Xt ) =
ηδ/ (1 − ηβd − ηβw − ηβm ).
λyt −λt
Pr(Yt = y|Ft−1 ) = e , y = 0, 1, . . . (12)
y!
such that λt+1 ∈ Ft . The average count over the next H periods is denoted,
H
1 X
ȲT,H = YT +h ,
H h=1
17
and we seek the conditional moments of ȲT,H given FT . This could in principle be
computed from the conditional distributions of (YT +1 , . . . , YT +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.
Table 4: ARP Model Estimation Results for Count data of CBOE VIX Jumps
ω β α E(λ) ℓ
Note: ML estimation results for the ARP model with robust standard errors in parentheses.
More generally, the analytical form of the conditional MGF of ȲT,H is given in the
following Proposition 3.
We estimate an ARP model for daily volatility jumps. Let Yt 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 Andersen et al. (2010) to identity the number
of daily volatility jumps.8 We estimated the ARP model with maximum likelihood.
8
A similar approach was used in Alitab et al. (2020) to determine the number of jumps in S&P 500
index.
18
The estimate parameters are presented in Table 3 along with their standard errors (in
parentheses). The average jump intensity is about 3.2 jumps per day and π̂ = α̂ + β̂ =
0.9517 shows that the jump intensity is persistent.
1 PH
Figure 6: Conditional moments, Et [(ȲT,H )r ] in ARG model, where of ȲT,H = H h=1 YT +h .
We present the four moments, r ∈ { 21 , 1, 32 , 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 2. Model parameters are taken from Table
1
4. The initial value of λT +1 is set as 10 E (λt ), where E (λt ) = ω/ (1 − β − α). The horizontal
axis indicates the calendar days.
19
4 Conclusion
In this paper, we introduced a novel method for computing moments, including frac-
tional 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 Hansen (1982)
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 appli-
cability of GMM to cover some problems that currently require simulated method of
moments (SMM), see McFadden (1989) and Duffie and Singleton (1993).
We found the new method to be very fast and highly accurate for computing mo-
ments of the normal-inverse Gaussian distribution. Moreover, the CMGF method is
especially useful in dynamic models where the MGF is known but derivatives are dif-
ficult 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 the MGF in closed form. An interesting extension is to apply the CMGF method
to multivariate settings, allowing for the computation of moments in multivariate dis-
tributions, including cross-moments that capture dependencies.
References
Alitab, D., Bormetti, G., Corsi, F., and Majewski, A. A. (2020). A jump and smile ride:
Jump and variance risk premia in option pricing. Journal of Financial Econometrics,
18:121–157.
Andersen, T. G., Bollerslev, T., Frederiksen, P., and Ørregaard Nielsen, M. (2010).
20
Continuous-time models, realized volatilities, and testable distributional implications
for daily stock returns. Journal of Applied Econometrics, 25:233–261.
Barndorff-Nielsen, O. E., Blæsid, P., Jensen, J. L., and Sørensen, M. (1985). The
fascination of sand. In Atkinson, A. and Feinberg, S., editors, A Celebration of
Statistics. Sprin, New York.
Barndorff-Nielsen, O. E., Hansen, P. R., Lunde, A., and Shephard, N. (2008). Designing
realised kernels to measure the ex-post variation of equity prices in the presence of
noise. Econometrica, 76:1481–536.
Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B. (2017). Julia: A fresh
approach to numerical computing. SIAM Review, 59:65–98.
Corsi, F., Fusari, N., and Vecchia, D. L. (2013). Realizing smiles: Options pricing with
realized volatility. Journal of Financial Economics, 107:284–304.
Cox, J. C., Ingersoll, J. E., and Ross, S. A. (1985). A theory of the term structure of
interest rates. Econometrica, 53:385–408.
Cressie, N. and Borkent, M. (1986). The moment generating function has its moments.
Journal of Statistical Planning and Inference, 13:337–344.
Cressie, N., Davis, A. S., Folks, J. L., and Policello, G. E. (1981). The moment-
generating function and negative integer moments. American Statistician, pages
148–150.
Engle, R. F. and Gallo, G. (2006). A multiple indicators model for volatility using
intra-daily data. Journal of Econometrics, 131:3–27.
21
Gourieroux, C. and Jasiak, J. (2006). Autoregressive gamma processes. Journal of
Forecasting, 25:129–152.
Hansen, P. R. and Tong, C. (2024). A globally valid integral expression for the reciprocal
gamma function. Working paper.
Jones, M. (1987). Inverse factorial moments. Statistics & Probability Letters, 6:37–42.
Laue, G. (1980). Remarks on the relation between fractional moments and fractional
derivatives of characteristic functions. Journal of Applied Probability, 17:456–466.
Majewski, A. A., Bormetti, G., and Corsi, F. (2015). Smile from the past: A general
option pricing framework with multiple volatility and leverage components. Journal
of Econometrics, 187:521–531.
Matsui, M. and Pawlas, Z. (2016). Fractional absolute moments of heavy tailed distri-
butions. Brazilian Journal of Probability and Statistics, 30:272–298.
Samko, S. G., Kilbas, A. A., and Marichev, O. I. (1993). Fractional integrals and
derivatives (theory and applications). Gordon and Breach, Switzerland.
Tomovski, Ž., Metzler, R., and Gerhold, S. (2022). Fractional characteristic functions,
and a fractional calculus approach for moments of random variables. Fractional
Calculus and Applied Analysis, 25:1307–1323.
22
A Appendix of Proofs
Some notation: We write C+ = {z ∈ C : Re[z] > 0} to denote the complex numbers
with positive real part.
In the proof of Theorem 1, we will need the following result (the Laplace transfor-
mation of xr ).
Proof. When z ∈ R+ (such that z = s, t = 0) the result follows from the definition
´∞ ´ +∞
of the gamma function, Γ(r + 1) = 0 y r e−y dy = 0 sr xr e−sx sdx, where we used
integration by substitution with y = sx.
For the general situation, z = s + it ∈ C+ , we consider
ˆ K ˆ sK+itK
r −zx 1
lim xe dx = y r e−y dy,
K→∞ 0 z r+1 0
Im
C
tK
D
λ
θ
Re
A B
sK
¸
Figure 7: The closed contour for the integral, y r e−y dy.
The function, f (y), is analytic over the closed contour for λ > 0, (i.e. it does not
23
include any singularity). Therefore, by Cauchy’s Integral Theorem, we have
˛ ˆ ˆ ˆ ˆ
0= f (y)dy = f (y)dy + f (y)dy + f (y)dy + f (y)dy,
AB BC CD DA
´
where the integral we seek is − CD f (y)dy as λ → 0 and K → ∞.
For the integral over the arc, DA, we have y = λeiθ = λ (cos θ + i sin θ), and dy =
iλeiθ dθ, such that
ˆ ˆ 0
r −y
|IDA | = y e dy = λr eiθr e−λ(cos θ+i sin θ) iλeiθ dθ
DA π/2
ˆ 0
|λ|
≤ λ r+1
e eiθr e−iλ sin θ ieiθ dθ
π/2
ˆ 0
≤ λr+1 e|λ| dθ = π2 λRe(r)+1 e|λ| → 0, as λ → 0 for Re (r) > −1.
π/2
because e−sK decays faster than K r+1 increases. So, as K → ∞ and λ → 0, we have
ˆ ˆ ! ˆ sK ˆ 0
!
0 = lim lim f (y)dy + f (y)dy = lim f (y)dy + f (y)dy .
K→∞ λ→0 AB CD K→∞ 0 sK+itK
24
Proof. We first consider the case x = 0. The result follows directly from,
ˆ +∞ ˆ +∞ +∞
ez·0 1 i 1
dt = dt = ,
−∞ z r+1 −∞ (s + it)r+1 r (s + it)r −∞
−iτ x
e
where f (τ ) = (s−iτ )r+1
. The seemingly odd change of variables, τ = −t, serves two
important purposes that we explain below.
We will use a closed contour integral in the complex plane with two part. The
first part is the line segment along the real axis, from −R to R, and the second part
is the semicircular arc, CR = {Reiθ , θ ∈ [0, π]}. The integral over the line segment,
´ +R
−R
f (τ )dτ , is the object of interest, whereas the semicircle defines an auxiliary integral
that is used to close the contour, such that we can apply Cauchy’s Integral Theorem.
The first reason for the change of variable, τ = −t, is to place the pole of f (ζ) outside
the closed contour for all R > 0 (the pole is at ζ = −is), and by Cauchy’s Integral
´R ´
Theorem, we have 0 = −R f (ζ)dζ + CR f (ζ)dζ. For the integral along the arc, we have
the expression, ˆ ˆ ˆ
f (ζ)dζ = e−iζx g(ζ)dζ = eiaζ g(ζ)dζ,
CR CR CR
1
where a = −x > 0 and g(ζ) = (s−iζ) r+1 is continuous and satisfies lim|ζ|→∞ g(ζ) = 0
+
for all r + 1 ∈ C . The second reason for the change of variable is that it enables us
to apply Jordan’s lemma, which requires a > 0 in the expression above. By Jordan’s
25
´
lemma we have that limR→∞ CR
e−ixζ g(ζ)dζ = 0, which implies that
ˆ +∞ ˆ +R ˆ
f (τ )dτ = lim f (τ )dτ = − lim f (ζ)dζ = 0.
−∞ R→∞ −R R→∞ CR
Proof of Lemma 1. We first focus on the case where x ≥ 0. With Re[z] = s > 0,
from Lemma A.1, we have the identity
ˆ ∞
Γ(r + 1)
xr e−zx dx = , for all r + 1 ∈ C+ .
0 z r+1
provided the integral is well defined. Using x > 0, we can rewrite the integral as
ˆ
r Γ(r + 1) σ+i∞ ezx
r
|x| = x = r+1
dz,
2πi σ−i∞ z
ˆ
Γ(r + 1) +∞ ezx
= r+1
dt, (A.1)
2π −∞ z
which is well-defined for all r + 1 ∈ C+ (because x > 0). For x = 0 it is well defined for
all r ∈ C+ . This finishes the proof of the third formula in Lemma 1. We proceed to get
ˆ ˆ
r Γ(r + 1) +∞ ezx
r Γ(r + 1) +∞ e−zx
|x| = x = r+1
dt + r+1
dt
2π −∞ z 2π −∞ z
ˆ
Γ(r + 1) +∞ ezx + e−zx
= dt, x ≥ 0 (A.2)
2π −∞ z r+1
26
Similarly, we proceed to get
ˆ ˆ
r r Γ(r + 1) +∞ e−zx Γ(r + 1) +∞ ezx
|x| = (−x) = r+1
dt + r+1
dt
2π −∞ z 2π −∞ z
ˆ
Γ(r + 1) +∞ e−zx + ezx
= dt, x < 0, (A.3)
2π −∞ z r+1
by using the similar steps in the proof of first formula in Lemma 1, we have
ˆ
k k! +∞
ezx + (−1)k e−zx
x = dt,
2π −∞ z k+1
"ˆ ∞
ezx
#
r Γ(r + 1) X
E[X ] = r+1
dt p(x)
2π x∈X −∞ z
ˆ
Γ(r + 1) ∞ X ezx
= r+1
p(x)dt
2π −∞ x∈X z
ˆ
Γ(r + 1) ∞ MX (z)
= r+1
dt,
2π −∞ z
27
arguments applies for continuous distribution, given by
ˆ ∞
"ˆ ∞ # ˆ ∞
r Γ(r+1) ezx Γ(r+1) MX (u)
E[X ] = 2π z r+1
dt f (x)dx = 2π z r+1
dt.
−∞ −∞ −∞
where p = 22. For T = T + H, assume that the conditional MGF for XT +H is given by
p
MT,H (z) = ET ezXT +H = exp A(H, z) +
X
Bj (H, z) XT +1−j .
j=1
ηz
A(1, z) = −δ log(1 − ηz), Bj (1, z) = βj , j = 1, . . . , p.
1 − ηz
28
The recursions (from H to H + 1) follow from
h i
ET (ezxT +H ) = ET ET +1 ezXT +H
p
X
= ET exp A(H − 1, z) + Bi (H − 1, z) XT +2−j
j=1
p
X
!
= exp A(h, z) + Bi (h, z) xt+2−i Et [exp (B1 (h, z) xt+1 )]
i=2
p−1
X
= exp A(h, z) + Bi+1 (h, z) xt+1−i
i=1
p !
ηB1 (h, z) X
× exp βi xt+1−i − δ log(1 − ηB1 (h, z))
1 − ηB1 (h, z) i=1
p
X
!
= exp A(h + 1, z) + Bi (h + 1, z) xt+1−i ,
i=1
where
A(h + 1, z) = A(h, z) − δ log(1 − ηB1 (h, z))
and
Bi+1 (h, z) + ηB1 (h,z)
β 1≤i<p
1−ηB1 (h,z) i
Bi (h + 1, z) =
ηB1 (h,z)
β i=p
1−ηB1 (h,z) i
H
1 X
ȲT,H = ST,H , where ST,H = YT +h .
H h=1
29
We will show, by induction, that MST,H (z) has an affine structure. For H = 1, we
have ST,1 = YT +1 , such that MST,1 (z) = Et [ezYT +1 ] = exp[λT +1 (ez − 1)], which can be
expressed in the affined form,
MST,1 (z) = exp Ã(1, z) + B̃(1, z)λT +1 , with Ã(1, z) ≡ 0 and B̃(1, z) ≡ ez − 1.
Next, suppose that MST,h (z) = exp(Ã(h, z) + B̃(h, z)λT +1 ) and consider
This proves that MST,h+1 (z) has an affine structure and, by an induction argument, so
does MST,h+2 (z), . . . , MST,H (z). Next, MȲT,H (z) = MST,H (z/H). So, if we set A(h, z) ≡
Ã(h, z/H) and B(h, z) = B(h, z/H), then MȲT,H (z) = exp[A(H, z) + B(H, z)]λT +1
has the affine structure, where A(h + 1, z) = A(k, z) + ωB(k, z) and B(h + 1, z) =
βB(h, z) + ez/H+αB(h,z) − 1 , with initial values A(1, z) = 0 and B(1, z) = (ez/H − 1).
□
30
B.2 The MGF in Heston-Nandi GARCH Model
We first have
h h i i
Et [exp (v1 rt+1 + v2 ht+2 )] = exp v1 r + v2 ω + v1 λ − 1
2
+ v2 β + αγ 2 ht+1
q
2
× Et v2 αzt+1 + (v1 − 2v2 αγ) ht+1 zt+1
h h i
= exp v1 r + v2 ω + v1 λ − 1
2
+ v2 β + αγ 2 ht+1
2 #
1 (v1 − 2v2 αγ)
− log (1 − 2v2 α) + ht+1
2 2 (1 − 2v2 α)
h
= exp v1 r + v2 ω − 21 log (1 − 2v2 α)
(v1 − 2v2 αγ)2
! #
1 2
+ v1 λ − 2
+ v2 β + αγ + ht+1 .
2 (1 − 2v2 α)
Next, we show by induction that the conditional MGF has the form,
H
!!
X
Et exp s rt+i = exp (A(H, s) + B(H, s)ht+1 ) .
i=1
For h = 1 we have
1 s2
Et (exp (srt+1 )) = exp sr + s λ − 2
+ 2
ht+1 = exp (A(1, s) + B(1, s)ht+1 )
s2
such that A(1, s) = sr and B(1, s) = s λ − 21 + 2
. Next, suppose that the MGF has
an affine structure for h, then for h + 1 we find
h+1 h+1
!! " !!#
X X
Et exp s rt+i = Et Et+1 exp s rt+i
i=1 i=1
h
" !!#
X
= Et exp (sRt+1 ) Et+1 exp s rt+1+i
i=1
31
where
1
A(h + 1, s) = A(h, s) + sr + B(h, s)ω − log (1 − 2B(h, s)α)
2
1 (s − 2B(h, s)αγ)2
B(h + 1, s) = s λ − + B(h, s) β + αγ 2 + ,
2 2 (1 − 2B(h, s)α)
1
for any sufficiently small s = Re(z) > 0, such that B(h, s) < 2α for all h = 1, . . . , H, to
ensure that log (1 − 2B(h, s)α) ∈ R. This completes the proof. □
32
Supplement
NIG Moments
The shape of the integrands provide some intuition for the advantage of the CMGF
method in the evaluation of moments of the NIG distribution with (ξ, χ) = ( 12 , − 31 ).
The CMGF integrand is non-negligible over a domain that is very similar for all values
of r and the range of the integrand is also similar across r. In fact, all CMGF integrands
have the same maximum value at zero. For the integrand, |x|r f (x), there is far more
variation across different moments, r, which may explain that this method is slower
than the CMGF method.
S.1
where K = ⌊ 2r ⌋and CK is a positive constant, see Kawata (1972, theorem 11.4.4), and
ˆ ∞ (k) (k)
φX (0) − φX (u)
r
E|X| = λ
Γ(1−λ)
[cos( πr
2
)]−1 du, (S.1)
0 u1+λ
where k = ⌊r⌋ and λ = r − k, see Schürger (2002), and another expression for non-
negative variables is
ˆ ∞
r 1 (k̃)
E[X ] = uλ̃−1 MX (−u)du, (X ≥ 0),
Γ(λ̃) 0
where k̃ = ⌈r⌉ > 0 and λ̃ = k̃ − r ∈ [0, 1), see Cressie and Borkent (1986). For strictly
positive variables we have the following expression for negative fractional moments,
r < 0, ˆ ∞
r 1 1
E[X ] = r+1
MX (−t)dt, (X > 0),
Γ(−r) 0 t
see Schürger (2002, theorem 1.1).
S.2
Some Computer Code
Moment of NIG distribution
S.3
Alpha = model (3) ;
Gamma = model (4) ;
Lam = model (5) ;
rf = model (6) ;
h1 = model (7) ;
uR = 5;
u = uR + s * complex (0 ,1) ;
v = -u ;
A1 = u * rf ;
B1 = u *( Lam -0.5) + u .^2/2;
A2 = v * rf ;
B2 = v *( Lam -0.5) + v .^2/2;
for m = 1: H -1
A1 = A1 + u * rf + B1 * Omega -0.5* log (1 -2* Alpha * B1 ) ;
A2 = A2 + v * rf + B2 * Omega -0.5* log (1 -2* Alpha * B2 ) ;
B1 = u *( Lam -0.5) + B1 *( Beta + Alpha * Gamma ^2) +0.5*( u -2* B1 * Alpha * Gamma ) .^2./( 1 -2* Alpha * B1 ) ;
B2 = v *( Lam -0.5) + B2 *( Beta + Alpha * Gamma ^2) +0.5*( v -2* B2 * Alpha * Gamma ) .^2./( 1 -2* Alpha * B2 ) ;
end
fs = real (1./( u .^( r +1) ) .*( exp ( A1 + B1 * hp1Q ) + exp ( A2 + B2 * hp1Q ) ) ) ;
end
S.4