Patricia Directrice

DR Reynaud-Bouret

Giulia, Julien, Stephano, Marie, etc. Cyrille

Lam Thank

Bien Thoa

Keywords: stochastic simulation, perfect simulation, Kalikow decomposition, Hawkes process, thinning

The goal of this thesis is to construct algorithms which are able to simulate the activity of a neural network. The activity of the neural network can be modeled by the spike train of each neuron, which are represented by a multivariate point processes. Most of the known approaches to simulate point processes encounter difficulties when the underlying network is large.

In this thesis, we propose new algorithms using a new type of Kalikow decomposition. In particular, we present an algorithm to simulate the behavior of one neuron embedded in an infinite neural network without simulating the whole network. We focus on mathematically proving that our algorithm returns the right point processes and on studying its stopping condition. Then, a constructive proof shows that this new decomposition holds for on various point processes.

Finally, we propose algorithms, that can be parallelized and that enables us to simulate a hundred of thousand neurons in a complete interaction graph, on a laptop computer. Most notably, the complexity of this algorithm seems linear with respect to the number of neurons on simulation.

Résumé

L' objectif de cette thèse est de construire des algorithmes capables de simuler l'activité d'un réseau de neurones. L'activité du réseau de neurones peut être modélisée par le train de spikes de chaque neurone, qui sont représentés par un processus ponctuel multivarié. La plupart des approches connues pour simuler des processus ponctuels rencontrent des difficultés lorsque le réseau sous-jacent est de grande taille.

Dans cette thèse, nous proposons de nouveaux algorithmes utilisant un nouveau type de décomposition de Kalikow. En particulier, nous présentons un algorithme permettant de simuler le comportement d'un neurone intégré dans un réseau neuronal infini sans simuler l' ensemble du réseau. Nous nous concentrons sur la preuve mathématique que notre algorithme renvoie les bons processus ponctuels et sur l' étude de sa condition d'arrêt. Ensuite, une preuve constructive montre que cette nouvelle décomposition est valable pour divers processus ponctuels.

Enfin, nous proposons des algorithmes, qui peuvent être parallélisés et qui permettent de simuler une centaine de milliers de neurones dans un graphe d'interaction complet, sur un ordinateur portable. Plus particulièrement, la complexité de cet algorithme semble linéaire par rapport au nombre de neurones à simuler.

Mots clef: simulation stochastique, simulation parfaite, décomposition de Kalikow, processus de Hawkes, algorithme de rejet, simulation à grande échelle, processus ponctuel, processus stochastique. v � Introduction

�.� Mathematical background

The main mathematical object underlying this thesis is the temporal point process. Let us start with some notation and definitions.

�.�.� Point process

A point process 𝑍 in R can be viewed as a (countable) collection of times in R, denoted by (𝑇 𝑛 ) 𝑛∈Z , with the convention that 𝑇 0 ≤ 0 < 𝑇 1 . We will call ℬ(R) the Borel subsets of R. For any bounded measurable set 𝐴 belonging to ℬ(R), we denote 𝑍(𝐴) the cardinal of 𝑍 ∩ 𝐴.

Throughout this thesis, we will only consider the simple point processes. A point process 𝑍 is simple if P(∀𝑡, 𝑍([𝑡]) = 0 or 1) = 1 with the convention that [𝑡] ∶= [𝑡, 𝑡].

The point process 𝑍 can also be defined by its associated discrete random measure 𝑑𝑍 𝑡 ∶= 𝑍(𝑑𝑡), i.e the random measure on (R, ℬ(R)) such that, for any positive measurable function 𝑓, we have ∫ R 𝑓 (𝑡)𝑑𝑍 𝑡 = ∑ 𝑘∈Z 𝑓 (𝑇 𝑘 ). Associated to the point process 𝑍, we also define a counting process 𝑍 𝑡 ∶= 𝑍((0, 𝑡]) if 𝑡 ≥ 0 and 𝑍 𝑡 ∶= -𝑍([𝑡, 0)) otherwise.

Moreover, associated to the process 𝑍, we consider the canonical filtration ℱ 𝑍 𝑡 defined by

ℱ 𝑍 𝑡 = 𝜎(𝑍 ∩ (-∞, 𝑡]).
If the process 𝑍 is such that, for any 𝑡 ℱ 𝑡 ⊃ ℱ 𝑍 𝑡 then ℱ 𝑡 is called a history of the process 𝑍. We will always be in this case in this manuscript.

It is more convenient to study point processes via their stochastic intensity, which is defined in the next section.

�.�.� Stochastic intensity

To facilitate the writing, in this section, we will only consider the point process 𝑍 in R + . The result can be easily extended in a more general framework.

The ℱ-predictable process (𝜙 𝑡 ) 𝑡≥0 is called ℱ-intensity of 𝑍 if 𝑍 𝑡 -∫ 𝑡 0 𝜙 𝑠 𝑑𝑠 is a ℱ 𝑡 martingale. This definition opened a new approach to study point process with abundant tools of martingale theory [2].

Introduction

In practice, to verify that 𝜙 𝑡 is the ℱ intensity of the process 𝑍, we often use the following equivalent definition of the stochastic intensity (Definition D7 of [2]).

Definition. Let 𝑍 𝑡 be a point process adapted to some history ℱ 𝑡 and let 𝜙 𝑡 be a nonnegative ℱ 𝑡 predictable process such that for all 𝑡 ≥ 0 ∫ 𝑡 0 𝜙 𝑠 𝑑𝑠 < ∞ 𝑎.𝑠 If for all nonegative ℱ 𝑡 predictable 𝐶 𝑡 we have

E (∫ ∞ 0 𝐶 𝑠 𝑑𝑍 𝑠 ) = E (∫ ∞ 0 𝐶 𝑠 𝜙 𝑠 𝑑𝑠 ) .
Then we say 𝑍 𝑡 admits the ℱ 𝑡 (predictable) intensity 𝜙 𝑡 .

Most of the models under consideration in this thesis come from neuroscience, that is the object of the next section.

�.� Neuroscience motivation

Neurons or nerve cells are the fundamental units of the brain and nervous system. A neuron receives thousands of synaptic inputs from other neurons, then sums them up (synaptic integration), if this sum is large enough (see Figure 1.2), it produces an action potential also called spike. The action potentials of a given neuron having roughly all the the time the same shape, therefore, we are only interested by the time of their emission. One way to present the neural activity is to record the occurrence of spikes, which we call spike train. That is why we use point processes to model the neural activity and 1. 3 Hawkes process in particular we use Hawkes process, that is presented in the upcoming section. The average number of spikes over a time unit is called firing rate.

�.� Hawkes process

Hawkes process is firstly introduced in the papers [12,13]. It then receives a lot of interest in the last decade and has numerous applications in various fields, to name but a few:

• in seismology [24] • in finance [1] • in genome analysis [START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF] • spread of infection disease [5] .

Throughout this thesis, the integral ∫ 𝑏 𝑎 stands for ∫ 𝑏- 𝑎 or ∫ [𝑎,𝑏) for the reason of predictability, if we do not mention it differently. Most of the time, we will work with multivariate and non-linear Hawkes process. Let us consider I the index set, I = {1, 2, … , 𝑁}. A multivariate Hawkes process (𝑍 1 , … , 𝑍 𝑁 ) is defined by its stochastic intensity for 𝑖 = 1, … , 𝑁, which has the following form in general:

𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜓 𝑖 ( ∑ 𝑗∈I ∫ 𝑡 -∞ ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑍 𝑗 𝑠 ) (1.3.1)
where 𝜓 𝑖 is called the associated intensity function and ℎ 𝑗𝑖 measures the influence of index 𝑗 on index 𝑖.

In neuroscience, the intensity 𝜙 𝑖 represents the firing rate of neuron 𝑖. Here, the firing rate of neuron 𝑖 depends on the spikes of all the neurons before time 𝑡. This formula imitates very well the synaptic integration procedure of a neural network. In addition, for mathematical purposes, we often assume that 𝜓 𝑖 is continuous function (or even Lipschitz for the reason of stability [3,4]). The continuity assumption here also ensures that the 𝜙 𝑖 is predictable [START_REF] Zhu | Nonlinear Hawkes processes[END_REF]. In addition, if we do not state otherwise, we always assume that ℎ 𝑗𝑖 is a 𝐿 1 function, i.e ‖ℎ 𝑗𝑖 ‖ 𝐿 1 = ∫ ∞ 0 |ℎ 𝑗𝑖 (𝑡)|𝑑𝑡 < ∞. We also suppose that 𝜓 𝑖 is locally integrable to avoid explosion [START_REF] Zhu | Nonlinear Hawkes processes[END_REF].

In the following, we will present several type of Hawkes processes. We begin with linear Hawkes process, which is obtained by considering the function 𝜓 𝑖 in Equation 1.3.1 has the following form 𝜓 𝑖 (𝑥) ∶= 𝜇 𝑖 +𝑥.

Then the conditional intensity is given by

𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜇 𝑖 + ∑ 𝑗∈I ∫ 𝑡 -∞ ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑍 𝑗 𝑠 (1.3.2)
where we keep all conditions as stated earlier, 𝜇 𝑖 is a positive number, which represents for the spontaneous rate of the neuron 𝑖 and ℎ 𝑗𝑖 is positive.

In neuroscience, it is reasonable to do not take into account the spikes which appear far from the time of consideration. It leads us to consider the linear Hawkes process with bounded support of interaction 𝐴, with intensity

𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜇 𝑖 + ∑ 𝑗∈I ∫ 𝑡 𝑡-𝐴
ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑍 𝑗 𝑠 (1.3.3) where 𝐴 is a given positive number.

As notice in the Figure 1.2, after a neuron emits a spike, the firing rate of this neuron rapidly decreases to zero. Then it follows by a refractory period, in which the neuron is unable to produce a new spike. This can be modelled by the nonlinear Hawkes process with hard refractory period [START_REF] Raad | Stability for Hawkes processes with inhibition[END_REF] with intensity

𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜓 𝑖 ( ∑ 𝑗∈I ∫ 𝑡 -∞ ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑍 𝑗 𝑠 )
1 𝑡-𝐿 𝑖 (𝑡)>𝛿 (1.3.4) where 𝛿 is a positive number, representing the length of the refractory period and 𝐿 𝑖 (𝑡) = sup{𝑇 ∈ 𝑍 𝑖 such that 𝑇 < 𝑡} the last spiking time before time 𝑡.

In this thesis, to simplify, we only focus on the positive interaction function ℎ 𝑗𝑖 , means we only consider the excitatory effects of neurons. However, our results are possible extended to the case of excitatory and inhibitory also. We left it for future works.

In mathematics, Hawkes processes have been extensively studied in probability theory, statistics and computer science:

1. stationary and stability [3] 2. long time behavior [START_REF] Zhu | Nonlinear Hawkes processes[END_REF] 3. mean field limits [4] 4. simulation [19,20,21] We focus in particular on the simulation of Hawkes processes. Simulation helps us to predict the pattern of point process in the future, for model checking, see the further discussions in [START_REF] Rasmussen | Temporal point processes: the conditional intensity function[END_REF]. In particular, brain simulation is received a lot of attentions recently. There exists several grand projects such as the Human Brain Project in Europe, the Brain Mapping in Japan and the Brain Initiative in the United States. In the next section, we present simulation algorithms of point processes.

�.� Simulation algorithms and Kalikow decomposition �.�.� Ogata's algorithm

We begin with the simulation algorithm proposed by Lewis and Shedler [18] in 1979, which is applied to nonhomogeneous Poisson processes. Inspired by this method, two years later, Ogata [23] introduced a thinning algorithm to simulate point processes having stochastic intensity. Throughout this thesis, we call classical Ogata's algorithm to refer to this method. In the following, we rewrite Proposition 1 of [23].

Say we want to simulate a multivariate point process 𝑍 = (𝑍 𝑖 ) 𝑖 with intensity (𝜙 𝑖 ) 𝑖∈I in the time interval [0, 𝑡 𝑚𝑎𝑥 ] Moreover, we assume that there is no point before time 0, i.e 𝑍 𝑖 0 = 0, ∀𝑖. Then Proposition 1 of [23] can be rewritten as follow.

Proposition. There exist a ℱ 𝑡 predictable process Λ 𝑡 which is defined path-wise such that Let us denote 𝑇 1 , 𝑇 2 , … , 𝑇 Π 𝑡 𝑚𝑎𝑥 be the points of process Π with the intensity Λ 𝑡 . We assign a mark 𝑝 = 1, … , 𝑁 + 1 to point 𝑇 𝑘 if

𝑝-1 ∑ 𝑖=1 𝜙 𝑖 (𝑇 𝑘 |ℱ 𝑇 𝑘 )/Λ 𝑇 𝑘 < 𝑈 𝑘 ≤ 𝑝 ∑ 𝑖=1 𝜙 𝑖 (𝑇 𝑘 |ℱ 𝑇 𝑘 )/Λ 𝑇 𝑘
where 𝑈 𝑘 is an uniform random variable on [0, 1]. Then the points with mark 𝑝 = 1, … , 𝑁 forms a multivariate point process with intensity (𝜙 𝑖 (𝑡|ℱ 𝑡 )) 𝑖 .

This method is quite generic, since it does not need any specific condition on the conditional intensity but it contains two main drawbacks. First, we can not obtain a stationary version of the Hawkes process since we start at time 0 without a point before 0. More importantly, at each point 𝑇 𝑘 , in order to assign a mark, we need to compute the cumulative sum of the intensities. However, this cumulative sum becomes extremely hard to compute when the number of neurons 𝑁 is big. This is one of the main motivations of the thesis, keeping in mind that the neural network of the humain brain is of size 𝑁 = 10 11 . Recently, there exists another approach by Mascart et al. [19], in which, they also aim to simulate such a huge neural network. Their algorithm is a variant of Ogata's algorithm. However, their network needs to be sparse, in the sense that one neuron is only connected to very few other neurons. In our thesis, we can consider a neural network which is complete, i.e each neuron connects to all the other neurons and we can even consider the infinite networks. In the next section, we present how we are able to do so.

�.�.� Kalikow decomposition

This section is devoted to introducing the Kalikow decomposition to improve the classical Ogata thinning algorithm. To begin, we define what a neighborhood is. Intuitively, the neighborhood of (𝑖, 𝑡) is a part of the past that is needed to compute the stochastic intensity 𝜙 𝑖 (𝑡|ℱ 𝑡 ). The rigorous mathematical definition can be found in Chapter 2 and 3. we need to look at either a neuron and a fraction of time before 𝑡 (the green line in Figure 1.3(a)), two neurons and two different fractions of times before 𝑡 (the blue lines in Figure 1.3(b)) or empty set (Figure 1.3(c)), i.e independent to the past.

The whole point of Kalikow decomposition is to pick at random a neighborhood in a family of spacetime neighborhoods (neighborhood family) according to a certain distribution. This can be done thanks to the Kalikow decomposition, which is presented in the following.

We say that a stochastic intensity 𝜙 𝑖 (𝑡|ℱ 𝑡 ) admits a Kalikow decomposition with respect to the neighborhood family V if there is a sequence of function 𝜙 𝑣 𝑘 𝑖 (𝑡|ℱ 𝑡 ) which only depends on the information in the neighborhood 𝑣 𝑘 and the probability distribution on V, 𝜆 𝑖,𝑡 (𝑣) such that 𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜆 𝑖,𝑡 (∅)𝜙 ∅ 𝑖 (𝑡) + ∑ 𝑘 𝜆 𝑖,𝑡 (𝑣 𝑘 )𝜙

𝑣 𝑘 𝑖 (𝑡|ℱ 𝑡 )
with 𝜆 𝑖,𝑡 (∅) + ∑ 𝑘 𝜆 𝑖,𝑡 (𝑣 𝑘 ) = 1.

The Kalikow decomposition can be interpreted as follows. The intensity function 𝜙 𝑖 (𝑡|ℱ 𝑡 ) in some cases might depend on the whole past. Instead of taking the whole information in the past, we draw a random neighborhood 𝑉, that is a finite part of the past according to the distribution 𝜆 𝑖,𝑡 (.). Once the neighborhood is chosen, assume that is is 𝑣 𝑘 , we then compute the intensity function according to 𝜙 𝑣 𝑘 𝑖 (𝑡|ℱ 𝑡 ) which depends only on the points in 𝑣 𝑘 . The information in 𝑣 𝑘 is much less than the whole past, which hopefully helps us to compute the intensity more efficiently.

Kalikow decomposition [16] was largely studied in the context of stochastic processes having long memory [6,10,11,[START_REF] Ost | Sparse space-time models: Concentration inequalities and Lasso[END_REF]. But the decomposition was applied to transition probability rather than to a stochastic intensity in our study. As a result, they consider only discrete-time processes. Despite numerous studies on discrete-time processes, up to our knowledge, there exists only one work on continuous processes by Hodara and Löcherbach [14]. It is worth noting that, in their paper, the intensities are assumed to be bounded by a constant, and 𝜆 𝑖,𝑡 (𝑣 𝑘 ) are not deterministic but depends on the realization of a dominating Poisson process, which in practice need to be simulated at prior to write the Kalikow decomposition.

To distinguish, we call such decomposition, the conditional Kalikow decomposition and the one defined with deterministic 𝜆 by unconditional Kalikow decomposition. A very natural question arises, does this unconditional Kalikow decomposition exists in different examples? If it is the case, how to write it? How we can apply this decomposition in the simulation? Remark 1. It is worth noting that, at this stage the distribution function 𝜆(.) in Kalikow decomposition can be either deterministic or random. We have not proved yet such a deterministic Kalikow decomposition exists in the continuous time process.

�.�.� Perfect simulation

In this section, we investigate several perfect simulation algorithms. By perfect simulation, we mean it can simulate the stationary distribution of a finite number of components in a potentially infinite network 1 Introduction within a given finite space-time window. To explain clearly and coherently, we first examine the algorithms for discrete-time processes then continuous-time processes. This thesis mostly deals with Hawkes processes, which are in particular a process with long memory. Therefore, we will focus to the perfect algorithm developed by Comets et al. [6] in the discrete-time processes and compare it with the famous perfect simulation, named Coupling From the Past (CFTP) of Propp and Wilson [START_REF] Propp | Exact sampling with coupled Markov chains and applications to statistical mechanics[END_REF]. Some other approaches to perfect simulations can be found at [8].

Perfect simulation is a widely discussed topic since the first seminal paper of Propp and Wilson [START_REF] Propp | Exact sampling with coupled Markov chains and applications to statistical mechanics[END_REF].

The goal of this paper is to simulate the invariant measure of Markov chains. Following Propp and Wilson [START_REF] Propp | Exact sampling with coupled Markov chains and applications to statistical mechanics[END_REF], fixed a time 0, for a fixed predeterminded time in the past, we simulate one trajectory per initial state. If the trajectories coalesce, the common value at time 0 is picked under the stationary distribution. If the trajectories do not coalesce, we extend their paths further in the past. Theorem 1 of [START_REF] Propp | Exact sampling with coupled Markov chains and applications to statistical mechanics[END_REF] showed that if we start from a far enough time in the past, the trajectories starting from all possible initial states coalesce.

Different from the Propp and Wilson's method, perfect simulation in the paper of Comets et al. [6] used to simulate the invariant measure of the processes with long memory, which is more general than Markov chains and Markov processes. This method is based on the regenerative construction of processes, which is (as the authors cited) introduced in the paper [9]. In their method, at each instant 𝑛, there is a (random) number 𝑘 𝑛 ≥ 0 such that the distribution of the move 𝑛 + 1 is the same for all histories agreeing the 𝑘 𝑛 preceding instants. This independence from 𝑘 𝑛 remote past yields to the times 𝜏 such that preceding past before 𝜏 are irrelevant for future moves. If there exist 𝜏 > -∞ such that 𝜏 ≤ 𝑛 -𝑘 𝑛 for all 𝑛 ≥ 𝜏, we call such 𝜏 the regeneration time. Moreover in Section 2 of [6], the authors presented explicitly the regenerative construction. This second approach by Comets et al. relies on the regenerative representation of transition probability.

Recently, by using the Kalikow decomposition's notation, Galves et al. provided a more general approach [10,11] with the same spirit of Comets's approach. In short, as discussed in the previous section, Kalikow decomposition allows us to (randomly) determine a space-time neighborhood, which then constitute a clan of ancestors of the interested point. Ancestor of a point 𝑇 is the point that might influence to 𝑇 (see Chapter 3 to more rigorous definition). By applying the Kalikow decomposition to these obtaining points in the clan of ancestor, we obtain the ancestors of ancestor of the interested point. Repeat this procedure until all the points either having empty neighborhood or visited neighborhood, means it only contains the visited points. By doing so, we create a branching processes containing all influent points to the interest point. Denote 𝜏 the point with minimum in time. We see that all history before time 𝜏 is irrelevant to determine the state of the interested point and that this property is similar to [6]. Back to our perfect simulation algorithm, once the clan is completely constructed, we process with the Forward algorithm.

In which, we begin with the smallest point. Then we continue with the point whose neighborhood is completely determined. We repeat this until the state of interested point is determined.

Most the existing papers focus on discrete time processes, in this thesis, we want to provide Perfect simulation for continuous time processes. Therefore, the approach of Comets need some adaptation. Indeed, we will need to decompose conditional intensities rather than transition probabilities. This leads to serious difficulties that usually prevent a more practical application of the Perfect Simulation algorithm.

Again, up to our knowledge, the only work dealing with continuous time counting processes, is the one by Hodara and Löcherbach [15]. Their decomposition is constructed under the assumption that there is a 1.5 Presentation of the thesis dominating Poisson process on each of the nodes, from which the points of the processes under interest can be thinned by rejection sampling (see also [22] for another use of thinning in simulation of counting processes). To prove the existence of a Kalikow decomposition and go back to a more classic discrete time setting, the authors need to freeze the dominating Poisson process, leading to a mixture, in the Kalikow decomposition, that depends on the realization of the dominating Poisson process. Such a mixture is not accessible in practice, and this prevents the use of their Perfect Simulation algorithm for more concrete purposes than mere existence.

There also exists other approaches to do perfect simulation for continuous point processes:

• by Møller and Rasmussen [21]: this method only applies for linear Hawkes processes which can be presented as Poisson cluster processes. Then the length distribution of the clusters are required to perform simulation in their algorithm. However as they discussed, these are not known even in simple cases. Therefore, the authors used the idea of dominated Coupling From The Past ( [17]), they replaced the (unkown) cumulative distribution function by the upper and lower bound of it.

By this method, they create coupled upper and lower process. These processes will converge to the true Hakwes processes.

• by Dassio and Zhao [7]: this method exploit some special properties of exponential Hawkes process.

�.� Presentation of the thesis

In Chapter 2, we introduce a new simulation algorithm in the stationary (time-homogeneous) system using the unconditional Kalikow decomposition. We apply this algorithm to a neuroscience question.

More precisely, we study the behavior of one part of a neural network with potentially an infinite number of neurons. We consider I the set of neurons, which might be infinite (countable). Throughout Chapter 2, we focus on the following conditional intensity

𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜇 𝑖 + ∑ 𝑗∈I 𝑤 𝑖𝑗 (𝑍 𝑗 ([𝑡 -𝐴, 𝑡)) ∧ 𝑀)
where 𝑀 is a given positive number and 𝑤 𝑖𝑗 ∈ [0, 1] is the proportion to the interaction between neuron 𝑖 and 𝑗 such that ∑ 𝑗 𝑤 𝑖𝑗 < 1 for all 𝑖. Clearly, within this formula, the conditional intensity is bounded.

We can apply Ogata's thinning algorithm [23]. However, to do so, we need to simulate the whole system in order to simulate one neuron. In addition, starting from time 0, it does not go backward in time, therefore we can not simulate a Hawkes process in a stationary regime. There exist a work of Moller and Rassmussen [21] that provides a perfect simulation but needs the branching structure of the process to do so. In this chapter, we present another approach that also overcomes this flaw. We define a new type of Kalikow decomposition in continuous time (Definition 2.4.1 of Chapter 2).

Definition. We say that the process admits a Kalikow decomposition with bound 𝑀 and neighborhood family V, if for any neuron 𝑖 ∈ I, for all 𝑣 ∈ V there exists a non negative 𝑀-bounded quantity 𝜙 𝑣 𝑖,𝑡 whose 𝑖 depends on the points in 𝑉 𝑇 , which are not known at this stage. To deal with this problem, we propose a modified version of the Perfect simulation, which we call Backward Forward algorithm (Algorithm 3). We prove that this algorithm ends in finite time almost surely in Proposition 1. Moreover, this algorithm also returns the right intensity process, which is proved in Proposition 2.

In Chapter 3, we show that unconditional Kalikow decomposition exists in various types of Hawkes processes. In this chapter, we also define for the first time an unconditional Kalikow decomposition in general case (Definition 3.3.3).

Definition. We say a process 𝑍 𝑖 for some 𝑖 ∈ I admits the Kalikow decomposition with respect to a neighborhood family V if for all 𝑡, for any 𝑣 ∈ V there exists a cylindrical function 𝜙 𝑣 𝑖 (𝑡|.) taking values in R + , whose value only depends on the points appearing in the neighborhood 𝑣, and a probability density function 𝜆 𝑖 (.) such that

𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜆 𝑖 (∅)𝜙 ∅ 𝑖 + ∑ 𝑣∈V,𝑣≠∅ 𝜆 𝑖 (𝑣)𝜙 𝑣 𝑖 (𝑡|ℱ 𝑡 ) (1.5.2) with 𝜆 𝑖 (∅) + ∑ 𝑣∈V,𝑣≠∅ 𝜆 𝑖 (𝑣) = 1.
Remark 3. Note that, in this definition of Kalikow decomposition, there is no restriction to 𝜙 𝑣 and the probability distribution 𝜆 is deterministic.

Next, very naturally, we show how to obtain the unconditional Kalikow decomposition for a various processes. This is a very constructive method and it also proved the existence of unconditional Kalikow decomposition, see Proposition 3 of Chapter 3.

On the other hand, for the nonlinear Hawkes process, in particular, we also propose another method using Taylor expansion in Proposition 5. We extend the definition of the random neighborhood in Chapter 2 and introduce the modified Perfect simulation with the general random neighborhood, see Backward-Foward algorithm in Section 3.5 of Chapter 3. Similar to the previous chapter, Proposition 7 gives the stopping condition for the Backward-Foward algorithm and Proposition 6 shows that it provides a right intensity process. In this chapter, we also study the complexity of the algorithm in terms of the number of simulated points (see Proposition 8 of Chapter 3). Based on these discussions, we show how to optimize the Kalikow decomposition in the last section of this chapter.

Presentation of the thesis

In Chapter 4, based on the theory developed in Chapter 3, we present several new simulation algorithms for point processes. We focus on two different hypotheses: either the conditional intensity is bounded by a deterministic number (Definition 7) or it is bounded by a predictable function (Definition 8).

In the first case, we rewrite Ogata's algorithm in our settings (see Algorithm 4). This version is sequential.

Notably, we construct a parallelized version of Ogata's algorithm (see Algorithm 5). Next, by adding Kalikow decomposition, we propose two new algorithms, one is purely sequential (Algorithm 6), the other is partially parallelized (Algorithm 7).

In the second case, when the intensity is bounded by a predictable function, we construct two sequential algorithms: one we call Ogata's algorithm (Algorithm 8), the second one is KalikowOgata algorithm (Algorithm 9).

In the deterministic bounded case, we design an algorithm for a particular situation, namely the Hawkes process with a refractory period (Section 4.3.5), whereas in the second case, we focus on proving mathematical the correctness of the algorithms (Section 4.4.3). In the last section, we provide several statistical tests and execution times of algorithms in Section 4.5. Notably, the Kalikow Ogata algorithm is showed to be significantly better than the classical Ogata algorithm and to have linear complexity even for a complete graph on simulation.

We left chapters 2, 3, 4 under their article form except for a few modifications in order to preserve consistency between the chapters. In particular, some definitions may be redundant with this introduction.

�.�.� Main contributions

First and foremost, my work is mostly inspired by practical questions. Therefore, I rather focus on the possibility of implementing the algorithms than the theoretical question about the existence of invariant measures like in [6,10,11,14]. For example, in the paper [14] (which is closest to our work and is also the only work in continuous time processes), the authors need to simulate the whole underlying network to obtain Kalikow decomposition and then do a perfect simulation. From a computational point of view, this approach is unfeasible (both for finite and infinite network). Our approach only simulates parts of the dominating Poisson processes which are needed (in the bounded case). More importantly, as we discussed earlier in Section 1.4.2, going from conditional Kalikow decomposition to the unconditional one is crucial in practice. In addition, many numerical results are presented in Chapter 2 and Chapter 4 together with numerous discussions.

Besides, this thesis also has some theoretical contributions. The previous works on Kalikow decomposition [6,9,10] mostly relied on the continuity assumption, which is not necessary in the thesis. The Linear Hawkes process which is indeed not satisfied the continuity assumption with the "nested" neighborhood, therefore the classical Kalikow decomposition can not be obtained. However, by considering the Linear Hawkes process with specific choice of neighborhood, we obtain Kalikow decomposition, see Chapter 3. Lastly, in Chapter 4 of the thesis, we also proposed an algorithm to simulate a point process with unbounded intensity. We extended the idea of Ogata [23] with Kalikow decomposition to introduce numerous new algorithms, see Chapter 4.

Introduction

�.�.� State of the works

Chapter 2 corresponds to the paper [START_REF] Phi | Event-scheduling algorithms with Kalikow decomposition for simulating potentially infinite neuronal networks[END_REF], in collaboration with Alexandre Muzy and Patricia Reynaud-Bouret. This paper is published in Springer Nature Computer Science.

Chapter 3 corresponds to the paper [START_REF] Phi | Kalikow decomposition for counting processes with stochastic intensity[END_REF]. This paper is submitted to Advances in Applied Probability, which is now in major revision.

�.� abstract

Event-scheduling algorithms can compute in continuous time the next occurrence of points (as events) of a counting process based on their current conditional intensity. In particular event-scheduling algorithms can be adapted to perform the simulation of finite neuronal networks activity. These algorithms are based on Ogata's thinning strategy [17], which always needs to simulate the whole network to access the behaviour of one particular neuron of the network. On the other hand, for discrete time models, theoretical algorithms based on Kalikow decomposition can pick at random influencing neurons and perform a perfect simulation (meaning without approximations) of the behaviour of one given neuron embedded in an infinite network, at every time step. These algorithms are currently not computationally tractable in continuous time. To solve this problem, an event-scheduling algorithm with Kalikow decomposition is proposed here for the sequential simulation of point processes neuronal models satisfying this decomposition. This new algorithm is applied to infinite neuronal networks whose finite time simulation is a prerequisite to realistic brain modeling.

�.� Introduction

Point processes in time are stochastic objects that model efficiently event occurrences with a huge variety of applications: time of deaths, earthquake occurrences, gene positions on DNA strand, etc. [1,20,22]).

