# Efficient Simulation of Hawkes Processes using their Affine Volterra Structure

Eduardo Abi Jaber, Elie Attal, Dimitri Sotnikov

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

## To cite this version:

Eduardo Abi Jaber, Elie Attal, Dimitri Sotnikov. Efficient Simulation of Hawkes Processes using their Affine Volterra Structure. 2025. hal-05369490

HAL Id: hal-05369490

https://hal.science/hal-05369490v1

Preprint submitted on 17 Nov 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient Simulation of Hawkes Processes using their Affine Volterra Structure

Eduardo Abi Jaber*¹, Elie Attal¹, and Dimitri Sotnikov†¹, ²

¹Ecole Polytechnique, CMAP

²Engie Global Markets

November 17, 2025

## Abstract

We introduce a novel and efficient simulation scheme for Hawkes processes on a fixed time grid, leveraging their affine Volterra structure. The key idea is to first simulate the integrated intensity and the counting process using Inverse Gaussian and Poisson distributions, from which the jump times can then be easily recovered. Unlike conventional exact algorithms based on sampling jump times first, which have random computational complexity and can be prohibitive in the presence of high activity or singular kernels, our scheme has deterministic complexity which enables efficient large-scale Monte Carlo simulations and facilitates vectorization. Our method applies to any nonnegative, locally integrable kernel, including singular and non-monotone ones. By reformulating the scheme as a stochastic Volterra equation with a measure-valued kernel, we establish weak convergence to the target Hawkes process in the Skorokhod $J_1$-topology. Numerical experiments confirm substantial computational gains while preserving high accuracy across a wide range of kernels, with remarkably improved performance for a variant of our scheme based on the resolvent of the kernel.

## 1 Introduction

Hawkes processes, originally introduced by Hawkes (1971), form a class of self-exciting counting processes where each event increases the likelihood of future events. This self-excitation is governed by a convolution kernel $K \in L^1([0, T], \mathbb{R}_+)$ . More precisely, the counting process $N = (N_t)_{t \le T}$ is characterized by its stochastic instantaneous intensity

$$ \lambda_t = g_0(t) + \int_0^t K(t-s) dN_s, \quad t \le T, \qquad (1.1) $$

where $g_0 \in L^1([0, T], \mathbb{R}_+)$ is the exogenous intensity, driving events that occur independently of the past trajectory. The second term, $K * dN$, is the endogenous component, which increases the intensity at each jump of $N$ and retains memory of these jumps through the kernel $K$. We refer to Laub, Lee, Pollett, and Taimre (2025) for a general review on Hawkes processes.

When $K(t) = c e^{-bt}$, with $c > 0, b \in \mathbb{R}$, the process $(N, \lambda)$ is Markovian. However, for general kernels $K$, the Volterra-type memory breaks the Markov property. This ability to encode self-excitation with memory makes Hawkes processes popular in diverse domains, including financial mathematics (Bacry, Mastromatteo, and Muzy, 2015; Hawkes, 2018; Embrechts, Liniger, and Lin, 2011; Rosenbaum and Tomas, 2021), insurance (Deaconu, Lejay, Lesage, Meira, Nichil, and State, 2022; Bessy-Roland, Boumezoued, and Hillairet, 2021; Baldwin, Gheyas, Ioannidis, Pym, and Williams, 2017; Hillairet, Réveillac, and Rosenbaum, 2023), seismology (Ogata, 1988; Baeumer, Davis, and Wang, 2024) and modeling of biological neurons (Hodara and Löcherbach, 2017; Goncalves and Gresland, 2022).

In the non-Markovian setting, computing key quantities, such as the distribution of $N_t$, its moments, or its quantiles, is generally intractable in closed form. Consequently, simulation for Monte Carlo methods plays a central role in applications.

<footnote>*EAJ is grateful for the financial support from the Chaires FiME-FDD and Financial Risks at Ecole Polytechnique.</footnote>

<footnote>†DS is grateful for the financial support provided by Engie Global Markets.</footnote>
One of the most widely used simulation methods is Ogata's thinning algorithm from Ogata (1981), based on the acceptance-rejection technique, which generates candidate jumps and discards some of them at random. It requires explicit upper bounds on the intensity $\lambda_t$, which are typically available when $K$ is bounded. However, the scheme is no longer exact when $K$ is unbounded (e.g., $K(0) = +\infty$), and the case of a non-monotone kernel requires special treatment.

Another possibility for simulation is to use the population representation (also called immigrant-birth or cluster representation, see Hawkes and Oakes (1974); Møller and Rasmussen (2005)), in which the Hawkes process is seen as a branching structure: each jump gives birth to an inhomogeneous Poisson process, which can be simulated directly. This approach only requires $K$ to be integrable.

Both simulation methods share the same important drawback: their computation cost is random, since jump times are simulated first, before building the processes $\lambda$ and $N$. Therefore, their complexity has order of magnitude $O(N_T^2)$, where the square comes from the non-Markovian Volterra structure. In regimes with strong self-excitation, $N_T$ can be very large, making even a single simulation costly. Furthermore, this makes vectorized sampling more difficult in the context of Monte Carlo methods. In addition, in many applications, the quantities of interest depend on the finite-dimensional distributions of $N$ (see, for instance, Deaconu, Lejay, Lesage, Meira, Nichil, and State (2022); Bessy-Roland, Boumezoued, and Hillairet (2021)), rather than on the event times themselves. In such cases, simulating the individual jump times, with their inherently random computational complexity, becomes undesirable and may substantially degrade performance. Even if one restricts to the Markovian case, where $K$ is an exponential kernel, an exact simulation method exists that achieves linear, but still stochastic, complexity of order $O(N_T)$ (see Dassios and Zhao (2013)).

This naturally leads to a central question:

*Can Hawkes processes be simulated efficiently with a deterministic complexity?*

To address this question, we make two key observations:

1. **Only integrated quantities matter:** Once the integrated intensity $\Lambda_{t_i, t_{i+1}} := \int_{t_i}^{t_{i+1}} \lambda_s ds$ has been simulated over a deterministic interval $[t_i, t_{i+1})$, the increment of the Hawkes process $N_{t_{i+1}} - N_{t_i}$, conditional on $\Lambda_{t_i, t_{i+1}}$, can be sampled using a Poisson distribution with parameter $\Lambda_{t_i, t_{i+1}}$.

2. **Hawkes processes are Affine Volterra processes with jumps:** as shown in Abi Jaber (2021, Example 2.3), Bondi, Livieri, and Pulido (2024); Cuchiero and Teichmann (2019). More precisely, based on the first observation, we focus on the dynamics of the (less conventional) integrated intensity $\Lambda$, rather than the instantaneous intensity $\lambda$. Integrating (1.1) and applying Fubini's theorem yields

$$ \Lambda_t := \int_0^t \lambda_s ds = \int_0^t g_0(s) ds + \int_0^t K(t-s)(\Lambda_s + Z_s) ds, \quad t \le T, \quad (1.2) $$

where $Z := N - \Lambda$ is a martingale with predictable quadratic variation equal to $\Lambda$, which, in other words, is affine in $\Lambda$. This affine Volterra structure translates into the knowledge of the characteristic function of the integrated intensity $\Lambda$ in semi-explicit form in terms of a non-linear deterministic Volterra equation.

Recently, in the context of continuous affine Volterra processes (see, e.g., Abi Jaber, Larsson, and Pulido (2019)), the affine structure of integrated quantities has been leveraged to develop an efficient simulation method, namely the iVi scheme of Abi Jaber (2025); Abi Jaber and Attal (2025), which relies on Inverse Gaussian increments.

Building on these ideas, we introduce a new paradigm for approximating $(\Lambda, N)$ directly on a fixed time grid $(t_i^n)_{0 \le i \le n} := (i \frac{T}{n})_{0 \le i \le n}$ with $n + 1$ points. Our method replaces the random $O(N_T^2)$ complexity of classical schemes with a deterministic $O(n^2)$, giving explicit control over the trade-off between computational cost and accuracy in Monte Carlo experiments. The scheme simulates the increments $\hat{\Lambda}_{i, i+1}^n \approx \Lambda_{t_{i+1}}^n - \Lambda_{t_i}^n$ through Inverse Gaussian sampling, where the parameters are derived from an approximation on the characteristic function of $\Lambda$, followed by Poisson sampling for $\hat{N}_{i, i+1}^n \approx N_{t_{i+1}}^n - N_{t_i}^n$. The approximation of the integrated intensity is then updated to reflect the dynamic (1.2). As an optional step, the jump times $T^n = (\tau_k^n)_{k=1, \dots, N_T^n}$ of the counting process $N$, where $N_T^n = \sum_{i=0}^{n-1} \hat{N}_{i, i+1}^n$ denotes the total number of jumps on $[0, T]$, can be approximated from $\hat{N}^n$. On each interval $[t_i, t_{i+1})$, conditional on the number of jumps $\hat{N}_{i, i+1}^n$ in that interval, we sample $\hat{N}_{i, i+1}^n$ jump times uniformly within the interval and then sort them in increasing order. The scheme
is detailed in Algorithm 1. We refer to our method as Hawkes iVi, for integrated Volterra implicit following the terminology of Abi Jaber (2025); Abi Jaber and Attal (2025).

<caption>**Algorithm 1 - The Hawkes iVi scheme:** Simulation of $\hat{\Lambda}, \hat{N}$</caption>

1: **Input:** $g_0, K \in L^1([0, T], \mathbb{R}_+)$, uniform partition of $[0, T]$, $t_i^n = iT/n$, with $i = 0, \dots, n$.

2: **Output:** $\hat{\Lambda}_{i,i+1}^n, \hat{N}_{i,i+1}^n$ for $i = 0, \dots, n-1$ and (optional) $\mathcal{T}^n = (\tau_k^n)_{k=1,\dots,N_T^n}$.

3: Compute

$$k_j^n = \int_{t_j^n}^{t_{j+1}^n} K(s) ds, \quad j = 0, \dots, n-1. \qquad (1.3)$$

4: **for** $i = 0$ to $n-1$ **do**

5: Compute the quantity:

$$\alpha_i^n = \int_{t_i^n}^{t_{i+1}^n} g_0(s) ds + \sum_{j=0}^{i-1} k_{i-j}^n \hat{N}_{j,j+1}^n, \qquad (1.4)$$

with $\alpha_0^n = \int_0^{t_1^n} g_0(s) ds$, for $i = 0$.

6: Simulate the increment:

$$\hat{N}_{i,i+1}^n \sim \mathcal{P}(\xi_i^n), \quad \text{where} \quad \xi_i^n \sim IG\left(\frac{\alpha_i^n}{1-k_0^n}, \left(\frac{\alpha_i^n}{k_0^n}\right)^2\right). \qquad (1.5)$$

7: Compute the increment:

$$\hat{\Lambda}_{i,i+1}^n = \alpha_i^n + k_0^n \hat{N}_{i,i+1}^n. \qquad (1.6)$$

8: Optional step (jump times): Simulate $\hat{N}_{i,i+1}^n$ i.i.d. uniform random variables on $[t_i, t_{i+1})$, sort, and append them to $\mathcal{T}^n$.

9: **end for**

$IG(\cdot, \cdot)$ and $\mathcal{P}(\cdot)$ denote respectively the Inverse Gaussian (see Appendix B) and the Poisson distributions.

The Inverse Gaussian sampling in (1.5) is justified by an approximation of the characteristic functional of $\Lambda_{t_{i+1}^n}^n - \Lambda_{t_i^n}^n$. Conditional on the latter, $N_{t_{i+1}^n} - N_{t_i^n}$ follows a Poisson distribution with parameter $\Lambda_{t_{i+1}^n}^n - \Lambda_{t_i^n}^n$. Since the scheme outputs approximations of the increments of $\Lambda$ and $N$ on the time grid, the piecewise constant càdlàg processes we are interested in are given by

$$\Lambda_t^n := \sum_{i=0}^{\lfloor nt/T \rfloor - 1} \hat{\Lambda}_{i,i+1}^n, \quad t \le T, \quad n \ge 1, \qquad (1.7)$$

$$N_t^n := \sum_{i=0}^{\lfloor nt/T \rfloor - 1} \hat{N}_{i,i+1}^n, \quad t \le T, \quad n \ge 1. \qquad (1.8)$$

Our main theoretical contribution, stated in Theorem 3.1, establishes the weak convergence of the pair of processes $(\Lambda^n, N^n)$, defined in (1.7)-(1.8), to the limiting pair $(\Lambda, N)$ in the Skorokhod $J_1$-topology. The general strategy follows the classical approach for weak convergence: proving tightness and identifying the limit. However, each of these steps is far from straightforward in our setting and requires a novel analysis. The key idea lies in recasting the discrete scheme as a stochastic Volterra equation with a measure-valued kernel (Lemma 5.3). This reformulation is crucial, as it reveals that the coupling between $N^n$ and $\Lambda^n$ naturally inherits a Volterra-type structure, and it enables refined moment estimates, leading to tightness of the approximating sequence and a clear identification of the limiting dynamics through stability results for stochastic Volterra equations (Lemma 5.10). Furthermore, a delicate point in our proof is to show that the limit process $N$ retains the point-process structure with unit jumps, a nontrivial property that is not automatic under weak convergence in $J_1$.
Beyond eliminating random complexity, our approach offers several other benefits:

* **Vectorization-friendly:** choosing $n$ before simulation enables efficient, parallel Monte Carlo sampling.

* **Robust kernel framework:** the method uses integrated kernel quantities only, accommodating any $K \in L^1([0, T], \mathbb{R}_+)$, including non-monotone and singular kernels.

* **Weak convergence guarantee:** we prove that the couple of processes $(\Lambda^n, N^n)$ defined in (1.7) and (1.8) converges weakly to $(\Lambda, N)$ in the Skorokhod topology, by reformulating the scheme as a stochastic Volterra equation with a measure-valued kernel (see Section 3).

* **Possible simulation of jump times:** although Algorithm 1 simulates $(\Lambda, N)$ directly, which is more convenient for Monte Carlo estimation, jump times can still be approximated satisfactorily (see Subsection 4.2).

* **Linear complexity in the Markovian case:** for exponential kernels $K(t) = c e^{-bt}$ with $c > 0, b \in \mathbb{R}$, a linear $O(n)$ variant exists. To see this, we remark that, in this case, (1.4) can be rewritten as

$$ \alpha_i^n = \int_{t_i^n}^{t_{i+1}^n} g_0(s) ds + e^{-bT/n} \left( \alpha_{i-1}^n - \int_{t_{i-1}^n}^{t_i^n} g_0(s) ds \right) + k_1^n \hat{N}_{i-1,i}^n, \quad 1 \le i \le n-1, \quad (1.9) $$

which can be computed in $O(1)$ time. This particular case is developed in Subsection 4.4.

* **Reformulation using resolvent kernels:** replacing $K$ by its resolvent of the second kind drastically improves performance empirically, although positivity of the scheme is no longer guaranteed. This alternative scheme is tested in Subsection 4.3.

