Fractional Moments by The Moment-Generating Function: Peter Reinhard Hansen Chen Tong

Download as pdf or txt
Download as pdf or txt
You are on page 1of 36

Fractional Moments by the Moment-Generating

Function∗
arXiv:2410.23587v1 [econ.EM] 31 Oct 2024

Peter Reinhard Hansen§ and Chen Tong‡

§
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

We introduce a novel method for obtaining a wide variety of moments of a ran-


dom 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 new moment 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 corre-
sponding densities and the derivatives of the MGF are either unavailable or very
difficult to obtain.

Keywords: Fractional Moments, Moment-Generating Function.


JEL Classification: C02, C40, C65


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.

2 A New Method for Computing Moments (CMGF)


We will use the following notation: N denotes the positive integers and N0 = {0} ∪ N.
We use N to denote a neighborhood about zero, for which (−δ, δ) ⊂ N for some δ > 0,
and the positive part of N is denoted N+ ≡ {s ∈ N : s > 0}. We use z = s + it
and ζ = σ + iτ to represent complex variables, where s = Re(z) and t = Im(z), and
similarly for ζ.
To avoid possible confusion, we use E[X k ] to denotes integer moments, k ∈ N,
whereas general moments (including fractional and complex moments) are represented
with E[X r ], r ∈ C.
Our results are applicable to any random variable with a well-defined MGF, as
characterized in the following assumption.

Assumption 1. X is a real random variable with

MX (s) = E[esX ] < ∞, for all s ∈ N ,

where N is a neighborhood about zero.

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

ψ(s, t) ≡ MX (s + it) = E[e(s+it)X ] ∈ C,

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.

Lemma 1. Let x ∈ R and z = s + it ∈ C, where s > 0 is an arbitrary positive constant.


ˆ +∞
r Γ(r + 1) ezx + e−zx
|x| = dt, Re(r) > −1,
2π −∞ z r+1

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 any x ̸= 0, and this identity also holds for x = 0 if k ∈ N.

The case, x = 0, arises as a special case in Lemma 1, requiring results to include a


qualification for distributions with positive mass at zero.
Our first result is the following expression for fractional absolute (possibly complex)
moments.

Theorem 1 (Absolute moments). Suppose X satisfies Assumption 1, ξ ∈ R, and


z = s + it with s ∈ N+ . Then
ˆ +∞
r Γ(r + 1) e−ξz MX (z) + eξz MX (−z)
E|X − ξ| = dt, (1)
2π −∞ z r+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.

Theorem 2 (Moments of non-negative random variables). Let X satisfies Assumption


1 and Pr(X ≥ ξ) = 1 for some ξ ∈ R. For any r ∈ C with Re(r) ≥ 0, we have,
ˆ +∞
r Γ(r + 1) e−ξz MX (z)
E[(X − ξ) ] = dt, (4)
2π −∞ z r+1

where z = s + it with s ∈ N+ . If X is strictly greater than ξ, i.e. Pr(X > ξ) = 1, then


the identity holds for Re(r) > −1.
If r is real, then (4) is equivalent to
ˆ +∞
e−ξz MX (z)
" #
Γ(r + 1)
E[(X − ξ)r ] = Re dt. (5)
π 0 z r+1

As an illustration, we apply (4) to an exponentially distributed random variable.


This example serves as an illustration because integer moments, k ∈ N, are easy to
(k)
obtain from the derivatives, MX (s), and other moments can be obtained by evaluating
´ ∞ r λx
0
x λe dx directly.

Example 1. For an exponentially distributed random variable, with parameter λ > 0,

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

Thus that E[X k ] = λ−k k! as expected. Deriving non-integer moments, r ∈


/ N, can
be derived using contour integrals and the Cauchy Integral Theorem. These are the
techniques used extensively to prove the main results in this paper.

For random variables whose support may include negative numbers, we have the
following result for integer moments.

Theorem 3. Let X satisfy Assumption 1 and let ξ ∈ R. For k ∈ N we have


ˆ ∞
e−ξz MX (z) + (−1)k eξz MX (−z)
" #
k! k
E[(X − ξ) ] = Re dt, (6)
π 0 z k+1

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

where z = s + it for any s ∈ N+ .

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<ξα } ] − ξα ,
α α

which is known as the expected shortfall.

2.1 Example: Moments of NIG Distribution