Most of the time, point processes are multivariate [6] in the sense that either several processes are considered at the same time, or in the sense that one process regroups together all the events of the different processes and marks them by their type. A typical example consists in considering either two processes, one counting the wedding events of a given person and one counting the children birth dates of the same person or only one marked process which regroups all the possible dates of birth or weddings independently and adds one mark per point, here wedding or birth.

Consider now a network of neurons each of them emitting action potentials (spikes). These spike trains can be modeled by a multivariate point process with a potentially infinite number of marks, each mark representing one given neuron. The main difference between classical models of multivariate point processes and the ones considered in particular for neuronal networks is the size of the network. A human brain consists in about 10 11 neurons whereas a cockroach contains already about 10 6 neurons. Therefore the simulation of the whole network is either impossible or a very difficult and computationally intensive task for which particular tricks depending on the shape of the network or the point processes have to be used [5,13,19].

Another point of view, which is the one considered here, is to simulate, not the whole network, but the events of one particular node or neuron, embedded in and interacting with the whole network. In this sense, one might consider an infinite network. This is the mathematical point of view considered in a series of papers [8,9,18] and based on Kalikow decomposition [11] coupled with perfect simulation theoretical algorithms [4,7]. However these works are suitable in discrete time and only provide a way to decide at each time step if the neuron is spiking or not. They cannot operate in continuous time, i.e. they cannot directly predict the next event (or spike). Up to our knowledge, there exists only one attempt of using such decomposition in continuous time [10], but the corresponding simulation algorithm is purely theoretical in the sense that the corresponding conditional Kalikow decomposition should exist given the whole infinite realization of a multivariate Poisson process, with an infinite number of marks, quantity which is impossible to simulate in practice.

The aim of the present work is to present an algorithm which • can operate in continuous time in the sense that it can predict the occurrence of the next event. In this sense, it is an event-scheduling algorithm;

• can simulate the behavior of one particular neuron embedded in a potentially infinite network without having to simulate the whole network;

• is based on an unconditional Kalikow decomposition and in this sense, can only work for point processes with this decomposition.

In Section 2.3, we specify the links between event-scheduling algorithms and the classical point process theory. In Section 2.4, we give the Kalikow decomposition. In Section 2.5, we present the backwardforward perfect simulation algorithm and prove why it almost surely ends under certain conditions. In Section 2.6, we provide simulation results and a conclusion is given in Section 2.7.

�.� Event-scheduling simulation of point processes

On the one hand, simulation algorithms of multivariate point processes [17] are quite well known in the statistical community but as far as we know quite confidential in the simulation (computer scientist) community. On the other hand, event-scheduling simulation first appeared in the mid-1960s [21] and was formalized as discrete event systems in the mid-1970s [24] to interpret very general simulation algorithms scheduling "next events". A basic event-scheduling algorithm "jumps" from one event occurring at a time stamp 𝑡 ∈ R + 0 to a next event occurring at a next time stamp 𝑡 ′ ∈ R + 0 , with 𝑡 ′ ≥ 𝑡. In a discrete event system, the state of the system is considered as changing at times 𝑡, 𝑡 ′ and conversely unchanging in between [23]. In [13], we have written the main equivalence between the point processes simulation algorithms and the discrete event simulation set-up, which led us to a significant improvement in terms of computational time when huge but finite networks are into play. Usual event-scheduling simulation algorithms have been developed considering independently the components (nodes) of a system. Our approach considers new algorithms for activity tracking simulation [16]. The event activity is tracked from active nodes to children (influencees).

Here we just recall the main ingredients that are useful for the sequel.

To define point processes, we need a filtration or history (ℱ 𝑡 ) 𝑡≥0 . Most of the time, and this will be the case here, this filtration (or history) is restricted to the internal history of the multivariate process (ℱ 𝑖𝑛𝑡 𝑡 ) 𝑡≥0 , which means that at time 𝑡-, i.e. just before time 𝑡, we only have access to the events that have already occurred in the past strictly before time 𝑡, in the whole network. The conditional intensity, 𝜙 𝑖 (𝑡|ℱ 𝑖𝑛𝑡 𝑡-), of the point process, representing neuron 𝑖 gives the instantaneous firing rate, that is the frequency of spikes, given the past contained in ℱ 𝑖𝑛𝑡 𝑡-. Let us just mention two very famous examples.

If 𝜙 𝑖 (𝑡|ℱ 𝑖𝑛𝑡 𝑡-) is a deterministic constant, say 𝑀, then the spikes of neuron 𝑖 form a homogeneous Poisson process with intensity 𝑀. The occurrence of spikes are completely independent from what occurs elsewhere in the network and from the previous occurrences of spikes of neuron 𝑖.

If we denote by I the set of neurons, we can also envision the following form for the conditional intensity:

𝜙 𝑖 (𝑡|ℱ 𝑖𝑛𝑡 𝑡-) = 𝜈 𝑖 + ∑ 𝑗∈ I 𝑤 𝑖𝑗 ( Nb 𝑗 [𝑡-𝐴,𝑡) ∧ 𝑀). (2.3.1)
This is a particular case of generalized Hawkes processes [3]. More precisely 𝜈 𝑖 is the spontaneous rate (assumed to be less than the deterministic upper bound 𝑀 > 1) of neuron 𝑖. Then every neuron in the network can excite neuron 𝑖: more precisely, one counts the number of spikes that have been produced by neuron 𝑗 just before 𝑡, in a window of length 𝐴, this is Nb 𝑗 [𝑡-𝐴,𝑡) ; we clip it by the upper bound 𝑀 and modulate its contribution to the intensity by the positive synaptic weight between neuron 𝑖 and neuron 𝑗, 𝑤 𝑖𝑗 . For instance, if there is only one spike in the whole network just before time 𝑡, and if this happens on neuron 𝑗, then the intensity for neuron 𝑖 becomes 𝜈 𝑖 + 𝑤 𝑖𝑗 . The sum over all neurons 𝑗 mimics the synaptic integration that takes place at neuron 𝑖. As a renormalization constraint, we assume that sup 𝑖∈ I ∑ 𝑗∈ I 𝑤 𝑖𝑗 < 1. This ensures in particular that such a process has always a conditional intensity bounded by 2𝑀.

Hence, starting for instance at time 𝑡, and given the whole past, one can compute the next event in the network by computing for each node of the network the next event in absence of other spike apparition.

To do so, remark that in absence of other spike apparition, the quantity 𝜙 𝑖 (𝑠|ℱ 𝑖𝑛𝑡 𝑠-) for 𝑠 > 𝑡 becomes for instance in the previous example

𝜙 𝑎𝑏𝑠 𝑖 (𝑠, 𝑡) = 𝜈 𝑖 + ∑ 𝑗∈ I 𝑤 𝑖𝑗 ( Nb 𝑗 [𝑠-𝐴,𝑡) ∧ 𝑀),
meaning that we do not count the spikes that may occur after 𝑡 but before 𝑠. This can be generalized to more general point processes. The main simulation loop is presented in Algorithm 1

𝑖 after 𝑡. It is a discrete event approach with the state corresponding to the function 𝜙 𝑎𝑏𝑠 𝑖 (., 𝑡). Ogata [17], inspired by Lewis' algorithm [12], added a thinning (also called rejection) step on top of this procedure because the integral ∫ 𝑇 (1) 𝑖 𝑡 𝜙 𝑎𝑏𝑠 𝑖 (𝑠, 𝑡)𝑑𝑠 can be very difficult to compute. To do so (and simplifying a bit), assume that 𝜙 𝑖 (𝑡|ℱ 𝑖𝑛𝑡 𝑡-) is upper bounded by a deterministic constant 𝑀. This means that the point process has always less points than a homogeneous Poisson process with intensity 𝑀. Therefore

Steps 5-6 of Algorithm 1 can be replaced by the generation of an exponential of parameter 𝑀, 𝐸 ′ 𝑖 and deciding whether we accept or reject the point with probability 𝜙 𝑎𝑏𝑠 𝑖 (𝑡 + 𝐸 ′ 𝑖 , 𝑡)/𝑀. There are a lot of variants of this procedure: Ogata's original one uses actually the fact that the minimum of exponential variables, is still an exponential variable. Therefore one can propose a next point for the whole system, then accept it for the whole system and then decide on which neuron of the network the event is actually appearing. More details on the multivariate Ogata's algorithm can be found in [13].

As we see here, Ogata's algorithm is very general but clearly needs to simulate the whole system to simulate only one neuron. Moreover starting at time 𝑡 0 , it does not go backward and therefore cannot simulate a Hawkes process in stationary regime. There has been specific algorithms based on clusters representation that aim at perfectly simulate particular univariate Hawkes processes [15]. The algorithm that we propose here, will also overcome this flaw.

�.� Kalikow decomposition

Kalikow decomposition relies on the concept of neighborhood, denoted by 𝑣, which are picked at random and which gives the portion of time and neuron subsets that we need to look at, to move forward. Typically, for a positive constant 𝐴, such a 𝑣 can be:

• {(𝑖, [-𝐴, 0))}, meaning we are interested only by the spikes of neuron 𝑖 in the window [-𝐴, 0);

• {(𝑖, [-2𝐴, 0)), (𝑗, [-2𝐴, -𝐴))}, that is, we need the spikes of neuron 𝑖 in the window [-2𝐴, 0) and the spikes of neuron 𝑗 in the window [-2𝐴, -𝐴);

• the emptyset ∅, meaning that we do not need to look at anything to pursue.

We need to also define 𝑙(𝑣) the total time length of the neighborhood 𝑣 whatever the neuron is. For instance, in the first case, we find 𝑙(𝑣) = 𝐴, in the second 𝑙(𝑣) = 3𝐴 and in the third 𝑙(𝑣) = 0.

We are only interested by stationary processes, for which the conditional intensity, 𝜙 𝑖 (𝑡 | ℱ 𝑖𝑛𝑡 𝑡 -), only depends on the intrinsic distance between the previous points and the time 𝑡 and not on the precise value of 𝑡 per se. In this sense the rule to compute the intensity may be only defined at time 0 and then shifted by 𝑡 to have the conditional intensity at time 𝑡. In the same way, the timeline of a neighborhood 𝑣 is defined as a subset of R * -so that information contained in the neighborhood is included in ℱ 𝑖𝑛𝑡 0-, and 𝑣 can be shifted (meaning its timeline is shifted) at position 𝑡 if need be. We assume that I the set of neurons is countable and that we have a countable set of possibilities for the neighborhoods 𝒱.

Then, we say that the process admits a Kalikow decomposition with bound 𝑀 and neighborhood family 𝒱, if for any neuron 𝑖 ∈ I, for all 𝑣 ∈ 𝒱 there exists a non negative 𝑀-bounded quantity 𝜙 𝑣 𝑖 , which is The interpretation of (2.4.1) is tricky and is not as straightforward as in the discrete case (see [18]). The best way to understand it is to give the theoretical algorithm for simulating the next event on neuron 𝑖 after time 𝑡 (cf. Algorithm 2).

Algorithm 2 Kalikow theoretical simulation algorithm

▷ With [𝑡 0 , 𝑡 1 ] the interval of simulation for neuron 𝑖 1: Initialize the family of points P = ∅ ▷ NB: since we are only interested by points on neuron 𝑖, 𝑗 𝑇 = 𝑖 is a useless mark here. Draw an exponential variable 𝐸 with parameter 𝑀, and compute 𝑇 = 𝑡 + 𝐸.

5:

Pick a random neighborhood according to the distribution 𝜆 𝑖 (.) given in the Kalikow decomposition and shift the neighborhood at time 𝑇: this is 𝑉 𝑇 . 𝑖 depends on the points in 𝑉 𝑇 , that are not known at this stage. That is why the efficient algorithm that we propose in the next section goes backward in time before moving forward.

Networks

�.� Backward Forward algorithm

Let us now describe the complete Backward Forward algorithm (cf. Algorithm 3). Note that to do so, the set of points P is not reduced, as in the two previous algorithms, to the set of points that we want to simulate but this contains all the points that need to be simulated to perform the task. Update 𝑉 𝑇 by drawing 𝑉 𝑇 according to 𝜆 𝑗 𝑇 shifted at time 𝑇. ▷ It is possible that the algorithm generated points before 𝑡 0 and they have to be removed Init.

Init.

Marks

Backward Forward At the difference with Algorithm 2, the main point is that in the backward part we pick at random all the points that may influence the thinning step. The fact that this loop ends comes from the following Proposition.

Networks

Proposition 1. If

sup 𝑖∈ I ∑ 𝑣∈𝒱 𝜆 𝑖 (𝑣)𝑙(𝑣)𝑀 < 1. (2.5.1)
then the backward part of Algorithm 3 ends almost surely in finite time.

The proof is postponed to Appendix 2.9. It is based on branching process arguments. Basically if in

Steps 8-16, we produce in average less than one point, either because we picked the empty set in 𝑉 𝑇 or because the simulation of the Poisson process ended up with a small amount of points, eventually none, then the loop ends almost surely because there is an extinction of the corresponding branching process.

In the backward part, one of the most delicate part consists in being sure that we add new points only if we have not visited this portion of time/neurons before (see ). If we do not make this verification, we may not have the same past depending on the neuron we are looking at and the procedure would not simulate the process we want.

In the forward part, because the backward algorithm stopped just before, we are first sure to have assess all 𝑉 𝑇 's. Since 𝜙

𝑉 𝑡 𝑗 (𝑡) is ℱ 𝑖𝑛𝑡 𝑡-measurable, for all 𝑡, 𝜙 𝑉 𝑇 𝑗 𝑇
(𝑇 ) only depends on the points in P with mark

𝑋 𝑇 = 1 inside 𝑉 𝑇 .
The problem in Algorithm 2, phrased differently, is that we do not know the marks 𝑋 𝑇 of the previous points when we have to compute 𝜙 𝑉 𝑇 𝑗 𝑇 (𝑇 ). But in the forward part of Algorithm 3, we are sure that the most backward point for which the thinning (𝑋 𝑇 =n.a.) has not taken place, satisfies

• either 𝑉 𝑇 = ∅

• or 𝑉 𝑇 ≠ ∅ but either there are no simulated points in the corresponding 𝑉 𝑇 or the points there come from previous rounds of the loop (Step 5). Therefore their marks 𝑋 𝑇 have been assigned.

Therefore, with the Backward Forward algorithm, and at the difference to Algorithm 2, we take the points in an order for which we are sure that we know the previous needed marks. the points that may influence the acceptation/rejection of point 𝑇 𝑛𝑒𝑥𝑡 . Notice that whereas usual activity tracking algorithms for point processes [13] automatically detect the active children (influencees), activity tracking in the backward steps detect the parents (influencers). The forward steps finally select the points.

�.� Illustration

To illustrate in practice the algorithm, we have simulated a Hawkes process as given in (2.3.1). Indeed such a process has a Kalikow decomposition (2.4.1) with bound 𝑀 and neighborhood family 𝒱 constituted of the 𝑣's of the form 𝑣 = {(𝑗, [-𝐴, 0))} for some neuron 𝑗 in I. To do that, we need the following choices:

𝜆 𝑖 (∅) = 1 -∑ 𝑗∈ I 𝑤 𝑖𝑗 and 𝜙 ∅ 𝑖 = 𝜈 𝑖 𝜆 𝑖 (∅)
and for 𝑣 of the form 𝑣 = {(𝑗, [-𝐴, 0))} for some neuron 𝑗 in I,

𝜆 𝑖 (𝑣) = 𝑤 𝑖𝑗 and 𝜙 𝑣 𝑖 = Nb 𝑗 [-𝐴,0) ∧ 𝑀.
We have taken I = Z2 and the 𝑤 𝑖𝑗 proportional to a discretized centred symmetric bivariate Gaussian distribution of standard deviation 𝜎. More precisely, once 𝜆 𝑖 (∅) = 𝜆 ∅ fixed, picking according to 𝜆 𝑖 consists in

• choosing whether 𝑉 is empty or not with probability

𝜆 ∅ • if 𝑉 ≠ ∅, choosing 𝑉 = {(𝑗, [-𝐴, 0))} with 𝑗 -𝑖 = round(𝑊 ) and 𝑊 obeys a bivariate 𝒩 (0, 𝜎 2 ).
The round function returns the integer number that is closest to 𝑊.

In all the sequel, the simulation is made for neuron 𝑖 = (0, 0) with 𝑡 0 = 0, 𝑡 1 = 100 (see Algorithm 3).

The parameters 𝑀, 𝜆 ∅ and 𝜎 vary. The parameters 𝜈 𝑖 = 𝜈 and 𝐴 are fixed accordingly by

𝜈 = 0.9𝑀𝜆 ∅ and 𝐴 = 0.9𝑀 -1 (1 -𝜆 ∅ ) -1 ,
to ensure that 𝜙 ∅ 𝑖 < 𝑀 and (2.5.1), which amounts here to (1 -𝜆 ∅ )𝐴𝑀 < 1. On Figure 2.4(a), with 𝑀 = 2, 𝜎 = 1 and 𝜆 ∅ small, we see the overall spread of the algorithm around the neuron to simulate (here (0, 0)). Because we chose a Gaussian variable with small variance for the 𝜆 𝑖 's, the spread is weak and the neurons very close to the neuron to simulate are requested a lot of time at Steps 9-11 of the algorithm. This is also where the algorithm spent the biggest amount of time to simulate Poisson processes. Note also that roughly to simulate 80 points, we need to simulate 10 times more points globally in the infinite network. Remark also on Figure 2.4(b), the avalanche phenomenon, typical of Hawkes processes: for instance the small cluster of black points on neuron with label 0 (i.e. (0,0)) around time 22, is likely to be due to an excitation coming for the spikes generated (and accepted) on neuron labeled 8 and self excitation. The beauty of the algorithm is that we do not need to have the whole time line of neuron 8 to trigger neuron 0, but only the small blue pieces: we just request them at random, depending on the Kalikow decomposition.

On Figure 2.5, we can first observe that when the parameter which governs the range of 𝜆 𝑖 's increase, the global spread of the algorithm increase. In particular, comparing the top left of Figure 2.5 to Figure 2.4 where the only parameter that changes is 𝜎, we see that the algorithm is going much further away and simulates much more points for a sensible equal number of points to generate (and accept) on neuron (0,0). Moreover we can observe that

Illustration

• From left to right, by increasing 𝜆 ∅ , it is more likely to pick an empty neighborhood and as a consequence, the spread of the algorithm is smaller. By increasing 𝜈 = 0.9𝑀𝜆 ∅ , this also increases the total number of points produced on neuron (0;0).

• From top to bottom, by increasing 𝑀, there are more points which are simulated in the Poisson processes (Step 12 of Algorithm 3) and there is also a stronger interaction (we do not truncate that much the number of points in 𝜙 𝑣 ). Therefore, the spread becomes larger and more uniform too, because there are globally more points that are making requests. Moreover, by having a basic rate 𝑀 which is 10 times bigger, we have to simulate roughly 10 times more points. For each neuron in Z 2 , that have been requested in Steps 9:11, except the neuron of interest (0, 0), have been counted the total number of requests, that is the number of time a 𝑉 𝑇 pointed towards this neuron (Steps 9 and 11) and the total time spent at this precise neuron simulating a homogeneous Poisson process (Step 12). Note that since the simulation is on [0, 100] the time spent at position (0, 0) is at least 100. On (a), the summary for one simulation with below the plot, number of points accepted at neuron (0, 0) and total number of points that have been simulated. Also annotated on (a), with labels between 0 and 8, the 9 neurons for which the same simulation in [20,40] is represented in more details on (b). More precisely on (b), in abscissa is time and the neuron labels in ordinate. A plain dot represents an accepted point, by the thinning step (Step 20 of Algorithm 3), and an empty round, a rejected point. The blue pieces of line represent the non empty neighborhoods that are in V.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 1 2 3 4 5 6
• • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • (b) Extract
2 Event-Scheduling Algorithms with Kalikow Decomposition for Simulating Potentially Infinite Neuronal Networks

• • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -20 0 20 -20 0 20
Accepted at (0,0): 86 Produced : 2001

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -30
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -20 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -20 0 20 -20 0 20
Accepted at (0,0): 969 Produced : 15378

�.� Conclusion

We derived a new algorithm for simulating the behavior of one neuron embedded in an infinite network. This is possible thanks to the Kalikow decomposition which allows picking at random the influencing neurons. As seen in the last section, it is computationnally tractable in practice to simulate open systems in the physical sense. A question that remains open for future work is whether we can prove that such a decomposition exists for a wide variety of processes, as it has been shown in discrete time (see [8,9,18]).

�.� Link between Algorithm � and the Kalikow decomposition

To prove that Algorithm 2 returns the desired processes, let us use some additional and more mathematical notation. Note that all the points simulated on neuron 𝑖 before being accepted or not can be seen as coming from a common Poisson process of intensity 𝑀, denoted Π 𝑖 . For any 𝑖 ∈ I, we denote the arrival times of Π 𝑖 , (𝑇 𝑖 𝑛 ) 𝑛∈Z , with 𝑇 𝑖 1 being the first positive time.

Link between Algorithm 2 and the Kalikow decomposition

As in Step 6 of Algorithm 2, we attach to each point of Π 𝑖 a stochastic mark 𝑋 given by,

𝑋 𝑖 𝑛 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if 𝑇 𝑖 𝑛 is accepted in the thinning procedure 0 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒. (2.8.1)
Let us also define 𝑉 𝑖 𝑛 the neighborhood choice of 𝑇 𝑖 𝑛 picked at random and independently of anything else according to 𝜆 𝑖 and shifted at time 𝑇 𝑖 𝑛 . In addition, for any 𝑖 ∈ I, define 𝑁 𝑖 = (𝑇 𝑖 𝑛 , 𝑋 𝑖 𝑛 ) 𝑛∈Z an 𝐸-marked point process with 𝐸 = {0; 1}. In particular, following the notation in Chapter VIII of [2], for any 𝑖 ∈ I, let