All in all, our numerical experiments in Section 4 confirm substantial computational gains while preserving high accuracy across a wide range of kernels, with remarkably improved performance for the variant of our scheme based on the resolvent of the kernel.

**Connection to the literature.** This work is inspired by the iVi scheme introduced by Abi Jaber (2025); Abi Jaber and Attal (2025) for continuous affine Volterra processes, and further enhanced by Grzelak and Zaugg (2025) for the lifted Heston model. Extending the approach to point processes requires several important adjustments. First, although both schemes rely on an Inverse Gaussian approximation, the sampling is performed differently in our jump setting to preserve the Volterra structure of the process $\Lambda$ while ensuring that $N$ remains integer-valued, a crucial distinction compared to the continuous case. Second, the proof of convergence necessitates new estimates, particularly to establish that the limiting process is a point process with unit jumps. From a numerical perspective, the iVi Hawkes scheme differs substantially from the simulation methods commonly used in the literature due to its deterministic complexity, and therefore requires a dedicated comparison to assess its computational advantages. Overall, while both iVi schemes share the same fundamental structure, each introduces improvements in its respective domain, demonstrating the effectiveness of the method for approximating affine Volterra processes in different contexts.

In the literature, weak convergence results towards Hawkes-type limits are rather scarce. To the best of our knowledge, Kirchner (2016, 2017) is one of the very few works studying such a convergence result using discrete-time INAR($p$) and INAR($\infty$) processes. One key difference with our approach is that we also establish the convergence of the intensity process $\Lambda^n$, whereas Kirchner (2016) focuses only on the point process itself. A numerical comparison of Algorithm 1 with a refined version of the approximation method of Kirchner (2016, 2017) is provided in Subsection 4.2. More general weak convergence and stability results for stochastic Volterra equations with jumps, such as those in Abi Jaber (2021); Abi Jaber, Cuchiero, Larsson, and Pulido (2021), cannot be directly applied in our context. These works rely on kernels $t \mapsto K(t)$ belonging to some $L^p$ space, while our setting is inherently discrete. After reformulation, our scheme can indeed be seen as a stochastic Volterra equation, but with a measure-valued kernel $K(ds)$, which requires a distinct treatment.

The most renowned simulation schemes for the Hawkes process, to which we compare our method in Section 4, are Ogata's algorithm (Ogata, 1981), the population approach (Møller and Rasmussen, 2005) and the exact exponential scheme for the Markovian case (Dassios and Zhao, 2013). Other methodologies can be found in the literature, often adapted to a specific class of kernels. Among these, Chen, Hawkes, and Scamas (2021) propose to adapt Ogata's algorithm by adding a small time-shift to handle singular Mittag-Leffler kernels, while Duarte, Löcherbach, and Ost (2019) introduce a thinning algorithm for sums of Erlang kernels.
**Outline of the paper.** In Section 2, we explain how the scheme can be derived from a joint approximation of the dynamics and of the characteristic functional of the integrated intensity, while Section 3 focuses on its weak convergence. Section 4 contains all numerical tests, and Section 5 is dedicated to proving the weak convergence result of Theorem 3.1. Finally, Appendix A defines the resolvent of kernel and gives their explicit form for some usual kernels, and Appendix B contains properties of the Inverse Gaussian distribution.

**Notation.** Throughout the paper, we use the symbol $*$ to denote convolution:

$$ (R * K)(t) := \int_{0}^{t} R(t-s)K(s) ds. $$

## 2 Deriving the Hawkes iVi scheme

The scheme is based on a joint approximation on the dynamics (1.2) and on the characteristic functional of the integrated intensity $\Lambda_t := \int_0^t \lambda_s ds$. We insist on the fact that the present paper deals with Hawkes process on $[0, T]$, as opposed to the whole nonnegative real line. Thus, we only require that $g_0, K \in L^1([0, T], \mathbb{R}_+)$.
The processes $(\Lambda_t, N_t)_{t \le T}$ can be interpreted as the restriction to $[0, T]$ of the Hawkes process on $\mathbb{R}_+$, with compactly supported exogenous intensity and memory kernel respectively given by

$$ \tilde{g}_0(t) := g_0(t) 1_{t \le T}, \quad \text{and} \quad \tilde{K}(t) := K(t) 1_{t \le T}, \quad t \ge 0. $$

Note, however, that Algorithm 1 is well-defined as long as $k_0^n := \int_0^{T/n} K(s) ds < 1$, due to (1.5), which is always true for $n$ large enough.

### 2.1 Approximating the dynamics

We recall that integrating (1.1) and applying Fubini's interchange leads to the dynamics of the integrated intensity

$$ \Lambda_t = \int_0^t g_0(s) ds + \int_0^t K(t-s) N_s ds, \quad t \le T. $$

Thus, similarly to Abi Jaber and Attal (2025, Proposition 1.1), we can decompose the increments $\Lambda_t - \Lambda_s$ into an $\mathcal{F}_s$-measurable part, corresponding to the non-Markovianity of the dynamics, and an increment between $s$ and $t$, as follows

$$ \Lambda_t - \Lambda_s = G_s(t) + \int_s^t K(t-r)(N_t - N_r) dr, \quad s \le t \le T, $$

where $G_s(t)$ is $\mathcal{F}_s$-measurable and given by

$$ G_s(t) := \int_s^t g_0(r) dr + \int_0^s \int_s^t K(u-r) du dN_r, \quad s \le t \le T. \quad (2.1) $$

We now want to approximate $\Lambda_{t^n_{i+1}} - \Lambda_{t^n_i}$. To do so, we derive a discretized proxy of $G_{t^n_i}(t^n_{i+1})$ using left-endpoint approximations of the integrals:

$$ \begin{align*} G_{t^n_i}(t^n_{i+1}) &= \int_{t^n_i}^{t^n_{i+1}} g_0(s) ds + \sum_{j=0}^{i-1} \int_{t^n_j}^{t^n_{j+1}} \int_{t^n_i}^{t^n_{i+1}} K(u-r) du dN_r \\ &\approx \int_{t^n_i}^{t^n_{i+1}} g_0(s) ds + \sum_{j=0}^{i-1} \int_{t^n_j}^{t^n_{j+1}} \int_{t^n_i}^{t^n_{i+1}} K(u-t^n_j) du dN_r \\ &\approx \int_{t^n_i}^{t^n_{i+1}} g_0(s) ds + \sum_{j=0}^{i-1} k^n_{i-j} \hat{N}^n_{j,j+1} =: \alpha^n_i, \quad i \le n-1, \quad n \ge 1, \end{align*} $$

where the $(k^n_i)_{i \le n-1}$ are defined in (1.3). Hence, we have done the approximation $G_{t^n_i}(t^n_{i+1}) \approx \alpha^n_i$, which explains the definition of $(\alpha^n_i)_{i \le n-1}$ in (1.4). Finally, to obtain (1.6), we also make a left-endpoint approximation in the second term of our decomposition of the increments of $\Lambda$, leading to

$$ \Lambda_{t^n_{i+1}} - \Lambda_{t^n_i} = G_{t^n_i}(t^n_{i+1}) + \int_{t^n_i}^{t^n_{i+1}} K(t^n_{i+1} - r) (N_{t^n_{i+1}} - N_r) dr \approx \alpha^n_i + k_0^n \hat{N}^n_{i,i+1} = \hat{\Lambda}^n_{i,i+1}. \quad (2.2) $$
We will see in Section 3 that this discretization amounts to replacing the memory kernel K by a measure-valued one, built on integrated quantities, hence preserving the Volterra structure of the process.

## 2.2 Deriving the Inverse Gaussian parameters

The Inverse Gaussian distribution comes from a short-time approximation of the characteristic functional of $\Lambda$. We remark that the exponent in the characteristic function of Inverse Gaussian random variables (see B.3) has the structure of a root of a second-order polynomial. It is a well-known fact that characteristic functionals of affine Volterra processes are given semi-explicitly, up to the solution to a Riccati–Volterra equation (see Abi Jaber (2021, Theorem 2.2) and Abi Jaber, Larsson, and Pulido (2019)). Consequently, using second-order Taylor expansions in an implicit Euler fashion, we can expect to obtain second-order polynomials as small-time approximations of such equations. We detail this methodology in the present subsection to identify the right Inverse Gaussian parameters for our scheme.

From Abi Jaber (2021, Example 2.3)<fnref n="1" />, we know that for $w \in \mathbb{C}$ with $\Re(w) \le 0$,

$$ \mathbb{E} \left[ \exp \left( w \left( \Lambda_{t_{i+1}^n} - \Lambda_{t_i^n} \right) \right) \mid \mathcal{F}_{t_i^n} \right] = \exp \left( \int_{t_i^n}^{t_{i+1}^n} F(\Psi(t_{i+1}^n - s)) dG_{t_i^n}(s) \right), \quad i \le n-1, \quad n \ge 1, $$

where $G_{t_i^n}$ is defined in (2.1). Here, the function F is given by

$$ F(u) = w + e^u - 1, \quad u \in \mathbb{C}, \tag{2.3} $$

and $\Psi$ is the continuous solution, with nonpositive real part, to the Volterra equation

$$ \Psi(t) = \int_{0}^{t} K(t-s) F(\Psi(s)) ds, \quad t \le T. $$

Using the notation $\Delta t = \frac{T}{n}$, we can make the left-endpoint approximations

$$ \int_{t_i^n}^{t_{i+1}^n} F(\Psi(t_{i+1}^n - s)) dG_{t_i^n}(s) \approx F(\Psi(\Delta t)) G_{t_i^n}(t_{i+1}^n) \approx F(\Psi(\Delta t)) \alpha_i^n, \tag{2.4} $$

along with

$$ \Psi(\Delta t) \approx \int_{0}^{\Delta t} K(s) ds F(\Psi(\Delta t)) = k_{0}^{n} F(\Psi(\Delta t)), \tag{2.5} $$

for small $\Delta t$. We remark that $\Psi(\Delta t)$ is here approximated by the solution to the equation $x = k_0^n F(x)$, which is not quadratic yet. Using the fact that $\Psi(0) = 0$ so that $\Psi(\Delta t)$ is small, we can approximate F by its second-order Taylor expansion, given by

$$ F(u) = w + e^u - 1 \approx w + u + \frac{u^2}{2}, \tag{2.6} $$

for small $u$. Combining all of our approximations, we obtain

$$ \mathbb{E} \left[ \exp \left( w \left( \Lambda_{t_{i+1}^n} - \Lambda_{t_i^n} \right) \right) \mid \mathcal{F}_{t_i^n} \right] \approx \exp \left( \hat{\Psi}(\Delta t) \frac{\alpha_i^n}{k_0^n} \right), $$

where $\hat{\Psi}(\Delta t)$ is the zero root with nonpositive real part of the second order polynomial equation

$$ \frac{k_0^n}{2} \left( \hat{\Psi}(\Delta t) \right)^2 - (1 - k_0^n) \hat{\Psi}(\Delta t) + w k_0^n = 0. $$

After computation, we obtain

$$ \hat{\Psi}(\Delta t) = \frac{1 - k_0^n}{k_0^n} \left( 1 - \sqrt{1 - 2 \left( \frac{k_0^n}{1 - k_0^n} \right)^2 w} \right), $$

<footnote n="1">¹The result is proved in a more general framework which requires assumptions of the memory kernel K. However, in the special case of Hawkes processes, it can be shown to hold for any $K \in L^1([0, T], \mathbb{R}_+)$ using similar arguments.</footnote>
which corresponds to the characteristic functional of an Inverse Gaussian random variable (see Appendix B) $\xi_i^n$ with parameters $(\frac{\alpha_i^n}{1-k_0^n}, (\frac{\alpha_i^n}{k_0^n})^2)$. Finally, since $N$ is a counting process with integrated intensity $\Lambda$, conditional on $\Lambda_{t_{i+1}}^n - \Lambda_t^n$, we have

$$N_{t_{i+1}}^n - N_{t_i^n} \sim \mathcal{P}(\Lambda_{t_{i+1}}^n - \Lambda_{t_i^n}) \approx \mathcal{P}(\xi_i^n),$$

where $\mathcal{P}(\lambda)$ denotes the Poisson distribution with parameter $\lambda$.

It is worth noting that we use equation (1.6) to compute the increment $\hat{\Lambda}_{i,i+1}^n$, rather than setting $\hat{\Lambda}_{i,i+1}^n = \xi_i^n$, as might be naturally suggested by the iVi scheme for Volterra processes introduced by Abi Jaber and Attal (2025). This adjustment is due to the discrete nature of point processes: since $N_{i,i+1}^n$ must be an integer, it can no longer be directly extracted from the equation

$$\alpha_i^n + k_0^n \hat{N}_{i,i+1}^n = \hat{\Lambda}_{i,i+1}^n,$$

which arises from the discretization (2.2). Instead, after simulating $\hat{N}_{i,i+1}^n$ as a Poisson random variable with the Inverse Gaussian parameter $\xi_i^n$, $\hat{\Lambda}_{i,i+1}^n$ is determined from the equation above, leading to (1.6) and ensuring the simulation of consistent trajectories for both the jump process $N$ and its integrated intensity $\Lambda$. By consistency, we understand the Volterra-type relation between the discretized processes $\Lambda^n$ and $N^n$ (see Lemma 5.3) that is a crucial ingredient of the proof of convergence. In particular, this helps us avoid situations in which the simulated counting process $N^n$ has no jumps (all sampled Poisson random variables are null) and $\Lambda^n > \int_0^T g_0(s) ds$, which should not be possible for a Hawkes process due to (1.2).

Combining the above, we arrive at Algorithm 1.

## 3 Weak convergence of the Hawkes iVi scheme

In this section, we formulate a convergence result and sketch the main steps of its proof. The technical details are postponed to Section 5. We assume that $g_0, K \in L^1([0, T], \mathbb{R}_+)$ and we denote by $n_0$ the smallest natural $n$, such that $k_0^n < 1$, so that the scheme in Algorithm 1 is well-defined. For the sake of notational simplicity only, we will assume that $n_0 = 1$, but we stress that all the results hold for general $n_0$. In particular, our proof of convergence does not require $\int_0^T K(s) ds < 1$.

We consider the piecewise constant non-decreasing processes $\Lambda^n$ and $N^n$ defined in (1.7) and (1.8) respectively. We assume to all our random variables are defined on the same complete probability space $(\Omega, \mathcal{F}, \mathbb{P})$ and we introduce the filtration

$$\mathcal{F}_t^n := \sigma(\hat{\mathcal{F}}_t^n \cup \mathcal{F}_0), \quad t \le T, \quad n \ge 1, \qquad (3.1)$$

where

$$\hat{\mathcal{F}}_t^n := \sigma(\hat{N}_{i,i+1}^n, 0 \le i \le \lfloor nt/T \rfloor - 1), \quad t \le T, \quad n \ge 1,$$

