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

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 $[0, \infty)$ .

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.

		$M_{X} (t_{1}, \dots, t_{n}) := E [e^{t_{1} X_{1} + \dots + t_{n} X_{n}}]$
		$= 1 + \sum k_{1}, \dots, k_{n} \geq 1 \frac{t_{1}^{k_{1}} \dots t_{n}^{k_{n}}}{n!} E [X_{1}^{k_{1}} \dots X_{n}^{k_{n}}],$

of a random vector

X = (X_{1}, \dots, X_{n})

, the joint cumulants of

(X_{1}, \dots, X_{n})

of orders

(l_{1}, \dots, l_{n})

are the coefficients

κ_{l_{1}, \dots, l_{n}} (X)

appearing in the log-MGF expansion

		$log M_{X} (t_{1}, \dots, t_{n}) = log (E [e^{t_{1} X_{1} + \dots + t_{n} X_{n}}])$		(1)
		$= \sum l_{1}, \dots, l_{n} \geq 1 \frac{t_{1}^{l_{1}} \dots t_{n}^{l_{n}}}{l_{1}! \dots l_{n}!} κ_{l_{1}, \dots, l_{n}} (X_{1}, \dots, X_{n}),$

for $(t_{1}, \dots, t_{n})$ in a neighborhood of zero in $R^{n}$ . In the sequel we let

κ (X_{1}, \dots, X_{n}) := κ_{1, \dots, 1} (X_{1}, \dots, X_{n}), n \geq 1,

and

κ^{(n)} (X) := κ_{1, \dots, 1} (X, \dots, X), n \geq 1.

The joint moments of $(X_{1}, \dots, X_{n})$ are then given by the joint moment-cumulant relation

E [X_{1} \dots X_{n}] = n \sum l = 1 \sum π_{1} \cup \dots \cup π_{l} = {1, \dots, n} l \prod j = 1 κ^{(| π_{j} |)} ((X_{i})_{i \in π_{j}}) .

where the sum runs over the partitions $π_{1}, \dots, π_{k}$ of the set ${1, \dots, n}$ , By the multivariate Faà di Bruno formula, (1) can be inverted as

	$κ (X_{1}, \dots, X_{n})$
	$= n \sum l = 1 (l - 1)! (- 1)^{l - 1} \sum π_{1} \cup \dots \cup π_{l} = {1, \dots, n} l \prod j = 1 E [\prod i \in π_{j} X_{i}] .$

In the univariate case, the moments $E [X^{n}]$ of a random variable $X$ are linked to its cumulants $(κ^{(n)} (X))_{n \geq 1}$ through the relation

	$E [X^{n}]$	$= B_{n} (κ^{(1)} (X), \dots, κ^{(n)} (X))$
		$= n \sum k = 1 B_{n, k} (κ^{(1)} (X) \dots κ^{(n - k + 1)} (X)),$

where

		$B_{n, k} (a_{1}, \dots, a_{n - k + 1}) = \frac{n!}{k!} \sum \frac{l_{1} + \dots + l_{k} = n}{l_{1} \geq 1, \dots, l_{k} \geq 1} \frac{a_{l_{1}}}{l_{1}!} \dots \frac{a_{l_{k}}}{l_{k}!}$
		$= \sum π_{1} \cup \dots \cup π_{k} = {1, \dots, n} a_{\| π_{1} \|} (X) \dots a_{\| π_{k} \|} (X),$		(2)

$1 \leq k \leq n$ , is the partial Bell polynomial of order $(n, k)$ , where the sum (2) holds on the integer compositions $(l_{1}, \dots, l_{k})$ of $n$ , see e.g. Relation (2.5) in [Luk55], and

B_{n} (a_{1}, \dots, a_{n})

= n \sum k = 1 B_{n, k} (a_{1}, \dots, a_{n - k + 1})

is the complete Bell polynomial of degree $n \geq 1$ . We also have the inversion relation

		$κ^{(n)} (X)$