𝑁 𝑖 𝑡 (𝑚𝑎𝑟𝑘) = ∑ 𝑛∈Z 1 𝑋 𝑖 𝑛 =𝑚𝑎𝑟𝑘 1 𝑇 𝑖 𝑛 ≤𝑡 for 𝑚𝑎𝑟𝑘 ∈ 𝐸 ℱ 𝑁 𝑡 = ⋁ 𝑖∈I 𝜎(𝑁 𝑖 𝑠 (0), 𝑁
E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝑑𝑁 𝑖 𝑡 (1)
⎞ ⎟ ⎟ ⎠ = ∞ ∑ 𝑛=1 E ( 𝐶 𝑇 𝑖 𝑛 1 𝑋 𝑖 𝑛 =1 )
Note that by Theorem T35 at Appendix A1 of [2], any point 𝑇 should be understood as a stopping time, and that by Theorem T30 at Appendix A2 of [2],

ℱ 𝑁 𝑇 -= ⋁ 𝑗 𝜎{𝑇 𝑗 𝑚 , 𝑋 𝑗 𝑚 such that 𝑇 𝑗 𝑚 < 𝑇 } So E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝑑𝑁 𝑖 𝑡 (1) ⎞ ⎟ ⎟ ⎠ = ∞ ∑ 𝑛=1 E ( 𝐶 𝑇 𝑖 𝑛 E(1 𝑋 𝑖 𝑛 =1 |ℱ 𝑁 𝑇 𝑖 𝑛 -, 𝑉 𝑖 𝑛 ) ) = ∞ ∑ 𝑛=1 E ⎛ ⎜ ⎜ ⎝ 𝐶 𝑇 𝑖 𝑛 𝜑 𝑉 𝑖 𝑛 𝑖 (𝑇 𝑖 𝑛 ) 𝑀 ⎞ ⎟ ⎟ ⎠ .
Let us now integrate with respect to the choice 𝑉 𝑖 𝑛 , which is independent of anything else.

E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝑑𝑁 𝑖 𝑡 (1)
⎞ ⎟ ⎟ ⎠ = ∞ ∑ 𝑛=1 E ⎛ ⎜ ⎜ ⎜ ⎝ 𝐶 𝑇 𝑖 𝑛 𝜆 𝑖 (∅)𝜑 ∅ 𝑖 + ∑ 𝑣∈𝒱 ,𝑣≠∅ 𝜆 𝑖 (𝑣) × 𝜑 𝑣 𝑖 (𝑇 𝑖 𝑛 ) 𝑀 ⎞ ⎟ ⎟ ⎟ ⎠ = E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝜑 𝑖 (𝑡) 𝑀 𝑑Π 𝑖 (𝑡) ⎞ ⎟ ⎟ ⎠ .
Since Π 𝑖 is a Poisson process with respect to (ℱ 𝑁 𝑡 ) 𝑡 with intensity 𝑀, and since 𝐶 𝑡 𝜑 𝑖 (𝑡)

𝑀 is ℱ 𝑁 𝑡-measurable, we finally have that

E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝑑𝑁 𝑖 𝑡 (1) ⎞ ⎟ ⎟ ⎠ = E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝜑 𝑖 (𝑡)𝑑𝑡 ⎞ ⎟ ⎟ ⎠
, which ends the proof.

�.� Proof of Proposition �

Proof. We do the proof for the backward part, starting with 𝑇 = 𝑇 𝑛𝑒𝑥𝑡 as the next point after 𝑡 0 (Step 4 of Algorithm 3), the proof being similar for the other 𝑇 𝑛𝑒𝑥𝑡 generated at Step 23. We construct a tree with root (𝑖, 𝑇 ). For each point (𝑗 𝑇 ′ , 𝑇 ′ ) in the tree, the points which are simulated in 𝑉 𝑇 ′ (Step 12 of Algorithm 3) define the children of (𝑗 𝑇 ′ , 𝑇 ′ ) in the tree. This forms the tree T.

Let us now build a tree C with root (𝑖, 𝑇 ) (that includes the previous tree) by mimicking the previous procedure in the backward part, except that we simulate on the whole neighborhood even if it has a part that intersects with previous neighborhoods (if they exist) (Step 11-12 of Algorithm 3). By doing so, we make the number of children at each node independent of anything else.

If the tree C goes extinct then so does the tree T and the backward part of the algorithm terminates.

But if one only counts the number of children in the tree C, we have a marked branching process whose reproduction distribution for the mark 𝑖 is given by • no children with probability 𝜆 𝑖 (∅)

• Poissonian number of children with parameter 𝑙(𝑣)𝑀 if 𝑣 is the chosen neighborhood with probability 𝜆 𝑖 (𝑣)

This gives that the average number of children issued from a node with the mark 𝑖 is

𝜁 𝑖 = 𝜆 𝑖 (∅) × 0 + ∑ 𝑣∈𝒱 ,𝑣≠∅
𝜆 𝑖 (𝑣)𝑙(𝑣)𝑀.

If we denote C

𝑘 as the collection of points in the tree C at generation 𝑘, and by 𝐾 𝑇 ′ the set of points generated independently as a Poisson process of rate 𝑀 inside 𝑉 𝑇 ′ , we see recursively that

C 𝑘+1 = ⋃ 𝑇 ′ ∈ C 𝑘 𝐾 𝑇 ′ But E(|𝐾 𝑇 ′ ||𝑇 ′ ) = 𝜁 𝑗 𝑇 ′ .
Therefore, if we denote the total number of sites in C 𝑘 by 𝑍 (𝑘) , we have

E(𝑍 (𝑘+1) | C 𝑘 ) ≤ 𝑍 (𝑘) sup 𝑖∈𝐼 𝜁 𝑖 .
2.9 Proof of Proposition 1

One can then conclude by recursion that,

E(𝑍 (𝑘) ) ≤ (sup 𝑖∈𝐼 𝜁 𝑖 ) 𝑘 < 1.
The last inequality use the sparsity neighborhood assumption. Then we deduce that, the mean number of children in each generation goes to 0 as 𝑘 tends to infinity. So by using classical branching techniques in [14], we conclude that the tree C will go extinct almost surely. This also implies that, the backward steps end a.s.

�.� abstract

We propose a new Kalikow decomposition and the corresponding Perfect Simulation algorithm for continuous time multivariate counting processes, on potentially infinite networks. We prove the existence of such a decomposition in various cases. This decomposition is not unique and we discuss the choice of the decomposition in terms of algorithmic efficiency. We apply these methods on several examples: linear

Hawkes process, age dependent Hawkes process, exponential Hawkes process.

�.� Introduction

Multivariate point (or counting) processes on networks have been used to model a large variety of situations : social networks [11], financial prices [2], genomics [22], etc. One of the most complex network models comes from neuroscience where the number of nodes can be as large as billions [16,19,23]. Several counting process models have been used to model such large networks: Hawkes processes [12,13],

Galves-Löcherbach models [10] etc. If the simulation of such large and potentially infinite networks is of fundamental importance in computational neuroscience [16,19], also the existence of such processes in stationary regime and within a potentially infinite network draws a lot of interest (see [10,14,18] in discrete or continuous time).

Kalikow decompositions [15] have been introduced and mainly used in discrete time. Such a decomposition provides a decomposition of the transition probabilities into a mixture of more elementary transitions. The whole idea is that even if the process is complex (infinite memory, infinite network), the elementary transitions look only at what happens in a finite neighborhood in time and space. Once the decomposition is proved for a given process, this can be used to write a Perfect Simulation algorithm [10,14,18]. Here the word "perfect" refers to the fact that it is possible in finite time to simulate what happens on one of the nodes of the potentially infinite network in a stationary regime. The existence of such an algorithm guarantees in particular the existence of the process in stationary regime. Most of the papers referring to Kalikow decomposition and Perfect Simulation are theoretical and aim at proving the existence of such processes on infinite networks in stationary regime [10,14].

In the present paper, we propose to go from discrete to continuous time. Therefore, we will decompose conditional intensities rather than transition probabilities. This leads to serious difficulties that usually prevent a more practical application of the Perfect Simulation algorithm. Indeed, up to our knowledge, the only work dealing with continuous time counting processes, is the one by Hodara and Löcherbach [14]. Their decomposition is constructed under the assumption that there is a dominating Poisson process on each of the nodes, from which the points of the processes under interest can be thinned by rejection sampling (see also [17] for another use of thinning in simulation of counting processes). To prove the existence of a Kalikow decomposition and go back to a more classic discrete time setting, the authors need to freeze the dominating Poisson process, leading to a mixture, in the Kalikow decomposition, that depends on the realization of the dominating Poisson process. Such a mixture is not accessible in practice, and this prevents the use of their Perfect Simulation algorithm for more concrete purposes than mere existence.

More recently, in a previous computational article [19], we have used another type of Kalikow decomposition, which does not depend on the dominating Poisson process. This leads to a Perfect Simulation algorithm, which can be used as a concrete way for Computational Neuroscience to simulate neuronal

networks as an open physical system, where we do not need to simulate the whole network to simulate what happens in a small part [19].

In the present work, we want to go further, by proposing an even more general Kalikow decomposition, which does not assume the existence of a dominating Poisson process at all. We also prove (and this is not done in [19]) that such decomposition exists on various interesting examples, even if these decompositions are not unique. Finally we propose a corresponding Perfect Simulation algorithm and we discuss its efficiency with respect to the decomposition that is used.

The paper is organized as follows. In Section 3.3, we introduce the basic notation and give the precise definition of a Kalikow decomposition. In Section 3.4, we present two methods to obtain a Kalikow decomposition for a counting process having stochastic intensity. The first method is very general and based on adequate choices of neighborhoods and weights. As an application, we study a classical example, the linear Hawkes process [12,13], and a very recent and promising process, the age dependent Hawkes process [20]. Though very simple cases of Hawkes processes are treated, it is worth to note that our method can be generalized very easily. The second method exists only for Hawkes processes and is based on Taylor expansion. Finally, in Section 3.5, we present a modified Perfect Simulation algorithm based on the Kalikow decomposition written in Section 3.4, and we discuss the efficiency of the algorithm with respect to the Kalikow decomposition.

�.� Notation and Kalikow decomposition �.�.� Notation and Definition

We start this section by recalling the definition of counting processes and stochastic intensity. We refer the reader to [3] and [7] for more complete statements.

Let I be a countable index set. A counting process 𝑍 𝑖 , 𝑖 ∈ I, can be described by its sequence of jump times in R, (𝑇 𝑖 𝑛 ) 𝑛∈Z [3]. Consider (ℱ 𝑡 ) 𝑡∈R the past filtration of the process 𝑍 = (𝑍 𝑖 ) 𝑖∈I :

ℱ 𝑡 = 𝜎(𝑍 𝑖 𝑠 , 𝑖 ∈ I, 𝑠 ≤ 𝑡).
Since a point process is fully characterized by its arrival times [3], we can denote by 𝒳 the canonical path space of 𝑍:

𝒳 = {({𝑡 𝑖 𝑛 } 𝑛∈Z ) 𝑖∈I such that ∀𝑛, 𝑖, 𝑡 𝑖 𝑛 < 𝑡 𝑖 𝑛+1 and 𝑡 𝑖 0 ≤ 0 < 𝑡 𝑖 1 },
where {𝑡 𝑖 𝑛 } 𝑛∈Z denotes a possible realization of (𝑇 𝑖 𝑛 ) 𝑛∈Z . Denote 𝒳 𝑡 the canonical path space of 𝑍 before time 𝑡:

𝒳 𝑡 = 𝒳 ∩ (-∞, 𝑡) 𝐼 .
A past configuration 𝑥 𝑡 is an element of 𝒳 𝑡 which is a realization of arrival times of 𝑍 before 𝑡.

Under suitable assumptions, the evolution of the point process 𝑍 𝑖 with respect to (ℱ 𝑡 ) 𝑡∈R is fully characterized by its stochastic intensity which depends on the past configuration, see Proposition 7.2.IV of [7]. Hence, in this paper, for any 𝑥 𝑡 ∈ 𝒳 𝑡 , given that the past before time 𝑡 is 𝑥 𝑡 , we denote by 𝜙 𝑖,𝑡 (𝑥 𝑡 ) the corresponding stochastic intensity of the process 𝑍 𝑖 at time 𝑡 for any 𝑖 ∈ I. More precisely, for any 𝑥 𝑡 ∈ 𝒳 𝑡 , we have

P ( 𝑍 𝑖 has jump in [𝑡, 𝑡 + 𝑑𝑡) | past before time 𝑡 = 𝑥 𝑡) = 𝜙 𝑖,𝑡 (𝑥 𝑡 )𝑑𝑡.
Together with the path space 𝒳 𝑡 , we denote V 𝑡 a countable collection of finite space-time neighborhoods 𝑣 𝑡 , in which each neighborhood 𝑣 𝑡 is a Borel subset of I × (-∞, 𝑡). More precisely, we call 𝑣 𝑡 a finite neighborhood if there exists a finite subset 𝐽 ⊂ I and a finite interval [𝑎, 𝑏] such that:

𝑣 𝑡 ⊂ 𝐽 × [𝑎, 𝑏].
For convenience, we denote 𝒳 the canonical path space of 𝑍 before time 0 instead of 𝒳 0 . In addition, a past configuration before 0 is denoted by 𝑥 and the intensity at time 0 is 𝜙 𝑖 (𝑥) instead of 𝜙 𝑖,0 (𝑥 0 ). Moreover, We say process 𝑍 = (𝑍 𝑖 ) 𝑖∈I satisfies a Kalikow decomposition w.r.t (V 𝑡 ) 𝑡∈R and (𝒴 𝑡 ) 𝑡∈R if for all 𝑖 ∈ I each process 𝑍 𝑖 satisfies a Kalikow decomposition w.r.t (V 𝑡 ) 𝑡∈R and (𝒴 𝑡 ) 𝑡∈R . For simplicity, all the conventions of (𝒳 𝑡 ) 𝑡∈R will be used for (𝒴 𝑡 ) 𝑡∈R , such as 𝒴 0 being replaced by 𝒴, etc.

if 𝑥 = ( {𝑡 𝑖 𝑛 } 𝑛∈Z -) 𝑖∈I ∈ 𝒳 ,
Remark 5. Note that the function 𝜆 𝑖,𝑡 (.) in Definition 3.3.3 is a deterministic function, that is why this decomposition is unconditional, whereas in [14], 𝜆 𝑖,𝑡 (.) was depending on the dominating Poisson processes (see Introduction). Secondly, we do not restrict ourself to a bounded intensity, which is a notable improvement compared to [19].

In the following subsection, we show that we can translate the decomposition at time 0 to any time 𝑡 when dealing with time homogeneous processes. �.�.� From the decomposition at time � to the decomposition at any time 𝑡.

In the context of time homogeneous process, for any 𝑡 ∈ R, we define

𝒴 𝑡 = 𝒴 →𝑡 = {𝑥 + 𝑡|𝑥 ∈ 𝒴 }.
For all 𝑖 ∈ I, assume that Equation (3.3.1) is satisfied at time 𝑡 = 0 for any 𝑥 ∈ 𝒳 ∩ 𝒴, we then prove that, for time homogeneous processes, this equation is achieved at any time 𝑡, for any 𝑥 𝑡 ∈ 𝒳 𝑡 ∩ 𝒴 𝑡 , for a particular choice of V 𝑡 . Hence, to show that a counting process 𝑍 𝑖 , 𝑖 ∈ I satisfies a Kalikow decomposition, it is then sufficient to write the Kalikow decomposition at time 0 only. Consider a neighborhood family at time 0, V 0 . Take a neighborhood 𝑣 ∈ V 0 and for 𝑡 ≥ 0, denote 𝑣 →𝑡 = {(𝑖, 𝑢 + 𝑡) ∶ (𝑖, 𝑢) ∈ 𝑣}, i.e, shift to the right the time component of the neighborhood 𝑣 by 𝑡 and define This shows that Equation (3.3.1) is satisfied at any time 𝑡. Thus, by definition, 𝑍 𝑖 admits a Kalikow decomposition w.r.t the neighborhood family (V 𝑡 ) 𝑡∈R where V 𝑡 = V →𝑡 . In addition, (V 𝑡 ) 𝑡∈R is completely determined once the neighborhood family at time 0, V 0 , is determined. From now on, it is sufficient to mention V 0 whenever we speak about the Kalikow decomposition.

V 𝑡 = V →𝑡 ∶= {𝑣 →𝑡 ∶ 𝑣 ∈ V 0 }. ( 3

�.�.� About the subspace 𝒴

The role of the subspace 𝒴 is to make the Kalikow decomposition achievable. There are many possible choices for such a subspace, depending on the model under consideration. In this paper, we focus for instance on the choice 𝒴 = 𝒳 >𝛿 , the subspace of 𝒳 where the distance between any two consecutive possible jumps is greater than 𝛿. More precisely,

𝒳 >𝛿 = {𝑥 = ({𝑡 𝑖 𝑛 } 𝑛∈Z -) 𝑖∈I such that ∀𝑛, 𝑖 𝑡 𝑖 𝑛+1 -𝑡 𝑖 𝑛 > 𝛿 and 𝑡 𝑖 0 ≤ 0}. (3.3.4)
Another possible choice of a subspace 𝒴 that we consider guarantees that the intensity is finite, by introducing 𝒴 = {𝑥 ∈ 𝒳 , ∀𝑖 𝜙 𝑖 (𝑥) < ∞}.

In practice, either the process is known, namely, the formula of the intensity is given, and we dispose of theoretical results that guarantee for instance that 𝑥 ∈ 𝒴 almost surely, or we handle the decomposition by creating a certain threshold. Before the process reaches this threshold, we are able to achieve the Kalikow decomposition by following one of the methods described in the following section. However, after having reached this threshold we can not say anything any more. This depends of course on the initial condition, and if we start for instance with no point on (-∞, 0), the intensity is in many cases likely to be finite at least at time 0. To proceed further with a Perfect Simulation algorithm is more tricky since we need to check that the intensity remains finite during all steps of the algorithm. This will be the main object of a future work.

�.� Main results

�.�.� The first method

In this section, we present step by step a general method to prove the existence of a Kalikow decomposition for the counting process 𝑍 𝑖 , 𝑖 ∈ I. We start by discussing the relevant family of neighborhoods.

For any subset 𝐽 ⊂ I we say a process 𝑍 𝑖 is locally dependent on a subprocess 𝑍 𝐽 ∶= (𝑍 𝑗 ) 𝑗∈𝐽 if the intensity 𝜙 𝑖 (𝑥) is a cylindrical function on 𝐽 × (-∞, 0). Roughly speaking, a process 𝑍 𝑖 is locally dependent on a subprocess 𝑍 𝐽 if the information of process 𝑍 𝑗 with 𝑗 ∈ 𝐽 is compulsory to compute the

Main results

intensity of process 𝑍 𝑖 . For short, we say 𝑖 is locally dependent on 𝐽. A more formal definition of local dependence can be found at [9]. Denote 𝒱 .→𝑖 ∶= {𝑗 ∈ I such that 𝑖 is locally dependent on 𝑗} (3.4.1)

and 𝒮 𝑖 ∶= 𝒱 .→𝑖 × (-∞, 0).

We denote by (P, ≤) a countable, ordered set. To simplify notation, in what follows, we consider P to be N, but depending on the concrete examples, different sets P will be considered, for example N × N, etc. Define a family of finite space-time neighborhoods at time 0 associated to the dynamic of 𝑍 𝑖 :

V 𝑖 = {(𝑣 𝑖 𝑘 ) 𝑘∈P ⊂ I × (-∞, 0) such that ∪ 𝑘 𝑣 𝑖 𝑘 = 𝒮 𝑖 } (3.4.2)
with the convention that 𝑣 𝑖 0 = ∅. In order to prove that 𝑍 𝑖 admits a Kalikow decomposition w.r.t the neighborhood family V 𝑖 and a convenient subspace 𝒴, we use the following condition. Once this condition is fulfilled, we will show in Proposition 3 below that we can derive the corresponding Kalikow decomposition. Assumption 1. There exists a non-negative sequence of functions (Δ 𝑖 𝑘 (𝑥)) 𝑘∈P which are cylindrical on 𝑣 𝑖 𝑘 such that

𝑛 ∑ 𝑘=0 Δ 𝑖 𝑘 (𝑥) → 𝜙 𝑖 (𝑥)
as 𝑛 → ∞, for every 𝑥 ∈ 𝒳 ∩ 𝒴. Now, we present a method to obtain a Kalikow decomposition under Assumption 1.

Step 1: For any 𝑘 ≥ 0, take Δ 𝑖 𝑘 (𝑥) from Assumption 1, it then guarantees that 𝜙 𝑖 (𝑥) can be written as follows:

𝜙 𝑖 (𝑥) = Δ 𝑖 0 + ∑ 𝑘≥1 Δ 𝑖 𝑘 (𝑥)
where Δ 𝑖 0 ∶= Δ 𝑖 0 (𝑥) is a constant, that is, it does not depend on 𝑥 since 𝑣 𝑖 0 = ∅. Note that, for 𝑘 ≥ 1, by definition, Δ 𝑖 𝑘 (𝑥) is cylindrical on 𝑣 𝑖 𝑘 and non negative.

Step 2: Define a deterministic nonnegative sequence (𝜂 𝑖 𝑘 ) 𝑘∈P such that

∑ 𝑘≥0 𝜂 𝑖 𝑘 = 1.
Obviously, 𝜙 𝑖 (𝑥) can be written as:

𝜙 𝑖 (𝑥) = 𝜂 𝑖 0 Δ 𝑖 0 𝜂 𝑖 0 + ∑ 𝑘≥1 𝜂 𝑖 𝑘 Δ 𝑖 𝑘 (𝑥) 𝜂 𝑖 𝑘
with the convention that 0/0 = 1. We can set 𝜂 𝑖 𝑘 to 0 whenever Δ 𝑖 𝑘 (𝑥) equals 0 for all 𝑥 and even discard the corresponding neighborhood from V 𝑖 . Proposition 3. Consider a point process 𝑍 𝑖 , 𝑖 ∈ I, with intensity at time 0, 𝜙 𝑖 (𝑥), for any 𝑥 ∈ 𝒳. For the neighborhood family V 𝑖 defined in (3.4.2), if 𝜙 𝑖 (𝑥) satisfies Assumption 1, then 𝑍 𝑖 has the following Kalikow decomposition with respect to V 𝑖 and the subspace 𝒴:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 𝜆 𝑖 (∅) = 𝜂 𝑖 0 𝜙 ∅ 𝑖 = Δ 𝑖 0 𝜂 𝑖 0 𝜆 𝑖 (𝑣 𝑖 𝑘 ) = 𝜂 𝑖 𝑘 𝜙 𝑣 𝑖 𝑘 𝑖 (𝑥) = Δ 𝑖 𝑘 (𝑥) 𝜂 𝑖 𝑘
for any choice of non negative weights (𝜂 𝑖 𝑘 ) 𝑘∈P such that

∑ 𝑘≥0 𝜂 𝑖 𝑘 = 1.

Remark 6. From

Step 2, it is clear that such a Kalikow decomposition is not unique. However, to perform the Perfect Simulation algorithm (see Section 3.5), (𝜂 𝑖 𝑘 ) 𝑘∈P can not be chosen arbitrarily. These weights need to satisfy several conditions, for example: the stopping condition, see Proposition 7 below, which enables the algorithm to end. On the other hand, for the simulation purposes, these weights should be chosen carefully to avoid the explosion of 𝜙 𝑣 𝑖 𝑘 𝑖 (𝑥). In addition, they also influence the mean number of total simulated points of this algorithm, and therefore the efficiency of the algorithm (see Section 3.5). To illustrate this point more clearly, one example is considered in Section 3.5.

�.�.� Examples of the first method

Linear Hawkes process.

In the following, we consider a linear Hawkes process [12,13], where the intensity at time 0, 𝜙 𝑖 (𝑥), is given as follows. For any 𝑥 ∈ 𝒳,

𝜙 𝑖 (𝑥) = 𝜇 𝑖 + ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ,
where ℎ 𝑗𝑖 (.) measures the local dependence of process 𝑍 𝑖 on 𝑍 𝑗 and 𝜇 𝑖 refers to the spontaneous rate of process 𝑍 𝑖 . We consider the following assumption for the method to work. Assumption 2. For all 𝑖, 𝑗 ∈ I, ℎ 𝑗𝑖 (.) is a non negative function.

Denote 𝜖 an arbitrary, fixed and positive number. Then we have

𝜙 𝑖 (𝑥) = 𝜇 𝑖 + ∑ 𝑗∈I ∑ 𝑛∈N * ∫ -𝑛𝜖+𝜖 -𝑛𝜖 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 .
Furthermore, in this example, we have

𝒱 .→𝑖 = {𝑗 ∈ I such that ℎ 𝑗𝑖 ≠ 0}.
Since I is at most countable, we may enumerate the indices in I by 𝑗 𝑙 for 𝑙 ∈ N * . In this example we put P = N * × N * , and we consider the neighborhood family V 𝑎𝑡𝑜𝑚 as follows: for 𝑙, 𝑛 ≥ 1 𝑣 𝑖 (𝑙,𝑛) ∶= {𝑗 𝑙 } × [-𝑛𝜖, -𝑛𝜖 + 𝜖).

We define, for all 𝑙, 𝑛 ≥ 1:

Δ 𝑖 (𝑙,𝑛) (𝑥) ∶= ∫ -𝑛𝜖+𝜖 -𝑛𝜖 ℎ 𝑗 𝑙 𝑖 (-𝑠)𝑑𝑥 𝑗 𝑙 𝑠 ≥ 0,
and Δ 𝑖 (0,0) ∶= Δ 𝑖 (0,0) (𝑥) ∶= 𝜇 𝑖 . Applying the Monotone Convergence theorem, we obtain

Δ 𝑖 (0,0) + (𝐿,𝑁) ∑ (𝑙,𝑛)=(1,1)
Δ 𝑖 (𝑙,𝑛) (𝑥) → 𝜙 𝑖 (𝑥),

when (𝐿, 𝑁) → (∞, ∞).
Notice that the intensity 𝜙 𝑖 (𝑥) is finite if and only if the series converges. However, what we do here is more general, since we allow the intensity to be infinite.

Thus, Assumption 1 is fulfilled. As a consequence, relying on Proposition 3, we conclude that 𝑍 𝑖 has a Kalikow decomposition with respect to V 𝑎𝑡𝑜𝑚 with 𝜂 𝑖 (𝑙,𝑛) = 1.

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 𝜆 𝑖 (∅) = 𝜂 (0,0) 𝜙 ∅ 𝑖 = 𝜇 𝑖 𝜂 (0,
Remark 7. Note that under mild conditions (see for instance [8]) it is well-known that the process exists in a stationary regime. On the other hand, for our purpose it is sufficient to consider the subspace 𝒴 = {𝑥 ∈ 𝒳 | ∀𝑖 𝜙 𝑖 (𝑥) < ∞} such that for all 𝑥 ∈ 𝒳 ∩ 𝒴, 𝜙 𝑖 (𝑥) and 𝜙 𝑣 𝑖 (𝑥) are finite.

Age dependent Hawkes process with hard refractory period.

In this section, we are interested in writing a Kalikow decomposition for Age dependent Hawkes processes with hard refractory period. Up to our knowledge, this process was first introduced in [6] and no Kalikow decomposition has been proved, even in a conditional framework, for this type of process. In our setting, the stochastic intensity of an Age dependent Hawkes process with hard refractory of length 𝛿 > 0 can be written as follows. For any 𝑖 ∈ I and 𝑥 ∈ 𝒳,

𝜙 𝑖 (𝑥) = 𝜓 𝑖 ( ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ) 1 𝑎 𝑖 0 (𝑥)>𝛿 (3.4.3)
with the convention that 𝑎 𝑖 0 (𝑥) = -𝐿 𝑖 0 (𝑥) and

𝐿 𝑖 0 (𝑥) = sup{𝑡 𝑖 𝑘 ∈ 𝑥 such that 𝑡 𝑖 𝑘 < 0} = 𝑡 𝑖 -1
is the last jump before 0 of process 𝑍 𝑖 , if it has at least one jump. Otherwise, we put 𝐿 𝑖 0 (𝑥) = -∞. Again, we have

𝒱 .→𝑖 = {𝑗 ∈ I such that ℎ 𝑗𝑖 ≠ 0} ∪ {𝑖}.
By definition of stochastic intensity (3.4.3), we observe that the distance of any two consecutive jumps have to be larger than 𝛿. This observation leads us to consider the subspace 𝒴 = 𝒳 >𝛿 that is introduced in (3.3.4). To prove a Kalikow decomposition, we consider the following assumptions. (ii) For every 𝑖, 𝜓 𝑖 (.) is an increasing, non negative continuous function.

On the other hand, we build a neighborhood family V 𝑛𝑒𝑠𝑡𝑒𝑑 by introducing a non decreasing sequence (𝑉 𝑖 (𝑘)) 𝑘≥0 of finite subsets of I such that 𝑉 𝑖 (0) = ∅, 𝑉 𝑖 (1) = {𝑖}, 𝑉 𝑖 (𝑘 -1) ⊂ 𝑉 𝑖 (𝑘) and ∪ 𝑘 𝑉 𝑖 (𝑘) = 𝒱 .→𝑖 ∪ {𝑖 }.

(3.4.4)

Consider P = N and set 𝑣 𝑖 𝑘 = 𝑉 𝑖 (𝑘) × [-𝑘𝛿, 0), then we obtain by construction an increasing, nested neighborhood sequence (𝑣 𝑖 𝑘 ) 𝑘∈N . By applying the methodology of Section 3.4.1, we prove the following Proposition.

Proposition 4. Consider 𝒳 >𝛿 defined as in (3.3.4). For any past configuration 𝑥 ∈ 𝒳, suppose that the intensity is of the form

𝜙 𝑖 (𝑥) = 𝜓 𝑖 ( ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ) 1 𝑎 𝑖 0 (𝑥)>𝛿 .
If Assumption 3 is fulfilled, then the corresponding Kalikow decomposition with respect to V 𝑛𝑒𝑠𝑡𝑒𝑑 and

𝒳 >𝛿 is for 𝑘 ≥ 1 ⎧ ⎪ ⎨ ⎪ ⎩ 𝜆 𝑖 (𝑣 𝑖 𝑘 ) = 𝜂 𝑖 𝑘 𝜙 𝑣 𝑖 𝑘 𝑖 (𝑥) = Δ 𝑖 𝑘 (𝑥) 𝜂 𝑖 𝑘 with (i) any positive sequence (𝜂 𝑖 𝑘 ) 𝑘≥1 such that ∑ 𝑘≥1 𝜂 𝑖 𝑘 = 1, (ii)
Δ 𝑖 𝑘 (𝑥) = 𝜓 𝑖 ⎛ ⎜ ⎜ ⎝ ∑ 𝑗∈𝑉 𝑖 (𝑘) 0 ∫ -𝑘𝛿 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ⎞ ⎟ ⎟ ⎠ 1 𝑎 𝑖 0 (𝑥)>𝛿 -𝜓 𝑖 ⎛ ⎜ ⎜ ⎝ ∑ 𝑗∈𝑉 𝑖 (𝑘-1) 0 ∫ -𝑘𝛿+𝛿 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ⎞ ⎟ ⎟ ⎠ 1 𝑎 𝑖 0 (𝑥)>𝛿
for 𝑘 ≥ 2 and Δ 1 𝑖 (𝑥) = 𝜓 𝑖 (0)1 𝑎 𝑖 0 (𝑥))>𝛿 .

Remark 8. Note that 𝜆 𝑖 (∅) does not appear in the above Kalikow decomposition. Amazingly, this does not cause any problems for the Perfect Simulation algorithm (see Section 3.5). We stress that this is one of the main differences with [14], the other one being that the decomposition does not depend on the dominating Poisson processes. By Assumption 3, (Δ 𝑖 𝑘 (𝑥)) 𝑘∈N is well-defined, non negative and cylindrical on 𝑣 𝑖 𝑘 . Moreover, we set Δ 𝑖 0 = 0 as well 𝜂 𝑖 0 = 0 and Δ 1 𝑖 (𝑥) = 𝜓 𝑖 (0)1 𝑎 𝑖 0 (𝑥))>𝛿 . Let

𝑟 [𝑛] 𝑖 (𝑥) ∶= 𝑛 ∑ 𝑘=1 Δ 𝑖 𝑘 (𝑥) = 𝜓 𝑖 ( ∑ 𝑗∈𝑉 𝑖 (𝑛) ∫ 0 -𝑛𝛿 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ) 1 𝑎 𝑖 0 (𝑥)>𝛿
and let us show that 𝑟

[𝑛] 𝑖 (𝑥) → 𝜙 𝑖 (𝑥) when 𝑛 → ∞. Consider the inner-term of the parenthesis,

∑ 𝑗∈𝑉 𝑖 (𝑛) ∫ 0 -𝑛𝛿 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 = ∫ 𝐼×R - ℎ 𝑗𝑖 (-𝑠)1 (𝑗,𝑠)∈𝑣 𝑖 𝑛 𝑑𝑥 𝑗 𝑠 𝑑𝜅 𝑗 ,
where we denote 𝑑𝜅 the counting measure on the discrete set I.

We have that ( ℎ 𝑗𝑖 (-𝑠)1 (𝑗,𝑠)∈𝑣 𝑖 𝑛 ) 𝑛∈Z is non negative and non decreasing sequence in 𝑛. In addition, it converges to ℎ 𝑗𝑖 (-𝑠)1 (𝑗,𝑠)∈𝒱 .→𝑖 ×(-∞,0) as 𝑛 → ∞. Moreover, since 𝜓 𝑖 (.) is a continuous and increasing function, the Monotone convergence theorem for Lebesgue Stieltjes measures implies that 𝑟

𝜓 𝑖 ( ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑧 𝑗 𝑠 ) 1 𝑎 𝑖 0 (𝑥)>𝛿 = 𝜓 𝑖 ( ∑ 𝑗∈𝑉 𝑖 (𝑘) ∫ 0 -𝑘𝛿 ℎ 𝑗𝑖 (-𝑠)𝑑𝑧 𝑗 𝑠 + + ∑ 𝑗∈𝑉 𝑖 (𝑘) ∫ -𝑘𝛿 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑧 𝑗 𝑠 + ∑ 𝑗∉𝑉 𝑖 (𝑘) ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑧 𝑗 𝑠 ) 1 𝑎 𝑖 0 (𝑥)>𝛿 .
By Assumption 3, we conclude that, inf

𝑧∈𝐷 𝑘 𝑖 (𝑥) 𝜓 𝑖 ( ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑧 𝑗 𝑠 ) 1 𝑎 𝑖 0 (𝑥)>𝛿 = 𝜓 𝑖 ( ∑ 𝑗∈𝑉 𝑖 (𝑘) ∫ 0 -𝑘𝛿 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ) 1 𝑎 𝑖 0 (𝑥)>𝛿 .
The equality is achieved when we consider the configuration 𝑧 having points only inside the neighborhood 𝑣 𝑖 𝑘 . Hence, for 𝑘 ≥ 2, we observe that

Δ 𝑖 𝑘 (𝑥) = inf 𝑧∈𝐷 𝑘 𝑖 (𝑥) 𝜓 𝑖 ( ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑧 𝑗 𝑠 ) 1 𝑎 𝑖 0 (𝑥)>𝛿 -inf 𝑧∈𝐷 𝑘-1 𝑖 (𝑥) 𝜓 𝑖 ( ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑧 𝑗 𝑠 ) 1 𝑎 𝑖 0 (𝑥)>𝛿 .
The above prescription corresponds to the classical method of obtaining a Kalikow decomposition in discrete time, discussed in [10,14].

�.�.� Another method for nonlinear Hawkes processes

In this section, we present a second method to prove the existence of a Kalikow decomposition for nonlinear Hawkes processes [4]. In the setting of this section, we suppose that the intensity of the nonlinear Hawkes process 𝑍 𝑖 , 𝑖 ∈ I, is given by

𝜙 𝑖 (𝑥) = 𝜓 𝑖 ( ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 )
for any 𝑥 ∈ 𝒳 .

We work under the following assumptions.

Assumption 4. (i) For any 𝑖, 𝑗 ∈ I, the function ℎ 𝑗𝑖 (.) is non negative and belongs to 𝐿 1 . Moreover, for every 𝑖, we have

∑ 𝑗∈I ‖ℎ 𝑗𝑖 ‖ 𝐿 1 < ∞.
(ii) For every 𝑖 ∈ I, 𝜓 𝑖 (.) is an analytic function on R with radius of convergence about 0 which is given by 𝐾, for some positive 𝐾. Moreover, its derivative of order 𝑛, 𝜓 (𝑛) 𝑖 (0), is non negative for all 𝑛 ≥ 1, and 𝜓 𝑖 (0) is non negative as well.

In this section, to develop our series using Taylor expansion, we choose

𝒴 𝐾 = { 𝑥 | sup 𝑖 ( ∑ 𝑗 ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ) < 𝐾 } .
We now describe our method. We fix 𝜖 an arbitrary positive parameter, and we put 𝐷 = I × N. An index in 𝐷 is 𝛼 = (𝑗, 𝑛) where 𝑗 ∈ I and 𝑛 ∈ N. We put |𝛼| ∶= |𝑗| + 𝑛. For any 𝑥 ∈ 𝒳 ∩ 𝒴 𝐾 , we introduce

𝑎 𝛼 (𝑥) ∶= ∫ -𝑛𝜖 -(𝑛+1)𝜖 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 .
Note that, whenever 𝑥 ∈ 𝒳 ∩ 𝒴 𝐾 , we have 𝑎 𝛼 (𝑥) converges to 0 when |𝛼| → ∞. Furthermore, we have

𝜙 𝑖 (𝑥) = 𝜓 𝑖 ( ∑ 𝑗∈I ∑ 𝑛∈N ∫ -𝑛𝜖 -(𝑛+1)𝜖 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ) = 𝜓 𝑖 ( ∑ 𝑗∈I ∑ 𝑛∈N 𝑎 𝑗𝑛 (𝑥) ) = 𝜓 𝑖 ( ∑ 𝛼∈𝐷 𝑎 𝛼 (𝑥) ) .
We suppose moreover that Assumption 5. For any 𝑖 ∈ I, there exist a deterministic sequence (𝜂 𝑖 𝛼 𝑘 ) 𝛼 𝑘 ∈𝐷 ∈ (0, 1) such that

∑ 𝑛≥1 ∑ 𝛼 1 ,…,𝛼 𝑛 ∈𝐷 𝜂 𝑖 𝛼 1 … 𝜂 𝑖 𝛼 𝑛 < 1.
Let us consider P = ∪ ∞ 𝑛=1 𝐷 𝑛 where 𝐷 𝑛 ∶= 𝐷 × 𝐷 × … × 𝐷 (𝑛 times). We construct the neighborhood family V 𝑇 𝑎𝑦𝑙𝑜𝑟 by defining for 𝛼 1∶𝑘 = (𝛼

𝜙 𝑖 (𝑥) = 𝜓 𝑖 ( ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ) .
Under Assumption 4 and 5, 𝑍 𝑖 has a Kalikow decomposition with respect to the neighborhood family V 𝑇 𝑎𝑦𝑙𝑜𝑟 defined by (3.4.6) and 𝒴 𝐾 as follows:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝜆 𝑖 (∅) = 1 -∑ 𝑘≥1 ∑ 𝛼 1 ,…,𝛼 𝑘 ∈𝐷 𝜂 𝑖 𝛼 1 … 𝜂 𝑖 𝛼 𝑘 𝜙 ∅ 𝑖 = 𝜓 𝑖 (0) 𝜆 𝑖 (∅) 𝜆 𝑖 (𝑣 𝑖 𝛼 1∶𝑘 ) = 𝜂 𝑖 𝛼 1 … 𝜂 𝑖 𝛼 𝑘 𝜙 𝑣 𝑖 𝛼 1∶𝑘 𝑖 (𝑥) = 𝜓 (𝑘) 𝑖 (0) 𝑘! 𝑎 𝛼 1 (𝑥) 𝜂 𝑖 𝛼 1 … 𝑎 𝛼 𝑘 (𝑥) 𝜂 𝑖 𝛼 𝑘 ,
for any weights (𝜂 𝑖 𝛼 𝑘 ) 𝛼 𝑘 ∈𝐷 satisfying Assumption 5 and for 𝛼 = (𝑗, 𝑛), 𝑎 𝛼 defined by

𝑎 𝛼 (𝑥) ∶= ∫ -𝑛𝜖 -(𝑛+1)𝜖 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 .
Proof. From Assumption 4, (ii), and noticing that for any 𝑥 ∈ 𝒳 ∩𝒴 𝐾 , we have 0 ≤ ∑ 𝑗∈I ∑ 𝑛∈N 𝑎 𝑗𝑛 (𝑥) = ∑ 𝛼∈𝐷 𝑎 𝛼 (𝑥) < 𝐾, we have a Taylor expansion of 𝜙 𝑖 (𝑥):

𝜙 𝑖 (𝑥) = 𝜓 𝑖 ( ∑ 𝑗∈I ∑ 𝑛∈N 𝑎 𝑗𝑛 (𝑥) ) = 𝜓 𝑖 ( ∑ 𝛼∈𝐷 𝑎 𝛼 (𝑥) ) = 𝜓 𝑖 (0) + ∑ 𝑘≥1 𝜓 (𝑘) 𝑖 (0) 𝑘! ( ∑ 𝛼=(𝑗,𝑛)∈𝐷 𝑎 𝛼 (𝑥) ) 𝑘 = 𝜓 𝑖 (0) + ∑ 𝑘≥1 ∑ 𝛼 1 …𝛼 𝑘 ∈𝐷 𝜓 (𝑘) 𝑖 (0) 𝑘! 𝑎 𝛼 1 (𝑥) … 𝑎 𝛼 𝑘 (𝑥).
Thus, to obtain the Kalikow decomposition, take (𝜂 𝑖 𝛼 𝑘 ) 𝛼 𝑘 ∈𝐷 from Assumption 5, then the conclusion follows.

Remark 10. The parameter (𝜂 𝑖 𝛼 𝑘 ) 𝛼 𝑘 ∈𝐷 plays the same role here as (𝜂 𝑖 𝑘 ) 𝑘∈N in the previous section. Again, the choice of this parameter is discussed in Section 3.5.

In the following, we give several examples in which we can apply the second method.

�.�.� Examples of second method

In this section, we always consider a deterministic sequence (𝜂 𝑖 𝛼 𝑘 ) 𝛼 𝑘 ∈𝐷 that satisfies Assumption 5.

Example 1. Exponential Hawkes process with 𝜓 𝑖 (.) ∶= exp(.). In this case we have 𝐾 = ∞, 𝜓 𝑖 (0) = 1 and 𝜓 (𝑘) 𝑖 (0) = 1 for 𝑘 ≥ 1. Therefore,

𝜙 𝑖 (𝑥) = exp ( ∑ 𝑗∈I ∫ 0 -∞ ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ) = 1 + ∑ 𝑘≥1 ∑ 𝛼 1 …𝛼 𝑘 ∈𝐷 1 𝑘! 𝑎 𝛼 1 (𝑥) … 𝑎 𝛼 𝑘 (𝑥).
Then, the Kalikow decomposition of 𝑍 𝑖 with respect to the neighborhood family V 𝑇 𝑎𝑦𝑙𝑜𝑟 in (3.4.6) and 𝒴 ∞ is given by

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 𝜆 𝑖 (∅) = 1 -∑ 𝑘≥1 ∑ 𝛼 1 ,…,𝛼 𝑘 ∈𝐷 𝜂 𝑖 𝛼 1 … 𝜂 𝑖 𝛼 𝑘 𝜙 𝑖 (∅) = 1 𝜆 𝑖 (∅) 𝜆 𝑖 (𝑣 𝑖 𝛼 1∶𝑘 ) = 𝜂 𝑖 𝛼 1 … 𝜂 𝑖 𝛼 𝑘 𝜙 𝑣 𝑖 𝛼 1∶𝑘 𝑖 (𝑥) = 1 𝑘! × 𝑎 𝛼 1 (𝑥) 𝜂 𝑖 𝛼 1 … 𝑎 𝛼 𝑘 (𝑥) 𝜂 𝑖 𝛼 𝑘 .
Example 2. We consider the non-linear Hawkes process with 𝜓 𝑖 (𝑢) = 𝑐ℎ(𝑢) where 𝑐ℎ(𝑢) = 𝑒 𝑢 + 𝑒 -𝑢 2 .

Then 𝐾 = ∞, 𝜓 𝑖 (0) = 1, and for 𝑘 ≥ 1, we have 𝜓 Therefore, 𝑍 𝑖 has the following Kalikow decomposition with respect to the neighborhood family V 𝑇 𝑎𝑦𝑙𝑜𝑟 in (3.4.6) and 𝒴 ∞ is given by

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝜆 𝑖 (∅) = 1 -∑ 𝑘≥1 ∑ 𝛼 1 ,…,𝛼 2𝑘-1 𝜂 𝑖 𝛼 1 … 𝜂 𝑖 𝛼 2𝑘-1 𝜙 𝑖 (∅) = 1 𝜆 𝑖 (∅) 𝜆 𝑖 (𝑣 𝑖 𝛼 1∶2𝑘-1 ) = 𝜂 𝑖 𝛼 1 … 𝜂 𝑖 𝛼 2𝑘-1 𝜙 𝑣 𝑖 𝛼 1∶2𝑘-1 𝑖 (𝑥) = 1 2(2𝑘 -1)! × 𝑎 𝛼 1 (𝑥) 𝜂 𝑖 𝛼 1 … 𝑎 𝛼 2𝑘-1 (𝑥) 𝜂 𝑖 𝛼 2𝑘-1 .
Remark 11. It is well-known that the Exponential Hawkes process may explode in finite time with probability non null [5] (that is, produce an infinite number of jumps within a finite time interval). To perform the Perfect Simulation, in this case, we shall consider the process having empty past before time 0, in other words, the process starting from 0. However, the difficulty lies on the fact that we need to check that the process stays in 𝒴 ∞ . To avoid this difficulty, we shall not consider the Exponential Hawkes process or other potentially explosive processes in the following section. The discussion of potentially explosive process will be the main object of upcoming work, as we mentioned earlier.

�.� Modified Perfect Simulation algorithm

In the following, we present a Perfect Simulation algorithm that simulates the process 𝑍 𝑖 on an interval [0, 𝑡𝑚𝑎𝑥], for some fixed 𝑡𝑚𝑎𝑥 > 0 in the stationary regime. Our algorithm is a modification of the method described in [19]. The procedure consists of backward and forward steps. In the backward steps, thanks to the Kalikow decomposition, we create a set of ancestors, which is a list of all the points that might influence the point under consideration. On the other hand, in the forward steps, where we go forward in time, by using the thinning method [17] we give the decision (keep or reject) to each visited point based on its neighborhood until the state of all considered points is decided. Further discussion can be found in [10,14,19].

For this algorithm to work, we introduce a subspace 𝒴 as follows:

𝒴 = {𝑥 ∈ 𝒳 such that ∀𝑣, 𝑖 ∶ 𝜙 𝑣 𝑖 (𝑥) ≤ Γ 𝑖 },
where Γ 𝑖 is a positive constant. We assume that the process generated by the algorithm stays in 𝒴 almost surely. For example, in the case of age dependent Hawkes processes in Section 3.4.2, we show in Proposition 9 below that Δ 𝑖 𝑘 (𝑥) are bounded. Thus, with an adequate choice of weights (𝜆 𝑖 (.)) 𝑖∈I and (Γ 𝑖 ) 𝑖∈I , the algorithm remains in 𝒴 almost surely. Moreover, for any 𝑥 ∈ 𝒴, we have

𝜙 𝑖 (𝑥) = 𝜆 𝑖 (∅)𝜙 ∅ 𝑖 + ∑ ∅≠𝑣∈V 𝑖 𝜆 𝑖 (𝑣)𝜙 𝑣 𝑖 (𝑥) ≤ Γ 𝑖 ,
which means that the intensity is bounded.

Furthermore, for any 𝑡 and 𝑥 𝑡 ∈ 𝒴 𝑡 , by the time homogeneity assumption,

𝜙 𝑖,𝑡 (𝑥 𝑡 ) = 𝜙 𝑖 (𝑥 ←𝑡 𝑡 ) ≤ Γ 𝑖 .
�.�.� Backward procedure Initial step. Fix 𝑖, set the initial time to be 0.

Step 1. Move to the first possible jump 𝑇 of 𝑍 𝑖 after 0 by taking

𝑇 ← 0 + 𝐸𝑥𝑝(Γ 𝑖 ).
Step 2. Recall that V 𝑖 is the neighborhood family associated to index 𝑖. Independently of anything else, pick a random neighborhood 𝑉 𝑖,𝑇 of (𝑖, 𝑇 ) according to the distribution (𝜆 𝑖 (𝑣)) Assume that 𝑉 𝑖,𝑇 = 𝑣 →𝑇 and define the projection to the second coordinate of 𝑣 by

𝜋 𝑗 (𝑣) ∶= {𝑡 ∈ R|(𝑗, 𝑡) ∈ 𝑣},
for any 𝑗 ∈ I. Notice that if for some 𝑗, (𝑗, 𝑡) ∉ 𝑣 for all 𝑡, then 𝜋 𝑗 (𝑣) = ∅.

Simulate Poisson processes in 𝑣 →𝑇 , that is for each 𝑗 ∈ I, we simulate a Poisson process Π 𝑗 with intensity Γ 𝑗 on 𝜋 𝑗 (𝑣 →𝑇 ). Put these points in 𝒞 1 𝑖,𝑇 , call it the first set of ancestors of (𝑖, 𝑇 ):

𝒞 1 𝑖,𝑇 = ⋃ 𝑗∈I { (𝑗, 𝑡)|𝑡 ∈ Π 𝑗( 𝜋 𝑗 (𝑣 →𝑇 ) )} .
Step 3. Recursively, we define the 𝑛 𝑡ℎ set of ancestors of (𝑖, 𝑇 ),

𝒞 𝑛 𝑖,𝑇 = ⋃ (𝑗,𝑠)∈𝒞 𝑛-1 𝑖,𝑇 𝒞 1 𝑗,𝑠 ⧵ ( 𝒞 1 𝑖,𝑇 ∪ … ∪ 𝒞 𝑛-1 𝑖,𝑇 ) ,
and for each (𝑗, 𝑡 ′ ) ∈ 𝒞 𝑛 𝑖,𝑇 , we perform Step 1 and Step 2. The backward scheme stops when 𝒞 𝑛 𝑖,𝑇 = ∅. We denote

𝑁 𝑖,𝑇 = inf{𝑛 ∶ 𝒞 𝑛 𝑖,𝑇 = ∅}.
The genealogy of (𝑖, 𝑇 ) is given by 𝒞 𝑖,𝑇 = ∪ 𝑁 𝑖,𝑇 𝑘=1 𝒞 𝑘 𝑖,𝑇 .

Remark 12.

Let us emphasize that 𝜆 𝑖 (∅) does not need to be strictly positive in the present Kalikow decomposition. This means that there is a chance that at every step of the Backward steps, we need to simulate a Poisson process in a non empty neighborhood. However, if there is no point simulated in these intervals, then we do nothing in the next step. Hence, the Backwards steps end. This is one of the main advantage of this Kalikow decomposition with respect to [10,14].

Remark 13. Notice that, comparing to the original Backward procedure [10,14], here, we start by simulating a point in the "future" by adding an exponential random variable to the initial time. Unfortunately, due to this additional step, the algorithm always requires that the intensity is bounded in 𝒴. This is one of the main drawback of this new algorithm.

�.�.� Forward procedure

We now attach to each point in 𝒞 𝑖,𝑇 a random variable 𝜒 whose value is either 0 or 1, where 1 means that this point is accepted.

For any point (𝑖, 𝑡) in the Backward steps, we know its neighborhood which contains all the points that might influence the state of (𝑖, 𝑡), but we do not know yet the state 𝜒 of the points in the neighborhood.

We start with the point (𝑗, 𝑠) ∈ 𝒞 𝑖,𝑇 which is the smallest in time, so that its associated neighborhood is either empty (𝑣 = ∅) or non empty but without any point of the Poisson process in it.

Step Remark 14. When implementing the algorithm on a computer, it is worth noticing that, whenever we simulate in an interval that intersects with previously simulated intervals, we will not simulate in the intersecting parts. In particular, no additional points is simulated in {𝑖} × (0, 𝑇 ) where (𝑖, 𝑇 ) is the first simulated point.

�.�.� Do we construct the right intensity?

For any 𝑖 ∈ I, we denote the arrival times of Π 𝑖 in Step 2 of the Backward steps by (𝜏 𝑖 𝑛 ) 𝑛∈Z , with 𝜏 𝑖 1 being the first positive time. Indeed even if we have simulated it on the randomly generated neighborhood, since the intensity is always the same, we can assume that all these points come only from one single Poisson process.

As in Step 1 of the Forward steps, we attach to each point of Π 𝑖 a stochastic mark 𝜒 given by,

𝜒 𝑖 𝑛 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if 𝜏 𝑖 𝑛 is accepted, 0 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒.
In addition, for any 𝑖 ∈ I, define 𝒵 𝑖 = (𝜏 𝑖 𝑛 , 𝜒 𝑖 𝑛 ) 𝑛∈Z an 𝐸-marked point process with 𝐸 = {0; 1}. In particular, following the notation in Chapter VIII of [3], for any 𝑖 ∈ I, let Moreover note that (𝒵 𝑖 𝑡 (1)) 𝑡∈R is the counting process associated to the accepted points of the algorithm. With these notations, we can prove the following result. This proof is already done in [19] but we add it here for sake of completeness.

𝑑𝒵 𝑖 𝑡 (𝑚𝑎𝑟𝑘) = ∑ 𝑛∈Z 1 𝜒 𝑖 𝑛 =𝑚𝑎𝑟𝑘 𝛿 𝜏 𝑖 𝑛 (𝑑𝑡) for 𝑚𝑎𝑟𝑘 ∈ 𝐸, ℱ 𝒵 𝑡 -= ⋁ 𝑖∈I 𝜎(𝒵 𝑖 𝑠 (0), 𝒵 𝑖 𝑠 (
E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝑑𝒵 𝑖 𝑡 (1)
⎞ ⎟ ⎟ ⎠ = ∞ ∑ 𝑛=1 E ( 𝐶 𝜏 𝑖 𝑛 1 𝜒 𝑖 𝑛 =1 ) .
Note that by Theorem T35 at Appendix A1 of [3], any point 𝑇 should be understood as a stopping time, and that by Theorem T30 at Appendix A2 of [3],

ℱ 𝒵 𝑇 -= ⋁ 𝑗∈I 𝜎{𝜏 𝑗 𝑚 , 𝜒 𝑗 𝑚 such that 𝜏 𝑗 𝑚 < 𝑇 , 𝑇 }.
Therefore,

E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝑑𝒵 𝑖 𝑡 (1) ⎞ ⎟ ⎟ ⎠ = ∞ ∑ 𝑛=1 E ( 𝐶 𝜏 𝑖 𝑛 E(1 𝜒 𝑖 𝑛 =1 |ℱ 𝒵 𝜏 𝑖 𝑛 -, 𝑉 𝑖 𝑛 ) ) = ∞ ∑ 𝑛=1 E ⎛ ⎜ ⎜ ⎜ ⎝ 𝐶 𝜏 𝑖 𝑛 𝜙 𝑉 𝑖 𝑛 𝑖,𝜏 𝑖 𝑛 (𝑥 𝜏 𝑖 𝑛 ) Γ 𝑖 ⎞ ⎟ ⎟ ⎟ ⎠ ,
where we denote 𝑉 𝑖 𝑛 for the neighborhood of 𝜏 𝑖 𝑛 . Let us now integrate with respect to the choice 𝑉 𝑖 𝑛 , which is independent of anything else.

E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝑑𝒵 𝑖 𝑡 (1) ⎞ ⎟ ⎟ ⎠ = ∞ ∑ 𝑛=1 E ⎛ ⎜ ⎜ ⎜ ⎝ 𝐶 𝜏 𝑖 𝑛 𝜆 𝑖 (∅)𝜙 ∅ 𝑖,𝜏 𝑖 𝑛 + ∑ 𝑣∈V →𝜏 𝑖 𝑛 ,𝑣≠∅ 𝜆 𝑖 (𝑣) × 𝜙 𝑣 𝑖,𝜏 𝑖 𝑛 (𝑥 𝜏 𝑖 𝑛 ) Γ 𝑖 ⎞ ⎟ ⎟ ⎟ ⎠ = E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝜙 𝑖,𝑡 (𝑥 𝑡 ) Γ 𝑖 𝑑Π 𝑖 (𝑡) ⎞ ⎟ ⎟ ⎠ .
Since Π 𝑖 is a Poisson process with respect to ℱ 𝒵 𝑡 with intensity Γ 𝑖 , and

𝐶 𝑡 𝜙 𝑖,𝑡 (𝑥 𝑡 ) Γ 𝑖 is ℱ 𝒵 𝑡-measurable, we finally have that E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝑑𝒵 𝑖 𝑡 (1) ⎞ ⎟ ⎟ ⎠ = E ⎛ ⎜ ⎜ ⎝ ∞ ∫ 0 𝐶 𝑡 𝜙 𝑖,𝑡 (𝑥 𝑡 )𝑑𝑡 ⎞ ⎟ ⎟ ⎠
, which ends the proof.

�.�.� Why does the Backward steps end?

We construct a tree with root (𝑖, 𝑇 ). For each point (𝑗 𝑇 ′ , 𝑇 ′ ) in the tree, the points which are simulated in 𝑉 𝑗 𝑇 ′,𝑇 ′ (Step 2 of the Backward steps) define the children of (𝑗 𝑇 ′ , 𝑇 ′ ) in the tree. This forms the tree T.

Let us now build a tree C with root (𝑖, 𝑇 ) (that includes the previous tree) by mimicking the procedure in the Backward steps, except that we simulate -independently on anything else -on the whole neighborhood even if it has a part that intersects with previous neighborhoods (if they exist) (Step 3 of the Backward steps). By doing so, we make the number of children larger but at each node, they are independent of anything else.

If the tree C goes extinct, then so does the tree T , and the backward part of the algorithm terminates.

To formulate a sufficient criterion that implies the almost sure extinction of C, let us denote the product measure 𝑃 on I × R, defined on the product sets that generate the Borel subsets of I × R, as follows.

𝑃 (𝐽 × 𝐴) ∶= ∑ 𝑗∈𝐽 Γ 𝑗 𝜇(𝐴)
for any 𝐽 ⊂ I, where 𝐴 is a Borel subset of R and where 𝜇 is the Lebesgue measure. The following proposition is already proved in [19], but without precisely defining 𝑃. This implies that ∑ 𝑘≥1 𝑃 (𝑣 𝑖 𝑘 )𝜆 𝑖 (𝑣 𝑖 𝑘 ) is the mean number of children issued from one point of type 𝑖. Then, the condition (3.5.1) says that the number children in each step should be less than one in average, for the tree to go extinct almost surely. That is a very classical result in Branching process [1].

�.�.� The complexity of the algorithm

In this section, we study the effect of (𝜆 𝑖 (𝑣 𝑖 𝑘 )) 𝑖,𝑘 on the number of points simulated by our algorithm. Until the end of this section, we suppose that Assumption 6. The index set I is finite, namely

|I| = 𝑁 < ∞.
Let us mention several notations that will be useful in the sequel. We denote 𝑒 𝑖 the 𝑖-th unit vector of R 𝑁 , 1 is the vector (1, 1, … , 1) 𝑇 and 𝜇 always stands for the Lebesgue measure. Finally, by a positive vector, we mean that all its components are positive.

Coming back to our problem, since the tree C dominates T, to give an upper bound of the number of points of T, it is sufficient to study the number of points of C. Recall that the root of the tree C is of type 𝑖 ∈ I. For 𝑛 ≥ 1, write 𝐾 𝑖 (𝑛) for the vector containing the numbers of ancestors in the 𝑛 𝑡ℎ set of ancestors of a single point of index 𝑖. We consider In addition, we denote 𝑋 𝑖 = (𝑋 1 𝑖 , 𝑋 2 𝑖 , … , 𝑋 𝑁 𝑖 ) 𝑇 and for any 𝜃 ∈ R 𝑁 , we consider the log-Laplace transform of 𝑋 𝑖 :

𝐾 𝑖 (𝑛) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ number of
𝜙 𝑖 (𝜃) ∶= log E 𝑖 ( 𝑒 𝜃 𝑇 𝑋 𝑖 )

where we denote P 𝑖 the law of a random tree C , whose root is of type 𝑖 and E 𝑖 the corresponding expectation.

For 𝑛 ≥ 2, given the population in the first set of ancestors 𝑋 𝑖 , we have the following relationship between the (𝑛 -1) 𝑡ℎ set of ancestors and the 𝑛 𝑡ℎ set of ancestors in C ,

𝐾 𝑖 (𝑛) = ∑ 𝑗∈I 𝑋 𝑗 𝑖 ∑ 𝑝=1 𝐾 (𝑝) 𝑗 (𝑛 -1),
where 𝐾 (𝑝) 𝑗 (𝑛 -1) is the vector of the number of ancestors in the (𝑛 -1) 𝑡ℎ set of ancestors issued from the 𝑝-th point of type 𝑗 in the first generation. In addition, for 𝑝 = 1, …, 𝑋 𝑗 𝑖 , we have that the 𝐾 (𝑝) 𝑗 (𝑛 -1)'s are independent copies of 𝐾 𝑗 (𝑛 -1). Note that this equation is trivial for 𝑛 = 1.

Moreover, we consider

𝑊 𝑖 (𝑛) = 𝑛 ∑ 𝑘=0 𝐾 𝑖 (𝑘)
the total number of ancestors in the first 𝑛 set of ancestors.

The log Laplace transform associated to the random vector 𝑊 𝑖 (𝑛) is given by:

Φ (𝑛) 𝑖 (𝜃) ∶= log E 𝑖 ( 𝑒 𝜃 𝑇 𝑊 𝑖 (𝑛)
) .

Therefore, the total number of points simulated of the tree C with initial point of type 𝑖 is

E ( 𝑊 𝑖 (𝑛) ) = E ( 𝑛 ∑ 𝑘=0 𝐾 𝑖 (𝑘) ) = ( 𝑛 ∑ 𝑘=0 𝑀 𝑘 ) 𝑒 𝑖 .
It leads us to conclude that,

E(𝑊 𝑖 ) = 1 𝑇 ( ∞ ∑ 𝑘=0 𝑀 𝑘 ) 𝑒 𝑖 ,
where the matrix 𝑀 is defined in (3.5.2).

To conclude, if we think the complexity of the algorithm in terms of the number of simulated points, 𝑊 𝑖 is a good bound, and by (3.5.3) we see that we can grasp the complexity of the algorithm by studying (3.5.3).

Proof. First, we prove that 𝜙 𝑖 (𝜃) is well defined for all 𝜃 ∈ 𝐵(0, 𝑟). Indeed, since (𝑋 𝑗 𝑖 ) 𝑗=1,…,𝑁 are independent, we have

𝜙 𝑖 (𝜃) = log E 𝑖 ( ∏ 𝑗 exp (𝜃 𝑗 𝑋 𝑗 𝑖 ) ) = ∑ 𝑗 log E 𝑖 ( exp (𝜃 𝑗 𝑋 𝑗 𝑖 ) ) .
Moreover, by conditioning on the first random neighborhood of type 𝑖, 𝑉 𝑖 , we have

𝑋 𝑗 𝑖 |𝑉 𝑖 = 𝑣 𝑖 𝑘 ∼ 𝒫 ( Γ 𝑗 𝜇(𝜋 𝑗 (𝑣 𝑖 𝑘 )) ) .
Therefore, we conclude that

E 𝑖 ( exp (𝜃 𝑗 𝑋 𝑗 𝑖 )|𝑉 𝑖 = 𝑣 𝑖 𝑘) = exp [ (𝑒 𝜃 𝑗 -1)Γ 𝑗 𝜇 ( 𝜋 𝑗 (𝑣 𝑖 𝑘 ) )] .
Finally, we obtain that

𝜙 𝑖 (𝜃) = ∑ 𝑗 log ( ∑ 𝑘 𝜆 𝑖 (𝑣 𝑖 𝑘 ) exp [ (𝑒 𝜃 𝑗 -1)Γ 𝑗 𝜇 ( 𝜋 𝑗 (𝑣 𝑖 𝑘 ) )] ) < ∞.
In the following, we prove that Φ (𝑛) 𝑖 (𝜃) satisfies the following recursion

Φ (𝑛) 𝑖 (𝜃) = 𝜃 𝑇 𝐾 𝑖 (0) + 𝜙 𝑖 (Φ (𝑛-1) (𝜃)).
Indeed, by definition of 𝑊 𝑖 (𝑛) we have

E 𝑖 ( 𝑒 𝜃 𝑇 𝑊 𝑖 (𝑛) ) = 𝑒 𝜃 𝑇 𝐾 𝑖 (0) E 𝑖 ( 𝑒 𝜃 𝑇 ∑ 𝑛 𝑘=1 𝐾 𝑖 (𝑘)
) .

In addition, from the definition of 𝐾 𝑖 (𝑘) we obtain,

E 𝑖 ( 𝑒 𝜃 𝑇 ∑ 𝑛 𝑘=1 𝐾 𝑖 (𝑘) ) = E 𝑖 ( 𝑒 𝜃 𝑇 ∑ 𝑛 𝑘=1 ∑ 𝑁 𝑗=1 ∑ 𝑋 𝑗 𝑖 𝑝=1 𝐾 (𝑝) 𝑗 (𝑘-1) ) = E 𝑖 ( 𝑁 ∏ 𝑗=1 𝑒 𝜃 𝑇 ∑ 𝑛 𝑘=1 ∑ 𝑋 𝑗 𝑖 𝑝=1 𝐾 (𝑝) 𝑗 (𝑘-1) ) = E 𝑖 [ E ( 𝑁 ∏ 𝑗=1 𝑒 𝜃 𝑇 ∑ 𝑛 𝑘=1 ∑ 𝑋 𝑗 𝑖 𝑝=1 𝐾 (𝑝) 𝑗 (𝑘-1) |𝑋 𝑖 )] .
Since conditioning on 𝑋 𝑖 , 𝐾

𝑗 (𝑘 -1)'s are independent, we have that

E 𝑖 ( 𝑒 𝜃 𝑇 ∑ 𝑛 𝑘=1 𝐾 𝑖 (𝑘) ) = E 𝑖 [ 𝑁 ∏ 𝑗=1 E ( 𝑒 𝜃 𝑇 ∑ 𝑛 𝑘=1 ∑ 𝑋 𝑗 𝑖 𝑝=1 𝐾 (𝑝) 𝑗 (𝑘-1) |𝑋 𝑗 𝑖 )] = E 𝑖 ⎡ ⎢ ⎢ ⎣ 𝑁 ∏ 𝑗=1 𝑋 𝑗 𝑖 ∏ 𝑝=1 E ( 𝑒 𝜃 𝑇 ∑ 𝑛 𝑘=1 𝐾 (𝑝) 𝑗 (𝑘-1) | 𝑋 𝑗 𝑖 ) ⎤ ⎥ ⎥ ⎦ = E 𝑖 ⎡ ⎢ ⎢ ⎣ 𝑁 ∏ 𝑗=1 𝑋 𝑗 𝑖 ∏ 𝑝=1 E ( 𝑒 𝜃 𝑇 ∑ 𝑛 𝑘=1 𝐾 (𝑝) 𝑗 (𝑘-1) ) ⎤ ⎥ ⎥ ⎦ = E 𝑖 [ 𝑁 ∏ 𝑗=1 ( E ( 𝑒 𝜃 𝑇 𝑊 𝑗 (𝑛-1)
))

𝑋 𝑗 𝑖 ] = E 𝑖 [ 𝑁 ∏ 𝑗=1 𝑒 Φ (𝑛-1) 𝑗 (𝜃)𝑋 𝑗 𝑖 ] = E 𝑖 [ 𝑒 Φ (𝑛-1) (𝜃) 𝑇 𝑋 𝑖 ] = 𝑒 𝜙 𝑖 (Φ (𝑛-1) (𝜃)) .
Finally, we conclude that

Φ (𝑛) 𝑖 (𝜃) = 𝜃 𝑇 𝐾 𝑖 (0) + 𝜙 𝑖 (Φ (𝑛-1) (𝜃)).
This equation holds for every 𝑖, therefore we have

Φ (𝑛) (𝜃) = 𝜃 + 𝜙(Φ (𝑛-1) (𝜃)). (3.5.4)
Let us consider the column sums of the matrix 𝑀:

∑ 𝑗∈I 𝑀 𝑗𝑖 = ∑ 𝑗∈I E 𝑖 (𝑋 𝑗 𝑖 ) = ∑ 𝑘≥1 𝑃 (𝑣 𝑖 𝑘 )𝜆 𝑖 (𝑣 𝑖 𝑘 ). Hence, ‖𝑀‖ 1 = sup ‖𝑥‖ 1 ≤1 {‖𝑀𝑥‖ 1 } = sup 𝑖 ∑ 𝑗∈I |𝑀 𝑗𝑖 | < 1
where ‖.‖ 1 is the induced norm for matrix on R 𝑁×𝑁 . Therefore, ‖𝐷𝜙(0)‖ 1 ≤ 𝐶 < 1. Moreover the norm is continuous and 𝐷𝜙(𝑠) is likewise, there is a 𝑟 > 0 such that, for

||𝑠|| 1 ≤ 𝑟, ‖𝐷𝜙(𝑠)‖ 1 ≤ 𝐶 < 1.
Hence, 𝜙(𝑠) is Lipschitz continuous in the ball 𝐵(0, 𝑟) and moreover 𝜙(0) = 0, which implies that

‖𝜙(𝑠)‖ 1 ≤ 𝐶‖𝑠‖ 1 for ‖𝑠‖ 1 ≤ 𝑟.
Moreover, take 𝜃 such that ‖𝜃‖ 1 1 -𝐶 ≤ 𝑟, hence ‖𝜃‖ 1 ≤ 𝑟. By induction we can show that

‖Φ (𝑛) (𝜃)‖ 1 ≤ ‖𝜃‖ 1 (1 + 𝐶 + … + 𝐶 𝑛 ) ≤ 𝑟 < ∞.
In addition, (𝑊 𝑖 (𝑛)) 𝑛≥1 is a positive, increasing vector sequence. Here, for any 𝑢, 𝑣 ∈ R 𝑁 , we say 𝑢 ≥ 𝑣 if (𝑢 -𝑣) is a positive vector.

For any 𝜃 such that ‖𝜃‖ 1 1 -𝐶 ≤ 𝑟, by using the theorem of monotone convergence, we have

Φ (𝑛) (𝜃) → Φ(𝜃)
when 𝑛 → ∞. Moreover, ‖Φ(𝜃)‖ 1 < ∞, and passing to the limit 𝑛 → ∞ in equation (3.5.4), we conclude that

Φ(𝜃) = 𝜃 + 𝜙(Φ(𝜃)).
In particular, choosing 𝜃 = 𝜗1 ∈ 𝐵(0, 𝑟) with 𝜗 ∈ R, we have

E 𝑖( 𝑒 𝜗𝑊 𝑖 ) < ∞,
where 𝑊 𝑖 = 1 𝑇 𝑊 𝑖 (∞) is the total number of point of C .

The last point is concluded by using Markov's inequality, that ends the proof.

�.�.� Efficiency of the algorithm and discussion of the choice of the weights on a particular example Finally, to illustrate the effect of 𝜂 to the complexity of Perfect Simulation algorithm, we consider the age dependent Hawkes process with hard refractory period of length 𝛿. We consider the setting of Proposition 4 and assume moreover that 𝜓 𝑖 (.) is an 𝐿-Lipschitz function, where 𝐿 is the Lipschitz constant. In this case we obtain an explicit upper bound for Δ 𝑖 𝑘 (𝑥), for any 𝑥 ∈ 𝒳 >𝛿 . Proposition 9. For any 𝑘 ≥ 1, we have

Δ 𝑖 𝑘 (𝑥) ≤ 𝐿 × [ ∑ 𝑗∈𝑉 𝑖 (𝑘)⧵𝑉 𝑖 (𝑘-1) ∑ 0≤𝑚<𝑘 ℎ 𝑗𝑖 (𝑚𝛿) + ∑ 𝑗∈𝑉 𝑖 (𝑘-1)
ℎ 𝑗𝑖 ((𝑘 -1)𝛿) ] .

To prove this proposition we use the following lemma, which is a particular case of Lemma 2.4 in [20].

Lemma 1. For any 0 ≤ 𝑘 < 𝑙, 𝑗 ∈ I, 𝑧 ∈ 𝒳 >𝛿 , 𝑡 ∈ R we have

∫ 𝑡-𝑘𝛿 𝑡-𝑙𝛿 ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑧 𝑗 𝑠 ≤ ∑ 𝑘≤𝑚<𝑙 ℎ 𝑗𝑖 (𝑚𝛿).
Proof of Lemma 1. For any 𝑗, fixed 𝜖 > 0, we have

∫ [𝑡-𝑙𝛿,𝑡-𝑘𝛿-𝜖] ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑧 𝑗 𝑠 = ∑ 𝑘≤𝑚<𝑙 ∫ [𝑡-(𝑚+1)𝛿-𝜖,𝑡-𝑚𝛿-𝜖] ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑧 𝑗 𝑠 ≤ ∑ 𝑘≤𝑚<𝑙 ℎ 𝑗𝑖 (𝑚𝛿 + 𝜖)
by Assumption 3 and the fact that there is at most one jump in the interval of length 𝛿. Therefore,

∫ [𝑡-𝑙𝛿,𝑡-𝑘𝛿) ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑧 𝑗 𝑠 = lim 𝜖↓0 ∫ [𝑡-𝑙𝛿,𝑡-𝑘𝛿-𝜖] ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑧 𝑗 𝑠 ≤ lim 𝜖↓0 ∑ 𝑘≤𝑚<𝑙 ℎ 𝑗𝑖 (𝑚𝛿 + 𝜖) ≤ ∑ 𝑘≤𝑚<𝑙 lim 𝜖↓0 ℎ 𝑗𝑖 (𝑚𝛿 + 𝜖) ≤ ∑ 𝑘≤𝑚<𝑙 ℎ 𝑗𝑖 (𝑚𝛿),
by using the theorem of monotone convergence and the fact that ℎ 𝑗𝑖 (.) is a decreasing function according to Assumption 3. This completes the proof of Lemma 1.

Proof of Proposition 9. By (ii) of Assumption 3, we have

Δ 𝑘 𝑖 (𝑥) ≤ 𝐿 × [ ∑ 𝑗∈𝑉 𝑖 (𝑘)⧵𝑉 𝑖 (𝑘-1) ∫ 0 -𝑘𝛿 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 + ∑ 𝑗∈𝑉 𝑖 (𝑘-1) ∫ -𝑘𝛿+𝛿 -𝑘𝛿 ℎ 𝑗𝑖 (-𝑠)𝑑𝑥 𝑗 𝑠 ] .
Applying Lemma 1 we conclude that

Δ 𝑖 𝑘 (𝑥) ≤ 𝐿 × [ ∑ 𝑗∈𝑉 𝑖 (𝑘)⧵𝑉 𝑖 (𝑘-1) ∑ 0≤𝑚<𝑘 ℎ 𝑗𝑖 (𝑚𝛿) + ∑ 𝑗∈𝑉 𝑖 (𝑘-1)
ℎ 𝑗𝑖 ((𝑘 -1)𝛿) ]

which ends the proof.

Remark 17. In addition, from (i) of Assumption 3, we have

∑ 𝑗∈𝑉 𝑖 (𝑘)⧵𝑉 𝑖 (𝑘-1) ∑ 0≤𝑚<𝑘 ℎ 𝑗𝑖 (𝑚𝛿) + ∑ 𝑗∈𝑉 𝑖 (𝑘-1) ℎ 𝑗𝑖 ((𝑘 -1)𝛿) → 0 as 𝑘 → ∞. Therefore Δ 𝑖 𝑘 (𝑥) → 0 when 𝑘 → ∞.
In what follows, to simplify the computation, we consider that I = Z, and for all 𝑖, we set 1. ℎ 𝑗𝑖 (𝑡) = 𝛽 𝑗𝑖 exp (-𝛼𝑡) where 𝛽 𝑗𝑖 , 𝛼 are positive constants for all 𝑗, 𝑖. In addition, we take 𝛼 = 1 𝛿 and 𝛽 𝑗𝑖 = 1 2|𝑗 -𝑖| 𝛾 for 𝑗 ≠ 𝑖 with a positive number 𝛾 and 𝛽 𝑖𝑖 = 1.

𝑉

𝑖 (0) = ∅, 𝑉 𝑖 (1) = {𝑖}, … , 𝑉 𝑖 (𝑘) = {𝑖 -𝑘 + 1, … , 𝑖, 𝑖 + 1, … , 𝑖 + 𝑘 -1} ∀𝑘 ≥ 2. 3. 𝜂 𝑖 𝑘 = 𝜂 𝑘 = 𝑐 𝜂 1 𝑘 𝑝 , ∀𝑘 ≥ 1,
where 𝑝 is a positive constant and 𝑐 𝜂 is a normalization constant.

From Proposition 9, we can choose

Γ 𝑖 = sup { 𝐿 𝑐 𝜂 , sup 𝑘≥2 [ 𝐿𝑘 𝑝 𝑐 𝜂 ( 1 -𝑒 -𝑘 1 -𝑒 -1 ∑ 𝑗∈{𝑖-𝑘+1,𝑖+𝑘-1} 𝛽 𝑗𝑖 + 𝑒 -𝑘+1 ∑ 𝑗∈𝑉 𝑖 (𝑘-1) 𝛽 𝑗𝑖 )]} = sup { 𝐿 𝑐 𝜂 , sup 𝑘≥2 [ 𝐿𝑘 𝑝 𝑐 𝜂 ( 1 -𝑒 -𝑘 (𝑘 -1) 𝛾 (1 -𝑒 -1 ) + 𝑒 -𝑘+1 ( 1 + 𝑘-2 ∑ 𝑚=1 1 𝑚 𝛾 ))]} ∶= Γ.
Since we define 𝜂 𝑖 𝑘 are independent of 𝑖, so we can consider that all points have the same type. Instead of comparing the Backward Steps to a multitype branching process, we may therefore put 𝐾(𝑘) ∶= 𝐾 𝑖 (𝑘) ∶= ∑ 𝑗 𝐾 𝑗 𝑖 (𝑘) and 𝑊 (𝑛) ∶= 𝑊 𝑖 (𝑛) ∶= ∑ 𝑗 𝑊 𝑗 𝑖 (𝑛) such that 𝐾(𝑘) is the total number of ancestors in the 𝑘 𝑡ℎ set of ancestors and 𝑊 (𝑛) is the total number of ancestors in the first 𝑛 set of ancestors. Hence, all the vector 𝐾(𝑘) ∶= 𝐾 𝑖 (𝑘) and 𝑊 (𝑛) = 𝑊 𝑖 (𝑛) in Proposition 8 now are numbers.

Moreover, 𝑣 𝑖 𝑘 = 𝑉 𝑖 (𝑘) × [-𝑘𝛿, 0) and we set

𝜁 = ∑ 𝑘 ′ ≥1 𝑃 (𝑣 𝑖 𝑘 ′ )𝜆(𝑣 𝑖 𝑘 ′ ) = ∑ 𝑘 ′ ≥1 𝜂 𝑘 ′ 𝑘 ′ 𝛿 ∑ 𝑗∈𝑉 𝑖 (𝑘 ′ ) Γ 𝑗 = ∑ 𝑘 ′ ≥1 𝑐 𝜂 1 𝑘 ′𝑝 𝑘 ′ 𝛿(2𝑘 ′ -1)Γ ∶= 𝑓 (𝑝).
By a classical well-known result in Branching processes [1] , we have for 𝑘 ≥ 1

E(𝐾(𝑘)) = 𝜁 𝑘 .
Therefore, the total expected number of simulated points 𝑊 ∶= 𝑊 (∞) is given by Once all these condition are fulfilled, we choose 𝑝 such that it minimizes 𝑓 (𝑝), so that the total number of simulated points is minimal.

E(𝑊 ) = 1 1 -𝜁 .
(𝑎) ∶ ∑ 𝑘 𝜂 𝑘 = 1 ⇒ 𝑝 > 1 (𝑏) ∶ Γ < ∞ ⇒ 𝑝 < 𝛾 (𝑐) ∶ ∑ ∫ ℎ 𝑗𝑖 < ∞ ⇒ 𝛾 > 1 (𝑑) ∶ 𝜁 < 1 ⇒ 𝑝 -2 > 1
Finally, we conclude that 3 < 𝑝 < 𝛾. Moreover, we want to choose 𝑝 such that 𝑓 (𝑝) to be smallest possible, which means that we need to take 𝑝 largest possible. Hence, we take 𝑝 close to 𝛾.

�.� Conclusion

Prior research has investigated the Perfect Simulation based on a conditional Kalikow decomposition in continuous time point processes [14]. However, it is impossible to implement this study in practice. In the present study, we continue our work in [19], to extend the class of point processes which satisfies a new Kalikow decomposition. This decomposition plays a vital role to build a tractable algorithm in practice.

In addition, we also improve the results in [14] and [19] on the existence of the Kalikow decomposition.

Most notably, this is the first study to our knowledge to establish a Kalikow decomposition for a stochastic intensity in general context. Our results provide a general method to write a Kalikow decomposition for large class of point processes, for example: linear Hawkes processes, exponential Hawkes processes and including very complex Hawkes processes: age dependent Hawkes processes. However, some limitations are worth noting. Although, we succeed to write the Kalikow decomposition for a variety of Hawkes processes, we still have to restrict ourselves to a bounded intensity to implement the Perfect Simulation. This is due to the simulation of the first jump in the Backward Steps. Future work should focus on removing totally the upper bound of the intensity and extending the Perfect Simulation to the unbounded intensities.

�.� Introduction

Simulating neural networks has received a lot of attention in recent years, for example, there exist several grand projects such as the Human Brain Project in Europe, the Brain Mapping in Japan and the Brain Initiative in the United States. A complete simulated model will help reduce the cost of performing biological experiments and create an important premise before the actual experiments. In practice, there are many approaches to simulate neural networks, to cite but a few [1,15,23]. Here, we will focus on the point process approach, especially the simulation algorithm using the classical Ogata's thinning method [16]. However, it has been found to be too expensive for the huge network, see [14]. Instead, Mascart et al [14] considered a sparse network, that is a node is only connected to only few other nodes. This algorithm [14] can not apply to a complete graph with a huge number of nodes.

Another way to handle this problem is to combine the classical Ogata algorithm and Kalikow decomposition [13]. There has been extensive research regarding this combination for theoretical purposes. For example, in probability theory, for discrete-time processes, the authors of [4,5,6] developed a Perfect simulation algorithm by using Kalikow decomposition (see also Galves et al. [7,8]). However, although the numerous studies in discrete time, little attention has been paid to the continuous-time processes.

Up to our knowledge, there only exists one theoretical result, which was developed in [12]. Despite the excellent result in theory, their algorithm is not feasible in practice (see the discussion in [19] and also Chapter 2). Another notable limitation of [12] is the boundedness hypothesis of the intensities. Recently, Phi et al [19](see also Chapter 2 of this thesis) showed an interesting result by using a new type of Kalikow decomposition. More precisely, we have successfully simulated the activity of one neuron embedded in an infinite network without simulating the whole network. However, it rather focuses on simulating the activity of a part of neural networks than the whole network.

In this chapter, we will focus on simulating the neural networks with a large but not infinite number of neurons and for mathematical interest, we focus on the case where no deterministic upper bound for intensities exists at prior. This chapter is organized as follows. In the first section, we give the mathematical definitions and several useful notation. In Section 4.3, we present four algorithms, two are modifications of the classical Ogata's algorithm, two are the combination of Ogata's algorithm and Kalikow decomposition.

All the proposed algorithms in this section are written for point processes having deterministically upper bounded intensity. In Section 4.4, we provide a modification of Ogata's algorithm and a combination of Ogata and Kalikow decomposition when such a deterministic bound does not exist. In Section 4.5, several numerical results are presented.

�.� Mathematical definitions and notation

We consider 𝑁 point processes in the time interval [0, 𝑡 𝑚𝑎𝑥 ] and assume that no events occurred before the • either a Poisson process with intensity 𝑀 𝑖 if Assumption 7 is satisfied,

• or a more intricate dominating process Π 𝑖 of stochastic intensity Λ 𝑖 if Assumption 8 is satisfied.

In both cases, these points are going to be accepted or rejected according to various marks or uniform variables that are drawn at each step. This is the reason why we define a 𝑝𝑜𝑖𝑛𝑡 (written in italic) as a vector whose components are given in the following table . The number of the components can be extended if needed. For example, in the second part, we need to draw several uniform random variables. Also, for the choice of the neighborhood in

A produced 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖

𝑇 is said to be known if the thinning mark 𝑋 𝑇 ≠ 𝑛.𝑎 . More precisely 𝑆 𝑖 𝑇 is said to be accepted if 𝑋 𝑇 = 1 and rejected if 𝑋 𝑇 = 0. For instance, in the Ogata's algorithm for the deterministic bounded intensity, 𝑋 𝑇 is defined as the following

𝑋 𝑇 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if 𝑈 𝑇 ≤ P(𝑋 𝑇 = 1|ℱ 𝑇 ), 0 otherwise. (4.2.3)
where P(𝑋 𝑇 = 1|ℱ 𝑇 ) is the ratio between the intensity 𝜙 𝑖 (𝑇 |ℱ 𝑇 ) and 𝑀 𝑖 .

From now on, to simplify writing and to use notation coming from neural networks, we say 𝑠𝑝𝑖𝑘𝑒 for an accepted 𝑝𝑜𝑖𝑛𝑡.

�.�.� Neuron structure

With the thinning method, several processes are taken into account to produce and accept a 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 . The analogy between a biological neuron and a point process reaches its limits. That is why we introduce a formal definition of a 𝑛𝑒𝑢𝑟𝑜𝑛 (written in italic). A 𝑛𝑒𝑢𝑟𝑜𝑛 is a mathematical structure that includes all the variables required to produce and accept the 𝑝𝑜𝑖𝑛𝑡𝑠 as well as the produced 𝑝𝑜𝑖𝑛𝑡𝑠 themselves (whether Remark 19. We define a relation " ⪯ " between two configurations of thinning marks on 𝒮 𝑇 : 𝒮 Compute the number of 𝑝𝑜𝑖𝑛𝑡𝑠 : 𝑛 𝑖 follows a 𝑃 𝑜𝑖𝑠𝑠𝑜𝑛(𝑀 𝑖 𝑡 𝑚𝑎𝑥 ). end for If the range of probability [P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|ℱ 𝑇 ), P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|ℱ 𝑇 )] is small enough, this allows us to determine the thinning mark.

𝑋 𝑇 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 if 𝑈 𝑇 ≤ P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|ℱ 𝑇 ) 0 if 𝑈 𝑇 ≥ P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|ℱ 𝑇 ) 𝑛.𝑎 otherwise. (4.3.12)

Overall description

To build a parallelized algorithm we introduce the structure 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 consisting on a subset of the set of the 𝑛𝑒𝑢𝑟𝑜𝑛𝑠. Let K be the set of index of 𝑁 𝐾 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠, K = {1, … , 𝑁 𝐾 }. The set of 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠 forms a partition of the set of 𝑛𝑒𝑢𝑟𝑜𝑛𝑠. This means that each 𝑛𝑒𝑢𝑟𝑜𝑛 belongs to exactly one 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛.

The algorithm computes sequentially each 𝑛𝑒𝑢𝑟𝑜𝑛 within the same population and the 𝑛𝑒𝑢𝑟𝑜𝑛𝑠 of two different 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠 are computed in parallel.

During the seeding step, the computation of each 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 is completely independent. Indeed this step consists in simulating independent Poisson processes for each 𝑛𝑒𝑢𝑟𝑜𝑛. 𝒮 𝑘 𝑝𝑜𝑝 , the set of produced 𝑝𝑜𝑖𝑛𝑡𝑠 for each 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑘 is sorted in an increasing 𝑇-order. Then to prepare for the next step, the complete set of produced 𝑝𝑜𝑖𝑛𝑡𝑠 is established as

𝒮 = ⋃ 𝑁 𝐾 1 𝒮 𝑘 𝑝𝑜𝑝 .
The thinning step is an alternation between parallel computation and sequential information sharing.

During a parallel phase, each 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑘 goes through the set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 𝑘 𝑝𝑜𝑝 trying to determine as many thinning marks as possible. Then, the 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠 synchronize, each one waiting for the others to finish their computation. During the sequential phase, the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 is updated with the thinning marks determined by each 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛. The alternation between the parallel and sequential phase is operated as long as there are still unknown 𝑝𝑜𝑖𝑛𝑡𝑠. for each 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑘 in parallel do for each 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑘 in sequence do

Pseudo-code

23:

Update the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 with the thinning marks Remark 24. The algorithm OtpP is only partially parallelized. In particular, there is an alternation between a parallel loop (line 15) and a sequential loop (line 22). There is also a sequential computation (line 13) that separates the two steps of the thinning.

Remark 25.

A crucial hypothesis that is needed for the algorithm OtpP works is the construction of P 𝑚𝑖𝑛 , P 𝑚𝑎𝑥 . In general, the algorithm OtpP works as soon the intensity is increasing with respect to the configuration of thinning marks, that is

𝜙(𝑇 |𝒮 𝑇 ) ≤ 𝜙(𝑇 |𝒮 ′ 𝑇 )
whenever 𝒮 𝑇 ⪯ 𝒮 ′ 𝑇 . In words, more accepted points more chances to emit a new spike.

The synchronization between the two steps seems to us unavoidable. On the other hand, it might be possible to avoid sequential computations within the second step. However, we have encountered difficulties with the use of the python multiprocessing package. When dynamic information was shared by parallel procedures, it generated synchronizations that we could not control. This is why, in the algorithm OtpP, the information managed by a 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 is always inaccessible to the information managed by another.

This of course has a cost, as we regularly interrupt parallel computation to consolidate and share information (line 22). Further discussion can be found in Section 4.5 of this chapter.

Remark 26.

A criticism that can be made about the algorithm partially parallelizing Ogata's thinning is the presence of the WHILE loop (line 14). As far as we know it is not possible to get rid of it.

However, it remains possible to reduce the number of iterations using Kalikow decomposition. Since the intensity using Kalikow decomposition depends on a smaller number of past event, the probability range [P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|ℱ 𝑇 ), P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|ℱ 𝑇 )] to accept 𝑝𝑜𝑖𝑛𝑡 (4.3.8) should be smaller in many cases.