and $\mathcal{F}_0$ is the set of all $\mathbb{P}$-null elements of $\mathcal{F}$. For each $n \ge 1$, the filtration $(\mathcal{F}_t^n)_{t \le T}$ is the one generated by the process $N^n$ and satisfies the usual conditions. Note that both processes $\Lambda^n$ and $N^n$ are adapted to the filtration.

We first mention two convenient properties of our scheme, that are pivotal in the proof of weak convergence. Firstly, for each $n \ge 1$, the process $N^n - \Lambda^n$ is shown to be a martingale in Lemma 5.2, which is enjoyable since $N - \Lambda$ is a martingale in the limit. However, it is important to note that $\Lambda^n$ is not a predictable process, since at each step $\hat{\Lambda}_{i,i+1}^n$ is set according to (1.6) after sampling $\hat{N}_{i,i+1}^n$. Second, the relation between $N^n$ and $\Lambda^n$ can also be rewritten as a stochastic Volterra equation, but with a measure-valued kernel (see Lemma 5.3). More precisely, introducing the finite non-negative measure $K^n := \sum_{i=0}^{n-1} k_i^n \delta_{t_i^n}$ on $[0, T]$, where the $(k_i^n)_{i \le n-1}$ are defined in (1.3) and $\delta_t$ corresponds to the Dirac measure in $t$, we have

$$\Lambda_t^n = \int_0^{\lfloor \frac{nt}{T} \rfloor} g_0(s) ds + \int_{[0,t]} N_{t-s}^n K^n(ds), \quad t \le T.$$
This relation explicits how our discretization procedure (see Subsection 2.1) transforms the Volterra structure of (1.2).

Combining these properties with estimates on the processes, we prove the weak convergence of the scheme, which is given by the following theorem.

**Theorem 3.1.** We have the weak convergence

$$ (\Lambda^n, N^n) \xrightarrow{n \to \infty} (\Lambda, N), $$

in the Skorokhod $J_1$-topology, where $N$ is a Hawkes process on $[0, T]$ with exogenous intensity $g_0$ and memory kernel $K$, and $\Lambda$ is its integrated intensity.

*Proof.* Lemma 5.7 gives the tightness of the couple of processes. Lemma 5.8 then characterizes any weak accumulation point $(\Lambda, N)$, showing that $N$ is a counting process (restricted to $[0, T]$) with integrated intensity $\Lambda = G_0 + K * N$. The uniqueness in law follows from the first lemma in Hawkes and Oakes (1974), with $H_0 = \emptyset$ in our case. Since all accumulation points have the same law, we obtain weak convergence of the sequence. $\square$

# 4 Numerical results

In this section, we compare the iVi Hawkes scheme with two popular approaches for simulating Hawkes processes, the population-based method (see, for instance, Møller and Rasmussen (2005)) and the thinning method introduced by Ogata (1981), both in terms of the marginal laws of the processes $N$ and $\Lambda$, and the distributions of the arrival times.

## 4.1 Population approach and Ogata's algorithm

We begin by recalling the ideas behind these two methods and highlighting their respective advantages and disadvantages.

**Population approach** The population approach (Møller and Rasmussen, 2005) models the Hawkes process as a superposition of non-homogeneous Poisson processes of two types:

1. “Migrants”, or arrival times generated by the Poisson component $g_0$ of the intensity.

2. “Descendants”, i.e., arrivals triggered by previous arrivals (of both types) corresponding to the terms $K(t - \tau)$, $t \ge \tau$ of the intensity for a given arrival time $\tau$.

This decomposition reveals the tree structure of the Hawkes process and reduces the simulation problem to simulating non-homogeneous Poisson processes applied sequentially to successive “generations” of arrival times. Simulation of a non-homogeneous Poisson process is straightforward when the inverse integrated intensity is available, thanks to the following result (cf. e.g., Brown and Nair (1988)).

**Theorem 4.1** (Random time change). Let $(\tau_1, \dots, \tau_{N_T})$ be a realization of a point process with integrated intensity $\Lambda$ over $[0, T]$ and $\Lambda(T) < \infty$ a.s. Then, the transformed arrivals $(\Lambda(\tau_1), \dots, \Lambda(\tau_{N_T}))$ form a standard Poisson process with unit rate on $[0, \Lambda(T)]$.

Hence, to simulate a non-homogeneous Poisson process with (deterministic) integrated intensity $\Lambda$, one simply simulates a standard Poisson process on $[0, \Lambda(T)]$ and then applies $\Lambda^{-1}$ to its arrival times. Algorithm 2 summarizes this procedure. We use the notation $\bar{K}(t) := \int_0^t K(s) ds$ for the integrated kernel.
Algorithm 2 Population approach

1: Input: $G_0, G_0^{-1}, \bar{K}, \bar{K}^{-1}$.

2: Output: event times $\mathcal{T} = (\tau_i)_{i=1,...,N_T}$.

3: Initialize empty sets: parent queue $\mathcal{P} \leftarrow \emptyset$, output set $\mathcal{T} \leftarrow \emptyset$.

4: Simulate “migrants” as an inhomogeneous Poisson process on $[0, T]$ with integrated intensity $G_0$, yielding arrivals $\tau_1^0, ..., \tau_{n_0}^0$, where $n_0$ denotes the total number of arrivals on $[0, T]$.

5: Add $\{\tau_1^0, ..., \tau_{n_0}^0\}$ to $\mathcal{P}$.

6: while $\mathcal{P} \neq \emptyset$ do

7: Remove the first element $\tau$ from $\mathcal{P}$.

8: Simulate arrivals $\tilde{\tau}_1^T, ..., \tilde{\tau}_{n_T}^T$ of a standard Poisson process on $[0, \bar{K}(T - \tau)]$.

9: Transform to descendant event times: $\tau_i^T \leftarrow \tau + \bar{K}^{-1}(\tilde{\tau}_i^T)$ for $i = 1, ..., n_T$.

10: Add $\{\tau_1^T, ..., \tau_{n_T}^T\}$ to $\mathcal{P}$.

11: Add $\tau$ to $\mathcal{T}$.

12: end while

13: Sort $\mathcal{T}$ in ascending order.

**Ogata's algorithm** The algorithm introduced by Ogata (1981) generalizes the standard thinning algorithm. The main idea is as follows: if the kernel $K$ is finite and decreasing, and $g_0$ is bounded, then for $t \in [0, T]$, given $\mathcal{F}_t$, one can bound the intensity of the Hawkes process between $t$ and the next arrival by

$$M = \max_{u \in [t, T]} \{g_0(u)\} + \sum_{\tau_i \le t} K(t - \tau_i),$$

and apply the standard thinning method (cf. Lewis and Shedler (1979)). Here, $(\tau_i)_{i=1,...,N_t}$ denotes the previously simulated arrivals. We summarize Ogata's method in Algorithm 3.

Algorithm 3 - Ogata's algorithm

1: Input: $g_0, K$ with $\max_{u \in [0, T]} g_0(u) < \infty$ and decreasing $K$, such that $K(0) < \infty$.

2: Output: $\mathcal{T} = (\tau_i)_{i=1,...,N_T}$.

3: Set $t \leftarrow 0$.

4: while $t < T$ do

5: Set $M \leftarrow \max_{u \in [t, T]} \{g_0(u)\} + \sum_{\tau_i \le t} K(t - \tau_i)$.

6: Simulate $\xi \sim \text{Exp}(1/M)$ and $U \sim \mathcal{U}(0, 1)$

7: Set $\tau \leftarrow t + \xi$.

8: If $\tau < T$ and $U < \frac{g_0(\tau) + \sum_{\tau_i \le t} K(\tau - \tau_i)}{M}$, set $\tau_{i+1} \leftarrow \tau$.

9: Set $t \leftarrow \tau$.

10: end while

Both methods, as presented, allow for exact simulation of the Hawkes process. However, the population approach requires knowledge of the inverse integrated kernel, which may be problematic for complex kernels like the Mittag-Leffler kernel used by Chen, Hawkes, and Scalas (2021). In addition, this kernel and the fractional kernel $K(t) = \frac{c}{\Gamma(\alpha)}t^{\alpha-1}$, $\alpha < 1$ violate the assumption $K(0) < \infty$ required by Ogata's algorithm. Moreover, the monotonicity assumption is violated by Erlang or, more generally, gamma kernels (see, e.g., Duarte et al. (2019)).

To handle singular kernels, Chen et al. (2021) propose shifting the kernel by a small $\epsilon > 0$ when computing $M$. In this case, Ogata's algorithm is no longer exact, though the error is almost negligible. Similarly, to overcome the monotonicity requirement, one may use any decreasing kernel $K^m \ge K$ with $K^m(0) < \infty$ to compute $M$. However, this increases the number of simulations in the acceptance-rejection procedure and thus reduces the efficiency of the algorithm.
## 4.2 Numerical experiments

For the numerical illustration, we will use the gamma kernel

$$K_{\Gamma}(t) = ce^{-bt} \frac{t^{\alpha-1}}{\Gamma(\alpha)}, \quad b > 0, \alpha > 0, c > 0, t \ge 0,$$

with the parameters $b = 3$, $\alpha = 2$, $c = 0.9b^{\alpha}$ and $\mu = 5$, $T = 1$. Note that when $\alpha$ is an integer, the gamma kernel reduces to the Erlang kernel widely used in the applications of the Hawkes processes.

Figure 1 shows the sample trajectories of $N$, $\Lambda$ and the instantaneous intensity $\lambda$ generated with the iVi scheme with time step 0.001. The instantaneous intensity is computed using the left-endpoint approximation

$$\lambda_{t_i} = g_0(t_i) + \int_0^{t_i} K(t_i - s) dN_s = g_0(t_i) + \sum_{j=0}^{i-1} \int_{t_j}^{t_{j+1}} K(t_i - s) dN_s \approx g_0(t_i) + \sum_{j=0}^{i-1} K(t_i - t_j) \hat{N}_{j,j+1}^n,$$

for $i = 0, \dots, n$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Sample trajectories simulated with the iVi scheme with time step 0.001.</figcaption></figure>

To illustrate the convergence of the Hawkes iVi scheme, we will examine:

1. The marginal laws of $N_T$ and $\Lambda_T$.

2. The distribution of the arrival times $\tau_1, \dots, \tau_{N_T}$.

**Laplace transforms.** To illustrate the convergence of the marginal laws, we compute the Laplace transforms $\phi_{N_T}(w) = \exp(wN_T)$ and $\phi_{\Lambda_T}(w) = \exp(w\Lambda_T)$ for different numbers of time steps. The convergence of the Monte Carlo estimator of the characteristic function, based on $10^6$ trajectories and with $w = -\frac{1}{\mathbb{E}[N_T]}$, is shown in Figure 2. The horizontal bar indicates three standard deviations of the estimator. Figure 3 reports the computation time required to simulate $10^5$ trajectories, together with the corresponding estimation error, for the iVi scheme under different discretization steps. The horizontal bars indicate the simulation time of the exact methods. These results demonstrate that the iVi scheme can provide substantial acceleration over Ogata's algorithm and the population-based approach while maintaining high precision.
<figure><img src="image_4.png" /><figcaption>Figure 2: Convergence of the Laplace transform with $w = -\frac{1}{\mathbb{E}[N_T]}$ for $N_T$ (left) and $\Lambda_T$ (right). The x-axis shows the number of time steps, and the y-axis shows the absolute error of the Monte Carlo estimator.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Simulation time for $10^5$ trajectories (y-axis) versus the absolute error of the Monte Carlo estimator of the Laplace transform (x-axis). Crosses represent estimators obtained with the iVi scheme under different discretization steps. Horizontal bars indicate the exact methods: Population (blue) and Ogata (purple).</figcaption></figure>

**Comparison of marginal laws.** The proximity of the exact marginal laws $\mathcal{L}(N_T)$ and $\mathcal{L}(\Lambda_T)$ to their counterparts generated by the iVi scheme is assessed by simulating $N = 10^4$ trajectories with each scheme and comparing the empirical distributions of $N_T$ and $\Lambda_T$ with those obtained from the exact population approach. In particular, we plot both the empirical distributions and the corresponding Q-Q plots against samples generated by the exact method. For the iVi scheme, the Kolmogorov-Smirnov test for $\Lambda_T$ yields a $p$-value equal to 0.92. The empirical distributions of $N_T$ and $\Lambda_T$ are shown in Figures 4 and 5, respectively.
<figure><img src="image_6.png" /><figcaption>Figure 4: Marginal law of $N_T$. Empirical densities (top left), distances to the empirical CDF of the population method (top right), Q-Q plots for Ogata algorithm (bottom left) and iVi scheme (bottom right). A time step 0.01 was taken for the iVi scheme.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Marginal law of $\Lambda_T$. Empirical densities (top left), distances to the empirical CDF of the population method (top right), Q-Q plots for Ogata algorithm (bottom left) and iVi scheme (bottom right). A time step 0.01 was taken for the iVi scheme.</figcaption></figure>

**Distribution of arrival times.** In the iVi scheme, the arrival times are simulated uniformly based on the simulated increments $\hat{N}_{i,i+1}^n$ over each discretization interval. Essentially, this means that the instantaneous
jump intensity $\lambda$ is approximated on each interval $[t_i, t_{i+1}]$ by

$$ \hat{\lambda}_i^n = \frac{\hat{\Lambda}_{i,i+1}^n}{t_{i+1}^n - t_i^n} $$

To assess the quality of this approximation, we perform standard time-change tests, see, e.g., Ogata (1988). Specifically, we simulate the arrival times $T^n = (\tau_i^n)_{i=1,...,N_T^n}$, compute

$$ \tau_i^{n,*} = \Lambda_{\tau_i^n}, \quad i = 1, \dots, N_T^n, $$

using (1.2), and verify that $(\tau_i^{n,*})_{i=1,...,N_T^n}$ follow the law of a standard Poisson process (by Theorem 4.1). To this end, we compare the empirical distribution of the interarrival times $(\tau_{i+1}^{n,*} - \tau_i^{n,*})_{i=1,...,N_T^n-1}$ with the exponential distribution Exp(1), corresponding to the waiting times of a Poisson process. In addition, we plot the pairs $(e^{-\tau_i^{n,*}}, e^{-\tau_{i+1}^{n,*}})_{i=1,...,N_T^n-1}$, which should be independent and uniformly distributed on $[0, 1]$.

For the iVi scheme with discretization step 0.1 and $T = 10$, the Kolmogorov-Smirnov test for the interarrival times $(\tau_{i+1}^{n,*} - \tau_i^{n,*})_{i=1,...,N_T^n-1}$ yields a p-value of 0.43. The empirical distribution of the transformed arrivals over $[0, 10]$ is shown in Figure 6.

