# Functional Laplace Transform of a Multivariate Hawkes Process, Subsequent Characteristics, and Numerical Approximations

Bartholomé Vieille, Rachid Senoussi, Samuel Soubeyrand

► To cite this version:

Bartholomé Vieille, Rachid Senoussi, Samuel Soubeyrand. Functional Laplace Transform of a Multivariate Hawkes Process, Subsequent Characteristics, and Numerical Approximations. 2025. hal-05168648

HAL Id: hal-05168648

https://hal.inrae.fr/hal-05168648v1

Preprint submitted on 18 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Functional Laplace Transform of a Multivariate Hawkes Process, Subsequent Characteristics, and Numerical Approximations

Bartholomé Vieille*†

Rachid Senoussi*‡

Samuel Soubeyrand*§¶

July 17, 2025

**Abstract.** Numerous studies grounded on Hawkes processes have been carried out in many fields including finance, biology and social network. Hawkes processes form a class of self-exciting simple point processes. In this article, we consider a general class of multivariate Hawkes processes envisioned to model dynamics of spatio-temporal epidemics. For this class, the igniting baseline intensity is time dependent and the exciting matrix function is a general one, making the model non-Markovian in most of the cases. In this article, we first provide the closed-form expression of the multivariate multi-temporal characteristic function of these Hawkes processes, extending in a natural way the classical single-time formula found in the Hawkes literature. Then, we use the infinitely divisible property of the Hawkes process to derive the equation system related to the probability distribution of counts at each single time, adapted to the general formulation of the Hawkes model considered in this article. Next, we provide closed-form formulas for the temporal structure of the two first moments of the process, which allows us to deduce an original expression of the multivariate covariance function at two distinct times, thereby extending existing results established for more restricted classes of Hawkes processes. Based on this expression, we analytically decompose the covariance at two distinct times into singular and continuous parts. We finish with brief numerical elements: We present a simple scheme for numerical approximations of the Laplace transform and the first two moments, and give examples of solutions of the different related integral equations. We also provides illustrative simulations of the multivariate Hawkes process for different model specifications.

**Keywords.** Hawkes process; Laplace transform; self-exciting process; spatio-temporal process; infinite divisibility property; covariance structure.

**MSC2020 subject classifications.** 60G55; 60E10.

<footnote>*INRAE, BioSP, 84914 Avignon, France</footnote>

<footnote>†bartholv@gmail.com</footnote>

<footnote>‡rachid.senoussi@inrae.fr</footnote>

<footnote>§samuel.soubeyrand@inrae.fr</footnote>

<footnote>¶https://orcid.org/0000-0003-2447-3067</footnote>
# 1 Introduction

A univariate simple counting process $N$ is defined by:

$$N(t) = \sum_{l=1}^{+\infty} 1_{\{T_l \le t\}},$$

where $T_1 < \dots < T_l < \dots$ is a sequence of positive increasing random variables, and $1_E$ equals 1 if event $E$ is true and 0 otherwise.

As usual, to describe the stochastic structure and the dynamics of $N$, we consider a complete probability space $(\Omega, \mathcal{F}, \mathbf{P})$ with the complete right-continuous natural filtration $\mathcal{F} = (\mathcal{F}_t)_{t \ge 0}$, with $\mathcal{F}_t = \sigma(N(s), s \le t)$.

$(N(t))_{t \ge 0}$ is said a *Hawkes process* driven by $(\lambda^0, \phi)$ if its $(\mathcal{F} - \mathbf{P})$ previsible compensator admits an intensity $\tilde{\lambda}^0(t)$ of the following form:

$$\tilde{\lambda}^0(t) = \lambda^0(t) + \sum_{l=1}^{+\infty} \phi(t - T_l) 1_{\{T_l \le t\}}. \qquad (1)$$

where $\lambda^0$, called the baseline intensity, and the function $\phi$, called the excitation function, are deterministic non-negative real functions defined on $\mathbb{R}_+$. This process is said self-exciting since each event $T_l$ induces a new point process with baseline intensity function $\phi$.

Hawkes processes were introduced by Alan G. Hawkes [16]. They are a generalization of Poisson processes, wherein multiple Poisson processes are added over time. Hawkes processes have been used in various fields such as finance [4, 17], genomics [32], epidemiology [34], insurance [37], neuroscience [15, 27], and seismology [5].

From the intensity expression (1), it is evident that Hawkes processes are usually non-Markovian and therefore challenging to investigate. Consequently, a significant portion of research on Hawkes processes focuses on specific models. For instance, extensive studies have been conducted on the case where the baseline intensity $\lambda^0$ is constant and the exciting function $\phi$ is exponential, in which case the analysis simplifies into the analysis of a two-dimensional Markovian process $(\tilde{\lambda}^0, N)$; see [14, 16, 30, 35].

The integral formula for the probability generating function restricted to a single time for the one-dimensional case was provided by [18]. This result was further extended to Hawkes processes with random jump sizes, as established by [13] who derived the expression for the Laplace Transform for the one-dimensional case. Moreover, [12] investigated the Laplace and Fourier transforms and the distribution of a Hawkes process in a specific case where the exciting function is exponential and with a particular baseline intensity. Furthermore, using martingale methods, [22] provided second-, third- and fourth-order moments at multiple times for one-dimensional Hawkes processes with constant baseline intensity, these results being used in the study of weak convergence of Hawkes processes toward Cox-Ingersoll-Ross models [23]; See also [20], who employed the Poisson embedding representation, Malliavin calculus, and the pseudo-chaotic expansion of the Hawkes process to derive expressions for second-order moments in the one-dimensional case with constant baseline intensity, as above.

A multivariate version of a Hawkes process was proposed by [9], and referred to as a marked temporal process $N$ with a finite mark set denoted by $\{1, \dots, d\}$. Subsequently, [11] obtained an integral formula for the characteristic function of the multivariate case, limited to a single time $t$; see also [26]. The extension of these formulas to multiple times is crucial for statistical inference if data include multiple observation times and for identifying process limits. Precisely, using a population representation, [24] obtained simple expressions of second- and third-order moments at multiple times for multivariate Hawkes processes that are stationary over $\mathbb{R}$, i.e., grounded on a constant baseline intensity and a short-time dependence assumption on the excitation
function. [31] provided implicit forms of functional moments generalizing the results of [24] to stationary spatial Hawkes processes. These forms are derived from the expression of the probability generating functional (Proposition 3.1 in [31]), which is analogous to our Theorem 2.3 obtained in the non-stationary case.

In this study, we consider a more general class of non-stationary, multidimensional Hawkes processes, characterized by a flexible and relatively broad excitation function. We characterize the multivariate multi-temporal Laplace transform and moments using integral equations. To achieve this, we use the immigrant-birth representation of a linear Hawkes Process, as described in [18] and [25] for marked Hawkes processes. This result not only characterizes the process but also provides intrinsic information about its behavior over time.

In Section 2, we present a Volterra-like system of equations expressing the multi-temporal characteristic function of our multivariate Hawkes process. A similar expression for the functional Laplace Transform of the process is also given, extending classical single-time expressions. In Section 3, we detail the precise form of the infinitely divisible property for our Hawkes process and the integral equations related to the probability distribution of the Hawkes process at any time. Section 4 is dedicated to expressing and calculating the first two moments of the process as well as the form of certain statistical functionals related to the general theory of spatial statistics. We obtain an original expression for the covariance function at two distinct times. Then, we propose in Section 5 some approaches to tackle the numerical aspects for the resolution of the system of Volterra-like equations mentioned above. We illustrate these numerical approaches for a specific Hawkes model, focusing on the mean and covariance functions. Section 6 gives a few simulated examples of the multi-dimensional Hawkes process to illustrate the effects of the excitation function and the interaction between the different dimensions (i.e., regions in a spatial setting). To conclude, Section 7 outlines the general framework of this work and mentions some avenues of potential applications for the results that we obtained, especially with regard to the development of statistical inference methods.

## Notations.

First, we convene that $j \in [r:k]$ holds for the enumeration $j = r, \dots, k$. In accordance with matrix calculus, a $d$-vector $w$ is always considered as a column vector while its transpose is the row vector $^\top w = (w_1, \dots, w_d)$, and for a $d \times n$-matrix $m = (m_i^j: i = 1: d, j = 1: n)$, its $j^{th}$ column is denoted $m_i^j$ and its $i^{th}$ row as $m_i^j$. We also set $\bar{w} = \sum_{i=1}^d w_i$ and adopt the useful notation of restriction $u_{[r:k]} = (u_r, \dots, u_k)$ for an $n$-row vector $u$ as well as $m_i^{[r:k]} = (m_i^j: 1 \le i \le d, r \le j \le k)$ for matrix $m$. We also convene that a $k$-vector having identical components, typically a scalar $u$, is written $u_{[k]}$, e.g. $0_{[d]}$ denotes the vector of zeros with dimension $d$. The left limit of a function $f$ is defined by $f^{-}(t) = \lim_{s \uparrow t} f(s)$ if it exists (as this is the case for non decreasing functions), and its $\tau$-translate is defined by $f_{\tau_{+}}(s) = f(\tau + s)$. Multiple summation over $q = (q_1, \dots, q_d) \in \mathbb{N}^d$ is simply written $\sum_q$. The unit vectors of $\mathbb{R}^d$ are denoted $e_j, j = 1: d$.

# 2 Multivariate Hawkes Process

A $d$-dimensional Hawkes process driven by $(\lambda^0, \phi)$ is a simple point process with $d$ marks

$$N(t) = ^\top(N_1, \dots, N_d)(t)$$

described by the collection of the occurrence times $\{T_l^j: j = 1: d, l \ge 1\}$ (indicating that an $l^{th}$ event of type $j$ occurs at $T_l^j$) and such that considering its natural filtration, $N$ has
a predictable compensator $\tilde{\Lambda}^0(t) = \top(\tilde{\Lambda}_1^0, \dots, \tilde{\Lambda}_d^0)(t)$ admitting a conditional intensity function $\tilde{\lambda}^0(t) = \top(\tilde{\lambda}_1^0, \dots, \tilde{\lambda}_d^0)(t)$ of the following form:

$$ \tilde{\lambda}_i^0(t) = \lambda_i^0(t) + \sum_{j=1}^{d} \sum_{l:T_l^j < t} \phi_i^j(t - T_l^j), \quad (2) $$

where $\lambda^0 : \mathbb{R}_+ \to (\mathbb{R}_+)^d$ and $\phi = (\phi_i^j) : \mathbb{R}_+ \to \mathcal{M}_d^+$, being the set of $d \times d$ matrices with non-negative coefficients.

This can be rewritten in integral-matrix forms as:

$$ \tilde{\lambda}^0(t) = \lambda^0(t) + \int_0^{t-} \phi(t-u) \, dN(u). \quad (3) $$

$N$ can represent various counting contexts involving interacting entities, such as populations (cities, countries, animal species, cells), finance and economic activities, or physical fundamental particles, types, or locations of earthquakes, etc. For the sake of clarity, in epidemiology, we assume that it describes the evolution of an infection in a spatial domain composed of $d$ distinct regions. Therefore, $N_i(t)$ represents the number of infected individuals in region $i$ at time $t$.

Equation (2) describes the dynamics of the process, indicating the presence of an initial unknown ancestor (labeled as $j = 0$) that triggers a $d$-dimensional Poisson process $N^0 = \top(N_1^0, \dots, N_d^0)$ with independent components and conditional intensity $\lambda^0 = \top(\lambda_1^0, \dots, \lambda_d^0)$. This process is described by occurrence times $T_l^{0,j}$, where $j = 1 : d$ and $l \ge 1$, representing individuals of the first generation. Each individual of the first generation of type $j$, occurring at time $T_l^{0,j}$, in turn, triggers a new $d$-dimensional Poisson process with conditional intensity $\phi_j^0 = \top(\phi_1^j, \dots, \phi_d^j)$, resulting in a second generation of individuals ($l' \ge 1$) of type $j'$ at times $(T_{l,l'}^{1,j,j'}} : j' = 1 : d, l' \ge 1$). This process continues recursively, generating individuals of the third generation and so on. The occurrence times $\{T_{l,l'}^{m,j,j'} : j, j' \in 1 : d; m \ge 0; l, l' \ge 1\}$ when ordered increasingly in a unique sequence $\{T_n^j : j = 1, \dots, d; n \ge 1\}$, actually form the Hawkes process $N$.

**Remark and Notation.** It is well known for multivariate point processes, see [21] for example, that when considering the natural filtration, the knowledge of the conditional intensity, parameterized here by $(\lambda^0, \phi)$, characterizes entirely the probability distribution of the process itself. Accordingly, when clarification is needed, we enhance this by using subscripts, eg $P_{(\lambda^0, \phi)}$, $E_{(\lambda^0, \phi)}$. Also, by replacing $\lambda^0$ with $\phi_j^{j'}$, $j' = 1 : d$, we can define $d$-new Hawkes processes $\tilde{N}^{j'}$, each directed by $(\phi_j^{j'}, \phi)$. These processes will play an important role in the sequel. Therefore convening that $\phi_j^0 = \lambda^0$ enables us to identify $N$ with $\tilde{N}^0$ and then to group together many formulas.

## 2.1 Population Representation of Multivariate Hawkes Processes

A useful representation of $N$ for analytic investigation is the so called population approach adopted in previous works like [18] and [11]. We mainly use here their notations to describe a Hawkes process $\tilde{N}^0(t)$, directed by $(\lambda^0, \phi)$, as the sum of an ancestor Poisson process $N^0$ and of all Hawkes processes ignited by each first generation individual. That is to say, at each event $T_l^{0,j}$, a new Hawkes process $\tilde{N}^{j,(l)}$ directed by $(\phi_j^j, \phi)$ starts (independently of index $l$), such that we have the representation

$$ \tilde{N}^0(t) = N^0(t) + \sum_{j=1}^{d} \sum_{l=1}^{N^0(t)} \tilde{N}^{j,(l)}(t - T_l^{0,j}). \quad (4) $$
Important remarks for the sequel

1. Actually from a point of view of filtration, the representation (4) is much richer than the proper filtration of $\tilde{N}^0$, since we cannot recover actually the ascendant of individuals (i.e., retrieve the $T_l^{m,j}$ from the knowledge of the $T_n^j$).

2. Conditionally to the knowledge of the sequence $T_l^{m,j}$, the Hawkes processes are independent of each others. Moreover for each index $j$, the $\tilde{N}^{j,(l)}$ are identically distributed, say as a reference Hawkes point process $\tilde{N}^j$.

3. Note that the sum of independent Hawkes processes $\tilde{N}^{0,(l)}$, $l = 1, \dots, L$, each directed by $(\lambda^{0,(l)}, \phi)$, is also a Hawkes process directed by $(\sum_{l=1}^L \lambda^{0,(l)}, \phi)$. Consequently, taking $\lambda^{0,(l)} = \lambda^0/L$, proves that Hawkes processes have infinitely divisible probability distributions.