In addition, the advantages are not limited to the case of parallel computing. As mentioned earlier (Remark 27), Kalikow decomposition could allow making linear the execution time with respect to the number of neurons, even in the case of sequential computation. This is why the first algorithm we present is sequential. We called it KOtS for Kalikow-Ogata's thinning Sequential. In a second sub-section, we present a partially parallelized version called KOtpP for Kalikow-Ogata's thinning partially parallelized.

These two algorithms are also based on the Ogata thinning method, which is a building of the 𝑠𝑝𝑖𝑘𝑒 trains in two steps. The first step produces the 𝑝𝑜𝑖𝑛𝑡𝑠 by simulation of Poisson processes (seeding step). The second step selects some of the 𝑝𝑜𝑖𝑛𝑡𝑠 to build the 𝑠𝑝𝑖𝑘𝑒 trains (thinning step).

�.�.� Sequential Kalikow-Ogata's thinning algorithm Kalikow decomposition

In this section, by adding Kalikow decomposition to the classical Ogata algorithm, we introduce two new algorithms in order to reduce the complexity of algorithms in previous sections. The Kalikow decomposition used here is unconditional see [18] and also Chapter 2 and Chapter 3.

Recall that V 𝑇 a family of the neighborhood is a set of all possible neighborhoods 𝑣 𝑇 (we refer to the formal definition of the neighborhood in Chapter 2 and Chapter 3). 𝑉 𝑇 is a random variable, taking values in V 𝑇 .