<figure><img src="image_8.png" /><figcaption>Figure 6: Transformed arrival times generated by the population approach (top), Ogata's algorithm (middle), and the iVi scheme (bottom). The columns correspond to: Q-Q plots against the exact law (left), empirical CDFs (middle), and scatter plots of $(e^{-\tau_i^{n,*}}, e^{-\tau_{i+1}^{n,*}})_{i=1,...,N_T^n-1}$ (right).</figcaption></figure>
**Importance of the implicit step.** We conclude this subsection by illustrating the relevance of the Inverse Gaussian simulation step (1.5). The motivation for introducing the random variable $\xi_i^n$ lies in the implicit discretization scheme for the characteristic function described in Subsection 2.2. A simpler alternative would be to use an explicit scheme: first approximate $U_{t_{i+1}} - U_{t_i} \approx \alpha_i^n$, then simulate $\hat{N}_{i,i+1} \sim P(\alpha_i^n)$, and finally update the integrated intensity using (1.6). We refer to this approach as the explicit integrated intensity scheme.

Note that this scheme can be viewed as a refined version of the discretization method presented in (Kirchner, 2017, Theorem 2.4), since we use the integrated kernel instead of discretizing the kernel itself. However, in high-activity regimes this scheme proves to be significantly less efficient, as we illustrate below for the fractional kernel

$$K_H(t) = \frac{ct^{H-0.5}}{\Gamma(H+0.5)}, \quad t \ge 0,$$

with $H = 0.1$, $c = 0.1$ and $\mu = 5$. We plot the convergence results for the Laplace transforms

$$\phi_{N_T}(w) = \exp(wN_T) \quad \text{and} \quad \phi_{\Lambda_T}(w) = \exp(w\Lambda_T),$$

with $w = -\frac{1}{\mathbb{E}[N_T]}$ and $T = 30$, in Figure 7. We observe that in this case, the iVi scheme is more accurate and exhibits a higher empirical rate of convergence than the explicit integrated intensity scheme, which justifies the use of the implicit approach we propose.

<figure><img src="image_9.png" /><figcaption>Figure 7: Convergence of the Laplace transform with $w = -\frac{1}{\mathbb{E}[N_T]}$ for $N_T$ (left) and $\Lambda_T$ (right). The x-axis shows the number of time steps, and the y-axis shows the absolute error of the Monte Carlo estimator.</figcaption></figure>

## 4.3 A resolvent version of the iVi scheme

In this subsection, we show how the numerical performance of the iVi scheme for Hawkes processes can be significantly improved by using the resolvent of its kernel $K$, i.e., the kernel $R \in L^1([0, T], \mathbb{R})$ satisfying

$$R * K = R - K.$$

We refer to Appendix A for the precise definition and a list of commonly used kernels along with their resolvents.

Recall that the dynamics of the integrated intensity are given by

$$\Lambda_t = \int_0^t g_0(s) ds + \int_0^t K(t-s) N_s ds, \quad t \le T. \qquad (4.1)$$

The following proposition allows us to rewrite (4.1) in terms of the resolvent $R$.

**Proposition 4.2.** The integrated intensity $\Lambda$ satisfies

$$\Lambda_t = G_0^R(t) + \int_0^t R(t-s) Z_s ds, \qquad (4.2)$$
where $Z_t := N_t - \Lambda_t$ and

$$G_0^R(t) := \int_0^t g_0(s) ds + \int_0^t R(t-s) \int_0^s g_0(u) du ds = G_0(t) + \int_0^t R(t-s) G_0(s) ds.$$

*Proof.* Note that equation (4.1) can be rewritten as a convolution:

$$\Lambda = G_0 + K * N,$$

where $G_0(t) = \int_0^t g_0(u) du$. Convolving both sides with $(\text{Id} + R)$ on the left and applying the resolvent equation yields:

$$\begin{aligned} \Lambda + R * \Lambda &= G_0 + R * G_0 + K * N + R * K * N \\ &= G_0 + R * G_0 + R * N. \end{aligned}$$

It remains to observe that $G_0^R = G_0 + R * G_0$, which concludes the proof. □

The equation (4.2) suggests the following modification of the iVi scheme.

**Algorithm 4 - The Hawkes Resolvent iVi scheme:** Simulation of $\hat{\Lambda}, \hat{N}$

1: **Input:** $G_0^R, R$, uniform partition of $[0, T]$, $t_i^n = iT/n$, with $i = 0, \dots, n$.

2: **Output:** $\hat{\Lambda}_{i,i+1}^n, \hat{N}_{i,i+1}^n$ for $i = 0, \dots, n-1$ and (optional) $T^n = (\tau_k^n)_{k=1,\dots,N_T^n}$.

3: Compute

$$r_j^n = \int_{t_j^n}^{t_{j+1}^n} R(s) ds, \quad j = 0, \dots, n-1. \qquad (4.3)$$

4: **for** $i = 0$ to $n-1$ **do**

5: Compute the increment:

$$\alpha_i^n = \max \left\{ (G_0^R(t_{i+1}^n) - G_0^R(t_i^n)) + \sum_{j=0}^{i-1} r_{i-j}^n \hat{Z}_{j,j+1}^n, 0 \right\}, \qquad (4.4)$$

with $\alpha_0^n = G_0^R(t_1^n)$, for $i = 0$.

6: Simulate the increment:

$$\hat{N}_{i,i+1}^n \sim \mathcal{P}(\xi_i^n), \quad \text{where} \quad \xi_i^n \sim \text{IG}\left(\alpha_i^n, \left(\frac{\alpha_i^n}{r_0^n}\right)^2\right). \qquad (4.5)$$

7: Compute the increments:

$$\hat{Z}_{i,i+1}^n = \hat{N}_{i,i+1}^n - \xi_i^n, \qquad (4.6)$$

$$\hat{\Lambda}_{i,i+1}^n = \frac{\alpha_i^n + r_0^n \hat{N}_{i,i+1}^n}{1 + r_0^n}. \qquad (4.7)$$

8: Optional step (jump times): Simulate $\hat{N}_{i,i+1}^n$ i.i.d. uniform random variables on $[t_i, t_{i+1})$, sort, and append them to $T^n$.

9: **end for**

The motivation for using the resolvent version of the scheme lies again in the characteristic function of $\Lambda$. Indeed, as in Section 2, we have, for $w \in \mathbb{C}$ with $\Re(w) \le 0$,

$$\mathbb{E} \left[ \exp \left( w \Lambda_{t_i^n, t_{i+1}^n} \right) \mid \mathcal{F}_{t_i^n} \right] = \exp \left( \int_{t_i^n}^{t_{i+1}^n} F^R \left( \Psi(t_{i+1}^n - s) \right) dG_{t_i^n}^R(s) \right), \quad i \le n-1, \quad n \ge 1,$$
where

$$F^R(u) = w + e^u - 1 - u, \quad u \in \mathbb{C}, \qquad (4.8)$$

$\Psi$ is the solution, with nonpositive real part, to the Volterra equation

$$\Psi(t) = \int_{0}^{t} R(t-s) F^{R}(\Psi(s)) ds, \quad t \le T,$$

and $G_{t_i^n}^R$ is given by

$$G_{t_i^n}^R(s) = G_0^R(s) - G_0^R(t_i^n) + \int_{t_i^n}^{s} \int_{0}^{t_i^n} R(u-r) d(N_r - \Lambda_r), \quad s \ge t_i^n, \quad i \le n-1, \quad n \ge 1.$$

Approximating $F^R$ by a quadratic function gives

$$F^R(u) = w + e^u - 1 - u \approx w + \frac{u^2}{2}, \quad u \in \mathbb{C}. \qquad (4.9)$$

Comparing (2.6) and (4.9) explains why the left-endpoint approximations

$$\int_{t_i^n}^{t_{i+1}^n} F^R(\Psi(t_{i+1}^n - s)) dG_{t_i^n}^R(s) \approx F^R(\Psi(\Delta t)) G_{t_i^n}^R(t_{i+1}^n),$$

and

$$\Psi(\Delta t) \approx \int_{0}^{\Delta t} R(s) ds F^R(\Psi(\Delta t)) = r_0^n F^R(\Psi(\Delta t)),$$

are more accurate than their counterparts in the iVi scheme given by (2.4) and (2.5). Indeed, while the function $F$ defined in (2.3) behaves like $F(u) \sim w + u$, the function $F^R$ is approximately quadratic: $F^R(u) \sim w + \frac{u^2}{2}$ as $u \to 0$, which increases the approximation order. We will show that this dramatically improves the speed of convergence.

The next steps of the proof are similar to those in Section 2 and lead to Algorithm 4.3.

Again, as in the scheme of Algorithm 1, the increment $\hat{\Lambda}_{i,i+1}^n$ is computed using (4.7), which follows from the relation

$$\Lambda_{t_{i+1}^n} - \Lambda_{t_i^n} \approx \alpha_i^n + r_0^n (\hat{N}_{i,i+1}^n - \hat{\Lambda}_{i,i+1}^n) = \hat{\Lambda}_{i,i+1}^n,$$

arising from a discretization similar to (2.2), see the discussion at the end of Subsection 2.2.

Importantly, for the computation of $\alpha_i^n$ in (4.4), one must use

$$\hat{Z}_{i,i+1}^n = \hat{N}_{i,i+1}^n - \xi_i^n,$$

and not

$$\hat{Z}_{i,i+1}^n = \hat{N}_{i,i+1}^n - \hat{\Lambda}_{i,i+1}^n,$$

since the latter choice significantly deteriorates the performance of the scheme. However, this leads to a potential loss of positivity in the discretization of (4.4). Although the integral $G_{t_i^n}^R(t_{i+1}^n)$ is positive, its discretization may not preserve positivity because the increments $\hat{Z}_{i,i+1}^n$ are no longer guaranteed to be positive and are difficult to control. While the capping in (4.4) does not affect the numerical results, it complicates the theoretical analysis of convergence.

Another inconvenience of the scheme is the need to compute the resolvent of the kernel $K$, which may be a challenging problem in some cases. However, when the resolvent is computable, the resolvent iVi scheme exhibits significantly better performance than the standard scheme given by Algorithm 1, as follows from the numerical experiments for the fractional kernel

$$K_H(t) = \frac{ct^{H-0.5}}{\Gamma(H+0.5)}, \quad t \ge 0,$$
with $H = 0.1$, $c = 0.1$ and $\mu = 5$. The convergence results for the Laplace transforms

$$ \phi_{N_T}(w) = \exp(wN_T) \quad \text{and} \quad \phi_{\Lambda_T}(w) = \exp(w\Lambda_T), $$

with $w = -\frac{1}{\mathbb{E}[N_T]}$ and $T = 30$, as well as the corresponding execution times, are shown in Figures 8 and 9.

For this example, the resolvent scheme achieves a precision of three standard deviations already with 80 time steps. In this case, the execution time for the iVi Resolvent scheme is approximately 1100 times faster than for the Population approach. Such a substantial difference occurs when the average number of jumps per trajectory becomes sufficiently large relative to the number of time steps needed for the iVi Resolvent scheme to reach convergence. Our experiments indicate that this is typically the case, with the iVi Resolvent scheme outperforming the exact methods across different choices of kernels and simulation horizons.

<figure><img src="image_10.png" /><figcaption>Figure 8: Convergence of the Laplace transform with $w = -\frac{1}{\mathbb{E}[N_T]}$ for $N_T$ (left) and $\Lambda_T$ (right). The x-axis shows the number of time steps, and the y-axis shows the absolute error of the Monte Carlo estimator.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Simulation time for $10^5$ trajectories (y-axis) versus the absolute error of the Monte Carlo estimator of the Laplace transform (x-axis). Crosses represent estimators obtained with the iVi scheme under different discretization steps. Horizontal bars indicate the exact methods: Population (blue) and Ogata (purple).</figcaption></figure>

## 4.4 Markovian version of the iVi scheme

In this subsection, we discuss a modification of Algorithms 1 and 4 that can be applied when the Hawkes process is Markovian.
For an arbitrary kernel $K$, the complexity of the simulation scheme is quadratic due to the computation of $\alpha_i^n$ given by (1.4) and (4.4). However, in the exponential case

$$K_{\exp}^{b,c}(t) = ce^{-bt}, \quad b > 0, c > 0, t \ge 0,$$

we have

$$k_{i+1}^n = e^{-bT/n} k_i^n, \quad i = 0, \dots, n-1, \qquad (4.10)$$

so that

$$\begin{aligned} \alpha_i^n &= \int_{t_i^n}^{t_{i+1}^n} g_0(s) ds + \sum_{j=0}^{i-1} k_{i-j}^n \hat{N}_{j,j+1}^n \\ &= \int_{t_i^n}^{t_{i+1}^n} g_0(s) ds + e^{-bT/n} \left( \alpha_{i-1}^n - \int_{t_{i-1}^n}^{t_i^n} g_0(s) ds \right) + k_1^n \hat{N}_{i-1,i}^n, \quad i = 1, \dots, n-1. \end{aligned}$$

This formula leads to a scheme with linear complexity, as all information about the past is summarized in the previous value $\alpha_{i-1}^n$. The same argument applies to the resolvent scheme, since the resolvent of the exponential kernel is itself an exponential kernel.

Similarly, a Markovian version of Ogata's Algorithm 3 can be easily derived by keeping track of the value

$$Y(t) = \sum_{\tau_i \le t} K_{\exp}^{b,c}(t - \tau_i),$$

which satisfies $Y(t+h) = e^{-bh}Y(t)$, $h > 0$. Moreover, in the Markovian case, there exists an exact algorithm to simulate the next arrival $\tau_{i+1}$ conditionally on $\mathcal{F}_{\tau_i}$, see Dassios and Zhao (2013).

We illustrate the convergence of the Laplace transform of $N_T$ and $\Lambda_T$ for $\mu = 10$ and $T = 2$, corresponding to the exponential kernel with $b = 5$ and $c = 4$. The results for the exact scheme, the Markovian version of Ogata's algorithm, and the iVi schemes are presented in Figures 10 and 11. Again, the iVi scheme achieves substantial acceleration (about 50 times faster than the exact method for the Resolvent iVi scheme with 200 time steps), even compared with the optimal exact method, thanks to time discretization and vectorization of the scheme.

<figure><img src="image_12.png" /><figcaption>Figure 10: Convergence of the Laplace transform with $w = -\frac{1}{E[N_T]}$ for $N_T$ (left) and $\Lambda_T$ (right). The x-axis shows the number of time steps, and the y-axis shows the absolute error of the Monte Carlo estimator.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11: Simulation time for $10^5$ trajectories (y-axis) versus the absolute error of the Monte Carlo estimator of the Laplace transform (x-axis). Crosses represent estimators obtained with the iVi scheme under different discretization steps. Horizontal bars indicate the exact methods: Population (blue) and Ogata (purple).</figcaption></figure>

**Remark 4.3.** A similar reasoning applies when the kernel is given by a sum of exponentials:

$$K(t) = \sum_{k=1}^{m} c_k e^{-b_k t} = \sum_{k=1}^{m} K_{\text{exp}}^{b_k, c_k}(t), \quad b_k > 0, c_k > 0, t \ge 0. \quad (4.11)$$

In this case, we can decompose $k_i^n$ as

