Anna Bonnet

Miguel Martinez Herrera

Maxime Sangnier

Maximum Likelihood Estimation for Hawkes Processes with self-excitation or inhibition

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The Hawkes model is a point process observed on the real line, which generally corresponds to the time, where any previously encountered event has a direct influence on the chances of future events occurring. This past-dependent mathematical model was introduced in [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF] and its first application was to model earthquakes occurrences [START_REF] Ogata | Statistical models for earthquake occurrences and residual analysis for point processes[END_REF][START_REF] Ogata | Space-time point-process models for earthquake occurrences[END_REF]. Since then, Hawkes processes have been widely used in various fields, for instance finance [START_REF] Bacry | Scaling limits for Hawkes processes and application to financial statistics[END_REF], social media [START_REF] Rizoiu | A tutorial on Hawkes processes for events in social media[END_REF][START_REF] Mishra | Feature driven and point process approaches for popularity prediction[END_REF], epidemiology [START_REF] Rizoiu | SIR-Hawkes: linking epidemic models and Hawkes processes to model diffusions in finite populations[END_REF], sociology [START_REF] Linderman | Discovering latent network structure in point process data[END_REF] and neuroscience [START_REF] Reynaud-Bouret | Goodness-of-fit tests and nonparametric adaptive estimation for spike train analysis[END_REF].

The main advantage of Hawkes processes is their ability to model different kinds of relationships between phenomena through an unknown kernel or transfer function. The Hawkes model was originally introduced as a self-exciting point process where the appearance of an event increases the chances of another one triggering. Several estimation procedures have been proposed for the kernel function, both in parametric [START_REF] Ogata | Statistical models for earthquake occurrences and residual analysis for point processes[END_REF][START_REF] Fonseca | Hawkes process: Fast calibration, application to trade clustering, and diffusive limit[END_REF][START_REF] Ozaki | Maximum likelihood estimation of Hawkes' self-exciting point processes[END_REF] and nonparametric [START_REF] Reynaud-Bouret | Goodness-of-fit tests and nonparametric adaptive estimation for spike train analysis[END_REF][START_REF] Bacry | Second order statistics characterization of Hawkes processes and nonparametric estimation[END_REF] frameworks.

However, the inhibition setting, where the presence of an event decreases the chance of another occurring, has drawn less attention in the literature, although it can be of great interest in several fields, in particular in neuroscience [START_REF] Reynaud-Bouret | Reconstructing the functional connectivity of multiple spike trains using Hawkes models[END_REF]. In this inhibition context, the cluster representation [START_REF] Hawkes | A cluster process representation of a self-exciting process[END_REF] on which is based the construction of a self-exciting Hawkes process, is no longer valid. While the existence and the construction of such nonlinear processes can be found in recent works for the univariate [START_REF] Costa | Renewal in Hawkes processes with self-excitation and inhibition[END_REF] and multivariate [START_REF] Chen | The multivariate Hawkes process in high dimensions: Beyond mutual excitation[END_REF] cases, statistical estimation of the kernel function has been hardly addressed. A first approach consists in computing an approximation of the likelihood as if the intensity function could take negative values, and optimizing it to get a maximum likelihood estimator [START_REF] Lemonnier | Nonparametric markovian learning of triggering kernels for mutually exciting and mutually inhibiting multivariate Hawkes processes[END_REF]. Alternatively, the type of interaction (excitation or inhibition) can be considered as a hidden variable, giving rise to a very practical estimation method [START_REF] Mei | The neural Hawkes process: A neurally self-modulating multivariate point process[END_REF].

In this paper, we propose a maximum likelihood procedure that can handle both excitation and inhibition scenarios for a univariate Hawkes process. Our approach is based on an explicit computation of the likelihood for any type of monotone kernel functions, which is facilitated by the introduction of the natural concept of restart points. The latter are the times when the intensity function, that can be null on some intervals, become strictly positive again. We show that these restart points have a closed-form expression when the kernel is exponential, which allows us to rewrite and maximize the likelihood without approximations that are proposed for instance in [START_REF] Lemonnier | Nonparametric markovian learning of triggering kernels for mutually exciting and mutually inhibiting multivariate Hawkes processes[END_REF]. Our estimator is implemented in Python (the code is freely available online 1 ). We also propose a numerical study which shows the good performance of our exact estimation procedure compared to approximated approaches, especially when the intensity function is frequently equal to zero.