We also assume that ∅ belongs to V 𝑇 for all 𝑇.

On the neighborhood family V 𝑇 , we define a probability distribution 𝜆(.) such that P(𝑉 𝑇 = 𝑣 𝑇 ) = 𝜆(𝑣 𝑇 ).

For ). It depends at most on only one pre-synaptic 𝑛𝑒𝑢𝑟𝑜𝑛 and on a very small range of the past time. This is the most important motivation for using the Kalikow decomposition. It can be used to upper bound the number of past events that must be known to determine the selection or rejection of a 𝑝𝑜𝑖𝑛𝑡. If this upper bound remains constant whatever the number of neurons is, this could allow making linear the execution time with respect to the number of neurons. Moreover, as the Kalikow decomposition allows to simplify the dependencies, it is an efficient tool to implement a parallelized algorithm. Finally, by assume the second hypothesis (ii.), it ensures that the thinning procedure is applicable for every neighborhood 𝑣 𝑇 .

𝜙 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) = ⎧ ⎪ ⎨ ⎪ ⎩ 𝜇 𝑖 1 𝑇 -𝐿 𝑍 𝑖 (𝑇 )>𝛿 𝜆 𝑖 (𝑣 𝑇 ) if 𝑗 = ∅ 1 𝑇 -𝐿 𝑍 𝑖 (𝑇 )>𝛿 𝜆 𝑖 (𝑣 𝑇 ) ∫ 𝑇 -𝑛𝛿 𝑇 -(𝑛+1)𝛿 ℎ 𝑗𝑖 (𝑇 -𝑡)𝑑𝑍 𝑗 𝑡 if 𝑗 ≠ ∅. ( 4
Remark 30. In our case, a neighborhood is written as the Cartesian product of two subsets (4.3.13), one for space and the other for a time.

Thus, the probability 𝜆(𝑣 𝑇 ) can be written as the product of two probabilities, one in space, one in time: In addition, a neighborhood 𝑣 𝑇 is randomly drawn for each 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 . As notice in Remark 30, 𝑣 𝑇 is the Cartesian product of two independent subset of the past information. We decompose the neighborhood in a neighborhood in space {𝑗} and a neighborhood in time [𝑇 -(𝑛+1)𝛿, 𝑇 -𝑛𝛿] that are drawn independently, respectively with probability 𝜆 𝑠𝑝𝑎𝑐𝑒 ({𝑗}) and 𝜆 𝑡𝑖𝑚𝑒 ({𝑛}). Then, the knowedge of the neighborhood 𝑣 𝑇 allows to perform a first thinning, called a pre-thinning. Indeed, as mentioned in Remark 29, the intensity 𝜙 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) is upper bounded by a deterministic constant 𝑀 𝑣 𝑇 smaller than 𝑀 𝑖 . This allows to eventually reject the 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 by setting its thinning mark as follows.

𝜆(𝑣 𝑇 ) = 𝜆 𝑠𝑝𝑎𝑐𝑒 ({𝑗}) × 𝜆 𝑡𝑖𝑚𝑒 ({𝑛}) (4
𝑋 𝑇 = ⎧ ⎪ ⎨ ⎪ ⎩ 0 if 𝑈 𝑇 > 𝑀 𝑣 𝑇 𝑀 𝑖 𝑛.𝑎 otherwise. (4.3.16)
In the case of our implementation 𝑀 𝑣 𝑇 only depends on the neighborhood in space {𝑗}. For this reason, the algorithm checks first the 𝑝𝑜𝑖𝑛𝑡𝑠 that are rejected then picks the neighborhoods in time only for still unknown 𝑝𝑜𝑖𝑛𝑡𝑠.

Then the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 is built as the unions of all sets 𝒮 𝑖 . The set 𝒮 is sorted to finally contain all the 𝑝𝑜𝑖𝑛𝑡𝑠 in an increasing 𝑇-order.

During the second step, for each unknown 𝑝𝑜𝑖𝑛𝑡 𝑆 Then the thinning mark 𝑋 𝑇 can be defined as a realization of a random variable Bernoulli with parameter P(𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 ).

(𝑇 |ℱ 𝑇 ) = 𝜓 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑣 𝑇 ) = ⎧ ⎪ ⎨ ⎪ ⎩ 𝜇 𝑖 𝜆 𝑖 (∅) if 𝑗 = ∅ 1 𝜆 𝑖 (𝑣 𝑇 ) ∑ 𝑆 𝑗 𝑇 ′ ∈𝒮 𝑣 𝑇 ℎ 𝑗𝑖 (𝑇 -𝑇 ′ )1 𝑋 𝑇 ′=1 if 𝑗 ≠ ∅. ( 4
Once every 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 get its final thinning mark, the 𝑠𝑝𝑖𝑘𝑒𝑠 train is obtained by keeping only the selected 𝑝𝑜𝑖𝑛𝑡𝑠. Compute the number of 𝑝𝑜𝑖𝑛𝑡𝑠 : 𝑛 𝑖 follows a 𝑃 𝑜𝑖𝑠𝑠𝑜𝑛(𝑀 𝑖 𝑡 𝑚𝑎𝑥 ). ) by discussing cases on the refractory period as we did in previous section:

Pseudo-code

P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 ) ≤ P(𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 ) ≤ P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 ). (4.3.24)
If this range of probability [P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 ), P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 )] is small enough it allows to determine the thinning mark.