$$k_i^n = \sum_{k=1}^{m} k_i^{n,(k)}, \quad i = 0, \dots, n-1,$$

where $k_i^{n,(k)}$ corresponds to the kernel $K_{\text{exp}}^{b_k, c_k}$ and satisfies the recursion (4.10).

The Markovian state is then given by the vector $y_i = (y_{i,1}, \dots, y_{i,m})$, with

$$y_{i,k} = \sum_{j=0}^{i-1} k_{i-j}^{n,(k)} \hat{N}_{j,j+1}^{n},$$

which evolves according to

$$y_{i+1,k} = e^{-b_k T/n} y_{i,k} + k_1^{n,(k)} \hat{N}_{i,i+1}^{n},$$

and the increment $\alpha_i^n$ is then given by

$$\alpha_i^n = \int_{t_i^n}^{t_{i+1}^n} g_0(s) ds + \sum_{k=1}^{m} y_{i,k}.$$

Sum of exponential kernels as in (4.11) can be used to approximate other kernels K, giving rise to so-called multifactor Markovian approximations of the Hawkes process with kernel K; see, for instance, (Abi Jaber, Cuchiero, Larsson, and Pulido, 2021, Section 7.2). This approach yields a computational complexity of order $O(mn)$; thus, when the number of factors m is smaller than the number of time steps n, it offers a significant computational advantage.

# 5 Proof of weak convergence

This section is dedicated to the proof of convergence of the Hawkes iVi scheme, as stated in Theorem 3.1.
## 5.1 Preliminary results

We begin with some basic properties on the moments of the random variables $(\alpha, \xi, \hat{N}, \hat{\Lambda})$ produced by Algorithm 1. We recall that the latter leads to the following structure at step $i \le n-1$:

* The quantity $\alpha_i^n$ is given by (1.4).

* Conditional on $\alpha_i^n$, $\xi_i^n$ follows an Inverse Gaussian distribution with parameters $\left(\frac{\alpha_i^n}{1-k_0^n}, \left(\frac{\alpha_i^n}{k_0^n}\right)^2\right)$.

* Conditional on $\xi_i^n$, $\hat{N}_{i,i+1}^n$ follows a Poisson distribution with parameter $\xi_i^n$.

* The increment $\hat{\Lambda}_{i,i+1}^n$ is defined as $\alpha_i^n + k_0^n \hat{N}_{i,i+1}^n$.

The filtrations $(\mathcal{F}_t^n)_{t \le T, n \ge 1}$ are defined in (3.1).

**Lemma 5.1.** Let $n \ge 1$. Then, the random variables $(\alpha_i^n)_{i \le n-1}$, $(\xi_i^n)_{i \le n-1}$, $(\hat{N}_{i,i+1}^n)_{i \le n-1}$ and $(\hat{\Lambda}_{i,i+1}^n)_{i \le n-1}$ have finite moments of order 2. Moreover, we have

$$ \mathbb{E} \left[ \hat{N}_{i,i+1}^n \mid \mathcal{F}_{t_j}^n \right] = \mathbb{E} \left[ \hat{\Lambda}_{i,i+1}^n \mid \mathcal{F}_{t_j}^n \right], \quad j \le i \le n-1. \quad (5.1) $$

$$ \mathbb{E} \left[ \left( \hat{N}_{i,i+1}^n - \hat{\Lambda}_{i,i+1}^n \right)^2 \mid \mathcal{F}_{t_j}^n \right] = \left( (k_0^n)^2 + (1-k_0^n)^2 \right) \mathbb{E} \left[ \hat{\Lambda}_{i,i+1}^n \mid \mathcal{F}_{t_j}^n \right], \quad j \le i \le n-1. \quad (5.2) $$

*Proof.* We prove the existence of second order moments by induction on $i \le n-1$. Firstly, $\hat{\alpha}_0^n = \int_0^n \frac{g_0(s)}{s} ds$ is deterministic. Moreover $\xi_0^n \sim IG\left(\frac{\alpha_0^n}{1-k_0^n}, \left(\frac{\alpha_0^n}{k_0^n}\right)^2\right)$ has finite moments of any order (see Appendix B). Then, we have

$$ \mathbb{E} \left[ \left( \hat{N}_{0,1}^n \right)^2 \right] = \mathbb{E} \left[ \mathbb{E} \left[ \left( \hat{N}_{0,1}^n \right)^2 \mid \xi_0^n \right] \right] = \mathbb{E} \left[ \xi_0^n + (\xi_0^n)^2 \right] < +\infty, $$

since $\hat{N}_{0,1}^n$ follows a Poisson distribution with parameter $\xi_0^n$, conditional on the latter. Finally, by linear combination,

$$ \hat{\Lambda}_{0,1}^n = \alpha_0^n + k_0^n \hat{N}_{0,1}^n $$

also admits a finite moment of order 2.

We next fix $i \le n-1$ and assume that $(\alpha_j^n)_{j<i}$, $(\xi_j^n)_{j<i}$, $(\hat{N}_{j,j+1}^n)_{j<i}$ and $(\hat{\Lambda}_{j,j+1}^n)_{j<i}$ all have finite moments of order 2. Then it is also the case for

$$ \alpha_i^n = \int_{t_i^n}^{t_{i+1}^n} g_0(s) ds + \sum_{j=0}^{i-1} k_{i-j}^n \hat{N}_{j,j+1}^n. $$

Using the formula (B.2) for the moments of the Inverse Gaussian distribution, we obtain

$$ \mathbb{E} \left[ (\xi_i^n)^2 \right] = \mathbb{E} \left[ \mathbb{E} \left[ (\xi_i^n)^2 \mid \alpha_i^n \right] \right] = \mathbb{E} \left[ \frac{(k_0^n)^2}{(1-k_0^n)^3} \alpha_i^n + \left( \frac{\alpha_i^n}{1-k_0^n} \right)^2 \right] < +\infty. $$

Finally, similarly to the case $i=0$, $\hat{N}_{i,i+1}^n$ and $\hat{\Lambda}_{i,i+1}^n$ have finite second-order moments.

For (5.1), we use the tower property of conditional expectation,

$$ \begin{aligned} \mathbb{E} \left[ \hat{N}_{i,i+1}^n \mid \mathcal{F}_{t_j}^n \right] &= \mathbb{E} \left[ \mathbb{E} \left[ \hat{N}_{i,i+1}^n \mid \mathcal{F}_{t_j}^n, \xi_i^n \right] \mid \mathcal{F}_{t_j}^n \right] \\ &= \mathbb{E} \left[ \xi_i^n \mid \mathcal{F}_{t_j}^n \right] \\ &= \mathbb{E} \left[ \mathbb{E} \left[ \xi_i^n \mid \mathcal{F}_{t_j}^n, \alpha_i^n \right] \mid \mathcal{F}_{t_j}^n \right] \\ &= \frac{1}{1-k_0^n} \mathbb{E} \left[ \alpha_i^n \mid \mathcal{F}_{t_j}^n \right], \quad j \le i \le n-1. \end{aligned} \quad (5.3) $$
Inserting this into (1.6), we obtain

$$
\begin{align*}
\mathbb{E}[\hat{\Lambda}_{i,i+1}^n | \mathcal{F}_{t_j^n}] &= \mathbb{E}\left[\left(\alpha_i^n \left(1 + \frac{k_0^n}{1-k_0^n}\right)\right) | \mathcal{F}_{t_j^n}\right] \\
&= \frac{1}{1-k_0^n} \mathbb{E}\left[\alpha_i^n | \mathcal{F}_{t_j^n}\right] \\
&= \mathbb{E}[\hat{N}_{i,i+1}^n | \mathcal{F}_{t_j^n}], \quad j \le i \le n-1,
\end{align*}
$$

thus proving (5.1).

We now move to the second equality. Using (1.6), we have

$$
\mathbb{E}\left[\left(\hat{N}_{i,i+1}^n - \hat{\Lambda}_{i,i+1}^n\right)^2 \middle| \mathcal{F}_{t_j^n}\right] = \mathbb{E}\left[\left((1-k_0^n) \hat{N}_{i,i+1}^n - \alpha_i^n\right)^2 \middle| \mathcal{F}_{t_j^n}\right], \quad j \le i \le n-1.
$$

An application of the tower property of conditional expectation, along with the expression of second-order moments of Poisson and Inverse Gaussian (see (B.2)) variables, yields that

$$
\begin{align*}
\mathbb{E}\left[\left(\hat{N}_{i,i+1}^n - \Lambda_{i,i+1}^n\right)^2 \middle| \mathcal{F}_{t_j^n}\right] &= \mathbb{E}\left[\mathbb{E}\left[(1-k_0^n)^2 \left(\hat{N}_{i,i+1}^n\right)^2 - 2(1-k_0^n) \hat{N}_{i,i+1}^n \alpha_i^n + (\alpha_i^n)^2 \middle| \mathcal{F}_{t_j^n}, \xi_i^n, \alpha_i^n\right] \middle| \mathcal{F}_{t_j^n}\right] \\
&= \mathbb{E}\left[(1-k_0^n)^2 \left(\left(\xi_i^n\right)^2 + \xi_i^n\right) - 2(1-k_0^n) \xi_i^n \alpha_i^n + (\alpha_i^n)^2 \middle| \mathcal{F}_{t_j^n}\right] \\
&= \mathbb{E}\left[\mathbb{E}\left[(1-k_0^n)^2 \left(\left(\xi_i^n\right)^2 + \xi_i^n\right) - 2(1-k_0^n) \xi_i^n \alpha_i^n + (\alpha_i^n)^2 \middle| \mathcal{F}_{t_j^n}, \alpha_i^n\right] \middle| \mathcal{F}_{t_j^n}\right] \\
&= \mathbb{E}\left[\frac{(k_0^n)^2}{1-k_0^n} \alpha_i^n + (\alpha_i^n)^2 + (1-k_0^n) \alpha_i^n - 2(\alpha_i^n)^2 + (\alpha_i^n)^2 \middle| \mathcal{F}_{t_j^n}\right] \\
&= \frac{(k_0^n)^2 + (1-k_0^n)^2}{1-k_0^n} \mathbb{E}\left[\alpha_i^n \middle| \mathcal{F}_{t_j^n}\right], \quad j \le i \le n-1,
\end{align*}
$$

which combined with (5.1) and (5.3) gives (5.2).

□

The following lemma identifies a martingale that will be pivotal in the proofs.

**Lemma 5.2.** For $n \ge 1$, the process $(Z^n)_{t \le T} := (N_t^n - \Lambda_t^n)_{t \le T}$ is a square-integrable $(\mathcal{F}_t^n)_{t \le T}$-martingale.
Its quadratic variation is given by

$$
[Z^n]_t = \sum_{i=0}^{\lfloor nt/T \rfloor - 1} (\hat{N}_{i,i+1}^n - \hat{\Lambda}_{i,i+1}^n)^2, \quad t \le T.
$$

Moreover, $([Z^n] - [(k_0^n)^2 + (1-k_0^n)^2] \Lambda^n)$ is a martingale, which implies

$$
\mathbb{E}[[Z^n]_t] = ((k_0^n)^2 + (1-k_0^n)^2) \mathbb{E}[\Lambda_t^n], \quad t \le T. \tag{5.4}
$$

Proof. We recall that square-integrability of all the quantities is provided by Lemma 5.1. Let $s \le t \le T$, we define $i := \lfloor nt/T \rfloor$ and $j := \lfloor ns/T \rfloor$. We have

$$
\begin{align*}
\mathbb{E}[N_t^n | \mathcal{F}_s^n] &= \mathbb{E} \left[ \sum_{l=0}^{i-1} \hat{N}_{l,l+1}^n \middle| \mathcal{F}_{t_j^n} \right] \\
&= N_s^n + \sum_{l=j}^{i-1} \mathbb{E} \left[ \hat{N}_{l,l+1}^n \middle| \mathcal{F}_{t_j^n} \right] \\
&= N_s^n + \sum_{l=j}^{i-1} \mathbb{E} \left[ \hat{\Lambda}_{l,l+1}^n \middle| \mathcal{F}_{t_j^n} \right] \\
&= N_s^n + \mathbb{E} \left[ \Lambda_t^n - \Lambda_s^n \middle| \mathcal{F}_s^n \right],
\end{align*}
$$
where we used (5.1) to compute the conditional expectations. This proves the martingality of $Z^n = N^n - \Lambda^n$. To identify its quadratic variation, we notice that it has bounded variation on $[0, T]$ (as it jumps a finite number of times) and is therefore a quadratic pure-jump martingale, see Protter (1990, Theorem II.26). Since it also starts from 0, we deduce that its quadratic variation is the sum of the squares of its jumps, i.e.

$$[Z^n]_t = [N^n - \Lambda^n]_t = \sum_{i=0}^{\lfloor nt/T \rfloor - 1} (\hat{N}_{i,i+1}^n - \hat{\Lambda}_{i,i+1}^n)^2, \quad t \le T.$$

By taking the expectation, applying (5.2) and recalling the definition of $\Lambda^n$ in (1.7), we obtain

$$\begin{aligned} \mathbb{E}[[Z^n]_t - [Z^n]_s | \mathcal{F}_s^n] &= \sum_{l=j}^{i-1} \mathbb{E} \left[ \left( \hat{N}_{l,l+1}^n - \hat{\Lambda}_{l,l+1}^n \right)^2 \middle| \mathcal{F}_{t_j}^n \right] \\ &= \left( (k_0^n)^2 + (1-k_0^n)^2 \right) \mathbb{E} [\Lambda_t^n - \Lambda_s^n | \mathcal{F}_s^n], \end{aligned}$$

proving that $([Z^n] - [(k_0^n)^2 + (1-k_0^n)^2] \Lambda^n)$ is a martingale.

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

A very convenient feature of our scheme is that it can be rewritten as a stochastic Volterra equation with a measure-valued kernel, which will be useful for deriving estimates on our processes. This is formulated in the following lemma.

**Lemma 5.3.** For $n \ge 1$, define the nonnegative discrete measure $K^n := \sum_{i=0}^{n-1} k_i^n \delta_{t_i^n}$ on $([0, T], \mathcal{B}_{[0, T]}).$
Then, we have

$$\Lambda_t^n = \int_0^{\lfloor \frac{n t}{T} \rfloor} g_0(s) ds + \int_{[0,t]} N_{t-s}^n K^n(ds), \quad t \le T. \quad (5.5)$$

*Proof.* We recall that combining (1.4) and (1.6) leads to

$$\hat{\Lambda}_{i,i+1}^n = \int_{t_i^n}^{t_{i+1}^n} g_0(s) ds + \sum_{j=0}^i k_{i-j}^n \hat{N}_{j,j+1}^n, \quad i \le n-1, \ n \ge 1.$$

Let $t \le T$ and set $i := \lfloor nt/T \rfloor$, then

