1 Introduction
Hawkes processes [Haw71] are self-exciting point processes that have found applications in fields such as neuroscience [CR10], genomics analysis [RBS10], as well as finance [ELL11], or social media [RLMX18].
The analysis of statistical properties of Hawkes processes is made difficult by their recursive nature, making the computation of moments difficult. In [JHR15]
, a tree-based method for the computation of cumulants has been introduced, with an explicit computation of third order cumulants. Ordinary differential equation (ODE) methods have been applied in
[DZ11]to the computation of the moment and probability generating functions of (generalized) Hawkes processes and their intensity, with the computation of first and second moments in the stationary case, see also
[EGG10], and [CHY20] and [DP20] for other ODE-based approaches.In [BDM12], stochastic calculus and martingale arguments have been applied to the computation of first and second order moments, however those approaches seem difficult to generalize to higher-orders moments. In [LRV20]
, cumulant recursion formulas have been obtained for general random variables using martingale brackets. Third-order cumulant expressions for Hawkes processes have been used in
[ABMR18] for the the analysis of order books in finance, and in [OJSBB17], [MMG20]for neuronal networks.
In [Pri21], the cumulants of Hawkes processes have been computed using using Bell polynomials, based on a recursive relation for the Probability Generating Functional (PGFl) of self-exciting point processes started from a single point. This provides a closed-form alternative to the tree-based approach of [JHR15].
In this note we apply the algorithm of [Pri21] to the recursive computation of joint moments of all orders of Hawkes processes, and present the corresponding codes written in Maple and Mathematica. The algorithm uses sums over partitions and Bell polynomials to compute joint cumulants in the case of an exponential branching intensity on .
We proceed as follows. After reviewing some combinatorial identities in Section 2, we will consider the computation of the joint cumulants of self-exciting Hawkes Poisson cluster processes in Section 3. Explicit computations for the time-dependent joint third and fourth cumulants of Hawkes processes with exponential kernels are presented in Section 4
, and are confirmed by Monte Carlo estimates.
2 Joint moments and cumulants
In this section we present background combinatorial results that will be needed in the sequel. Given the Moment Generating Function (MGF) Thiele, T.N.
of a random vector
, the joint cumulants of of orders are the coefficients appearing in the log-MGF expansion(1) | ||||
for in a neighborhood of zero in . In the sequel we let
and
The joint moments of are then given by the joint moment-cumulant relation
where the sum runs over the partitions of the set , By the multivariate Faà di Bruno formula, (1) can be inverted as
In the univariate case, the moments of a random variable are linked to its cumulants through the relation
where
(2) |
, is the partial Bell polynomial of order , where the sum (2) holds on the integer compositions of , see e.g. Relation (2.5) in [Luk55], and
is the complete Bell polynomial of degree . We also have the inversion relation
, see e.g. Theorem 1 of [Luk55], and also [LS59], Relations (2.8)-(2.9) in [McC87], or Corollary 5.1.6 in [Sta99].
As an example we consider the recursive computation of Borel cumulants as an example. Let be a branching process started at
with Poisson distributed offspring count
of parameter , and let denote the total count of offsprings generated by It is known, see [PS98] and § 3.2 of [CF06] that has the Borel distribution Consul, P.C. Famoye, F. Pólya, G. Szegö, G.We have and the induction relation
(3) |
, see § 8.4.3 in [CF06] and Proposition 2.1 in [Pri21]. The recursion (3) is implemented in the following Maple code.
In particular, the command in Maple yields the second cumulant , and by the commands and we find , and , see also (8.85) page 159 of [CF06]. Those results can be recovered from the command , and using the following Mathematica code.
3 Joint Hawkes cumulants
In the cluster process framework of [HO74], we consider a real-valued self-exciting point process on , with Poisson offspring intensity and Poisson immigrant intensity on , built on the space
of locally finite configurations on , whose elements are identified with the Radon point measures , where denotes the Dirac measure at . In particular, any initial immigrant point branches into a Poisson random sample denoted by and centered at , with intensity measure on . Figure 1 presents a graph of the point measure followed by the corresponding sample paths of the self-exciting counting process and its stochastic intensity , , in the exponential kernel example of the next section.