To illustrate the new method, we use Theorem 1 to compute fractional moments of the
normal-inverse Gaussian (NIG) distribution. The CMGF method is significantly faster
and more accurate than both simulation-based methods and direct integration using
the NIG density.
The NIG distribution was introduced in Barndorff-Nielsen (1978) and has four pa-
rameters, λ (location), δ (scale), α (tail heaviness), and β (asymmetry). It has density
 q 
f (x) = √ αδ
K1 α δ2 + (x − λ)2 eδγ+β(x−λ) , x∈R
π δ 2 +(x−λ)2


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

r CMGF Integate with Simulations


density (N = 106 )

-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 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 100.3 318,081

Note: Computation time (in microseconds) for evaluating E|X| r


´ ∞ using three methods: the
CMGF method of Theorem 1, by numerical integration of −∞ |x|r f (x)dx, and by simu-
lations, where we generate X1 , . . . , XN independent and identically NIG distributed, with
N =1,000,000 and take the average of |Xi |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 Chen and Revels (2016). Computation were done with
Julia v1.11.0, see Bezanson et al. (2017), on a MacBook Pro M1 Max with 32 GM memory.

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.

3 Applications to Dynamic Models


We consider three applications that demonstrate the usefulness of the new CMGF
expressions. In the first application, we use Theorems 1 and 3 to obtain moments
for cumulative returns from a Heston-Nandi GARCH model. The second application
uses Theorem 2 to compute fractional moments of realized volatility that follow an
Autoregressive Gamma process. The third application also applies Theorem 2 to obtain
fractional moments, but in this case for a discrete distribution that has probability mass
at zero.

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 .

3.1 Moments of Cumulative Returns in a GARCH Model


The Heston-Nandi GARCH (HNG) model by Heston and Nandi (2000) is given by
  q
1
rt+1 = rf + λ − 2
ht+1 + ht+1 zt+1 , with zt ∼ iidN (0, 1),
 q 2 (7)
ht+1 = ω + βht + α zt − γ ht ,

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.

Proposition 1. Let rt , t = 1, 2, . . ., be given by (7) and define cumulative returns,


RT,H = rT +1 + · · · + rt+H . Then the conditional MFG has the affine form,

MRT,H |T (z) = exp (A(H, z) + B(H, z)hT +1 ) , (8)

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).

Table 2: Heston-Nandi GARCH Estimation Results for Daily Log-returns of SPX

λ ω β α γ ℓ

1.9781 1.15 ×10−14 0.7593 5.67×10−6 185.5 17,963


(0.3166) (0.85×10−14 ) (0.0190) (1.18×10−6 ) (23.5)

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.

We can now illustrate Theorem 1 by computing the fractional absolute moments of


cumulative returns, RT,H using the MGF we derived in Proposition 1. Figure 3 presents
absolute moments, E|RT,H |r , for a range of r ∈ (−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 3 in this application, by computing the conditional
integer moments, ET [(RT,H )k ], k ∈ {1, 2, 3, 4}. From these moments, we compute

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 ,

ET [RT3 +H ] − µ3T,H µt,H


SkewT,H = 3
−3 ,
σT,H σT,H
ET [RT4 +H ] − µ4T,H µT,H µ2T,H
KurtT,H = 4
−4 SkewT,H −6 2 .
σT,H σT,H σT,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.

3.2 Moments in Autoregressive Gamma Model


In this application, Xt represents the daily realized variance, which is computed from
the intraday transaction data. We follow Corsi et al. (2013) and Majewski et al. (2015)
and adopt a Heterogeneous Autoregressive Gamma (HARG) model for Xt . The HARG

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

(d) (w) (m)


θt = βd Xt + βw Xt + βm Xt , (9)

(d) (w) (m)


= 41 5i=2 Xt−i , and Xt 1 22
P P
where the variables, Xt = Xt−1 , Xt = 17 i=6 Xt−i
represent lagged daily, “weekly”, and “monthly” averages. Observe that θt is Ft−1 -
measurable. The HARG model implies a restricted AR(22) structure,

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.

Proposition 2. Suppose Xt follows a HARG(p) process. Then the conditional MGF


for XT +H , given Ft , is given by
 
  p
MXT +H |T (z) ≡ Et ezXt+H = exp A(H, z) +
X
Bj (H, z)XT +1−j  ,
j=1

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

and the recursions,

A(h + 1, z) = A(h, z) − δ log[1 − ηB1 (h, z)],


B1 (h,z)
Bj (h + 1, z) = ϕ
1−ηB1 (h,z) j
+ 1{j<p} Bj+1 (h, z), j = 1, . . . , p.

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.

Table 3: HARG Estimation Results for Daily Realized Variance

β̃d β̃w β̃m η δ ℓ

0.4896 0.2789 0.0357 0.0053 0.9644 17,298


(0.0279) (0.0333) (0.0120) (0.0003) (0.0226)

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.

We adopt the realized (Parzen) kernel estimator, RKt by Barndorff-Nielsen et al.


(2008), as our realized measure of the variance. The daily RKt is estimated for the S&P
500 index over the sample period from January 1, 2000 to November 30, 2021 (5,490
trading days).7 The realized measures are annualized by the scaling, Xt = 252 × RKt .
Table 3 presents the maximum-likelihood estimates for the HARG models along with
their corresponding standard errors (in parentheses).
 
Under Proposition 2, we compute the term structure of moments Et XTr +H for
H = 1, . . . , 180 and r ∈ {− 21 , 12 , 32 , 2}. Figure 5 plots the term structure of these four
1
moments when the initial value of lagged Xt components are all set as 10 E[Xt ], where
E[Xt ] = ηδ/ (1 − ηβd − ηβw − ηβm ) . We include both the simulated value (solid line)
from one million Monte-Carlo simulations and the numerical value (dashed line) from
Proposition 2. We can find the numerical values fit the simulated values very well.
7
The data of realized kernel was obtained from the Realized Library at the Oxford-Man Institute,
which was discontinued in 2022.

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 ).