$$\begin{aligned} \Lambda_t^n &= \sum_{j=0}^{i-1} \hat{\Lambda}_{j,j+1}^n \\ &= \int_0^{t_i^n} g_0(s) ds + \sum_{j=0}^{i-1} \sum_{l=0}^j k_{j-l}^n \hat{N}_{l,l+1}^n \\ &= \int_0^{t_i^n} g_0(s) ds + \sum_{j=0}^{i-1} \sum_{l=0}^j k_l^n \hat{N}_{j-l,j-l+1}^n \\ &= \int_0^{t_i^n} g_0(s) ds + \sum_{l=0}^{i-1} k_l^n \sum_{j=l}^{i-1} \hat{N}_{j-l,j-l+1}^n \\ &= \int_0^{t_i^n} g_0(s) ds + \sum_{l=0}^{i-1} k_l^n N_{t_i^n-l}^n. \end{aligned}$$

From the definition of $N^n$ in (1.8), and the fact that our partition of $[0, T]$ is uniform, we have $N_{t_i^n-l}^n = N_{t-t_i^n}^n$. Since $N_0^n = 0$, we deduce that

$$\sum_{l=0}^{i-1} k_l^n N_{t_i^n-l}^n = \int_{[0,t]} N_{t-s}^n K^n(ds),$$

which concludes the proof of (5.5).

<figure><img src="image_15.png" /></figure>
## 5.2 Estimates

We are now ready to state the estimates on the process $\Lambda^n$ that will be central in proving the tightness of our processes. We make use of the notation $G_0 = \int_0 g_0(s) ds$.

**Lemma 5.4.** We have the following moment bounds,

$$C_{\Lambda} := \sup_{n \ge 1} \mathbb{E} \left[ \sup_{0 \le t \le T} \Lambda_t^n \right] < +\infty, \qquad (5.6)$$

$$C_{\Lambda^2} := \sup_{n \ge 1} \mathbb{E} \left[ \sup_{0 \le t \le T} (\Lambda_t^n)^2 \right] < +\infty. \qquad (5.7)$$

Moreover, for any $j < i \le n$,

$$\Lambda_{t_i^n}^n - \Lambda_{t_j^n}^n \le G_0(t_i^n) - G_0(t_j^n) + \|N^n\|_{\infty} \left( \int_0^{t_{j+1}^n} |K(s+t_i^n-t_j^n) - K(s)| ds + \int_0^{t_i^n-t_j^n} K(s) ds \right). \quad (5.8)$$

*Proof.* We first prove (5.7), which implies (5.6) using $|x| \le 1+x^2$. Fix $t \le T$ and $n \ge 1$, we recall that thanks to the Volterra equation (5.5),

$$\Lambda_t^n \le G_0(T) + \int_{[0,t]} N_{t-s}^n K^n(ds).$$

Taking squares and applying Jensen's inequality on the measure $\frac{K^n(ds)}{K^n([0,t])}$, we have

$$(\Lambda_t^n)^2 \le 2 G_0^2(T) + 2 \int_0^T K(s) ds \int_{[0,t]} (N_{t-s}^n)^2 K^n(ds),$$

where we used the fact that $K^n([0,t]) \le K^n([0,T]) = \int_0^T K(s) ds$. We can then leverage the martingale property of $Z^n = N^n - \Lambda^n$ from Lemma 5.2 to apply the BDG inequality. We obtain

$$\begin{aligned} \mathbb{E}[(N_t^n)^2] &\le 2 \mathbb{E}[(Z_t^n)^2] + 2 \mathbb{E}[(\Lambda_t^n)^2] \\ &\le 8 \mathbb{E}[[Z^n]_t] + 2 \mathbb{E}[(\Lambda_t^n)^2] \\ &\le 8 \left( (k_0^n)^2 + (1-k_0^n)^2 \right) \mathbb{E}[\Lambda_t^n] + 2 \mathbb{E}[(\Lambda_t^n)^2] \\ &\le 8 \left( (k_0^n)^2 + (1-k_0^n)^2 \right) \mathbb{E}[1+(\Lambda_t^n)^2] + 2 \mathbb{E}[(\Lambda_t^n)^2]. \end{aligned}$$

This, along with the fact that $k_0^n \to 0$ as $n \to \infty$, proves that we can find a constant $C > 0$, independent of $t$ and $n$, such that

$$\mathbb{E}[(\Lambda_t^n)^2] \le C + C \int_{[0,t]} \mathbb{E}[(\Lambda_{t-s}^n)^2] K^n(ds).$$

We now choose $\theta > 0$ such that

$$C \int_0^T e^{-\theta s} K(s) ds \le \frac{1}{2}.$$

Such a number exists by the dominated convergence theorem. It allows us to bound

$$\begin{aligned} \int_{[0,T]} e^{-\theta s} K^n(ds) &= \sum_{i=0}^{n-1} e^{-\theta t_i^n} \int_{t_i^n}^{t_{i+1}^n} K(s) ds \\ &= \sum_{i=0}^{n-1} \int_{t_i^n}^{t_{i+1}^n} e^{\theta(s-t_i^n)} e^{-\theta s} K(s) ds \\ &\le e^{\theta T/n} \int_0^T e^{-\theta s} K(s) ds \\ &\le \frac{1}{2C} e^{\theta T/n}, \quad n \ge 1. \end{aligned}$$
Thus, there exists a rank $n_1$, such that

$$ \int_{[0,T]} e^{-\theta s} K^n(ds) \leq \frac{3}{4C}, \quad n \geq n_1. $$

Hence, we obtain

$$ \begin{aligned} \sup_{0 \le t \le T} e^{-\theta t} \mathbb{E}[(\Lambda_t^n)^2] & \le C + C \sup_{0 \le t \le T} \int_{[0,t]} e^{-\theta(t-s)} \mathbb{E}[(\Lambda_{t-s}^n)^2] e^{-\theta s} K^n(ds) \\ & \le C + C \left( \sup_{0 \le t \le T} e^{-\theta t} \mathbb{E}[(\Lambda_t^n)^2] \right) \int_{[0,T]} e^{-\theta s} K^n(ds) \\ & \le C + \frac{3}{4} \sup_{0 \le t \le T} e^{-\theta t} \mathbb{E}[(\Lambda_t^n)^2], \quad n \ge n_1. \end{aligned} $$

For each $n$, the quantity $\sup_{0 \le t \le T} e^{-\theta t} \mathbb{E}[(\Lambda_t^n)^2]$ is finite, as it is bounded by $\mathbb{E}[(\Lambda_T^n)^2]$ which is finite thanks to Lemma 5.1. Therefore, we obtain,

$$ \sup_{0 \le t \le T} e^{-\theta t} \mathbb{E}[(\Lambda_t^n)^2] \le 4C, \quad n \ge n_1, $$

and finally,

$$ \mathbb{E}[(\Lambda_T^n)^2] \le e^{\theta T} \sup_{0 \le t \le T} e^{-\theta t} \mathbb{E}[(\Lambda_t^n)^2] \le 4C e^{\theta T}, \quad n \ge n_1. $$

Since $\Lambda^n$ is non-decreasing, this proves (5.7).

For (5.8), we fix $n \ge 1$ and $j < i \le n$. We have

$$ \begin{aligned} \Lambda_{t_i^n}^n - \Lambda_{t_j^n}^n &= \int_{t_j^n}^{t_i^n} g_0(s) ds + \sum_{l=0}^{i} N_{t_i^n-t_l^n}^n k_l^n - \sum_{l=0}^{j} N_{t_j^n-t_l^n}^n k_l^n \\ &= \int_{t_j^n}^{t_i^n} g_0(s) ds + \sum_{l=0}^{i-j-1} N_{t_i^n-t_l^n}^n k_l^n + \sum_{l=i-j}^{i} N_{t_i^n-t_l^n}^n k_l^n - \sum_{l=0}^{j} N_{t_j^n-t_l^n}^n k_l^n \\ &= \int_{t_j^n}^{t_i^n} g_0(s) ds + \sum_{l=0}^{i-j-1} N_{t_i^n-t_l^n}^n k_l^n + \sum_{l=0}^{j} N_{t_j^n-t_l^n}^n (k_{l+i-j}^n - k_l^n) \\ &\le \int_{t_j^n}^{t_i^n} g_0(s) ds + \|N^n\|_{\infty} \left( \sum_{l=0}^{i-j-1} k_l^n + \sum_{l=0}^{j} |k_{l+i-j}^n - k_l^n| \right). \end{aligned} $$

We conclude by computing

$$ \sum_{l=0}^{i-j-1} k_l^n = \int_0^{t_i^n - t_j^n} K(s) ds, $$

and

$$ \sum_{l=0}^{j} |k_{l+i-j}^{n} - k_{l}^{n}| = \sum_{l=0}^{j} \left| \int_{t_l^n}^{t_{l+1}^n} (K(s+t_i^n-t_j^n) - K(s)) ds \right| \le \int_{0}^{t_{j+1}^n} |K(s+t_i^n-t_j^n) - K(s)| ds. $$

**Remark 5.5.** For any $n \ge 1$, $Z^n = N^n - \Lambda^n$ is a martingale, therefore by the BDG inequality,

$$ \mathbb{E}\left[\sup_{0 \le t \le T} |Z_t^n|^2\right] \le 4 \mathbb{E}\left[[Z^n]_T\right] \le 4 \left((k_0^n)^2 + (1-k_0^n)^2\right) \mathbb{E}[\Lambda_T^n], $$

where we used (5.4). Combining this with the bound (5.6), we obtain

$$ C_Z := \sup_{n \ge 1} \mathbb{E} \left[ \sup_{0 \le t \le T} |Z_t^n|^2 \right] < +\infty, \qquad (5.9) $$

since $k_0^n$ goes to 0 in the limit $n \to \infty$.
We will need one additional property to characterize the limiting processes. Namely, we want to ensure the limit of $(N^n)_{n \ge 1}$ is a counting process, and therefore only has unit jumps.

**Lemma 5.6.** We have

$$ \mathbb{P} \left( \sup_{0 \le t \le T} \Delta N_t^n > 1 \right) \xrightarrow{n \to \infty} 0. $$

*Proof.* Let $n \ge 1$, we have

$$ \mathbb{P} \left( \sup_{0 \le t \le T} \Delta N_t^n > 1 \right) = \mathbb{P} \left( \bigcup_{i=0}^{n-1} \{ \hat{N}_{i,i+1}^n \ge 2 \} \right) \le \sum_{i=0}^{n-1} \mathbb{P} \left( \hat{N}_{i,i+1}^n \ge 2 \right) $$

Then, since $\hat{N}_{i,i+1}^n \sim \mathcal{P}(\xi_i^n)$ conditional on $\xi_i^n$, we have

$$ \mathbb{P}(\hat{N}_{i,i+1}^n \ge 2) = \mathbb{E}[1 - e^{-\xi_i^n} - \xi_i^n e^{-\xi_i^n}], \quad i \le n-1. $$

Using the bound $1 - e^{-x} - xe^{-x} \le x^2$ for $x \ge 0$, we obtain

$$ \mathbb{P}(\hat{N}_{i,i+1}^n \ge 2) \le \mathbb{E}[(\xi_i^n)^2] = \mathbb{E}[\mathbb{E}[(\xi_i^n)^2 | \alpha_i^n]] = \frac{(k_0^n)^2}{(1-k_0^n)^3} \mathbb{E}[\alpha_i^n] + \mathbb{E}[(\alpha_i^n)^2], \quad i \le n-1, $$

where we used (B.2) for computing the second order moment of an Inverse Gaussian variable. We now use the fact that $\alpha_i^n$ is nonnegative (see (1.4)) and $\alpha_i^n \le \hat{\Lambda}_{i,i+1}^n$ thanks to (1.6). Therefore,

$$ \begin{aligned} \mathbb{P}(\hat{N}_{i,i+1}^n \ge 2) & \le \frac{(k_0^n)^2}{(1-k_0^n)^3} \mathbb{E}[\hat{\Lambda}_{i,i+1}^n] + \mathbb{E}[(\hat{\Lambda}_{i,i+1}^n)^2] \\ & \le \frac{(k_0^n)^2}{(1-k_0^n)^3} \mathbb{E}[\hat{\Lambda}_{i,i+1}^n] + \mathbb{E}[\hat{\Lambda}_{i,i+1}^n \max_{0 \le j \le n-1} \hat{\Lambda}_{j,j+1}^n], \quad i \le n-1. \end{aligned} $$

Summing this inequality, we get

$$ \mathbb{P} \left( \sup_{0 \le t \le T} \Delta N_t^n > 1 \right) \le \mathbb{E} \left[ \Lambda_T^n \left( \frac{(k_0^n)^2}{(1-k_0^n)^3} + \max_{0 \le j \le n-1} \hat{\Lambda}_{j,j+1}^n \right) \right]. $$

We treat the two terms on the right-hand side separately. First,

$$ \frac{(k_0^n)^2}{(1-k_0^n)^3} \mathbb{E}[\Lambda_T^n] \le \frac{(k_0^n)^2}{(1-k_0^n)^3} C_\Lambda \xrightarrow{n \to \infty} 0, $$

from (5.6) and the fact that $k_0^n \to 0$ in the limit. Secondly,

$$ \mathbb{E} \left[ \Lambda_T^n \max_{0 \le j \le n-1} \hat{\Lambda}_{j,j+1}^n \right] \le w \left( G_0, \frac{T}{n} \right) \mathbb{E} [\Lambda_T^n] + \mathbb{E} [\Lambda_T^n \|N^n\|_\infty] \left( \int_0^{T-\frac{T}{n}} \left| K \left( t + \frac{T}{n} \right) - K(t) \right| dt + k_0^n \right) $$

from (5.8). Here, $w$ stands for the classical modulus of continuity. The right-hand side vanishes in the limit since

$$ w \left( G_0, \frac{T}{n} \right), \quad k_0^n, \quad \text{and} \quad \int_0^{T-\frac{T}{n}} \left| K \left( t + \frac{T}{n} \right) - K(t) \right| dt $$

all converge to 0 (see Lemma 5.9), and

$$ \sup_{n \ge 1} \mathbb{E} [\Lambda_T^n \|N^n\|_\infty] < +\infty $$

by Lemma 5.4 and Remark 5.5. Thus, we can conclude that

$$ \mathbb{P} \left( \sup_{0 \le t \le T} \Delta N_t^n > 1 \right) \xrightarrow{n \to \infty} 0. $$

□
## 5.3 Tightness and characterization of the limiting processes

We now have all the elements to prove the tightness of our processes in the Skorokhod topology. By C-tightness, we mean tightness with almost surely continuous limits, see Jacod and Shiryaev (2013, Definition VI.3.25).

**Lemma 5.7.** The sequence $(\Lambda^n)_{n \ge 1}$ is C-tight, and $(N^n)_{n \ge 1}$ is tight.

*Proof.* We start by proving C-tightness of $(\Lambda^n)_{n \ge 1}$ using criterion (ii) from Jacod and Shiryaev (2013, Proposition VI.3.26). Firstly, (5.6) implies