To outline the paper, besides a quick introduction to self-regulating Hawkes processes (also referred to as self-correcting Hawkes processes or Hawkes processes with inhibition), Section 2 introduces the concepts of underlying intensity function and restart points. General results concerning the compensator and the exact maximum likelihood estimation procedure are described in Section 3. At last, after a brief discussion about goodness-of-fit in Section 4, Section 5 concludes with a numerical study of the estimation error.

The Hawkes process

Let N be a point process on R * + , where R * + = {x > 0 : x ∈ R}, and (T k ) k≥1 its associated event times (with convention T 0 = 0). For any t ≥ 0, let us note N (t) = k≥1 1 T k ≤t the number of events in [0, t] (where 1 • stands for the indicator function), and λ its conditional intensity function [START_REF] Daley | An introduction to the theory of point processes[END_REF]:

λ(t) = lim h→0 P(N (t + h) -N (t) > 0) h .
A univariate Hawkes process is a point process defined by the conditional intensity function:

λ(t) = λ 0 + t 0 h(t -s) dN (s) + =   λ 0 + T k ≤t h(t -T k )   + , (1)
where x + = max(0, x) denotes the positive part of any real value x, λ 0 ∈ R * + is the baseline intensity and h : R + → R is the kernel, which is assumed to be a monotone measurable function with lim t→∞ h(t) = 0. The kernel function h is the key component of a Hawkes process: it translates the influence (generally assumed to fade away over time) of a past event over the process. Here, h is allowed to take negative values, meaning that it can model both self-exciting and self-regulating Hawkes processes.

Working with such Hawkes processes may prove to be difficult as the positive part function is non-linear. In particular, while computing the compensator function [START_REF] Daley | An introduction to the theory of point processes[END_REF]

Λ(t) = t 0 λ(t) dt, ∀t ≥ 0, (2)
is very easy in the self-exciting case (by linearity of the intensity), it becomes more challenging for the self-regulating Hawkes process. As it is the keystone to derive the likelihood function (and then to obtain a parametric estimation method), our first contribution is to provide an exact expression of the compensator.

For this purpose, let us first introduce the underlying intensity function and the restart time, two quantities which will allow us to derive the computation of the likelihood of a monotone Hawkes process, in a framework unifying self-correcting and self-exciting Hawkes processes. Definition 2.1. Let the underlying intensity function of N be:

λ (t) = λ 0 + t 0 h(t -s) dN (s).
In addition, let the restart time T k be, for any positive integer k:

T k = inf {t ≥ T k | λ(t) > 0}, along with its corresponding cooldown interval τ * k = T * k -T k .
T 1 T 1 T 2 T 2 T 3 T 3 T 4 T 4

Underlying intensity Conditional intensity

Figure 1: Example of the intensity (red curve) and underlying intensity (blue curve) for a self-regulating Hawkes process, with the associated restart times. We only see the negative values of the blue curve since they precisely correspond to the values for which the two intensity functions are not equal.

As illustrated in Figure 1, λ corresponds to the intensity λ as if it were allowed to take negative values. Moreover, as the kernel is assumed to be monotone, the restart time associated to one occurrence can be interpreted as the first moment after this occurrence from which λ and λ become equal (in particular, the restart time and the occurrence time coincide if the intensity function is nonnegative at this time, see Figure 1):

T k = inf {t ≥ T i | ∀t ∈ (T * k , T k+1 ), λ(t) = λ (t)}.

Maximum likelihood estimation and the exponential model

Assume a parametric model P = {λ θ , θ ∈ Θ} for the conditional intensity function λ, where θ contains unknown quantities such as the baseline λ 0 and the kernel h. Then, with convention log(t) = -∞ for t ≤ 0, the log-likelihood t of any θ ∈ Θ with respect to the observations T 1 , . . . , T N (t) in the time interval [0, t] is [19, Proposition 7.2.III.], [START_REF] Ozaki | Maximum likelihood estimation of Hawkes' self-exciting point processes[END_REF]:

