Etienne Chevalier
email: etienne.chevalier@univ-evry.fr

Yadh Hafsi
email: yadh.hafsi@universite-paris-saclay.fr

Vathana Ly Vath
email: vathana.lyvath@ensiie.fr

Uncovering Market Disorder and Liquidity Trends Detection

Keywords: Liquidity Risk, Quickest Detection, Change-point Detection, Minimax Optimality, Marked Hawkes Processes, Limit Order Book

The primary objective of this paper is to conceive and develop a new methodology to detect notable changes in liquidity within an order-driven market. We study a market liquidity model which allows us to dynamically quantify the level of liquidity of a traded asset using its limit order book data. The proposed metric holds potential for enhancing the aggressiveness of optimal execution algorithms, minimizing market impact and transaction costs, and serving as a reliable indicator of market liquidity for market makers. As part of our approach, we employ Marked Hawkes processes to model trades-through which constitute our liquidity proxy. Subsequently, our focus lies in accurately identifying the moment when a significant increase or decrease in its intensity takes place. We consider the minimax quickest detection problem of unobservable changes in the intensity of a doubly-stochastic Poisson process. The goal is to develop a stopping rule that minimizes the robust Lorden criterion, measured in terms of the number of events until detection, for both worst-case delay and false alarm constraint. We prove our procedure's optimality in the case of a Cox process with simultaneous jumps, while considering a finite time horizon. Finally, this novel approach is empirically validated by means of real market data analyses.

Introduction

Assets liquidity is an important factor in ensuring the efficient functioning of a market. [START_REF] Glosten | Estimating the components of the bid/ask spread[END_REF] define liquidity as the ability of an asset to be traded rapidly, in significant volumes and with minimal price impact. Measuring liquidity, therefore, involves three aspects of the trading process: time, volume and price. This is reflected in Kyle's description of liquidity as a measure of the tightness of the bid-ask spread, the depth of the limit order book and its resilience (see [START_REF] Kyle | Continuous Auctions and Insider Trading[END_REF]). Hence, as highlighted in [START_REF] Lybek | Measuring Liquidity in Financial Markets[END_REF] and in [START_REF] Bińkowski | Endogeneous Dynamics of Intraday Liquidity[END_REF], it is essential for liquidity-driven variables to not only capture the transaction cost associated with a relatively small quantity of shares but also to assess the depth accessible to large market participants.

1

Additionally, measures related to the pace or speed of the market are considered since trading slowly can help mitigate implicit transaction costs. Typical indicators in this category include value traded and volatility over a reference time interval.

Several studies have been devoted to identifying different market regimes rapidly and in an automated fashion. [START_REF] Hamilton | Rational-expectations econometric analysis of changes in regime: An investigation of the term structure of interest rates[END_REF] proposes the initial notion of regime switches, wherein he establishes a relationship between cycles of economic activity and business cycle regimes. The idea remains relevant, with ongoing exploration of novel approaches. [START_REF] Horvath | Clustering Market Regimes using the Wasserstein Distance[END_REF] suggest an unsupervised learning algorithm for partitioning financial time series into different market regimes based on the Wassertein k-means algorithm for example. [START_REF] Bucci | Market regime detection via realized covariances[END_REF] employ covariance matrices to identify market regimes and identify shifts in volatile states.

This paper introduces a novel liquidity regime change detection methodology aimed at assessing the resilience of an order book using tick-by-tick market data. More precisely, our objective is to study methods derived from the theory of disorder detection and to develop a liquidity proxy that effectively captures its inherent characteristics. This will enable us to thoroughly examine how the distribution of the liquidity proxy evolves over time. By employing these methods, we seek to gain a detailed understanding of the dynamics involved in liquidity changes and their impact within the distribution patterns of our chosen proxy. To accomplish this, our methodology centers around the introduction of a liquidity proxy known as "trades-through" (see Definition 2.1). Trades-through possess a high informational content, which is of interest in this study. Specifically, trades-through serve as an indicator of the order book's resilience, allowing the examination of activity at various levels of depth. They can also provide insightful information on the volatility of the financial instrument under examination. Additionally, trades-through can indicate whether an order book has been depleted or not in comparison to its previous states and to which extent. [START_REF] Pomponio | Multiple-limit trades: Empirical facts and application to lead-lag measures[END_REF] investigate several stylized facts related to trades-through, which are similar to other microstructural features observed in the literature. Their research highlights the statistical robustness of tradesthrough concerning the size of orders placed. This implies that trades-though are not solely a result of low quantities available on the best limits and that the information they provide is of significant importance. Another noteworthy result they observe is the presence of self-excitement and clustering patterns in trades-through. In other words, a trade-through is more likely to occur soon after another trade-through than after any other trade. This result supports the idea of modelling trades-through using Hawkes processes as in [START_REF] Toke | Modelling trades-through in a limit order book using Hawkes processes[END_REF]. This approach is particularly interesting given that Hawkes processes fall into the class of branched processes. Their dynamics can therefore be illustrated by a representation of immigration and birth that allows interpretation. Indeed, the intensity, at which events take place, is formed by an exogenous term representing the arrival of "immigrants" and a self-reflexive endogenous term representing the fertility of "immigrants" and their "descendants". We may refer to [START_REF] Brémaud | HAWKES BRANCHING POINT PROCESSES WITHOUT ANCESTORS[END_REF] for more comprehensive details.

Our focus here will be to use the Marked-Hawkes process to model trades-through rather than standard Hawkes processes in order to capture the impact of the orders' volumes on the intensity of the counting processes. Our research centers on employing the Marked-Hawkes process as a modelling technique for trades-through, as opposed to conventional Hawkes processes. The rationale behind this choice lies in our intention to account for the influence of order volumes on the intensity of the counting processes.

Upon the conclusion of the modelling phase, we will utilize this proxy to identify intraday liquidity regimes 1 . In order to achieve this, our model involves detecting the times at which the distribution of liquidity (represented by a counting process) undergoes changes as fast as possible. This entails comparing the intensities of two doubly stochastic Poisson processes. This type of problem is commonly known in literature as "Quickest change-point detection" or "Disorder detection". Disorder detection has been studied through two research paths, namely the Bayesian approach and the minimax approach. The Bayesian perspective, originally proposed by [START_REF] Girshick | A Bayes Approach to a Quality Control Model[END_REF], typically assumes that the change point is a random variable and provides prior knowledge about its distribution. As an example, [START_REF] Dayanik | Compound Poisson Disorder Problem[END_REF] conduct their study on compound Poisson processes assuming an exponential prior distribution for the disorder time. In their work, [START_REF] Bayraktar | Poisson Disorder Problem with Exponential Penalty for Delay[END_REF] address a Poisson disorder problem with an exponential penalty for delays, assuming an ex-ponential distribution for the delay. However, the minimax approach does not include any prior on the distribution of the disorder times, see for instance, [START_REF] Page | Continuous inspection schemes[END_REF] and [START_REF] Liniger | Procedures for Reacting to a Change in Distribution[END_REF]. Due to the scarcity of relevant literature and data pertaining to disorder detection problems, accurately estimating the prior distribution in the context of market microstructure poses a challenging task. Hence, our preference lies in adopting a non-Bayesian min-max approach. [START_REF] Karoui | Minimax optimality in robust detection of a disorder time in doubly-stochastic poisson processes[END_REF] prove the optimality of the CUSUM 2 procedure 3.1 for ρ < 1 and ρ > 1 for doubly stochastic poisson processes which extends the proof proposed by G. [START_REF] Moustakides | Sequential change detection revisited[END_REF] and [START_REF] Poor | Quickest Detection[END_REF] for ρ < 1. However, the proposed solution is limited to cases where point processes do not exhibit simultaneous jumps. This presents a problem in our case since definition 2.1 suggests that each time a n-limit trade-through event is observed, n -1 other trade-through events occur simultaneously. Hence, we need to demonstrate the optimality of these results by considering scenarios in which point processes exhibit a finite number of simultaneous arrival times which is precisely what we obtained. We are also able to produce a tractable formula of the average run delay. More importantly, our findings encompass the intricacies inherent to modelling limit order books and investigating market microstructure. We use our results to identify changes in liquidity regimes within an order book. To do this, we apply our detection procedure to real market data, enabling us to achieve convincing detection performance. To the best of our knowledge, we believe that this is the first attempt to quantify market liquidity using a methodology combining the notion of trades-through, disorder detection theory, and marked-hawkes processes.

The article is structured in the following manner. In Section 2, the primary aim is to establish a model for trades-through using Marked Hawkes Processes (MHP) and subsequently provide the relevant mathematical framework. In Section 3, various results related to the optimality of the CUSUM procedure within the framework of sequential test analysis for a simultaneous jump Cox process will be presented. Finally, in Section 4, we provide an overview of the goodness-of-fit results obtained from applying a Marked Hawkes Process to analyze trades-through data in the Limit Order Book of BNP Paribas stock. Additionally, we present the outcomes of utilizing our disorder detection methodology to identify various liquidity regimes.

Trades-Through modelling and Hawkes Processes

The aim of this section is to build a model for trades-through by means of Marked Hawkes Processes (MHP) and to present the corresponding mathematical framework. We begin by defining the concept of trades-through.

Definition 2. . Formally, a trade-through can be defined by the vector (type, depth, volume) ∈ {-1, 1} × N × R + where type equals 1 if the trade-through occurs on the bid side of the book and -1 otherwise, depth represents the number of limits consumed by the market event relative to the trade-through, and volume represents the corresponding traded volume.

In brief, an n-limit trade-through is defined as an event that exhausts the available n-th limit in the order book. It is considered that a trade-through at limit m is also a trade-through of limit n if m is greater than n [START_REF] Pomponio | Multiple-limit trades: Empirical facts and application to lead-lag measures[END_REF]. In this context, the term 'limit n' refers to any order that is positioned n ticks away from either the best bid price or the best ask price. This definition is broad as it encompasses different types of orders, including market orders, cancellations, and hidden orders which underlines its significant information content. Let (Ω, F, F = {F t } t≥0 , P) be a complete filtered probability space endowed with right continuous filtration {F t } t≥0 . Going forward, whenever we mention equalities between random variables, we mean that they hold true almost surely under the probability measure P. We may not explicitly indicate this notion (P-a.s) in most cases. Let 0

< τ A,i 0 < τ A,i 1 < τ A,i 2 < . . . (resp. 0 < τ B,i 0 < τ B,i 1 < τ B,i 2 < . . . )
form a sequence of increasing F-measurable stopping times and (v A,i k ) k≥0 (resp. (v B,i k ) k≥0 ) be a sequence of i.i.d R + -valued random variables for each i ∈ {1, . . . , M }. Here, (τ A,i k ) k≥0 (resp. (τ B,i k ) k≥0 ) represent the arrival times of trade-throughs at limit i on the ask (resp. bid) side of the limit order book while (v A,i k ) k≥0 (resp. (v B,i k ) k≥0 ) represent the sizes of volume jumps of trades-though for these arrival times.

Remark 2.1. The (τ A,i k ) k≥0 (resp. (τ B,i k ) k≥0 ) values are arranged in a specific order for a given i, but it is essential to note that this order may not be consistently maintained when comparing them to the (τ A,j k ) k≥0 (resp. (τ B,j k ) k≥0 ) values for a different j. We proceed by gathering the arrival times of trade-through orders for all the limits on the bid (resp. ask) side. Subsequently, these arrival times are organized into sequences (τ A k ) k≥0 and (τ B k ) k≥0 respectively that are sorted in ascending order.

Our goal is to design a model that takes into account the volumes of the trades-through for all pertinent limits and that replicates the clustering phenomenon observed in the corresponding stylized facts (See [START_REF] Pomponio | Multiple-limit trades: Empirical facts and application to lead-lag measures[END_REF]). To achieve this, we suggest using a modified version of Cox processes3 , specifically Marked Hawkes Processes. We define a 2-dimensional marked point process N A×B = N A , N B where N A (resp. N B ) is the counting measure on the measurable space (R

+ × R + , B ⊗ L) associated to the set of points (τ A k , v A k ); k ≥ 0 (resp. (τ B k , v B k ); k ≥ 0 ), B being the Borel σ-field on R + while L is the Borel σ-field on R + . This means that ∀C ∈ B ⊗ L : N A (C) = k≥0 1 C (τ A k , v A k ) , N B (C) = k≥0 1 C (τ B k , v B k ) = 1≤i≤M k≥0 1 C (τ A,i k , v A,i k ) = 1≤i≤M k≥0 1 C (τ B,i k , v B,i k )
Let F A×B = F A×B t t≥0 be a sub-σ-field of the complete σ-field F where