$$ \lim_{R \to \infty} \limsup_{n \to \infty} \mathbb{P} \left( \sup_{0 \le t \le T} |\Lambda_t^n| \ge R \right) = 0. $$

Then, using (5.8), and noticing that we can rewrite the modulus of continuity of $\Lambda^n$ as

$$ w(\Lambda^n, \delta) = \max_{\substack{0 \le j \le i \le n \\ t_i^n - t_j^n \le \delta}} |\Lambda_{t_i^n}^n - \Lambda_{t_j^n}^n|, \quad \delta > 0, $$

we get

$$ w(\Lambda^n, \delta) \le w(G_0, \delta) + \|N^n\|_{\infty} \left( \sup_{0 \le \eta \le \delta} \int_0^{T-\eta} |K(s+\eta) - K(s)| ds + \int_0^{\delta} K(s) ds \right). $$

Since $g_0, K \in L^1([0, T], \mathbb{R}_+)$, the dominated convergence theorem shows that $w(G_0, \delta)$ and $\int_0^\delta K(s) ds$ go to 0 as $\delta \to 0^+$, and Lemma 5.9 proves that it is also the case for

$$ \sup_{0 \le \eta \le \delta} \int_0^{T-\eta} |K(s+\eta) - K(s)| ds. $$

Using the bounds (5.6)-(5.9), we get

$$ \sup_{n \ge 1} \mathbb{E} [\|N^n\|_{\infty}] < +\infty, $$

so that finally

$$ \lim_{\delta \to 0^+} \limsup_{n \to \infty} \mathbb{E}[w(\Lambda^n, \delta)] = 0. \quad (5.10) $$

Thus, $(\Lambda^n)_{n \ge 1}$ is C-tight.

Finally, we prove that $(N^n)_{n \ge 1}$ is tight. We use Aldous' criterion on $(Z^n)_{n \ge 1}$, see Jacod and Shiryaev (2013, Theorem VI.4.5). Let $n \ge 1$, $\delta > 0$, and $0 \le S_1^n \le S_2^n \le T$ two $(\mathcal{F}_t^n)_{t \le T}$-stopping times, such that $S_2^n - S_1^n \le \delta$ almost surely. Then, using the martingality of $[Z^n] - ((k_0^n)^2 + (1-k_0^n)^2) \Lambda^n$ from Lemma 5.2, we have

$$ \begin{aligned} \mathbb{E} \left[ \left( Z_{S_2^n}^n - Z_{S_1^n}^n \right)^2 \right] &= \mathbb{E} \left[ [Z^n]_{S_2^n} - [Z^n]_{S_1^n} \right] \\ &= \left( (k_0^n)^2 + (1-k_0^n)^2 \right) \mathbb{E} \left[ \Lambda_{S_2^n}^n - \Lambda_{S_1^n}^n \right] \\ &\le \left( (k_0^n)^2 + (1-k_0^n)^2 \right) \mathbb{E} [w(\Lambda^n, \delta)]. \end{aligned} $$

Therefore, defining the sets

$$ \mathcal{T}_{\delta}^{n} := \left\{ S_1, S_2 \left( \mathcal{F}_{t}^{n} \right)_{t \le T} \text{stopping times} : 0 \le S_1 \le S_2 \le \min(T, S_1 + \delta) \right\}, \quad \delta > 0, \ n \ge 1, $$

we have

$$ \limsup_{n \ge 1} \sup_{(S_1^n, S_2^n) \in \mathcal{T}_{\delta}^n} \mathbb{P} \left( \left| Z_{S_2^n}^n - Z_{S_1^n}^n \right| > \varepsilon \right) \le \limsup_{n \ge 1} \frac{(k_0^n)^2 + (1-k_0^n)^2}{\varepsilon^2} \mathbb{E}[w(\Lambda^n, \delta)], \quad \delta > 0, \ \varepsilon > 0, $$

which vanishes in the limit $\delta \to 0^+$, from (5.10). Thus, $(Z^n)_{n \ge 1}$ is tight for the Skorokhod topology. We conclude by using the fact that $(\Lambda^n)_{n \ge 1}$ is C-tight, so that $(N^n)_{n \ge 1} = (Z^n + \Lambda^n)_{n \ge 1}$ is tight, thanks to Jacod and Shiryaev (2013, Corollary VI.3.33). $\square$
From Lemma 5.7, $(\Lambda^n, N^n)_{n \ge 1}$ is tight. By Prokhorov's theorem, it admits a weakly convergent subsequence. The following lemma characterizes the accumulation points.

**Lemma 5.8.** Let $(\Lambda, N)$ be a weak accumulation point of $(\Lambda^n, N^n)_{n \ge 1}$. We have the following:

(i) $\Lambda$ is a continuous, non-decreasing and nonnegative process, starting from 0.

(ii) $N$ is a piecewise constant càdlàg process with unit jumps, starting from 0.

(iii) $\Lambda = \int_0^t \lambda_s ds$ with $\lambda$ satisfying

$$ \lambda_t = g_0(t) + \int_0^t K(t-s) dN_s, \quad t \le T. $$

(iv) $N - \Lambda$ is a martingale with respect to the filtration generated by $N$.

*Proof.* We first prove (i), (ii) by standard arguments.

The first point comes directly from the definition of C-tightness and the fact that each $\Lambda^n$ is non-decreasing, nonnegative, and starts from 0.

For (ii), $N_0 = 0$ almost surely as it is also the case for each $N^n$. Furthermore, we know that $\mathbb{P}(N_t^n \in \mathbb{N}, t \le T, n \ge 1) = 1$. Since $\{f \in D([0, T]): f(t) \in \mathbb{N}, t \le T\}$ is a closed set for the Skorokhod topology. We deduce by Portmanteau's theorem that $N$ takes its values in $\mathbb{N}$, almost surely. It remains to prove that its jumps can only be equal to 1. From Lemma 5.6, we have

$$ \mathbb{P} \left( \sup_{0 \le t \le T} \Delta N_t^n > 1 \right) \xrightarrow{n \to \infty} 0. $$

Denoting by $\mathbb{P}'_n$ the law of $N^n$, this can be rewritten as

$$ \mathbb{P}'_n (S_{\Delta}^{-1}((1, +\infty))) \xrightarrow{n \to \infty} 0, $$

where $S_\Delta: f \mapsto \sup_{0 \le t \le T} |\Delta f(t)|$ is shown to be continuous for the Skorokhod topology in Lemma 5.11, so that $S_\Delta^{-1}((1, +\infty))$ is an open set, and by applying Portmanteau's theorem, we obtain

$$ \Delta N_t \le 1, \quad t \le T, \quad \text{a.s.} $$

which concludes the proof of (ii).

For the proof of (iii), we make use of Skorokhod's representation theorem. Since we have a subsequence $(\Lambda^{n_k}, N^{n_k})_{k \ge 0}$ converging weakly to $(\Lambda, N)$, there exists càdlàg processes $(\tilde{\Lambda}^k, \tilde{N}^k)_{k \ge 0}$ and $(\tilde{\Lambda}, \tilde{N})$ all defined on the same probability space $(\tilde{\Omega}, \tilde{\mathcal{F}}, \tilde{\mathbb{P}})$, such that:

* For any $k \ge 0$, $(\tilde{\Lambda}^k, \tilde{N}^k) \sim (\Lambda^{n_k}, N^{n_k})$, and $(\tilde{\Lambda}, \tilde{N}) \sim (\Lambda, N)$.

* For any $\omega \in \tilde{\Omega}$, $(\tilde{\Lambda}^k(\omega), \tilde{N}^k(\omega)) \xrightarrow{k \to \infty} (\tilde{\Lambda}(\omega), \tilde{N}(\omega))$ in the Skorokhod topology.

Using the same reasoning as above, $\tilde{\Lambda}$ is almost surely continuous. We use the notation $G_0^n(t) := \int_0^{\lfloor n \frac{t}{T} \rfloor \frac{T}{t}} g_0(s) ds$, and we define the following subsets of $(D([0, T])^2)$,