Similarly, any subdivision of $\mathbb{R}_+ = \cup_{l=1}^L [t_{l-1}, t_l[$ shows that $\tilde{N}^0$ directed by $(\lambda^0, \phi)$ is the sum of $L$ independent copies of Hawkes processes directed by $(\lambda^0 \chi_{[t_{l-1},t_l[}, \phi), l = 1, \dots, L$, where $\chi_{[t_{l-1},t_m[}$ is the indicator function that gives 1 if its argument is in $[t_{l-1}, t_l[$ and 0 otherwise, that is $\chi_{[t_{l-1},t_l[}(t) = \mathbf{1}_{t \in [t_{l-1},t_l[}$.

Probabilities of Zero and One Events

As a consequence of the population representation (4), we have the following result:

**Lemma 2.1.** For any $t \ge 0$ and $j' = 0 : d$, we have for $0_{[d]} = (0, \dots, 0) \in \mathbb{R}^d$

$$ \tilde{p}_0^{j'}(t) = P(\tilde{N}^{j'}(t) = 0_{[d]}) = \exp \left( - \sum_{j=1}^{d} \int_{0}^{t} \phi_j^{j'}(s) ds \right) \quad (5) $$

$$ \tilde{p}_{e_i}^{j'}(t) = P(\tilde{N}^{j'}(t) = e_i) = \left( \int_0^t \phi_i^{j'}(s) ds \right) \exp \left( - \sum_{j=1}^d \int_0^t \phi_j^{j'}(s) ds \right). \quad (6) $$

*Proof.* For $j' = 0$, Equation (4) tells us that the events "$\tilde{N}_i^0(t) = 0$" (resp. "$\tilde{N}_i^0(t) = 1$") correspond exactly to the independent Poissonian events "$N_i^0(t) = 0$" (resp. "$N_i^0(t) = 0$") directed by $\phi_i^0 = \lambda_i^0$, for all $i = 1 : d$. Clearly, this also holds true for $j' = 1 : d$. $\square$

## 2.2 Multi-Time Functional Laplace Transform of a Multivariate Hawkes Process

Let $a = (a_i^j)$ be a $d \times n$ matrix with non negative values and $t = (t_1, \dots, t_n)$ a set of ordered times $0 \le t_1 \le \dots \le t_n$.

The $n$-multi-temporal Laplace transform of the $d$-multivariate point process $\tilde{N}^{j'}$, directed by $(\phi_{\bullet}^{j'}, \phi)$, $j' = 0 : d$, is defined as:

$$ \tilde{L}^{j', (n)}(a, t) := \mathbf{E}_{(\phi_{\bullet}^{j'}, \phi)} \left[ \exp \left( - \sum_{r=1}^{n} \top a_{\bullet}^{r} \tilde{N}^{j'}(t_r) \right) \right]. $$

### 2.2.1 Classical Laplace Transforms of Hawkes Processes

Let $\tilde{N}^{j'}$, $j' = 0 : d$, be Hawkes processes respectively driven by functional parameters $(\phi_{\bullet}^{j'}, \phi)$ with respective Laplace transforms $\tilde{L}^{j', (k)}$, $j' = 0 : d, k = 1 : n$. Then we have the following set of integral equations relating their multi-temporal Laplace transforms:
**Theorem 2.2.** With previous notations $t_{[k,m]}$, $a_{\bullet}^{[k,m]}$, $\overline{a_j^{[k,m]}} = \sum_{r=k}^m a_j^r$ and putting $t_0 = 0$, we have for every $j' = 0 : d$ and $m = 1 : n$,

$$ \tilde{L}^{j', (m)}(a, t) = \tilde{p}_0^{j'}(t_m) \prod_{k=1}^{m} \exp \left( \sum_{j=1}^{d} \int_{t_{k-1}}^{t_k} e^{-\overline{a_j^{[k,m]}}} \tilde{L}^{j, (m+1-k)}(a_{\bullet}^{[k,m]}, (t-u)_{[k,m]}) \phi_j^{j'}(u) du \right). \quad (7) $$

**Remark.** This theorem mainly tells us that the Laplace transform related to any $\lambda^0$ can be derived through the solution of the system of $d$ basic Laplace transforms $\tilde{L}^{j', (k)}$, $j' = 1 : d$. The system of non-linear integral equations (7) will yield Volterra equations for the first and second order moments of the process (see Section 4). Hence, we call this system a Volterra-like system of equations.

*Proof.* We thoroughly use the representation (4) and the associated remarks, where $N^0$ is a multivariate Poisson process with independent components. For sake of simplicity we denote for $j = 1 : d$, the set of occurrence times $\mathcal{T}_k^j = \{T_l^{0,j} : t_{k-1} \le T_l^{0,j} < t_k\}$ and the corresponding set of indices $E_k^j = \{l : T_l^{0,j} \in \mathcal{T}_k^j\}$ both of random size $Q_k^j$.

For each $m = 1, \dots, n$, the expansion

$$ \tilde{N}^0(t_m) = \sum_{k=1}^{m} \left( [N^0(t_k) - N^0(t_{k-1})] + \sum_{j,l \in E_k^j} \tilde{N}^{j,(l)}(t_m - T_l^{0,j}) \right), $$

shows that $\tilde{N}^0(t_m)$ is the sum of the contributions of $m$ independent Hawkes processes $\tilde{N}^{0,(k)}(t_m)$, each directed by $(\lambda^0 \chi_{[t_{k-1}, t_k[}, \phi)$.

Therefore, inverting summation order and regrouping similar terms yields:

$$ \sum_{m=1}^{n} \top a_{\bullet}^{m} \tilde{N}^{0}(t_{m}) = \sum_{k=1}^{n} \left( \top \overline{a_{\bullet}^{[k,n]}}[N^{0}(t_{k}) - N^{0}(t_{k-1})] + \sum_{m=k}^{n} \sum_{j,l \in E_{k}^{j}} \top a_{\bullet}^{m} \tilde{N}^{j,(l)}(t_{m} - T_{l}^{0,j}) \right), $$

which reads as the sum of the contributions of $n$ independent Hawkes processes $\tilde{N}^{0,(k)}$.

Next, for the Poisson process $N^0$ and $q = \top(q_1, \dots, q_d) \in \mathbb{N}^d$, we recall that

$$ \begin{aligned} p_k^0(q) &= \mathbf{P}(N^0(t_k) - N^0(t_{k-1}) = q) \\ &= \prod_{j=1}^{d} \exp \left( \Lambda_j^0(t_{k-1}) - \Lambda_j^0(t_k) \right) \left( \Lambda_j^0(t_k) - \Lambda_j^0(t_{k-1}) \right)^{q_j} / q_j!. \end{aligned} $$

We also recall that for each $k$, conditionally to the event $\{N^0(t_k) - N^0(t_{k-1}) = q\}$, the $d$ sets of occurrence times $\mathcal{T}_k^j$, $j = 1 : d$, are independent, that for each given $j$, the set $\mathcal{T}_k^j = \{T_l^{0,j} : t_{k-1} \le T_l^{0,j} < t_k\}$ of size $Q_k^j$ is a sequence of i.i.d. random variables having probability density $\lambda_j^0(u) / (\Lambda_j^0(t_k) - \Lambda_j^0(t_{k-1}))$ with $u \in [t_{k-1}, t_k[ [10]$. Consequently the set $\{\tilde{N}^{j,(l)}(\cdot - T_l^{0,j}) : t_{k-1} \le T_l^{0,j} < t_k\}$ of size $Q_k^j$ also forms a sequence of i.i.d. Hawkes processes.

So, we have
$$
\begin{align*}
\tilde{L}^{0,(n)}(a,t) &= \prod_{k=1}^{n} \left( \sum_{q} p_k^0(q) \prod_{j=1}^{d} e^{-a_j^{[k,n]}} q_j \mathbf{E} \left[ \prod_{l \in E_k^j} e^{(-\sum_{m=k}^n \top a_m^m \tilde{N}^{j,(l)} (t_m - T_l^{0,j}))} |Q_k^j = q_j, \mathcal{T}_k^j \right] \right) \\
&= \prod_{k=1}^{n} \left( \sum_{q} p_k^0(q) \prod_{j=1}^{d} e^{-a_j^{[k,n]}} q_j \mathbf{E} \left[ \prod_{l \in E_k^j} \tilde{L}^{j,(n+1-k)} (a_j^{[k,n]}, (t-T_l^{0,j})_{[k,n]}) |Q_k^j = q_j \right] \right) \\
&= \prod_{k=1}^{n} \left( \sum_{q} p_k^0(q) \prod_{j=1}^{d} e^{-a_j^{[k,n]}} q_j \left( \int_{t_{k-1}}^{t_k} \tilde{L}^{j,(n+1-k)} (a_j^{[k,n]}, (t-T_l^{0,j})_{[k,n]}) \frac{\lambda_j^0(u) du}{\Lambda^0(t_k) - \Lambda(t_{k-1})} \right)^{q_j} \right) \\
&= \prod_{k=1}^{n} \exp \left( \sum_{j=1}^{d} \int_{t_{k-1}}^{t_k} \left( e^{-a_j^{[k,n]}} \tilde{L}^{j,(n+1-k)} (a_j^{[k,n]}, (t-u)_{[k,n]}) - 1 \right) \lambda_j^0(u) du \right) \\
&= \exp \left( -\sum_{j=1}^{d} \int_{0}^{t_n} \lambda_j^0(u) du \right) \prod_{k=1}^{n} \exp \left( \sum_{j=1}^{d} \int_{t_{k-1}}^{t_k} e^{-a_j^{[k,n]}} \tilde{L}^{j,(n+1-k)} (a_j^{[k,n]}, (t-u)_{[k,n]}) \lambda_j^0(u) du \right).
\end{align*}
$$

We notice that the computation of any characteristic function $\tilde{L}^{0,(n)}$ is uniquely determined by the knowledge of the Laplace transforms $\tilde{L}^{j',(k)}$, $j' = 1: d$, $k = 1: n$. Consequently, replacing $\lambda^0$ with $\phi_j^{j'}$, $j' = 1: d$, yields a fundamental set of $d$ dependent integral equations:

$$
\tilde{L}^{j', (n)}(a,t) = \tilde{p}_0^{j'}(t_n) \prod_{k=1}^{n} \exp \left( \sum_{j=1}^{d} \int_{t_{k-1}}^{t_k} e^{-a_j^{[k,n]}} \tilde{L}^{j,(n+1-k)}(t-u)_{[k,n]} \phi_j^{j'}(u) du \right). \quad (8)
$$

<figure><img src="image_2.png" /></figure>

Equation (8) actually yields a recursive formula of Laplace transforms.

**Corollary 2.2.1.** Enhancing the dependence of Laplace transforms with respect to $\lambda^0$, we get

$$
\tilde{L}_{\lambda^0}^{0,(n)}(a,t) = \tilde{L}_{\lambda^0 \chi_{[0,t_n]}}^{0,(n)}(a,t) = \prod_{k=1}^{n} \tilde{L}_{\lambda^0 \chi_{[t_{k-1},t_k]}}^{0,(n+1-k)}(a_{\bullet}^{[k,n]}, (t-t_{k-1})_{[k,n]}).
$$

The same recurrence holds for the equations satisfied by $\tilde{L}_{\lambda^0}^{j', (n)}(a,t)$, $j' = 1: d$.

2.3 Functional Laplace Transforms

Let $t \in \mathbb{R}$ and $\psi: \mathbb{R}_+ \to \mathbb{R}_+^d$ be a bounded measurable vector function. We now state a somewhat more general theorem concerning the functional Laplace transform of multivariate Hawkes processes $\tilde{N}^{j'}$ driven by $(\phi_j^{j'}, \phi)$, $j' = 0: d$, where the functional Laplace transform is defined as:

$$
\tilde{S}^{j'}(\psi, t) := \mathbf{E}_{(\phi_j^{j'}, \phi)} \left[ \exp \left( - \int_0^t \top \psi(s) \tilde{N}^{j'}(ds) \right) \right].
$$

We have the following result, whose proof is very similar to that of Theorem 2.2
**Theorem 2.3.** For $t \ge 0$ and any bounded measurable $\psi : \mathbb{R}_+ \to \mathbb{R}_+^d$, we have for $j' = 0 : d$,

$$ \tilde{S}^{j'}(\psi, t) = \tilde{p}_0^{j'}(t) \exp \left( - \sum_{j=1}^{d} \int_{0}^{t} e^{\psi_j(s)} \tilde{S}^{j}(\psi_{s+}, t-s) \phi_j^{j'}(s) ds \right). $$

*Proof.* We start dealing with the case $j' = 0$ driven by $(\lambda^0, \phi)$. Once again, like in the proof of Theorem 2.2, let us denote for $j = 1 : d$ the set of occurrence times $\mathcal{T}^j = \{T_l^{0,j} : T_l^{0,j} \le t\}$ and its associate set of indices $E^j = \{l : T_l^{0,j} \in \mathcal{T}^j\}$ both of random size $N_j^0(t)$.

From the process representation (4), we have

$$ \begin{aligned} \int_0^t \top \psi(s) \tilde{N}^0(ds) &= \int_0^t \top \psi(s) N^0(s) + \sum_{j=1}^d \sum_{l=1}^{N_j^0(t)} \left( \int_0^{t-T_l^{0,j}} \top \psi(T_l^{0,j} + s) \tilde{N}^{j,l}(ds) \right) \\ &= \sum_{j=1}^d \sum_{l=1}^{N_j^0(t)} \left( \psi_j(T_l^{0,j}) + \int_0^{t-T_l^{0,j}} \top \psi(T_l^{0,j} + s) \tilde{N}^{j,l}(ds) \right). \end{aligned} $$

Hence, considering the independence of the point processes $N_j^0$, $j = 1 : d$, and that conditionally to $N_j^0(t) = n_j^0$, the set $\mathcal{T}^j$ forms an i.i.d. sample of size $n_j^0$ of random variables with probability distribution function $\lambda_j^0(s)/\Lambda_j^0(t)$ concentrated on the interval $[0, t]$ for each $j = 1, \dots, d$, entails that the set of random variables

$$ \{ \tilde{S}^j(\psi_{T_l^{0,j,+}}, t - T_l^{0,j}) = E_{(\phi_j^j, \phi)} \left( \exp \left( - \int_0^{t-T_l^{0,j}} \top \psi(T_l^{0,j} + s) \tilde{N}^{j,(l)}(ds) \right) \mid l = 1, \dots, n_j^0 \} $$

constitutes also a sample of i.i.d. random variables.

Therefore, we have

$$ \begin{aligned} \tilde{S}^0(\psi, t) &= \mathbf{E}_{(\lambda^0, \phi)} \left[ \exp \left( - \sum_{j=1}^d \sum_{l=1}^{N_j^0(t)} \left( \psi_j(T_l^{0,j}) + \int_0^{t-T_l^{0,j}} \top \psi(T_l^{0,j} + s) \tilde{N}^{j,l}(ds) \right) \right) \right] \\ &= \prod_{j=1}^d \mathbf{E}_{(\lambda^0, \phi)} \left( \mathbf{E} \left( \prod_{l=1}^{N_j^0(t)} e^{-\psi_j(T_l^{0,j})} \mathbf{E}_{(\phi_j^j, \phi)} \left( \exp \left( - \int_0^{t-T_l^{0,j}} \top \psi(T_l^{0,j} + s) \tilde{N}^{j,l}(ds) \right) \mid N_j^0(t) \right) \right) \right) \\ &= \prod_{j=1}^d \mathbf{E}_{(\lambda^0, \phi)} \left( \left[ \mathbf{E} \left( e^{-\psi_j(T_l^{0,j})} \tilde{S}^j(\psi_{T_l^{0,j,+}}, t - T_l^{0,j}) \right) \right]^{N_j^0(t)} \right) \\ &= \prod_{j=1}^d \mathbf{E}_{(\lambda^0, \phi)} \left( \left[ \int_0^t e^{-\psi_j(s)} \tilde{S}^j(\psi_{s+}, t-s) \frac{\lambda_j^0(s)}{\int_0^t \lambda_j^0(v) dv} ds \right]^{N_j^0(t)} \right). \end{aligned} $$

Since $N_j^0(t)$ are Poisson distributed with respective mean parameter $\Lambda_j^0(t) = \int_0^t \lambda_j^0(v) dv$, we finally get

$$ \begin{aligned} \tilde{S}^0(\psi, t) &= \prod_{j=1}^d \left( \sum_{k=0}^{\infty} \frac{e^{-\Lambda_j^0(t)}}{k!} (\Lambda_j^0(t))^k \left[ \int_0^t e^{-\psi_j(s)} \tilde{S}^j(\psi_{s+}, t-s) \frac{\lambda_j^0(s) ds}{\Lambda_j^0(t)} \right]^k \right) \\ &= e^{-\sum_{j=1}^d \Lambda_j^0(t)} \exp \left( \sum_{j=1}^d \int_0^t e^{-\psi_j(s)} \tilde{S}^j(\psi_{s+}, t-s) \lambda_j^0(s) ds \right). \end{aligned} $$
Similarly, replacing $\lambda^0$ with $\phi_j^{j'}$; $j' = 1:d$, yields

$$ \tilde{S}^{j'}(\psi, t) = e^{-\sum_{j=1}^{d} \int_{0}^{t} \phi_j^{j'}(s) ds} \exp \left( \sum_{j=1}^{d} \int_{0}^{t} e^{-\psi_j(s)} \tilde{S}^j(\psi_{s+}, t-s) \phi_j^{j'}(s) ds \right). $$

□

## Remarks

1. Actually Theorem 2.2 can be deduced from Theorem 2.3 by considering the following piece-wise constant function:

$$ \psi(s) = \sum_{k=1}^{n} \overline{a^{[k,n]} 1_{[t_{k-1}, t_k]}}(s) $$

and by observing that for $s \in [t_{k-1}, t_k]$ and $j = 1:d$, we have $\psi_j(s) = \overline{a_j^{[k,n]}}$ and $\tilde{S}^j(\psi_{s+}, t-s) = \tilde{L}^{j,(n+1-k)}(a_j^{[k,n]}, (t-s)_{[k,n]}).$

2. Moreover, emphasizing the dependence upon $\lambda^0$ via the notation $\tilde{S}^0 = \tilde{S}_{\lambda^0}$, we obtain a similar recurrence equation: for $0 < t_1 < \dots < t_n < t$,

$$ \tilde{S}_{\lambda^0}(\psi, t) = \prod_{k=1}^{n} \tilde{S}_{\lambda^0 \chi_{[t_{k-1}, t_k]}}(\psi_{t_{k-1}^{+}}, t - t_{k-1}). $$

3. Functional Laplace transforms may be of some importance in applications. For example, suppose that we are interested in the assessment of the medical or financial costs of an epidemic in $d$ regions, while the counting process is only observed over deterministic region-specific time intervals $[C_j^{2k}, C_j^{2k+1}]$, $j = 1:d$, $k \ge 0$. In this case, the relevant Laplace transform would concern the functions $\psi_j(s) = a_j \theta_j(s) \sum_{k \ge 0} \chi_{[C_j^{2k}, C_j^{2k+1}]}(s)$, where $\theta_j(s)$ corresponds to the cost induced by the occurrence of an event in region $j$ at time $s$ and the vector $a = (a_1, \dots, a_d)$ holds for variables of the Laplace transform $L_C^j(a, t)$ of censored observations, that is $\tilde{S}^{j'}(\psi, t) = L_C^{j'}(a, t)$.

# 3 Infinite Divisibility and Distribution Probabilities of Counts

## 3.1 Preliminaries

We detail here the precise form of the infinitely divisible property for Hawkes processes and the integral equations related to the probability distribution of its counts. First we start with the following analytic result.

**Proposition 3.1.** Let $f(z) = \sum_{l \in \mathbb{N}^d} \alpha_l z^l$ be a $d$-multivariate analytic function on the complex polydisc $D^d$ with $z = (z_1, \dots, z_d)$, $l = (l_1, \dots, l_d)$ and $z^l = z_1^{l_1} \times \dots \times z_d^{l_d}$. Then $g(z) = e^{f(z)}$ admits the representation $g(z) = e^{\alpha_0 (1 + \sum_{k \neq 0} \beta_k z^k)}$ where for $k = (k_1, \dots, k_d) \neq 0$ and $|k| = k_1 + \dots + k_d$

$$ \beta_k = \sum_{r=1}^{|k|} \frac{1}{r!} \sum_{\substack{l^{(1)} \neq 0, \dots, l^{(r)} \neq 0 \\ l^{(1)} + \dots + l^{(r)} = k}} \alpha_{l^{(1)}} \times \dots \times \alpha_{l^{(r)}}. \qquad (9) $$
Conversely, we have for $l \neq 0$

$$ \alpha_l = \sum_{r=1}^{|l|} \frac{(-1)^{r-1}}{r} \sum_{\substack{k^{(1)} \neq 0, \dots, k^{(r)} \neq 0 \\ k^{(1)} + \dots + k^{(r)} = l}} \beta_{k^{(1)}} \times \dots \times \beta_{k^{(r)}}. \quad (10) $$

*Proof.* Simply apply the series expansion of exp and log(1 + ·) functions about point 0. □

**Corollary 3.1.1.** Let $X$ be an $\mathbb{N}^d$-valued random variable with probability generating function $g(z) = \mathbf{E}(z^X) = \sum_{k \in \mathbb{N}^d} p_k z^k$ and $p_k = P(X = k)$. If $p_0 > 0$, then there exists a unique multivariate analytic function $f(z) = \sum_{l \in \mathbb{N}^d} \alpha_l z^l$ such that $g(z) = e^{f(z)}$, with $\alpha_0 = \log(p_0)$ and for all $l \neq 0$, $\alpha_l$ satisfies Equation (10) with $\beta_k = p_k/p_0$.

**Remark.** Notice that the summation in the previous formulas describes the number of distinct paths in the lattice $\mathbb{N}^d$ starting from point $(0, \dots, 0)$ and leading to point $k = (k_1, \dots, k_d)$ in $r$ steps. Moreover the calculation of $\beta_k$ (or $\alpha_l$) is recursive/causal since it only requires the knowledge of coefficients $\alpha_{l(p)}$ with indices $l(p) \le k$, that is to say $1 \le l_j^{(p)} \le k_j$, $j = 1, \dots, d$.

**Proposition 3.2.** Let $X$ be a multivariate random variable with probability distribution function $F(dx)$ on $\mathbb{R}_+^d$, Laplace transform $L(a) = \mathbf{E}(e^{-\top a X}), a \in \mathbb{R}_+^d$, and such that $0 < p_0 = F(\{0\}) < 1$. Let $\bar{F}^{(*)r}(dx)$ denote the $r^{th}$ convolution of $\bar{F} = F(dx)/(1-p_0)$ probability measure on $(\mathbb{R}_+^d)_*$. Then there exists a unique finite signed measure $\gamma(dx)$ on $(\mathbb{R}_+^d)_*$ such that $\gamma((\mathbb{R}_+^d)_*) = -\log(p_0)$ and

$$ L(a) = \exp \left( \log(p_0) + \int_{(\mathbb{R}_+^d)_*} e^{-\top a x} \gamma(dx) \right) \quad (11) $$

with

$$ \gamma(dx) = \sum_{r=1}^{\infty} \frac{(-1)^{r-1}}{r} \left( \frac{1-p_0}{p_0} \right)^r \bar{F}^{*r}(dx). \quad (12) $$

*Proof.* Since $\bar{F}(dx) = F(dx)/(1-p_0)$ is a probability distribution on $(\mathbb{R}_+^d)_*$ with Laplace transform $\bar{L}(a)$, the convolution rule applies and yields:

$$ \bar{L}^r(a) = \left( \int_{(\mathbb{R}_+^d)_*} e^{-\top a x} \bar{F}(dx) \right)^r = \int_{(\mathbb{R}_+^d)_*} e^{-\top a x} \bar{F}^{(*)r}(dx), $$

such that

$$ \log(L(a)) = \log \left( p_0 + \int_{(\mathbb{R}_+^d)_*} e^{-\top a x} F(dx) \right) = \log (1 + (1 - p_0) (\bar{L}(a) - 1)). \quad (13) $$

Since $|(1 - p_0)(\bar{L}(a) - 1)| < 1$, the series expansion yields

$$ \begin{aligned} \log(L(a)) &= \sum_{r=1}^{\infty} \frac{(-1)^{r-1}}{r} (1-p_0)^r \left( C_0^r (-1)^r + \sum_{j=1}^{r} C_j^r (-1)^{r-j} \bar{L}^j(a) \right) \\ &= -\sum_{r=1}^{\infty} \frac{1}{r} (1-p_0)^r + \sum_{j=1}^{\infty} \left( \sum_{r=j}^{\infty} \frac{1}{r} (1-p_0)^r C_j^r \right) (-1)^{j-1} \bar{L}^j(a). \end{aligned} $$
Taking into account the equality $(1-x)^{-(j+1)} = \sum_{r=j}^{\infty} C_j^r x^{r-j}$, we get

$$ 
\begin{aligned}
\log(L(a)) &= \log(p_0) + \sum_{j=1}^{\infty} \frac{(-1)^{j-1}}{j} \left(\frac{1-p_0}{p_0}\right)^j \bar{L}^j(a) \\
&= \log(p_0) + \int_{(\mathbb{R}_+^d)_*} e^{-\tau ax} \left( \sum_{j=1}^{\infty} \frac{(-1)^{j-1}}{j} \left(\frac{1-p_0}{p_0}\right)^j \bar{F}^{(*)j}(dx) \right).
\end{aligned}
$$

Now, we observe that

$$ \gamma(dx) = \sum_{j=1}^{\infty} \frac{(-1)^{j-1}}{j} \left(\frac{1-p_0}{p_0}\right)^j \bar{F}^{(*)j}(dx) $$

can be written as the difference of two non-negative measures:

$$ 
\begin{aligned}
\gamma_1(dx) &= \sum_{j=1}^{\infty} \frac{1}{2j-1} \left(\frac{1-p_0}{p_0}\right)^{2j-1} \bar{F}^{(*)}(2j-1)(dx) \\
\gamma_2(dx) &= \sum_{j=1}^{\infty} \frac{1}{2j} \left(\frac{1-p_0}{p_0}\right)^{2j} \bar{F}^{(*)2j}(dx).
\end{aligned}
$$

Since $\bar{F}^{(*)j}$ is for each $j \ge 1$ a probability measure on $(\mathbb{R}_+^d)_*$, we can deduce that $\gamma_1$ and $\gamma_2$ are non negative possibly infinite measures.

It remains to show that $\gamma_1$ and $\gamma_2$ are finite. In this aim, it is sufficient to consider any Borelian interval product $B = \prod_{j=1}^d [0, b_j]$ with $b_j \ge 0$ and to observe that $\bar{F}^{(*)r}(B) \le \bar{F}^r(B)$. Therefore, $\gamma_1(B) \le \sum_{j=1}^{\infty} \frac{1}{2j-1} \left(\frac{1-p_0}{p_0}\right)^{2j-1}$ and so $\gamma$ is a $\sigma$-additive signed measure on $(\mathbb{R}_+^d)_*$. Moreover $\gamma((\mathbb{R}_+^d)_*) = \sum_{r=1}^{\infty} \frac{(-1)^{r+1}}{r} \left(\frac{1-p_0}{p_0}\right)^r = \log(1 + \frac{1-p_0}{p_0}) = -\log(p_0)$. $\square$

## 3.2 Infinite Divisibility Property

First, we recall the most achieved formulation of the Lévy Kinchine representation theorem for infinitely divisible random vector $X$ with non-negative components, i.e. $P(X \in \mathbb{R}_+^d) = 1$, see [28].

**Theorem 3.3.** Let $X$ be a multivariate random variable with values on $\mathbb{R}_+^d$ and Laplace transform $L(a) = \mathbf{E}(e^{-\tau a X}), a \in \mathbb{R}_+^d$. Then $X$ is infinitely divisible if and only if there exist a unique $\beta \in \mathbb{R}_+^d$ and a $\sigma$-finite measure $\gamma$ on $\mathbb{R}_+^d$ such that $\gamma(\{0\}) = 0$, $\int_{\mathbb{R}_+^d} (\|x\| \wedge 1)\gamma(dx) < \infty$ and

$$ L(a) = \exp \left( - \tau a \beta + \int_{(\mathbb{R}_+^d)_*} (e^{-\tau a x} - 1) \gamma(dx) \right). \quad (14) $$

Nevertheless, one can say a little more for non negative integer valued random vectors.

**Theorem 3.4.** Let $X$ be a multivariate random variable with values in $\mathbb{N}^d$ with $P(X = 0) = p_0 > 0$ and Laplace transform $L(a) = \mathbf{E}(e^{-\tau a X}), a \in \mathbb{R}_+^d$. Then $X$ is infinitely divisible if and only if there exists a (unique) finite measure $\gamma(dx) = \sum_{l \in (\mathbb{N}^d)_*} \alpha_l \delta_l(dx)$ on $(\mathbb{N}^d)_*$, satisfying $\sum_{l \in (\mathbb{N}^d)_*} \alpha_l = -\log(p_0)$ and such that

$$ L(a) = p_0 \exp \left( \sum_{l \neq 0} e^{-\tau a l} \alpha_l \right), $$

where the coefficients $\alpha_l, l \neq 0$ satisfy Equation (10) with $\beta_k = p_k/p_0$.
Proof. Since $L(a)$ can be written as $g(z) = p_0 + \sum_{k \neq 0} z^k p_k$ with $z = (e^{-a_1}, \dots, e^{-a_d})$ and $p_0 > 0$, then according to Corollary 3.1.1, we have a unique function $f(z) = \sum_{l \in \mathbb{N}^d} \alpha_l z^l$ satisfying $g(z) = e^{f(z)}$ with $\alpha_0 = \log(p_0)$ and the $\alpha_l$ satisfying Equation (10) with $\beta_k = p_k/p_0$.

We now prove that in the representation (14), $\beta = 0$ and $\gamma$ is finite.

Indeed, if we take $a = (u, \dots, u) \in \mathbb{R}_+^d$, then for all $k \in (\mathbb{N}^d)_*$, we have $0 \le e^{-\top ak} = e^{-u|k|} \le e^{-u}$ and then $p_0 \le L(a) \le p_0 + e^{-u}(1-p_0) \underset{u \to +\infty}{\longrightarrow} p_0$.

On the other hand, for $x \neq 0$, we have $0 \le 1 - e^{-u|x|} \underset{u \to +\infty}{\uparrow} 1$; hence, by Beppo Levi's lemma we have, (14) yields $L(a) \underset{u \to +\infty}{\longrightarrow} \exp(-\infty|\beta| - \gamma((\mathbb{R}_+^d)_*)) = p_0 > 0$. But this holds true only for $\beta = 0$ and $-\log(p_0) = \gamma((\mathbb{R}_+^d)_*) < \infty$. This also means that $\gamma$ is a finite measure.

Moreover, since $f$ can be written $f((e^{-a_1}, \dots, e^{-a_d})) = \log(p_0) + \int_{(\mathbb{N}^d)_*} e^{-\top ax} \alpha(dx)$, with $\alpha(dx) = \sum_{l \neq 0} \alpha_l \delta_l(dx)$ and that the measure $\gamma(dx)$ is unique, we simply have $\gamma(dx) = \alpha(dx)$, whose support is $(\mathbb{N}^d)_*$. $\square$

**Remark.** If $X$ is indefinitely divisible and $\mathbb{N}^d$ valued, we always have $E(x) = \sum_{k \neq 0} k \alpha_k$ and $Cov(X) = \sum_k k^\top k \alpha_k$ for both expressions, either finite or infinite.

### 3.3 Distribution Probabilities of Hawkes Process Countings

Without loss of generality, let us assume for the indefinitely divisible Hawkes processes $\tilde{N}^{j'}(t)$, $j' = 0: d$, that $\lambda_0(s)$ and $\phi$ are integrable over any interval $[0, t]$, and satisfy $\tilde{p}_0^{j'}(t) > 0$ and $\tilde{p}_0^{j'}(t) > 0$ for all $t > 0$.

**Proposition 3.5.** Let $\tilde{N}^{j'}$, $j' = 0: d$, be multivariate Hawkes processes driven by integrable functions $(\phi_j^{j'}, \phi)$ and let $\tilde{p}_k^{j'}(t) = P_{(\phi_j^{j'}, \phi)}(\tilde{N}^{j'}(t) = k), k \in \mathbb{N}^d$, the probability distribution functions of $\tilde{N}^{j'}(t)$, $j' = 0: d, t \in \mathbb{R}_+$. Then we have

$$ \frac{\tilde{p}_l^{j'}(t)}{\tilde{p}_0^{j'}} - \sum_{r=2}^{|l|} \frac{(-1)^r}{r} \sum_{\substack{k^{(1)} \neq 0, \dots, k^{(r)} \neq 0 \\ k^{(1)} + \dots + k^{(r)} = l}} \frac{\tilde{p}_{k^{(1)}}^{j'}(t)}{\tilde{p}_0^{j'}} \times \dots \times \frac{\tilde{p}_{k^{(r)}}^{j'}(t)}{\tilde{p}_0^{j'}} = \sum_{j=1}^d \int_0^t (\tilde{p}_{l-e_j}^j(t-s)) \phi_j^{j'}(s) ds. \quad (15) $$

Proof. First, let us deal with $j' = 0$ and then observe that equations of Theorem 2.2 reduce to the following forms for a single time $t > 0$ and $j' \in \{1, \dots, d\}$:

$$ \tilde{L}_{\phi_j^{j'}}^{j'(1)}(a, t) = \tilde{p}_0^{j'}(t) \exp \left( \sum_{j=1}^{d} \int_{0}^{t} e^{-a_j} \tilde{L}^{j,(1)}(a, (t-s)) \phi_j^{j'}(s) ds \right). \quad (16) $$

Since $\tilde{N}^j$, $j = 0: d$, are infinitely divisible, then according to Theorem 3.4, there exist finite
measures $\gamma^{j'}(t, dx)) = \sum_{k \in (\mathbb{N}^d)_*} \alpha_k^{j'}(t) \delta_k(dx)$, for all $t > 0$ and $j' = 0 : d$, such that

$$
\begin{align*}
\sum_{l \neq 0} e^{-\top al} \alpha_l^{j'}(t) &= \sum_{j=1}^{d} \int_{0}^{t} e^{-a_j} \tilde{L}^{j,(1)}(a, (t-s)) \phi_j^{j'}(s) ds \\
&= \sum_{j=1}^{d} \int_{0}^{t} e^{-a_j} \left( \sum_{l \in \mathbb{N}^d} e^{-\top al} \tilde{p}_l^j(t-s) \right) \phi_j^{j'}(s) ds \\
&= \sum_{l \in \mathbb{N}^d} e^{-\top al} \sum_{j=1}^{d} \int_{0}^{t} e^{-a_j} \left( \tilde{p}_l^j(t-s) \right) \phi_j^{j'}(s) ds \\
&= \sum_{l \neq 0} e^{-\top al} \left( \sum_{j=1}^{d} \int_{0}^{t} \left( \tilde{p}_{l-e_j}^j(t-s) \right) \phi_j^{j'}(s) ds \right)
\end{align*}
$$

with $e_j$ being the $j^{th}$ unit vector of $\mathbb{R}^d$. Therefore we have $\alpha_l^{j'}(t) = \sum_{j=1}^d \int_0^t \tilde{p}_{l-e_j}^j(t-s) \phi_j^{j'}(s) ds$. Then, according to Equation (10), we obtain the following recursive formula for $l \neq 0$ and $j' = 0 : d$:

$$
\frac{\tilde{p}_l^{j'}}{\tilde{p}_0^{j'}}(t) - \sum_{r=2}^{|l|} \frac{(-1)^r}{r} \sum_{\substack{k^{(1)} \neq 0, \dots, k^{(r)} \neq 0 \\ k^{(1)} + \dots + k^{(r)} = l}} \frac{\tilde{p}_{k^{(1)}}^{j'}}{\tilde{p}_0^{j'}}(t) \times \dots \times \frac{\tilde{p}_{k^{(r)}}^{j'}}{\tilde{p}_0^{j'}}(t) = \sum_{j=1}^d \int_0^t (\tilde{p}_{l-e_j}^j(t-u)) \phi_j^{j'}(u) du.
$$

□

**Example.** Let us consider the case $d = 1$ with $\phi^0 = \lambda^0$ and $\phi = \phi^1$. From the population representation (14), we first have for $l$ equals 0 and 1: $\tilde{p}_0^{j'}(t) = e^{-\Phi^{j'}(t)}$ and $\tilde{p}_1^{j'}(t) = e^{-\Phi^{j'}(t)}\Phi^{j'}(t)$. Then for $l = 2$, we get:

$$
\tilde{p}_2^{j'}(t) e^{\Phi^{j'}(t)} = \frac{1}{2} (\Phi^{j'}(t))^2 + \int_0^t e^{-\Phi^1(t-u)} \Phi^1(t-u) \phi^{j'}(u) du,
$$

and so on...

**Remark.** Clearly, Proposition 3.5 can be generalized to multi-dimensional multi-temporal counts $\top \tilde{N}^{j'} = (\top \tilde{N}^{j'}(t_1), \dots, \top \tilde{N}^{j'}(t_n))$ of dimension $d \times n$. Indeed, the probability distribution of the process $\tilde{N}^{j'}$ is infinitely divisible and one can adapt Equation (15) to $(\mathbb{N}^d)^n$ (instead of $\mathbb{N}^d$), providing summation formulas over increasing paths passing through the different points $\tilde{N}^{j'}(t_r); r = 1 : n$.

# 4 Integral Equations of the First Two Multi-temporal Moments

The Laplace transform is an important tool that characterizes all properties of probability distributions of non negative random variables. This is also the case of the multivariate multi-time distributions of a Hawkes process with Laplace transform $\tilde{L}^{j',(n)}(a, t)$ satisfying Equation (7). Hence, if $a = (a_j^m : m = 1 : n, j = 1 : d)$ and $t = (t_1, \dots, t_n)$, for any non negative integer valued set of triplets $Q = \{(q_1, m_1, j_1), \dots, (q_r, m_r, j_r)\}$ such that $q_k \ge 1, 1 \le m_k \le n, 1 \le j_k \le d), k = 1 : r$ and $|q| = \sum_{k=1}^r q_k$, we have the following formula for the $Q^{th}$ moment of
$\tilde{N}^{j'}, j' = 0 : d:$

$$
\begin{align*}
M^{j', (Q)}(t) &= \mathbf{E}_{(\phi_{\bullet}^{j'}, \phi)} \left( \prod_{k=1}^{r} [\tilde{N}_{jk}^{j'}(t_{mk})]^{q_k} \right) \\
&= (-1)^{|q|} \partial_Q \tilde{L}^{j', (n)}(a, t) \Big|_{a=0},
\end{align*}
$$

with

$$
\partial_Q = \partial_{(a_{j_1}^{m_1}, q_1 \dots (a_{j_r}^{m_r}, q_r)}^{|q|}.
$$

Here, $\mathcal{Q}$ represents the set $\mathcal{Q}$ of size $|q|$ containing $q_1$ copies of $a_{j_1}^{m_1}$, ..., $q_r$ copies of $a_{j_r}^{m_r}$; hence, $\partial_Q$ can be replaced by $\partial_Q$ without any risk of confusion.

Therefore using the Faà di Bruno formula for partial derivatives of composed functions, namely $\tilde{L}^{j', (n)}(a, t) = \exp(\theta^{j'}(a, t))$ satisfying $\tilde{L}^{(j', n)}(0_{[d]}, t) \equiv 1$ for all $t, j', n$, we get

$$
M^{j', (Q)}(t) = (-1)^{|q|} \left. \sum_{\pi \in \mathcal{P}} \prod_{\tilde{Q}_l \in \pi} \partial_{\tilde{Q}_l} \theta^{j'}(a, t) \right|_{a=0}, \quad (17)
$$

with $\pi$ running through the set $\mathcal{P}$ of all partitions $\{\tilde{Q}_1, \dots, \tilde{Q}_p\}$ of the set $\mathcal{Q}$. For example in the case of $\partial_{(a_1^1)^2 a_2^2}^3$ we have $\mathcal{Q} = \{a_1^1, a_1^1, a_2^2\}$. Moreover, as the components $\tilde{Q}_l$ of the partition $\pi$ may contain several copies of a same element $a_j^m$, Equation (17) simplifies a lot in the case where some of the $q_m$ are strictly greater than 1, such that the number of distinct terms then reduces to $\prod_{k=1}^r q_k! (k!)^{q_k}$.

In what follows, we focus on the derivation of formal and/or implicit equations related to the
first and second order moments.

4.1 Equations of First Order Moments

**Proposition 4.1.** Consider the set of vector means $\tilde{M}^{j', (1)}(t) = E_{(\phi_{\bullet}^{j'}, \phi)} (\tilde{N}^{j'}(t))$, $j' = 0 : d$, and the $d \times d$ matrix of base vector means $\tilde{M}^{(1)}(t) = [\tilde{M}^{1,(1)}, \dots, \tilde{M}^{d,(1)}]$. Let $\Phi(t) = \int_0^t \phi(s)ds$ and $\Lambda^0(t) = \int_0^t \lambda^0(s)ds$. Then, $\tilde{M}^{(1)}$ is solution of the following matrix integral equation

$$
\tilde{M}^{(1)}(t) = \Phi(t) + \tilde{M}^{(1)} * \phi(t), \quad (18)
$$

that corresponds to

$$
\tilde{M}_{l}^{j', (1)}(t) = \Phi_{l}^{j'}(t) + \sum_{j=1}^{d} \int_{0}^{t} \tilde{M}_{l}^{j, (1)}(t-s) \phi_{j}^{j'}(s) ds, \quad (19)
$$

for $1 \le j', l \le d$. Moreover,

$$
\tilde{M}^{0,(1)}(t) = \Lambda^0(t) + \tilde{M}^{(1)} * \lambda^0(t). \quad (20)
$$

*Proof.* Let $k = 1$, $t \in \mathbb{R}_{+*}$ and $\top a = (a_1, \dots, a_d) \in (\mathbb{R}_+)^d$, then Equation (8) yields for $j' = 0 : d$

$$
\tilde{L}^{j', (1)}(a, t) = \tilde{p}_0^{j'}(t) \exp \left( \sum_{j=1}^{d} e^{-a_j} \tilde{L}^{j, (1)}(a, \cdot) * \phi_j^{j'}(t) \right). \quad (21)
$$

Then, applying Equation (17), we get

$$
\tilde{M}^{j', (1)}(t) = -\partial_{a_l}^1 \tilde{L}^{j', (1)}(0_{[d]}, t).
$$
Hence, for each $j' = 1 : d$,

$$ \partial_{a_l}^1 \tilde{L}^{j', (1)} (a, t) = \tilde{L}^{j', (1)} (a, t) \left( \sum_{j=1}^d \left( (-\mathbf{1}_{l=j}) e^{-a_j} \tilde{L}^{j,(1)}(a,.) + e^{-a_j} \partial_{a_l}^1 \tilde{L}^{j,(1)}(a,.) \right) * \phi_j^{j'} \right) (t), \quad (22) $$

and, consequently,

$$ \begin{aligned} \partial_{a_l}^1 \tilde{L}^{j', (1)} (0_{[d]}, t) &= \tilde{L}^{j', (1)} (0_{[d]}, t) \sum_{j=1}^d \left( (-\mathbf{1}_{l=j}) \tilde{L}^{j,(1)}(0_{[d]},.) + \mathbf{1} \times \partial_{a_l}^1 \tilde{L}^{(1)}(0_{[d]},.) \right) * \phi_j^{j'}(t) \\ &= \sum_{j=1}^d \left( (-\mathbf{1}_{l=j}) - \tilde{M}_l^{j,(1)} \right) * \phi_j^{j'}(t). \end{aligned} $$

since $\tilde{L}^{j', (1)}(0_{[d]}, t) \equiv 1$, for $j' = 0 : d$. Therefore, for all $j' = 1 : d$ and $l = 1 : d$, we obtain

$$ \tilde{M}_l^{j', (1)}(t) = \Phi_l^{j'}(t) + \sum_{j=1}^{d} \left( \tilde{M}_l^{j,(1)} * \phi_j^{j'} \right) (t). \quad (23) $$

The latter formula can be rewritten under a matrix form and replacing $\phi_j^{j'}$ with $\lambda^0$ yields the second part of the proposition. $\square$

### 4.1.1 Explicit Formulas for First Order Moments

Actually, the previous implicit/integral equations for the first order moments can be developed explicitly.

Remind that the $r^{th}$-convolution of non-negative $d \times d$-matrix function $\phi(t) = (\phi_i^j)(t)$ on $\mathbb{R}_+$ is well defined recursively as $\phi^{(*)r} = \phi^{(*)(r-1)} * \phi$, e.g. $(\phi^{(*)2})_{j'}^j(t) = \sum_{r=1}^d \phi_j^r * \phi_{r'}^{j'}(t) = \sum_{r=1}^d \phi_r^{j'} * \phi_j^r(t)$. We also convene that for $r = 0$, we have

$$ \phi^{(*)0}(t) = \mathrm{diag}(\delta_0(\mathrm{d}t), \dots, \delta_0(\mathrm{d}t)) = I_d \delta_0(\mathrm{d}t), $$

where $\delta_0(\mathrm{d}t)$ is the Dirac distribution (or generalized function) supported by $0 \in \mathbb{R}_+$.

**Proposition 4.2.** Finite or infinite, the mean function $\tilde{M}^{(1)}(t)$ is written

$$ \tilde{M}^{(1)}(t) = \Phi * \left( \sum_{r=0}^{\infty} \phi^{(*)r} \right) (t) = \left( \sum_{r=0}^{\infty} \phi^{(*)r} \right) * \Phi(t) \quad (24) $$

with derivative

$$ \tilde{m}^{(1)}(t) = \sum_{r=1}^{\infty} \phi^{(*)r}(t). \quad (25) $$

Similarly, we have

$$ \tilde{M}^{0,(1)}(t) = \left( \sum_{r=0}^{\infty} \phi^{(*)r} \right) * \Lambda^0(t), \quad \text{with derivative} \quad \tilde{m}^{0,(1)}(t) = \left( \sum_{r=0}^{\infty} \phi^{(*)r} \right) * \lambda^0(t). \quad (26) $$

Moreover, if $\|\phi\|_{L^1} = \sup(\|\phi_i^j\|_{L^1} : 1 \le i, j \le d) \le \alpha$ with $d\alpha < 1$, then $\tilde{M}^{(1)}(t)$ is finite for all $t \ge 0$. Additionally, if $\lambda^0 \in L^1$, then $\tilde{m}^{0,(1)} \in L^1$ with primitive $\tilde{M}^{0,(1)}(t)$ finite for all $t \ge 0$.
Proof. Since all functions $\phi_j^{j'}$ are non negative, one can expand by iteration Equation (18) and get the following increasing converging series of matrix functions (whose components may be finite or infinite):

$$ \tilde{M}^{(1)}(t) = \Phi * \left( \sum_{r=0}^{\infty} \phi^{(*)r} \right) (t). $$

Moreover, observe that using differential rules for convolution yields $\Phi_l^r * \phi_j^l = \phi_l^r * \Phi_j^l$, which implies that commutativity of matrix convolution $\Phi * \phi = \phi * \Phi$ holds in this specific case. Iterating the process proves the commuting property $\Phi * \phi^{(*)r} = \phi^{(*)r} * \Phi$ for all $r \ge 0$.

Next, if $\tilde{M}^{(1)}$ is finite on an interval $I$ of $\mathbb{R}_+$, then it is differentiable with derivative

$$ \tilde{m}^{(1)}(t) = \sum_{r=1}^{\infty} \phi^{(*)r}(t). $$

and, due to Equation (20), if $\tilde{M}^{0,(1)}$ is finite, then it is differentiable with derivative:

$$ \tilde{m}^{0,(1)}(t) = \lambda^0(t) + \tilde{m}^{(1)} * \lambda^0(t) = \left( \sum_{r=0}^{\infty} \phi^{(*)r} \right) * \lambda^0(t). $$

We now prove that the given conditions are sufficient to ensure finite values of the means and their derivatives. For that purpose, let us recall the Young's inequalities for convolution in Lebesgue spaces $L^p$: For any $f \in L^p$, $g \in L^q$ and $1/p + 1/q = 1 + 1/r$, we have $f * g \in L^r$ and $\|f * g\|_{L^r} \le \|f\|_{L^p} \|f\|_{L^q}$.

Thus, in the particular case where $r = p = q = 1$ and $\|\phi_i^j\|_{L^1} \le \alpha, 1 \le i, j \le d$, we get $\|\sum_{l=1}^d \phi_i^l * \phi_l^j\|_{L^1} \le d\alpha^2 = (d\alpha)^2/d$.

Iterating the procedure, we get $\|\phi^{(*)r}\|_{L^1} \le (d\alpha)^r/d$ for $r \ge 1$. For $d\alpha < 1$, we finally get $\|\sum_{r=1}^\infty \phi^{(*)r}\|_{L^1} \le \sum_{r=1}^\infty (d\alpha)^r/d = \alpha/(1 - d\alpha) < \infty$.

Consequently the series $\sum_{r=1}^\infty \phi^{(*)r}$ converges in $L^1$ to $\tilde{m}^{(1)}$ and, hence, all its components are almost everywhere finite (and integrable).

Clearly, $\Phi$ the primitive of $\phi$ is therefore bounded and differentiable on any interval, the mean $\tilde{M}^{(1)}$ is also differentiable and finite for all $t \ge 0$. Similar results hold true for $M^{0,(1)}$ if $\lambda^0$ is in $L^1$. $\square$

### 4.1.2 Examples and Remarks

1. We start with the classical example of constant $\phi$, that is $\phi(t) = W\chi(t)$ where $W = (w_i^j)$ is a constant matrix and $\chi(t) = \mathbf{1}_{[0,\infty]}(t)$. Hence $\phi \in L^\infty$ and condition of Proposition 4.2 holds true. Then, we observe that $\chi^{(*)r}(t) = t^{r-1}/(r-1)!$ and, therefore, we get $\tilde{m}^{(1)}(t) = \sum_{r=1}^\infty \phi^{(*)r}(t) = W(\sum_{k=0}^\infty W^k t^k/k!) = We^{tW}$.

Consequently, since $\Phi(t) = tW$, we also have $\tilde{M}^{(1)}(t) = \Phi * (\sum_{r=0}^\infty \phi^{(*)r})(t) = tW + W^2(\sum_{r=0}^\infty W^r t^{r+2}/(r+2)!) = e^{tW} - I_d$.

It is worth noticing that using the Cauchy formula related to $r$-repeated integrals namely $I^{(r)}[g](t) = \int_0^t \int_0^{u_r} \dots \int_0^{u_1} g(u_1)du_1 \dots du_r = \int_0^t g(u)(t-u)^{(r-1)/(r-1)!}du$, we get for an initial igniting intensity $\lambda^0$, the formulas $\tilde{m}^{0,(1)}(t) = \sum_{r=0}^\infty W^r I^{(r)}[\lambda^0](t)$ and $\tilde{M}^{0,(1)}(t) = \sum_{r=0}^\infty W^r I^{(r+1)}[\lambda^0](t)$.

The latter expression can be useful in the case of the existence of recursive formulas for the calculus of repeated integrals of function such as exponential, sine, etc.

2. The condition of Proposition 4.2 holds for $\phi_i^j(t) = \alpha_i^j e^{-\beta_i^j t}$, if $d\alpha_i^j < \beta_i^j$, $1 \le i, j \le d$ and $\lambda^0(t)$ is a constant vector.
3. One can also consider the interesting class of $\Gamma$ intensities admitting many different shapes. For example if $\phi_i^j(t) = w_i^j f(t; \kappa, \theta) = w_i^j \frac{t^{\kappa-1}e^{-t/\theta}}{\Gamma(\kappa)\theta^{\kappa}}$, then using the stability of the $\Gamma$ probability distribution family for convolution, we get for $W = (w_i^j)_{i,j}$ the expression

$$ \phi^{(*)r}(t) = W^r f(t; r\kappa, \theta). $$

Therefore, we obtain $\tilde{m}^{(1)}(t) = \sum_{r=1}^{\infty} \phi^{(*)r}(t) = \sum_{r=1}^{\infty} \frac{W^r}{\theta\Gamma(r\kappa)} \left(\frac{t}{\theta}\right)^{r\kappa-1} e^{-t/\theta}$ and $\tilde{M}^{(1)}(t) = \sum_{r=1}^{\infty} \frac{W^r}{\Gamma(r\kappa)} \gamma(r\kappa, t/\theta)$, with $\gamma$ being the incomplete gamma function.

* In particular for $\kappa = 1$ we have

$$ \tilde{m}^{(1)}(t) = \frac{e^{-t/\theta}}{t} \left( e^{W \frac{t}{\theta}} - I_d \right) \text{ and } \tilde{M}^{(1)}(t) = \sum_{r=1}^{\infty} \frac{W^r}{r!} \gamma(r, t/\theta) $$

* If $\kappa = 2$ and $W$ has a square root say $W^{1/2}$ (e.g if $W$ is positive definite), then $\sum_{r=1}^{\infty} \phi^{(*)r}(t) = \frac{e^{-t/\theta}}{t} (\cosh(W^{1/2} \frac{t}{\theta}) - I_d)$.

In any case, one has to use the Jordan form of $W$ and $W^{1/2}$ to make easier the computation (or approximation) of the series.

4. One can arguably consider that $\sum_{r=0}^{\infty} \phi^{(*)r}(t)$ is simply the inverse (in the both senses of matrix product and convolution product) of the matrix valued generalized function $I_d\delta_0(ds) - \phi(s)ds$.

## 4.2 Second Order Moments

Taking $k = 2, a = (a_\bullet^1, a_\bullet^2)$ and $t = (t_1, t_2)$ in Equation (8) yields for $j' = 0 : d$

$$ \tilde{L}^{j', (2)}(a, t) = \tilde{p}_0^{j'}(t_2) \exp \left( \sum_{j=1}^{d} \int_{0}^{t_1} e^{-(a_j^1+a_j^2)} \tilde{L}^{j', (2)}(a, (t_1-u, t_2-u)) \phi_j^{j'}(u) du \right) \quad (27) $$

$$ \times \exp \left( \sum_{j=1}^{d} \int_{t_1}^{t_2} e^{-a_j^2} \tilde{L}^{j', (1)}(a_\bullet^2, (t_2-u)) \phi_j^{j'}(u) du \right). \quad (28) $$

Thus, according to moment formula (17) with $Q = ((1, m_1, k), (1, m_2, l))$ and for $j' = 0 : d$, we obtain:

$$ \tilde{M}_{k,l}^{j', (2)}(t_{m_1}, t_{m_2}) = \mathbf{E}_{(\phi_\bullet^j, \phi)} \left( \tilde{N}_k^{j'}(t_{m_1}) \tilde{N}_l^{j'}(t_{m_2}) \right) = \partial_{a_k^{m_1} a_l^{m_2}}^2 \tilde{L}^{j', (2)}(a, t) \Big|_{a=0} $$

We can also write expressions of the covariance function satisfying:

$$ \tilde{C}_{k,l}^{j'}(t_{m_1}, t_{m_2}) = \tilde{M}_{k,l}^{j', (2)}(t_{m_1}, t_{m_2}) - \tilde{M}_k^{j', (1)}(t_{m_1}) \tilde{M}_l^{j', (1)}(t_{m_2}). $$

Below, we consider two cases: the case corresponding to a single time where we are interested in the covariance function of the $d$-vector $\tilde{N}^j(t_1)$, and the case of distinct times $t_1 < t_2$ where we consider the covariance function of the $2d$-vector $\top(\top \tilde{N}^j(t_1), \top \tilde{N}^j(t_2))$ with covariance matrix partitioned as follows

$$ \tilde{C}^{j'}(t_1, t_2) = \begin{pmatrix} \tilde{C}^{j'}(t_1, t_1) & \tilde{C}^{j'}(t_1, t_2) \\ \tilde{C}^{j'}(t_2, t_1) & \tilde{C}^{j'}(t_2, t_2) \end{pmatrix}. $$

In this partitioned matrix, the diagonal blocks correspond to single times $t_1$ and $t_2$.
## Remarks

* The case $m_1 = m_2$ means that we must only deal with a single time $t_1$ (resp. $t_2$) and hence to consider $a_\bullet^1$ (resp. $a_\bullet^1 + a_\bullet^2$) instead of $a = (a_\bullet^1, a_\bullet^2)$, and then we must simply deal with $\tilde{L}^{j',(1)}$ instead of $\tilde{L}^{j',(2)}$. This finally leads to twice differentiate $\tilde{L}^{j',(1)}(a, s)$ with $a = a_\bullet^1$ and $s = t_1$ (resp. $a = a_\bullet^1 + a_\bullet^2 u$ and $s = t_2$).

* The case $m_1 \neq m_2$, say $m_1 = 1, m_2 = 2$, implies that $t_1 \neq t_2$ and in this case $\partial_{a_{j_1}^1 a_{j_2}^2}^2 \tilde{L}^{j',(2)}(a, t)$ simplifies a lot since the third term in Equation (27) only depends on $a_\bullet^2$.

* Finally, note that $\tilde{C}^{j'}(t_2, t_1) = \tilde{C}^{j'}(t_1, t_2)$.

### 4.2.1 Case of a Single Time

**Proposition 4.3.** For a single time $t \ge 0$ and $j' = 0 : d$, we have

$$ \tilde{C}^{j'}(t) = \sum_{j=1}^{d} \left[ (e_j + \tilde{M}^{j,(1)})^{\top} (e_j + \tilde{M}^{j,(1)}) + \tilde{C}^{j'} \right] * \phi_j^*(t). \quad (29) $$

Moreover for each $j', k, l \in 1 : d$, the vector of distinct covariance components $\tilde{C}_{k,l}^\bullet(t) = \top(\tilde{C}_{k,l}^1, \dots, \tilde{C}_{k,l}^d)(t)$ satisfies the linear integral equation

$$ \tilde{C}_{k,l}^{\bullet}(t) = R_{k,l}^{\bullet}(t) + \top \phi * \tilde{C}_{k,l}^{\bullet}(t) \quad (30) $$

whose solution is

$$ \tilde{C}_{k,l}^{\bullet}(t) = \top \left( \sum_{r=1}^{d} \phi^{(*)r} \right) * \tilde{R}_{k,l}^{\bullet}(t) \quad (31) $$

with components $\tilde{R}_{k,l}^j(u) = \top e_k (e_j + \tilde{M}^{j,(1)}(u))^{\top} (e_j + \tilde{M}^{j,(1)}(u)) e_l$ and $R_{k,l}^\bullet = \top \phi * \tilde{R}_{k,l}^\bullet$.

*Proof.* For a single time $t$ and a vector $a$, we have for $j' = 0 : d$

$$ \mathbf{E}_{(\phi^{j'}, \phi)} (\tilde{N}_k^{j'}(t) \tilde{N}_l^{j'}(t)) = \tilde{M}_{k,l}^{j',(2)}(t) = \partial_{a_k, a_l}^2 \tilde{L}^{j',(1)}(a, t) \Big|_{a=0} $$

Using Formula (4.2),

$$ \begin{aligned} \partial_{a_l}^1 \tilde{L}^{j',(1)}(a, t) &= \tilde{L}^{j',(1)}(a, t) \times U_l^{j'}(a, t) \\ \partial_{a_k, a_l}^2 \tilde{L}^{j',(1)}(a, t) &= \tilde{L}^{j',(1)}(a, t) [U_k^{j'}(a, t) U_l^{j'}(a, t) + \partial_{a_k}^1 U_l^{j'}(a, t)], \end{aligned} $$

with $U_l^{j'}(a, t) = \sum_{j=1}^d ((-1)_{l=j} e^{-a_j} \tilde{L}^{j,(1)}(a,.) + e^{-a_j} \partial_{a_l}^1 \tilde{L}^{j,(1)}(a,.) ) * \phi_j^{j'}(t)$. So, we have

$$ \begin{aligned} \partial_{a_k}^1 U_l^{j'}(a, t) &= \sum_{j=1}^d \left[ \partial_{a_k}^1 \left( (-1)_{l=j} e^{-a_j} \tilde{L}^{j,(1)}(a,.) + e^{-a_j} \partial_{a_l}^1 \tilde{L}^{j,(1)}(a,.) \right) \right] * \phi_j^{j'}(t) \\ &= \sum_{j=1}^d (-1)_{l=j} \left( (-1)_{k=j} e^{-a_j} \tilde{L}^{j,(1)}(a,.) + e^{-a_j} \partial_{a_k}^1 \tilde{L}^{j,(1)}(a,.) \right) * \phi_j^{j'}(t) \\ &\quad + \sum_{j=1}^d ((-1)_{k=j} e^{-a_j} \partial_{a_l}^1 \tilde{L}^{j,(1)}(a,.) + e^{-a_j} \partial_{a_k a_l}^2 \tilde{L}^{j,(1)}(a,.) ) * \phi_j^{j'}(t). \end{aligned} $$
Setting $a = 0$ and recalling that for all $s \ge 0$ and $j' = 0 : d$, we have

$$ \tilde{L}^{j', (1)}(0, s) \equiv 1, \partial_{a_k}^1 \tilde{L}^{j', (1)}(0, s) = U_k^{j'}(0, s) = -\tilde{M}_k^{j', (1)}(s), \text{ and } \partial_{a_k a_l}^2 \tilde{L}^{j', (1)}(0, s) = \tilde{M}_{k,l}^{j', (2)}(s), $$

we get

$$ \partial_{a_k}^1 U_l^{j'}(0, t) = \sum_{j=1}^d \left( 1_{l=j} 1_{k=j} + 1_{l=j} \tilde{M}_k^{j,(1)}(.) + 1_{k=j} \tilde{M}_l^{j,(1)}(.) + \tilde{M}_{kl}^{j,(2)}(.) \right) * \phi_j^{j'}(t) $$

and we finally obtain for $1 \le k, l \le d$

$$ \tilde{M}_{k,l}^{j', (2)}(t) = \left( \tilde{M}_k^{j', (1)} \tilde{M}_l^{j', (1)} + 1_{k=l} \Phi_k^{j'} + \tilde{M}_k^{l,(1)} * \phi_l^{j'} + \tilde{M}_l^{k,(1)} * \phi_k^{j'} + \sum_{j=1}^d \tilde{M}_{k,l}^{j,(2)} * \phi_j^{j'} \right) (t). $$

Next, let us consider the $d \times d$ components of the matrix covariance function at a single time $t$

$$ \tilde{C}_{k,l}^{j'}(t) = \tilde{M}_{k,l}^{j', (2)}(t) - \tilde{M}_k^{j', (1)}(t) \tilde{M}_l^{j', (1)}(t). $$

For $j' = 0 : d$, the latter expression can thus be written

$$ \tilde{C}_{k,l}^{j'}(t) = \sum_{j=1}^{d} \left[ (1_{j=k} + \tilde{M}_k^{j,(1)}) (1_{j=l} + \tilde{M}_l^{j,(1)}) + \tilde{C}_{k,l}^{j} \right] * \phi_j^{j'}(t). \quad (32) $$

Next, if we restrain index $j'$ to $\{1, ..., d\}$, we also observe vector integral equations involving $d$ functions, that is for each pair $(k, l)$, the component $\tilde{C}_{k,l}^{\bullet}(t) = \top(\tilde{C}_{k,l}^{1}, ..., \tilde{C}_{k,l}^{d})(t)$, satisfies an integral equation similar to that of the vector of means, that is

$$ \tilde{C}_{k,l}^{\bullet}(t) = R_{k,l}^{\bullet}(t) + \top \phi * \tilde{C}_{k,l}^{\bullet}(t). \quad (33) $$

Since $(\top \phi) * (\top \phi) = \top(\phi * \phi)$ and all functions are non negative, the solution, finite or infinite, corresponds to the series:

$$ \tilde{C}_{k,l}^{\bullet}(t) = \top \left( \sum_{r=0}^{\infty} \phi^{(*)r} \right) * R_{k,l}^{\bullet}(t). \quad (34) $$

Furthermore, we notice that Equation (32) can also be written

$$ \tilde{C}_{k,l}^{j'}(t) = \sum_{j=1}^{d} \left[ \top e_k (e_j + \tilde{M}^{j,(1)}) \top (e_j + \tilde{M}^{j,(1)}) e_l + \top e_k \tilde{C}^j e_l \right] * \phi_j^{j'}(t), \quad (35) $$

which takes the matrix form,

$$ \tilde{C}^{j'}(t) = \sum_{j=1}^{d} \left[ (e_j + \tilde{M}^{j,(1)})^{\top} (e_j + \tilde{M}^{j,(1)}) + \tilde{C}^{j} \right] * \phi_j^{j'}(t). \quad (36) $$

To end the proof, let us observe that $R_{k,l}^{\bullet} = \top \phi * \tilde{R}_{k,l}^{\bullet}$ with $\tilde{R}_{k,l}^{\bullet}$ having components $\tilde{R}_{k,l}^{j}(u) = \top e_k (e_j + \tilde{M}^{j,(1)}(s))^{\top} (e_j + \tilde{M}^{j,(1)}(s)) e_l$.

□
### 4.2.2 Covariance for Distinct Times

We now deal with the off diagonal covariance structures related to distinct times for indices $j' = 0 : d$

$$ \tilde{C}^{j'}(t_1, t_2) = \mathbf{E}_{(\phi^{j'}, \phi)} (\tilde{N}^{j'}(t_1)^{\top} \tilde{N}^{j'}(t_2)) - \tilde{M}^{j', (1)}(t_1)^{\top} \tilde{M}^{j', (1)}(t_2). $$

**Proposition 4.4.** For $t = (t_1, t_2)$ with $t_1 < t_2$ and convening again that $u_{[2]} = (u, u)$ if $u \in \mathbb{R}_+$, then for $j' = 0 : d$, the sub-covariance structure $\tilde{C}^{j'}(t)$ satisfies the matrix integral equation:

$$ \tilde{C}^{j'}(t) = \sum_{j=1}^{d} \int_{0}^{t_1} \left[ (e_j + \tilde{M}^{j,(1)}(t_1 - u))^{\top} (e_j + \tilde{M}^{j,(1)}(t_2 - u)) + \tilde{C}^{j}(t - u_{[2]}) \right] \phi_j^{j'}(u) du. \quad (37) $$

Moreover, if $\tilde{C}_{k,l}^{\bullet} = (\tilde{C}_{k,l}^{1}, \dots, \tilde{C}_{k,l}^{d})$ and $R_{k,l}^{\bullet} = (R_{k,l}^{1}, \dots, R_{k,l}^{d})$ with

$$ R_{k,l}^{j'}(t) = \sum_{j=1}^{d} \int_{0}^{t_1} \top e_k (e_j + \tilde{M}^{j,(1)}(t_1 - u))^{\top} (e_j + \tilde{M}^{j,(1)}(t_2 - u)) e_l \phi_j^{j'}(u) du, $$

then each index pair $(k, l)$ gives rise to independent vector integral equations

$$ \tilde{C}_{k,l}^{\bullet}(t) = R_{k,l}^{\bullet}(t) + \int_{0}^{t_1} \top \phi(u) \tilde{C}_{k,l}^{\bullet}(t - u_{[2]}) du. \quad (38) $$

*Proof.* We now consider the case $Q = \{(1, 1, k), (1, 2, l)\}$, that is

$$ \tilde{M}_{k,l}^{j', (2)}(t) = \partial_{a_k^1 a_l^2}^2 \tilde{L}^{j', (2)}(a, t) \Big|_{a=0} $$

with the Laplace transform $\tilde{L}^{j', (2)}(a, t)$ given by Equation (27). We have

$$ \partial_{a_l^2}^1 \partial_{a_k^1}^1 \tilde{L}^{j', (2)}(a, t) = \tilde{L}^{j', (2)}(a, t) \left( V_l^{j'}(a, t) U_k^{j'}(a, t) + \partial_{a_l^2}^1 U_k^{j'}(a, t) \right) $$

where

$$ U_k^{j'}(a, t) = \sum_{j=1}^{d} e^{-(a_j^1 + a_j^2)} \int_{0}^{t_1} (-1_{k=j} + \partial_{a_k^1}^1) \tilde{L}^{j', (2)}(a, (t - u_{[2]})) \phi_j^{j'}(u) du $$

satisfies when $a = 0$

$$ U_k^{j'}(0, t) = \sum_{j=1}^{d} \int_{0}^{t_1} (-1_{k=j} - \tilde{M}_k^{j,(1)}(t_1 - u)) \phi_j^{j'}(u) du = -\tilde{M}^{j,(1)}(t_1). $$

Moreover,

$$ \partial_{a_l^2}^1 U_k^{j'}(a, t) = \sum_{j=1}^{d} e^{-(a_j^1 + a_j^2)} \int_{0}^{t_1} \left( 1_{l=j} 1_{k=j} - 1_{k=j} \partial_{a_l^2}^1 - 1_{l=j} \partial_{a_k^1}^1 + \partial_{a_l^2 a_k^1}^2 \right) \tilde{L}^{j', (2)}(a, t - u_{[2]}) \\ \times \phi_j^{j'}(u) du $$

satisfies when $a = 0$

$$ \partial_{a_l^2}^1 U_k^{j'}(0, t) = \sum_{j=1}^{d} \int_{0}^{t_1} \left( 1_{l=j} 1_{k=j} + 1_{k=j} \tilde{M}_l^{j,(1)}(t_2 - u) + 1_{l=j} \tilde{M}_k^{j,(1)}(t_1 - u) + \tilde{M}_{k,l}^{j,(2)}(t - u_{[2]}) \right) \\ \times \phi_j^{j'}(u) du. $$
We also have

$$
\begin{aligned}
V_l^{j'}(a,t) = & \sum_{j=1}^{d} e^{-(a_j^1+a_j^2)} \int_0^{t_1} (-1_{l=j} + \partial_{a_l^1}^1) \tilde{L}^{j,(2)}(a, (t-u_{[2]})) \phi_j^{j'}(u) du \\
& + \sum_{j=1}^{d} e^{-a_j^2} \int_{t_1}^{t_2} (-1_{l=j} + \partial_{a_l^2}^1) \tilde{L}^{j,(1)}(a_{\bullet}^2, t_2-u) \phi_j^{j'}(u) du
\end{aligned}
$$

with

$$
\begin{aligned}
V_l^{j'}(0,t) &= \sum_{j=1}^{d} \left[ \int_{0}^{t_1} (-1_{l=j} - \tilde{M}_l^{j,(1)}(t_2-u)) \phi_j^{j'}(u) du \right. \\
&\qquad \left. + \int_{t_1}^{t_2} (-1_{l=j} - \tilde{M}_l^{j,(1)}(t_2-u)) \phi_j^{j'}(u) du \right] \\
&= -\tilde{M}_l^{j',(1)}(t_2).
\end{aligned}
$$

All in all, we finally get for $j' = 0 : d$

$$
\begin{aligned}
\tilde{M}_{k,l}^{j',(2)}(t) &= \tilde{M}_k^{j',(1)}(t_1) \tilde{M}_l^{j',(1)}(t_2) + \sum_{j=1}^{d} \int_{0}^{t_1} \left( \mathbf{1}_{l=j} \mathbf{1}_{k=j} + \mathbf{1}_{k=j} \tilde{M}_l^{(1),j}(t_2-u) \right. \\
&\qquad \left. + \mathbf{1}_{l=j} \tilde{M}_k^{j,(1)}(t_1-u) + \tilde{M}_{k,l}^{j,(2)}(t-u_{[2]}) \right) \phi_j^{j'}(u) du \\
&= \tilde{M}_k^{j',(1)}(t_1) \tilde{M}_l^{j',(1)}(t_2) + \sum_{j=1}^{d} \int_{0}^{t_1} {}^{\top}e_k \left( e_j {}^{\top}e_j + e_j {}^{\top}\tilde{M}^{j,(1)}(t_2-u) + \tilde{M}^{j,(1)}(t_1-u) {}^{\top}e_j \right. \\
&\qquad \left. + \tilde{M}^{j,(2)}(t-u_{[2]}) \right) e_l \phi_j^{j'}(u) du
\end{aligned}
$$

which can be written under the matrix form

$$
\tilde{M}^{j', (2)}(t) = \tilde{M}^{j', (1)}(t_1)^{\top} \tilde{M}^{j', (1)}(t_2) \\
+ \sum_{j=1}^{d} \int_{0}^{t_1} \left( e_j^{\top} e_j + e_j^{\top} \tilde{M}^{j', (1)}(t_2 - u) + \tilde{M}^{j', (1)}(t_1 - u)^{\top} e_j + \tilde{M}^{j', (2)}(t - u_{[2]}) \right) \phi_j^{j'}(u) du.
$$

Similarly, the covariance function can be written

$$
\tilde{C}^{j'}(t) = \sum_{j=1}^{d} \int_{0}^{t_1} \left[ (e_j + \tilde{M}^{j,(1)}(t_1 - u))^{\top} (e_j + \tilde{M}^{j,(1)}(t_2 - u)) + \tilde{C}^{j'}(t - u_{[2]}) \right] \phi_j^{j'}(u) du \quad (39)
$$

whose components satisfies for all $1 \le k, l \le d$ and $j' = 0 : d$

$$
\tilde{C}_{k,l}^{j'}(t) = R_{k,l}^{j'}(t) + \sum_{j=1}^{d} \int_{0}^{t_1} \tilde{C}_{k,l}^{j}(t - u_{[2]}) \phi_j^{j'}(u) du \quad (40)
$$

with

$$
R_{k,l}^{j'}(t) = \sum_{j=1}^{d} \int_{0}^{t_1} {}^{\top}e_k \left( e_j + \tilde{M}^{j,(1)}(t_1 - u) \right)^{\top} \left( e_j + \tilde{M}^{j,(1)}(t_2 - u) \right) e_l \phi_j^{j'}(u) du. \quad (41)
$$

Once again, the restriction of index $j'$ to 1: $d$, yields a set of $d$ matrix integral equations (39) related to the set of $d$ base matrices.
For each pair of indices $k,l$, let $\tilde{C}_{k,l}^{\bullet} = \top(\tilde{C}_{k,l}^{1}, \dots, \tilde{C}_{k,l}^{d})$ and $R_{k,l}^{\bullet} = \top(R_{k,l}^{1}, \dots, R_{k,l}^{d})$, then Equation (40) can be summed up to the following vector integral equation

$$ \tilde{C}_{k,l}^{\bullet}(t) = R_{k,l}^{\bullet}(t) + \int_{0}^{t_1} \top \phi(u) \tilde{C}_{k,l}^{\bullet}(t - u_{[2]}) du. \quad (42) $$

<figure><img src="image_3.png" /></figure>

### 4.2.3 Explicit Solution

Here, we provide an explicit solution of Equations (42). In this aim, we recourse again to generalized function calculus. Before formulating the proposition, we highlight that a real valued locally integrable function $\psi$ considered as a measure $\psi(u)du$ on $\mathbb{R}_+$, can be extended to a measure $\tilde{\psi}(du, dv) = \psi(u)\delta_u(dv)du$ supported by the diagonal $\Delta$ of $\mathbb{R}_+ \times \mathbb{R}_+$ and defined for any measurable function $f$ on $\mathbb{R}_+^2$ as follows

$$ \tilde{\psi}(f) = \int_{0}^{\infty} \int_{0}^{\infty} f(u, v) \psi(u) \delta_{u}(dv) du = \int_{0}^{\infty} f(u, u) \psi(u) du. $$

Moreover, by extension, we will consider the $d \times d$ matrix of positive measures $\tilde{\phi}(du, dv) = (\tilde{\phi}_{j}^{j'}(du, dv))_{j',j \in 1:d}$ supported by the diagonal $\Delta$ of $\mathbb{R}_+ \times \mathbb{R}_+$.

**Proposition 4.5.** For $t_1 < t_2$ and $k$ and $l$ in $1: d$, the vector $\tilde{C}_{k,l}^{\bullet}(t) = \top(\tilde{C}_{k,l}^{1}, \dots, \tilde{C}_{k,l}^{d})(t)$ satisfying Equation (42) has for solution

$$ \tilde{C}_{k,l}^{\bullet}(t_1, t_2) = \top \left( \sum_{r=1}^{\infty} \phi^{(*)r} \right) * \tilde{R}_{k,l}^{\bullet}(t_1, t_2) \quad (43) $$

where, for $1 \le j, k, l \le d$, the vector bivariate function $\tilde{R}_{k,l}^{\bullet}$ has components

$$ \tilde{R}_{k,l}^{j}(u, v) = \top e_k \left( e_j + \tilde{M}^{j,(1)}(u) \right) \top \left( e_j + \tilde{M}^{j,(1)}(v) \right) e_l. $$

Moreover, for $j' = 0$, that is to say $\phi^0 = \lambda^0$, we have

$$ \tilde{C}_{k,l}^{0}(t_1, t_2) = \top \left( \left( \sum_{r=0}^{\infty} \phi^{(*)r} \right) * \tilde{\lambda}^{0} \right) * \tilde{R}_{k,l}^{\bullet}(t_1, t_2). \quad (44) $$

*Proof.* With the previous convention related to convolution, we set

$$ \tilde{\phi}^{(*)0}(du, dv) = \text{diag}(\delta_{0|2}, \dots, \delta_{0|2}) = I_d \delta_{0|2}(du, dv), $$

which is of dimension $d \times d$, and $\tilde{\phi}^{(*)1} = \tilde{\phi}$. For higher orders, we obtain according to the definition of convolution of generalized functions:

$$ \begin{aligned} \tilde{\phi}_{j}^{j'} * \tilde{\phi}_{i}^{j'}(f) &= \int_{\mathbb{R}_{+}^{2}} \int_{\mathbb{R}_{+}^{2}} f(u_{1} + u_{2}, v_{1} + v_{2}) \phi_{j}^{j'}(u_{1}) \delta_{u_{1}}(dv_{1}) \phi_{i}^{i'}(u_{2}) \delta_{u_{2}}(dv_{2}) du_{1} du_{2} \\ &= \int_{\mathbb{R}_{+}^{2}} f(u_{1} + u_{2}, u_{1} + u_{2}) \phi_{j}^{j'}(u_{1}) \phi_{i}^{i'}(u_{2}) du_{1} du_{2} \\ &= \int_{0}^{\infty} f(w, w) \left( \int_{0}^{\infty} \phi_{j}^{j'}(u_{1}) \phi_{i}^{i'}(w - u_{1}) du_{1} \right) dw \\ &= \int_{0}^{\infty} f(w, w) \left( \int_{0}^{w} \phi_{j}^{j'}(u_{1}) \phi_{i}^{i'}(w - u_{1}) du_{1} \right) dw \\ &= \tilde{\phi}_{j}^{j'} * \phi_{i}^{i'}(f). \end{aligned} $$
Similarly, combining matrix product properties and iterating the latter result, we get $(\tilde{\phi})^{(*)r} = (\tilde{\phi}^{(*)r})$ for all $r \ge 2$.

Then, let us observe that

$$ \begin{aligned} \int_0^{t_1} \tilde{C}_{k,l}^j (t - u_{[2]}) \phi_j^{j'}(u) du &= \int_0^{t_1} \int_0^{t_2} \tilde{C}_{k,l}^j (t - (u, v)) \phi_j^{j'}(u) \delta_u(dv) du \\ &= \tilde{C}_{k,l}^j * \tilde{\phi}_j^{j'}(t_1, t_2). \end{aligned} $$

Similarly, using tensor product $\otimes$ for functions, we get from Equation (41)

$$ \begin{aligned} R_{k,l}^{j'}(t_1, t_2) &= \sum_{j=1}^{d} \int_{0}^{t_1} \int_{0}^{t_2} {}^{\top}e_k \left( e_j + \tilde{M}^{j,(1)}(t_1 - u) \right) {}^{\top} \left( e_j + \tilde{M}^{j,(1)}(t_2 - v) \right) e_l \phi_j^{j'}(u) \delta_u(dv) du \\ &= \sum_{j=1}^{d} \tilde{R}_{k,l}^j * \tilde{\phi}_j^{j'}(t_1, t_2) \end{aligned} $$

with,

$$ \tilde{R}_{k,l}^j(u, v) = {}^{\top}e_k \left( e_j + \tilde{M}^{j,(1)}(u) \right) \otimes {}^{\top} \left( e_j + \tilde{M}^{j,(1)}(v) \right) e_l. $$

To sum up, the vector $\tilde{C}_{k,l}^\bullet = (\tilde{C}_{k,l}^1, \dots, \tilde{C}_{k,l}^d)$ satisfies

$$ \tilde{C}_{k,l}^{\bullet}(t_1, t_2) = R_{k,l}^{\bullet}(t_1, t_2) + {}^{\top}\tilde{\phi} * \tilde{C}_{k,l}^{\bullet}(t_1, t_2) $$

and has for solution:

$$ \tilde{C}_{k,l}^{\bullet} = {}^{\top}\left(\sum_{r=0}^{\infty} \tilde{\phi}^{(*)r}\right) * R_{k,l}^{\bullet} = {}^{\top}\left(\sum_{r=1}^{\infty} \tilde{\phi}^{(*)r}\right) * \tilde{R}_{k,l}^{\bullet}, $$

with $\tilde{R}_{k,l}^\bullet = {}^{\top}(\tilde{R}_{k,l}^1, \dots, \tilde{R}_{k,l}^d)$ and $R_{k,l}^\bullet = {}^{\top}\tilde{\phi} * \tilde{R}_{k,l}^\bullet$.

Finally, for $j' = 0$, that is letting $\phi^0 = \lambda^0$, we get from Equation (39)

$$ \begin{aligned} \tilde{C}_{k,l}^0 &= {}^{\top}\tilde{\lambda}^0 * \tilde{R}_{k,l}^\bullet + {}^{\top}\tilde{\lambda}^0 * {}^{\top}\left(\sum_{r=1}^{\infty} \tilde{\phi}^{(*)r}\right) * \tilde{R}_{k,l}^\bullet \\ &= {}^{\top}\tilde{\lambda}^0 * {}^{\top}\left(\sum_{r=0}^{\infty} \tilde{\phi}^{(*)r}\right) * \tilde{R}_{k,l}^\bullet. \end{aligned} $$

□

## Remarks

* Each component $\tilde{C}_{k,l}^{j'}$ or $\tilde{M}_{k,l}^{(2)j'}$ satisfies an integral equation only including components of indices $k$ and $l$, such that this equation can be solved independently from the other equations corresponding to other pairs of indices.

* We can consider that $(\sum_{r=0}^{\infty} \tilde{\phi}^{(*)r})$ is the inverse (with respect to both matrix product and convolution product of distribution functions) of the signed matrix valued measure supported by the diagonal $\Delta$ of $\mathbb{R}_+^2$, namely $I_d\delta_{0[2]} - \tilde{\phi}$.

* The covariance measures $C_{k,l}^{j'}$, $j' = 0: d$ and $1 \le k, l \le d$, are all positive measures and, hence, components $N_k^{j'}(t_1)$ and $N_l^{j'}(t_2)$ are always non negatively correlated.

* Note finally that for $t_1 = t_2$, Solution (43) reduces to Solution (31) given for a single time.
### 4.2.4 Example

As in the example related to the mean equations and described in Section 4.1.2, let us take the $\Gamma$ intensities with $\phi(t) = Wf(t; \kappa, \theta)$ whose primitive function is $\Phi(t) = W\gamma(\kappa, t/\theta)/\Gamma(\kappa)$. For $\kappa = 1$, we readily get

$$ \top\left(\sum_{r=1}^{\infty} \widetilde{\phi}^{(*)r}(u, v)\right)(du, dv) = \frac{e^{-u/\theta}}{u} \left(e^{Wu/\theta} - I_d\right) \delta_u(dv)du. $$

On the other hand, with $\tilde{R}_{k,l}^j(u, v) = \top e_k(I_d + \tilde{M}^{(1)}(u)) e_j \top e_j^\top(I_d + \tilde{M}^{(1)}(v)) e_l$ and $\tilde{M}^{(1)}(u) = \sum_{r=1}^\infty \frac{W^r}{r!} \gamma(r, u/\theta)$, we get after transposition operations

$$ \tilde{C}_{k,l}^{j'}(t_1, t_2) = \sum_{j=1}^{d} \int_{0}^{t_1} \frac{e^{-u/\theta}}{u} \top e_{j'}(e^{Wu/\theta} - I_d) e_j \tilde{R}_{k,l}^j(t_1 - u, t_2 - u) du $$

which has no known analytic expression but which can be numerically approximated.

## 4.3 Lebesgue Decomposition of the Covariance Measure

The covariance functions $C$ of a point process $N$ on a measurable space $(D, \mathcal{D})$ defined by $C(A \times B) = \mathbf{E}(N(A)N(B)) - \mathbf{E}(N(A))\mathbf{E}(N(B))$, $A, B \in \mathcal{D}$, is advantageously considered from a statistical point of view as a signed moment measure of order 2, that is $C(du, dv)$ on $D^2$ or simply $C(t_1, t_2) = C([0, t_1] \times [0, t_2])$. It is well known that it has at least a singular part $^S C$ whose support is the diagonal $\Delta = \{u|_2 = (u, u), u \in D\}$ of $D^2$. The singular and continuous parts may take different forms depending on the probability distribution of $N$. We develop this issue for our multivariate Hawkes process and prove that the singular part is also supported by $\Delta$ but has an additional part to the usual first order moment. The absolutely continuous part is also determined. For that purpose we recourse once again to a measure and generalized function context.

**Proposition 4.6.** Let the matrix valued function be continuous and non negative. Then the covariance measure $\tilde{C}_{k,l}^{j'}$ is locally finite, it has for all $j' = 0: d$ and $1 \le k, l \le d$, an absolute continuous part, namely

$$ \begin{aligned} {}^{AC}\tilde{C}_{k,l}^{j'}(du, dv) = & \left(\tilde{m}_k^{j',(1)}(u) \tilde{m}_l^{k,(1)}(v-u)\right. \\ & \left. + \sum_{j=1}^{d} \int_{0}^{u\wedge v} \tilde{m}_j^{j',(1)}(w) \tilde{m}_k^{j,(1)}(u-w) \tilde{m}_l^{j,(1)}(v-w) dw\right) du dv. \end{aligned} $$

The singular part is supported by the diagonal $\Delta$ of $\mathbb{R}_+^2$ and is written:

$$ {}^S\tilde{C}_{k,l}^{j'}(du, dv) = \left( \delta_{k=l} \tilde{m}_k^{j',(1)} + \tilde{m}_l^{j',(1)} * \tilde{m}_k^{l,(1)} \right) (u) \delta_u(dv) du. $$

*Proof.* We first recall some results on differential calculus as well as some notations related to generalized functions with support on $\mathbb{R}_+^2$.

The Heaviside function $H(u, v) = 1_{\mathbb{R}_+^2}(u, v)$ considered as a distribution (precisely the restricted Lebesgue measure) satisfies the equation $\partial_u \partial_v H(du, dv) = \delta_0(du)\delta_0(dv) = \delta_{0|2}(du, dv)$.

Let us assume firstly that $j' \ne 0$. The convergent series $h(u) = (\sum_{r=1}^\infty \phi^{(*)r})(u) = \tilde{m}^{(1)}(u)$ if finite, is necessarily continuous since it is written as the convolution of a generalized function
with a continuous one, that is to say $h = (\sum_{r=0}^{\infty} \phi^{(*)r}) * \phi$, inheriting therefore the continuity property of $\phi$.

Then, emphasizing that the constant 1 is actually a function $1(u) = 1_{[0,\infty]}(u)$, we can develop $\tilde{R}_{k,l}^j$ as follows

$$ \tilde{R}_{k,l}^j(u, v) = \delta_{k=l} 1(u) 1(v) + \delta_{j=l} \tilde{M}_k^{l,(1)}(u) 1(v) + \delta_{j=k} 1(u) \tilde{M}_l^{k,(1)}(v) + \tilde{M}_k^{j,(1)}(u) \tilde{M}_l^{j,(1)}(v). $$

Moreover, according to Proposition 4.5 with $\tilde{h}(du, dv) = h(u)\delta_u(dv)du$, we have

$$ \tilde{C}_{k,l}^{j'}(t_1, t_2) = \left( \sum_{j=1}^{d} \tilde{h}_j^{j'} * \tilde{R}_{k,l}^{j} \right) (t_1, t_2). $$

Now, we separately compute each term of $\tilde{h}_j^{j'} * \tilde{R}_{k,l}^j$ using the expression of $\tilde{R}_{k,l}^j$ provided above. We recall that $\tilde{m}^{(1)}$ has support $\mathbb{R}_+$, and we first get

$$ \begin{aligned} \left( \tilde{h}_j^{j', (1)} * [1(\cdot)1(\cdot)] \right) (t_1, t_2) &= \int_0^{t_1} \tilde{m}_j^{j', (1)}(u) du \\ \left( \tilde{h}_j^{j'} * [\tilde{M}_k^{l,(1)}(\cdot)1(\cdot)] \right) (t_1, t_2) &= \int_0^{t_1} h_j^{j'}(u) \tilde{M}_k^{l,(1)}(t_1 - u) du \\ &= \int_0^{t_1} \int_0^{t_1} \tilde{m}_j^{j', (1)}(u) \tilde{m}_k^{l,(1)}(w-u) du dw \\ &= \int_0^{t_1} \left( \int_0^w \tilde{m}_j^{j', (1)}(u) \tilde{m}_k^{l,(1)}(w-u) du \right) dw \\ &= \int_0^{t_1} \tilde{m}_j^{j', (1)} * \tilde{m}_k^{l,(1)}(w) dw. \end{aligned} $$

Let us observe that both components depend only on the first variable $t_1$ meaning that both measures have the diagonal $\Delta$ as support.

For the third term, we obtain

$$ \begin{aligned} \left( \tilde{h}_j^{j'} * [1(\cdot)\tilde{M}_l^{k,(1)}(\cdot)] \right) (t_1, t_2) &= \int_0^{t_1} h_j^{j'}(u) \tilde{M}_l^{k,(1)}(t_2 - u) du \\ &= \int_0^{t_1} \int_0^{t_2} \tilde{m}_j^{j', (1)}(u) \tilde{m}_l^{k,(1)}(v - u) du dv. \end{aligned} $$

This means that this measure component is absolutely continuous with respect to the Lebesgue measure on $\mathbb{R}_+^2$ with density $\tilde{m}_j^{j', (1)}(u) \tilde{m}_l^{k,(1)}(v-u)$.

For the last term, we get

$$ \begin{aligned} \left( \tilde{h}_j^{j'} * [\tilde{M}_k^{j,(1)}(\cdot) \tilde{M}_l^{j,(1)}(\cdot)] \right) (t_1, t_2) &= \int_0^{t_1} h_j^{j'}(u) \tilde{M}_k^{j,(1)}(t_1 - u) \tilde{M}_l^{j,(1)}(t_2 - u) du \\ &= \int_0^{t_1} \int_0^{t_1} \int_0^{t_2} \tilde{m}_j^{j', (1)}(u) \tilde{m}_k^{j,(1)}(w-u) \\ &\qquad \times \tilde{m}_l^{j,(1)}(v-u) du dw dv \\ &= \int_0^{t_1} \int_0^{t_2} \left( \int_0^{v \wedge w} \tilde{m}_j^{j', (1)}(u) \tilde{m}_k^{j,(1)}(w-u) \right. \\ &\qquad \left. \times \tilde{m}_l^{j,(1)}(v-u) du \right) dw dv. \end{aligned} $$
Therefore this last component is absolutely continuous with density

$$ \int_{0}^{v \wedge w} \tilde{m}_{j}^{j', (1)}(u) \tilde{m}_{k}^{j, (1)}(w-u) \tilde{m}_{l}^{j, (1)}(v-u) du. $$

To sum up, $\tilde{C}^{j'}(du, dv)$ has a singular component supported by the diagonal $\Delta$, namely

$$ S\tilde{C}_{k,l}^{j'}(du, dv) = \left( \delta_{k=l} \tilde{m}_{k}^{j'} + \tilde{m}_{l}^{j'} * \tilde{m}_{k}^{l, (1)} \right) (u) \delta_{u}(dv) du, $$

and the absolute continuous component is written

$$ \begin{aligned} AC\tilde{C}_{k,l}^{j'}(du, dv) = & \left( \tilde{m}_{k}^{j', (1)}(u) \tilde{m}_{l}^{k, (1)}(v-u) \right. \\ & \left. + \sum_{j=1}^{d} \int_{0}^{u \wedge v} \tilde{m}_{j}^{j', (1)}(w) \tilde{m}_{k}^{j, (1)}(u-w) \tilde{m}_{l}^{j, (1)}(v-w) dw \right) du dv. \end{aligned} $$

If $j' = 0$, we get the same formulas by defining $h = (\sum_{r=0}^{+\infty} \phi^{(*)r}) * \lambda^0$ which is still continuous as soon as $\lambda^0$ is continuous. $\square$

## Remarks

* It is worth noticing that the singular and continuous parts satisfy the same integral equation when considering the associated components of $\tilde{R}_{k,l}^j$.

* Fundamentally, the first and second moments roughly satisfy the same type of integral equation, save that the space dimension increases. It is sensibly expected that the moments of higher order involving several times, would satisfy similar integral equations with explicit solutions having more complex singular parts supported by one dimensional diagonal and semi diagonal hyperplanes, etc.

# 5 Schemes for Numerical Approximation

In this section, we describe some algorithms for computing the numerical solutions of the equations of the first two moments and the Laplace transforms developed in the previous sections. Since these functions are time-dependent, we use the following discretization scheme: For $T > 0$, we consider the regular partition of the time interval $[0, T]$ consisting of successive sub-intervals with length $\tau = 1/M > 0$, $M \in \mathbb{N}_*$, that is $0 < \tau < 2\tau < ... < M\tau = T$, and we call it the grid denoted by $G = \{0, ..., M\}$.

## 5.1 Time Approximation of Laplace Transforms

Since the solution of the Laplace transform of any Hawkes process driven by $(\lambda^0, \phi)$, is simply expressible as the solution of the basic system of Volterra-like equations of Theorem 2.2 (see the corresponding remark), we propose the algorithm described below to approach this basic system on any time interval $[0, T]$. We further assume that the multi-time Laplace transforms $\tilde{L}^{j', (k)}(a, t)$ are continuous with respect to the time variable $t$, such that the limits of their discrete time approximation remain meaningful. We also assume that the matrix parameter $a$ is fixed and, hence, we often omit it in the sequel.
### 5.1.1 Single-time Laplace Transform Approximation

For any fixed vector parameter $\top a = (a_1^1, \dots, a_d^1)$, we propose to approximate the values of $\tilde{L}^{j',(1)}(a_\bullet^1, \cdot)$, $j' = 1: d$, on the grid $\mathcal{G}$, that is $\bar{L}_a^{j',(1)}(m) \approx \tilde{L}^{j',(1)}(a, m\tau)$. From Equation (7) in Theorem 2.2 restricted to the case $k = 1$, and after a time inversion, we get, for $j' = 0: d$,

$$ \tilde{L}^{j',(1)}(a, t) = \exp \left( \sum_{j=1}^{d} \int_{0}^{t} \left( e^{-a_j} \tilde{L}^{j,(1)}(a, u) - 1 \right) \phi_j^{j'}(t-u) \, du \right). $$

Hence, we trivially set, for all $j' = 1: d$,

$$ \bar{L}_a^{j',(1)}(0) = 1 $$

$$ \bar{L}_a^{j',(1)}(1) = \exp \left( \tau \sum_{j=1}^{d} \left( e^{-a_j} \bar{L}_a^{j,(1)}(0) - 1 \right) \phi_j^{j'}(\tau) \right) $$

$$ \bar{L}_a^{j',(1)}(2) = \exp \left( \tau \sum_{j=1}^{d} \left( e^{-a_j} \bar{L}_a^{j,(1)}(0) - 1 \right) \phi_j^{j'}(2\tau) + \left( e^{-a_j} \bar{L}_a^{j,(1)}(1) - 1 \right) \phi_j^{j'}(\tau) \right), $$

and for $m = 3: M$,

$$ \bar{L}_a^{j',(1)}(m) = \exp \left( \tau \sum_{j=1}^{d} \sum_{r=0}^{m-1} \left( e^{a_j} \bar{L}_a^{j,(1)}(r) - 1 \right) \phi_j^{j'}((m-r)\tau) \right). \quad (45) $$

### 5.1.2 Two-times Laplace Transform Approximation

In this case, the general formula in Theorem 2.2 with the $d \times 2$-matrix $a = (a_\bullet^1, a_\bullet^2)$ and the time vector $t = (t_1, t_2)$ becomes, for $j' = 0: d$,

$$ \tilde{L}^{j',(2)}(a, t) = \exp \left( \sum_{j=1}^{d} \left( \int_{0}^{t_1} \left[ e^{-(a_j^1+a_j^2)} \tilde{L}^{j,(2)}(a, (t_1-u, t_2-u)) - 1 \right] \phi_j^{j'}(u) du \right. \right. \\ \left. + \int_{t_1}^{t_2} \left[ e^{-a_j^2} \tilde{L}^{j,(1)}(a_\bullet^2, (t_2-u)) - 1 \right] \phi_j^{j'}(u) du \right) \right). $$

In the same way as in the previous subsection, for any fixed $a$, we propose to approximate $\tilde{L}^{j',(2)}(a, \cdot)$ on the nodes $(m_1\tau, m_2\tau)$, $0 \le m_1 < m_2 < M$, of the bi-dimensional grid $\mathcal{G}^2$ by $\bar{L}_a^{j',(2)}(m_1, m_2)$.

Considering the previous equation, we have essentially to approximate $\bar{L}_a^{j',(2)}(m_1, m_2)$ recursively only along secondary diagonals $(m_1 - r, m_2 - r)$, $r = 0: m_1$. Clearly, we have for all $j' = 1: d$ and all $m = 0: M$,

$$ \bar{L}_a^{j',(2)}(0, m) = \bar{L}_{a^2}^{j',(1)}(m). $$

Hence, after the time inversion $v = t_1 - u$, we propose the following recursive formulas:

$$ \bar{L}_a^{j',(2)}(1, m+1) = \exp \left( \tau \sum_{j=1}^{d} \left( \left[ e^{-(a_j^1+a_j^2)} \bar{L}_a^{j,(2)}(0, m) - 1 \right] \phi_j^{j'}(\tau) \right. \right. \\ \left. + \sum_{r'=1}^{m} \left[ e^{-a_j^2} \bar{L}_{a^2}^{j,(1)}(r') - 1 \right] \phi_j^{j'}((m+1-r')\tau) \right) \right), $$
and for $1 < r \le m$

$$ \bar{L}_{a}^{j', (2)}(r, m+r) = \exp \left( \tau \sum_{j=1}^{d} \left( \sum_{r'=0}^{r-1} \left[ e^{-(a_j^1+a_j^2)} \bar{L}_{a}^{j', (2)}(r', m+r') - 1 \right] \phi_j^{j'}((r-r')\tau) \right. \right. \\ \left. \left. + \sum_{r'=r}^{m+r} \left[ e^{-a_j^2} \bar{L}_{a\bullet}^{j', (1)}(r') - 1 \right] \phi_j^{j'}((m+r-r')\tau) \right) \right). \quad (46) $$

Thus, using the latter formula, for any $0 < m_1 < m_2 \le M$ and $j' = 1 : d$, we finally get:

$$ \bar{L}_{a}^{j', (2)}(m_1, m_2) = \bar{L}_{a}^{j', (2)}(m_1, (m_2 - m_1) + m_1). $$

### 5.1.3 Remarks

1. One may consider that the previous approximations $\bar{L}_{a}^{j', (k)}$, $k = 1 : 2$, simply constitute the first step to initialize an iterative process using the implicit integral equation of Theorem 2.2 to yield successive approximations, using the well known fixed-point theorem.

2. Using the same fixed-point theorem, one may as well initialize the iteration processes with the well known analytic form of the Laplace transforms of the d multi-type Poisson processes each driven by density $\phi_j^{j'}$, $j' = 1 : d$, and compare the two approaches.

3. For a general Hawkes process driven by $(\lambda^0, \phi)$, we have only to substitute $\phi_j^{j'}$ for $\lambda_j^0$ in Equations (45) and (46)

4. Even tedious, similar procedures of numerical approximations can be applied to higher-order multi-time Laplace transforms $\tilde{L}^{j', (k)}(a, (t_1, \dots, t_n))$, $2 < k \le n$.

## 5.2 Approximation of the First Two Moment Functions

With regard to the calculation of the mean and covariance functions of Hawkes processes, we have to essentially deal with the expression of the fundamental matrix function $\sum_{r=0}^{\infty} \phi^{(*)r}$. For specific models, the elements of the matrix function $\phi = (\phi_j^{j'})$ can be easily convoluted analytically with each other at any order, and even yield an explicit expression for the series; see Examples 4.1.2 and 4.2.4. Otherwise, the different mean and covariance formulas in Equations (26), (43) and (44) can be simply approximated on any given grid $\mathcal{G}$ by limiting the fundamental matrix series to $\sum_{r=0}^{K} \phi^{(*)r}$ with a large enough integer $K$.

Indeed, plain numerical methods can be implemented for convolution product of functions. For example, if $f_1$ and $f_2$ defined on $\mathbb{R}^+$ have values $f_1(m_1\tau)$ and $f_2(m_2\tau)$ on the grid $\mathcal{G}$ then the values of $f_3 = f_1 * f_2$ can be approximated by $\tilde{f}_3(m\tau) = \tau \sum_{r=1}^{m} f_1(r\tau)f_2((m+1-r)\tau)$ on the same grid.

This algorithm extends easily to convolution products of matrix functions. Note also that the expression $(\sum_{r=0}^{\infty} \phi^{(*)r})(t_1, t_2)$ related to the covariance solution in Proposition 4.5 mainly consists in restraining integration over the diagonal of $\mathbb{R}_+^2$.

Even if all calculation procedures presented above deal with generic multidimensional Hawkes process, what follows illustrate these procedures for a specific model of a 2D-Hawkes process, for which we give some numerical results about its first two moments.

### 5.2.1 Specification of a 2D-Hawkes process

Thus, let us consider the following 2D-Hawkes process driven by the functional parameter $(\lambda^0, \phi)$, where the baseline intensity $\lambda^0$ (that may also be called the ignition or immigration intensity) is periodic, and each component of the basic matrix function $\phi$ is a beta-like distribution,
more precisely:

$$ 
\begin{aligned}
\lambda_i^0(t) &= a_i + b_i \sin(c_i t), \\
\phi_j^{j'}(t) &= \alpha_j^{j'} t^{\beta_j^{j'}} [\gamma_j^{j'} - t]^{\rho_j^{j'}} \quad \text{if } t \le \gamma_j^{j'} \text{ and } 0 \text{ otherwise,}
\end{aligned} 
\qquad (47) 
$$

where all parameters are real numbers, $i = 1:2$, $a_i > 0$, $a_i > b_i$ and $1 \le j', j \le 2$. Observe that $(\lambda^0, \phi)$ depends on a finite number of scalar parameters and thus is considered as a classical parametric model for the sequel. For illustration purpose, parameters $a, b, c, \alpha, \beta, \gamma, \rho$ were randomly drawn, yielding:

<table><thead><tr><th>$i$ or $j$</th><th>$a$</th><th>$b$</th><th>$c$</th><th>$\alpha_{\bullet}^{1}$</th><th>$\alpha_{\bullet}^{2}$</th><th>$\beta_{\bullet}^{1}$</th><th>$\beta_{\bullet}^{2}$</th><th>$\gamma_{\bullet}^{1}$</th><th>$\gamma_{\bullet}^{2}$</th><th>$\rho_{\bullet}^{1}$</th><th>$\rho_{\bullet}^{2}$</th></tr></thead><tbody><tr><td>1</td><td>1.057</td><td>0.031</td><td>0.845</td><td>0.073</td><td>0.046</td><td>0.060</td><td>1.254</td><td>1.576</td><td>1.831</td><td>0.598</td><td>0.897</td></tr><tr><td>2</td><td>1.061</td><td>0.093</td><td>0.817</td><td>0.050</td><td>0.096</td><td>1.897</td><td>1.923</td><td>0.369</td><td>0.182</td><td>0.789</td><td>0.713</td></tr></tbody></table>

The whole time interval is defined with $T = 10$, the mesh length is $\tau = 0.005$ and the infinite sum is stopped at $r = 500$.

Figure 1 shows the shapes of the functional parameters constituting $\lambda^0$ and $\phi$. The periodic shape of $\lambda^0$ was chosen to picture a potential seasonality effect in the ignition of epidemics. The different beta-like shapes of $\phi$ components presented here may represent differences in the characteristics of the epidemic dynamics in and between different regions (e.g., heterogeneity in latency, virulence, infection duration, spatial connectivity, properties of pathogen strains, social and care conditions, etc.). Figure 2 shows the behavior of the fundamental series, and illustrates in this case the rapid extinction of the epidemic if no ignition / immigration occurs.

<figure><img src="image_4.png" /><figcaption>Figure 1: Functional parameters constituting $\lambda^0$ and $\phi$ specified in Equation (47), and drawn over the time interval [0, 10]. Left panel: Ignition intensities $\lambda_j^0(t) \equiv \phi_j^0(t)$, $j = 1, 2$. Right panel: Basic 2×2-matrix functions $\phi_j^{j'}(t)$, $1 \le j', j \le 2$. Black and grey colors correspond to column $j' = 1$, red and orange to column $j' = 2$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 2: Fundamental series $\sum_{r=1}^{\infty} \phi^{(*)r}(t)$, for $t \in [0, 10]$, approximated up to order $r = 500$. Black and grey colors correspond to column $j' = 1$, red and orange to column $j' = 2$.</figcaption></figure>

### 5.2.2 Calculation of the Mean functions

The left panel of Figure 3 gives the basic intensity ratios for $1 \le j, j' \le 2$, that is to say: $\tilde{m}_j^{j', (1)}(t) / (\tilde{m}_1^{j', (1)}(t) + \tilde{m}_2^{j', (1)}(t))$. We observe the fast convergence of the component ratios towards two limits. This phenomenon is also observed for all models we have tried out and is probably related to Perron theorem on positive matrix. The right panel shows the basic mean functions $\tilde{M}_j^{j', (1)}(t)$. We observe in this example finite asymptotic limits in the absence of an ignition / immigration process.

Figure 4 shows, on the left, a subtle effect of the periodicity of the baseline intensity on the behavior of the intensities $\tilde{m}_j^{0,(1)}$, $j = 1:2$, and, on the right, the increasing mean components $\tilde{M}_j^{0,(1)}$, $j = 1:2$, of the original Hawkes process driven by $(\lambda^0, \phi)$.

<figure><img src="image_6.png" /><figcaption>Figure 3: Ratios of basic intensities (left panel) and corresponding cumulative intensities, i.e., the mean functions (right panel), over the time interval $[0, 10]$. Black and grey colors hold for $j' = 1$ while red and orange hold for $j' = 2$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: First order moments of the 2D-Hawkes process driven by the functional parameter $(\lambda^0, \phi)$, over the time interval $[0, 10]$. Left panel: The two components of the intensities $\tilde{m}_j^{0,(1)}$, $j = 1:2$. Right panel: The correspondent cumulative intensities, that is the mean functions $\tilde{M}_j^{0,(1)}$, $j = 1:2$.</figcaption></figure>

### 5.2.3 Calculation of the Covariance Structure

From the algorithmic point of view, the computation of the covariance structure is very similar to that of the mean. It is however much longer since the covariance function depends on two variates. Figure 5 shows two functional elements of the fundamental covariance structure, namely the inner-region variance $\tilde{C}_{1,1}^{j'}(t, t)$ (left) and the inter-regions covariance $\tilde{C}_{1,2}^{j'}(t, t)$ (right). Once again, due to the beta like form of $\phi$, we observe asymptotic finite values pointing out the extinction in the absence of immigration. Figure 6 represents images with level contours of two components of the $2 \times 2$ matrix function of correlations, which depends on 2 variables.

<figure><img src="image_8.png" /><figcaption>Figure 5: Example of basic covariance structure over the time interval $[0, 10]$, for $j' = 1$. Left panel: Variance function $t \mapsto \tilde{C}_{1,1}^{j'}(t, t)$. Right panel: Covariance function $t \mapsto \tilde{C}_{1,2}^{j'}(t, t)$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6: Continued example. Basic correlation structure of a point process driven by $(\phi_{\bullet}^{1}, \phi)$. Left panel: Correlation within the region $j' = 1$ obtained by standardization of $\tilde{C}_{1,1}^{1}(t_1, t_2)$, $0 \le t_1, t_2 \le 10$. Right panel: Correlation between regions $j' = 1$ and $j = 2$ obtained by standardization of $\tilde{C}_{1,2}^{1}(t_1, t_2)$, $0 \le t_1, t_2 \le 10$.</figcaption></figure>

Figure 7 shows similar graphs related to the original Hawkes process (i.e., when the ignition/immigration process is accounted for). The effect of the periodic trend of $\lambda^0$ can be slightly observed through the inner-region correlation structure (left; see the subtle fluctuations of the isolines). The inter-regions covariance structure is displayed in the right panel.

<figure><img src="image_10.png" /><figcaption>Figure 7: Correlation structures of the Hawkes process driven by $(\lambda^0, \phi)$, with periodic $\lambda^0$. Left panel: Correlation within the region $j' = 2$, obtained from $\tilde{C}_{2,2}^{0}(t_1, t_2)$, $0 \le t_1, t_2 \le 10$. Right panel: Correlation between regions $j' = 1$ and $j = 2$ obtained from $\tilde{C}_{2,1}^{0}(t_1, t_2)$, $0 \le t_1, t_2 \le 10$.</figcaption></figure>

# 6 Trajectories of Multivariate Hawkes Processes

In this section, we present simulated trajectories for different multivariate Hawkes processes in order to illustrate the various behaviors that the Hawkes Process can embrace.
First, we look at the effect of the excitation function by considering exponential, gamma, constant and beta excitation functions with constant baseline intensity.

**Case 1: Exponential Excitation Function with Constant Baseline Intensity.** In this case, the excitation function is exponential, and the baseline intensity is constant. The intensity function $\lambda(t)$ at time $t$ is given by:

$$ \lambda(t) = \mu + \sum_{i=1}^{n} \alpha_i e^{-\beta(t-t_i)^{\gamma}} $$

where $\mu$ is the constant baseline intensity, $\alpha_i$ is the excitation parameter for the $i$-th process, $\beta$ is the decay rate of the exponential excitation, and $(t-t_i)^{\gamma}$ denotes the positive part of $(t-t_i)$, which is zero for $t < t_i$ and $t - t_i$ for $t \ge t_i$. Figure 8 displays two contrasted realizations of trajectories for this setup. The 3-dimensional counting process $N$ is displayed as well as the evolution of the 3-dimensional intensity function.

### Realization 1

<figure><img src="image_11.png" /></figure>

### Realization 2

<figure><img src="image_12.png" /><figcaption>Figure 8: Trajectory of the counting process and evolution of the intensity function for the exponential excitation function with a constant baseline intensity (Case 1).</figcaption></figure>

**Case 2: Gamma Excitation Function with Constant Baseline Intensity.** In this case, the excitation function corresponds to a gamma distribution, and the baseline intensity remains constant. The intensity function $\lambda(t)$ is given by:

$$ \lambda(t) = \mu + \sum_{i=1}^{n} \alpha_i \left( \frac{(t - t_i)^{\gamma-1}}{\Gamma(\gamma)} \right) e^{-\delta(t - t_i)^{\gamma}} $$
where $\mu$ is the baseline intensity, $\alpha_i$ is the excitation parameter for the $i$-th process, $\gamma$ is the shape parameter of the gamma distribution, $\delta$ is the rate parameter, and $(t - t_i)^+$ is the positive part function as described previously. An example of trajectories is shown in Figure 9.

<figure><img src="image_13.png" /><figcaption>Figure 9: Trajectory of the counting process and evolution of the intensity function for the gamma excitation function with a constant baseline intensity (Case 2).</figcaption></figure>

**Case 3: Constant Excitation Function with Constant Baseline Intensity.** In this case, the excitation function is constant, and the baseline intensity is also constant. The intensity function $\lambda(t)$ is given by:

$$ \lambda(t) = \mu + \sum_{i=1}^{n} \alpha_i $$

where $\mu$ is the constant baseline intensity, and $\alpha_i$ is the constant excitation parameter for the $i$-th process. Figure 10 gives a realization of trajectories in this case.

<figure><img src="image_14.png" /><figcaption>Figure 10: Trajectory of the counting process and evolution of the intensity function for the constant excitation function with a constant baseline intensity (Case 3).</figcaption></figure>

**Case 4: Beta Excitation Function with Constant Baseline Intensity.** In this case, the excitation function corresponds to a Beta distribution defined on the interval [0, 1] with two shape parameters labeled $\alpha$ and $\beta$ such that the probability density function of the Beta distribution is given by:

$$ f(x; \alpha, \beta) = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha, \beta)}, $$
where $B(\alpha, \beta)$ is the Beta function, defined as:

$$B(\alpha, \beta) = \int_{0}^{1} x^{\alpha-1}(1-x)^{\beta-1}dx.$$

The intensity function $\lambda(t)$ in this case is given by:

$$\lambda(t) = \mu + \sum_{i=1}^{n} \alpha_i \left( \frac{(t - t_i)}{(t - t_i + \delta)} \right)^{\gamma - 1} \left( 1 - \frac{(t - t_i)}{(t - t_i + \delta)} \right)^{\beta - 1},$$

where $\mu$ is the baseline intensity, $\alpha_i$ is the excitation parameter for the i-th process, and the term involving $(t - t_i)$ represents the Beta-shaped excitation function. Figure 11 shows a realization of the trajectory of the multivariate Hawkes Process in this case.

<figure><img src="image_15.png" /><figcaption>Figure 11: Trajectory of the counting process and evolution of the intensity function for the Beta excitation function with a constant baseline intensity (Case 4).</figcaption></figure>

Second, we simulated 4-dimensional realizations of the Hawkes process for varying interaction matrices between four entities in order to illustrate the effect of the interactions on the trajectory of the multi-dimensional counting process. We considered an intensity function with the form $\lambda(t) = \mu + W \text{diag}(\phi)$ with exponential $(\phi_i^i)_{i=1,...,d}$. Figure 12 shows different behaviors obtained for various W matrices and various ordering in parameter values. We clearly observe contrasted curve shapes consistent with input parameters.

# 7 Conclusion

In this work, we proposed and characterized a self-exciting point process taking the form of a Hawkes process defined in a spatial context. The space may correspond to distinct geographic regions or any other entities, whose connections drive how much the process that is going on in a given entity generates offspring points in any other entity. Hence, this process may be applied to represent self- and inter-exciting dynamics at the nodes of any weighted-oriented network like those referred in [7, 33].

In this article, we detailed the precise form of the infinitely divisible property associated with such a Hawkes process. We computed its multi-dimensional and multi-temporal characteristic function, which allows a thorough comprehension of the process dynamics across both time and space. The generalized Laplace Transform of the process was also described as well as its first two moment functions (mean and covariance structures). We were able to obtain closed-form formula for the characteristic and moment functions, and we proposed a numerical scheme for solving the equation that we obtained.
$$W = \begin{pmatrix} a & a & 0 & 0 \\ a & a & 0 & 0 \\ 0 & 0 & b & 0 \\ 0 & 0 & 0 & b \end{pmatrix}$$

$$\mu = (c, c, c, c)$$
$$b < a$$

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

$$W = \begin{pmatrix} a & a & a & a \\ 0 & b & 0 & 0 \\ 0 & 0 & b & 0 \\ 0 & 0 & 0 & b \end{pmatrix}$$

$$\mu = (c, C, C, C)$$
$$b < a; c < C$$

$$W = \begin{pmatrix} a & a & a & a \\ 0 & b & 0 & 0 \\ 0 & 0 & b & 0 \\ 0 & 0 & b & b \end{pmatrix}$$

$$\mu = (c, C, C, C)$$
$$b < a; c < C$$

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

$$W = \begin{pmatrix} a & a & a & a \\ 0 & b & b & 0 \\ 0 & b & b & 0 \\ 0 & 0 & 0 & b \end{pmatrix}$$

$$\mu = (c, c, c, c)$$
$$b < a$$

<caption>Figure 12: Trajectory of the counting process for different interaction matrices and parameter ordering.</caption>

As pointed out in the introduction, the main results in our article are the formulas concerning the covariance structure of the process and more generally the formulas that concern multiple times for multivariate, non-stationary Hawkes processes. These formulas and all the accompany-
ing material provided in this article will allow, in further studies, the development of estimators
for the parameters of our process and the description of its behavior.

In terms of estimation, our spatial Hawkes process could be fitted to data collected during epidemics such as those caused by the phytopathogenic bacterium *Xylella fastidiosa*. The spatial Hawkes process that we proposed could allow us to draw inferences not simply based on a deterministic propagation model with a unique introduction of the disease (such assumptions were made in [1, 2]), but based on a stochastic model with multiple introductions driven by the baseline intensity function $\lambda^0$, as suggested by [3]. The most trivial ways to make estimation for our spatial Hawkes process may consist in using either the moment functions at multiple times or the probability distribution of counts in the framework of the minimum contrast method [8, 36].

In terms of description of the model behavior, several objectives may be considered. Suppose that the model is used to describe epidemics in multiple regions, then properties linking the baseline intensity $\lambda^0$ with the overall epidemic size (typically the cumulative number of points at a given time) may be derived to select the $d' < d$ regions in which the baseline intensity should be reduced to minimize the overall epidemic size, given $d'$ and the reduction factor. Similarly, properties linking inter-region migration terms in $\phi$ with the overall epidemic size may be derived to select the $d' < d(d-1)$ migration terms should be reduced to minimize the overall epidemic size, given $d'$ and the reduction factor. Beyond the control of the epidemic size, another interesting issue concerns the surveillance of the epidemics, with two typical objectives: early detection [29] and disease delimitation [6]. For example, further studies may focus on deriving the properties of the model allowing us to determine how the sampling strategies for early detection or disease delimitation should be adapted to the baseline intensity function $\lambda^0$ or the excitation function $\phi$ (including the migration terms).

Finally, it would be interesting to explore how the approach and tools proposed here could be mobilized to further characterize multidimensional self-exciting processes with dependencies (MSPD), recently introduced by [19]. MSPDs generalize Hawkes processes by allowing the excitation function to depend on a mark associated with each point event. In an epidemiological context, such dependencies could reflect the influence of pathogen variants with differing transmission profiles or of infected hosts exhibiting heterogeneous contact patterns—including the effect of super-spreaders.

**Acknowledgements.** This work was supported by an ‘INRAE MathNum – Région PACA’ PhD grant, the BEYOND Project funded by ANR (grant ANR-20-PCPA-0002) and the BeXyl Project funded by the HORIZON.2.6 programme (grant 101060593).

## References

[1] Candy Abboud, Olivier Bonnefon, Eric Parent, and Samuel Soubeyrand, *Dating and local-izing an invasion from post-introduction data and a coupled reaction-diffusion-absorption model*, Journal of Mathematical Biology **79** (2019), 765–789.

[2] Candy Abboud, Eric Parent, Olivier Bonnefon, and Samuel Soubeyrand, *Forecasting pathogen dynamics with Bayesian model-averaging: Application to Xylella fastidiosa*, *Bulletin of Mathematical Biology* **85** (2023), 67.

[3] Candy Abboud, Rachid Senoussi, and Samuel Soubeyrand, *Piecewise-deterministic Markov processes for spatio-temporal population dynamics*, Statistical Inference for Piecewise-deterministic Markov Processes, 2018, pp. 209–255.
[4] Emmanuel Bacry, Iacopo Mastromatteo, and Jean-François Muzy, *Hawkes processes in finance*, Market Microstructure and Liquidity **1** (2015), 1550005.

[5] Loubna Ben Allal, Antoine Lejay, and Radu Stefan Stoica, *Hawkes point processes based inference applied to seismic data analysis*, 2020 RING Meeting, ASGA, hal-02928408, 2020.

[6] M Charras-Garrido, Lamiae Azizi, Florence Forbes, Senan Doyle, Nathalie Peyrard, and David Abrial, *On the difficulty to delimit disease risk hot spots*, International Journal of Applied Earth Observation and Geoinformation **22** (2013), 99-105.

[7] Maria Choufany, Davide Martinetti, Rachid Senoussi, Cindy E Morris, and Samuel Soubeyrand, *Spatiotemporal large-scale networks shaped by air mass movements*, Frontiers in Applied Mathematics and Statistics **6** (2021), 602621.

[8] D. Dacunha-Castelle and M. Duflo, *Probabilités et Statistiques: Problèmes à Temps Mobile*, vol. 2, Masson, Paris, 1982.

[9] Daryl J Daley, David Vere-Jones, et al., *An Introduction to the Theory of Point Processes, Volume I: Elementary Theory and Methods*, Springer, New York, 2003.

[10] H.A. David and H.N. Nagaraja, *Order statistics*, Wiley, Hoboken, 2004.

[11] Omar El Euch and Mathieu Rosenbaum, *The characteristic function of rough Heston models*, Mathematical Finance **29** (2019), 3-38.

[12] Eymen Errais, Kay Giesecke, and Lisa R Goldberg, *Affine point processes and portfolio credit risk*, SIAM Journal on Financial Mathematics **1** (2010), 642-665.

[13] Xuefeng Gao, Xiang Zhou, and Lingjiong Zhu, *Transform analysis for Hawkes processes with applications in dark pool trading*, Quantitative Finance **18** (2018), 265-282.

[14] Xuefeng Gao and Lingjiong Zhu, *Functional central limit theorems for stationary Hawkes processes and application to infinite-server queues*, Queueing Systems **90** (2018), 161-206.

[15] László Gerencsér, György Perczel, Loránd Eross, Dániel Fabó, and Zsuzsanna Vágó, *Hawkes processes: Some key ideas, links to neuroscience and system identification*, Communications in Information and Systems **21** (2021), 385-413.

[16] Alan G Hawkes, *Spectra of some self-exciting and mutually exciting point processes*, Biometrika **58** (1971), 83-90.

[17] ———, *Hawkes processes and their applications to finance: A review*, Quantitative Finance **18** (2018), 193-198.

[18] Alan G Hawkes and David Oakes, *A cluster process representation of a self-exciting process*, Journal of Applied Probability **11** (1974), 493-503.

[19] Caroline Hillairet, Thomas Peyrat, and Anthony Réveillac, *Multivariate self-exciting processes with dependencies*, arXiv preprint arXiv:2503.15958 (2025).

[20] Caroline Hillairet and Anthony Réveillac, *Explicit correlations for the Hawkes processes*, arXiv preprint arXiv:2304.02376 (2023).

[21] Jean Jacod, *Multivariate Point Processes: Predictable projection, Radon-Nycodym derivatives, representation of martingales*, Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete **31** (1975), 235-253.
[22] Thibault Jaisson and Mathieu Rosenbaum, *Limit theorems for nearly unstable Hawkes processes: Version with technical appendix*, 2014.

[23] __________, *Limit theorems for nearly unstable Hawkes processes*, The Annals of Applied Probability **25** (2015), 600–631.

[24] Stojan Jovanović, John Hertz, and Stefan Rotter, *Cumulants of Hawkes point processes*, Physical Review E **91** (2015), 042802.

[25] Dmytro Karabash and Lingjiong Zhu, *Limit theorems for marked Hawkes processes with application to a risk model*, Stochastic Models **31** (2015), 433–451.

[26] Raviar Karim, Roger J A Laeven, and Michel Mandjes, *Exact and asymptotic analysis of general multivariate hawkes processes and induced population processes*, arXiv preprint arXiv:2106.03560 (2021).

[27] Régis C Lambert, Christine Tuleau-Malot, Thomas Bessaih, Vincent Rivoirard, Yann Bouret, Nathalie Leresche, and Patricia Reynaud-Bouret, *Reconstructing the functional connectivity of multiple spike trains using Hawkes models*, Journal of Neuroscience Methods **297** (2018), 9–21.

[28] Michel Loève, *Elementary probability theory*, Springer, 1977.

[29] Davide Martinetti and Samuel Soubeyrand, *Identifying lookouts for epidemio-surveillance: Application to the emergence of Xylella fastidiosa in France*, Phytopathology **109** (2019), 265–276.

[30] David Oakes, *The Markovian self-exciting process*, Journal of Applied Probability **12** (1975), 69–77.

[31] Nicolas Privault, *Recursive computation of the Hawkes cumulants*, Statistics & Probability Letters **177** (2021), 1091–1096.

[32] Patricia Reynaud-Bouret and Sophie Schbath, *Adaptive estimation for Hawkes processes: Application to genome analysis*, The Annals of Statistics **38** (2010), 2781–2822.

[33] Hervé Richard, Davide Martinetti, Daphné Lercier, Yann Fouillat, Buyung Hadi, Maged Elkahky, Jianqi Ding, Lucie Michel, Cindy E Morris, Karine Berthier, et al., *Computing geographical networks generated by air-mass movement*, GeoHealth **7** (2023), e2023GH000885.

[34] Marian-Andrei Rizoiu, Swapnil Mishra, Quyu Kong, Mark Carman, and Lexing Xie, *SIR-Hawkes: linking epidemic models and Hawkes processes to model diffusions in finite populations*, Proceedings of the 2018 World Wide Web Conference, 2018, pp. 419–428.

[35] Youngsoo Seol, *Limit theorems for an inverse Markovian Hawkes process*, Statistics & Probability Letters **155** (2019), 1085–1089.

[36] Samuel Soubeyrand, Florence Carpentier, Nicolas Desassis, and Joël Chadœuf, *Inference with a contrast-based posterior distribution and application in spatial statistics*, Statistical Methodology **6** (2009), no. 5, 466–477.

[37] Anatoliy Swishchuk, Rudi Zagst, and Gabriela Zeller, *Hawkes processes in insurance: Risk model, application to empirical data and optimal investment*, Insurance: Mathematics and Economics **101** (2021), 107–124.