An algorithm for the computation of joint Hawkes moments with exponential kernel

The purpose of this paper is to present a recursive algorithm and its implementation in Maple and Mathematica for the computation of joint moments and cumulants of Hawkes processes with exponential kernels. Numerical results and computation times are also discussed. Obtaining closed form expressions can be computationally intensive, as joint fifth cumulant and moment formulas can be respectively expanded into up to 3,288 and 27,116 summands.

READ FULL TEXT VIEW PDF

Authors

page 1

page 2

page 3

page 4

03/13/2018

On moments of exponential functionals of additive processes

Let X = (X t) t>0 be a real-valued additive process, i.e., a process wit...
03/01/2018

Fast and accurate computation of orthogonal moments for texture analysis

In this work we propose a fast and stable algorithm for the computation ...
12/14/2020

Recursive computation of the Hawkes cumulants

We propose a recursive method for the computation of the cumulants of se...
07/02/2021

Moments of Subsets of General Equiangular Tight Frames

This note outlines the steps for proving that the moments of a randomly-...
06/28/2021

Exact simulation of extrinsic stress-release processes

We present a new and straightforward algorithm that simulates exact samp...
03/13/2018

On moments of integral exponential functionals of additive processes

Let X = (X t) t>0 be a real-valued additive process, i.e., a process wit...
06/26/2019

Correlators of Polynomial Processes

A process is polynomial if its extended generator maps any polynomial to...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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.

method=escape,ActualText= 1 c := proc(n, mu) local tmp, k, z1; option remember; if n = 1 then return 1/(1 - mu); end if;
method=escape,ActualText= 2   tmp := 0; z1 := []; for k from n by -1 to 2 do z1 := [op(z1), c(n - k + 1, mu)]; tmp := tmp + IncompleteBellB(n, k, op(z1)); end do;
method=escape,ActualText= 3   return mu*tmp/(1 - mu); end proc;
method=escape,ActualText= 4 m := proc(n, mu) local tmp, z, k; option remember; if n = 0 then return 1; end if;
method=escape,ActualText= 5   tmp := 0; z := []; for k from n by -1 to 1 do z := [op(z), c(n - k + 1, mu)]; tmp := tmp + IncompleteBellB(n, k, op(z)); end do;
method=escape,ActualText= 6   return tmp; end proc;

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.

method=escape,ActualText= 1 c[n_, mu_] := c[n, mu] = (Module[{tmp, k}, If[n == 1, Return[1/(1 - mu)]]; tmp = 0; z1 = {};
method=escape,ActualText= 2   For[k = n, k >= 2, k--, z1 = Append[z1, Block[{i = n - k + 1}, c[i, mu]]]; tmp += BellY[n, k, z1]];
method=escape,ActualText= 3   Simplify[mu*tmp/(1 - mu)]]);
method=escape,ActualText= 4 m[n_, mu_] := (Module[{tmp, z, k}, tmp = 0; If[n == 0, Return[1]];
method=escape,ActualText= 5    z = {}; For[k = n, k >= 1, k--, z = Append[z, Block[{i = n - k + 1}, c[i, mu]]]; tmp += BellY[n, k, z]]; Simplify[tmp]])

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.

Figure 1: Sample paths of and of the intensity .

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 .

The joint cumulants of can be obtained as a consequence of Proposition 3.1, by the combinatorial summation

(5)

see Corollary 3.4 and Proposition 3.5 in [Pri21]. Joint moments can then be recovered by the joint moment-cumulant relation

(6)

which can be inverted as

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 .

Figure 2: 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.

Figure 3: Second joint moment.

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

Figure 4: Third joint moments.
Figure 5: Fourth joint moments.

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
Table 1: Counts of summands and cumulant computation times in Maple.

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.

Figures 6 and 7 show the numerical evaluation of fifth and sixth joint moments. and .

Figure 6: Fifth joint moments.
Figure 7: Sixth joint moments.

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
Table 2: Counts of summands and moment computation times in Maple.