t (θ) = N (t) k=1 log (λ θ (T - k )) -Λ θ (t), (3)
where the compensator Λ θ is defined as in Equation ( 2) and 3) reveals the importance of being able to compute the compensator Λ (equivalently Λ θ ) in order to provide a practical implementation of the maximum likelihood estimator of λ. Thus, a first contribution of this paper lies in Proposition 3.1, which establishes a decomposition of the compensator Λ using the underlying intensity function λ and the restart times T 1 , . . . , T N (t) . Proposition 3.1. For any t > 0, the compensator Λ can be expressed as:

λ θ (T - k ) = lim t→T - k λ θ (t). Equation (
Λ(t) =        λ 0 t if t < T 1 λ 0 T 1 + N (t)-1 k=1 T k+1 T k λ (u) du + t T N (t) λ (u) du if t ≥ T 1 , (4)
with the conventions that the sum is equal to 0 if N (t) = 1 and the last integral is equal to

0 if t < T * N (t) .
Proof. This comes directly from splitting the integral of Λ(t) = t 0 λ(t) dt on the intervals [T k , T k+1 ) (k ∈ {0, . . . , N (t) -1}) and [T N (t) , t], and by remarking that, since

h is monotone, ∀t ∈ [T k , T k+1 ), λ(t) = λ (t)1 [T * k ,T k+1 ) (t).
In order to give an explicit computation of the quantity

T k+1 T k λ (u) du (equivalently t T N (t) λ (u) du)
which appears in Proposition 3.1, we focus on the classical scenario where we consider an exponential kernel h(t) = αe -βt , for some α ∈ R and β ∈ R * + . Let us notice that α can be either positive or negative, meaning that the process may be either self-exciting or self-regulating.

Then, the underlying intensity function can be written as:

λ (t) = λ 0 + t 0 αe -β(t-s) dN (s). (5)
The forthcoming proposition steps forward in computing the compensator for an exponential kernel.

Proposition 3.2 (Compensator for exponential kernel). Let t > 0 and k ∈ {1, . . . , N (t)}. The restart times read:

T k = T k + β -1 log λ 0 -λ (T k ) λ 0 1 λ (T k )<0 ,
and the compensator is expressed as in Equation (4), with, for any τ ∈ [T * k , T k+1 ]:

τ T k λ (u) du = λ 0 (τ -T k ) + β -1 (λ (T k ) -λ 0 )(e -β(T k -T k ) -e -β(τ -T k ) ).
Proof. The proof is in A.

Corollary 3.1 (Log-likelihood for exponential kernel). Let

P = λ θ = λ0 + t 0 ᾱe -β(t-s) dN (s) : θ = ( λ0 , ᾱ, β) ∈ Θ , (6)
be a parametric exponential model for the conditional intensity function λ with Θ = R * + × R × R * + , along with the candidate compensator Λ θ , the underlying intensity function λ θ and the restart times T θ,1 , . . . , T θ,N (t) associated to λ θ (see Equation (2) and Definition 2.1).

For any θ = ( λ0 , ᾱ, β) ∈ Θ, by denoting

Λ θ,k = λ0 (T k -T θ,k-1 ) + β-1 (λ θ (T k-1 ) -λ0 )(e -β(T θ,k-1 -T k-1 ) -e -β(T k -T k-1 ) ),
the log-likelihood reads (with convention log(x) = -∞ for x ≤ 0):

t (θ) = log λ0 -λ0 T 1 + N (t) k=2 log λ0 + (λ θ (T k-1 ) -λ0 )e -β(T k -T k-1 ) -Λ θ,k -λ0 (t -T θ,N (t) ) + β-1 (λ θ (T N (t) ) -λ0 ) e -β(T θ,N (t) -T N (t) ) -e -β(t-T N (t) ) 1 t>T θ,N (t) . (7)
Proof. By Equation [START_REF] Linderman | Discovering latent network structure in point process data[END_REF] in the proof of Proposition 3.2,

λ θ (T - k ) = λ0 if k = 1, λ0 + (λ θ (T k-1 ) -λ0 )e -β(T k -T k-1 ) if k ≥ 2.
Combining this expression with Propositions 3.1 and 3.2 leads to the result.

Corollary 3.1 exhibits that the log-likelihood for self-regulating Hawkes processes with an exponential kernel can be evaluated in O(N (t)) operations (by computing iteratively the quantities T θ,k and Λ θ,k appearing in the summation of Equation ( 7)), as already known for self-exciting exponential Hawkes processes [START_REF] Laub | Hawkes processes: Simulation, estimation, and validation[END_REF]Chapter 4.2]. For other monotone kernels without the Markov property, evaluating the log-likelihood with the method proposed here requires O(N (t) 2 ) operations, similarly to existing approaches for self-exciting Hawkes processes.