F A×B t = σ N A×B (s, v); s ∈ [0, t], v ∈ R + .
Simply put, F A×B t contains all the information about the process N A×B up to time t. The collection of all such sub-σ-fields, denoted by F A×B t , t ≥ 0 , is called the internal history of the process N A×B . One of the core concepts underlying the dynamics of N A×B is the notion of conditional intensity. Specifically, the dynamics of N A×B are characterized by the non-negative, F A×B -progressively measurable processes λ A and λ B , which model the conditional intensities. These processes satisfy the following conditions:

E K N A (t, v) -N A (s, v)dv F A×B s = E [s,t[×K λ A (u, v)dudv F A×B s , 0 ≤ s ≤ t E K N B (t, v) -N B (s, v)dv F A×B s = E [s,t[×K λ B (u, v)dudv F A×B s , 0 ≤ s ≤ t (1)
Moving forward, the point process N A×B will be represented as a bivariate Hawkes process, with its intensity wherein its intensity is dependent on the marks (v A k ) k≥0 and (v B k ) k≥0 . By viewing the marked point processes N A and N B as a processes on the product space R + × R + , we can derive the marginal processes N A g (.) = N A (., R + ) and N B g (.) = N B (., R + ), which refer to as ground process. In other worlds, for t ∈ R + ,

N A g (t) = [0,t[×R+ N A ( du × dv A ), N B g (t) = [0,t[×R+ N B ( du × dv B ) (2)
These ground processes therefore describe solely the arrival times of the trades-through and can be understood as being Hawkes processes as well. In instances where N A×B exhibit sufficient regularity (See The definition 6.4.III and section 7.3 in [START_REF] Daley | An Introduction to the Theory of Point Processes[END_REF]), it is possible to factorize its conditional intensities λ A and λ B without presupposing stationarity and under the assumption that the mark's conditional distributions at time t given

F A×B t - have densities of R + → R + : v → f A (v F A×B t - ) and R + → R + : v → f B (v F A×B t -
). The conditional intensity factors out as :

λ A (t, v) = λ A g (t)f A (v F A×B t - ), λ B (t, v) = λ B g (t)f B (v F A×B t - ), (t, v) ∈ R + × R + (3)
where λ g (t) is referred to as the F A×B t - -intensity of the ground intensity. In a heuristic manner, one can summarize equation 3 for (t, v) ∈ R + × R + as being in the form of :

λ A (t, v)dtdv ≈ E N A (dt × dv) F A×B t - ≈ λ A g (t)f A (v F A×B t - )dtdv λ B (t, v)dtdv ≈ E N B (dt × dv) F A×B t - ≈ λ B g (t)f B (v F A×B t - )dtdv
Remark 2.2. It is noteworthy that the conditional intensity of the ground processes considered in this context is evaluated with respect to the filtration

F A×B = (F A×B t ) t≥0 rather than F Ng = (F Ng t ) t>≥0 with F Ng t = σ (N g (s), s ∈ [0, t]
). This is attributed to the fact that unlike F Ng , the filtration F A×B accounts for the information pertaining to the marks associated with the process N A×B .

We assume that the conditional intensities λ A g and λ B g have the following form :

λ i g (t) = µ i (t) + j=A,B k≥0 γ ij (t -τ j k , v j k ) = µ i (t) + j=A,B [0,t[×R+ γ ij (t -s, v)N j ( ds × dv), i = A, B (4)
where µ i : R + → R + and γ ij : R + × R + → R + are non-negative measurable functions. We choose a standard factorized form for the ground intensity kernel

γ ij (t -u, v) = α ij e βij (t-u) × g j (v)
to describe the weights of the marks, g j being a measurable function defined on R + that characterizes the mark's impact on the intensity. The exponential shape of the kernel is derived from the findings of [START_REF] Toke | Modelling trades-through in a limit order book using Hawkes processes[END_REF]. The properties of functions g A and g B will be examined in Section 4.

Remark 2.3. As we progress to proposition 4.2, it will become apparent that the factorization of the latter kernel presents a way to represent the conditional intensities of these processes in a Markovian manner.

Thus, the conditional ground intensities of this model have the following integral representation :

λ A g (t) = µ A (t)+ [0,t[×R+ α AA e -β AA (t-u) g A (v) N A (du × dv) + [0,t[×R+ α AB e -β AB (t-u) g B (v) N B (du × dv), t ≥ 0 (5) λ B g (t) = µ B (t)+ [0,t[×R+ α BA e -β BA (t-u) g A (v) N A (du × dv) + [0,t[×R+ α BB e -β BB (t-u) g B (v) N B (du × dv), t ≥ 0 (6)
where (α ij , β ij ) (i,j)∈{A,B} 2 are non-negative reals.

Conventionally, the impact functions g i have to satisfy the following normalizing condition which enhances the overall stability of the numerical results (see section 1.3 of Liniger (2009) for further details) :

+∞ 0 g i (v)f i (v)dv = 1, i = A, B. (7)
This means that the impact functions of the Marked-Hawkes based model would be equal to

g i (v) = v η i E(v η i ) .
The following stability conditions of our model are based on Theroem 8 of [START_REF] Brémaud | HAWKES BRANCHING POINT PROCESSES WITHOUT ANCESTORS[END_REF].

Proposition 2.1 (Existence and uniqueness). Suppose that the following conditions hold :

1. The spectral radius of the branching matrix (also referred to as branching ratio) satisfies :

sup λ∈ρ(∥Γ∥)
|λ| < 1 2. The decay functions satisfy

+∞ 0 tγ ij (t)dt < +∞, ∀i, j ∈ {A, B}
where ∥Γ∥ = γ ij (i,j)∈{A,B} 2 , γ ij (t) = α ij e -βij t g j (v(t)) and ρ (∥Γ∥) represents the set of all eigenvalues of ∥Γ(t)∥, then there exists a unique point process N A×B with associated intensity process λ A×B .

Remark 2.4. Existence in proposition 2.1 means that we can find a probability space (Ω, F, P) which is rich enough to support such the process N A×B . Uniqueness means that any two processes complying with the above conditions have the same distribution. In our case, the conditions mentioned above are described as follows :

1. 1 2 α AA β AA + α BB β BB + α AA β AA -α BB β BB 2 + 4 α AB β AB α BA β BA < 1 2. α ij β ij 2 < +∞, ∀i, j ∈ {A, B}
3 Sequential Change-Point Detection: CUSUM-based optimal stopping scheme

We now shift our focus towards developing a methodology that can distinguish between distinct liquidity regimes. This can be done by detecting changes or disruptions in the distribution of a given liquidity proxy, which, in our case, is represented by the number of trades-through. Here, we will focus on identifying disruptions that affect the intensity of the marked multivariate Hawkes process discussed in the previous section. This will enable us to compare the resilience of our order book to that of previous days over a finite horizon T . The sequential detection methodology involves comparing the distribution of the observations to a predefined target distribution. The objective is to detect changes in the state of the observations as rapidly as possible while minimizing the number of false alarms which correspond to sudden changes in the arrival rate in the case of Poisson processes with stochastic intensity. To be specific, let us consider a general point process N with a known rate λ, where the arrivals of certain events are associated to this process. At a certain point in time θ, the rate of occurrence of events of the process N undergoes an abrupt shift from λ to ρλ. However, the value of the disorder time θ is unobservable. The goal is to identify a stopping time τ that relies exclusively on past and current observations of the point process N and can detect the occurrence of the disorder time θ swiftly. Thus, we will focus in the following on solving the sequential hypothesis testing problems in the aforementioned general form :

H 0 : λ(t), H 1 : λθ1 (t) = λ(t)1 {t<θ1} + ρ 1 λ(t)1 {t≥θ1} , ρ 1 > 1 (8)
and,

H 0 : λ(t), H 2 : λθ2 (t) = λ(t)1 {t<θ2} + ρ 2 λ(t)1 {t≥θ2} , ρ 2 < 1 (9)
We propose the introduce a change-point detection procedure that is applicable to our model. The goal is to study the distribution of the process N A×B process. The problem at hand can be tackled through various approaches. Here, we opt for a non-Bayesian approach, specifically the min-max approach introduced by PAGE (1954), which assumes that the timing of the change is unknown and non-random. Previous research has shown the optimality of the CUSUM procedure in solving this problem (see [START_REF] Poor | Quickest Detection[END_REF]G. V. [START_REF] Moustakides | Optimality of the CUSUM procedure in continuous time[END_REF] [START_REF] Karoui | Minimax optimality in robust detection of a disorder time in doubly-stochastic poisson processes[END_REF]). However, these demonstrations are confined to scenarios where the process involves Cox jumps of size 1. This constraint poses significant limitations in our case, as the definition 2.1 of trade-throughs implies for of simultaneous jumps. Therefore, the objective of this chapter is to demonstrate the optimality of the CUSUM procedure in a broader context that considers multiple jumps and to formulate the relevant detection delay.

The identification of different liquidity regimes in an order book can be done in various ways. One can apply our disorder detection methodology separately on processes N A and N B to study liquidity either on the bid or on the ask, or consider process N A + N B and take both into account at the same time. To enhance the clarity and maintain the general applicability of our results, we consider the point process N as a general Cox process that decomposes into the form D i=1 N i where 0 < D < +∞ and N i is a finite point process defined on the same probability space introduced in the previous section and whose arrivals are non-overlapping. We denote λ i the F-intensity of N i and by t → Λ i (t) = t 0 λ i (s)ds it's compensator for i = 1, . . . , D.

We denote the probability measure for the case with no disorder by P ∞ , where θ = +∞. The probability measure for the case with a disorder that starts from the beginning, where θ = 0, is denoted by P 0 . If there is a disorder at the instant θ, represented by the value of the intensity of the counting process equal to λθ (t), then the probability measure is denoted by P θ . The constructed measure satisfies the following conditions :

P θ = P 0 , if θ = 0 P ∞ , if θ = +∞ Remark 3.1.
In what follows, we use the notation E to refer to the expectations that are evaluated under the probability measure P ∞ , while E θ represents the expectations evaluated under the probability measure P θ .

We initiate the discourse by introducing a few notations that will be instrumental in presenting our results. Let U i = N i -β(ρ)Λ i and i ∈ {1, . . . , D}. The process D i=1 U i is referred to as the Log Sequential Probability Ratio (LSPR) between the two probability measure P ∞ and P 0 where β(ρ) = (ρ -1)/ log ρ and Λ i (t) = t 0 λ i (s)ds for all 0 ≤ t ≤ T . The probability measure P 0 is defined as the measure equivalent to P ∞ with density :

dP 0 dP ∞ Ft = ρ D i=1 U i (t) , 0 ≤ t ≤ T (10)
Or, more generally :

dP θ dP ∞ Ft = D i=1 ρ U i (t) ρ U i (θ) , 0 ≤ θ ≤ t = ρ D i=1 U i (t)-U i (θ) (11)
According to Girsanov's theorem, these densities are defined if [START_REF] Øksendal | Applied Stochastic Control of Jump Diffusions[END_REF]). Hence,

E ρ D i=1 U i (t)-U i (θ) = 1 which is the case since ρ D i=1 Ui(t) = e log(ρ) D i=1 N i (t)-(ρ-1) D i=1 Λ i (t) is a local martingale if D i=1 Λ i (t) is càdlàg (see
D i=1 N i (t) -ρΛ i (t) is a P θ -martingale iff ρ D i=1 U i (t)-U i (θ)
is an P ∞ -martingale. We can therefore now define the log-likelihood ratio :

ℓ t,θ : = D i=1 t θ log λi θ (s) λ i (s) dN i (s) - D i=1 t θ λi θ (s) -λ i (s) ds = log (ρ) D i=1 N i (t) -N i (θ) - ρ -1 log (ρ) Λ i (t) -Λ i (θ) = log (ρ) D i=1 U i (t) -U i (θ) (12) with λi θ (t) = λ i (t), 0 ≤ t ≤ θ ρλ i (t), t > θ
which is the intensity of the node i if the disorder takes place at the time θ.

In statistical problems, it is common to observe two distinct components of loss. The first component is typically associated with the costs of conducting the experiment, or any delay in reaching a final decision that may incur additional losses. The second component is related to the accuracy of the analysis. The min-max problem (see [START_REF] Page | Continuous inspection schemes[END_REF]) can be formulated as finding a stopping rule that minimizes the worst-case expected cost under the Lorden criterion. The objective is therefore to find the stopping time that minimizes the following cost :

inf τ ∈Tπ C(τ ) ( 13
)
Where

C(τ ) = sup θ∈[0,+∞] ess sup E θ (τ -θ) + F θ , T π is the class of F-stopping times that statisfies E(τ ) ≥ π and π > 0 a constant.
As explained by [START_REF] Liniger | Procedures for Reacting to a Change in Distribution[END_REF], the formulation of the problem under criterion C(τ ) is robust in the sense that it seeks to estimate the worst possible delay before the disorder time. The false alarm rate can be controlled by controlling E(τ ) through π. It has been observed that Lorden's optimality criterion is characterized by a high level of stringency as it necessitates the optimization of the maximum average detection delay over all feasible observation paths leading up to and following the changepoint. This rigorous approach often produces conservative outcomes. However, it is worth noting that the motivation behind the use of min-max optimization can be linked to the Kullback-Leibler divergence. This has notably been the case in the works of G. V. [START_REF] Moustakides | Optimality of the CUSUM procedure in continuous time[END_REF] for detecting disorder occurring in the drift of Ito processes. In our case, for τ ∈ T π and 0 ≤ θ ≤ τ , equation 12 yields :

E θ log dP θ dP ∞ F θ = log (ρ) E θ D i=1 N i (τ ) -N i (θ) - ρ -1 log (ρ) Λ i (τ ) -Λ i (θ) F θ = (log (ρ) -ρ + 1) D i=1 E θ N i (τ ) -N i (θ) + F θ (14)
Hence, the following minimax formulation of the detection problem is proposed :

inf τ C(τ ) = inf τ sup θ∈[0,+∞] ess sup D i=1 E θ N i (τ ) -N i (θ) + F θ with the constraint D i=1 E N i (τ ) ≥ π (15)
Karoui, Loisel, and Salhi (2017) suggests another justification of this relaxation based on the Timerescaling theorem (see [START_REF] Daley | An Introduction to the Theory of Point Processes[END_REF]) which allows us to transform a Cox process into a Poisson process verifying :

E θ D i=1 N i (τ ) -N i (θ) + F θ = ρ D i=1 E θ τ τ ∧θ λ i (s)ds F θ (16)
Heuristically, this is equivalent to the formulation referenced in 15. Indeed,

D i=1 E θ τ τ ∧θ λ i (s)ds F θ = D i=1 λ i,cst E θ (τ -θ)
+ F θ in the case of a Poisson process with constant intensity.

Remark 3.2. The formulation 15 is particularly useful in finance trading applications, where it is more advantageous to focus solely on transactions/events that occurred, rather than the state of the market at every moment.

Remark 3.3. Another alternative is to consider the criterion

E θ sup 1≤i≤D N i (τ ) -N i (θ) + F θ ,
which would be to minimize the worst possible detection delay among all processes (N i (t)) 1≤i≤D . Since the Lorden criterion is already a conservative measure, adding another criteria function could result in excessive restrictions that may cause significant delays in the detection process. Therefore, it is preferable to refrain from using additional criteria functions. Notably, the Lorden criterion produces larger values compared to the criteria proposed by [START_REF] Pollak | Optimality and Almost Optimality of Mixture Stopping Rules[END_REF] or [START_REF] Shiryaev | On stochastic models and optimal methods in the quickest detection problems[END_REF], which further supports its sufficiency in detecting changes (see Albert [START_REF] Shiryaev | Stochastic Disorder Problems[END_REF] p. 9-15).

Given that the likelihood ratio test is considered the most powerful test for comparing simple hypotheses according to the Neyman-Pearson lemma, we will utilize ρ D i=1 U i (t)-U i (θ) to assess the shift in this distribution. The log-likelihood ratio D i=1 U i (t) -U i (θ) usually exhibits a negative trend before a change and a positive trend after a change (see section 8.2.1 of Tartakovsky, Nikiforov, and Basseville (2014) for more details). Consequently, the key factor in detecting a change is the difference between the current value of D i=1 U i (t) -U i (θ) and its minimum value. As a result, we are prompted to introduce the CUSUM process.

Definition 3.1. Let ρ ̸ = 1, m > 0, 0 ≤ t ≤ T and U = D i=1 N i -βΛ i = D i=1 U i .
The CUSUM processes are defined as the reflected processes of U :

Û (t) = U (t) -inf 0≤s≤t U (s), if ρ > 1 U (t) = sup 0≤s≤t U (s) -U (t), if ρ < 1 (17)
The stopping times for the CUSUM procedures are given by :

TC = inf t : Û (t) > m , T C = inf t : U (t) > m
where the infimum of the empty set is +∞.

Remark 3.4. We note that ρβ(1/ρ) = β(ρ). Therefore, for t > 0,

sup θ<t ℓ t,θ ≤ ℓ C t,θ :=    log (ρ) D i=1 U i (t) -inf θ<t U i (θ) if ρ > 1 log 1 ρ D i=1 sup θ<t U i (θ) -U i (t) if ρ < 1 and t : sup θ<t ℓ t,θ > m ⊆ t : ℓ C t,θ > m
Consequently, initial lower bounds can be derived for the stopping times TC and T C based on the reflected CUSUM processes associated with each component N i .

Remark 3.5. Considering the increasing event times {τ i , i = 1, 2, . . .} of the process D i=1 N i , and for any fixed t ≥ 0 and k ≥ 1 where t > t k , the following relationships hold (see Lemma 4.1 in Wang et al. ( 2021)) :

sup τ k <s≤min{τ k+1 ,t} ℓ t,s = lim s→τ + k ℓ t,s = ℓ t,τ + k
, and sup 0≤s≤τ1 ℓ t,s = ℓ t,0 . This implies for k > 0 and t > τ k that :

sup θ<t ℓ t,θ = max{ sup 0<s≤τ1 ℓ t,s , sup τ1<s≤τ2 ℓ t,s , . . . , sup τ k <s≤min{τ k+1 ,t} ℓ t,s } = max{ℓ t,0 , ℓ t,τ + 1 , . . . , ℓ t,τ + k }
This allows us to extract a more time-efficient method for calculating the optimal stopping times T C and TC . To do this, we formulate the log-likelihood ratio recursively :

ℓ t,τ + n+1 = D i=1 U i (t) -U i (τ + n+1 ) = D i=1 U i (t) -U i (τ + n+1 ) + U i (τ + n ) -U i (τ + n ) = ℓ t,τ + n - D i=1 N i (τ + n+1 ) -N i (τ + n ) + β (ρ) D i=1 τ + n+1 τ + n λ i (s)ds = ℓ t,τ + n - D i=1 N i (τ + n+1 ) -N i (τ + n ) + β (ρ) D i=1 Λ i (τ + n+1 ) -Λ i (τ + n ) One noteworthy characteristic of D i=1 N i ( T C ) and D i=1 N i ( TC )
is that their conditional expectations can be explicitly calculated based on their initial values and the parameter m. As a result, the theorems 3.1 and 3.2 offer a valuable tool for controlling the ARL constraint through the parameter m, i.e, expressing

D i=1 E N i ( T C ) and D i=1 E N i ( TC ) as tractable formulas for ρ > 1 and ρ < 1. Theorem 3.1. Let T C the CUSUM stopping time defined by T C = inf t ≤ T : U (t) > m for ρ < 1.
Assume that the intensities λ i of the processes N i are càdlàg, ∀i ∈ {1, . . . , D}. The constraint on the Average Run Length (ARL) in T C is equal to :

E y D i=1 N i ( T C ) = m y 1 β(ρ) ⌊x⌋ k=0 (-1) k k! ((x -k)/β(ρ)) k exp((x -k)/β(ρ))dx (18)
with P y = P(. U (0) = y).

Proof. Proof postponed to appendix.

Theorem 3.2. Let TC the CUSUM stopping time defined by TC = inf t ≤ T : Û (t) > m for ρ > 1.

Assume that the intensities λ i of the processes N i are càdlàg, ∀i ∈ {1, . . . , D}. The constraint on the Average Run Length (ARL) in TC is equal to :

E v D i=1 N i ( TC ) = W (m -v) W (m) W ′ (m) - m-v 0 W (y)dy, E m - D i=1 N i ( TC ) = W (m) βW ′ (m) (19)
with

P v = P(. Û (0) = v), W (x) = 1 β(ρ) ⌊x⌋ k=0 (-1) k k! ((x -k)/β(ρ)) k exp((x -k)/β(ρ)) and x 0 W (y)dy = ⌊x⌋ k=0 e (x-k)/β(ρ) k i=0 (-1) j j! ((x -k)/β(ρ)) j -1 .
Proof. Proof postponed to appendix. Two subsequent theorems 3.3 and 3.4 substantiate the optimality of the CUSUM procedure by separately considering the cases where ρ < 1 and ρ > 1. These theorems establish that the CUSUM stopping time attains the lower bound of the Lorden criterion C(τ ) = sup θ∈[0,+∞[ ess sup E θ (τ -θ) + F θ in each case, thereby demonstrating its optimality in detecting changes or shifts in sequential observations for both scenarios.

Theorem 3.3. Let T C = inf t ≤ T : U (t) > m denote the CUSUM stopping time. Then, T C is
proven to be the optimal solution for 15 for ρ < 1 where

D i=1 E N i (T ) = D i=1 E N i ( T C ) .
Proof. Proof postponed to appendix.

Theorem 3.4. Let TC = inf t ≤ T : Û (t) > m denote the CUSUM stopping time. Then, TC is proven to be the optimal solution for 15 for ρ > 1 where

D i=1 E N i (T ) = D i=1 E N i ( TC ) .
Proof. Proof postponed to appendix.

Experimental Results

We proceed to presenting the experimental findings related to our methodology. To begin with, we provide a comprehensive outline of the final form of the kernel in the Marked Hawkes Process. Subsequently, we conduct an assessment of its goodness of fit to the BNP Paribas data. Following that, we perform a thorough analysis to evaluate the methodology's effectiveness in detecting changes in liquidity regimes.

Data Description

We use Refinitiv limit order book tick-by-tick data (timestamp to the nanosecond) for BNP-Paribas (BNP.PA) stock from 01/01/2022 to 31/12/2022 to fit our Marked Hawkes models and carry out our research on disorder detection. This data give us access to the full depth of the order book and include prices and volumes on the bid and ask. Throughout the study, we only consider trades flagged as 'normal' trades. In particular, we did not consider any block-trade or off-book trade in the following. We restrict our data time-frame to the period of the day not impacted by auction phases (from 9:30 PM to 17 PM, Paris time). Only trades ranging from limit one to four are considered here. The selection is motivated by the relatively rare instances of trade-through events when the limit exceeds four. Indeed, the comprehensive analysis of transactions involving the BNP stock reveals that approximately 4.89% of the overall transaction count led to at least one trade-through event. While these trade-through events bear significance, it is pertinent to note that trades-through with limit surpassing four constitute a low percentage of merely 0.0126% within the total transaction count during that specific timeframe. These findings align with the research conducted by [START_REF] Pomponio | Multiple-limit trades: Empirical facts and application to lead-lag measures[END_REF]. The VSTOXX® Volatility Index serves as a tool for assessing the price dynamics and liquidity of stocks comprising the EURO STOXX 50® Index throughout the year 2022. This index will facilitate the identification of appropriate days for conducting our comparative analysis. Figure 3 shows that day 05/01/2022 corresponds to the day with the lowest volatility, while 04/03/2022 is the day with the highest volatility. We will consider 31/08/2022 as the reference day where the volatility index is set equal to the average of the VSTOXX® (index vol ≈ 27) over the entire year and where the volatility index on the previous and following days does not deviate greatly from the latter.

Estimation

The current state of the art regarding the estimation of Multivariate Hawkes processes revolves around methods : Method of moments, Maximum likelihood estimation (MLE) and Least squares estimation.

A number of papers have reviewed each method. For exemple, Cartea, Cohen, and Labyad (2021) construct an adaptive stratified sampling parametric estimator of the gradient of the least squares estimator. [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF] used the techniques developed by Bartlett to analyze the spectra of point processes that where which were later retrieved by Bacry, Dayri, and J. Muzy (2011) to propose a nonparametric estimation method for stationary Multivariate Hawkes processes with symmetric kernels. Bacry and J.-F. Muzy (2016) relaxed these assumptions, except for stationarity, and showed that the Multivariate Hawkes process parameters solve a system of Wiener-Hopf equations. In their work, Daley and Vere-Jones (2007) put forth the idea of maximizing the log-likelihood of the sample path.

We adopt this approach due to the efficient computation of the likelihood function in our case, which consequently eases the estimation of the model's parameters.

The subsequent analysis involves parameter estimation within our model, utilizing the maximum likelihood method. We define an observation period denoted as [0, T ], which corresponds to the time span during which empirical data was collected. To construct the likelihood function 20, we initiate by defining the probability densities f A and f B of the marks. Figure 4 The observations exhibit a good fit with an exponential distribution. We will characterize the impact of volume on trades-through using a normalized power decay function 4 . This means that g i (v) = c i × v ηi with η i as a non-negative real and c i as a normalizing constant, i = A, B. By normalizing the impact function using its first moment (as specified in 7), we derive the expression:

g i (v) = v ηi E (v ηi ) = β ηi i Γ(1 + η i ) v ηi
In this equation, β i > 0 denotes the parameter of the exponential distribution of volumes, defined on the half-line ]0, +∞[, with i = A, B. Furthermore, the baselines µ A and µ B are defined as a non-homogeneous piecewise linear continuous functions t → Remark 4.2. The model under consideration aligns with the one proposed by [START_REF] Toke | Modelling trades-through in a limit order book using Hawkes processes[END_REF], specifically when g i ≡ 1 where i = A, B.

We proceed by computing the likelihood function L ({N (t)} t≤T ), which represents the joint probability of the observed data as a function of the underlying parameters.

Proposition 4.1 (Likelihood). Consider a marked Hawkes process N A×B = N A , N B with intensity λ A×B = λ A , λ B and compensator Λ A×B = Λ A , Λ B . We denote the probability densities of the marks by (f A , f B ). The log-likelihood of {N A×B (t)} t≤T relative to the unit rate Poisson process ℓ is equal to :

ℓ = T - i∈{A,B}    T 0 µ i (u) du + j∈{A,B} τ j k ′ <T α ij β ij β ηj j Γ(1 + η j ) v ηj k ′ 1 -e -βij (T -τ j k ′ )    + i∈{A,B} τ i k ≤T ln   µ i τ i k + j∈{A,B} τ j k ′ <τ i k α ij β ηj j Γ(1 + η j ) v ηj k ′ e -βij (τ i k -τ j k ′ )    + i∈{A,B} τ i k ≤T ln f i v i k ( 20
)
4 Also referred to as Omori's Law in the context of earthquake modelling.

where

f i (v) = β i e -βiv , i = A, B.
Proof. The computation of the log-likelihood of a multidimensional Hawkes process involves summing up the likelihood of each coordinate i ∈ {A, B} :

ℓ := ln L ({N (t)} t≤T ) L 0 = T + i∈{A,B} ln L N i (t) t≤T (21)
According to Theorem 7.3.III of [START_REF] Daley | An Introduction to the Theory of Point Processes[END_REF], one can express the likelihood of a realization

τ i 1 , v i 1 ,. . . , (τ i N i g (T ) , v i N i g (T )
) of N i (t) in the following form :

L {N i (t)} t≤T =   N i g (T ) k=1 λ i g τ i k     N i g (T ) k=1 f i v i k   exp - T 0 λ i g (u)du ( 22
)
The likelihood of the process N A×B relative to the unit rate Poisson process is thus expressed as follows :

ℓ := ln L {N A×B (t)} t≤T L 0 = T + i∈{A,B} ln L N i (t) t≤T = T - T 0 λ A g (u) du + T 0 ln λ A g (u) N A (du × dv) + T 0 ln f A (v) N A (du × dv) - T 0 λ B g (u)du + T 0 ln λ B g (u)N B (du × dv) + T 0 ln f B (v) N B (du × dv) = T - i∈{A,B}    T 0 µ i (u) du + j∈{A,B} τ j k ′ <T α ij β ij β ηj j Γ(1 + η j ) v ηj k ′ 1 -e -βij (T -τ j k ′ )    + i∈{A,B} τ i k ≤T ln   µ i τ i k + j∈{A,B} τ j k ′ <τ i k α ij β ηj j Γ(1 + η j ) v ηj k ′ e -βij (τ i k -τ j k ′ )    + i∈{A,B} τ i k ≤T ln f i v i k = T + ℓ λg + ℓ v ,
where f A (resp. f B ) is the density function relative to the distribution of the volumes and the pre-factor c j is the normalizing constant.

Remark 4.3. Equation 2 allows us to decompose the log-likelihood into the sum of two terms. The first term being generated by the ground intensity and the second by the conditional distribution of marks. Note that the absence of any shared parameter between ℓ λg and ℓ v , i.e, between f i and λ i g allows for the independent maximization of each term.

Goodness-of-fit analysis

We now intend to analyze the quality of fit of our model. We begin by examining the stability of the obtained parameters. Figure 5 displays the boxplots related to the baseline intensity parameters The exogenous part of the intensity function effectively captures the market activity profile with a noticeable U-shaped curve. To gain deeper insights, we further explore the remainder of the kernel. In Figure 6, we present boxplots for the values of the parameters (α ij ) i,j∈{A,B} , (β ij ) i,j∈{A,B} , (η i ) i∈{A,B} , and (β i ) i∈{A,B} related to the endogenous aspect of the kernel, providing valuable information about its self-excitation and mutual-excitation properties. We can observe that the values of β A and β B range between 0.009 and 0.013, which confirms our initial guess shown in Figure 4. Additionally, we notice that mutual excitation effects are weak compared to self-excitation effects, consistent with the findings reported by [START_REF] Toke | Modelling trades-through in a limit order book using Hawkes processes[END_REF]. Furthermore, we observe that the parameters have low variance, and the conditions specified in 2.4 are consistently met with an average branching ratio of 0.83. This indicates the stability of our model and suggests that we are operating in a sub-critical regime.

(µ A i ) 1≤i≤14 and (µ B i ) 1≤i≤14 . B 1 B 2 B 3 B 4 B 5 B 6 B 7 B 8 B 9 B 10 B 11 B 12 B 13 B 14 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 A 1 A 2 A 3 A 4 A 5 A 6 A 7 A 8 A 9 A 10 A 11 A 12 A
To assess the validity of our fit, we perform a residual analysis of our counting processes using the Time-Rescaling theorem (see [START_REF] Daley | An Introduction to the Theory of Point Processes[END_REF]). Proposition 4.2 enables us to recursively compute the distances Λ i τ i k -Λ i τ i k-1 , reducing the complexity of the goodness-of-fit algorithm from O(N 2 ) to O(N ), as outlined in [START_REF] Ozaki | Maximum likelihood estimation of Hawkes self-exciting point processes[END_REF].

Proposition 4.2. Consider a marked D-dimensional Hawkes process with exponential decays and an impact function of power law type, denoted as (N i ) 1≤i≤D . Let τ i k k≥0 represent the event times of the process N i , and N i g denote its related ground process for all i ∈ 1, . . . , D. Then,

V i k := Λ i g τ i k -Λ i g τ i k-1 = τ i k τ i k-1 µ i (s)ds + D j=1 α ij β ij A ij (k -1) × 1 -e -βij (τ i k -τ i k-1 ) + τ i k-1 ≤τ j k ′ <τ i k g j (v j k ′ ) 1 -e -βij (τ i k -τ j k ′ ) (23)
where A ij (k) is defined recursively as

A ij (k) = e -βij (τ i k -τ i k-1 ) A ij (k -1) + τ i k-1 ≤τ j k ′ <τ i k g j v j k ′ e -βij (τ i k -τ j k ′ )
with initial condition:

A ij (0) = 0 and V 1 k k≥0 , V 2 k k≥0 , . . . , V D k k≥0
are D independent sequences of independent identically distributed exponential random variables with unit rate.

We compare the deviation between the distribution of the transformed process and that of a stationary Poisson process of intensity equal to 1 through the Kolmogorov-Smirnov test and verify the independence of its observations through the Ljung-Box test. The Q-Q plot in Figure 7 The outcomes of our calibration for the month of May 2022 indicate that our model successfully passes the Kolmogorov-Smirnov test 50.2% of the time at a significance level of 5%, 60.8% of the time at a significance level of 2.5%, and 82.5% at a significance level of 1%. Moreover, our model consistently passes the Ljung-Box test up to the twentieth term at a significance level of 5%. Taking these factors into account, along with the notably convincing Q-Q plots and the presence of low-variance model parameters that meet stability criteria, we can assert that our model fits the market data effectively.

Detection of shifts in liquidity regimes

In this section, the results of the CUSUM procedure studied in section 3 will be used to make a dayto-day comparative analysis of the liquidity state. Our approach involves selecting a specific day and treating the intensity of the Hawkes processes that model the associated trades-through as the baseline intensity of the financial instrument being analyzed. Sequential hypothesis testing is then performed based on this established intensity. In other words, this amounts to comparing the intensities of the trades-through processes related to two distinct days and thus to comparing the state of liquidity between two given days.

We will use the VSTOXX index and bid-ask spread as a marker/proxy of liquidity to assess the quality of our detection. The primary challenge of implementing our sequential hypothesis testing procedure through backtesting is the unavailability of data on disorder periods. Nonetheless, it is noticeable that the commencement of CME's open-outcry trading phase for major equity index futures at 2:30 PM Paris time, along with significant US macroeconomic news releases such as the ISM Manufacturing Index at 4 PM Paris time, leads to amplified market volatility in Europe. The impact of these events on BNP Paribas' stock prices will be utilized as an indicator to evaluate the detection's delay. We start by fixing the parameters m and ρ as defined in section 3 in order to define the hypotheses to be tested and to determine the average detection delays

D i=1 E N i ( T C ) and D i=1 E N i ( TC )
which is relative to them. As the Average Run Length is an increasing function of m for all ρ > 0 and as Lorden's criterion is pessimistic (in the sense that it seeks to minimise the worst detection delay), it is preferable to choose m small enough not to have a very large detection delay (i.e., second type error). This said, a too small m will result in an increase of false alarms (i.e. first type error). A compromise must therefore be found between the two. We therefore propose to carry out our hypothesis tests for ρ down = 0.5 and ρ up = 1.5 and m = 5. Figures 2a and2b show that these parameters yield It can be seen that comparing two days where the difference in volatility is quite large implies that the difference in liquidity between two days is detected quite quickly. However, this does not allow us to assess these variations from a more immediate point of view. Another approach would be to compare two relatively close days in order to monitor the evolution of the state of liquidity. This would allow us to better assess the variations in intra-day liquidity. A comparison between 10/05/2022 and 11/05/2022 is therefore suggested. It can be seen visually that there is a fairly consistent correlation between a low spread and a state where we are under the ρ <1 assumption and vice versa (i.e., a high intra-day liquidity). The quality of the detection can also be checked by directly inspecting the number of trades-through that have taken place during the day as shown in figure 11. Remark 4.4. It is worth emphasizing that our methodology involves comparing intensities among counting processes. Consequently, selecting the appropriate reference day becomes imperative when conducting such comparisons.

D i=1 E N i ( T C ) = 75.

Conclusion and discussions

In conclusion, this paper introduces a novel methodology to systematically quantify liquidity within a limit order book and detecting significant changes in its distribution. The use of Marked Hawkes Processes and a CUSUM-based optimal stopping scheme has led to the development of a reliable liquidity proxy with promising applications in financial markets.

The proposed metric holds great potential in enhancing optimal execution algorithms, reducing market impact, and serving as a reliable market liquidity indicator for market makers. The optimality results of the procedure, particularly in the case of a Cox process with simultaneous jumps and a finite time horizon, highlight the robustness of our approach. The tractable formulation of the average detection delay allows for convenient manipulation of the trade-off between the occurrence of false detections and the speed of identifying moments of liquidity disorder.

Moreover, the empirical validation using real market data further reinforces the effectiveness of our methodology, affirming its practicality and relevance in real-world trading scenarios. Future research can build upon these findings to explore its applications. D i=1 N i with a process (see definition A.1 for ρ < 1 and definition A.2 for ρ < 1) with jumps of size 1, in order to extend these results to cases where there are simultaneous jumps. This will allow us to compute the Average Run Length (ARL) and explicit the lower bounds of the Expected Detection Delay (EDD) of the CUSUM procedure. Subsequently, we provide proof that the CUSUM procedure attains this lower bound, establishing its optimality.

Definition A.1. Let τ i,ϵ k k≥1 = τ i k + i D ϵ inf 1≤l,l ′ ≤D inf 0≤τ l j <τ l ′ j ′ ≤τ i k | τ l j -τ l ′ j ′ | k≥1 = τ i k + ϵ i k k≥1
where τ i 0 = 0 and τ i k k≥1 are the event times of the process N i for all i ranging from 1 to D. In this context, we define N i,ϵ as the counting process whose arrival times are τ i,ϵ k k≥1

, i.e, N i,ϵ (t) = k≥1 1 {τ i,ϵ k ≤t} , ∀0 ≤ t ≤ T . We also define F ϵ = (F ϵ t ) t≥0 as the natural filtration associated to the process

D i=1 N i,ϵ .
Remark A.1. Definition A.1 implies that each element N i,ϵ is F-adapted and that, ∀0 ≤ t ≤ T ,

N i,ϵ (t) = k≥1 1 {τ i,ϵ k ≤t} ≤ N i t - 1 D ϵ inf 1≤l,l ′ ≤D inf 0≤τ l j <τ l ′ j ′ ≤T | τ l j -τ l ′ j ′ | ≤ N i (t), a.s
One can clearly observe that N i -N i,ϵ is an F-sub-martingale. According to the Doob-Meyer decomposition, it can be inferred that Λ i -Λ i,ϵ is a predictable process that is non-decreasing and starts from zero. Consequently, ∀0 ≤ t ≤ T , Λ i,ϵ (t) ≤ Λ i (t), a.s Notation A.1. In the subsequent discussion, we will refer to the process t → sup 0≤s≤t U (s) (resp.

t → sup 0≤s≤t U ϵ (s)) as U (resp. U ϵ ) where U ϵ (t) := D i=1 N i,ϵ (t) -β(ρ)Λ i,ϵ (t), ∀0 ≤ t ≤ T . We will also use the notation U ϵ to represent t → U ϵ (t) := sup 0≤s≤t U ϵ (s) -U ϵ (t).
D i=1 N i are noteworthy. Indeed, the set A i,ϵ := ω ∈ Ω : ∃k ∈ N/τ i k+1 (ω) ≤ τ i,ϵ k (ω)
may not be empty. This can be problematic as it can influence the order of arrival times in τ i,ϵ k k≥1

.

However, it can be established through analysis that as ϵ tends to 0 + , the limit of A i,ϵ has a measure of zero. In fact,

ω ∈ A i,ϵ ⇒ ∃k ∈ N : τ i k+1 (ω) ≤ τ i,ϵ k (ω) ⇒ ∃k ∈ N : τ i k+1 (ω) ≤ τ i k (ω) + ϵT Consequently, A i,ϵ ⊂ {ω ∈ Ω : 0 ≤ k ≤ N i (ω, T ) -1/τ i k+1 (ω) -τ i k (ω) ≤ ϵT }
This entails that,

P A i,ϵ ≤ E   k≤N i (T ) P τ i k+1 -τ i k ≤ ϵT | N i (T )   ≤ E   k≤N i (T ) τ i k +ϵT τ i k λ i (s)ds e - τ i k +ϵT τ i k λ i (s)ds  
Once more, utilizing the theorem of dominated convergence yields the following conclusion and the fact that N i is almost surely finite, we can conclude that :

lim ϵ→0 + E   k≤N i (T ) τ i k +ϵT τ i k λ i (s)ds e - τ i k +ϵT τ i k λ i (s)ds   = E   k≤N i (T ) lim ϵ→0 + τ i k +ϵT τ i k λ i (s)ds e - τ i k +ϵT τ i k λ i (s)ds   = 0 Hence, lim ϵ→0 + P A i,ϵ = 0 (24)
We now introduce an intermediate result that will be useful later in the proof of Theorem 3.1.

Lemma 1. For all t ≥ 0 and i ∈ {1, . . . , D},

lim ϵ→0 + E | N i (t) -N i,ϵ (t) | = lim ϵ→0 + E | Λ i (t) -Λ i,ϵ (t) | = 0
Proof. Let 0 ≤ t ≤ T , i ∈ {1, . . . , D} and τ i kt define the maximum arrival time of N i occurring before t. From Definition A.1, it follows that, if τ i,ϵ kt < τ i kt+1 ,

N i (t) -N i,ϵ (t) = 0 if τ i,ϵ kt ≤ t < τ i kt+1 1 if τ i kt ≤ t < τ i,ϵ kt
Hence, we get that

E (N i (t) -N i,ϵ (t))1 A ϵ,c ≤ E 1 {τ i k t ≤t<τ i,ϵ k t } ≤ P 0 ≤ t -τ i kt ≤ ϵ i kt ≤ P 0 ≤ t -τ i kt ≤ i D ϵT
Based on the definition of a doubly stochastic Poisson process (see chapter 7 of [START_REF] Snyder | Doubly Stochastic Poisson-Processes[END_REF]), we have :

P 0 ≤ t -τ i kt ≤ i D ϵT = E P t - i D ϵT ≤ τ i kt ≤ t | λ i (s) : t - i D ϵT ≤ s ≤ t = E   1 k t ! t t-i D ϵT λ i (s)ds kt e -t t-i D ϵT λ i (s)ds  
Consequently, by employing the theorem of dominated convergence, it follows that :

lim ϵ→0 + P 0 ≤ t -τ i kt ≤ i D ϵT = 0
and that :

lim ϵ→0 + E (N i (t) -N i,ϵ (t))1 A ϵ,c = 0 (25)
Since N i,ϵ ≤ N i (see Remark A.1) and using the fact that lim ϵ→0

+ P A i,ϵ = 0 (see Remark A.2), we get that E | N i (t) -N i,ϵ (t) | = E N i (t) -N i,ϵ (t) and that : lim ϵ→0 + E | N i (t) -N i,ϵ (t) | = lim ϵ→0 + E (N i (t) -N i,ϵ (t))1 A ϵ + E (N i (t) -N i,ϵ (t))1 A ϵ,c = 0 Hence, lim ϵ→0 + E | N i (t) -N i,ϵ (t) | = lim ϵ→0 + E | Λ i (t) -Λ i,ϵ (t) | = 0 Lemma 2. Let t ∈ [0, T ] and U (resp. U ϵ ) the CUSUM process of D i=1 N i (resp. D i=1 N i,ϵ ) for ρ < 1. Then, lim ϵ→0 + E y | U ϵ (t) -U (t) | = 0 with E y = E(. | U 0 = y).
Proof. Let (τ k ) k≥1 denote the distinct and ordered jump times of D i=1 N i and τ ϵ k define the maximum arrival time of D i=1 N i,ϵ occurring before τ k+1 for all k ∈ N * . The goal here is to demonstrate the convergence L 1 by bounding | U ϵ -U | by a process that converges to 0 under the L 1 norm. One approach is to evaluate the difference between the processes U and U ϵ at the arrival times of τ k and τ ϵ k , as U and U ϵ are continuous on the intervals ]τ k , τ ϵ k [ and ]τ ϵ k , τ k+1 [ and only exhibit jumps at those times. To initiate the proof, we start by employing an inductive argument to establish boundedness. As discussed in Remark A.2, the Definition A.1 is not completely satisfactory. In order to address the raised concerns, the inequalities established in the subsequent induction hold for events within the set A ϵ,c := Ω\ D i=1 A i,ϵ . Subsequently, we will utilize the aforementioned result in Remark A.2 to facilitate the passage to the limit. For k = 1, since D i=1 N i (resp. D i=1 N i,ϵ ) jumps at the moment τ 1 (resp. τ ϵ 1 ), two cases can arise :

Case 1. sup 0≤u≤τ1

D i=1 N i (u) -β(ρ)Λ i (u) > D i=1 N i (τ 1 ) -β(ρ)Λ i (τ 1 ). This means that D i=1 β(ρ)Λ i (τ 1 ) > D i=1 N i (τ 1 ). Since D i=1 N i,ϵ (τ 1 ) = 0 by definition, it follows that U (τ 1 ) = U ϵ (τ 1 ) = 0 and that : | U (τ 1 ) -U ϵ (τ 1 ) |= 0 Moreover, | U (τ ϵ 1 ) -U ϵ (τ ϵ 1 ) | =| sup u≤τ ϵ 1 D i=1 N i (u) -β(ρ)Λ i (u) -sup u≤τ ϵ 1 D i=1 N i,ϵ (u) -β(ρ)Λ i,ϵ (u) | = D i=1 N i,ϵ (τ ϵ 1 ) -β(ρ)Λ i,ϵ (τ ϵ 1 ) if U ϵ (τ ϵ 1 ) = 0 0 else ≤ β(ρ) D i=1 Λ i (τ ϵ 1 ) -Λ i,ϵ (τ ϵ 1 ) if U ϵ (τ ϵ 1 ) = 0 0 else Case 2. sup 0≤u≤τ1 D i=1 N i (u) -β(ρ)Λ i (u) = D i=1 N i (τ 1 ) -β(ρ)Λ i (τ 1 ).
The evaluation of the difference in τ 1 is straightforward :

| U (τ 1 ) -U ϵ (τ 1 ) | = D i=1 N i (τ 1 ) -β(ρ)Λ i (τ 1 ) ≤ D i=1 N i (τ 1 ) D i=1 β(ρ)Λ i,ϵ (τ ϵ 1 ) > D i=1 N i,ϵ (τ ϵ 1 ). Given that D i=1 N i (τ 1 ) = D i=1 N i,ϵ (τ ϵ 1 ), we can deduce that D i=1 β(ρ)Λ i,ϵ (τ ϵ 1 ) > D i=1 N i (τ 1
) in this case. Thus :

| U (τ ϵ 1 ) -U ϵ (τ ϵ 1 ) | =| sup u≤τ ϵ 1 D i=1 N i (u) -β(ρ)Λ i (u) -sup u≤τ ϵ 1 D i=1 N i,ϵ (u) -β(ρ)Λ i,ϵ (u) | = β(ρ) D i=1 | Λ i (τ 1 ) -Λ i,ϵ (τ ϵ 1 ) | if U ϵ (τ ϵ 1 ) = 0 D i=1 N i (τ 1 ) -β(ρ)Λ i (τ 1 ) else ≤ β(ρ) D i=1 | Λ i (τ 1 ) -Λ i,ϵ (τ ϵ 1 ) | if U ϵ (τ ϵ 1 ) = 0 β(ρ) D i=1 Λ i,ϵ (τ ϵ 1 ) -Λ i (τ 1 )
else Consequently, we have the following bounds for k = 1 :

| U (τ ϵ 1 ) -U ϵ (τ ϵ 1 ) | ≤ β(ρ) D i=1 | Λ i,ϵ (τ ϵ 1 ) -Λ i (τ 1 ) | + | Λ i,ϵ (τ ϵ 1 ) -Λ i (τ ϵ 1 ) | (26) | U (τ 1 ) -U ϵ (τ 1 ) | = D i=1 N i (τ 1 ) -N i,ϵ (τ 1 ) (27) Let k ≥ 1, such that 0 ≤ τ k ≤ T . Assume that: | U (τ ϵ k ) -U ϵ (τ ϵ k ) | ≤ β(ρ) j≤k D i=1 | Λ i,ϵ (τ ϵ j ) -Λ i (τ j ) | + | Λ i,ϵ (τ ϵ j ) -Λ i (τ ϵ j ) | = L ϵ k (28) | U (τ k ) -U ϵ (τ k ) | ≤ j≤k D i=1 N i (τ j ) -N i,ϵ (τ j ) + L ϵ k-1 = M ϵ k (29)
where L ϵ 0 = 0. We will now demonstrate the validity of these bounds for k +1. The proof follows a similar approach to what we have previously seen for the case of k = 1. We can categorize the analysis into two scenarios, depending on whether U attains a supremum at τ k+1 or not.

Case 1.

sup u≤τ k+1 D i=1 N i (u) -β(ρ)Λ i (u) > D i=1 N i (τ k+1 ) -β(ρ)Λ i (τ k+1
). In this situation, the maximum value in U is attained prior to the (k + 1)-th jump. Therefore, there exists k 0 ≤ k such that U (τ k+1 ) = U (τ k0 ). In case U ϵ (τ ϵ k+1 ) = 0, we have that :

U ϵ (τ ϵ k+1 ) ≥ D i=1 N i,ϵ (τ ϵ k0 ) -β(ρ)Λ i,ϵ (τ ϵ k0 ) = U (τ ϵ k+1 ) + β(ρ) D i=1 Λ i (τ k0 ) -Λ i,ϵ (τ ϵ k0 ) U ϵ (τ ϵ k+1 ) -U (τ ϵ k+1 ) ≥ β(ρ) D i=1 Λ i,ϵ (τ ϵ k0 ) -Λ i (τ k0 ) (30) Moreover, U ϵ (τ ϵ k+1 ) = D i=1 N i,ϵ (τ ϵ k+1 ) -β(ρ)Λ i,ϵ (τ ϵ k+1 ) = D i=1 N i (τ k+1 ) -β(ρ)Λ i (τ k+1 ) + β(ρ)Λ i (τ k+1 ) -β(ρ)Λ i,ϵ (τ ϵ k+1 ) ≤ U (τ ϵ k+1 ) + β(ρ) D i=1 Λ i (τ k+1 ) -Λ i,ϵ (τ ϵ k+1 ) i.e, U ϵ (τ ϵ k+1 ) -U (τ ϵ k+1 ) ≤ β(ρ) D i=1 Λ i (τ k+1 ) -Λ i,ϵ (τ ϵ k+1 ) (31)
Therefore, according to equations 30 and 31,

| U (τ ϵ k+1 ) -U ϵ (τ ϵ k+1 ) | ≤ β(ρ) D i=1 | Λ i,ϵ (τ ϵ k0 ) -Λ i (τ k0 ) | + | Λ i (τ k+1 ) -Λ i,ϵ (τ ϵ k+1 ) | if U ϵ (τ ϵ k+1 ) = 0 | U (τ ϵ k ) -U ϵ (τ ϵ k ) | else ≤ L ϵ k+1 if U ϵ (τ ϵ k+1 ) = 0 L ϵ k else ≤ L ϵ k+1
Furthermore, by making use of assumption 28,

| U (τ k+1 ) -U ϵ (τ k+1 ) | =| U (τ ϵ k ) -U ϵ (τ ϵ k ) | ≤ L ϵ k ≤ M ϵ k+1
Thus, the induction has been established for the first case.

Case 2.

sup 0≤u≤τ k+1 D i=1 N i (u) -β(ρ)Λ i (u) = D i=1 N i (τ k+1 ) -β(ρ)Λ i (τ k+1 ). Suppose that U ϵ (τ ϵ k+1 ) > 0, then there exists k 0 ≤ k such that U ϵ (τ ϵ k+1 ) = U ϵ (τ ϵ k0
) and :

U (τ ϵ k+1 ) ≥ D i=1 N i (τ k0 ) -β(ρ)Λ i (τ k0 ) = D i=1 N i,ϵ (τ ϵ k0 ) -β(ρ)Λ i,ϵ (τ ϵ k0 ) + β(ρ)Λ i,ϵ (τ ϵ k0 ) -β(ρ)Λ i (τ k0 ) = U ϵ (τ ϵ k+1 ) + β(ρ) D i=1 Λ i,ϵ (τ ϵ k0 ) -Λ i (τ k0 ) i.e, U (τ ϵ k+1 ) -U ϵ (τ ϵ k+1 ) ≥ β(ρ) D i=1 Λ i (τ k0 ) -Λ i,ϵ (τ ϵ k0 ) (32)
In addition,

U (τ ϵ k+1 ) = D i=1 N i (τ ϵ k+1 ) -β(ρ)Λ i (τ ϵ k+1 ) = D i=1 N i,ϵ (τ ϵ k+1 ) -β(ρ)Λ i,ϵ (τ ϵ k+1 ) + β(ρ)Λ i,ϵ (τ ϵ k+1 ) -β(ρ)Λ i (τ ϵ k+1 ) ≤ U ϵ (τ ϵ k+1 ) + β(ρ) D i=1 Λ i,ϵ (τ ϵ k+1 ) -Λ i (τ ϵ k+1 ) i.e, U (τ ϵ k+1 ) -U ϵ (τ ϵ k+1 ) ≤ β(ρ) D i=1 Λ i,ϵ (τ ϵ k+1 ) -Λ i (τ ϵ k+1 ) (33)
Therefore, based on equations 32 and 33,

| U (τ ϵ k+1 ) -U ϵ (τ ϵ k+1 ) | = | D i=1 N i (τ k+1 ) -N i,ϵ (τ ϵ k+1 ) -β(ρ)Λ i (τ k+1 ) + β(ρ)Λ i,ϵ (τ ϵ k+1 ) | if U ϵ (τ ϵ k+1 ) = 0 β(ρ) D i=1 | Λ i,ϵ (τ ϵ k+1 ) -Λ i (τ ϵ k+1 ) | + | Λ i (τ k0 ) -Λ i,ϵ (τ ϵ k0 ) | else ≤ β(ρ) D i=1 | Λ i,ϵ (τ ϵ k+1 ) -Λ i (τ k+1 ) | if U ϵ (τ ϵ k+1 ) = 0 L ϵ k+1 else ≤ L ϵ k+1 Moreover, U (τ k+1 ) -U ϵ (τ k+1 ) ≤ D i=1 N i (τ k+1 ) -β(ρ)Λ i (τ k+1 ) -N i,ϵ (τ k+1 ) + β(ρ)Λ i,ϵ (τ k+1 ) ≤ D i=1 N i (τ k+1 ) -N i,ϵ (τ k+1 )
and

U ϵ (τ k+1 ) -U (τ k+1 ) ≤ D i=1 N i,ϵ (τ ϵ k0 ) -β(ρ)Λ i,ϵ (τ ϵ k0 ) -N i (τ ϵ k0 ) + β(ρ)Λ i (τ ϵ k0 ) ≤ β(ρ) D i=1 Λ i (τ ϵ k0 ) -Λ i,ϵ (τ ϵ k0 ) Thus, | U (τ k+1 ) -U ϵ (τ k+1 ) | ≤ D i=1 N i (τ k+1 ) -N i,ϵ (τ k+1 ) + β(ρ) D i=1 Λ i (τ ϵ k0 ) -Λ i,ϵ (τ ϵ k0 ) ≤ M ϵ k+1
which completes the recurrence. We proceed to complete the proof by employing a convergence argument. Let t ∈ [0, T ] and τ kt be the largest arrival time of D i=1 N i (t) + N i,ϵ (t) less than or equal to t. Therefore, we have,

| U ϵ (t) -U (t) | =| U (τ kt ) -U ϵ (τ kt ) + D i=1 N i,ϵ (τ kt ) -N i (τ kt ) -β(ρ)Λ i,ϵ (t) + β(ρ)Λ i (t) | ≤ L ϵ kt + M ϵ kt + D i=1 N i,ϵ (τ kt ) -N i (τ kt ) -β(ρ)Λ i,ϵ (t) + β(ρ)Λ i (t)
In addition, according to Lemma 1,

lim ϵ→0 + E y | D i=1 N i (t) -N i,ϵ (t) | 1 A ϵ,c = lim ϵ→0 + E y | D i=1 Λ i (t) -Λ i,ϵ (t) | 1 A ϵ,c = 0 Since D i=1 N i (t) < +∞ (i.e, E y (card ({k ∈ N; τ k ≤ t})) < +∞), we conclude that lim ϵ→0 + E y L ϵ kt + M ϵ kt + D i=1 N i,ϵ (τ kt ) -N i (τ kt ) -β(ρ)Λ i,ϵ (t) + β(ρ)Λ i (t) 1 A ϵ,c = 0
and that lim

ϵ→0 + E y | U ϵ (t) -U (t) | 1 A ϵ,c = 0
As stated in Remark A.2, lim ϵ→0 + P D i=1 A i,ϵ = 0. By utilizing the Cauchy Schwartz inequality, we

get that lim ϵ→0 + E y | U ϵ (t) -U (t) | 1 A ϵ = 0. Hence, lim ϵ→0 + E y | U ϵ (t) -U (t) | = lim ϵ→0 + E y | U ϵ (t) -U (t) | 1 A ϵ + E y | U ϵ (t) -U (t) | 1 A ϵ,c
= 0

We now introduce the functions g m and h m (resp. gm and hm for ρ = 1 ρ ) solutions of the delayed differential equations define in 34 and 35. This will allow us to compute the average delay of the disorder detection5 E D i=1 N i ( T C ) as a function of the performance/threshold criterion m that defines the CUSUM statistic and to determine the Expected Detection Delay

E θ D i=1 N i (T ) -N i (θ) + | F θ .
Proposition A.1 [START_REF] Karoui | Minimax optimality in robust detection of a disorder time in doubly-stochastic poisson processes[END_REF]). Let functions g m and h m the regular finite variations solutions of delayed differential equations (DDE), denoted DDE(β), restricted to the interval [0, m] :

βg ′ m (x) = g m (x) -g m (x -1) + -1 (34)
with g m (0) = 0 (Cauchy problem),

βh ′ m (x) = h m (x + 1) -h m (x) + 1 (35)
with h m (0) = 0 and h ′ m (0) = 0 (Neumann problem). The same properties hold true for the functions gm and hm , solutions of the same system where β is replaced by β, with ρ = 1/ρ, and

β = β(ρ) = β(1/ρ) = β(ρ)/ρ. Then, g m (y) = m y W (z)dz, with g m (0) = m 0 W (z)dz (36) h m (x) = W (m -x) W (m) W ′ (m) - m-x 0 W (y)dy, h m (m -) = W (m) βW ′ (m) (37) gm (y) = ρ m y ρ z W (z)dz, ( 38
) hm (x) = W (m -x) W (m) W ′ (m) - m-x 0 W (z)dz, hm (m -) = W (m) β W ′ (m) (39)
where

W (x) = 1 β ⌊x⌋ k=0 (-1) k k! ((x -k)/β) k exp((x -k)/β) x 0 W (y)dy = ⌊x⌋ k=0 e (x-k)/β k i=0 (-1) j j! ((x -k)/β) j -1 if β > 1 and W (x) = ρρ x W (x) if β < 1.
Proof. Refer to sections 4 and 6 in the work of [START_REF] Karoui | Minimax optimality in robust detection of a disorder time in doubly-stochastic poisson processes[END_REF].

Proof of Theorem 3.1. Let ρ < 1, m ≥ 0, 0 ≤ t ≤ T , U (t) = D i=1 N i (t) -β(ρ)Λ i (t), U (t) = sup 0≤s≤t U (s), U (t) = sup 0≤s≤t U (s) -U (t)
and g m the continuous solution of the equation 34. Since N i -β(ρ)Λ i decreases between two event times of N i , the process U changes only when N i jumps at time t such that U (t) = 0. The jumps of U are therefore negative with sizes less than or equal to D. This means that, ∀0 ≤ t ≤ T ,

U (t) -U (t -) = 1 { U (t)=0} D i=1 N i (t) -β(ρ)Λ i (t) -U (t -) - D i=1 N i (t -) + β(ρ)Λ i (t -)
Because (Λ i ) 1≤i≤D are continuous, we have that

D i=1 Λ i (t) = D i=1 Λ i (t -) and that U (t) -U (t -) = 1 { U (t)=0} D i=1 N i (t) -N i (t -) -U (t -
) . Thus, U and U are solutions of the following stochastic differential equations :

dU (t) = D i=1 1 { U (t)=0} 1 - U (t -) D j=1 N j (t) -N j (t -) dN i (t) (40) d U (t) = - D i=1 1 ∧ U (t -) D j=1 N j (t) -N j (t -) dN i (t) + β(ρ) D i=1 dΛ i (t) (41)
The transformation described in A.1 is formulated in a manner that ensures the jumps of U ϵ depend on one direction only when considering events that are elements of set A ϵ,c , i.e,

d U ϵ (t) = - D i=1 1 ∧ U ϵ (t -) dN i,ϵ (t) + β(ρ) D i=1 dΛ i,ϵ (t)
It should be noted that, for clarity in notation, we have omitted the fact that the expectations evaluated before taking the limit should include

1 A ϵ,c . Let dN i,m, U ϵ (t) = 1 [0,m] U ϵ (t -) dN i,ϵ (t), ∀i ∈ {1, . . . , D}. Applying the itô formula to G m,ϵ (t) := g m ( U ϵ (t)) -g m ( U (0)) + D i=1 N i,m, U ϵ (t) yields the following decomposition : dG m,ϵ (t) = dg m ( U ϵ (t)) + D i=1 dN i,m, U ϵ (t) = β(ρ)g ′ m ( U ϵ (t -)) D i=1 dΛ i,ϵ (t) + g m ( U ϵ (t)) -g m ( U ϵ (t -)) + D i=1 dN i,m, U ϵ (t) = β(ρ)g ′ m ( U ϵ (t -)) D i=1 dΛ i,ϵ (t) - D i=1 g m ( U ϵ (t -)) -g m (( U ϵ (t -) -1) + ) dN i,m, U ϵ (t) + D i=1 dN i,m, U ϵ (t) Since βg ′ m (x) = g m (x) -g m ((x -1) + ) -1 on [0, m] and g m (x) = g ′ m (x) = 0 if x > m or x < 0 : dG m,ϵ (t) = β(ρ)g ′ m ( U ϵ (t -)) D i=1 d Λ i,ϵ (t) -N i,ϵ (t)
Thus, G m,ϵ is a F-martingale. According to Doob's stopping theorem :

E y g m ( U ϵ ( T C )) -g m ( U (0)) + D i=1 N i,m, U ϵ ( T C ) = E y g m ( U ϵ (0)) -E y g m ( U (0)) = 0 (42)
where T C is the following F-stopping time :

T C = inf t ≤ T : U (t) > m 29
The stopping time

T C ensures that U ϵ ( T C ) converges to U ( T C ) such that 1 [0,m] U ϵ ( T - C
) does not nullify. Indeed, by employing a straightforward induction argument, we can demonstrate that U ϵ ≤ U holds over the intervals [τ ϵ k , τ k+1 [ for all k ∈ 1, . . . , D i=1 N i (T ) . Leveraging the fact that there does not exist any k ∈ N and ω ∈ Ω such that T C (ω) = τ k (ω), as a jump in D i=1 N i corresponds to a decrease in U , we get that lim ϵ→0 + P (B ϵ ) = 0, where B

ϵ := ω ∈ Ω : ∃k ∈ N/ T C (ω) ∈]τ k (ω), τ ϵ k (ω)] .
Hence, for sufficiently small ϵ, it follows that lim ϵ→0 + P τ ϵ k ≤ T C < τ k+1 = 1. Drawing upon the results established in Lemma 2, we have lim ϵ→0

+ U ϵ ( T - C ) = U ( T - C ) = m, under L 1 . Consequently, lim ϵ→0 + E y D i=1 N i,m, U ϵ ( T C ) = lim ϵ→0 + E y 1 B ϵ,c D i=1 N i,m, U ϵ ( T C ) = E y D i=1 N i ( T C ) .
Out of continuity of g m and using the theorem of dominated convergence, we have

lim ϵ→0 + E y g m U ϵ ( T C ) = E y g m ( lim ϵ→0 + U ϵ ( T C )) = g m (m) = 0
Taking ϵ to 0 in equation 42, yields the following result :

E y D i=1 N i ( T C ) = g m (y), ∀m ≥ 0
In the previous proof, we introduced the processes (N i,ϵ ) 1≤i≤D because we looked for a transformation U ϵ of the process U with non-simultaneous jumps and that is inferior to the original reflected process U (t) outside of the intervals [τ k , τ ϵ k [ in order to converge leftwards to U at T C with k ∈ 1, . . . , D i=1 N i (T ) . This allowed us to apply the results cited in [START_REF] Karoui | Minimax optimality in robust detection of a disorder time in doubly-stochastic poisson processes[END_REF] and [START_REF] Poor | Quickest Detection[END_REF] to the D i=1 N i,ϵ and U ϵ processes and then go to the limit to obtain the desired results on the original processes. In what follows, we propose to introduce new processes N i,ϵ and Û ϵ to adapt this line of reasoning to the case where ρ > 1.

Definition A.2. Let ( N i,ϵ ) 1≤i≤D be the counting process whose arrival times are :

( τ i,ϵ k ) k≥1 = τ i k -i D ϵ inf 1≤l,l ′ ≤D inf 0≤τ l j <τ l ′ j ′ ≤τ i k | τ l j -τ l ′ j ′ | k≥1 = τ i k -ϵ i k k≥1
where τ i 0 = 0, (τ i k ) k≥1 are the event times of the process N i , ∀i ∈ {1, . . . , D}. In this context, we define N i,ϵ as the counting process whose arrival times are τ i,ϵ k k≥1

, i.e, N i,ϵ (t) = k≥1 1 { τ i,ϵ k ≤t} , ∀0 ≤ t ≤ T . We also define Fϵ = Fϵ t t≥0

as the natural filtration associated to

D i=1 N i,ϵ with Fϵ t = σ( N i,ϵ (s), 1 ≤ i ≤ D; s ≤ t). Remark A.3. It is worth noting that F ϵ t ⊂ F t-ϵT while Fϵ t ⊂ F t+ϵT , for all 0 ≤ t ≤ T . Since ∀1 ≤ i ≤ D, N i is cadlàg, we have lim ϵ→0 + N i,ϵ (t) = lim ϵ→0 + N i (t + ϵ i ) = N i (t + ) = N i (t), a.s.
Therefore, we find almost-sure convergence for the process N i,ϵ and we limit ourselves to convergence in L 1 for the process N i,ϵ .

Remark A.4. This means that N i,ϵ is Fϵ -adapted and that, ∀0 ≤ t ≤ T ,

N i,ϵ (t) = k≥1 1 {τ i,ϵ k ≤t} ≥ N i (t + 1 D ϵ inf 1≤l,l ′ ≤D inf 0≤τ l j <τ l ′ j ′ ≤T | τ l j -τ l ′ j ′ |) ≥ N i (t), a.s
Note that unlike N i -N i,ϵ (see A.1), N i -N i,ϵ is an F ϵ -sub-martingale. According to the Doob-Meyer decomposition, it can be inferred that Λ i -Λ i,ϵ is a predictable process that is decreasing and starts from zero. Consequently, ∀0 ≤ t ≤ T , Λ i (t) ≤ Λ i,ϵ (t), a.s Notation A.2. In the subsequent discussion, we shall refer to the process inf 0≤s≤t U (s) (resp. inf 0≤s≤t U ϵ (s)) as U (t) (resp. U ϵ (t)) where, in the case ρ > 1, U ϵ (t) := D i=1 N i,ϵ (t)-β(ρ) Λ i,ϵ (t), ∀0 ≤ t ≤ T . We will also use the notation

Û ϵ to represent Û ϵ (t) := D i=1 N i,ϵ (t)+β(ρ) Λi,ϵ (t)-inf 0≤s≤t D i=1 N i,ϵ (s) -β(ρ) Λi,ϵ (s) .
i=1 N i (t) (resp. D i=1 N i,ϵ (t)) for ρ > 1. Thus, lim ϵ→0 + | Û ϵ (t) -Û (t) |= 0, a.s.
Proof. The proof is similar to that of the Lemma 2. Indeed, ( U (t), U ϵ (t)) and ( Û ϵ (t), Û (t)) play symmetric roles.

Proof of Theorem 3.

2. Let ρ > 1, 0 ≤ t ≤ T , U (t) = D i=1 N i (t) -β(ρ)Λ i (t)
and h m a continuous solution of the equation 35. We adopt the same notations as in the prior part and introduce the process Û ϵ (t). Note that Û ϵ (t) only increases if a jump takes place and that U ϵ (t) = inf 0≤s≤t D i=1 N i,ϵ (s)β(ρ) Λi,ϵ (s) is continuous. Applying the Itô formula to Û ϵ (t) µ+1 , we get :

d Û ϵ (t) µ+1 = D i=1 Û ϵ (t -) + 1 µ+1 -Û ϵ (t -) µ+1 d N i,ϵ (t) -(1 + µ) Û ϵ (t -) µ β(ρ)d D i=1 Λi,ϵ (t) -dU ϵ (t) Since U ϵ only decreases when Û ϵ (t) = Û ϵ (t -) = 0, then Û ϵ (t -) µ β(ρ) D i=1 d Λi,ϵ (t) -dU ϵ (t) = Û ϵ (t -) µ β(ρ)d D i=1
Λi,ϵ (t) Going to the limit µ → 0, we get :

d Û ϵ (t) = D i=1 d N i,ϵ (t) -1 { Û ϵ (t)>0} β(ρ)d Λi,ϵ (t), a.s Let J m, Û ϵ (t) := lim µ →0 µ≤s≤t 1 { Û ϵ (s-µ)>m> Û ϵ (
E v h m Û ϵ ( T ϵ C ) -h m ( Û (0)) -h m (m -) + D i=1 N i,m, Û ϵ T ϵ C = E v h m ( Û ϵ (0)) -E v h m ( Û (0))
where TC and T ϵ C are the following F-stopping times : TC = inf t ≤ T : Û (t) > m and T ϵ C = TC + ϵT

The stopping time T ϵ C is introduced to ensure that Û ϵ ( T ϵ C ) converges to Û ( TC ) (see Lemma 3) such that the process 1 [0,m] Û ϵ ( T ϵ- C ) does not nullify. This is achieved by maintaining Û ϵ ≤ Û over the intervals [τ ϵ k , τ k+1 [ for all k ∈ 1, . . . , D i=1 N i (T ) . In order to achieve this, it is sufficient to show that lim ϵ→0 + P (C ϵ ) = 0 where C ϵ := ω ∈ Ω : ∃k ∈ N/τ k (ω) < T C (ω) ≤ τ ϵ k+1 (ω) ≤ T ϵ C (ω) . Given that the demonstration bears resemblance to that which was undertaken in Lemma 1, we will not detail it here. Since lim ϵ→0 + Û ϵ ( T ϵ C ) = Û ( TC ) -= m -, under L 1 , by continuity of h m , we have lim ϵ→0 + h m ( Û ϵ ( T ϵ C )) = h m (lim ϵ→0 + Û ϵ ( T ϵ C )) = h m (m -). Moreover, since N i is càdlàg and finite almost-surely, then lim ϵ→0 + N i,ϵ (t) = N i (t) and

lim ϵ→0 + E v D i=1 N i,m, Û ϵ (t) ( T ϵ C ) = lim ϵ→0 + E v D i=1 1 C ϵ,c N i,m, Û ϵ (t) ( T ϵ C ) = E v D i=1 N i ( TC )
Taking ϵ to 0, yields the following result :

E v D i=1 N i ( TC ) = h m (v), ∀m ≥ 0
The remaining part of the proof is derived from the results established previously in this section on the Average Run Length (ARL). The underlying approach involves initially establishing lower bounds on the Lorden criterion for both ρ < 1 and ρ > 1. Subsequently, it is demonstrated that the CUSUM stopping time attains these lower bounds, thereby establishing its optimality.

Proof of Theorem 3.3. Let D θ,ϵ = ω ∈ Ω : ∃k ∈ N/τ i k (ω) ≤ θ < τ ϵ k (ω) and D θ,ϵ,c = Ω\D θ,ϵ , for all θ ∈ R + . We also define the modified criterion C ϵ , for all F-stopping times τ , as : We start with calculations under the probability P and take the primitive of the conditional criterion Γ τ t = E θ 1 D t,ϵ,c τ t d D i=1 N i,ϵ (s) F t with respect to the non-decreasing process ρ -Ū ϵ = ρ Ū ϵ . An integration by parts yields :

E θ τ t Γ τ s dρ -Ū ϵ (s) F t = E θ τ t 1 D s,ϵ,c ]s,τ ] d D i=1 N i,ϵ (u) dρ -Ū ϵ (s) F t = E θ τ t d D i=1 N i,ϵ (u) ρ -Ū ϵ (u -) -ρ -Ū ϵ (t) 1 D u,ϵ,c F t
Since Ĉϵ (τ ) ≥ Γ τ t , we have :

C ϵ (τ )E θ ρ -Ū ϵ (τ ) F t ≥ E θ τ t D i=1 dN i,ϵ (u)ρ -Ū ϵ (u -) 1 D u,ϵ,c F t + E θ C ϵ (τ )ρ -Ū ϵ (t) - τ t D i=1 dN i,ϵ (u)ρ -Ū ϵ (t) 1 D u,ϵ,c F t ≥ E θ τ t D i=1 dN i,ϵ (u)ρ -Ū ϵ (u -) 1 D u,ϵ,c F t + E θ ρ -Ū ϵ (t) C ϵ (τ ) - τ t d D i=1 N i,ϵ (u) F t ≥ E θ τ t D i=1 dN i,ϵ (u)ρ -Ū ϵ (u -) 1 D u,ϵ,c F t
This yields the following inequalities for all F-stopping times τ :

ρE θ τ t D i=1 dN i,ϵ (u)ρ -U ϵ (u -) 1 D u,ϵ,c F t ≤ C ϵ (τ ) E θ ρ -U ϵ (τ ) F t
Taking the expectation on both sides and using the fact that ρ U is the martingale density of P with respect to P, we establish the following lower bound :

ρE τ t D i=1 dN i,ϵ (u)ρ -U ϵ (u -) ρ U ϵ (u -)-U (u -) 1 D u -,ϵ,c ≤ C ϵ (τ ) E ρ -U ϵ (τ ) ρ U ϵ (τ )-U (τ ) (45)
Since, ∀τ, θ ≥ 0, D i=1 E θ N i,ϵ (τ ) -N i,ϵ (θ)

+ 1 D θ,ϵ,c F θ = D i=1 E θ N i,ϵ (τ ) -N i (θ) + 1 D θ,ϵ,c F θ ≤ D i=1 E θ N i (τ ) -N i (θ) + 1 D θ,ϵ,c F θ ≤ D i=1 E θ N i (τ ) -N i (θ) + F θ Hence, C ϵ (τ ) ≤ C(τ )
Note that we can establish the proof for lim ϵ→0 + P D θ,ϵ = 1 for all θ ≥ 0 using a similar approach as shown in 1. Drawing upon the insights presented in equation 45 and leveraging the convergence findings established in Lemma 2, we can derive the ensuing expression:

ρE τ t D i=1 dN i (u)ρ -U (u -) ≤ C(τ )E ρ -U (τ ) (46)
Applying (i) of Theorem 7 in [START_REF] Karoui | Minimax optimality in robust detection of a disorder time in doubly-stochastic poisson processes[END_REF] to the process D i=1 N i,ϵ , we get that :

gm (0)E ρ -U ϵ (τ ) ≤ ρE τ 0 D i=1 dN i,ϵ (u)ρ -U ϵ (u -) (47)
By taking the limit and utilizing equation 46, we obtain : gm (0) ≤ C(τ )

In accordance with Theorem 3.1, we can conclude that the CUSUM stopping time T C achieves the lower bound of the Lorden criterion C. This establishes the optimality of the CUSUM stopping time.

Proof of Theorem 3.4. The proof of this theorem is very similar to proof 3.3.

A.2 Proof of section 4

Proof of Proposition 4.2.

Λ i g τ i k -Λ i g τ i k-1 = τ i k τ i k-1 λ i g (s)ds = τ i k τ i k-1 µ i (s)ds + D j=1 τ j k ′ <τ i k-1 α ij β ij g j v j k ′ e -βij (τ i k-1 -τ j k ′ ) -e -βij (τ i k -τ j k ′ ) + D j=1 τ i k-1 ≤τ j k ′ <τ i k α ij β ij g j v j k ′ 1 -e -βij (τ i k -τ j k ′ ) = τ i k τ i k-1 µ i (s)ds + D j=1 α ij β ij 1 -e -βij (τ i k -τ i k-1 ) τ j k ′ <τ i k-1 g j v j k ′ e -βij (τ i k-1 -τ j k ′ ) + D j=1 τ i k-1 ≤τ j k ′ <τ i k α ij β ij g j v j k ′ 1 -e -βij (τ i k -τ j k ′ )
34

Following [START_REF] Ozaki | Maximum likelihood estimation of Hawkes self-exciting point processes[END_REF], we have that:

A ij (k) = τ j k ′ <τ i k g j v j k ′ e -βij (τ i k -τ j k ′ ) = e -βij (τ i k -τ i k-1 ) A ij (k -1) + τ i k-1 ≤τ j k ′ <τ i k g j v j k ′ e -βij (τ i k -τ j k ′ )
Hence, ∀k ≥ 2 :

Λ i g τ i k -Λ i g τ i k-1 = τ i k τ i k-1 µ i (s)ds + D j=1 α ij β ij A ij (k -1) × 1 -e -βij (τ i k -τ i k-1 ) + τ i k-1 ≤τ j k ′ <τ i k g j v j k ′ 1 -e -βij (τ i k -τ j k ′ )
where A ij (0) = 0, ∀i, j ∈ {1, . . . , D}.

By virtue of the Time-Rescaling theorem (Daley and Vere-Jones ( 2007)), we can draw the conclusion that the elements V 1 k k≥0 , V 2 k k≥0 , . . . , V D k k≥0 are D sequences of independent identically distributed exponential random variables with unit rate.

Figure 1 :

 1 Figure1: Illustration of a trade-through event of limit 2 on the bid side.2 The abbreviation CUSUM stands for CUmulative SUM.

Figures

  Figures 2a and 2b give a visual image of the values taken by Average Run Length (ARL) when varying the ratio ρ and the parameter m.

  Average Run Length for ρ ∈ [0.6, 1[ and m ∈ [1, 15] Average Run Length for ρ ∈ ]1, 1.4]and m ∈[START_REF] Wang | Sequential change-point detection for mutually exciting point processes over networks[END_REF] 15]

  provides an illustrative representation of the distribution of volume of trades-through on the bid and ask sides of the order book for BNP Paribas stock during May 2022.

Figure 4 :

 4 Figure 4: Log-log scale Complementary Cumulative Density Function (CCDF) of the volume of trades-through on the bid (left figure) and ask (right figure) sides from May 2nd, 2022, 09:30:00, to May 31, 2022, 17:00:00. The dashed red lines correspond to the CCDF of a Pareto distribution with parameter η = 0.45, the dashed black line corresponds to the CCDF of a Gaussian distribution, and the dashed green line corresponds to the CCDF of an Exponential distribution with a coefficient equal to 0.01.

  i 14 T over a subdivision of the time interval [0, T ] into half-hour intervals with µ A i > 0 and µ B i > 0 for all i ∈ {1, . . . , 14}. This adaptive approach allows the model to effectively accommodate fluctuations in intraday market activity.

Figure 5 :

 5 Figure 5: Tukey Boxplot of the values of the baseline intensity parameters for bid side (left figure) and ask side (right figure).

Figure 6 :

 6 Figure 6: Tukey Boxplot of the values of the parameters (α ij ) i,j∈{A,B} , (β ij ) i,j∈{A,B} , (η i ) i∈{A,B} and (β i ) i∈{A,B} .

  shows a typical visual example of the fit of our model on the 01/02/2022 trades-through data.

Figure 7 :

 7 Figure 7: A depiction of a Q-Q plot showcasing the relationship between the theoretical quantiles, derived from an exponential distribution with a parameter of 1, and the sample quantiles generated through a Marked Hawkes model applied to Bid (or Ask) trades-through data for the BNP Paribas stock on 01/02/2022.

  97 and D i=1 E N i ( TC ) = 60.4. We initially compare the intensities of the trades through between the reference day 31/08/2022 and the day with the highest volatility 04/03/2022.

Figure 8 :Figure 9 :

 89 Figure 8: Performing hypothesis testing on the trade-through intensities between 31/08/2022 and 04/03/2022, with parameters ρ up = 1.5, ρ down = 0.5, and m = 5. The first figure (upside) presents the reflected processes U and Û . The second figure (downside) provides a comparison of the detection results against the spread.

Figure 10 :

 10 Figure 10: Performing hypothesis testing of the trade-through intensities between 10/05/2022 and 11/05/2022, with parameters ρ up = 1.5, ρ down = 0.5 and m = 5. The first figure (upside) presents the reflected processes U and Û . The second figure (downside) provides a comparison of the detection results against the spread.

Figure 11 :

 11 Figure 11: Number of cumulative trades-through on the first 4 limits during 11/05/2022.

Figures

  Figures 12a and 12b showcase a comparative simulation of the processes D i=1 N i,ϵ and U ϵ with respect to D i=1 N i and U .

Figures

  Figures 13a and 13b showcase a comparative simulation of the processes D i=1 N i,ϵ and Û ϵ with respect to D i=1 N i and Û .

D

  s+µ), Û ϵ (s)= Û ϵ (s -)=m} 6 , ∀0 ≤ t ≤ T . The Itô-Tanaka-Meyer formula applied to H m,ϵ (t) := h m Û ϵ (t) -h m ( Û (0)) + D i=1 N i,m, Û ϵ (t) + h m (m -)J m, Û ϵ (t)6 Also referred to as discontinuous local time.results in the following decomposition :dH m,ϵ (t) = dh m ( Û ϵ (t)) + D i=1 d N i,m, Û ϵ (t) + h m (m -)dJ m, Û ϵ (t) = β(ρ)h ′ m ( Û ϵ (t -))1 { Û ϵ (t -)>0} D i=1 d Λi,ϵ (t) + h m ( Û ϵ (t)) -h m ( Û ϵ (t -ϵ (t -))1 { Û ϵ (t -Û ϵ (t -) + 1) -h m ( Û ϵ (t -)) d N i,ϵ (t) + D i=1 d N i,m, Û ϵ (t) where d N i,m, Û ϵ (t) (t) = 1 [0,m] Û ϵ (t -) d N i,ϵ (t) and m ≥ 0. Since βh ′ m (x) = h m (x + 1) -h m (x) + 1 if x ∈ [0, m], h m (0) = h ′ m (0) = 0 and h m (x) = 0 if x ≥ m, thus : dH m,ϵ (t) = β(ρ)h ′ m ( Û ϵ (t -))1 { Û ϵ (t)>0} D i=1 d Λi,ϵ (t) -N i,ϵ (t)(43)Hence, H m,ϵ is an F-martingale, with Ft t≥0 = Ft+ϵT t≥0 . Therefore, according to Doob's stopping theorem :

E

  θ N i,ϵ (τ ) -N i,ϵ (θ) + 1 D θ,ϵ,c F θ (44)

Moments of transition between liquidity regimes will be defined as instances where the distribution of trade-throughs experiences a change.

Also referred to as doubly stochastic Poisson process.

The same comparison is undertaken only this time between the reference day 31/08/2022 and the day with the lowest volatility index 05/01/2022.

Often referred to as the Average Run Length (ARL).