3.3 Moments in Autoregressive Poisson Model


The autoregressive Poisson (ARP) model is given by

λyt −λt
Pr(Yt = y|Ft−1 ) = e , y = 0, 1, . . . (12)
y!

where the dynamic intensity parameter evolves according to

λt+1 = ω + βλt + αYt , t = 1, 2, . . . , (13)

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(λ) ℓ

0.1548 0.7473 0.2043 3.2024 -9722


(0.0386) (0.0349) (0.0243)

Note: ML estimation results for the ARP model with robust standard errors in parentheses.

For H = 1, the conditional MGF is

MYT +1 |T (z) = exp (λT +1 (ez − 1)) . (14)

More generally, the analytical form of the conditional MGF of ȲT,H is given in the
following Proposition 3.

Proposition 3. Let Yt be given by (12) and (13). Then

MȲT,H |T (z) = exp (A(H, z) + B(H, z)λT +1 ) ,

where A(H, z) and B(H, z) are given from

A(h + 1, z) = A(h, z) + ωB(h, z),


B(h + 1, z) = βB(h, z) + (ez/H+αB(h,z) − 1),

with initial value A(1, z) = 0, and B(1, z) = ez/H − 1.

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.

Under Proposition 3, we compute the term structure of moments ET [(ȲT,H )r ] for


r ∈ { 12 , 1, 32 , 2} and H ranging from 1 day to six months. Figure 6 plots the term
1
structure of these four moments when the initial value of λt+1 set as 10 E (λt ). We
include both the simulated value (solid line) from one million Monte-Carlo simulations
and the numerical value (dashed line) from Proposition 3. We can find the numerical
values fit the simulated values very well.

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. (1978). Hyperbolic distributions and distributions on hyper-


bolae. Scandinavian Journal of Statistics, 5:151–157.

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.

Chen, J. and Revels, J. (2016). Robust benchmarking in noisy environments. arXiv


e-prints.

Corsi, F. (2009). A Simple Approximate Long-Memory Model of Realized Volatility.


Journal of Financial Econometrics, 7:174–196.

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.

Duffie, D. and Singleton, K. J. (1993). Simulated moments estimation of markov-models


of asset prices. Econometrica, 61:929–952.

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, L. P. (1982). Large Sample Properties of Generalized Method of Moments


Estimators. Econometrica, 50(4):1029–54.

Hansen, P. R. and Tong, C. (2024). A globally valid integral expression for the reciprocal
gamma function. Working paper.

Heston, S. L. and Nandi, S. (2000). A closed-form GARCH option valuation model.


Review of Financial Studies, 13:585–625.

Jones, M. (1987). Inverse factorial moments. Statistics & Probability Letters, 6:37–42.

Kawata, T. (1972). Fourier analysis in probability theory. Academic Press.

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.

McFadden, D. (1989). A method of simulated moments for estimation of discrete


response models without numerical integration. Econometrica, pages 995–1026.