𝑋 𝑇 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 if 𝑈 𝑇 ≤ P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 ) 0 if 𝑈 𝑇 > P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 ) 𝑛.𝑎 otherwise. ( 4

�.�.� Applications

In this section, we explain how we concretely apply the algorithms to a multivariate linear Hawkes process with a hard refractory period 𝛿, where the intensity of index 𝑖 is recalled below. with 𝛽 𝑗𝑖 as follows

𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜓 𝑖 (𝑡|ℱ 𝑡 )1 𝑡-𝐿 𝑍 𝑖(𝑡)>𝛿 = ( 𝜇 𝑖 + ∑ 𝑗∈I ∫ 𝑡 0 ℎ 𝑗𝑖 (𝑡 -𝑠)𝑑𝑍 𝑗 𝑠 ) 1 𝑡-𝐿 𝑍 𝑖(𝑡)>𝛿 . ( 4
𝛽 𝑗𝑖 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 2|𝑗 -𝑖| 𝑝 if 𝑗 ≠ 𝑖 1 if 𝑗 = 𝑖. (4.3.28)

Upper bounds intensities for Ogata's thinning algorithm

The computation of the upper bound 𝑀 𝑖 is the same as pointed out in Lemma 1 of Section 4.6 of [18] or Chapter 3. We recall that 𝑀 𝑖 = 𝜇 𝑖 + ∑ 𝑗 𝛽 𝑗𝑖 1exp (-𝛼𝛿) the dominating intensity for index 𝑖.

Upper bounds intensities for Kalikow-Ogata's thinning algorithm

We prove the existence of 𝑀 𝑣 and 𝑀 𝑖 for a multivariate linear Hawkes process with the hard refractory period. The following calculations and proofs are mostly based on [18] or Chapter 3. Since it is sufficient to do the calculus at a time 𝑇 arbitrary, hence to simplify, we omit the subscript 𝑇 at the notations. We Following the method in [18], it admits a Kalikow decomposition at time 𝑇, for 𝑣

= {𝑗}×[𝑇 -(𝑛+1)𝛿, 𝑇 - 𝑛𝛿) ⎧ ⎪ ⎨ ⎪ ⎩ 𝜆 𝑖 (𝑣) = 𝜆 𝑗𝑖 > 0 𝜙 𝑣 𝑖 (𝑇 |ℱ 𝑇 ) = ∫ 𝑇 -𝑛𝛿 𝑇 -(𝑛+1)𝛿 ℎ 𝑗𝑖 (𝑇 -𝑠)𝑑𝑍 𝑗 𝑠 𝜆 𝑗𝑖 1 𝑇 -𝐿 𝑍 𝑖(𝑇 )>𝛿 . (4.3.29) Otherwise, if 𝑣 = ∅, we have ⎧ ⎪ ⎨ ⎪ ⎩ 𝜆 𝑖 (𝑣) = 𝜆 𝑖 (∅) > 0 𝜙 𝑣 𝑖 (𝑇 |ℱ 𝑇 ) = 𝜇 𝑖 𝜆 𝑖 (∅) 1 𝑇 -𝐿 𝑍 𝑖(𝑇 )>𝛿 (4.3.30)
with any sequence (𝜆 𝑗𝑖 ) 𝑗 such that ∑ 𝑗 𝜆 𝑗𝑖 = 1 -𝜆 𝑖 (∅).

Moreover, we also choose

𝜆 𝑗𝑖 = 𝐶 2|𝑗 -𝑖| 𝑝 = 𝐶𝛽 𝑗𝑖
with 𝐶 is a renormalization constant. Also from Lemma 1 in Section 4.6 of [18] we can give the dominating intensities 𝑀 𝑖 and 𝑀 𝑣 as follow:

for 𝑣 = {𝑗} × [𝑇 -(𝑛 + 1)𝛿, 𝑇 -𝑛𝛿), we have

𝑀 𝑣 = 𝑒 -𝑛𝛿 𝐶 .
Otherwise, if 𝑣 = ∅ we have 𝑀 ∅ = 𝜇 𝑖 𝜆 𝑖 (∅) .

Therefore, we can choose

𝑀 𝑖 = sup 𝑣 𝑀 𝑣 = max { 𝜇 𝑖 𝜆 𝑖 (∅)
, 𝑒 -𝛿 𝐶 } .

�.� Second part: Stochastically bounded intensities

Most of the time, finding a constant that dominates the intensity of a neuron for the whole time of simulation is impossible, for example, the Linear Hawkes process [10], etc. However, as stated earlier, it remains possible to construct a dominating intensity, which changes over time and dominates the intensity of every neuron as we go along, see Assumption 8. This idea is not new, it can be found in [16]. This dominating intensity is updated point after point, so we construct it gradually.

For the simulation purposes, recall that the point processes we consider have empty past, i.e no point before time 0. To simplify the writing, throughout this section, we set 𝑡𝑚𝑎𝑥 = ∞ and the goal is to simulate the first 𝑛 0 (deterministic) spikes. We focus on the two following versions of Linear Hawkes process. For the classical version of the Linear Hawkes process, we show later that in some cases, the dominating intensity Λ 𝑖 (𝑡) might be a non-decreasing simple predictable function [22]. In other words, the dominating intensity may explode in finite time. This explains why we propose to stop the algorithm after obtaining 𝑛 0 spikes.

For the Linear Hawkes process with bounded support, the dominating intensity Λ 𝑖 (𝑡) now can be decreased after a delay of time 𝐴.

With stochastically bounded intensity, the Ogata's algorithm [16] can be rewritten as follow i. Seeding step: We construct the dominating intensity Λ 𝑖 (𝑡) gradually. For each point 𝑆 𝑖 𝑇 , we update the dominating intensity to find a candidate point.

ii. Thinning step: It is selected or rejected via an exact computation of the intensity at the corresponding time.

�.�.� Sequential Ogata's thinning algorithm Overall description

In the following, recall that the dominating process Π 𝑖 correspond to index 𝑖 and has intensity Λ 𝑖 (𝑡).

Moreover, recall that Π = (Π 𝑖 ) 𝑖 and denote Π(1) for the list of accepted points after the thinning steps.

We consider the filtration F such that Λ 𝑖 (𝑡) is adapted to F𝑡 for every 𝑖: Intuitively, as soon as no point is produced, whatever the picked neighborhood 𝑣 is, the intensity function restricted to this neighborhood is bounded by a predictable function. That allows us, at the same time, to construct the dominating intensity and do a thinning step.

Remark 41. Assumption 10 is stronger than Assumption 8. Indeed, from Assumption 10, we conclude that 𝜙 𝑖 (𝑡|ℱ 𝑡 )1 𝑍(𝑇 ,𝑡)=0 ≤ 𝑀 𝑖 (𝑇 ) ∀𝑖 ∈ I, ∀𝑡 > 𝑇 , ∀𝑇 .

Analogous to Ogata's algorithm, we start with the seeding step, followed by a Kalikow step where we pick a random neighborhood. Finally, we do a thinning step to determine the marks. More precisely, 1. At the seeding step, we create a candidate to find the next point after time 𝑇. At this stage, we compute the dominating intensity 𝑀(𝑇 ) = ∑ 𝑁 𝑖=1 𝑀 𝑖 (𝑇 ) at time 𝑇, that creates a candidate 𝑇 ′ by adding an exponential of parameter 𝑀(𝑇 ). We attribute the neuron mark to this received candidate 𝑇 ′ , this candidate belongs to neuron 𝑖 with probability 𝑀 𝑖 (𝑇 )/𝑀(𝑇 ). This is exactly step 1 of the previous algorithm 2. At the Kalikow step, by using the Kalikow decomposition, we pick a random neighborhood of 𝑇 ′ . Without loss of generality, assume that 𝑉 𝑇 ′ = 𝑣 𝑇 ′ . Compute the dominating intensity of the candidate 𝑀(𝑇 ) = ∑ 𝑖 𝑀 𝑖 (𝑇 ), where 𝑀 𝑖 (𝑇 ) is defined in Assumption 10. Assume that, this candidate belongs to neuron 𝑖. Proof. The proof is similar to Proposition 2. Take 𝐶 𝑡 a positive, predictable function with respect to F𝑡 .

Let us denote 𝑇 with 𝑗 𝑇 = 𝑖 by 𝑇 𝑖 𝑘 for some 𝑘. We have This holds for any 𝐶 𝑡 that is F𝑡 predictable. Therefore, it implies that 𝜙 𝑖 (𝑡| F𝑡 ) is the F𝑡 intensity of Π 𝑖 (1). In other words, the obtaining points from Algorithm 9 admits 𝜙 𝑖 (𝑡| F𝑡 ) as stochastic intensity. This completes our proof.

E
Remark 42. To prove Proposition 10, it is then sufficient to consider the Kalikow decomposition in a particular case, where the neighborhood 𝑉 𝑇 takes the whole past before time 𝑇 with probability 1. Indeed, Ogata's algorithm is a particular case of the KalikowOgata algorithm by considering the neighborhood family V 𝑇 having only one element: the whole past before time 𝑇.

�.�.� Applications

To illustrate our simulation algorithms, in this section, we consider Linear Hawkes processes [10,11] and the Linear Hawkes processes with a bounded support [9]. Recall that the intensities are given as follow: where 𝐴 is a positive number.

Settings

We consider the parameters in very general settings, that is

• 𝜇 𝑖 a positive number, provides the intensity of the Hawkes process 𝑖 when it does not interact with any process in the network.

• Each interaction function (ℎ 𝑗𝑖 ) ∶ R + ↦ R + continuous and non increasing function, that represents the strength of the synaptic connection between the pre-synaptic process 𝑗 and the postsynaptic process 𝑖 [21].

Stochastic upper bounds for the intensities in Ogata's thinning algorithm In the following of this example (and only in this example), we will write ℎ 𝑗𝑖 as ℎ 𝑖 𝑗 to remove any unexpected confusion to the notations. In which, we recall that 𝑍((𝑎, 𝑏)) counts the number of spikes of the process 𝑍 in the interval (𝑎, 𝑏). Clearly, for each 𝑇 ′ rejected and until we meet a spike, the dominating intensity will decrease after delays of length 𝐴. We set 𝑀 𝑖 (𝑇 ) = sup ( sup 𝑘 ( 1 𝜆 𝑘𝑖 ( 𝑍 𝑘 ((0, 𝑇 )) + 1 ) ℎ 𝑘𝑖 (0) ) , 𝜇 𝑖 𝜆 0𝑖 ) .

For the Linear Hawkes process with bounded support, the result is very similar. We have Kalikow On Figure 4.4, we can observe that the vertical lines stay in the blue shaded region. This means that there is no significant correlation between the delays. simulating 400 000 𝑛𝑒𝑢𝑟𝑜𝑛𝑠. In general, in the huge network (more than 20 000), the KOtpP performs significantly better than OtS, and better than KOtS.

However, with the number of 𝑛𝑒𝑢𝑟𝑜𝑛𝑠 smaller than 20 000, the parallelization lost performances. Two reasons are noticed. The first one is the existence of a fixed cost in time with parallel computation. We assume that this is the synchronization time, it would probably be reduced with fewer alternations between sequential and parallel phases of computations. The other reason is the fact that for a small number of neurons it becomes more efficient to sort the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠. Indeed we noticed that the smaller

  The human brain has about 10 11 neurons that fire electrochemical signals to communicate with each other via specialized connections called synapses. These signals are transmitted from the cell body of a neuron through the axon (then synapse) to the next neuron (see Figure 1.1).

Figure 1

 1 Figure 1.1: Neuron Source: tutorix.com Figure 1.2: Action potential Source: teachmephysiology.com

𝜙

  𝑖 (𝑡|ℱ 𝑡 ) ≤ Λ 𝑡 and we set 𝜙 𝑁+1 (𝑡|ℱ 𝑡 ) = Λ 𝑡 -𝑁 ∑ 𝑖=1 𝜙 𝑖 (𝑡|ℱ 𝑡 ).

1 . 3 ,

 13 Figure 1.3: A neighborhood

Algorithm 3 1 : 6 :

 316 Backward Forward Algorithm ▷ With [𝑡 0 , 𝑡 1 ] the interval of simulation for neuron 𝑖 ∈ I Initialize the family V of non empty neighborhoods with {(𝑖, [𝑡 0 , 𝑡 1 ])} 2: Initialize the family of points P = ∅ ▷ Each point is a time 𝑇 with 3 marks: 𝑗 𝑇 is the neuron on which 𝑇 appears, 𝑉 𝑇 for the choice of neighborhood, 𝑋 𝑇 for the thinning step (accepted/rejected) 3: Draw 𝐸 an exponential variable with parameter 𝑀 4: Schedule 𝑇 𝑛𝑒𝑥𝑡 = 𝑡 0 + 𝐸 5: while 𝑇 𝑛𝑒𝑥𝑡 < 𝑡 1 do Append to P, the point 𝑇 𝑛𝑒𝑥𝑡 , with 3 marks: 𝑗 𝑇 𝑛𝑒𝑥𝑡 = 𝑖, 𝑉 𝑇 𝑛𝑒𝑥𝑡 = n.a. and 𝑋 𝑇 𝑛𝑒𝑥𝑡 = n.a. (n.a. stand for not assigned yet) 7: while There are points 𝑇 in P with 𝑉 𝑇 = n.a. do 8: for each point 𝑇 in P with 𝑉 𝑇 = n.a. do 9:

10 : 12 :

 1012 if 𝑉 𝑇 ≠ ∅ then11: Find the portion of time/neurons in 𝑉 𝑇 which does not intersect the existing non empty neighborhoods in V Simulate on it a Poisson process with rate 𝑀 13: Append the simulated points, 𝑇 ′ , if any, to P with their neuron 𝑗 𝑇 ′ and with 𝑉 𝑇 ′ = 𝑋 𝑇 ′ =n.a the 𝑇's in P with 𝑋 𝑇 =n.a. in increasing order 19: for each of them starting with the most backward do 20: Draw 𝑋 𝑇 as a Bernoulli variable with parameter 𝜙 𝐸 ′ another exponential variable with parameter 𝑀 23: 𝑇 𝑛𝑒𝑥𝑡 ← 𝑇 𝑛𝑒𝑥𝑡 + 𝐸 ′ 24: end while 25: The desired points are the points in P with marks 𝑗 𝑇 = 𝑖, 𝑋 𝑇 = 1 and that appear in [𝑡 0 , 𝑡 1 ]

Figure 4

 4 Figure 4 describes an example to go step by step through Algorithm 3. The backward steps determine all

2Remark 4 .

 4 Event-Scheduling Algorithms with Kalikow Decomposition for Simulating Potentially Infinite Neuronal Networks Main flow example for Algorithm 3, with backward steps (cf. Algorithm 3, Steps 7-17) and forward steps (cf. Algorithm 3, Steps 18-25). Following circled numbers: (1) The next point 𝑇 𝑛𝑒𝑥𝑡 = 𝑡 0 +𝐸 (cf. Algorithm 3, Step 4) is scheduled, (2) The neighborhood 𝑉 𝑇 𝑛𝑒𝑥𝑡 is selected in the first backward step, a first generation of three points (𝑎, 𝑏 on neuron 𝑘 and 𝑐 on neuron ℓ) is drawn (cf. Algorithm 3, Step 9), thanks to a Poisson process, (cf. Algorithm 3, Steps 11-12) and appended to P (cf. Algorithm 3, Step 13), (3) at the second generation, a non empty neighborhood is found, i.e. 𝑉 𝑏 ≠ ∅ (cf. Algorithm 3, Steps 9-1), but the Poisson process simulation does not give any point in it (cf. Algorithm 3, Step 12), (4) at the second generation, the neighborhood 𝑉 𝑎 is picked, it is not empty and overlap the neighborhood of the first generation (cf. Algorithm 3, Steps 9-11): therefore there is no new simulation in the overlap (𝑐 is kept and belongs to 𝑉 𝑏 as well as 𝑉 𝑎 ) but there is a new simulation thanks to a Poisson process outside of the overlap leading to a new point 𝑑 (cf. Algorithm 3, Step 12)(5) at the second generation, for point 𝑐, one pick the empty neighborhood, i.e. 𝑉 𝑐 = ∅ (cf. Algorithm 3, Step 9) and therefore we do not simulate any Poisson process, (6) at third generation, similarly no point and no interval are generated, i.e. 𝑉 𝑑 = ∅ (cf. Algorithm 3, Step 9). This is the end of the backward steps and the beginning of the forward ones, (7) the point 𝑑 is not selected, acceptation/selection taking place with probability 𝜙 ∅ ℓ 𝑀 (cf. Algorithm 3, Step 20), (8) the point 𝑐 is accepted, here again with probability 𝜙 ∅ ℓ 𝑀 (cf. Algorithm 3, Step 20), (9) the point 𝑏 is not selected, acceptation taking place, here, with probability 𝜙 𝑉 𝑏 𝑘 (𝑏) 𝑀 (cf. Algorithm 3, Step 20), (10) the point 𝑎 is selected, acceptation taking place, here, with probability 𝜙 𝑉 𝑎 𝑘 (𝑎) 𝑀 (cf. Algorithm 3, Step 20), (11) The neighborhood of neuron 𝑖 contains two points, one on neuron 𝑘 and one on neuron ℓ and one selects 𝑇 𝑛𝑒𝑥𝑡 with probability 𝜙 𝑉 𝑇 𝑛𝑒𝑥𝑡 𝑖 (𝑇 𝑛𝑒𝑥𝑡 ) 𝑀 .

of the time simulation for neurons 0 to 9 Figure 2 . 4 :

 924 Figure 2.4: Simulation for 𝑀 = 2, 𝜎 = 1, 𝜆 ∅ = 0.25. For each neuron in Z 2 , that have been requested in Steps 9:11,except the neuron of interest (0, 0), have been counted the total number of requests, that is the number of time a 𝑉 𝑇 pointed towards this neuron (Steps 9 and 11) and the total time spent at this precise neuron simulating a homogeneous Poisson process (Step 12). Note that since the simulation is on [0, 100] the time spent at position (0, 0) is at least 100. On (a), the summary for one simulation with below the plot, number of points accepted at neuron (0, 0) and total number of points that have been simulated. Also annotated on (a), with labels between 0 and 8, the 9 neurons for which the same simulation in[20,40] is represented in more details on (b). More precisely on (b), in abscissa is time and the neuron labels in ordinate. A plain dot represents an accepted point, by the thinning step (Step 20 of Algorithm 3), and an empty round, a rejected point. The blue pieces of line represent the non empty neighborhoods that are in V.

Figure 2 . 5 :

 25 Figure 2.5: Simulation for 4 other sets of parameters, all of them with 𝜎 = 3. Summaries as explained in Figure 2.4. On top, 𝑀 = 2; on bottom, 𝑀 = 20. On the left part, 𝜆 ∅ = 0.25, on the right part, 𝜆 ∅ = 0.5.

ℎ

  0) 𝜆 𝑖 ( 𝑣 𝑖 (𝑙,𝑛)) = 𝜂 (𝑙,𝑛) 𝑗 𝑙 𝑖 (-𝑠)𝑑𝑥 𝑗 𝑙 𝑠 𝜂 (𝑙,𝑛) , for any choice of non negative weights ( 𝜂 𝑖 (𝑙,𝑛)) (𝑙,𝑛)∈N * ×N *

Assumption 3 .

 3 (i) For all 𝑖, 𝑗 ∈ I, ℎ 𝑗𝑖 (.) is a non negative, non increasing 𝐿 1 function. Moreover, for every 𝑖∑ 𝑗∈I ‖ℎ 𝑗𝑖 ‖ 𝐿 1 < ∞.

Proof.

  For any 𝑥 ∈ 𝒳 >𝛿 , 𝑖 ∈ I and 𝑘 ≥ 2, consider Δ 𝑖 𝑘 (𝑥) = 𝜓 𝑖 ( ∑ 𝑗∈𝑉 𝑖 (𝑘)

  [𝑛] 𝑖 (𝑥) → 𝜙 𝑖 (𝑥) as 𝑛 → ∞. As a consequence, Assumption 1 is fulfilled, and by Proposition 3, the conclusion follows.

Remark 9 .

 9 Denote 𝐷 𝑘 𝑖 (𝑥) = {𝑧 ∈ 𝒳 >𝛿 ∶ 𝑧 𝑣 𝑖 𝑘 = 𝑥}. Clearly,

  1); 𝑠 < 𝑡) and ℱ

Proposition 7 .Remark 15 .

 715 If sup 𝑖∈I ∑ 𝑘≥1 𝑃 (𝑣 𝑖 𝑘 )𝜆 𝑖 (𝑣 𝑖 𝑘 ) < 1 (3.5.1) then the Backward steps in the Perfect Simulation algorithm terminate almost surely in finite time. For any neighborhood 𝑣, we have ∑ 𝑗∈𝑣 E ( Π 𝑗( 𝜋 𝑗 (𝑣) )) = 𝑃 (𝑣).

  ancestors in the 𝑛 𝑡ℎ set of ancestors of type 1 number of ancestors in the 𝑛 𝑡ℎ set of ancestors of type 2 … number of ancestors in the 𝑛 𝑡ℎ set of ancestors of type 𝑁 that 𝐾 𝑖 (0) = 𝑒 𝑖 for every 𝑖. For every 𝑗 ∈ I, denote 𝑋 𝑗 𝑖 ∶= 𝐾 𝑗 𝑖 (1) the number of ancestors of type 𝑗 in the first set of ancestors with initial point of type 𝑖. Moreover, as we discussed earlier, 𝑋 𝑗 𝑖 is the cardinal of the points that a Poisson process of intensity Γ 𝑗 puts on 𝜋 𝑗 (𝑉 𝑖 ), where 𝑉 𝑖 is the random neighborhood from a point of type 𝑖. In other words, if we denote 𝒫 the Poisson distribution, conditioning on 𝑉 𝑖 = 𝑣, we have 𝑋 𝑗 𝑖 ∼ 𝒫 ( Γ 𝑗 𝜇(𝜋 𝑗 (𝑣)) ) .

  Now we are looking for values of 𝑝 such that: (a) (𝜂 𝑘 ) 𝑘 defines a probability, (b) the algorithm stays in 𝒴, (c) the Kalikow decomposition exists (Proposition 4), (d) the branching process goes extinct in finite time almost surely.

time 0 .

 0 Let I be the set of indices of these processes, I = {1, 2, … , 𝑁}. Denote 𝑍 𝑖 the counting process associated with index 𝑖. For any 0 ≤ 𝑠 < 𝑡 < 𝑡 𝑚𝑎𝑥 , 𝑍 𝑖 ((𝑠, 𝑡]) records the number of points produced by the process 𝑍 𝑖 during the interval (𝑠, 𝑡]. Denote ℱ a history of the process 𝑍 = (𝑍 𝑖 ) 𝑖∈𝐼 (see Introduction of the thesis, Chapter 1). Let us recall that (𝜙 𝑖 (𝑡|ℱ 𝑡 )) 𝑡 is the intensity function of process 𝑍 𝑖 .

3 : 4 : 5 : 6 :

 3456 for each 𝑝𝑜𝑖𝑛𝑡 indexed by 𝑘 ≤ 𝑛 𝑖 do Simulate the time : 𝑇 𝑖 𝑘 follows a 𝑈 ([0, 𝑡 𝑚𝑎𝑥 ]) distribution Define the uniform : 𝑈 𝑇 𝑖 𝑘 follows a 𝑈 ([0, 1]) distribution Set the 𝑝𝑜𝑖𝑛𝑡 as unknown : the thinning mark 𝑋 𝑇 𝑖 𝑘 = 𝑛.𝑎 7:

Algorithm 5 2 : 3 :

 523 Algorithm OtpP (for Ogata's thinning partially parallelized) ▷ Step 1: Seeding step 1: for each 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑘 in parallel do for each 𝑛𝑒𝑢𝑟𝑜𝑛 𝑖 ∈ population do Compute the number of 𝑝𝑜𝑖𝑛𝑡𝑠 : 𝑛 𝑖 follows a 𝑃 𝑜𝑖𝑠𝑠𝑜𝑛(𝑀 𝑖 𝑡 𝑚𝑎𝑥 ).

4 : 5 : 6 : 7 :

 4567 for each 𝑝𝑜𝑖𝑛𝑡 indexed by 𝑙 ≤ 𝑛 𝑖 do Simulate the time : 𝑇 𝑖 𝑙 follows a 𝑈 ([0, 𝑡 𝑚𝑎𝑥 ]) distribution Define the uniform : 𝑈 𝑇 𝑖 𝑙 follows a 𝑈 ([0, 1]) distribution Set the 𝑝𝑜𝑖𝑛𝑡 as unknown : the thinning mark 𝑋 𝑇 𝑖 𝑙 𝒮 𝑘 𝑝𝑜𝑝 = ⋃ 𝒮 𝑖 the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠 belonging to the 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 11: Sort this set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 𝑘 𝑝𝑜𝑝 = ⋃ 𝒮 𝑖 in an increasing 𝑇-order 12: end for 13: Define the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠 : 𝒮 = ⋃ 𝒮 𝑘 𝑝𝑜𝑝 ▷ Step 2: Thinning 14: while At least one 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 ∈ 𝒮 is unknown do 15:

Remark 22 .Remark 23 .

 2223 end while ▷ Final step : format the spikes trains26: Delete the points that have a thinning mark 𝑋 𝑇 = 0 The algorithm OtpP stops after a finite number of steps almost surely. Indeed only a finite number of 𝑝𝑜𝑖𝑛𝑡𝑠 are produced after the seeding step. Then at each WHILE loop, at least one 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 accepts or rejects a 𝑝𝑜𝑖𝑛𝑡. Indeed the 𝑝𝑜𝑖𝑛𝑡 with the smallest time can always be accepted or rejected. The algorithm OtpP is more general than the algorithm OtS (the sequential version). Indeed the algorithm OtpP with a single 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 is completely equivalent to the Algorithm OtS.

Algorithm 6

 6 Algorithm KOtS (for Kalikow-Ogata's thinning Sequential) ▷ Step 1: Seeding step 1: for each 𝑛𝑒𝑢𝑟𝑜𝑛 𝑖 do 2:

  .3.26) Settings An interaction function ℎ 𝑗𝑖 is an exponentially decreasing function multiplied by a constant 𝛽 𝑗𝑖 which depends on the distance |𝑖 -𝑗| between the two 𝑛𝑒𝑢𝑟𝑜𝑛𝑠. ℎ 𝑗𝑖 (𝑡) = 𝛽 𝑗𝑖 exp(-𝛼𝑡) (4.3.27)

4. 3

 3 First part: Deterministically bounded intensities and Parallelization consider the family of neighborhood V with 𝑣 = {𝑗} × [𝑇 -(𝑛 + 1)𝛿, 𝑇 -𝑛𝛿) for 𝑗 ∈ I and 𝑣 ∅ = ∅.

[

  Linear Hawkes process] 𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜇 𝑖 + ∑ process with bounded support 𝐴 > 0] 𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜇 𝑖 + ∑

3 ) 10 ( 6 ) 2 .

 31062 F𝑡 = ⋁ 𝑗∈I {(𝑇 , 𝑋 𝑇 ) ∶ 𝑇 ∈ Π 𝑗 ∩ [0, 𝑡)} (4.4.Assumption Stronger version of the stochastically bounded intensity). There exists a dominating point process Π = (Π 𝑖 ) 𝑖 such that 1. At any 𝑇 of the dominating process Π, there exists 𝑀 𝑖 (𝑇 ) a F𝑇 measurable random variable such that sup 𝑡>𝑇 sup 𝑣 𝑡 ∈V 𝑡 𝜙 𝑣 𝑡 𝑖 (𝑡|ℱ 𝑡 )1 𝑍(𝑇 ,𝑡)=0 ≤ 𝑀 𝑖 (𝑇 ) < ∞. (4.4.For each index 𝑖, Π 𝑖 has F intensity Λ 𝑖 (𝑡) = ∑ 𝑇 𝑀 𝑖 (𝑇 )1 𝑇 <𝑡≤𝑛 Π (𝑇 ) (4.4.7) Note that, Λ 𝑖 (𝑡) ≥ 𝜙 𝑖 (𝑡|ℱ 𝑡 ) almost surely and that therefore one might deduce points of 𝑍 by thinning points of Π and recall that F is the extended filtration of ℱ (4.4.3).

3 .Algorithm 9

 39 At the thinning step, we then assign to it a thinning mark with Bernoulli distribution of parameterP(𝑋 𝑇 ′ = 1|𝑉 𝑇 ′ = 𝑣 𝑇 ′ , F𝑇 ′ ) = 𝜙 𝑣 𝑇 ′ 𝑖 (𝑇 ′ | F𝑇 ′ )/𝑀 𝑖 (𝑇 ) .The thinning mark at the candidate 𝑇 ′ is determined by𝑋 𝑇 ′ = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if 𝑈 𝑇 ′ ≤ P(𝑋 𝑇 ′ = 1|𝑉 𝑇 ′ = 𝑣 𝑇 ′ , F𝑇 ′ ) AlgorithmKOtSS for KalikowOgata thinning sequential with stochastic bounded intensity 1: Initialize 𝑇 = 0. 2: Initialize number of spikes 𝑛 = 0. 3: while 𝑛 ≤ 𝑛 0 do ▷ Step 1: Seeding 4:

5 : 6 :

 56 Draw a candidate 𝑇 ′ ← 𝑇 + 𝐸𝑥𝑝(𝑀(𝑇 ))Distribute to 𝑇 ′ a mark 𝑖 with probability 𝑀 𝑖 (𝑇 )/𝑀(𝑇 ) ▷ Step 2: Kalikow 7:

8 :Proposition 11 .

 811 Pick a neighborhood 𝑉 𝑇 ′ = 𝑣 𝑘 with probability 𝜆 𝑖,𝑇 ′ (𝑣 𝑘 ) in Kalikow decomposition. To simplify set 𝑣 𝑇 ′ = 𝑣 𝑘 . ▷ Step 3: Thinning 9: Compute the thinning rate P(𝑋 𝑇 ′ = 1|𝑉 𝑇 ′ = 𝑣 𝑇 ′ , F𝑇 ′ ) 10: Determine the thinning mark 𝑋 𝑇 ′ by Equation (4.4.8). If accepted, increase 𝑛 by 1. 11: Update 𝑇 ← 𝑇 ′ 12: end while ▷ Final step : format the spikes trains 13: Delete the points that have a thinning mark 𝑋 𝑇 = 0By using Proposition 1 in[16], we can prove as before that the simulated points after Step 1 is a multivariate point process Π = (Π 𝑖 ) 𝑖 having conditional intensity (Λ 𝑖 (𝑡)) 𝑖 respectively up to time 𝜏. We can prove that Conditioning on the event {𝜏 < ∞}, the point process Π 𝑖 (1) obtained by Algorithm 9 admits 𝜙 𝑖 (𝑡| F𝑡 ) as F𝑡 -predictable intensity up to the stopping time 𝜏.

[

  Linear Hawkes process] 𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜇 𝑖 + ∑ process with a bounded support] 𝜙 𝑖 (𝑡|ℱ 𝑡 ) = 𝜇 𝑖 + ∑

For

  any 𝑇 and 𝑖, we recall the definition of upper stochastical bound 𝑀 𝑖 (𝑇 ) in Assumption 8, sup 𝑡>𝑇 𝜙 𝑖 (𝑡|ℱ 𝑡 )1 𝑍(𝑇 ,𝑡)=0 ≤ 𝑀 𝑖 (𝑇 ).

  With the Linear Hawkes process (4.4.1), we split the integral into two partssup 𝑡>𝑇 𝜙 𝑖 (𝑡|ℱ 𝑡 )1 𝑍(𝑇 ,𝑡)=0 = 𝜇 𝑖 + sup 𝑡>𝑇 ( ∑ 𝑗 ∫ 𝑇 0 ℎ 𝑖 𝑗 (𝑡 -𝑠)𝑑𝑍 𝑗 𝑠 + ℎ 𝑖 𝑗 𝑇 (𝑡 -𝑇 )1 𝑋 𝑇 =1 ) .Finally, since ℎ 𝑖 𝑗 is non increasing function, we set𝑀 𝑖 (𝑇 ) = 𝜇 𝑖 + 𝑍((0, 𝑇 )) max 𝑗 ℎ 𝑖 𝑗 (0) + ℎ 𝑖 𝑗 𝑇 (0)1 𝑋 𝑇 =1 .