$$ \mathcal{V}^n := \left\{ (f,g) \in (D([0,T])^2 : f(t) = G_0^n(t) + \int_{[0,t]} g(t-s) K^n(ds), t \le T \right\}, \quad n \ge 1, $$

$$ \mathcal{V} := \left\{ (f,g) \in (D([0,T])^2 : f(t) = G_0(t) + \int_0^t K(t-s) g(s) ds, t \le T \right\}. $$

From the Volterra equation of Lemma 5.3, we know that

$$ \mathbb{P}((\Lambda^n, N^n) \in \mathcal{V}^n) = 1, \quad n \ge 1. $$
By equality in law of the couples,

$$ \tilde{P}((\tilde{\Lambda}^k, \tilde{N}^k) \in \mathcal{V}^{n_k}) = 1, \quad k \ge 0, $$

leading to

$$ \tilde{P}((\tilde{\Lambda}^k, \tilde{N}^k) \in \mathcal{V}^{n_k}, k \ge 0) = 1 $$

by countability. Thus, the set

$$ \mathcal{A} := \{ \omega \in \tilde{\Omega} : (\tilde{\Lambda}^k(\omega), \tilde{N}^k(\omega)) \in \mathcal{V}^{n_k}, k \ge 0 \} \cap \{ \tilde{\Lambda} \text{ continuous} \}, $$

has unit measure, and for any $\omega \in \mathcal{A}$, we can apply Lemma 5.10 with $(f_k)_{k \ge 0} := (\tilde{N}^k(\omega))_{k \ge 0}$ and $(F_k)_{k \ge 0} := (\tilde{\Lambda}^k(\omega) - G_0^{n_k})_{k \ge 0}$. This leads to $\tilde{P}((\tilde{\Lambda}, \tilde{N}) \in \mathcal{V}) = 1$, and therefore,

$$ \Lambda_t = G_0(t) + \int_0^t K(t-s) N_s ds, \quad t \le T, \quad \text{a.s.} $$

We conclude by using the fact that $N$ has finite variation on $[0, T]$ since it is non-decreasing, which allows us to rewrite the previous equation as

$$ \Lambda_t = \int_0^t \left( g_0(s) + \int_0^s K(s-r) dN_r \right) ds, \quad t \le T, \quad \text{a.s.,} $$

identifying $\lambda$ explicitly.

Finally, for (iv), we have the bound

$$ \sup_{n \ge 1} \sup_{0 \le t \le T} \mathbb{E} \left[ (N_t^n - \Lambda_t^n)^2 \right] \le C_Z $$

from (5.9), which shows that $(N_t^n - \Lambda_t^n)_{t \le T, n \ge 1}$ is uniformly integrable. For each $n \ge 1$, $N^n - \Lambda^n$ is a martingale with respect to the filtration $(\mathcal{F}^n)_{t \le T}$ generated by $N^n$ (see (3.1)), which is also the filtration generated by $(\Lambda^n, N^n - \Lambda^n)$. Thus, applying Whitt (2007, Theorem 5.3) shows that $N - \Lambda$ is a martingale with respect to the filtration generated by $(\Lambda, N - \Lambda)$. This filtration is identical to the one generated by $(\Lambda, N)$, which is also the one generated by $N$ alone since we have shown that $\Lambda = G_0 + K * N$ almost surely. This concludes the proof. $\square$

## 5.4 Additional results

We regroup here some technical results that are used in the previous subsections.

The following lemma is useful for proving the tightness of our processes.

**Lemma 5.9.** Let $K \in L^1([0, T], \mathbb{R})$, we have

$$ \lim_{\delta \to 0^+} \sup_{0 \le \eta \le \delta} \int_0^{T-\eta} |K(s+\eta) - K(s)| ds = 0. $$

*Proof.* Let $\varepsilon > 0$. By density of continuous functions in $L^1([0, T], \mathbb{R})$, there is $f : [0, T] \to \mathbb{R}$ continuous such that $\|K - f\|_{L^1} \le \varepsilon$. We can thus bound

$$ \int_0^{T-\eta} |K(s+\eta) - K(s)| ds \le 2 \|K - f\|_{L^1} + \int_0^{T-\eta} |f(s+\eta) - f(s)| ds \\ \le 2\varepsilon + T w(f, \delta), \quad \eta \le \delta, $$

which gives, by continuity of $f$,

$$ \lim_{\delta \to 0^+} \sup_{0 \le \eta \le \delta} \int_0^{T-\eta} |K(s+\eta) - K(s)| ds \le 2\varepsilon $$

and concludes the proof. $\square$
We also need to ensure some the stability of our Volterra equations in a deterministic setting.

**Lemma 5.10.** Let $(f_n)_{n \ge 1}$ be a sequence of càdlàg functions on $[0, T]$ converging to $f$, càdlàg with bounded variation, in $L^1([0, T])$. We define

$$F_n(t) := \int_{[0,t]} f_n(t-s) K^n(ds), \quad t \le T, \quad n \ge 1,$$

where the measures $(K^n)_{n \ge 1}$ are defined in Lemma 5.3. If $F_n$ converges to $F$ continuous, in $L^1([0, T])$, then

$$F(t) = \int_{0}^{t} K(t-s) f(s) ds, \quad t \le T.$$

*Proof.* For any $n \ge 1$, we can write

$$\begin{aligned} \int_{0}^{T} \left| F(t) - \int_{0}^{t} K(t-s) f(s) ds \right| dt &\le \int_{0}^{T} \left| F(t) - F_n(t) \right| dt + \int_{0}^{T} \left| \int_{[0,t]} (f_n(t-s) - f(t-s)) K^n(ds) \right| dt \\ &\quad + \int_{0}^{T} \left| \int_{[0,t]} f(t-s) K^n(ds) - \int_{0}^{t} f(t-s) K(s) ds \right| dt \\ &=: \mathcal{I}_1^n + \mathcal{I}_2^n + \int_{0}^{T} \left| \mathcal{I}^n(t) - \mathcal{I}(t) \right| dt, \quad n \ge 1. \end{aligned}$$

By $L^1$-convergence, $\mathcal{I}_1^n$ goes to 0 as $n$ goes to infinity. Similarly,

$$\mathcal{I}_2^n \le K^n([0, T]) \int_0^T |f_n(t) - f(t)| dt = \int_0^T K(s) ds \int_0^T |f_n(t) - f(t)| dt, \quad n \ge 1,$$

which vanishes in the limit.

For the remaining term, we use $K^n([0, T]) = \int_0^T K(s) ds$, and for $0 \le t < T$,

$$K^n([0, t]) = \sum_{j=0}^{\lfloor nt/T \rfloor} \int_{t_j^n}^{t_{j+1}^n} K(s) ds = \int_0^T \mathbb{1}_{[\frac{x}{n}, (\lfloor \frac{nt}{T} \rfloor + 1)\frac{x}{n}]}(s) K(s) ds \xrightarrow{n \to \infty} \int_0^t K(s) ds,$$

by the dominated convergence theorem since $K \in L^1([0, T])$. We then rewrite $\mathcal{I}^n(t)$ as

$$\begin{aligned} \mathcal{I}^n(t) &= \int_{[0,t]} f(t-s) K^n(ds) \\ &= f(0) K^n([0,t]) + \int_{[0,t]} \int_0^{t-s} df(u) K^n(ds) \\ &= f(0) K^n([0,t]) + \int_0^t K^n([0,t-u]) df(u), \quad n \ge 0, \end{aligned}$$

which, thanks to the domination $K^n([0, t]) \le K^n([0, T]) = \int_0^T K(s) ds$, converges to

$$f(0) \int_{0}^{t} K(s) ds + \int_{0}^{t} \int_{0}^{t-u} K(s) ds df(u) = \int_{0}^{t} f(t-s) K(s) ds = \mathcal{I}(t).$$

Finally, we can apply the dominated convergence theorem with the bound

$$|\mathcal{I}^n(t) - \mathcal{I}(t)| \le 2 \|f\|_{\infty} \int_0^T K(s) ds, \quad t \le T, \quad n \ge 1,$$

to show that $\int_0^T |\mathcal{I}^n(t) - \mathcal{I}(t)| dt$ goes to 0 as $n \to \infty$. We conclude that

$$\int_{0}^{T} \left| F(t) - \int_{0}^{t} K(t-s) f(s) ds \right| dt = 0,$$
so that, by continuity,

$$F(t) = \int_{0}^{t} K(t-s) f(s) ds, \quad t \leq T. \quad \square$$

Finally, the continuity of the following functional is used to prove that the limit of $(N^n)_{n \geq 1}$ is a counting process.

**Lemma 5.11.** On $D([0, T])$, the function

$$S_{\Delta}: f \mapsto \sup_{0 \leq t \leq T} |\Delta f(t)|, \quad f \in D([0, T]),$$

is continuous for the Skorokhod topology.

*Proof.* Let $f_n \to f$ in the Skorokhod topology. According to Jacod and Shiryaev (2013, Theorem VI.1.14), there exists continuous increasing bijections $(\lambda_n)_{n \geq 0}$ from $[0, T]$ onto itself such that

$$\|f_n \circ \lambda_n - f\|_{\infty} \to 0, \quad \|\lambda_n - \text{Id}\|_{\infty} \to 0,$$

as $n$ goes to infinity. Since

$$\begin{aligned} \left|\left|f_n(\lambda_n(t)) - f_n(\lambda_n(s))\right| - \left|f(t) - f(s)\right|\right| &\leq \left|f_n(\lambda_n(t)) - f(t) + f(s) - f_n(\lambda_n(s))\right| \\ &\leq 2 \|f_n \circ \lambda_n - f\|_{\infty}, \quad s \leq t \leq T, \quad n \geq 0, \end{aligned}$$

we deduce that

$$|S_{\Delta}(f_n \circ \lambda_n) - S_{\Delta}(f)| \leq 2 \|f_n \circ \lambda_n - f\|_{\infty}, \quad n \geq 0,$$

showing that $S_{\Delta}(f_n \circ \lambda_n) \xrightarrow{n \to \infty} S_{\Delta}(f)$. We conclude by noticing that since for each $n \geq 0$ $\lambda_n$ is a continuous bijection from $[0, T]$ onto itself, $S_{\Delta}(f_n \circ \lambda_n) = S_{\Delta}(f_n)$. $\square$

# A Kernels and resolvents

Consider a kernel $K \in L_{\text{loc}}^{1}(\mathbb{R}_{+}, \mathbb{R}_{+})$. There exists a unique kernel $R \in L_{\text{loc}}^{1}(\mathbb{R}_{+}, \mathbb{R}_{+})$, called the *resolvent of the second kind*, or simply the *resolvent*, such that

$$R * K = K * R = R - K.$$

Existence and uniqueness are guaranteed by Gripenberg, Londen, and Staffans (1990, Theorems 2.3.1 and 2.3.5).

Resolvents are particularly useful for solving linear Volterra equations. More precisely, the solution to

$$f = g + K * f$$

is given by

$$f = g + R * g.$$

Below, we provide a table listing some frequently used kernels, along with their integrals $\bar{K}(t) := \int_{0}^{t} K(s) ds$ and the corresponding resolvents. We use $c \geq 0$, $b \in \mathbb{R}$, $\alpha > 0$, and we denote by $E_{\alpha, \beta}$ the Mittag-Leffler function defined by

$$E_{\alpha, \beta}(z) := \sum_{n=0}^{+\infty} \frac{z^n}{\Gamma(\alpha n + \beta)}, \quad z \in \mathbb{C}, \quad \alpha, \beta > 0,$$

and by $\gamma$ the lower incomplete gamma function

$$\gamma(\alpha, x) := \frac{1}{\Gamma(\alpha)} \int_{0}^{x} t^{\alpha-1} e^{-t} dt, \quad x \geq 0, \quad \alpha > 0.$$
<table><thead><tr><th>K(t)</th><th>K̄(t)</th><th>R(t)</th></tr></thead><tbody><tr><td>ce<sup>-bt</sup></td><td>c/b(1 - e<sup>-bt</sup>)</td><td>ce<sup>- (b-c)t</sup></td></tr><tr><td>ct<sup>α-1</sup>/Γ(α)</td><td>ct<sup>α</sup>/Γ(α + 1)</td><td>ct<sup>α-1</sup> E<sub>α,α</sub>(ct<sup>α</sup>)</td></tr><tr><td>ce<sup>-bt</sup>t<sup>α-1</sup>/Γ(α)</td><td>c/b<sup>α</sup>γ(α, bt)</td><td>ce<sup>-bt</sup>t<sup>α-1</sup> E<sub>α,α</sub>(ct<sup>α</sup>)</td></tr><tr><td>ct<sup>α-1</sup>E<sub>α,α</sub>(ct<sup>α</sup>)</td><td>ct<sup>α</sup>E<sub>α,α+1</sub>(ct<sup>α</sup>)</td><td>ct<sup>α-1</sup>E<sub>α,α</sub>(2ct<sup>α</sup>)</td></tr><tr><td>ce<sup>-bt</sup>t<sup>α-1</sup>E<sub>α,α</sub>(ct<sup>α</sup>)</td><td>Σ<sub>n=1</sub><sup>+∞</sup> (c/b<sup>α</sup>)<sup>n</sup>γ(αn, bt)</td><td>ce<sup>-bt</sup>t<sup>α-1</sup>E<sub>α,α</sub>(2ct<sup>α</sup>)</td></tr></tbody></table>

## B The Inverse Gaussian distribution

The Inverse Gaussian distribution, is a probability distribution on $\mathbb{R}_+$ given by the following density

$$f_{\mu,\lambda}(x) = \sqrt{\frac{\lambda}{2\pi x^3}} \exp\left(-\frac{\lambda(x-\mu)^2}{2\mu^2 x}\right), \quad x > 0, \quad \mu, \lambda > 0, \qquad (B.1)$$

If $X \sim IG(\mu, \lambda)$, it admits finite moments of any order, that can be computed using the recursive formula

$$\mathbb{E}[X] = \mu, \quad \text{and} \quad \mathbb{E}[X^n] = (2n-3) \frac{\mu^2}{\lambda} \mathbb{E}[X^{n-1}] + \mu^2 \mathbb{E}[X^{n-2}], \quad n \ge 2. \qquad (B.2)$$

Finally, its characteristic function is given by

$$\mathbb{E}[\exp(wX)] = \exp\left(\frac{\lambda}{\mu}\left(1 - \sqrt{1 - 2\frac{\mu^2}{\lambda}w}\right)\right), \quad w \in \mathbb{C}, \quad \Re(w) \le 0. \qquad (B.3)$$

Inverse Gaussian random variables can be simulated with one standard normal and one uniform variable, using the acceptance-rejection method from Michael, Schucany, and Haas (1976), as follows.

**Algorithm 5** Sampling from the Inverse Gaussian Distribution

1: **Input:** Parameters $\mu > 0$, $\lambda > 0$.

2: **Output:** Sample IG from the Inverse Gaussian distribution.

3: Generate $\xi \sim N(0, 1)$ and compute $Y = \xi^2$.

4: Compute the candidate value $X$:

$$X = \mu + \frac{\mu^2 Y}{2\lambda} - \frac{\mu}{2\lambda} \sqrt{4\mu\lambda Y + \mu^2 Y^2}.$$

5: Generate a uniform random variable: Sample $\eta \sim \text{Uniform}(0, 1)$.

6: Select the output:

7: **if** $\eta \le \frac{\mu}{\mu+X}$ **then**

8: Set the output $IG = X$.

9: **else**

10: Set $IG = \frac{\mu^2}{X}$.

11: **end if**

## References

Eduardo Abi Jaber. Weak existence and uniqueness for affine stochastic Volterra equations with $L^1$-kernels. *Bernoulli*, 27(3):1583–1615, 2021.
Eduardo Abi Jaber. Simulation of square-root processes made simple: applications to the Heston model. *Risk Magazine (Cutting Edge Section)*, June, 2025.

Eduardo Abi Jaber and Elie Attal. Simulating integrated Volterra square-root processes and Volterra Heston models via Inverse Gaussian. *arXiv preprint arXiv:2504.19885*, 2025.

Eduardo Abi Jaber, Martin Larsson, and Sergio Pulido. Affine Volterra processes. *The Annals of Applied Probability*, 29(5), 2019.

Eduardo Abi Jaber, Christa Cuchiero, Martin Larsson, and Sergio Pulido. A weak solution theory for stochastic Volterra equations of convolution type. *The Annals of Applied Probability*, 31(6):2924-2952, 2021.

Emmanuel Bacry, Iacopo Mastromatteo, and Jean-François Muzy. Hawkes processes in finance. *Market Microstructure and Liquidity*, 1(01):1550005, 2015.

Boris Baeumer, Louis Davis, and Ting Wang. A fractional Hawkes process model for earthquake aftershock sequences. *Journal of the Royal Statistical Society Series C: Applied Statistics*, 73(5):1185-1202, 2024.

Adrian Baldwin, Iffat Gheyas, Christos Ioannidis, David Pym, and Julian Williams. Contagion in cyber security attacks. *Journal of the Operational Research Society*, 68(7):780-791, 2017.

Yannick Bessy-Roland, Alexandre Boumezoued, and Caroline Hillairet. Multivariate Hawkes process for cyber insurance. *Annals of Actuarial Science*, 15(1):14-39, 2021.

Alessandro Bondi, Giulia Livieri, and Sergio Pulido. Affine Volterra processes with jumps. *Stochastic Processes and their Applications*, 168:104264, 2024.

Timothy C. Brown and M. Gopalan Nair. A simple proof of the multivariate random time change theorem for point processes. *Journal of Applied Probability*, 25(1):210-214, 1988.

Jing Chen, Alan G. Hawkes, and Enrico Scalas. *A Fractional Hawkes Process*, pages 121-131. Springer International Publishing, Cham, 2021.

Christa Cuchiero and Josef Teichmann. Markovian lifts of positive semidefinite affine Volterra-type processes. *Decisions in Economics and Finance*, 42(2):407-448, 2019.

Angelos Dassios and Hongbiao Zhao. Exact simulation of Hawkes process with exponentially decaying intensity. *Electronic Communications in Probability*, 18:1-13, 2013.

Madalina Deaconu, Antoine Lejay, Laurent Lesage, Jorge A. Meira, Geoffrey Nichil, and Radu State. Hawkes processes framework with a Gamma density as excitation function: application to natural disasters for insurance. *Methodology and Computing in Applied Probability*, 24(4):2509-2537, 2022.

Aline Duarte, Eva Löcherbach, and Guilherme Ost. Stability, convergence to equilibrium and simulation of non-linear Hawkes processes with memory kernels given by the sum of Erlang kernels. *ESAIM: PS*, 23: 770-796, 2019.

Paul Embrechts, Thomas Liniger, and Lu Lin. Multivariate Hawkes processes: an application to financial data. *Journal of Applied Probability*, 48(A):367-378, 2011.

Branda Goncalves and Paul Gresland. Perfect simulation for interacting Hawkes processes with variable length memory. *arXiv preprint arXiv:2209.09143*, 2022.

Gustaf Gripenberg, Stig-Olof Londen, and Olof Staffans. *Volterra Integral and Functional Equations*. Encyclopedia of Mathematics and its Applications. Cambridge University Press, 1990.

Lech A. Grzelak and Nicola F. Zaugg. Lifted Heston Model: efficient Monte Carlo simulation with large time steps. *arXiv preprint arXiv:2510.08805*, 2025.

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

Alan G. Hawkes. Hawkes processes and their applications to finance: a review. *Quantitative Finance*, 18(2): 193-198, 2018.
Alan G. Hawkes and David Oakes. A cluster process representation of a self-exciting process. *Journal of applied probability*, 11(3):493-503, 1974.

Caroline Hillairet, Anthony Réveillac, and Mathieu Rosenbaum. An expansion formula for Hawkes processes and application to cyber-insurance derivatives. *Stochastic Processes and their Applications*, 160:89-119, 2023.

Pierre Hodara and Eva Löcherbach. Hawkes processes with variable length memory and an infinite number of components. *Advances in Applied Probability*, 49(1):84-107, 2017.

Jean Jacod and Albert N. Shiryaev. *Limit theorems for stochastic processes*, volume 288 of *Grundlehren der mathematischen Wissenschaften*. Springer, 2nd edition, 2013.

Matthias Kirchner. Hawkes and INAR($\infty$) processes. *Stochastic Processes and their Applications*, 126(8): 2494-2525, 2016.

Matthias Kirchner. An estimation procedure for the Hawkes process. *Quantitative Finance*, 17(4):571-595, 2017.

Patrick J. Laub, Young Lee, Philip K. Pollett, and Thomas Taimre. Hawkes models and their applications. *Annual Review of Statistics and Its Application*, 12(2025):233-258, 2025.

Peter A. Lewis and Gerald S. Shedler. Simulation of nonhomogeneous Poisson processes by thinning. *Naval Research Logistics Quarterly*, 26(3):403-413, 1979.

John R. Michael, William R. Schucany, and Roy W. Haas. Generating random variates using transformations with multiple roots. *The American Statistician*, 30(2):88-90, 1976.

Jesper Møller and Jakob G. Rasmussen. Perfect simulation of Hawkes processes. *Advances in Applied Probability*, 37(3):629-646, 2005.

Yosihiko Ogata. On Lewis’ simulation method for point processes. *IEEE Transactions on Information Theory*, 27(1):23-31, 1981.

Yosihiko Ogata. Statistical models for earthquake occurrences and residual analysis for point processes. *Journal of the American Statistical association*, 83(401):9-27, 1988.

Philip Protter. *Stochastic Integration and Differential Equations*. Springer Berlin, Heidelberg, 1990.

Mathieu Rosenbaum and Mehdi Tomas. From microscopic price dynamics to multidimensional rough volatility models. *Advances in Applied Probability*, 53(2):425-462, 2021.

Ward Whitt. Proofs of the martingale FCLT. *Probability Surveys*, 4:268-302, 2007.