One-variable examples are computed using the following codes which use Bell polynomials instead of sums over partitions.

method=escape,ActualText= 1 kz := proc(z, n, a, b, t) local tmp, z1, k; option remember;
method=escape,ActualText= 2   if n = 1 then if a = b then return 1 + a*(t - z); else return b/(b - a) + a*exp((a - b)*(t - z))/(a - b); end if; end if;
method=escape,ActualText= 3   tmp := 0; z1 := []; for k from n by -1 to 2 do z1 := [op(z1), kz(y + z, n - k + 1, a, b, t)]; tmp := tmp + IncompleteBellB(n, k, op(z1)); end do;
method=escape,ActualText= 4   return int(a*exp((a - b)*y)*tmp, y = 0 .. t - z); end proc;
method=escape,ActualText= 5 c := proc(a, b, t, n) local y, k, z, temp; option remember; temp := 0; y := [];
method=escape,ActualText= 6   for k from n by -1 to 1 do y := [op(y), kz(z, n - k + 1, a, b, t)]; temp := temp + IncompleteBellB(n, k, op(y)); end do;
method=escape,ActualText= 7   return int(temp, z = 0 .. t); end proc;
method=escape,ActualText= 8 m := proc(a, b, t, n) local tmp, z, k; option remember;
method=escape,ActualText= 9 tmp := 0; if n = 0 then return 1; end if; z := []; for k from n by -1 to 1 do z := [op(z), c(a, b, t, n - k + 1)]; tmp := tmp + IncompleteBellB(n, k, op(z)); end do;
method=escape,ActualText= 10   return tmp; end proc;

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.

method=escape,ActualText= 1 kz[z_, n_, a_, b_, t_] :=
method=escape,ActualText= 2   kz[z, n, a, b, t] = (Module[{tmp, k, y, z1},
method=escape,ActualText= 3      If[n == 1, If[a === b, Return[1 + a*(t - z)],
method=escape,ActualText= 4        Return[b/(b - a) + a*E^((a - b)*(t - z))/(a - b)]]]; tmp = 0;
method=escape,ActualText= 5      z1 = {}; For[k = n, k >= 2, k--, z1 = Append[z1, Block[{i = n - k + 1, u = y + z}, kz[u, i, a, b, t]]];
method=escape,ActualText= 6       tmp += BellY[n, k, z1];]; a*Integrate[E^((a - b)*y)*tmp, {y, 0, t - z}]]);
method=escape,ActualText= 7 c[a_, b_, t_, n_] := c[a, b, t, n] = (Module[{y, k, z, temp}, temp = 0; y = {}; For[k = n, k >= 1, k--,
method=escape,ActualText= 8       y = Append[y, Block[{i = n - k + 1, u = z}, kz[u, i, a, b, t]]]; temp += BellY[n, k, y]]; Return[Integrate[temp, {z, 0, t}]]]);
method=escape,ActualText= 9 m[a_, b_, t_, n_] := m[a, b, t, n] = (Module[{tmp, z, k}, tmp = 0; If[n == 0, Return[1]]; z = {}; For[k = n, k >= 1, k--, z = Append[z, c[a, b, t, n - k + 1]]; tmp += BellY[n, k, z]]; tmp])

Appendix A - Maples codes