Samko, S. G., Kilbas, A. A., and Marichev, O. I. (1993). Fractional integrals and
derivatives (theory and applications). Gordon and Breach, Switzerland.

Schürger, K. (2002). Laplace transforms and suprema of stochastic processes. In Ad-


vances in Finance and Stochastics, pages 285–294. Springer, Berlin, Heidelberg.

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.

Wolfe, S. (1975). On moments of probability distribution functions. In Fractional


Calculus and Its Applications, pages 306–316.

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 ).

Lemma A.1. Let r + 1 ∈ C+ and z ∈ C+ . Then


ˆ ∞
Γ(r + 1)
xr e−zx dx = .
0 z r+1

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

where we used the substitution: y = zx ∈ C. The function in the last integral,


f (y) ≡ y r e−y , has a pole at y = 0 if r < 0, and is a multi-valued function when r is a
non-integer. We will adopt the principal value by restricting the angle of y = λeiθ to
θ ∈ (−π, π], and construct a closed contour integral that avoids the pole at y = 0. The
contour is illustrated in Figure 7.

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

For the BC-integral, we have y = sK + iuK and dy = iKdu, such that


ˆ sK+itK ˆ t
|IBC | = r −y
y e dy = iK (sK + iuK)r e−(sK+iuK) du
sK 0
ˆ t
r+1 −sK
= iK e (s + iu)r e−iuK du
0
ˆ t
≤K Re(r)+1 −sK
e |(s + iu)r | du → 0, as K → ∞,
0

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

The integral we seek is therefore given by,


ˆ sK+itK ˆ sK
lim f (y)dy = lim f (y)dy = Γ(r + 1).
K→∞ 0 K→∞ 0

This completes the proof.

Lemma A.2. For x < 0, z = s + it ∈ C+ , and r + 1 ∈ C+ we have


ˆ +∞
ezx
dt = 0.
−∞ z r+1

For x = 0, the identity holds for r ∈ C+ .

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 −∞

because the last term is zero for r with Re (r) > 0.


Next, for x < 0 we have
ˆ +∞ ˆ +∞
ezx e(s+it)x
dt = r+1 dt
−∞ z r+1 −∞ (s + it)
ˆ −∞
e(s−iτ )x
=− r+1 dτ, τ = −t
+∞ (s − iτ )
ˆ +∞
e(s−iτ )x
= r+1 dτ
−∞ (s − iτ )
ˆ +∞
sx e−iτ x
=e r+1 dτ
−∞ (s − iτ )
ˆ +R
sx
= lim e f (τ )dτ,
R→∞ −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

This completes the proof.

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

This is the Laplace transform of xr and the inverse Laplace transform is


ˆ σ+i∞
r 1 Γ(r + 1) zx
x = e dz, x > 0,
2πi σ−i∞ 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

where the second term is zero from Lemma A.2.


Then, we focus on the case that x < 0. From the identity A.1, we directly have
ˆ +∞
r Γ(r + 1) r e−zx
|x| = (−x) = dt, x<0
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

where the second term is zero from Lemma A.2.


Combine the two formulas (A.2) and (A.3), we have the first formula in Lemma 1
ˆ +∞
r Γ(r + 1) ezx + e−zx
|x| = dt, x ∈ R.
2π −∞ z r+1

The second formula in Lemma 1, from expression (A.1), we have


ˆ
k Γ(k + 1) +∞ ezx
x = k+1
dt, x > 0,
2π −∞ z
ˆ +∞ −zx
k k Γ(k + 1) e
x = (−1) k+1
dt, x < 0,
2π −∞ z

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

where we use Γ(k + 1) = k! for k ∈ N. This completes the proof. □


Proof of Theorem 1. The proof is the same as in the proof the Theorem 2 by using
the first formula in Lemma 1. □
Proof of Theorem 2. We prove the result for discrete and continuous distributions
separately. Consider first a non-negative discrete random variable, X, with support X .
Rewriting E[X r ] = x∈X xr p(x) using the third formula in Lemma 1 yields
P

"ˆ ∞
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

for r + 1 ∈ C+ if Pr(X = 0) = 0 and r ∈ C+ if Pr(X = 0) > 0. The same line of

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.
−∞ −∞ −∞

This completes the proof. □


Proof of Theorem 3. The proof is the same as in the proof the Theorem 1 by using
the second formula in Lemma 1. □

B Result for the Dynamic Models