Goodness-of-fit

Even though computing the compensator Λ (equivalently Λ θ ) was clearly motivated by maximum likelihood estimation, it turns out that it is of great benefit to assess goodness-of-fit, and in particular to check the validity of a maximum likelihood estimation. This is possible thanks to the Time Change Theorem, a result originally stated for inhomogeneous Poisson processes. Let us note that we can find applications of Theorem 4.1 to self-exciting Hawkes processes in the literature [START_REF] Laub | Hawkes processes: Simulation, estimation, and validation[END_REF]Chapter 5]. Since for self-regulating Hawkes processes Λ is still monotone, this result can also be applied in our case.

To be more precise, let us consider θ ∈ Θ and the null hypothesis: "(U k ) k≥1 is a realization of an exponential Hawkes process with parameter θ". This hypothesis can be tested by applying a Kolmogorov-Smirnov test between the empirical distribution of (Λ θ (U k+1 ) -Λ θ (U k )) k≥1 and an exponential distribution with parameter 1. This procedure is illustrated in Table 1, Section 5.

Numerical Results

This section is aimed at assessing the maximum likelihood estimation method for self-regulating Hawkes processes, based on the exact computation of the compensator Λ θ in the exponential model (6) (Corollary 3.1). This procedure is compared to the approximated maximum likelihood estimation proposed in [START_REF] Lemonnier | Nonparametric markovian learning of triggering kernels for mutually exciting and mutually inhibiting multivariate Hawkes processes[END_REF], which consists in approximating Λ θ by:

Λ LM θ (t) = t 0 λ θ (u) du.
This optimization procedure is performed with the L-BFGS-B algorithm from the Scipy package (with (1, 0, 1) as a starting guess and a bounds argument such that λ 0 ≥ 0, α ∈ R, β ≥ 0). In other words, estimators are:

θ ∈ arg max θ∈Θ T Nmax (θ) = Nmax k=1 log (λ θ (T - k )) -Λ θ (T Nmax ) ,
where N max = 200 is the total number of jumps and Λ θ can be replaced by Λ LM θ to obtain the approximated likelihood proposed in [START_REF] Lemonnier | Nonparametric markovian learning of triggering kernels for mutually exciting and mutually inhibiting multivariate Hawkes processes[END_REF].

The comparison between the exact and the approximated estimation procedure is based on simulated data sets coming from self-correcting Hawkes processes of the form (6) with 6 different values of θ = ( λ0 , ᾱ, β) ∈ Θ (see Table 1) which have been chosen in order to explore different scenarios, in particular depending on whether the intensity function is frequently null or not. Observations are sets of time jumps generated with a sampling algorithm (see the algorithm in B and Python implementation online), which is a particular case of Ogata's thinning simulation method [START_REF] Ogata | On Lewis' simulation method for point processes[END_REF] that can handle Hawkes processes with either self-excitation or inhibition. Figure 2 represents the relative absolute errors of estimations θ = ( λ0 , α, β) for each of the 6 simulated models. We observe that the exact approach provides more accurate estimations than the approximated procedure (as illustrated in the boxplots of Figure 2 and by the p-values of the goodnessof-fit tests in Table 1). As expected, the more time the conditional intensity equals 0 (from left to right in Figure 2), the greater the differences between the two procedures. Furthermore, the leftmost boxplot confirms that when the underlying intensity is nonnegative both methods are mostly identical. Let us note that in this case the estimation of ᾱ is rather wrong (the estimation of β is impacted consequently) probably because its value is close to 0 compared to the magnitude of λ0 .

Discussion

In this paper we proposed a maximum likelihood approach for Hawkes processes that can handle both self-exciting and self-regulating scenarios, the first case being already covered in the literature and the latter being our main contribution. For this purpose, we define the concepts of underlying intensity function and restart times when working with monotone kernel functions. In particular we obtain exact expressions of the compensator for the exponential Hawkes process which is the key step of the estimation procedure. We present numerical results on synthetic data that show the efficiency of our procedure, with a substantial improvement compared to approximated approaches when the intensity function is frequently null.