method=escape,ActualText= 1 kz := proc(z, a, b, t::list) local pm, pp2, p, pp, tmp, k, y, h, i, j, ii, u, n, zz, c; option remember; n := nops(t);
method=escape,ActualText= 2   if n = 1 then if a = b then return 1 + a*(t[1] - z); else return b/(b - a) + a*exp((a - b)*(t[1] - z))/(a - b); end if; end if;
method=escape,ActualText= 3   tmp := 0; pp2 := Iterator:-SetPartitions(n); for pp in pp2 do p := pp2:-ToSets(pp);
method=escape,ActualText= 4   if 2 <= nops(p) then c := 1; for i to nops(p) do c := c*kz(y + z, a, b, map(op, convert(p[i], list), t)); end do;
method=escape,ActualText= 5   tmp := tmp + c; end if; end do; return a*int(exp((a - b)*y)*tmp, y = 0 .. t[1] - z); end proc;
method=escape,ActualText= 1 c := proc(a, b, t::list) option remember; local y, e, k, pm, tmp, p2, pp, p, c, i, zz, j, u, ii, n; n := nops(t);
method=escape,ActualText= 2   tmp := kz(y, a, b, t); if 2 <= n then pm := Iterator:-SetPartitions(n); for pp in pm do p := pm:-ToSets(pp);
method=escape,ActualText= 3   if 2 <= nops(p) then e := 1; for i to nops(p) do e := e*kz(y, a, b, map(op, convert(p[i], list), t)); end do; tmp := tmp + e;
method=escape,ActualText= 4   end if; end do; end if; return int(tmp, y = 0 .. t[1]); end proc;
method=escape,ActualText= 1 m := proc(a, b, t::list) option remember; local y, e, k, u, ii, pm, tmp, p2, pp, p, i, zz, j, n; n := nops(t); tmp := c(a, b, t);
method=escape,ActualText= 2   if 2 <= n then pm := Iterator:-SetPartitions(n); for pp in pm do p := pm:-ToSets(pp); if 2 <= nops(p) then e := 1; for i to nops(p) do e := e*c(a, b, map(op, convert(p[i], list), t)); end do; tmp := tmp + e;
method=escape,ActualText= 3   end if; end do; end if; return tmp; end proc;

Appendix B - Mathematica codes

method=escape,ActualText= 1 Needs["Combinatorica`"];
method=escape,ActualText= 2 kz[z_, a_, b_, t__] :=
method=escape,ActualText= 3   kz[z, a, b, t] = (Module[{tmp, y, i, c, n}, n = Length[t];
method=escape,ActualText= 4      If[n == 1, If[a === b, Return[1 + a*(t[[1]] - z)],
method=escape,ActualText= 5        Return[b/(b - a) + a*E^((a - b)*(t[[1]] - z))/(a - b)]]];
method=escape,ActualText= 6      tmp = 0; Do[c = 1; If[Length[p] >= 2, For[i = 1, i <= Length[p], i++,
method=escape,ActualText= 7         c = c*Block[{u = y + z, v = t[[p[[i]]]]}, kz[u, a, b, v]]];
method=escape,ActualText= 8        tmp += c], {p, SetPartitions[n]}];
method=escape,ActualText= 9      Return[a*Integrate[E^((a - b)*y)*tmp, {y, 0, t[[1]] - z}]]]);
method=escape,ActualText= 1 c[a_, b_, t__] :=
method=escape,ActualText= 2   c[a, b, t] = (Module[{y, e, tmp, n, i}, n = Length[t]; tmp = 0;
method=escape,ActualText= 3      Do[e = 1; For[i = 1, i <= Length[p], i++,
method=escape,ActualText= 4        e = e*Block[{u = y, v = t[[p[[i]]]]}, kz[u, a, b, v]]];
method=escape,ActualText= 5        tmp += Flatten[{e}][[1]],
method=escape,ActualText= 6        {p, SetPartitions[n]}];
method=escape,ActualText= 7      Return[Integrate[tmp, {y, 0, t[[1]]}]]]);
method=escape,ActualText= 1 m[a_, b_, t__] := (Module[{n, e, i, tmp}, tmp = 0; n = Length[t]; If[n == 0, Return[1]];
method=escape,ActualText= 2    Do[e = 1; For[i = 1, i <= Length[p], i++, e = e*c[a, b, t[[p[[i]]]]]];
method=escape,ActualText= 3     tmp += e, {p, SetPartitions[n]}]; Flatten[{tmp}][[1]]])

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