B.1 MGF for ARG model
Proof of Proposition 1. First, the dynamic process of θt can be rewritten as


 βd j=1
p


X 
θT = βj XT −j , βi = βw /4 2≤j≤5

j=1 

β /17 6 ≤ j ≤ 22

m

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

For the case, H = 1, we can use the MGF (11),


 
    p
ηz
ET ezXT +1 = exp δ log 1−ηz
X
1
+ θ
1−ηz T +1
= exp A(1, z) + Bj (1, z)XT +1−j  ,
j=1

for s < η −1 where

η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

Note that we need the s to be such that

B1 (h, s) < η −1 for ∀ h ≤ T − t.

That completes the proof of Proposition 2. □

Conditional MGF for Autoregressive Poisson


Proof of Proposition 2. The conditional MGF for a Poisson distributed random
 
variable, Yt , with intensity λt , is Et−1 ezYt = exp (λt (ez − 1)). So for the ARP process,

Yt |Ft−1 ∼ Poisson(λt ), with λt = ω + βλt−1 + αnt−1 .

The average number over the next H periods, is

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

MST,h+1 (z) = ET [exp(zST,h+1 )] = ET [exp(zYT +1 + zST +1,h )]


= ET [ET +1 {exp(zYT +1 + zST +1,h )}]
= ET [exp{zYT +1 + Ã(h, z) + B̃(h, z)λT +2 }]
= ET [exp{zYT +1 + Ã(h, z) + B̃(h, z)(ω + βλT +1 + αYT +1 )}]
= exp{Ã(h, z) + B̃(h, z)(ω + βλT +1 )}ET [exp{(z + αB̃(h, z))YT +1 )}]
= exp{Ã(h, z) + B̃(h, z)(ω + βλT +1 )} exp{λT +1 (ez+αB̃(h,z) − 1)}
= exp{Ã(h, z) + B̃(h, z)(ω + βλT +1 ) + λT +1 (ez+B̃(h,z)α − 1)}
= exp{[Ã(h, z) + ω B̃(h, z)] + [β B̃(h, z) + (ez+B̃(h,z)α − 1)]λT +1 }.
| {z } | {z }
=Ã(h+1,z) B̃(h+1,z)

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

= Et [exp (sRt+1 + A(h, s) + B(h, s)ht+2 )]


1
 
= exp A(h, s) + sr + B(h, s)ω − log (1 − 2B(h, s)α)
2
(s − 2B(h, s)αγ)2
! !
1
   
2
× exp s λ − + B(h, s) β + αγ + ht+1
2 2 (1 − 2B(h, s)α)
= exp (A(h + 1, s) + B(h + 1, s)ht+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.

Figure S.1: Integrands for some selected NIG moments.

Existing Methods for Fractional Moments


Existing expressions for positive fractional moments include,
ˆ ∞
" K
#
r −(1+r)
X u2k (2k)
E|X| = CK u −Re[φX (u)] + φX (0) du,
0 k=0 (2k)!

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 Laue (1980).


Closely relate to (S.1) is the following expression for non-negative random variables,
ˆ ∞ (k) (k)
r k λ MX (0) − MX (−u)
E[X ] = (−1) Γ(1−λ)
du, (X ≥ 0),
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

Moments in HNG Model


% MATLAB ( R2024a version )
function Y = S i m u l a t i o n M o m e n t s (N ,H ,r , model )
Omega = model (1) ;
Beta = model (2) ;
Alpha = model (3) ;
Gamma = model (4) ;
Lam = model (5) ;
rf = model (6) ;
h1 = model (7) ;
Z = randn (N , H ) ;
h = zeros (N , H ) ;
h (: ,1) = h1 ;
for t =1: H -1
h (: , t +1) = Omega + Beta * h (: , t ) + Alpha *( Z (: , t ) - Gamma * sqrt ( h (: , t ) ) ) .^2;
end
Rmat = H * rf +( Lam -1/2) * sum (h ,2) + sum ( sqrt ( h ) .* Z ,2) ;
Y = mean ( abs ( Rmat ) .^ r ) ;
end

function Y = CMGF (H ,r , model , Tol )


Y = gamma ( r +1) / pi * integral ( @ ( s ) IntFunction (s ,H ,r , model ) ,0 , Inf , ’ AbsTol ’ ,Tol , ’ ArrayValued ’ , false )
;
end

function fs = IntFunction (s ,H ,r , model )


Omega = model (1) ;
Beta = model (2) ;

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

You might also like