Moreover, Λ 1

 Λ1 𝑖 (𝑡) = ∑ 𝑇 𝑀 𝑖 (𝑇 )1 𝑇 <𝑡≤𝑛 Π (𝑇 ) is ℱ Π 𝑡 predictable, non decreasing function where recall that 𝑛 Π (𝑇 ) is the next point after time 𝑇 of the process Π. Indeed, denote 𝑇 ′ = 𝑛 Π (𝑇 ), we have𝑀 𝑖 (𝑇 ) = 𝜇 𝑖 +𝑍((0, 𝑇 )) max 𝑗 ℎ 𝑖 𝑗 (0)+ℎ 𝑖 𝑗 𝑇 (0)1 𝑋 𝑇 =1 ≤ 𝜇 𝑖 +𝑍((0, 𝑇 ′ )) max 𝑗 ℎ 𝑖 𝑗 (0)+ℎ 𝑖 𝑗 𝑇 ′ (0)1 𝑋 𝑇 ′=1 = 𝑀 𝑖 (𝑇 ′ ).We also can construct a more precise dominating intensity that decreases after each rejected point by considering𝑀 𝑖 (𝑇 ) = 𝜇 𝑖 + ∑ 𝑗 ∫ 1 𝑋 𝑇 =1 .That is also ℱ Π 𝑇 measurable. Moreover, if 𝑇 ′ = 𝑛 Π (𝑇 ) is the candidate point after 𝑇 and 𝑋 𝑇 ′ = 0, we have𝑀 𝑖 (𝑇 ′ ) = 𝜇 𝑖 + ∑ 𝑗 ∫ 𝑋 𝑇 =1 < 𝑀 𝑖 (𝑇 )since ℎ 𝑖 𝑗 is a non-increasing function. We obtain a similar result with the Linear Hawkes process with bounded support (4.4.10), 𝑀 𝑖 (𝑇 ) = 𝜇 𝑖 + 𝑍([𝑇 -𝐴, 𝑇 )) max 𝑗 ℎ 𝑖 𝑗 (0) + ℎ 𝑖 𝑗 𝑇 (0)1 𝑋 𝑇 =1 .

4. 4

 4 Second part: Stochastically bounded intensities Stochastic upper bounds for the intensities in KalikowOgata's thinning algorithm Firstly, we begin to write Kalikow decomposition for the intensity function 𝜙 𝑖 (𝑡|ℱ 𝑡 ). Let us consider the neighborhood family V 𝑡 with 𝑣 𝑘 = {𝑘}×[0, 𝑡) for 𝑘 ∈ I. Follow the method in [18], the intensity of Linear Hawkes process admits a Kalikow decomposition at time 𝑡 for 𝑘 > 0: each 𝑇, by Assumption 10, sup 𝑡>𝑇 sup 𝑘 𝜙 𝑣 𝑘 𝑖 (𝑡|ℱ 𝑡 )1 𝑍((𝑇 ,𝑡))=0 ≤ 𝑀 𝑖 (𝑇 ).

  decomposition at time 𝑡 for 𝑘 > 0: 𝑖( 𝑣 𝑘) = 𝜆 𝑘𝑖 > 0 𝜙 𝑣 𝑘 𝑖 (𝑡|ℱ 𝑡 ) = ∫ 𝑡 𝑡-𝐴 ℎ 𝑘𝑖 (𝑡 -𝑠)𝑑𝑍 𝑘 𝑠 𝑘 ((𝑇 -𝐴, 𝑇 )) + 1 ) ℎ 𝑘𝑖 (0) ) , 𝜇 𝑖 𝜆 0𝑖 ) .

Figure 4 . 2 :

 42 Figure 4.2: On each subfigure, the x-axis displays the number of lags and the y-axis displays the auto-correlation at that number of lags. By default, the plot starts at lag = 0 and the auto-correlation will always be 1 at lag = 0. The blue shaded region is the confidence interval with a default value of 𝛼 = 0.05.

Figure 4 . 4 :

 44 Figure 4.4: On each subfigure, the x-axis displays the number of lags and the y-axis displays the auto-correlation at that number of lags. By default, the plot starts at lag = 0 and the auto-correlation will always be 1 at lag = 0. The blue shaded region is the confidence interval with a default value of 𝛼 = 0.05.

Figure 4 . 6 :

 46 Figure 4.6: On each subfigure, the x-axis displays the number of lags and the y-axis displays the auto-correlation at that number of lags. By default, the plot starts at lag = 0 and the auto-correlation will always be 1 at lag = 0. The blue shaded region is the confidence interval with a default value of 𝛼 = 0.05.

Figure 4 . 7 : 2 .

 472 Figure 4.7: Execution time of the algorithms sequential Ogata's thinning, sequential Kalikow-Ogata's thinning, and partially parallel Kalikow-Ogata's thinning with 2 and 4 𝑝𝑜𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠: For each number of neurons, the algorithms are run 3 times. They simulate the process during 5s. The median values are used to perform a polynomial regression of degree 2 .

1

  Introduction value only depends on the points appearing in the neighborhood 𝑣, and a probability density function 𝜆 𝑖 (.) such that 𝜙 𝑖 (𝑡 | ℱ 𝑡 ) = 𝜆 𝑖 (∅)𝜙 ∅ 𝑖 (𝑡) + ∑ Note that, at this stage, we only define the Kalikow decomposition for a bounded intensity and also we need to restrict our self to bounded 𝜙 𝑣 . Within this new type of Kalikow decomposition, we propose a theoretical algorithm that combines the idea of Ogata algorithm and Kalikow decomposition, see Algorithm 2 of Chapter 2. But this algorithm is purely theoretical since the computation of 𝜙

		𝜆 𝑖 (𝑣) × 𝜙 𝑣 𝑖 (𝑡|ℱ 𝑡 )	(1.5.1)
		𝑣∈V,𝑣≠∅
	with 𝜆 𝑖 (∅) + ∑ 𝑣∈V,𝑣≠∅	𝜆 𝑖 (𝑣) = 1.
	Remark 2. 𝑉 𝑇

  Compute the time 𝑇 of the next spike of the system after 𝑡, and the neuron where the spike occurs by 𝑇 ← min 𝑖∈ I 𝑇 𝑖 , with 𝑗 𝑇 ← arg min 𝑖∈ I 𝑇 𝑖

	2 Event-Scheduling Algorithms with Kalikow Decomposition for Simulating Potentially Infinite Neuronal
	Networks
	Algorithm 1 Classical point process simulation algorithm Note that the quantity 𝜙 𝑎𝑏𝑠 𝑖 (𝑠, 𝑡) can be also seen as the hazard rate of the next potential point 𝑇
		▷ With [𝑡 0 , 𝑡 1 ] the interval of simulation
	1: Initialize the family of points P = ∅
		▷ Each point is a time 𝑇 with a mark, 𝑗 𝑇 , which is the neuron on which 𝑇 appears
	2: Initialize 𝑡 ← 𝑡 0
	3: repeat
	4:	for each neuron 𝑖 ∈ I do
	5:	Draw independently an exponential variable 𝐸 𝑖 with parameter 1
	6:	Apply the inverse transformation, that is, find 𝑇 𝑖 such that
		∫ 𝑡	𝑇 𝑖	𝜙 𝑎𝑏𝑠 𝑖 (𝑠, 𝑡)𝑑𝑠 = 𝐸 𝑖 .
	7:	end for
	8:
	9:	if 𝑇 ≤ 𝑡 1 then
	10:	append 𝑇 with mark 𝑗 𝑇 to P
	11:	end if
	12:	𝑡 ← 𝑇
	13: until 𝑡 > 𝑡 1

6: Draw 𝑋 𝑇 a Bernoulli variable with parameter 𝜙 𝑉 𝑇 𝑖 (𝑇 ) 𝑀 7: if 𝑋 𝑇 = 1 and 𝑇 ≤ 𝑡 1 then

	8:	append 𝑇 to P
	9:	end if
	10:	𝑡 ← 𝑇
	11: until 𝑡 > 𝑡 1
	This Algorithm is close to Algorithm 1 but adds a neighborhood choice (Step 5) with a thinning step
	(Steps 6-9).
	In Appendix 2.8, we prove that this algorithm indeed provides a point process with an intensity given
	by (2.4.1) shifted at time 𝑡.
	The previous algorithm cannot be put into practice because the computation of 𝜙 𝑉 𝑇

  𝑡∈R is the counting process associated to the point processs P simulated by Algorithm 2. Let us denote by 𝜑 𝑖 (𝑡), the formula given by (2.4.1) and shifted at time 𝑡. Note that since the 𝜙 𝑣 Proof. Following the technique in Chapter 2 of[2], let us take 𝐶 𝑡 a non negative predictable function with respect to (w.r.t) ℱ

		𝑖 𝑠 (1); 𝑠 ≤ 𝑡) and ℱ 𝑡 𝑁(1)	𝑖∈I = ⋁	𝜎(𝑁 𝑖 𝑠 (1); 𝑠 ≤ 𝑡).
	Moreover note that (𝑁 𝑖 𝑡 (1)) 𝑁(1)
		𝑡 -
		𝑁(1) 𝑡	-predictable intensity.
	𝑁 𝑖 (1) 𝑡	𝑁(1) 𝑡-thatis ℱ

𝑖 's are ℱ 𝑖𝑛𝑡 0-= ℱ 𝑁(1) 0-, 𝜑 𝑖 (𝑡) is ℱ measurable. We also denote 𝜑 𝑣 𝑖 (𝑡) the formula of 𝜙 𝑣 𝑖 shifted at time 𝑡.

With this notation, we can prove the following.

Proposition 2. The process (𝑁 𝑖

𝑡 (1)) 𝑡∈R admits 𝜑 𝑖 (𝑡) as ℱ measurable and therefore ℱ 𝑁 𝑡-measurable . We have, for any 𝑖 ∈ I,

Definition 3.3.1. For

  then for any 𝑖 ∈ I, we denote the point measure associated to index 𝑖 by Let 𝑥 𝑡 = ( {𝑡 𝑖 𝑛 } 𝑛∈Z) 𝑖∈I ∈ 𝒳 𝑡 and denote 𝑥 ←𝑡 𝑡 ∶= ( {𝑡 𝑖 𝑛 -𝑡} 𝑛) 𝑖 ∈ 𝒳 the shifted configuration. In this paper, if we do not mention otherwise, we always consider time homogeneous point processes, which in our setting can be defined as follows. a given 𝑖 ∈ I, a counting process 𝑍 𝑖 with stochastic intensity (𝜙 𝑖,𝑡 (𝑥 𝑡 )) 𝑡∈R is said to be time homogeneous if 𝜙 𝑖,𝑡 (𝑥 𝑡 ) = 𝜙 𝑖 (𝑥 ←𝑡 𝑡 ) for all 𝑡 ∈ R and 𝑥 𝑡 ∈ 𝒳 𝑡 . Before giving the definition of a Kalikow decomposition, we introduce the definition of a cylindrical function as follows. For any neighborhood 𝑣 𝑡 ∈ V 𝑡 and 𝑥 𝑡 , 𝑦 𝑡 ∈ 𝒳 𝑡 , we say 𝑥 𝑡 𝑣 𝑡 = 𝑦 𝑡 whenever 𝑥 𝑡 = 𝑦 𝑡 in 𝑣 𝑡 . This means that, for all 𝑖 ∈ I, 𝑛 ∈ Z, such that 𝑡 𝑖 𝑛 ∈ 𝑥 𝑡 and (𝑖, 𝑡 𝑖 𝑛 ) ∈ 𝑣 𝑡 , we have 𝑡 𝑖 𝑛 ∈ 𝑦 𝑡 and vice-versa. A real valued function 𝑓 is called cylindrical in 𝑣 𝑡 if 𝑓 (𝑥 𝑡 ) = 𝑓 (𝑦 𝑡 ) for any 𝑥 𝑡 𝑣 𝑡 = 𝑦 𝑡 , and we usually stress the dependence in 𝑣 𝑡 by denoting 𝑓 𝑣 𝑡 (𝑥 𝑡 ). what follows, we give the definition of the Kalikow decomposition for a counting process 𝑍 𝑖 with stochastic intensity 𝜙 𝑖,𝑡 (𝑥 𝑡 ). We say a time homogeneous process 𝑍 𝑖 for some 𝑖 ∈ I admits the Kalikow decomposition with respect to (w.r.t) a neighborhood family (V 𝑡 ) 𝑡∈R and a sequence of subspaces (𝒴 𝑡 ) 𝑡∈R of 𝒳 ∞ , if for all 𝑡, the intensity 𝜙 𝑖,𝑡 (𝑥 𝑡 ) admits a convex decomposition for any past configuration 𝑥 𝑡 ∈ 𝒳 𝑡 ∩ 𝒴 𝑡 , that is, for any 𝑣 𝑡 ∈ V 𝑡 there exists a cylindrical function 𝜙 𝑣 𝑡 𝑖,𝑡 (.) on 𝑣 𝑡 taking values in R + and a probability density function 𝜆 𝑖,𝑡 (.) such that ∀𝑥 𝑡 ∈ 𝒳 𝑡 ∩ 𝒴 𝑡 , 𝜙 𝑖,𝑡 (𝑥 𝑡 ) = 𝜆 𝑖,𝑡 (∅)𝜙 ∅ 𝑖,𝑡 + ∑ 𝑣 𝑡 ∈V 𝑡 ,𝑣 𝑡 ≠∅ 𝜆 𝑖,𝑡 (𝑣 𝑡 )𝜙 𝑣 𝑡 ∈V 𝑡 ,𝑣 𝑡 ≠∅ 𝜆 𝑖,𝑡 (𝑣 𝑡 ) = 1.

	𝑑𝑥 𝑖 𝑠 = ∑ 𝑚∈Z Definition 3.3.2. In Definition 3.3.3. 𝑣 𝑡 𝑖,𝑡 (𝑥 𝑡 ) with 𝜆 𝑖,𝑡 (∅) + ∑	(3.3.1)

- 𝛿 𝑡 𝑖 𝑚 (𝑑𝑠) in which 𝛿 𝑡 (.) is a Dirac measure at 𝑡. Throughout this article, without further mentioning, the integral ∫ 𝑏 𝑎 stands for ∫ [𝑎,𝑏) with 𝑎, 𝑏 ∈ R.

  Since 𝑍 𝑖 is a time homogeneous process, by Definition 3.3.1, we have that 𝜙 𝑖,𝑡 (𝑥 𝑡 ) = 𝜆 𝑖,𝑡 (∅)𝜑 ∅ 𝑖,𝑡 + ∑ 𝑣 →𝑡 ∈V 𝑡 ,𝑣 →𝑡 ≠∅ 𝜆 𝑖,𝑡 (𝑣 →𝑡 )𝜑 𝑣 →𝑡 𝑖,𝑡 (𝑥 𝑡 )

			(3.3.3)
	with
	𝜆 𝑖,𝑡 (∅) +	∑
			.3.2)
	Since Equation (3.3.1) is satisfied at time 0 for any 𝑥 ∈ 𝒳 ∩ 𝒴 and since 𝑥 ←𝑡 𝑡 ∈ 𝒳 ∩ 𝒴, this implies that
	𝜙 𝑖 (𝑥 ←𝑡 𝑡 ) = 𝜆 𝑖,0 (∅)𝜙 ∅ 𝑖,0 + ∑ 𝑣∈V 0 ,𝑣≠∅	𝜆 𝑖,0 (𝑣)𝜙 𝑣 𝑖,0 (𝑥 ←𝑡 𝑡 ).

Define the cylindrical function 𝜑

𝑣 →𝑡 𝑖,𝑡 , that is cylindrical on 𝑣 →𝑡 such that 𝜑 𝑣 →𝑡 𝑖,𝑡 (𝑥 𝑡 ) = 𝜙 𝑣 𝑖,0 (𝑥 ←𝑡 𝑡 ) and 𝜑 ∅ 𝑖,𝑡 ∶= 𝜙 ∅ 𝑖,0 . Moreover, for any nonempty neighborhood 𝑣 in V 0 , we consider 𝜆 𝑖,𝑡 (𝑣 →𝑡 ) = 𝜆 𝑖,0 (𝑣). 𝑣 →𝑡 ∈V 𝑡 ,𝑣 →𝑡 ≠∅ 𝜆 𝑖,𝑡 (𝑣 →𝑡 ) = 1.

  1 , … , 𝛼 𝑘 ),

	𝑣 𝑖 𝛼 1∶𝑘 = 𝑣 𝑖 (𝛼 1 ,…,𝛼 𝑘 ) ∶= ∪ 𝑘 𝑙=1 𝑤 𝛼 𝑙 ,	(3.4.6)

where 𝑤 𝛼 𝑙 ∶= {𝑗} × [-(𝑛 + 1)𝜖, 𝑛𝜖) if 𝛼 𝑙 = (𝑗, 𝑛), and 𝑣 𝑖 0 ∶= ∅.

Proposition 5. Consider a non linear Hawkes process 𝑍 𝑖 , 𝑖 ∈ I, with intensity 𝜙 𝑖 (𝑥) given by

  𝑣∈V 𝑖 , which corresponds to the Kalikow decomposition of intensity 𝜙 𝑖,𝑇 (𝑥 𝑇 ) in (3.3.1) where 𝑥 𝑇 ∈ 𝒴 𝑇 .

More precisely, this means that we attach to (𝑖, 𝑇 ) a random variable 𝑉 𝑖,𝑇 with values in (V 𝑖 ) →𝑇 . For any 𝑣 ∈ V 𝑖 , P(𝑉 𝑖,𝑇 = 𝑣 →𝑇 ) = 𝜆 𝑖 (𝑣).

Proposition 6 .

 6 If we suppose that the process (𝒵 𝑖 𝑡 (1)) 𝑡∈R stays in (𝒴 𝑡 ) 𝑡∈R almost surely then it admits 𝜙 𝑖,𝑡 (𝑥 𝑡 ) as ℱ Proof. Take 𝐶 𝑡 a non negative predictable process with respect to ℱ

	𝒵 (1) 𝑡 -	-predictable intensity.
		𝒵 𝑖 (1) 𝑡	, that is ℱ 𝑡-𝒵 (1)	measurable and
	therefore ℱ 𝒵 𝑡-measurable . We have, for any 𝑖 ∈ I,

  Kalikow decomposition, it is more convenient to add 𝑉 𝑇 , the neighborhood picked at time 𝑇.

	Component	Terminology	Meaning
	𝑇 ∈ [0, 𝑡 𝑚𝑎𝑥 ]	Time definition	When the 𝑝𝑜𝑖𝑛𝑡 occurred
	𝑗 𝑇 ∈ I	Index definition	Where the 𝑝𝑜𝑖𝑛𝑡 occurred.
	𝑋 𝑇 ∈ {0, 𝑛.𝑎, 1}	Thinning mark	indicates if the 𝑝𝑜𝑖𝑛𝑡 is rejected (𝑋 𝑇 = 0), selected (𝑋 𝑇 = 1) or still unknown (𝑋 𝑇 = 𝑛.𝑎).
	𝑈 𝑇 ∈ [0, 1]	Uniform	Uniform random variable in [0, 1].
	𝑉 𝑇	Neighborhood	Random neighborhood associated to the 𝑝𝑜𝑖𝑛𝑡, only needed if we use Kalikow decomposition.

Table 4 .

 4 1: Definition of a 𝑝𝑜𝑖𝑛𝑡 A 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 is of index 𝑖 and at time 𝑇. Saying that a 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 is produced means that the simulation of the Poisson process Π 𝑖 has an event at time 𝑇. As soon as 𝑆 𝑖 𝑇 is produced, it is said to be unknown (𝑋 𝑇 = 𝑛.𝑎) and 𝑈 𝑇 is randomly set as an Uniform random variable on [0, 1].

  1 𝑇 and 𝒮 2 𝑇 . We say 𝒮 1 𝑇 ⪯ 𝒮 2 𝑇 if all the accepted point in 𝒮 1 𝑇 is also included and accepted in 𝒮 2 𝑇 , that is ∀𝑆 𝑖 𝑇 ∈ 𝒮 1 𝑇 , if 𝑋 𝑇 = 1 then 𝑆 𝑖 𝑇 ∈ 𝒮 2 𝑇 It is worth noting that 𝜓 𝑖 (𝑇 |𝒮 𝑇 ) is increasing with respect to the configuration of thinning marks on 𝒮 𝑇 , whenever 𝒮 𝑇 ⪯ 𝒮 ′ 𝑇 . It also implies that the intensity function 𝜙 𝑖 (𝑇 |𝒮 𝑇 ) is increasing with respect to the configuration of thinning marks on 𝒮 𝑇 . Knowing 𝜙 𝑖 (𝑇 |ℱ 𝑇 ) allows computing the probability to select the 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 using the thinning formula Then the thinning mark 𝑋 𝑇 can be defined as a realization of a random variable Bernoulli with parameter P(𝑋 𝑇 = 1|ℱ 𝑇 ) (formula 4.2.3).Once every 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 gets its final thinning mark, the 𝑠𝑝𝑖𝑘𝑒𝑠 train is obtained by keeping only the selected points.

	Pseudo-code
	Algorithm 4 Algorithm OtS (for Ogata's thinning Sequential)
	▷ Step 1: Seeding step
	1: for each 𝑛𝑒𝑢𝑟𝑜𝑛 𝑖 do
	2:
	that is
	𝜓(𝑇 |𝒮 𝑇 ) ≤ 𝜓(𝑇 |𝒮 ′ 𝑇 )
	P(𝑋 𝑇 = 1|ℱ 𝑇 ) =	𝜙 𝑖 (𝑇 |ℱ 𝑇 ) 𝑀 𝑖	.	(4.3.4)

P

  : Define the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠 : 𝒮 = ⋃ 𝒮 𝑖 10: Sort this set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 in an increasing 𝑇-order 11: for each point 𝑆 𝑖 𝑇 ∈ 𝒮 do The algorithm OtS cannot be parallelized without an adaptation. Indeed, at line 12 of Algorithm 4 to compute 𝜙 𝑖 (𝑇 |ℱ 𝑇 ), we need the past before time 𝑇, ℱ 𝑇 . This past is available only if the 𝑝𝑜𝑖𝑛𝑡𝑠 in 𝒮 are sorted in time and we go through point by point. In the case of parallel computation, the 𝑝𝑜𝑖𝑛𝑡𝑠 are at best partially sorted. Hence when a 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑇 is computed it may exist 𝑆 𝑇 ′ with 𝑇 ′ < 𝑇 where 𝑋 𝑇 ′ = 𝑛.𝑎. To deal with the problem in Remark 21, let us note that even if the intensity 𝜙 𝑖 (𝑇 |ℱ 𝑇 ) cannot be determined, it remains possible to provide nontrivial bounds for 𝜓 𝑖 (𝑇 |𝒮 𝑇 ), the part without the refractory period. Indeed, we already note (Remark 19) that 𝜓 𝑖 (𝑇 |𝒮 𝑇 ) is increasing with respect to any thinning mark of the 𝑝𝑜𝑖𝑛𝑡𝑠 produced before 𝑇. Hence it should be possible to build two sets 𝒮 𝑚𝑖𝑛 𝑇 and 𝒮 𝑚𝑎𝑥 𝑇 of 𝑝𝑜𝑖𝑛𝑡𝑠 all known (that is none with the thinning mark n.a) such that 𝜓 𝑖 (𝑇 |𝒮 𝑚𝑖𝑛 𝑇 ) ≤ 𝜓 𝑖 (𝑇 |𝒮 𝑇 ) ≤ 𝜓 𝑖 (𝑇 |𝒮 𝑚𝑎𝑥 𝑇 ). (4.3.5) For any 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑇 ′ ∈ 𝒮 𝑇 , we define the 𝑝𝑜𝑖𝑛𝑡𝑠 𝑆 𝑚𝑖𝑛 𝑇 ′ and 𝑆 𝑚𝑎𝑥 𝑇 ′ as the 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑇 ′ with the respective thinning mark 𝑋 𝑚𝑖𝑛 𝑇 ′ and 𝑋 𝑚𝑎𝑥 𝑇 ′ such that Our goal is to use the bounds for 𝜓 𝑖 (𝑇 |𝒮 𝑇 ) to determine bounds for the acceptance probability of the 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 . We are looking for P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|ℱ 𝑇 ) and P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|ℱ 𝑇 ) such thatP 𝑚𝑖𝑛 (𝑋 𝑇 = 1|ℱ 𝑇 ) ≤ P(𝑋 𝑇 = 1|ℱ 𝑇 ) ≤ P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|ℱ 𝑇 ).Since the relation (4.3.5) provides bounds for 𝜓 𝑖 (𝑇 |ℱ 𝑇 ) and not for 𝜙 𝑖 (𝑇 |ℱ 𝑇 ) we have to discuss three possible cases concerning the refractory period.i. Case 1 : The refractory period is not respected if ∃𝑆 𝑖 𝑇 ′ with 𝑇 -𝛿 ≤ 𝑇 ′ < 𝑇 and 𝑋 𝑇 ′ = 1. In this case 𝑆 𝑖 𝑇 is surely rejected (𝑋 𝑇 = 0):P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|ℱ 𝑇 ) = P(𝑋 𝑇 = 1|ℱ 𝑇 ) = P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|ℱ 𝑇 ) = 0.The respect of refractory period cannot be determined if ∃𝑆 𝑖 𝑇 ′ with 𝑇 -𝛿 ≤ 𝑇 ′ < 𝑇 and 𝑋 𝑇 ′ = 𝑛.𝑎. In this case : 𝑚𝑖𝑛 (𝑋 𝑇 = 1|ℱ 𝑇 ) = 0 P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|ℱ 𝑇 ) = The refractory period is respected if ∀𝑆 𝑖 𝑇 ′ with 𝑇 -𝛿 ≤ 𝑇 ′ < 𝑇 the thinning mark 𝑋 𝑇 ′ = 0. Hence 𝜙 𝑖 (𝑇 |ℱ 𝑇 ) = 𝜓 𝑖 (𝑇 |𝒮 𝑇 ) then

	8: end for ▷ Step 2: Thinning 12: Compute the intensity 𝜙 𝑖 (𝑇 |ℱ 𝑇 ) with formula (4.3.2) 13: Compute the probability P(𝑋 𝑇 = 1|ℱ 𝑇 ) = 𝜙 𝑖 (𝑇 |ℱ 𝑇 ) 𝑀 𝑖 14: Determine the thinning mark 𝑋 𝑇 of the 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 (4.2.3) to accept 𝑆 𝑖 𝑇 15: end for ▷ Final step : format the spikes trains 16: Delete the points that have a thinning mark 𝑋 𝑇 = 0 �.�.� Partially parallelized Ogata's thinning algorithm 𝑇 ′ = { 0 if 𝑋 𝑇 ′ = 𝑛.𝑎 𝑋 𝑇 ′ otherwise Remark 21. 𝑋 𝑚𝑖𝑛 ii. Case 3 : ⎧ ⎪ ⎨ ⎪ ⎩ 𝜓 𝑖 (𝑇 |𝒮 𝑚𝑖𝑛 𝑇 ) P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|ℱ 𝑇 ) = 𝑀 𝑖 P 𝑚𝑎𝑥 (𝑋 𝑇 = 1|ℱ 𝑇 ) = 𝜓 𝑖 (𝑇 |𝒮 𝑚𝑎𝑥 𝑀 𝑖 . 𝑇 )	(4.3.6) (4.3.11)
	Remark 20. The algorithm OtS stops after a finite number of steps almost surely. Indeed only a finite 𝑋 𝑚𝑎𝑥 𝑇 ′ = { 1 if 𝑋 𝑇 ′ = 𝑛.𝑎 𝑋 𝑇 ′ otherwise (4.3.7)
	number of 𝑝𝑜𝑖𝑛𝑡𝑠 are produced after the seeding step. The sets 𝒮 𝑚𝑖𝑛 𝑇 = ⋃ 𝑆 𝑚𝑖𝑛 𝑇 ′ and 𝒮 𝑚𝑎𝑥 𝑇 = ⋃ 𝑆 𝑚𝑎𝑥 𝑇 ′ satisfy the relation (4.3.5). (4.3.8)
				(4.3.9)
	ii. Case 2 : ⎧ ⎪ ⎨ ⎪ ⎩	𝜓 𝑖 (𝑇 |𝒮 𝑚𝑎𝑥 𝑀 𝑖 𝑇 )	.	(4.3.10)

9

  Compute the upper bound probabilityP 𝑚𝑎𝑥 (𝑋 𝑇 = 1|ℱ 𝑇 ) (4.3.8)Compute the lower bound probability P 𝑚𝑖𝑛 (𝑋 𝑇 = 1|ℱ 𝑇 ) (4.3.8)

	16:	for each unknown 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 ∈ 𝒮 𝑘 𝑝𝑜𝑝 do
	17:
	18:
	19:	Determine the thinning mark 𝑋 𝑇 if it is possible (4.3.12)
	20:	end for
	21:	end for
	22:

  any neuron 𝑖, suppose that the intensity 𝜙 𝑖 (𝑡|ℱ 𝑡 ) obeys a Kalikow decomposition for all 𝑡 ∈ [0, 𝑡 𝑚𝑎𝑥 ]. This means that, at each time 𝑡, there exists a neighborhood family V 𝑡 , in which, we define a probability distribution 𝜆 𝑖 (𝑣 𝑡 ) and functions 𝜙 𝑣 𝑡 𝑖 (𝑡|ℱ 𝑡 ) that depend only on the spikes in neighborhood 𝑣 𝑡 such that 𝜙 𝑖 (𝑡|ℱ 𝑡 ) = ∑ 𝑣 𝑡 ∈V 𝑡 𝜆 𝑖 (𝑣 𝑡 )𝜙 𝐼 ⋃{∅} and 𝑛 ∈ N. With this concrete definition, the associated intensity 𝜙

	𝑣 𝑡 𝑖 (𝑡|ℱ 𝑡 ).
	To illustrate our result, we focus on a very specific choice of neighborhood 𝑣 𝑇	1 , which is given by:
	𝑣 𝑇 = {𝑗} × [𝑇 -(𝑛 + 1)𝛿, 𝑇 -𝑛𝛿)	(4.3.13)
	with 𝑗 ∈ 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) is as
	follows

Remark 28 .

 28 𝑀 𝑖 , the upper bound of the intensity 𝜙 𝑖 (𝑇 |ℱ 𝑇 ) is not an upper bound of the intensity 𝜙 𝑉 𝑇 𝑖 (𝑇 |ℱ 𝑇 ). However, as mentioned in Chapter 2 in the case of the Hawkes process with refractory period, it remains possible to provide an upper bound of 𝜙 𝑉 𝑇 𝑖 (𝑇 |ℱ 𝑇 ). (𝑇 |ℱ 𝑇 ) depends on the interaction functions and the chosen decomposition. For the simulation purpose, let us assume that the upper bound exist. For any 𝑖, 𝑇 and neighborhood 𝑣 𝑇 ∈ V 𝑇 , there exists a deterministic positive number 𝑀 𝑣 𝑇 such that 𝜙 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) ≤ 𝑀 𝑣 𝑇 .

	The upper bound of 𝜙	𝑉 𝑇
	Assumption 9. i. ii. Grant the first hypothesis (i.), there exists moreover a deterministic positive number 𝑀 𝑖 such that
		sup	sup
		𝑇	𝑣 𝑇 ∈V 𝑇