From a theoretical point of view, future work will consist in adapting analytical results to study the convergence of our estimator in the self-regulating case. Regarding modeling, it would be of great interest to consider kernel functions outside the classical exponential scenario. Another important step is the extension of our concepts and algorithms to the multivariate version of the process, which is not straightforward since in the multivariate setting the expression of the restart times are no longer explicit. This last point is essential in order to target real-world datasets since in many applications, being limited to the univariate case will lead to detect self-excitation. However, a model that accounts for potential inhibition effects is of great interest when considering interactions between events of different natures, which will typically be modeled by a multivariate process. This multidimensional extension is the object of a future work, with a further perspective to use our procedure in neuroscience applications in order to detect attraction and repulsion effects between neurons.

Theorem 4 . 1 (

 41 [START_REF] Daley | An introduction to the theory of point processes[END_REF] Theorem 7.4.IV]). Assume that Λ is continuous, monotone and Λ(t) ---→ t→∞ ∞ a.s. Then a.s., a sequence of event times (U k ) k≥1 is a realization of N if and only if (Λ(U k )) k≥1 is a realization of a homogeneous Poisson process with unit intensity.

Figure 2 :

 2 Figure 2: Top panel: relative absolute errors of estimations θ = ( λ0 , α, β). Bottom panel: example of simulated intensities for each set of values θ = ( λ0 , ᾱ, β) with the corresponding average percentage of time when the intensities are equal to zero.

Table 1 :

 1 × 10 7 -1.15 × 10 10 0.37 5.12 × 10 -06 Quantitative assessment of the numerical study: sets of true parameters (left), average estimations over 100 repetitions (middle) and average p-values for the test of Section 4.

			Parameters			Estimations
		λ0	ᾱ	β	λ0	α	β	p-value
	Exact Approx	0.5	-0.001	0.4	0.53 0.54	0.05 0.05	4.25 4.23	0.78 0.78
	Exact Approx	0.5	-0.2	0.4	0.52 0.52	-0.21 -0.22	0.42 0.44	0.72 0.70
	Exact Approx	1.05	-0.75	0.8	1.06 1.14	-0.76 -0.88	0.80 0.82	0.69 0.55
	Exact Approx	2.43	-0.98	0.4	2.55 2.83	-1.01 -1.22	0.39 0.42	0.73 0.51
	Exact Approx	2.85	-2.5	1.8	2.86 8.44 × 10 3 -8.15 × 10 6 2.66 -2.58 1.84	0.73 0.29
	Exact Approx	1.6	-0.75	0.1	1.61 1.36	-0.75	0.11	0.70

https://github.com/migmtz/hawkes-inhibition-expon

Acknowledgments

We thank the Associate Editor and the referees for their valuable comments helping to improve greatly the overall quality of this letter.

A Proof of Proposition 3.2

Let us begin by expressing the underlying intensity function between two event times. First, λ (t) = λ 0 for t ∈ [0, T 1 ). Then, for any k ∈ N, for all t ∈ [T k , T k+1 ), the underlying intensity is differentiable in t and (λ ) (t) = -β(λ (t) -λ 0 ), with the left condition: λ k := λ (T k ). Solving this differential equation leads to

Now, by definition of the restart time

, we obtain T k by solving for t: λ (t) = 0. Thus, by Equation ( 8):

Gathering both situations, we obtain the first part of Proposition 3.2:

which is the second part of Proposition 3.2.

B Simulation algorithm

Algorithm 1 builds upon Ogata's thinning simulation method [21, Proposition 1] in order to handle Hawkes processes with either self-excitation or inhibition.

Algorithm 1: Thinning algorithm for monotone Hawkes process.

Input Parameters λ 0 , h a monotone function, and a stopping criteria (end-time T or maximal number of jumps N max ); Initialization Initialize λ k = λ 0 , t k = 0 and list of times T = ∅; while Stopping criteria not fulfilled do

Set

λmax then Add t to sequence of times T ; end Set t k = t cand ; end return the sequence of jumps T .