In the sequel, we assume that and consider the integral operator defined as
and the inverse operator given by
with
. The first cumulant of given that is started from a single point at is given by for . The next proposition provides a way to compute the higher order cumulants of given that is started from a single point at by an induction relation based on set partitions, see Proposition 3.5 in [Pri21].
Proposition 3.1
For , the joint cumulants of given that is started from a single point at are given by the induction relation
(4) | ||||
, where the above sum is over set partitions , , and denotes the cardinality of the set .
4 Joint Hawkes moments with exponential kernel
In this section we consider the exponential kernel , , and constant Poisson intensity , . In this case,
defines the self-exciting Hawkes process with stochastic intensity
In this case, the integral operator satisfies
and the recursive calculation of joint moments and cumulants will be performed by evaluating in Proposition 3.1 on the family of functions of the form , , , as in the next lemma.
Lemma 4.1
For in the linear span generated by the functions , , , the operator is given by
.
Proof. For all we have the equality
which follows from the fact that the sum of exponential random variables with parameter
has a gamma distribution with shape parameter
and scaling parameter .Using Lemma 4.1, we can rewrite (4) for as
(7) | ||||
with
if , and
, if . The recursive computation of in (7) is implemented in the following Mathematica code using Lemma 4.1.
The computation of joint Hawkes cumulants by the recursive relation (5) is then implemented in the following code.
Finally, joint moments are computed from the joint moment-cumulant relation (6) which is implemented in the following code. The joint moments of are obtained from the above code using the command in Maple or in Mathematica. Figures 2 to 7 are plotted with , , , , and one million Monte Carlo samples, and Figure 2 presents the first moment .

For example, the second joint moment of is obtained by the command in Maple or in Mathematica presented in appendix, which yields
see Figure 3.

Figures 4 and 5 show the numerical evaluation of third and fourth joint moments, obtained from and .


The following tables count the (approximate) numbers of summands appearing in joint cumulant and moment expressions when expanded as a sum of the form
excluding factorizations and simplifications of expressions.
One variable | ||||
Cumulant | Time | Count | All variables | |
Sixth | 64s | 671 | Time | Count |
Fifth | 11s | 226 | 2403s | 3288 |
Fourth | 1.7s | 81 | 31s | 536 |
Third | 0.5s | 35 | 1.6s | 91 |
Second | 0.2s | 12 | 0.3s | 14 |
First | 0.06s | 4 |
The tables also display the corresponding computation times on a 8-core laptop computer with 8Gb RAM. Symbolic computation appears faster with Maple, although computation times become similar at the order six and above.


Computation times are presented in seconds for symbolic calculations using the variables , and can be significantly shorter when the variables are set to specific numerical values. Moment computations times in Table 2 are similar to those of Table 1, and can be sped up if cumulant functions are memoized after repeated calls.
One variable | ||||
Moment | Time | Count | All variables | |
Sixth | 66s | 2159 | Time | Count |
Fifth | 11s | 762 | 2544s | 27116 |
Fourth | 1.9s | 265 | 35s | 2236 |
Third | 0.5s | 88 | 1.8s | 266 |
Second | 0.2s | 22 | 0.5s | 29 |
First | 0.05s | 4 |
One-variable examples are computed using the following codes which use Bell polynomials instead of sums over partitions.
For example, the second moment of is obtained in Maple by the command , which yields
The same result can be obtained in Mathematica from the command using the code below.
Appendix A - Maples codes
Appendix B - Mathematica codes
References
- [ABMR18] M. Achab, E. Bacry, J. F. Muzy, and M. Rambaldi. Analysis of order book flows using a non-parametric estimation of the branching ratio matrix. Quant. Finance, 18(2):199–212, 2018.
- [BDM12] E. Bacry, K. Dayri, and J.F. Muzy. Non-parametric kernel estimation for symmetric Hawkes processes. Application to high frequency financial data. Eur. Phys. J. B, 85:157–168, 2012.
-
[CF06]
P.C. Consul and F. Famoye.
Lagrangian probability distributions
. Birkhäuser Boston, Inc., Boston, MA, 2006. - [CHY20] L. Cui, A. Hawkes, and H. Yi. An elementary derivation of moments of Hawkes processes. Adv. in Appl. Probab., 52:102–137, 2020.
- [CR10] S. Cardanobile and S. Rotter. Multiplicatively interacting point processes and applications to neural modeling. Journal of Computational Neuroscience, 28:267–284, 2010.
- [DP20] A. Daw and J. Pender. Matrix calculations for moments of Markov processes. Preprint arXiv:1909.03320, 2020.
- [DZ11] A. Dassios and H. Zhao. A dynamic contagion process. Adv. in Appl. Probab.