𝑖

𝑀 𝑣 𝑇 ≤ 𝑀 𝑖 .

Remark 29.

Here for each neighborhood 𝑣 𝑇 , we suppose that there exists a deterministic upper bound 𝑀 𝑣 𝑇 for 𝜙 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ). This hypothesis eventually guarantees to do thinning on each neighborhood 𝑣 𝑇 .

  .3.15) Remark 31. Note that, in our settings in Section 4.3.5, 𝜆 𝑠𝑝𝑎𝑐𝑒 ({𝑗}) and 𝜆 𝑡𝑖𝑚𝑒 ({𝑛}) can be computed independently of 𝑇.As this algorithm is based the Ogata thinning method (Section 4.3.1), it is also written in two main steps.During the first step, as for the algorithm OtS, each 𝑛𝑒𝑢𝑟𝑜𝑛 𝑖 ∈ I produced the set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 𝑖 . In particular , the vector (𝑇 𝑖 𝑘 ) 𝑘∈[|1,𝑛 𝑖 |] corresponds to the time arrivals of a homogeneous Poisson process with intensity 𝑀 𝑖 .

	Overall description

  𝑖 𝑇 in the set 𝒮, the intensity 𝜙 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) is computed using the formula defined in the section about Kalikow decomposition(4.3.14). The algorithm split this formula to isolate the Hawkes process intensity 𝜓 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) from the refractory period constraint: 𝑇 -𝐿 𝑍 𝑖(𝑇 )>𝛿 = 0, the refractory period is not respected, 𝜙 (𝑇 |ℱ 𝑇 ). In the opposite case, if 1 𝑇 -𝐿 𝑍 𝑖(𝑇 )>𝛿 = 1 then 𝜙 (𝑇 |ℱ 𝑇 ) only depends on 𝒮 𝑣 𝑇 , the set of the 𝑝𝑜𝑖𝑛𝑡𝑠 produced inside the neighborhood 𝑣 𝑇 .

	with 𝜓	𝑣 𝑇 𝑖 defined as follows.
			𝜓	𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) =	⎧ ⎪ ⎨ ⎪ ⎩	𝜓 ∅ 𝑖 (𝑇 |ℱ 𝑇 ) = 1 𝜆 𝑖 (𝑣 𝑇 ) ∫ 𝑇 -𝑛𝛿 𝑇 -(𝑛+1)𝛿 ℎ 𝑗𝑖 (𝑇 -𝑡)𝑑𝑍 𝜇 𝑖 𝜆 𝑖 (𝑣 𝑇 ) if 𝑗 = ∅	𝑗 𝑡 if 𝑗 ≠ ∅.	(4.3.18)
	If 1 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) = 0 without the need to compute 𝜓 𝑣 𝑇 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) = 𝜓 𝑣 𝑇
		𝜓	𝑣 𝑇 𝑖
				𝜙	𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) = 𝜓	𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 )1 𝑇 -𝐿 𝑍 𝑖(𝑇 )>𝛿	(4.3.17)

𝑖 𝑖

  .3.19) Remark 32. Again, similarly to Remark 19, we observe that 𝜓 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑣 𝑇 ) is increasing with respect to the configuration of thinning marks in 𝑣 𝑇 . In the case where 𝑗 ≠ ∅, since the time neighborhood [𝑇 -(𝑛 + 1)𝛿, 𝑇 -𝑛𝛿] has a length equals to a refractory period 𝛿, only one term of the sum is a non zero term. In other words, 𝜓 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑣 𝑇 ) depends on at most one previous event. It was mentioned within Remark 27 as a condition to make the execution time linear with respect to the number of 𝑛𝑒𝑢𝑟𝑜𝑛𝑠. The algorithm can then compute the probability to select the 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 using the thinning formula with the intensity 𝜙 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) : P(𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 ) = 𝜙

	Remark 33. 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) 𝑀 𝑖	.	(4.3.20)

: Kalikow step and pre-thinning

	3: 4:	for each 𝑝𝑜𝑖𝑛𝑡 indexed by 𝑘 ≤ 𝑛 𝑖 do Simulate the time : 𝑇 𝑖 𝑘 follows a 𝑈 ([0, 𝑡 𝑚𝑎𝑥 ]) distribution
	5:	Define the uniform : 𝑈 𝑇 𝑖 𝑘	follows a 𝑈 ([0, 1]) distribution
	▷ Step 26: Pick neighborhood in space : {𝑗} wp 𝜆 𝑠𝑝𝑎𝑐𝑒 𝑖	({𝑗})
	7: 8: 9:	which only depends on 𝑗 Reject the 𝑝𝑜𝑖𝑛𝑡 : the thinning mark 𝑋 𝑇 𝑖 Compute the upper bound 𝑀 𝑣 𝑇 𝑖 𝑘 if 𝑈 𝑇 𝑖 𝑘 > 𝑀 𝑣 𝑇 𝑖 𝑘 then 𝑀 𝑖 𝑘 = 0
	10:	else
	11: 12:	Set the 𝑝𝑜𝑖𝑛𝑡 as unknown : the thinning mark 𝑋 𝑇 𝑖 𝑘 Pick neighborhood in time : [𝑇 𝑖 𝑘 -(𝑛 + 1)𝛿, 𝑇 𝑖 𝑘 -𝑛𝛿] wp 𝜆 𝑡𝑖𝑚𝑒 = 𝑛.𝑎 𝑖	({𝑛})
	13:	end if
	14:	end for
	15: end for
		▷ Step 3 : Thinning
	19: 21:	𝑣 𝑇 Determine the thinning mark 𝑋 𝑇 of the 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 Compute the intensity 𝜙 𝑇 (4.2.3) 𝑣 𝑇 𝑖 (𝑇 |ℱ 𝑇 ) 𝑀 𝑖	to accept 𝑆 𝑖 𝑇
	22:

16: Define the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠 : 𝒮 = ⋃ 𝒮 𝑖 17: Sort this set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 in an increasing 𝑇-order 18: for each point 𝑆 𝑖 𝑇 ∈ 𝒮 do 𝑖 (𝑇 |ℱ 𝑇 ) with formula (4.3.17) 20:

Compute the probability P(𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇 ) = 𝜙

end for ▷ Final step : format the spikes trains 23

  : Delete the points that have a thinning mark 𝑋 𝑇 = 0 The problem is very similar to the case of algorithms without Kalikow decomposition.We follow the same pattern as in Section 4.3.2 in which we proposed a parallelized version of Ogata's thinning algorithm. The main issue was to define bounds for the intensity if it cannot be exactly computed.By Remark 33 we know that a unique 𝑝𝑜𝑖𝑛𝑡 of the neighborhood will be selected at the end. This remark simplifies the procedure. First, if a 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑗 𝑇 ′ has already been selected we know that it is the only one and that the other 𝑝𝑜𝑖𝑛𝑡𝑠 will be rejected. Hence 𝜓 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑣 𝑇 ) can be exactly computed as follows Otherwise there is no accepted 𝑝𝑜𝑖𝑛𝑡 in 𝒮 𝑣 𝑇 , in this case, the lower bound of 𝜓 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑇 ) is 0. It remains to find the upper bound. We define a set 𝒮 With the same remark as before we can build 𝒮 𝑣 𝑇 ,𝑚𝑎𝑥 𝑇 as a set of a unique 𝑝𝑜𝑖𝑛𝑡, the one not rejected that maximize the interaction function ℎ 𝑗𝑖 , with 𝑗 fixed by 𝑣 𝑇 . Hence the bounds are as follows

	Remark 34. For a general choice of neighborhood, it suffices to replace Step 2 by choosing a neighborhood 𝑣 𝑇 𝑖 𝑘 with a general Kalikow decomposition and do pre-thinning for each 𝑇 𝑖 𝑘 . �.�.� Partially parallelized Kalikow-Ogata's thinning algorithm 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑇 ) = ℎ 𝑗𝑖 (𝑇 -𝑇 ′ ) 𝜆 𝑖 (𝑣 𝑇 ) . (4.3.21) 𝑣 𝑇 ,𝑚𝑎𝑥 𝑇 such that 0 ≤ 𝜓 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑇 ) ≤ 𝜓 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑣 𝑇 ,𝑚𝑎𝑥 𝑇 ). (4.3.22) 0 ≤ 𝜓 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑇 ) ≤ 1 𝜆 𝑖 (𝑣 𝑇 ) max 𝑆 𝑗 𝑇 ′ ∈𝒮 𝑣 𝑇 ℎ 𝑗𝑖 (𝑇 -𝑇 ′ )1 𝑋 𝑇 ′≠0 . (4.3.23) Remark 35. 𝜓 Then we use the bounds for 𝜓 𝑣 𝑇

Let us start with the intensity without considering the refractory period. We are looking for bounds for 𝜓 𝑣 𝑇 𝑖 (𝑇 |𝒮 𝑣 𝑇 ) if there exist a 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑗 𝑇 ′ ∈ 𝒮 𝑣 𝑇 with the thinning mark 𝑋 𝑇 ′ = 𝑛.𝑎. 𝑖 (𝑇 |𝒮 𝑇 ) to build bounds for P(𝑋 𝑇 = 1|𝑉 𝑇 = 𝑣 𝑇 , ℱ 𝑇

  Note that, in our setting, P 𝑚𝑖𝑛 stay to be 0 until one point in 𝑣 𝑇 is said to be accepted. In this case, P 𝑚𝑖𝑛 = P 𝑚𝑎𝑥 = P. With more general choices of neighborhood, it suffices to accept/reject all unknown points in neighborhood 𝑣 𝑇 to construct P 𝑚𝑖𝑛 , P 𝑚𝑎𝑥 , as long as the intensities are monotone with respect to 𝒮 𝑇 . Again, the crucial hypothesis in this section is the existence of P 𝑚𝑖𝑛 and P 𝑚𝑎𝑥 . This is guaranteed as soon the intensity restricted to any neighborhood 𝑣 𝑇 is increasing with respect to the configuration of thinning marks in neighborhood 𝑣 𝑇 , that is𝜙 𝑣 𝑇 (𝑇 |𝒮

	Remark 37. 1. 2. In addition, with more general choices of neighborhood, it suffices to replace the Step 2 by choosing
	a neighborhood 𝑣 𝑇 𝑖 𝑙	by Kalikow decomposition and do pre-thinning for each 𝑇 𝑖 𝑙 .
	3. 𝑣 𝑇 𝑇 ) ≤ 𝜙 𝑣 𝑇 (𝑇 |𝒮 ′ 𝑇	𝑣 𝑇 )
	whenever 𝒮	𝑣 𝑇 𝑇 ⪯ 𝒮 ′ 𝑇	𝑣 𝑇 .

.3.25) Remark 36.

  Note that 𝜙 𝑣 (.| F𝑡 ) is the cylindrical function that only depends on the points on the neighborhood 𝑣.Let us now integrate with respect to the choice 𝑉 𝑇 𝑖 𝑘 , which is independent of anything else. Since for any 𝑇, Λ 𝑖 (𝑇 ) = 𝑀 𝑖 (𝐿 Π (𝑇 )), we conclude that by definition of Kalikow decomposition for 𝜙 𝑖 (𝑡| F𝑡 ) and the fact that Π 𝑖 is empty outside the interval [0, 𝜏]. Moreover, since Π 𝑖 𝑡 is constructed as point process with F𝑡 intensity Λ 𝑖 (𝑡), it leads us to conclude that 𝐶 𝑡 𝜙 𝑖 (𝑡| F𝑡 )𝑑𝑡 ) .

	E (∫ 0	𝜏	𝐶 𝑡 𝑑Π 𝑖 𝑡 (1) )	= ∑ 𝑘	E	⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝐶 𝑇 𝑖 𝑘	∑ 𝑣 𝑇 𝑖 𝑘 ∈V 𝑇 𝑖 𝑘	𝜆 𝑖 (𝑣 𝑇 𝑖 𝑘 Λ 𝑖 (𝑇 𝑖 )𝜙 𝑣 𝑇 𝑖 𝑘 𝑖 (𝑇 𝑖 𝑘 | F𝑇 𝑖 𝑘 𝑘 )	)	1 𝑇 𝑖 𝑘 ≤𝜏	⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
											= E (∫ 0	𝜏	𝐶 𝑡	𝜙 𝑖 (𝑡| F𝑡 ) Λ 𝑖 (𝑡)	𝑑Π 𝑖 𝑡 )	.
							E (∫ 0	𝜏	𝐶 𝑡 𝑑Π 𝑖 𝑡 (1) )	= E (∫ 0	𝜏
						(∫ 0	𝜏	𝐶 𝑡 𝑑Π 𝑖 𝑡 (1) )	= ∑ 𝑘∈N	E (	𝐶 𝑇 𝑖 𝑘	1 𝑋 𝑇 𝑖 𝑘	=1 1 𝑇 𝑖 𝑘 ≤𝜏 )	.
	Conditioning on {𝑉 𝑇 𝑖 𝑘	, F𝑇 𝑖 𝑘	}, we have
	E (∫ 0	𝜏	𝐶 𝑡 𝑑Π 𝑖 𝑡 (1) )	= ∑ 𝑘	E (	𝐶 𝑇 𝑖 𝑘	E ( 1 𝑋 𝑇 𝑖 𝑘	=1 1 𝑇 𝑖 𝑘 ≤𝜏 |𝑉 𝑇 𝑖 𝑘	, F𝑇 𝑖 𝑘 ))
												= ∑ 𝑘	E ⎛ ⎜ ⎜ ⎜ ⎝ 𝐶 𝑇 𝑖 𝑘	𝜙 𝑀 𝑖 (𝐿 Π (𝑇 𝑖 𝑉 𝑇 𝑖 𝑘 𝑖 (𝑇 𝑖 𝑘 | F𝑇 𝑖 𝑘 𝑘 )) )	1 𝑇 𝑖 𝑘 ≤𝜏	⎞ ⎟ ⎠ ⎟ ⎟ .

Event-Scheduling Algorithms with Kalikow Decomposition for Simulating Potentially Infinite Neuronal Networks

Note that, we write this neighborhood notation here informally. In fact, 𝑣 𝑇 in the previous chapter equal to 𝑣 𝑇 ∪{𝑖}×(𝑇 -𝛿, 𝑇 )

New methods for simulating point processes.Pseudo-code

Acknowledgements

When I write these lines, I know that it marks the end of a long and memorable journey in my life . I am grateful to many people for the help and support during the various stage of writing my PhD thesis. First and foremost, I would like to thank you the defense committee: Philippe Robert, Nancy Lopes Garcia, Etienne Tanre, and Alexandre Muzy for crucial feedbacks which greatly improve my thesis. Manon and Roland played a crucial role in my committee de these over 3 years, providing endless assistance and for this, I am grateful to them.

Next and certainly, most importantly, I turn to Patricia and Eva. Surely a page, even a book, is not enough for me to express my gratitude to them. No words are enough to express my gratitude. I always consider them as my second mothers. They give me advice, both in math and in real life. I have known Patricia for 5 years. Perhaps no one understands my abilities better than Patricia. Patricia taught me so much and spent a lot of time talking to me in her limited time. In Patricia, I found constant creativity, enthusiasm for science, and daring to take on challenges. She also taught me lessons about interaction with other colleagues. Although I don't have much time to work directly with Eva, in return, every time I get to meet Eva, she always takes precious time to talk to me. Her extensive knowledge has helped me a lot to get to where I am today.

I would like to thank the friends who have helped me a lot along this way. First, I would like to express Acknowledgements This work was supported by the French government, through the UCA 𝐽 𝑒𝑑𝑖 Investissements d' Avenir managed by the National Research Agency (ANR-15-IDEX-01) and by the interdisciplinary Institute for Modeling in Neuroscience and Cognition (NeuroMod) of the Université Côte d' Azur. The authors would like to thank Professor E.Löcherbach from Paris 1 for great discussions about Kalikow decomposition and Forward Backward Algorithm. Acknowledgements This paper is dedicated to the Institute of Mathematics Hanoi, where the author was honored to work for more than three years. The author would like to express his very great appreciation to his supervisors Patricia Reynaud-Bouret and Eva Löcherbach for their patient guidance, enthusiastic encouragement and useful critiques on this research. This work was supported by the French government, through the UCAJedi and 3IA Côte d' Azur Investissements d' Avenir managed by the National Research Agency (ANR-15-IDEX-01 and ANR-19-P3IA-0002), by the CNRS through the "Mission pour les Initiatives Transverses et Interdisciplinaires" (Projet DYNAMO, "APP Modélisation du Vivant"), by the interdisciplinary Institute for Modeling in Neuroscience and Cognition (NeuroMod) of the Université Côte d' Azur, and directly by the National Research Agency (ANR-19-CE40-0024) with the ChaMaNe project.

Introduce 𝑊 𝑖 (∞) ∶= lim 𝑛→∞ 𝑊 𝑖 (𝑛), if it exists, and let 𝑊 𝑖 = 1 𝑇 𝑊 𝑖 (∞) be the total number of points of C . The log-Laplace transform of 𝑊 𝑖 (∞) is given by

) ,

and we write again Φ(𝜃) = (Φ 1 (𝜃), Φ 2 (𝜃), … , Φ 𝑁 (𝜃)) 𝑇 .

With this approach, we can prove the following exponential inequality. This result is inspired by Lemma 1 of [21]. Define ‖.‖ 1 a 1norm on R 𝑁 and the ball with respect to this norm on R 𝑁 by 𝐵(.).

Proposition 8. Grant Assumption 6 and assume that the weights (𝜆 𝑖 (.)) 𝑖∈I satisfy the following condition:

Suppose moreover that there exists a positive number 𝑟 that depends on the matrix 𝑀 such that, for all positive vectors 𝜃 belonging to 𝐵(0, 𝑟), we have

Then for all 𝜃 ∈ 𝐵(0, 𝑟), we conclude that Φ 𝑖 (𝜃) < ∞ and moreover Φ(𝜃) = 𝜃 + 𝜙(Φ(𝜃)).

In particular, for 𝜗 ∈ R + such that 𝜃 = 𝜗1 ∈ 𝐵(0, 𝑟), there exists a constant 𝑐 0 that depends on 𝑀 and 𝑖 such that

for all 𝑥 > 0.

Remark 16. From Chapter V of [1] , we have

� New methods for simulating point processes.

This chapter is written in collaboration with Paul Gresland during his master internship, under supervision of Patricia Reynaud-Bouret and Alexandre Muzy.

Status: working progress.

Mathematical definitions and notation

�.�.� Deterministic and stochastic upper bound of intensities

and we also consider that 𝐿 𝑍 (𝑡) = sup{𝑇 ∈ 𝑍 such that 𝑇 < 𝑡} the last arrival time before time 𝑡 of the process 𝑍, with the convention that 0 is the last arrival time before the first event of any process 𝑍.

The second assumption consists in saying that the intensity now is bounded by a particular stochastic predictable function. In this case, we construct a dominating intensity point after point, gradually starting from a initial point 𝑇 = 0. I remove the measuring in the following assumption since we have not constructed the extended filtration. We have not constructed Π and don't have the notations 𝑋.

Assumption 8 (Stochastically bounded intensity

�.� First part: Deterministically bounded intensities and Parallelization

This section is devoted to present four algorithms: two based on Ogata's thinning methods and the two other methods being a combination of Ogata's algorithm and Kalikow decomposition. To simplify, we write our algorithms in a particular example: the Linear Hawkes processes with refractory period [3,20].

The general algorithms can be constructed very similarly, we give more details in the remarks at the end of each algorithm (if needed). Recall that the intensity of Linear Hawkes process with hard refractory period is given by

where recall that 𝐿 𝑍 𝑖 (𝑡) = 𝑠𝑢𝑝{𝑇 ∈ 𝑍 𝑖 such that 𝑇 < 𝑡}, 𝜇 𝑖 's positive real numbers and ℎ 𝑗𝑖 the non negative functions with support in R + .

We begin with two modifications of Ogata's algorithms. The first one called OtS, is the classical Ogata's thinning algorithm in [16]. The second, called OtpP, is a partially parallelized version of OtS.

�.�.� Sequential Ogata's thinning algorithm

Ogata's thinning algorithm

The modificated Ogata's thinning method aims to simulate all the processes 𝑖 ∈ I in two steps: seeding step and thinning step:

i. Seeding step : For each index 𝑖 ∈ I, a list of candidate points is built by simulating Π 𝑖 a homogeneous

Poisson process with intensity 𝑀 𝑖 defined in Assumption 7.

ii. Thinning step: Each of the candidate points are selected or rejected via an exact computation of the intensity at the corresponding time.

The algorithm presented in this section is a rewriting of Ogata's algorithm [16] in accordance with our programming method.

Remark 18. This method is indeed a modification of Ogata's algorithm [16] since in the original paper, Ogata constructs a dominating intensity process point after point before performing the thinning step.

Throughout this chapter, whenever we say Ogata's algorithm, it refers to this modification of Ogata's algorithm, and not the original one in [16].

Overall description

As this algorithm is essentially the Ogata thinning method (Section 4.3.1), it is also written in two main steps.

During the first step each 𝑛𝑒𝑢𝑟𝑜𝑛 𝑖 ∈ I produced the set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 𝑖 . In particular, the vector (𝑇 𝑖 𝑘 ) 𝑘∈[|1,𝑛 𝑖 |] corresponds to the times arrival of a homogeneous Poisson process with intensity 𝑀 𝑖 . Then the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 is built as the unions of all sets 𝒮 𝑖 . The set 𝒮 is sorted to finally contain all the 𝑝𝑜𝑖𝑛𝑡𝑠 in an increasing 𝑇-order, that is 𝑆 𝑖 𝑇 is "smaller" than 𝑆 𝑗 𝑇 ′ if and only if 𝑇 < 𝑇 ′ .

During the second step, for each 𝑝𝑜𝑖𝑛𝑡 𝑆 𝑖 𝑇 in the set 𝒮, the intensity 𝜙 𝑖 (𝑇 |ℱ 𝑇 ) is computed using the Hawkes formula with refractory period (4.3.26). The algorithm splits this formula to isolate the Hawkes process intensity 𝜓 𝑖 (𝑇 |ℱ 𝑇 ) from the refractory period constraint:

where for each 𝑛𝑒𝑢𝑟𝑜𝑛 𝑖 ∈ population do 3:

Compute the number of 𝑝𝑜𝑖𝑛𝑡𝑠 : 𝑛 𝑖 follows a 𝑃 𝑜𝑖𝑠𝑠𝑜𝑛(𝑀 𝑖 𝑡 𝑚𝑎𝑥 ).

4:

for each 𝑝𝑜𝑖𝑛𝑡 indexed by 𝑘 ≤ 𝑛 𝑖 do for each 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑘 in sequence do 30:

Update the whole set of 𝑝𝑜𝑖𝑛𝑡𝑠 𝒮 with the thinning marks 31:

end for 32: end while ▷ Final step : format the spikes trains 33: Delete the points that have a thinning mark 𝑋 𝑇 = 0 with 𝑋 𝑇 is the thinning random variable corresponds to 𝑇. Because F includes the accepted points, one can see F as an extension of ℱ. In addition, with this extended filtration, the intensity of Linear Hawkes process can be written as

which is F𝑡 predictable. Similarly for Linear Hawkes process with bounded support defined in (4.4.2). In addition, we define the first time where the system has 𝑛 0 points.

Note that the dominating process Π will be constructed on [0, 𝜏], outside of this interval, its value is 0.

Similarly to the deterministically bounded intensity, the proposed algorithm consists of two steps: seeding step and thinning step.

1. At the seeding step, we create a candidate to find the next point after time 𝑇. At this stage, we compute the dominating intensity 𝑀(𝑇 ) = ∑ 𝑁 𝑖=1 𝑀 𝑖 (𝑇 ) at time 𝑇, that creates a candidate 𝑇 ′ by adding an exponential random variable of parameter 𝑀(𝑇 ). We attribute a mark to this received candidate. This candidate 𝑇 ′ belongs to neuron 𝑖 with probability 𝑀 𝑖 (𝑇 )/𝑀(𝑇 ).

2. At the thinning step, assuming that it belongs to neuron 𝑖, we then assign to it a thinning mark with

Bernoulli distribution of parameter P(𝑋 𝑇

Remark 38. Note that, in order to find the next spike after time 𝑇, we use 𝑀 𝑖 (𝑇 )/𝑀(𝑇 ) and 𝜙 𝑖 (𝑇 ′ | F𝑇 ′ )/𝑀 𝑖 (𝑇 ).

They are F𝑇 measurable, in other words, we only need the information up to time 𝑇 and 𝑇 included.

The thinning rate at the candidate point of time 𝑇 ′ is given as follow:

Remark 39.

1. In comparison to previous section, where we start by simulating a Poisson process on each neuron, here the dominating intensity needs to be constructed gradually, point after point.

This prevents us to do parallel computation.

2. Note that, this algorithm is a modification of the original paper of Ogata [16]. In the original paper, Ogata do labelling task. More precisely, at each candidate point, he assigns a mark from 0 to 𝑁 where 0 means rejection. To do this task, we need to compute the convolutional sum of the intensities. This idea also appeared in [14]. However, when the number of neuron is huge, compute the convolutional sum of the intensities is very expensive. Whereas, in our algorithm, we only need to compute 𝑀 𝑖 (𝑇 ) for all 𝑖.

Pseudo-code

Algorithm 8 Algorithm OtSS for Ogata thinning sequential with stochastically bounded intensity 1: Initialize 𝑇 = 0.

2: Initialize the number of spikes 𝑛 = 0. Assume that, this candidate belongs to neuron 𝑖.

8:

Compute the thinning rate P(𝑋

Determine the thinning mark 𝑋 𝑇 ′ by Equation (4.4.5). If accepted, increase 𝑛 by 1.

It means

is well defined as a probability.

By using Proposition 1 in [16], we can prove that the simulated points after Step 1 is a multivariate point process Π = (Π 𝑖 ) 𝑖 having conditional intensity (Λ 𝑖 (𝑡)) 𝑖 up to stopping time 𝜏. We can prove that Proposition 10. Conditioning on the event {𝜏 < ∞}, the point process Π 𝑖 (1) obtained points by Algorithm 9 admits 𝜙 𝑖 (𝑡| F𝑡 ) as F𝑡 -predictable intensity up to the stopping time 𝜏.

Proof. This proposition is a particular case of Proposition 11. We accept it for a moment.

�.�.� Sequential KalikowOgata algorithm Overall description

For any neuron 𝑖, suppose that the intensity 𝜙 𝑖 (𝑡|ℱ 𝑡 ) satisfies a Kalikow decomposition for all 𝑡. This means that, at each time 𝑡, there exists a neighborhood family V 𝑡 , in which, we define a probability distribution 𝜆 𝑖 (𝑣 𝑡 ) and cylindrical functions 𝜙

�.� Numerical results

In this section, we prove statistically that the algorithms in Section 4.3 return Hawkes processes with the right intensity and assert that the new algorithms outperform the classical Ogata's algorithm [16]. The simulation is performed with a laptop computer with a processor Intel(R) Core(TM) i7-1165G7, base frequency 2,80 GHz, RAM 16 Go, and 4 hearts. The algorithms were implemented in Python programming language (version 3.8.3). The algorithms were run in parallel by using package multiprocessing in Python (version 3.8.3). The plots and statistical analysis were obtained by using Python (version 3.8.3), part of it using library statsmodels (v0.13.2).

�.�.� Statistical test

In 1988, Ogata [17] introduced a standard method to assess if the data obey a point process with a given intensity, in particular Hawkes processes. This method is based on the time rescaling theorem (see for instance [2]), which says that if 𝜆(.) is the conditional intensity of the point process 𝑍 and Λ(𝑡) = ∫ 𝑡 0 𝜆 𝑠 𝑑𝑠, then the point process Z = {Λ(𝑇 ), 𝑇 ∈ 𝑍} is the homogeneous Poisson process with rate 1. Ogata then derived a method to test a point process 𝑍 has a given intensity 𝜆(.) as follows:

• Applying the transform formula Λ to each arrival time of 𝑍 to obtain Z .

• Test 1: Test that the consecutive delay between points of Z obeys the exponential distribution of rate 1, by using the Kolmogorov-Smirnov test.

• Test 2: Test that the points of Z are uniformly distributed by using the Kolmogorov-Smirnov test.

• Test 3: Test that the delays between points of Z are independent, for instance, we perform an autocorrelation test between delays with a certain lag, here we choose up to lag 9.

We perform a test with the 200 neurons that are simulated during 5s (i.e 𝑡 𝑚𝑎𝑥 = 5). We consider the Hawkes process with a hard refractory period and bounded support with intensity is given by:

with 𝜇 𝑖 = 1∀𝑖, 𝐴 = 0. Remark 43. The OtpP algorithm is extremely slow, even worse than OtS, hence it is out of interest in practice. Therefore, we will not present its result in this section. On Figure 4.2, we can observe that the vertical lines stay on the blue shaded region. This means that there is no significant correlation between the delays. the number of 𝑛𝑒𝑢𝑟𝑜𝑛𝑠, the larger the number of 𝑝𝑜𝑖𝑛𝑡𝑠 unknown after the first iteration. This can be explained by the fact that when the number of 𝑛𝑒𝑢𝑟𝑜𝑛𝑠 is small the 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠 are small. Therefore 𝑛𝑒𝑢𝑟𝑜𝑛𝑠 are more likely to have important interactions with 𝑛𝑒𝑢𝑟𝑜𝑛𝑠 from other 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠. In other words, the larger a 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛, the more independent it is. We can see it well for 4 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠 and several 𝑛𝑒𝑢𝑟𝑜𝑛𝑠 between 10 and 100. In this interval, the execution time decreases with respect to the number of 𝑛𝑒𝑢𝑟𝑜𝑛𝑠.

Remark 44. However, it is worth noting that, Python is not suitable for parallel computing. There are 2 paradigms of parallel computing: multithreading and multiprocessing. There exist a mechanism in Python called "Global Interpreter lock" or GIL that prevents two or more thread from executing simultaneously. Therefore it is impossible to benefit from multithreading in Python. For multiprocessing, it is possible to run two or more Python processes in parallel, as each process has its interpreter. However, each process exists in its memory stack and is thus unable to exchange data unless an external synchronization mechanism (for example, MPI) is employed. To conclude, Paul Gresland and the engineers in NeuroMod and Inria (Sophia Antipolis) are working together to recode the algorithms in C++ to benefit the parallel computing.