$n \geq 1$ , 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 $(X_{n})_{n \geq 0}$ be a branching process started at $X_{0} = 1$

with Poisson distributed offspring count

N

of parameter

μ \in (0, 1)

, and let

X

denote the total count of offsprings generated by

(X_{n})_{n \geq 0}

It is known, see [PS98] and § 3.2 of [CF06] that

X

has the Borel distribution Consul, P.C. Famoye, F. Pólya, G. Szegö, G.

P (X = n) = e^{- μ n} \frac{(μ n)^{n - 1}}{n!}, n \geq 1.

We have $κ^{(1)} (X) = 1 / (1 - μ)$ and the induction relation

		$κ^{(n)} (X) = \frac{μ}{1 - μ} (B_{n} (κ^{(1)} (X), \dots, κ^{(n)} (X)) - κ^{(n)} (X))$
		$= \frac{μ}{1 - μ} n \sum k = 2 B_{n, k} (κ^{(1)} (X), \dots, κ^{(n - k + 1)} (X)),$		(3)

$n \geq 2$ , 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 $c (2, μ)$ in Maple yields the second cumulant $κ^{(2)} (X) = μ / (1 - μ)^{3}$ , and by the commands $c (3, μ)$ and $c (4, μ)$ we find $κ^{(3)} (X) = μ (1 + 2 μ) / (1 - μ)^{5}$ , and $κ^{(4)} (X) = μ (1 + 8 μ + 6 μ^{2}) / (1 - μ)^{7}$ , see also (8.85) page 159 of [CF06]. Those results can be recovered from the command $c [2, μ]$ , $c [3, μ]$ and $c [4, μ]$ 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 $[0, \infty)$ , with Poisson offspring intensity $γ (d x)$ and Poisson immigrant intensity $ν (d x)$ on $[0, \infty)$ , built on the space

	$Ω = {ξ = {x_{i}}_{i \in I} \subset [0, \infty) :$

of locally finite configurations on $[0, \infty)$ , whose elements $ξ \in Ω$ are identified with the Radon point measures $ξ (d z) = \sum x \in ξ ϵ_{x} (d z)$ , where $ϵ_{x}$ denotes the Dirac measure at $x \in R_{+}$ . In particular, any initial immigrant point $y \in R_{+}$ branches into a Poisson random sample denoted by $ξ_{γ} (y + d z) = \sum x \in ξ ϵ_{x + y} (d z)$ and centered at $y$ , with intensity measure $γ (y + d z)$ on $[0, \infty)$ . Figure 1 presents a graph of the point measure $ξ (d z)$ followed by the corresponding sample paths of the self-exciting counting process $X_{t} (ξ) := ξ ([0, t]) = \sum x \in ξ 1_{[0, t]} (x)$ and its stochastic intensity $λ_{t}$ , $t \in [0, 10]$ , in the exponential kernel example of the next section.

Figure 1: Sample paths of $X_{t}$ and of the intensity $λ (t)$ .

In the sequel, we assume that $γ ([0, \infty)) < 1$ and consider the integral operator $Γ$ defined as

Γ f (z) = \int_{0}^{\infty} f (z + y) γ (d y), z \in R_{+},

and the inverse operator $(I_{d} - Γ)^{- 1}$ given by

	$(I_{d} - Γ)^{- 1} f (z) = f (z)$
	$+ \infty \sum m = 1 \int_{R_{+}^{m}} f (z + x_{1} + \dots + x_{m}) γ (d x_{1}) \dots γ (d x_{m}),$

with

	$(I_{d} - Γ)^{- 1} Γ f (z) = (I_{d} - Γ)^{- 1} f (z) - f (z)$
	$= \infty \sum m = 1 \int_{R_{+}^{m}} f (z + x_{1} + \dots + x_{m}) γ (d x_{1}) \dots γ (d x_{m}),$

$z \in R_{+}$ . The first cumulant $κ_{z}^{(1)} (f)$ of $\sum x \in ξ f (x)$ given that $ξ$ is started from a single point at $z \in R_{+}$ is given by $κ_{z}^{(1)} (f) = (I_{d} - Γ)^{- 1} f (z)$ for $n = 1$ . The next proposition provides a way to compute the higher order cumulants $κ_{z}^{(n)} (f)$ of $\sum x \in ξ f (x)$ given that $ξ$ is started from a single point at $z \in R_{+}$ by an induction relation based on set partitions, see Proposition 3.5 in [Pri21].

Proposition 3.1

For $n \geq 2$ , the joint cumulants $κ_{z}^{(n)} (f_{1}, \dots, f_{n})$ of $\sum x \in ξ f_{1} (x), \dots, \sum x \in ξ f_{n} (x)$ given that $ξ$ is started from a single point at $z \in R_{+}$ are given by the induction relation

		$κ_{z}^{(n)} (f_{1}, \dots, f_{n})$		(4)
		$= n \sum k = 2 \sum π_{1} \cup \dots \cup π_{k} = {1, \dots, n} (I_{d} - Γ)^{- 1} Γ κ_{z}^{(\| π_{j} \|)} ((f_{i})_{i \in π_{j}}),$

$n \geq 2$ , where the above sum is over set partitions $π_{1} \cup \dots \cup π_{k} = {1, \dots, n}$ , $k = 2, \dots, n$ , and $| π_{i} |$ denotes the cardinality of the set $π_{i} \subset {1, \dots, n}$ .

The joint cumulants $κ^{(n)} (f_{1}, \dots, f_{n})$ of $\sum x \in ξ f_{1} (x), \dots, \sum x \in ξ f_{n} (x)$ can be obtained as a consequence of Proposition 3.1, by the combinatorial summation

		$κ^{(n)} (f_{1}, \dots, f_{n})$		(5)
		$= n \sum k = 1 \sum π_{1} \cup \dots \cup π_{k} = {1, \dots, n} \int_{0}^{\infty} k \prod j = 1 κ_{z}^{(\| π_{j} \|)} ((f_{i})_{i \in π_{j}}) ν (d z),$

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

		$E [\sum x \in ξ f_{1} (x) \dots \sum x \in ξ f_{n} (x)]$		(6)
		$= n \sum l = 1 \sum π_{1} \cup \dots \cup π_{l} = {1, \dots, n} l \prod j = 1 κ^{(\| π_{j} \|)} ((f_{i})_{i \in π_{j}}),$

which can be inverted as

	$κ^{(n)} (f_{1}, \dots, f_{n})$
	$= n \sum l = 1 (l - 1)! (- 1)^{l - 1} \sum π_{1} \cup \dots \cup π_{l} = {1, \dots, n} l \prod j = 1 E [\prod i \in π_{j} \sum x \in ξ f_{i} (x)] .$

4 Joint Hawkes moments with exponential kernel

In this section we consider the exponential kernel $γ (d x) = a 1_{[0, \infty)} (x) e^{- b x} d x$ , $0 < a < b$ , and constant Poisson intensity $ν (d z) = ν d z$ , $ν > 0$ . In this case,

X_{t} (ξ) := ξ ([0, t]) = \sum x \in ξ 1_{[0, t]} (x), t \in R_{+},

defines the self-exciting Hawkes process with stochastic intensity

λ_{t} := ν + a \int_{0}^{t} e^{- b (t - s)} d X_{s}, t \in R_{+} .

In this case, the integral operator $Γ$ satisfies

Γ f (z) = a \int_{0}^{\infty} f (z + y) e^{- b y} d y, z \in R_{+},

and the recursive calculation of joint moments and cumulants will be performed by evaluating $(I_{d} - Γ)^{- 1} Γ$ in Proposition 3.1 on the family of functions $e_{p, η, t}$ of the form $e_{p, η, t} (x) := x^{p} e^{η x} 1_{[0, t]} (x)$ , $η < b$ , $p \geq 0$ , as in the next lemma.

Lemma 4.1

For $f$ in the linear span generated by the functions $e_{p, η, t}$ , $p \geq 0$ , $η \in R$ , the operator $(I_{d} - Γ)^{- 1} Γ$ is given by

(I_{d} - Γ)^{- 1} Γ f (z) = a \int_{0}^{t - z} f (z + y) e^{(a - b) y} d y,

$z \in [0, t]$ .

Proof. For all $p, η \geq 0$ we have the equality

	$(I_{d} - Γ)^{- 1} Γ e_{p, η, t} (z)$
	$= \infty \sum n = 1 \int_{[0, t]^{n}} e_{p, η, t} (z + x_{1} + \dots + x_{n}) γ (d x_{1}) \dots γ (d x_{n})$
	$= \infty \sum n = 1 \frac{a^{n}}{(n - 1)!} \int_{0}^{t - z} (z + y)^{p} e^{η (z + y)} y^{n - 1} e^{- b y} d y$
	$= a e^{η z} \int_{0}^{t - z} (z + y)^{p} e^{(η + a - b) y} d y, z \in [0, t],$

which follows from the fact that the sum $τ_{1} + \dots + τ_{n}$ of $n$ exponential random variables with parameter $b > 0$

has a gamma distribution with shape parameter

n \geq 1

and scaling parameter

b > 0

□

Using Lemma 4.1, we can rewrite (4) for $t_{1} < \dots < t_{n}$ as

		$κ_{z}^{(n)} (1_{[0, t_{1}]}, \dots, 1_{[0, t_{n}]})$		(7)
		$= n \sum k = 2 \sum π_{1} \cup \dots \cup π_{k} = {1, \dots, n} \int_{0}^{t_{1} - z} a e^{(a - b) y} κ_{z + y}^{(\| π_{j} \|)} ((1_{[0, t_{i}]})_{i \in π_{j}}) d y,$

with

	$κ_{z}^{(1)} (1_{[0, t]})$	$= (I_{d} - Γ)^{- 1} 1_{[0, t]} (z)$
		$= \frac{b}{b - a} + \frac{a}{a - b} e^{(a - b) (t - z)}, z \in [0, t],$

if $a \neq b$ , and

κ_{z}^{(1)} (1_{[0, t]}) = (I_{d} - Γ)^{- 1} 1_{[0, t]} (z) = 1 + a (t - z),

$z \in [0, t]$ , if $a = b$ . The recursive computation of $κ_{z}^{(n)} (1_{[0, t_{1}]}, \dots, 1_{[0, t_{n}]})$ 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 $E [X_{t_{1}} \dots X_{t_{n}}]$ of $X_{t_{1}}, \dots, X_{t_{n}}$ are obtained from the above code using the command $m (a, b, [t_{1}, \dots, t_{n}])$ in Maple or $m [a, b, {t_{1}, \dots, t_{n}}]$ in Mathematica. Figures 2 to 7 are plotted with $ν = 1$ , $a = 0.5$ , $b = 1$ , $T = 2$ , and one million Monte Carlo samples, and Figure 2 presents the first moment $m_{1} (t) = m (a, b, [t]) = m [a, b, {t}]$ .

For example, the second joint moment of $(X_{t}, X_{T})$ is obtained by the command $m_{2} (t, T) = m (a, b, [t, T])$ in Maple or $m_{2} (t, T) = m [a, b, {t, T}]$ in Mathematica presented in appendix, which yields

	$E [X_{t} X_{T}] =$
	$\frac{1}{(a - b)^{4}} (\frac{1}{2} e^{- a t - b (t + T)} (2 b^{4} t e^{a t + b (t + T)}$
	$+ a^{2} b (- e^{a (2 t + T)} + e^{2 b t + a T} + 2 b t e^{2 a t + b T}$

	$- e^{b t + a (t + T)} (1 + b t)) - 2 a b^{2} (e^{2 b t + a T} - 2 e^{2 a t + b T}$
	$- e^{b t + a (t + T)} + e^{a t + b (t + T)} (2 + b t)))$
	$+ (b^{2} t + a (e^{(a - b) t} - b t - 1))$
	$\times (b^{2} T + a (e^{(a - b) T} - b T - 1))),$

see Figure 3.

Figures 4 and 5 show the numerical evaluation of third and fourth joint moments, obtained from $m_{3} (t_{1}, t_{2}, t_{3}) = m (a, b, [t_{1}, t_{2}, t_{3}]) = m [a, b, {t_{1}, t_{2}, t_{3}}]$ and $m_{4} (t_{1}, t_{2}, t_{3}, t_{4}) = m (a, b, [t_{1}, t_{2}, t_{3}, t_{4}]) = m [a, b, {t_{1}, t_{2}, t_{3}, t_{4}}]$ .

The following tables count the (approximate) numbers of summands appearing in joint cumulant and moment expressions when expanded as a sum of the form

\sum \frac{k, l, p_{1}, \dots, p_{n} \geq 0}{q_{1}, \dots, q_{n}, r_{1}, \dots, r_{n} \geq 0} a^{k} b^{l} t_{1}^{p_{1}} \dots t_{n}^{p_{n}} e^{q_{1} a t_{1} + \dots + q_{n} a t_{n} + r_{1} b t_{1} + \dots + r_{n} b t_{n}},

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. $m_{5} (t_{1}, t_{2}, t_{3}, t_{4}, t_{5})$ and $m_{6} (t_{1}, t_{2}, t_{3}, t_{4}, t_{5}, t_{6})$ .

Computation times are presented in seconds for symbolic calculations using the variables $a, b, t_{1}, \dots, t_{n}$ , 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 $X_{t}$ is obtained in Maple by the command $m (a, b, t, 2)$ , which yields

	$E [X_{t}^{2}] =$	$\frac{(b^{2} t + a (- 1 + e^{(a - b) t} - b t))^{2}}{(a - b)^{4}}$
		$+ \frac{1}{2 (a - b)^{4}} (e^{- 2 b t} (2 b^{4} e^{2 b t} t$
		$+ a^{2} b (- e^{2 a t} + e^{2 b t} + 4 b e^{(a + b) t} t)$
		$- 2 a b^{2} (- 3 e^{(a + b) t} + e^{2 b t} (3 + b t))$
		$+ 2 a^{3} (e^{2 a t} - e^{(a + b) t} (1 + 2 b t)))) .$

The same result can be obtained in Mathematica from the command $m [a, b, t, 2]$ 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.</

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

Authors

On moments of exponential functionals of additive processes

Fast and accurate computation of orthogonal moments for texture analysis

Recursive computation of the Hawkes cumulants

Moments of Subsets of General Equiangular Tight Frames

Exact simulation of extrinsic stress-release processes

On moments of integral exponential functionals of additive processes

Correlators of Polynomial Processes

1 Introduction

2 Joint moments and cumulants

3 Joint Hawkes cumulants

Proposition 3.1

4 Joint Hawkes moments with exponential kernel

Lemma 4.1

Appendix A - Maples codes

Appendix B - Mathematica codes

References

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

Authors

Related Research

On moments of exponential functionals of additive processes

Fast and accurate computation of orthogonal moments for texture analysis

Recursive computation of the Hawkes cumulants

Moments of Subsets of General Equiangular Tight Frames

Exact simulation of extrinsic stress-release processes

On moments of integral exponential functionals of additive processes

Correlators of Polynomial Processes

This week in AI

1 Introduction

2 Joint moments and cumulants

3 Joint Hawkes cumulants

Proposition 3.1

4 Joint Hawkes moments with exponential kernel

Lemma 4.1

Appendix A - Maples codes

Appendix B - Mathematica codes

References