Bernard Lapeyre

Mathieu Rosenbaum

Emmanuel Bacry

Jean-Philippe Bouchaud

Jim Gatheral

Fabrizio Lillo

Sylvie Mon

tient aussi aux conditions de travail excellentes offertes par le laboratoire MAS.

Introduction et Principaux Résultats 1. Contexte et Motivations

Les économies modernes doivent une partie de leur richesse au bon fonctionnement de leurs marchés financiers. L'industrie financière a vécu trois changements majeurs depuis les années 1930 et il est pertinent de les rappeler ici [START_REF] Kirilenko | Moore's law vs. Murphy's law: algorithmic trading and its discontents[END_REF].

Le premier a trait à la complexité accrue du système financier lui-même. Cela est un corollaire de l'évolution économique et démographique ainsi que la mondialisation des échanges.

Le deuxième développement est l'émergence de la finance quantitative. Les progrès en économie financière ont été nombreux, pour n'en citer que quelques-uns : la théorie d'optimisation du portefeuille de Markowitz, le modèle d'évaluation des actifs financiers (CAPM) de Sharpe, le modèle multi-facteurs du risque (BARRA) de Rosenberg, la formule de Black, Scholes et Merton pour la valorisation des options et son principe de couverture dynamique, ainsi que la théorie de valorisation par martingales de Harrison et Pliska. Ces travaux et d'autres ont fourni le fondement sur lequel repose un pan de l'industrie financière.

Le troisième bouleversement vient des progrès en technologie informatique, fût-ce au niveau matériel (figure 1), logiciel, collection et stockage des données, ou encore la connectivité et les réseaux. La puissance de calcul a rendu possible la résolution de problèmes auparavant insolubles-le problème d'optimisation du portefeuille en est un des premiers exemples. Les avancées logicielles ont diminué les couts des opérations bancaires, et la connectivité rendu possible l' "électronification" des marchés.

Une autre conséquence de la révolution informatique est qu'elle permet l'enregistrement de toutes les transactions financières, voire tous les ordres passés, offrant au scientifique une quantité astronomique de données et rendant possible une approche quasi-expérimentale des marchés. Au centre des transactions électroniques se trouve le carnet d'ordres. C'est là où l'offre et la demande se rencontrent. Le carnet d'ordres répertorie à chaque instant tous les ordres d'achat ou de vente non-exécutés sur un titre, avec priorité selon le prix et le temps de soumission. Cet objet, dont les règles de fonctionnement sont assez simples, est, regardé de près, très complexe. Le but de cette thèse est d'explorer cet objet sous un angle mathématique.

Plus précisément, notre but est triple:

(1) Proposer quelques fondements théoriques sous des hypothèses simples.

(2) Étudier quelques propriétés empiriques des carnets d'ordres, sur des données dites tick by tick.

(3) Montrer le potentiel applicatif de certain modèles de carnet basés sur une classe de processus ponctuel, les processus de Hawkes, en en soulignant les succès et les limites.

Contributions

Nous résumons dans cette section les points-clés de cette thèse, chapitre par chapitre.

Chapitre 2. Le but de ce chapitre est d'explorer les liens entre la description microscopique de la formation des prix (modélisation multi-agents) et l'approche par équations différentielles stochastiques utilisée classiquement pour décrire l'évolution des prix à des échelles de temps macroscopiques.

Nous présentons une étude mathématique du carnet d'ordres comme une chaîne de Markov multidimensionnelle en temps continu et prouvons plusieurs résultats dans le cas de temps d'arrivée poissonniens indépendants.

Nous montrons que la structure des annulations est un facteur important pour l'existence d'une distribution stationnaire pour le carnet d'ordres et la convergence exponentielle envers elle. Nous démontrons aussi, par l'intermédiaire du théorème central limite fonctionnel, que la limite à grande échelle du processus de prix est un mouvement Brownien.

Chapitre 3. Nous illustrons l'analyse théorique du deuxième chapitre par simulation numérique et comparons les résultats aux données de marché. C'est aussi l'occasion pour nous de décrire les données tick by tick et leur traitement.

Chapitre 4. Depuis leur introduction dans [Haw71b], les processus de Hawkes ont été appliqués dans un large éventail de domaines de recherche allant de de la sismologie (à l'origine), au risque de crédit, la contagion financière, et plus récemment la modélisation de la microstructure des marchés. En microstructure des marchés, et en particulier la modélisation de carnet d'ordres, la pertinence de ces processus provient au moins de deux propriétés empiriques du flux des ordres à l'échelle microscopique:

(1) Time-clustering: les ordres ont tendance à arriver par grappes.

(2) Dépendance mutuelle: les flux d'ordres présentent des dépendances croisées non négligeables. Par exemple, les ordres au marché "excitent" les ordres limites et vice versa.

Dans ce chapitre, nous posons un modèle de carnet d'ordres basé sur le processus de Hawkes dans un cadre markovien, et en utilisant des techniques de la théorie des chaînes de Markov et la stabilité stochastique, montrons que le carnet d'ordres est stable et conduit à un prix diffusif limite à 1. INTRODUCTION ET PRINCIPAUX R ÉSULTATS de grandes échelles de temps, généralisant ainsi les résultats du deuxième chapitre.

Chapitre 5. Ce chapitre est axé sur les aspects pratiques de la modélisation stochastique des carnets d'ordres, à savoir l'estimation d'un modèle réaliste de carnet basé sur les processus de Hawkes et son utilisation pour l'évaluation de stratégies de trading algorithmique.

Introduction and Background

The emergence of electronic trading as a major means of trading financial assets makes the study of the order book central to understanding the mechanisms of price formation. In order-driven markets, buy and sell orders are matched continuously subject to price and time priority. The order book is the list of all buy and sell limit orders, with their corresponding price and size, at a given instant of time. Essentially, three types of orders can be submitted:

• Limit order: Specify a price (also called "quote") at which one is willing to buy or sell a certain number of shares; • Market order: Immediately buy or sell a certain number of shares at the best available opposite quote; • Cancellation order: Cancel an existing limit order.

In the literature, "agents" who submit exclusively limit orders are referred to as liquidity providers. Those who submit market orders are referred to as liquidity takers.

Limit orders are stored in the order book until they are either executed against an incoming market order or canceled. The ask price P A (or simply the ask) is the price of the best (i.e. lowest) limit sell order. The bid price P B is the price of the best (i.e. highest) limit buy order. The gap between the bid and the ask

S := P A -P B , (1.1)
is always positive and is called the spread. Prices are not continuous, but rather have a discrete resolution ∆P, the tick, which represents the smallest quantity by which they can change. We define the mid-price as the average between the bid and the ask

P := P A + P B 2 . (1.2)
The price dynamics is the result of the interplay between the incoming order flow and the order book [START_REF] Bouchaud | Statistical properties of stock order books: empirical results and models[END_REF]. Figure 1 is a schematic illustration of this process [START_REF] Ferraris | Equity market impact models[END_REF]. Note that we chose to represent quantities on the bid side of the book by non-positive numbers.

Although in reality orders can have any size, we shall assume in most of the chapter that all orders have a fixed unit size q. This assumption is convenient to carry out our analysis and is, for now, of secondary importance to the problem we are interested in 1 . Throughout this chapter, we may refer to three different "times":

• Physical time (or clock time) in seconds,

• Event time (or tick time): the time counter is incremented by 1 every time an event (i.e. market, limit or cancellation order) occurs, • Trade time (or transaction time): the time counter is incremented every trade (i.e. every market order).

1.1. Related Literature. Order book modelling has been an area of intense research activity in the last decade. The remarkable interest in this area is due to two factors: . Order book schematic illustration: a buy market order arrives and removes liquidity from the ask side, then sell limit orders are submitted and liquidity is restored.

• Widespread use of algorithmic trading in which the order book is the place where offer and demand meet, • Availability of tick by tick data that record every change in the order book and allow precise analysis of the price formation process at the microscopic level.

Schematically, two modelling approaches have been successful in capturing key properties of the order book-at least partially. The first approach, led by economists, models the interactions between rational agents who act strategically: they choose their trading decisions as solutions to individual utility maximization problems (See e.g. [START_REF] Parlour | Limit order markets: a survey[END_REF] and references therein).

In the second approach, proposed by econophysicists, agents are assumed to act randomly. This is sometime referred to as zero-intelligence order book modelling, in the sense that order arrivals and placements are entirely stochastic. The focus here is more on the "mechanistic" aspects of the continuous double auction rather than the strategic interactions between agents. Despite this apparent limitation, zero-intelligence (or statistical) order book models do capture many salient features of real markets (See [DFG + 03, FPZ05] and references therein). Two notable developments in this strand of research are [START_REF] Maslov | Simple model of a limit order-driven market[END_REF] who proposed one of the earliest stochastic order book models, and [CS01] who added the possibility to cancel existing limit orders.

In their seminal paper [START_REF] Smith | Statistical theory of the continuous double auction[END_REF], Smith et al. develop a dynamical statistical order book model under the assumption of independent Poissonian order flows. They provide a thorough analysis of the model using simulation, dimensional analysis and mean field approximation. They study key characteristics of the model, namely:

(1) Price diffusion.

(2) Liquidity characteristics: average depth profile, bid-ask spread, price impact and time and probability to fill a limit order.

One of the most important messages of their analysis is that zero-intelligence order book models are able to produce reasonable market dynamics and liquidity characteristics. Our focus here is on the first point, that is, the convergence of the price process, which is a jump process at the microscopic level, to a diffusive process 2 at macroscopic time scales. The authors in [START_REF] Smith | Statistical theory of the continuous double auction[END_REF] suggest that a diffusive regime is reached. Their argument relies on a mean field approximation. Essentially, this amounts to neglecting the dependence between order fluctuations at adjacent price levels.

Another important paper of interest to us is [START_REF] Cont | A stochastic model for order book dynamics[END_REF]. Cont et al. propose to model the order book dynamics from the vantage point of queuing systems. They remarkably succeed in deriving many conditional probabilities of practical importance such as the probability of an increase in the midprice, of the execution of an order at the bid before the ask quote moves, and of "making the spread". To our knowledge, they are the first to clearly set the problem of stochastic order book modelling in the context of Markov chains, which is a very powerful and well-studied mathematical concept.

1.2. Outline. In this chapter, we build on the models of [START_REF] Cont | A stochastic model for order book dynamics[END_REF] and [START_REF] Smith | Statistical theory of the continuous double auction[END_REF] to present a stylized description of the order book, and derive 2 We mean (abusively) by "diffusive process" or simply "diffusion" the mathematical concept of Brownian motion. several mathematical results in the case of independent Poissonian arrival times. In particular, we show that the cancellation structure is an important factor ensuring the existence of a stationary distribution for the order book and the exponential convergence towards it. We also prove, by means of the functional central limit theorem (FCLT), that the rescaled-centered price process converges to a Brownian motion, which is a new result.

The remainder of the chapter is organized as follows. In section 2, we motivate our approach using an elementary example where the spread is kept constant ("perfect market making"). In sections 3 trough 5, we compute the infinitesimal generator associated with the order book in a general setting, and link the price dynamics to the instantaneous state of the order book. In section 6, we prove that the order book is ergodic-in particular it has a stationary distribution-that it converges to its stationary state exponentially fast, and that the large-scale limit of the price process is a Brownian motion. Our proofs rely on the theory of infinitesimal generators and Foster-Lyapunov stability criteria for Markov chains. We outline an order book simulation algorithm in section 1 and provide a numerical illustration. Finally, section 4 summarizes our results and contains critiques of Markovian order book models.

An Elementary Approximation: Perfect Market Making

We start with the simplest agent-based market model:

• The order book starts in a full state: All limits above P A (0) and below P B (0) are filled with one limit order of unit size q. The spread starts equal to 1 tick; • The flow of market orders is modeled by two independent Poisson processes M + (t) (buy orders) and M -(t) (sell orders) with constant arrival rates (or intensities) λ + and λ -; • There is one liquidity provider, who reacts immediately after a market order arrives so as to maintain the spread constantly equal to 1 tick. He places a limit order on the same side as the market order (i.e. a buy limit order after a buy market order and vice versa) with probability u and on the opposite side with probability 1u.

The mid-price dynamics can be written in the following form

dP(t) = ∆P (dM + (t) -dM -(t))Z, (2.1)
where Z is a Bernoulli random variable

Z = 0 with probability (1 -u), (2.2)
and Z = 1 with probability u.

(2.

3)

The infinitesimal generator 3 L associated with this dynamics is

L f (P) = u λ + ( f (P + ∆P) -f ) + λ -( f (P -∆P) -f ) , (2.5)
where f denotes a test function. It is well known that a continuous limit is obtained under suitable assumptions on the intensity and tick size. Noting that (2.5) can be rewritten as

L f (P) = 1 2 u (λ + + λ -)(∆P) 2 f (P + ∆P) -2 f + f (P -∆P) (∆P) 2 + u (λ + -λ -)∆P f (P + ∆P) -f (P -∆P) 2∆P , (2.6)
and under the following assumptions

u (λ + + λ -)(∆P) 2 -→σ 2 as ∆P → 0, (2.7)
and u (λ +λ -)∆P-→µ as ∆P → 0, (2.8) the generator converges to the classical diffusion operator

σ 2 2 ∂ 2 f ∂P 2 + µ ∂ f ∂P , (2.9)
corresponding to a Brownian motion with drift. This simple case is worked out as an example of the type of limit theorems that we will be interested in in the sequel. One should also note that a more classical approach using the Functional Central limit Theorem (FCLT) as in [START_REF] Billingsley | Convergence of Probability Measures[END_REF] or [START_REF] Whitt | Stochastic-Process Limits[END_REF] yields similar results ; For given fixed values of λ + , λ -and ∆P, the rescaledcentred price process

P(nt) -nµt √ nσ (2.10) 3
The infinitesimal generator of a time-homogeneous Markov process (X(t)) t≥0 is the operator L, if exists, defined to act on sufficiently regular functions f : R n → R, by

L f (x) := lim t↓0 E[ f (X(t))|X(0) = x] -f (x) t .
(2.4)

It provides an analytical tool to study (X(t)).

converges as n → ∞, to a standard Brownian motion (B(t)) where

σ = ∆P (λ + + λ -)u, (2.11)
and µ = ∆P(λ +λ -)u.

(2.12)

Let us also mention that one can easily achieve more complex diffusive limits such as a local volatility model by imposing that the limit is a function of P and t u (λ + + λ -)(∆P) 2 → σ 2 (P, t), (2.13) and u (λ +λ -)∆P → µ(P, t).

(2.14) This is the case if the original intensities are functions of P and t themselves.

Order Book Dynamics

3.1. Model Setup: Poissonian Arrivals, Reference Frame and Boundary Conditions. We now consider the dynamics of a general order book under the assumption of Poissonian arrival times for market orders, limit orders and cancellations. We shall assume that each side of the order book is fully described by a finite number of limits K, ranging from 1 to K ticks away from the best available opposite quote. We will use the notation 4 X(t) := (a(t); b(t)) := (a 1 (t), . . . , a K (t); b 1 (t), . . . , b K (t)) ,

(3.1)

where a := (a 1 , . . . , a K ) designates the ask side of the order book and a i the number of shares available i ticks away from the best opposite quote, and b := (b 1 , . . . , b K ) designates the bid side of the book. By doing so, we adopt the representation described in [START_REF] Cont | A stochastic model for order book dynamics[END_REF] or [START_REF] Smith | Statistical theory of the continuous double auction[END_REF] 5 , but depart slightly from it by adopting a finite moving frame, as we think it is realistic and more convenient when scaling in tick size will be addressed.

Let us now recall the events that may happen:

• arrival of a new market order; • arrival of a new limit order;

• cancellation of an already existing limit order.

These events are described by independent Poisson processes:

4

In what follows, bold notation indicates vector quantities.

5

See also [START_REF] Gatheral | Zero-intelligence realized variance estimation[END_REF] for an interesting discussion.

• M ± (t): arrival of new market order, with intensity λ M + I(a 0) and λ M -I(b 0); • L ± i (t): arrival of a limit order at level i, with intensity λ L ± i ; • C ± i (t): cancellation of a limit order at level i, with intensity λ

C + i a i and λ C - i |b i |.
q is the size of any new incoming order, and the superscript "+" (respectively "-") refers to the ask (respectively bid) side of the book. Note that the intensity of the cancellation process at level i is proportional to the available quantity at that level. That is to say, each order at level i has a lifetime drawn from an exponential distribution with intensity λ C ± i . Note also that buy limit orders L - i (t) arrive below the ask price P A (t), and sell limit orders L + i (t) arrive above the bid price P B (t). We impose constant boundary conditions outside the moving frame of size 2K: Every time the moving frame leaves a price level, the number of shares at that level is set to a ∞ (or b ∞ depending on the side of the book). Our choice of a finite moving frame and constant6 boundary conditions has three motivations. Firstly, it assures that the order book does not empty and that P A , P B are always well defined. Secondly, it keeps the spread S and the increments of P A , P B and P = (P A + P B )/2 bounded-This will be important when addressing the scaling limit of the price. Thirdly, it makes the model Markovian as we do not keep track of the price levels that have been visited (then left) by the moving frame at some prior time. Figure 1 is a representation of the order book using the above notations.

Comparison to Previous

Results and Models. Before we proceed, we would like to recall some results already present in the literature and highlight their differences with respect to our analysis. Smith et al. have already investigated in [START_REF] Smith | Statistical theory of the continuous double auction[END_REF] the scaling properties of some liquidity and price characteristics in a stochastic order book model. These results are summarized in table 1. In the model of Smith et al. [START_REF] Smith | Statistical theory of the continuous double auction[END_REF], orders arrive on an infinite price grid (This is consistent as limit orders arrival rate per price level is finite). Moreover, the arrival rates are independent of the price level, which has the advantage of enabling the analytical predictions summarized in table 1.

a ∞ = 4, b ∞ = -4.
The shape of the order book is such that a(t) = (0, 0, 0, 0, 1, 3, 5, 4, 2) and b(t) = (0, 0, 0, 0, -1, 0, -4, -5, -3). The spread S (t) = 5 ticks. Assume that at time t ′ > t a sell market order dM -(t ′ ) arrives, then a(t ′ ) = (0, 0, 0, 0, 0, 0, 1, 3, 5), b(t ′ ) = (0, 0, 0, 0, 0, 0, -4, -5, -3) and S (t ′ ) = 7. Assume instead that at t ′ > t a buy limit order dL - 1 (t ′ ) arrives one tick away from the best opposite quote, then a(t ′ ) = (1, 3, 5, 4, 2, 4, 4, 4, 4), b(t ′ ) = (-1, 0, 0, 0, -1, 0, -4, -5, -3) and S (t ′ ) = 1.

Quantity

Scaling relation Average asymptotic depth

λ L /λ C Average spread λ M /λ L f (ǫ, ∆P/p c ) Slope of average depth profile (λ L ) 2 /λ M λ C g(ǫ, ∆P/p c ) Price "diffusion" parameter at short time scales (λ M ) 2 λ C /λ L ǫ -0.5 Price "diffusion" parameter at long time scales (λ M ) 2 λ C /λ L ǫ 0.5 Table 1. Results of Smith et al. ǫ := q/(λ M /2λ C
) is a "granularity" parameter that characterizes the effect of discreteness in order sizes, p c := λ M /2λ L is a characteristic price interval, and f and g are slowly varying functions.

We stress that, to our understanding, these results are obtained by meanfield approximations, which assume that the fluctuations at adjacent price levels are independent. This allows fruitful simplifications of the complex dynamics of the order book. In addition, the authors do not characterize the convergence of the coarse-grained price process in the sense of Stochastic Process Limits, nor do they show that the limiting process is precisely a Brownian motion (theorem 4.2).

In the model of Cont el al. [START_REF] Cont | A stochastic model for order book dynamics[END_REF], arrival rates are indexed by the distance to the best opposite quote, which is more realistic. The order book is constrained to a finite price grid [1, P max ] that facilitates the analysis of the Markov chain. Here, we use a combination of the two models in that the arrival rates are not uniformly distributed across prices, and the reference frame is finite but moving. Cont et al. [START_REF] Cont | A stochastic model for order book dynamics[END_REF] have considered the question of the ergodicity of their order book model. We also address this question following a different route, and more importantly to our analysis, exhibit the rate of convergence to the stationary state, which turns out to be the key of the proof of theorem 4.2.

3.3. Evolution of the Order Book. We can write the following coupled SDEs for the quantities of outstanding limit orders in each side of the order book: 7

da i (t) = -        q - i-1 k=1 a k        + dM + (t) + qdL + i (t) -qdC + i (t) + (J M -(a) -a) i dM -(t) + K i=1 (J L - i (a) -a) i dL - i (t) + K i=1 (J C - i (a) -a) i dC - i (t), (3.2) 7
Remember that, by convention, the b i 's are non-positive.

and

db i (t) =        q - i-1 k=1 |b k |        + dM -(t) -qdL - i (t) + qdC - i (t) + (J M + (b) -b) i dM + (t) + K i=1 (J L + i (b) -b) i dL + i (t) + K i=1 (J C + i (b) -b) i dC + i (t), (3.3)
where the J's are shift operators corresponding to the renumbering of the ask side following an event affecting the bid side of the book and vice versa.

For instance the shift operator corresponding to the arrival of a sell market order dM -(t) of size q is 8

J M -(a) =          0, 0, . . . , 0 k times , a 1 , a 2 , . . . , a K-k          , (3.4)
with k := inf{p :

p j=1 |b j | > q} -inf{p : |b p | > 0}. (3.5)
Similar expressions can be derived for the other events affecting the order book.

In the next sections, we will study some general properties of the order book, starting with the generator associated with this 2K-dimensional continuous-time Markov chain.

8

For notational simplicity, we write J M -(a) instead of J M -(a; b) etc. for the shift operators.

Infinitesimal Generator

Let us work out the infinitesimal generator associated with the jump process described above. We have

L f (a; b) = λ M + ( f [a i -(q -A(i -1)) + ] + ; J M + (b) -f ) + K i=1 λ L + i ( f a i + q; J L + i (b) -f ) + K i=1 λ C + i a i ( f a i -q; J C + i (b) -f ) + λ M -f J M -(a); [b i + (q -B(i -1)) + ] --f + K i=1 λ L - i ( f J L - i (a); b i -q -f ) + K i=1 λ C - i |b i |( f J C - i (a); b i + q -f ), (4.1)
where, to ease the notations, we note f (a i ; b) instead of f (a 1 , . . . , a i , . . . , a K ; b) etc. and

x + := max(x, 0), x -:= min(x, 0), x ∈ R.

(4.2)

The operator above, although cumbersome to put in writing, is simple to decipher: a series of standard difference operators corresponding to the "deposition-evaporation" of orders at each limit, combined with the shift operators expressing the moves in the best limits and therefore, in the origins of the frames for the two sides of the order book. Note the coupling of the two sides: the shifts on the a's depend on the b's, and vice versa. More precisely the shifts depend on the profile of the order book on the other side, namely the cumulative depth up to level i defined by

A(i) := i k=1 a k , (4.3)
and

B(i) := i k=1 |b k |, (4.4)
and the generalized inverse functions thereof

A -1 (q ′ ) := inf{p : p j=1 a j > q ′ }, (4.5)
and

B -1 (q ′ ) := inf{p : p j=1 |b j | > q ′ }, (4.6)
where q ′ designates a certain quantity of shares 9 .

Remark 4.1. The index corresponding to the best opposite quote equals the spread S in ticks, that is

i A := A -1 (0) = inf{p : p j=1 a j > 0} = S ∆P := i S , (4.7
)

and i B := B -1 (0) = inf{p : p j=1 |b j | > 0} = S ∆P := i S = i A .
(4.8)

Price Dynamics

We now focus on the dynamics of the best ask and bid prices, denoted by P A (t) and P B (t). One can easily see that they satisfy the following SDEs:

dP A (t) = ∆P[(A -1 (q) -A -1 (0))dM + (t) - K i=1 (A -1 (0) -i) + dL + i (t) + (A -1 (q) -A -1 (0))dC + i A (t)], (5.1)
and

dP B (t) = -∆P[(B -1 (q) -B -1 (0))dM -(t) - K i=1 (B -1 (0) -i) + dL - i (t) + (B -1 (q) -B -1 (0))dC - i B (t)], (5.2)
which describe the various events that affect them: change due to a market order, change due to limit orders inside the spread, and change due to the 9 Note that a more rigorous notation would be A(i, a(t)) and A -1 (q ′ , a(t))

for the depth and inverse depth functions respectively. We drop the dependence on the last variable as it is clear from the context. cancellation of a limit order at the best price. Equivalently, the respective dynamics of the mid-price and the spread are:

dP(t) = ∆P 2 (A -1 (q) -A -1 (0))dM + (t) -(B -1 (q) -B -1 (0))dM -(t) - K i=1 (A -1 (0) -i) + dL + i (t) + K i=1 (B -1 (0) -i) + dL - i (t) + (A -1 (q) -A -1 (0))dC + i A (t) -(B -1 (q) -B -1 (0))dC - i B (t) , (5.3
)

dS (t) = ∆P (A -1 (q) -A -1 (0))dM + (t) + (B -1 (q) -B -1 (0))dM -(t) - K i=1 (A -1 (0) -i) + dL + i (t) - K i=1 (B -1 (0) -i) + dL - i (t) + (A -1 (q) -A -1 (0))dC + i A (t) + (B -1 (q) -B -1 (0))dC - i B (t) .
(5.4)

The equations above are interesting in that they relate in an explicit way the profile of the order book to the size of an increment of the mid-price or the spread, therefore linking the price dynamics to the order flow. For instance the infinitesimal drifts of the mid-price and the spread, conditional on the shape of the order book at time t, are given by:

E [dP(t)|(a; b)] = ∆P 2 (A -1 (q) -A -1 (0))λ M + -(B -1 (q) -B -1 (0))λ M - - K i=1 (A -1 (0) -i) + λ L + i + K i=1 (B -1 (0) -i) + λ L - i + (A -1 (q) -A -1 (0))λ C + i A a i A -(B -1 (q) -B -1 (0))λ C - i B |b i B | dt, ( 5.5)
and

E [dS (t)|(a; b)] = ∆P (A -1 (q) -A -1 (0))λ M + + (B -1 (q) -B -1 (0))λ M - - K i=1 (A -1 (0) -i) + λ L + i - K i=1 (A -1 (0) -i) + λ L - i + (A -1 (q) -A -1 (0))λ C + i A a i A + (B -1 (q) -B -1 (0))λ C - i B |b i B | dt.
(5.6)

Ergodicty and Diffusive Limit

In this section, our interest lies in the following questions:

(1) Is the order book model defined above stable?

(2) What is the stochastic-process limit of the price at large time scales?

The notions of "stability" and "large-scale limit" will be made precise below. We first need some useful definitions from the theory of Markov chains and stochastic stability. Let (Q t ) t≥0 be the Markov transition probability function of the order book at time t, that is

Q t (x, E) := P [X(t) ∈ E|X(0) = x] , t ∈ R + , x ∈ S, E ⊂ S, (6.1)
where S ⊂ Z 2K is the state space of the order book. We recall that a (aperiodic, irreducible) Markov process is ergodic if an invariant probability measure π exists and

lim t→∞ ||Q t (x, .) -π(.)|| = 0, ∀x ∈ S, (6.2)
where ||.|| designates for a signed measure ν the total variation norm 10 defined as ||ν|| := sup

f :| f |<1 |ν( f )| = sup E∈B(S) ν(E) -inf E∈B(S)
ν(E). (6.4)

In (2.32), B(S) is the Borel σ-field generated by S, and for a measurable function f on S, ν( f ) := S f dν.

V-uniform ergodicity. A Markov process is said V-uniformly ergodic if there exists a coercive 11 function V > 1, an invariant distribution π, and constants 0 < r < 1, and R < ∞ such that

||Q t (x, .) -π(.)|| ≤ Rr t V(x), x ∈ S, t ∈ R + .
(6.5)

V-uniform ergodicity can be characterized in terms of the infinitesimal generator of the Markov process. Indeed, it is shown in [START_REF]Markov Chains and Stochastic Stability[END_REF][START_REF]Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes[END_REF] that it is equivalent to the existence of a coercive function V (the "Lyapunov test function") such that LV(x) ≤ -βV(x) + γ, (Geometric drift condition.) (6.6) 10

The convergence in total variation norm implies the more familiar pointwise convergence lim t→∞ |Q t (x, y)π(y)| = 0, x, y ∈ S.

(6.3) Note that since the state space S is countable, one can formulate the results without the need of a "measure-theoretic" framework. We prefer to use this setting as it is more flexible, and can accommodate possible generalizations of these results.

11

That is, a function such that V(x) → ∞ as |x| → ∞.

for some positive constants β and γ. (Theorems 6.1 and 7.1 in [START_REF]Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes[END_REF].) Intuitively, condition (2.34) says that the larger V(X(t)) the stronger X is pulled back towards the center of the state space S. A similar drift condition is available for discrete-time Markov processes (X n ) n∈N and reads

DV(x) ≤ -βV(x) + γI C (x), (6.7)
where D is the drift operator

DV(x) := E[V(X n+1 ) -V(X n )|X n = x]. (6.8)
and C ⊂ S a finite set. (Theorem 16.0.1 in [START_REF]Markov Chains and Stochastic Stability[END_REF].) We refer to [START_REF]Markov Chains and Stochastic Stability[END_REF] for further details.

6.1. Ergodicity of the Order Book and Rate of Convergence to the Stationary State. Of utmost interest is the behavior of the order book in its stationary state. We have the following result:

Theorem 6.1. If λ C = min 1≤i≤K {λ C ± i } > 0, then (X(t)) t≥0 = (a(t); b(t)
) t≥0 is an ergodic Markov process. In particular (X(t)) has a stationary distribution π. Moreover, the rate of convergence of the order book to its stationary state is exponential. That is, there exist r < 1 and R < ∞ such that

||Q t (x, .) -π(.)|| ≤ Rr t V(x), t ∈ R + , x ∈ S.
(6.9)

Proof. Let

V(x) := V(a; b) := K i=1 a i + K i=1 |b i | + q (6.10)
be the total number of shares in the book (+q shares). Using the expression of the infinitesimal generator (4.1) we have

LV(x) ≤ -(λ M + + λ M -)q + K i=1 (λ L + i + λ L - i )q - K i=1 (λ C + i a i + λ C - i |b i |)q + K i=1 λ L + i (i S -i) + a ∞ + K i=1 λ L + i (i S -i) + |b ∞ | (6.11) ≤ -(λ M + + λ M -)q + (Λ L -+ Λ L + )q -λ C qV(x) + K(Λ L -a ∞ + Λ L + |b ∞ |), (6.12)
where

Λ L ± := K i=1 λ L ± i and λ C := min 1≤i≤K {λ C ± i } > 0. (6.13)
The first three terms in the right hand side of inequality (6.11) correspond respectively to the arrival of a market, limit or cancellation order-ignoring the effect of the shift operators. The last two terms are due to shifts occurring after the arrival of a limit order inside the spread. The terms due to shifts occurring after market or cancellation orders (which we do not put in the r.h.s. of (6.11)) are negative, hence the inequality. To obtain inequality (6.12), we used the fact that the spread i S is bounded by K + 1-a consequence of the boundary conditions we impose-and hence (i Si) + is bounded by K.

The drift condition (6.12) can be rewritten as

LV(x) ≤ -βV(x) + γ, (6.14)
for some positive constants β, γ. Inequality (6.14) together with theorem 7.1 in [START_REF]Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes[END_REF] let us assert that (X(t)) is V-uniformly ergodic, hence (6.9).

Corollary 6.1. The spread S (t) = A -1 (0, a(t))∆P = S (X(t)) has a well-defined stationary distribution-This is expected as by construction the spread is bounded by K + 1.

6.2. The Embedded Markov Chain. Let (X n ) denote the embedded Markov chain associated with (X(t)). In event time, the probabilities of each event are "normalized" by the quantity

Λ(x) := λ M + + λ M -+ Λ L + + Λ L -+ K i=1 λ C + i a i + K i=1 λ C - i |b i |. (6.15)
For instance, the probability of a buy market order when the order book is in state x, is

P["Buy market order at time n"|X n-1 = x] := p M + (x) = λ M + Λ(x) . (6.16)
The choice of the test function V(x) = i a i + i b i + q does not yield a geometric drift condition, and more care should be taken to obtain a suitable test function. Let z > 1 be a fixed real number and consider the function 12

V(x) := z i a i + i |b i | := z ϕ(x) . (6.17)

To save notations, we always use the letter V for the test function.

We have Theorem 6.2. (X n ) is V-uniformly ergodic. Hence, there exist r 2 < 1 and R 2 < ∞ such that

||U n (x, .) -ν(.)|| ≤ R 2 r n 2 V(x), n ∈ N, , x ∈ S. (6.18)
where (U n ) n∈N is the transition probability function of (X n ) n∈N and ν its stationary distribution.

Proof.

DV(x) ≤ λ M + Λ(x) (z i a i -q+ i |b i | -V(x)) + j λ L + j Λ(x) (z i a i +q+ i |b i |+K|b ∞ | -V(x)) + j λ C + j a j Λ(x) (z i a i -q+ i |b i | -V(x)) + λ M - Λ(x) (z i a i + i |b i |-q -V(x)) + j λ L - j Λ(x) (z i a i +Ka ∞ + i |b i |+q -V(x)) + j λ C - j |b j | Λ(x) (z i a i + i |b i |-q -V(x)). (6.19)
If we factor out V(x) = z a i + b i in the r.h.s of (6.19), we get

DV(x) V(x) ≤ λ M + + λ M - Λ(x) (z -q -1) + Λ L -+ Λ L - Λ(x) (z q+Kd ∞ -1) + j λ C + j a j + j λ C - j |b j | Λ(x) (z -q -1), (6.20)
where

d ∞ := max{a ∞ , |b ∞ |}. (6.21) Then DV(x) V(x) ≤ λ M + + λ M - λ M + + λ M -+ Λ L + + Λ L -+ λ C ϕ(x) (z -q -1) + Λ L + + Λ L - λ M + + λ M -+ Λ L + + Λ L -+ λ C ϕ(x) (z q+Kd ∞ -1) + λ C ϕ(x) λ M + + λ M -+ Λ L + + Λ L -+ λ C ϕ(x) (z -q -1), (6.22)
with the usual notations

λ C := min λ C ± i and λ C := max λ C ± i . (6.23) Denote the r.h.s of (6.22) B(x). Clearly lim ϕ(x)→∞ B(x) = λ C (z -q -1) λ C < 0, (6.24)
hence there exists A > 0 such that for x ∈ S and ϕ(x

) > A DV(x) V(x) ≤ λ C (z -q -1) 2λ C := -β < 0. (6.25)
Let C denote the finite set

C = {x ∈ S : ϕ(x) = i a i + i b i ≤ A}. (6.26)
We have DV(x) ≤ -βV(x) + γI C (x), (6.27) with γ := max x∈C DV(x). (6.28) Therefore (X n ) n≥0 is V-uniformly ergodic, by theorem 16.0.1 in [START_REF]Markov Chains and Stochastic Stability[END_REF].

The Case of Constant Cancellation

Rates. The proof above can be applied to the case where the cancellation rates do not depend on the state of the order book X ′ (t)-We shall denote the order book X ′ (t) in order to highlight that the assumption of proportional cancellation rates is relaxed. The probability of a cancellation

dC ± i (t) in [t, t + δt] is now P[C ± i (t + δt) -C ± i (t) = 1|X ′ (t) = x ′ ] = λ C ± i δt + o(δt), (6.29) instead of P[C + i (t + δt) -C + i (t) = 1|X ′ (t) = x ′ ] = λ C + i a i (t)δt + o(δt), (6.30) where lim δt→0 o(δt)/δt = 0. Since Λ = λ M + + λ M -+ Λ L + + Λ L -+ K i=1 λ C + i + K i=1 λ C - i
does not depend on x ′ , the analysis of the stability of the continuoustime process (X ′ (t)) and its discrete-time counterpart (X ′ n ) are essentially the same.

We have the following result:

Theorem 6.3. Set Λ C ± := K i=1 λ C ± i and Λ L ± := K i=1 λ L ± i . (6.31)
Under the condition

λ M + + λ M -+ Λ C + + Λ C -> (Λ L + + Λ L -)(1 + Kd ∞ ), (6.32) (X ′ n ) is V-uniformly ergodic. There exist r 3 < 1 and R 3 < ∞ such that ||U ′n (x, .) -ν ′ (.)|| ≤ R 3 r n 3 V(x), n ∈ N, x ∈ S. (6.33)
The same is true for (X ′ (t)).

Proof. Let us prove the result for (X ′ n ). Inequality (6.20) is still valid by the same arguments, but this time the arrival rates are independent of

x ′ DV(x ′ ) V(x ′ ) ≤ λ M + + λ M - Λ (z -q -1) + Λ L + + Λ L - Λ (z q+Kd ∞ -1) + Λ C + + Λ C - Λ (z -q -1). (6.34) Set z =: 1 + ǫ > 1. (6.35) A Taylor expansion in ǫ gives Λ DV(x) V(x) ≤ (λ M + + λ M -)(-qǫ) + (Λ L + + Λ L -)(q + Kd ∞ )ǫ + (Λ C + + Λ C -)(-qǫ) + o(ǫ). (6.36)
For ǫ > 0 small enough, the sign of (6.36) is determined by the quantity

-(λ M + + λ M -) + (Λ L + + Λ L -)(1 + Kd ∞ ) -(Λ C + + Λ C -). (6.37)
Hence, if (6.32) holds DV(x) ≤ -βV(x) for some β > 0, (6.38) and a geometric drift condition is obtained for X ′ .

If for concreteness we set q = 1 share, and all the arrival rates are symmetric and do not depend on i, then condition (6.32) can be rewritten as

λ M + Kλ C > Kλ L (1 + Kd ∞ ). (6.39)
where K is the size of the order book and d ∞ is the depth far away from the mid-price. Note that the above is a sufficient condition for (V-uniform) stability.

6.4. Large-scale Limit of the Price Process. We are now able to answer the main question of this chapter. Let us define the process e(t) ∈ {1, . . . , 2(2K + 1)} which indicates the last event

{M ± , L ± i , C ± i } i∈{1,...,K} ,
that has occurred before time t.

Lemma 6.1. If we append e(t) to the order book (X(t)), we get a Markov process

Y(t) := (X(t), e(t)) (6.40)
which still satisfies the drift condition (2.34).

Proof. Since e(t) takes its values in a finite set, the arguments of the previous sections are valid with minor modifications, and with the test functions

V(y) := q + a i + |b i | + e, (continuous-time setting) (6.41) V(y) := e a i + |b i |+e . (discrete-time setting) (6.42)
The V-uniform ergodicity of (Y(t)) and (Y n ) follows.

Given the state X n-1 of the order book at time n -1 and the event e n , the price increment at time n can be determined. (See equation (5.3).) We define the sequence of random variables

η n := Ψ(X n-1 , e n ) := Φ(Y n , Y n-1 ), (6.43)
as the price increment at time n. Ψ is a deterministic function giving the elementary "price-impact" of event e n on the order book at state X n-1 . Let µ be the stationary distribution of (Y n ), and M its transition probability function.

We are interested in the random sums

P n := n k=1 η n = n k=1 Φ(Y k , Y k-1 ), (6.44)
where

η k := η k -E µ [η k ] = Φ k = Φ k -E µ [Φ k ], (6.45)
and the asymptotic behavior of the rescaled-centered price process

P (n) (t) := P ⌊nt⌋ √ n , (6.46)
as n goes to infinity.

Theorem 6.4. The series

σ 2 = E µ [η 2 0 ] + 2 ∞ n=1 E µ [η 0 η n ] (6.47)
converges absolutely, and the rescaled-centered price process is a Brownian motion in the limit of n going to infinity. That is

P (n) (t) n→∞ -→ σB(t), (6.48)
where (B(t)) is a standard Brownian motion.

Proof. The idea is to apply the functional central limit theorem for (stationary and ergodic) sequences of weakly dependent random variables with finite variance. Firstly, we note that the variance of the price increments η n is finite since it is bounded by K + 1. Secondly, the V-uniform ergodicity of (Y n ) is equivalent to

||M n (x, .) -µ(.)|| ≤ Rρ n V(x), n ∈ N, (6.49)
for some R < ∞ and ρ < 1. This implies thanks to theorem 16.1.5 in [START_REF]Markov Chains and Stochastic Stability[END_REF] 13 that for any g 2 , h 2 ≤ V, k, n ∈ N, and any initial condition y

|E y [g(Y k )h(Y k+n )] -E y [g(Y k )]E y [h(Y k+n )]| ≤ Rρ n [1 + ρ k V(y)], (6.50) where E y [.] means E[.|Y 0 = y]. This in turn implies |E y [h(Y k )g(Y k+n )]| ≤ R 1 ρ n [1 + ρ k V(y)] (6.51) for some R 1 < ∞, where h = h -E µ [h], g = g -E µ [g]
. By taking the expectation over µ on both sides of (4.17) and noting that E µ [V(Y 0 )] is finite by theorem 14.3.7 in [START_REF]Markov Chains and Stochastic Stability[END_REF], we get

|E µ [h(Y k )g(Y k+n )]| ≤ R 2 ρ n =: ρ(n), k, n ∈ N. (6.52)
Hence the stationary version of (Y n ) satisfies a geometric mixing condition, and in particular

n ρ(n) < ∞. (6.53)
Theorems 19.2 and 19.3 in [START_REF] Billingsley | Convergence of Probability Measures[END_REF] on functions of mixing processes 14 let us conclude that

σ 2 := E µ [η 2 0 ] + 2 ∞ n=1 E µ [η 0 η n ] (6.54)
is well-defined-the series in (4.20) converges absolutely-and coincides with the asymptotic variance

lim n→∞ 1 n E µ        n k=1 (η k ) 2        = σ 2 . (6.55) Moreover P (n) (t) n→∞ -→ σB(t), (6.56)
where (B(t)) is a standard Brownian motion. The convergence in (4.22) happens in D[0, ∞), the space of R-valued càdlàg functions, equipped with the Skorohod topology.

Remark 6.1. In the large-scale limit, the mid-price P, the ask price P A = P + S 2 , and the bid price P B = P -S 2 converge to the same process (σB(t)).

13

We refer to §16.1.2 "V-geometric mixing and V-uniform ergodicity" in [MT09] for more details. 14 See also theorem 4.4.1 in [START_REF] Whitt | Stochastic-Process Limits[END_REF] and discussion therein. Remark 6.2. Theorem 4.2 is also true with constant cancellation rates under condition (6.32). In this case the result holds both in event time and physical time. Indeed, let (N(t)) t∈R + denote a Poisson process with intensity

Λ = λ M ± + Λ L ± + K i=1 λ C ± i .
The price process in physical time (P c (t)) t∈R + can be linked to the price in event time (P n ) n∈N by P c (t) = P N(t) .

(6.57)

Then P ⌊kt⌋ √ k k→∞ -→ σB(t)
as in theorem 4.2, (6.58)

and since N(u) Λu u→∞ -→ 1 a.s., P c (kt) √ k = P N(kt) √ k k→∞ ∼ → P ⌊Λkt⌋ √ k k→∞ -→ √ ΛσB(t).
(6.59) Remark 6.3. Yet another specification of the cancellation process. Another interesting specification of the cancellation process (C i (t)) is to assume that the arrival rate is constant (for each i) but the canceled volume is proportional to the queue size |X i |. In this case, the treatments of the continuous time chain and its embedded discrete-time counterpart are equivalent, and theorems 6.1-4.2 can be obtained in an analogous manner to the proofs in this section.

CHAPTER 3

Markovian Order Book Models II: Numerical Example

Numerical Example

In order to gain a better intuitive understanding of the "mechanics" of the model, we sketch in Algorithm 1 below the simulation procedure in pseudo-code (See also [START_REF] Gatheral | Zero-intelligence realized variance estimation[END_REF] for a similar description.) For simplicity, we take a symmetric order book model. We also let (usual notations):

λ L := λ L 1 , . . . , λ L K , (1.1)
Λ L := K i=1 λ L i , (1.2)
λ C (a) := λ C 1 a 1 , . . . , λ C K a K , (1.3) Λ C (a) := K i=1 λ C i a i , (1.4) λ C (b) := λ C 1 |b 1 |, . . . , λ C K |b K | , (1.5) Λ C (b) := K i=1 λ C i |b i |, (1.6) Λ(a, b) := 2(λ M + Λ L ) + Λ C (a) + Λ C (b). (1.7)
In order to put the simulation results and the data on the same footing, we relax the assumption of constant order sizes; we draw the order volumes from lognormal distributions.

The parameters of the model are estimated from tick by tick data as detailed in section 2. For concreteness1 , we use the parameters of the stock SCHN.PA (Schneider Electric) in March 2011 for the plots. They are summarized in tables 3 and 4.

Algorithm 1 Order book simulation.

Require: Model parameters-Arrival rates: λ M , {λ L i } i∈{1,...K} , {λ C i } i∈{1,...K} , order book size: K, reservoirs: a ∞ , b ∞ , volume distribution parameters: (v M , s M ), (v L , s L ), (v C , s C ). Simulation Parameters-Number of time steps: N. Initialization-t ← 0, X(0) ← X init . 1: for time step n = 1, . . . , N, do 2:
Compute the best bid p B and best ask p A .

3:

Compute Λ C (b) = K i=1 λ C i |b i |, i.e
. the weighted sum of shares at price levels from

p A -K to p A -1. 4: Compute Λ C (a) = K i=1 λ C i a i .

5:

Draw the waiting time τ for the next event from an exponential distribution with parameter

Λ(a, b) = 2(λ M + Λ L ) + Λ C (a) + Λ C (b).

6:

Draw a new event according to the probability vector

λ M , λ M , Λ L , Λ L , Λ C (a), Λ C (b) /Λ(a, b).
These probabilities correspond respectively to a buy market order, a sell market order, a buy limit order, a sell limit order, a cancellation of an existing sell order and a cancellation of an existing buy order. 7:

Depending on the event type, draw the order volume from a lognormal distribution with parameters (v M , s M ), (v L , s L ) or (v C , s C ). 8:

If the selected event is a limit order, select the relative price level from {1, 2, . . . , K} according to the probability vector

λ L 1 , . . . , λ L K /Λ L .

9:

If the selected event is a cancellation, select the relative price level at which to cancel an order from {1, 2, . . . , K} according to the probability vector

λ C 1 a 1 , . . . , λ C K a K /Λ C (a). (or λ C (b)/Λ C (b) for the bid side.) 10:
Update the order book state according to the selected event. 11:

Enforce the boundary conditions:

a i = a ∞ , i ≥ K + 1, b i = b ∞ , i ≥ K + 1.
12: Increment the event time n by 1 and the physical time t by τ. 13: end for Remark 1.1. For the practical implementation, it is easier to work with an "absolute" price frame ∆p × {1 . . . L} where L ≫ K.

Figure 1 represents the average depth profile, that is, the average number of outstanding shares at a distance of i ticks from the best opposite price. The agreement between the simulation and the data is fairly good (See panel (a) of figure 10 for a cross-sectional view on CAC 40 stocks.) We also plot the distribution of the spread in figure 2. Note that the simulated distribution is tighter than the actual one (this is systematic and is documented in panel (b) of figure 10.) Figure 3 shows the fast decay of the autocorrelation function of the price increments. Note the high negative autocorrelation of simulated trade prices relatively to the data. In accordance with the theoretical analysis, figures 4-6 illustrate the asymptotic normality of price increments.

The signature plot of the price time series is defined as the variance of price increments at lag h normalized by the lag, that is

σ 2 h := V [P(t + h) -P(t)] h . (1.8)
This function measures the variance of price increments per time unit. It is interesting in that it shows the transition from the variance at small time scales where micro-structure effects dominate, to the long-term variance. By theorem 4.2 2 lim h→∞ σ 2 h = σ 2 , for some fixed value σ.

(1.9)

We verify this numerically in figure 7. Two remarks are in order regarding the signature plot:

Long-term variance-The simulated long-term variance is systematically lower than the variance computed from the data (This is documented in panel (c) of figure 10.) Intuitively, the depth of the order book is expected to increase from the best price towards the center of the book. In the absence of autocorrelation in trade signs, this would cause prices to wander less often far away from the current best as they hit a higher "resistance". We also suspect that actual prices exhibit locally more "drifting phases" than in a (symmetric) Markovian model where the expected price drift is null at all times. An interesting analysis of a simple order book model that allows time-varying arrival rates can be found in [CS03].

Short-term variance-

The signature plot predicted by the model is too high at short time scales relative to the asymptotic variance, especially for 2 Strictly spreaking, we proved the result in event-time.

traded prices. This is classically known to be due to bid-ask bounce. It is however remarkable that the signature plot of actual trade prices looks much flatter compared to the simulation (See figure 7.) This was discovered and discussed in detail by Bouchaud et al. in [START_REF] Bouchaud | Fluctuations and response in financial markets: The subtle nature of random price changes[END_REF], and Lillo and Farmer in [START_REF] Lillo | The long memory of the efficient market[END_REF] (See also [START_REF] Farmer | Market efficiency and the longmemory of supply and demand: Is price impact variable and permanent or fixed and temporary?[END_REF] and [START_REF] Bouchaud | How markets slowly digest changes in supply and demand[END_REF].) They note that actual order signs exhibit positive long-ranged correlations. They also note that actual prices are diffusive-the signature plot is flat-even at small time scales. They solve this apparent paradox by showing that diffusivity results from two opposite effects: autocorrelation in trade signs induces persistence in the prices, just at the exact amount to counterbalance the mean reversion induced by the liquidity stored in the order book. This subtle equilibrium between liquidity takers and liquidity providers which guarantees price diffusivity at short lags, is not accounted for by the bare Markovian order book model we study, and one can speak about anomalous diffusion at short time scales for Markovian order book models [START_REF] Smith | Statistical theory of the continuous double auction[END_REF]. Because of the absence of positive autocorrelation in trade signs in the model, this effect is magnified when one looks at trades. The next paragraph elaborates on this point.

Anomalous Diffusion at Short Time Scales.

A qualitative understanding of the discrepancy between the model and the data signature plots at short time scales can be gleaned with the following heuristic argument. In what follows, we reason in trade time t. Denote by P T r (t) the price of the trade at time t, and α(t) its sign: We assume that the two signs are equally probable (symmetric model). But to make the argument valid for both the model (for which successive trade signs are independent) and the data (for which trade signs exhibit long memory) we do not assume independence of successive trade signs. Let also for a quantity Z ∆Z(t) := Z(t + 1) -Z(t).

α(t) = 1,
(1.12)

We have by definition where P(t -) and S (t -) are respectively the prevailing mid-price and spread just before the trade. From equation (1.13)

P T r (t) = P(t -) + 1 2 α(t)S (t -), (1.13)
σ T r 1 2 := V[∆P T r (t)] = E ∆P T r (t) 2 = E ∆P(t -) 2 + E ∆P(t -)∆(α(t)S (t -)) + 1 4 E ∆(α(t)S (t -)) 2 . (1.14)
The first term in the r.h.s. is the variance of mid-price increments σ 2 1 . The second term represents the covariance of mid-price increments and the trade sign (weighted by the spread) and we assume it is negligible 3 . Let us focus 3 This amounts to neglecting the correlation between trade signs and mid-quote movements, which is justified by the dominance of cancellations and limit orders in comparison to market orders in order book data.

∆(α(t)S (t -)) = α(t + 1)∆S (t -) + S (t -)∆α(t).
(1.15)

Then E ∆(α(t)S (t -)) 2 = E (∆α(t)) 2 E S (t -) 2 + 2E α(t + 1)∆S (t -)S (t -)∆α(t) + E α(t + 1) 2 E (∆S (t -)) 2 . (1.16)
Again, we neglect the cross term 4 in the r.h.s. and we are left with

E ∆(α(t)S (t -)) 2 ≈ E (∆α(t)) 2 E S (t -) 2 + E (∆S (t -)) 2 . (1.17) 4
This time, we are neglecting the correlation between trade signs and spread movements.

E (∆α(t)) 2 = E α(t + 1) 2 + E α(t) 2 -2E [α(t)α(t + 1)] = 2 (1 -ρ 1 (α)) , (1.18)
where ρ 1 (α) is the autocorrelation of trade signs at the first lag.

Finally 5 :

σ T r 1 2 ≈ σ 1 2 + 1 2 (1 -ρ 1 (α)) E S (t -) 2 + 1 4 E (∆S (t -)) 2 .
(1.20)

Two effects are clear from equation (1.20):

5 More generally, after n trades: (1) The trade price variance at short time scales is larger than the midprice variance, (2) Autocorrelation in trade signs dampens this discrepancy. This partially explains 6 why the trades signature plot obtained from the data is flatter than the model predictions: ρ 1 (α) model = 0, while ρ 1 (α) data ≈ 0.6.

σ T r n 2 ≈ σ n 2 + 1 2n (1 -ρ n (α)) E S (t -) 2 . (1.19)
From a modeling perspective, a possible solution to recover the diffusivity even at very short time scales, is to incorporate long-ranged correlation in the order flow. Toth et al. [TLD + 11] have investigated numerically this route using a "ǫ-intelligence" order book model. In this model, market orders signs are long-ranged correlated, that is, in trade time

ρ n (α) = E [α(t + n)α(t)] ∝ n -γ , γ ∈]0, 1[. (1.21) 6
Interestingly, although the arguments that led to (1.20) are rather qualitative, a back of the envelope calculation with E S 2 ∈ [1, 9], gives a difference σ T r 2σ 2 in the range [0.5, 4.5]; which has the same order of magnitude of the values obtained by simulation. And the size of incoming market orders is a fraction f of the volume displayed at the best opposite quote, with f drawn from the distribution

P ξ ( f ) = ξ(1 -f ) ξ-1 , (1.22)
They show that, by fine tuning the additional parameters γ and ξ, one can ensure the diffusive behavior of the price both at mesoscopic (≈ a few trades) and macroscopic (≈ hundred trades) time scales 7 .

7 Note that Toth. el al. [TLD + 11] model the "latent order book", not the actual observable order book. The former represents the intended volume at each price level p, that is, the volume that would be revealed should the price come close to p. So that the interpretation of their parameters, in particular the expected lifetime τ life of an order, does not strictly match ours.

Model Parameters Estimation

2.1. Description of the Data. For reproducibility, we summarize in tables 3 and 4 the parameters used to obtain figures 1-7. These correspond to estimating the model for the stock SCHN.PA (Schneider Electric). Our dataset consists of TRTH8 data for the CAC 40 index constituents in March 2011 (23 trading days). We have tick by tick order book data up to 10 price levels, and trades. A snapshot of these files is given in tables 1 and 2. In order to avoid the diurnal seasonality in trading activity (and the impact of the US market open on European stocks), we somehow arbitrarily restrict our attention to the time window [9 : 30-14 : 00] CET.

Trades and Tick by Tick Data Processing.

As one cannot distinguish market orders from cancellations in tick by tick data, and since the timestamps of the trades and tick by tick data files are asynchronous, we use h := V [P(t + h) -P(t)]/h. y axis unit is tick 2 per trade for panel (a) and tick 2 .second -1 for panel (b). We used a 1,000,000 event simulation run for the model signature plots. Data signature plots are computed separately for each trading day [9 : 30-14 : 00] then averaged across 23 days. For calendar time signature plots, prices are sampled every second using the last tick rule. The inset is a zoom-in. a matching procedure to reconstruct the order book events. In a nutshell, we proceed as follows for each stock and each trading day:

(1) Parse the tick by tick data file to compute order book state variations:

• If the variation is positive (volume at one or more price levels has increased), then label the event as a limit order. • If the variation is negative (volume at one or more price levels has decreased), then label the event as a "likely market order".

• If no variation-this happens when there is just a renumbering in the field "Level" that does not affect the state of the bookdo not count an event.

(2) Parse the trades file and for each trade:

(a) Compare the trade price and volume to likely market orders whose timestamps are in [t T r -∆t, t T r + ∆t], where t T r is the trade timestamp and ∆t is a predefined time window 9 . (b) Match the trade to the first likely market order with the same price and volume and label the corresponding event as a market order-making sure the change in order book state happens at the best price limits. (c) Remaining negative variations are labeled as cancellations.

Doing so, we have an average matching rate of around 85% for CAC 40 stocks. As a byproduct, one gets the sign of each matched trade, that is, whether it is buyer or seller initiated.

Parameters Estimation.

If T be the trading duration of interest each day (T = 4.5 hours-[9 : 30-14 : 00]-in our case.) Then

λ M := #trades 2T , (2.1)
and

λ L i := 1 2T .
(#buy limit orders arriving i tick away from the best opposite quote + #sell lim. orders etc.) . (2.2)

For cancellations, we need to normalize the count by the average number of shares X i at distance i from the best opposite quote:

λ C i := 1 X i 1 2T .
(#cancellation orders in the bid side arriving i tick away from the best opposite quote + #cancellation orders in the ask side etc.) , (2.3) 9

We set ∆t = 3 s for CAC 40 stocks. We found that the median reporting delay for trades is -900 ms: on average, trades are reported 900 milliseconds before the change is recorded in tick by tick data.

K 30 a ∞ 250 b ∞ 250 (v M , s M ) (4.00, 1.19) (v L , s L ) (4.47, 0.83) (v C , s C ) (4.48, 0.82)
λ M ± 0.1237 Table 3. Model parameters for the stock SCHN.PA (Schneider Electric) in March 2011 (23 trading days). Figures 8 and9 are graphical representation of these parameters.

We then average λ M , λ L i and λ L i across 23 trading days to get the final estimates. As for the volumes, we estimate by maximum likelihood the parameters ( v, s) of a lognormal distribution separately for each order type. We depict the parameters in figures 8 and 9.

A Typical Parmaters Set.

λ M = 0.1, (order / second) λ L = 1, (order / second / tick) λ C = 0.2, (order / second / tick) K = 10, (ticks)

q = 1, (share) a ∞ = b ∞ = 5. (shares)
(2.4)

Results for CAC 40 Stocks

In order to get a cross-sectional view of the performance of the model on all CAC 40 stocks, we estimate the parameters separately for each stock and run a 100, 000 event simulation for each parameter set. We then compare in figure 10 the average depth, average spread and the long-term "volatility" measured directly from the data, to those obtained from the simulations. Dashed line is the identity function-It would correspond to a perfect match between model predictions and the data. Solid line is a linear regression Note that despite the good agreement between the average depth profiles (panel (a)), and although the model successfully predicts the relative magnitudes of the long-term variance σ 2 ∞ and the average spread S for different stocks, it tends to systematically underestimate σ 2 ∞ and S . This may be related to the absence of autocorrelation in order signs in the model and the presence of more drifting phases in actual prices than in those obtained by simulation.

z data = b 1 + b 2 z model for each quantity of interest z. i (ticks) X i (shares) λ L ± i 10 3 .λ C ± i 1 276 0.

Conclusions

In the previous two chapters, we analyzed a simple Markovian order book model, in which elementary changes in the price and spread processes are explicitly linked to the instantaneous shape of the order book and the order flow parameters.

Two basic properties were investigated: the ergodicity of the order book and the large-scale limit of the price process. The first property, which we answered positively, is desirable in that it assures the stability of the order book in the long run, and gives a theoretical underpinning to statistical measurements on order book data. The scaling limit of the price process is, as anticipated, a Brownian motion. A key ingredient in this result is the convergence of the order book to its stationary state at an exponential rate, a property equivalent to a geometric mixing condition satisfied by the stationary version of the order book. This short memory effect, plus a constraint on the variance of price increments guarantee a diffusive limit at large time scales. Our assumptions are independent Poissonian order flows, proportional cancellation rates, and the presence of two reservoirs of liquidity K ticks away from the best quotes to guarantee that the spread does not diverge. 10 10 We believe this assumption can be relaxed under a balance condition on the arrival rates. One has however to consider an order book model with finite but unbounded support, and control not only the stability of the spread but also of all the gaps in the book. We believe the results hold for a wide class of Markovian order book models: In general, one can state that price increments in a stable Markovian order book model are aggregationally Gaussian 11 .

In a sense, this could offer a mathematical justification to the Bachelier model of asset prices, from a market microstructure perspective. In reality, the picture is however more subtle: even if the price process is asymptotically diffusive, at short time scales, the model produces stronger anticorrelation in traded prices than what is actually observed in the data. At those time scales, price diffusivity is arguably the result of a balance between persistent liquidity taking and anti-persistent liquidity providing.

We believe however that the approach presented here is interesting for clearly identifying conditions under which the asymptotic normality of price increments holds; and more importantly, for introducing a set of mathematical tools for further investigating the price dynamics in more sophisticated stochastic order book models. Indeed, using the same techniques, we are studying extensions of our results to the case of mutually exciting-and therefore dependent-order flows (point 1 below). This will be published elsewhere.

Our work can naturally be extended in several ways. In the following lines, we suggest some possible avenues to explore.

Firstly, actual order flows exhibit non-negligible cross dependences. As documented in [MT11b], market orders excite limit orders and vice versa. A possible solution for endogenously incorporating these dependences is the use of mutually exciting processes:

λ M (t) = λ M (0) + t 0 ϕ MM (t -s)dN M (s) + t 0 ϕ LM (t -s)dN L (s), (4.1)
and,

λ L (t) = λ L (0) + t 0 ϕ LL (t -s)dN L (s) + t 0 ϕ ML (t -s)dN M (s), (4.2) 11
Rigorously, the convergence to the stationary state has to happen fast enough. That is, with an integrable convergence rate ρ(n) as in (4.19).

This model has the additional advantage of capturing clustering in order arrivals (due to the self-excitation terms ϕ MM and ϕ LL ), and for exponentially decaying kernels12 can be cast into a Markovian setting.

Besides, long-ranged correlation in order signs is a very important feature of the data, as discussed in section 1. Analyzing this mathematically is more difficult since the model is no longer Markovian. Moreover, it is natural to add another source of randomness on the rates themselves, for instance

dλ(t) = θ(λ(t) -λ(t))dt + ν λ(t)dW(t), (4.3)
where λ is a (deterministic) background intensity to account for the Ushaped daily trading activity and θ, ν are the parameters of a CIR process. Such stochastic arrival rates would lead to stochastic volatility in the prices.

Although we argued that the simple Markovian order book model we study is stable and asymptotically diffusive, markets do show signs of fragility quite often and large jumps do occur in actual prices. Understanding how these macroscopic jumps (or departure from equilibrium) arise from events at the order book level, for instance via sudden evaporation of liquidity in one side of the book is much needed.

Finally, richer price dynamics (e.g. fat-tailed return distributions) can be obtained using feedback loops between the arrival rates and the price (or its volatility) as in [START_REF] Preis | Multi-agent-based order book model for financial markets[END_REF].

These extensions may, however, render the model less amenable to mathematical analysis, and we leave the investigation of such interesting (but sometime difficult) questions for future research. In market microstructure, and particularly order book modelling, the relevance of these processes comes at least from two empirical properties of (market, limit and cancellation) order flows at the microscopic level:

(1) Time clustering: order arrivals are highly clustered in time.

(2) Mutual dependence: order flow exhibit non-negligible cross dependences. For instance, as documented in [MT11b], market orders excite limit orders and vice versa.

At the microscopic level, point process-based microstructure models capture by construction the intrinsic discreteness of prices and volumes. A question of interest in this context is the microscopic to macroscopic transition in the price dynamics. This strand of research has attracted a lot of interest of late [AJ13, BDHM13a, BDHM13b, BM13, CdL13, CdL12, HP13, ZRA13].

In this chapter, we cast a Hawkes process-based order book model into a Markovian setting, and using techniques from the theory of Markov chains and stochastic stability [START_REF]Markov Chains and Stochastic Stability[END_REF], show that the order book is stable and leads to a diffusive price limit at large time scales.

1.1. Outline. Section 2 is a distillation of some mathematical results about Hawkes processes and Markov chains stochastic stability. Section 3 63 contains three auxiliary stability results which, apart from their own interest, are useful to prove the stability of the order book. Section 4 is an application to a particular order book model and is the main contribution of the chapter.

1.2. Notations. The following notations appear frequently throughout this chapter, and we recall them here for reference:

• (X n ): discrete-time process, • (X(t)): continuous-time process, • |x| = p i=1 |x i |, • 1, p = {1, 2, . . . p}.

Preliminary Remarks

We collect in this section several definitions and results that are useful for the rest of this chapter. The presentation is rather informal.

Point Processes.

Definition 2.1 (Point process). A point process is an increasing sequence (T n ) n∈N of positive random variables defined on a measurable space (Ω, F , P).

We will restrict our attention to processes that are nonexplosive, that is, for which lim n→∞ T n = ∞. To each realization (T n ) corresponds a counting function (N(t)) t∈R + defined by

N(t) = n if t ∈ [T n , T n+1 [, n ≥ 0.
(2.1) (N(t)) is a right continuous step function with jumps of size 1 and carries the same information as the sequence (T n ), so that (N(t)) is also called a point process.

Definition 2.2 (Multivariate point process). A multivariate point process (or marked point process) is a point process (T n ) for which a random variable X n is associated to each T n . The variables X n take their values in a measurable space (E, E).

We will restrict our attention to the case where E = {1, . . . , M}, M ∈ N * . For each m ∈ {1, . . . , M}, we can define the counting processes

N m (t) = n≥1 I(T n ≤ t)I(X n = i).
(2.2)

We also call the process

N(t) = (N 1 (t), . . . , N M (t))
a multivariate point process.

Definition 2.3 (Intensity of a point process). A point process (N(t)) t∈R + can be completely characterized by its (conditional) intensity function, λ(t), defined as

λ(t) = lim u→0 P [N(t + u) -N(u) = 1|F t ] u , (2.3)
where F t is the history of the process up to time t, that is, the specification of all points in [0, t]. Intuitively

P [N(t + u) -N(u) = 1|F t ] = λ(t) u + o(u),
(2.4)

P [N(t + u) -N(u) = 0|F t ] = 1 -λ(t)u + o(u), (2.5) P [N(t + u) -N(u) > 1|F t ] = o(u). (2.6)
This is naturally extended to the multivariate case by setting for each m ∈ {1, . . . , M}

λ m (t) = lim u→0 P [N m (t + u) -N m (u) = 1|F t ] u .
(2.7) 2.2. Hawkes Processes.

Hawkes Process.

Definition 2.4. A Hawkes process (N(t)) t∈R + is a point process whose intensity is specified by

λ(t) = µ + α t 0 e -β(t-s) dN(s) = µ + α 0≤s i ≤t e -β(t-s i ) ,
(2.8)

for a triplet (µ, α, β) of positive real numbers 1 .

The process thus defined is self-excited: it has a base intensity µ, plus exponentially decaying shocks due to previous jumps. The parameter α characterizes the scale of the excitation and β its decay in time.

Proposition 2.1. The process X(t) = (N(t), λ(t)) is Markov.

λ(t) = µ + t 0 ϕ(t -s)ds,
(2.9) with an unspecified kernel ϕ > 0. But we only consider exponentially decaying kernels in this chapter.

Proof. From a straightforward calculation, we have for any

t 2 > t 1 λ(t 2 ) = µ + α t 2 0 e -β(t 2 -s) dN(s) = µ + α t 1 0 e -β(t 2 -s) dN(s) + α t 2 t 1 e -β(t 2 -s) dN(s) = µ + e -β(t 2 -t 1 ) (λ(t 1 ) -µ) + t 2 t 1 e -β(t 2 -s) dN(s).
(2.10) So that in order to compute λ(t 2 ), we only need to know λ(t 1 ) and {N(t) :

t 1 ≤ t ≤ t 2 }- the information contained in {N(t), λ(t) : 0 ≤ t < t 1 } is irrelevant. Hence P [(N(t 2 ), λ(t 2 )) ∈ A|{N(t), λ(t) : t ∈ [0, t 1 ]}] = P [(N(t 2 ), λ(t 2 )) ∈ A|N(t 1 ), λ(t 1 )] ,
(2.11) for any measurable set A ⊂ N × R + , and X is Markov.

Multivariate Hawkes Process.

Definition 2.5. We say that N = (N 1 , . . . , N M ) is a multivariate Hawkes process when

λ m (t) = µ m + M j=1 α m j t 0
e -β m j (t-s) dN j (s).

(2.12) Proposition 2.2. Let Y i j (t) = α i j t 0 e -β i j (t-s) dN j (s), 1 ≤ i, j ≤ M, and

Y(t) = {Y i j (t)} 1≤i, j≤M . The process X(t) = (N(t), Y(t)) is Markov. Proof. Let t 2 > t 1 . Since Y m j (t 2 ) = e -β m j (t 2 -t 1 ) Y m j (t 1 ) + t 2 t 1 e -β m j (t 2 -s) dN j (s), ( 2
A i j = α i j β i j , 1 ≤ i, j ≤ M. (2.16) If ρ(A) < 1
then there exists a (unique) stationary multivariate point process N(t) = (N 1 (t), . . . , N m (t)) whose intensity is specified as in definition 2.5. ρ(A) is the spectral radius of the matrix A, that is, its largest eigenvalue.

The Embedded Discrete-time Hawkes

Process. Throughout this chapter, we will mostly work with processes sampled in discrete time. We show in this section how to construct a discrete-time version (X n ) n∈N out of a multivariate Hawkes process X(t) = (N(t), Y(t)) t∈R + , where Y is defined by

Y i j (t) = α i j t 0 e -β i j (t-s) dN j (s), 1 ≤ i, j ≤ M, (2.17)
as in proposition 2.2.

First denote (T n ) n≥1 the jump times of the process (and set T 0 = 0), and

X n = X(T n ) = (N(T n ), Y(T n )).
(2.18)

We define E n = E(T n ) ∈ {1, . . . , M} as the mark of the process. The value of E n indicates which component of N(t) has jumped at time T n . We also define the waiting times (τ n ) n≥1 between two successive jumps as

τ n = T n+1 -T n . (2.19) Given that (N n , Y n ) = (ξ, y), (N n+1 , Y n+1
) is generated as follows: Set

τ n+1 = min(τ 1 n+1 , . . . , τ M n+1 ), (2.20)
where conditional on (N n , Y n ) = (ξ, y), the distribution of τ 1 n+1 , . . . , τ M n+1 is that of independent positive random variables whose marginal distributions are determined by hazard rates

h m (t) := µ m + M j=1 y m j e -β i j t , t ≥ 0, 1 ≤ m ≤ M.
(2.21)

Then set

E n+1 = argmin 1≤m≤M τ m n+1 , (2.22
)

N n+1 = (ξ 1 , . . . , ξ E n + 1, . . . , ξ M ), (2.23)
and Y m j n+1 = y m j e -β m j τ n+1 + α m j I(E n+1 = j).

(2.24)

2.4. Drift of a Discrte-time Markov Process.

Definition 2.7. The drift operator D is defined to act on any nonnegative measurable function V by

DV(x) = E[V(X n+1 ) -V(X n )|X n = x].
(2.25)

We will also use the notation

PV(x) = E[V(X n+1 )|X n = x], (2.26) hence DV(x) = PV(x) -V(x).
(2.27)

As will be clear in the next section, the importance of this operator stems from the existence of criteria based on the drift to establish properties of the process. It can be interpreted as the analogous for a process to the derivative for a function2 .

A Digression on Stochastic Stability. Let (X n

) n∈N be a Markov process on a sate space S and (Q n ) n∈N * its transition probability function, that is

Q n (x, A) = P [X n ∈ A|X 0 = x] ,
(2.28) for x ∈ S and A a measurable subset of S.

Ergodicity of a Markov process.

Ergodicity is a strong form of "stability": To rephrase [START_REF]Markov Chains and Stochastic Stability[END_REF], it means that "there is an invariant regime described by a measure π such that if the process starts in this regime (that is, if X 0 has distribution π) then it remains in the regime. And moreover if the process starts in some other regime, then it converges in a strong probabilistic sense with π as a limiting distribution."

Formally, a (aperiodic, irreducible) Markov process is ergodic if an invariant 3 probability measure π exists and

lim n→∞ ||Q n (x, .) -π(.)|| = 0, ∀x ∈ S, (2.30)
where ||.|| designates for a signed measure ν the total variation norm 4 defined as ||ν|| := sup

f :| f |<1 |ν( f )| = sup A∈B(S) ν(A) -inf A∈B(S)
ν(A).

(2.32)

In (2.32), B(S) is the Borel σ-field generated by S, and for a measurable function f on S, ν( f ) := S f dν.

2.5.2. V-uniform ergodicity. We say that a Markov process is V-uniformly ergodic if there exists a coercive 5 function V > 1, an invariant distribution π, and constants 0 < r < 1, and R < ∞ such that

||Q n (x, .) -π(.)|| ≤ Rr n V(x), x ∈ S.
(2.33)

That is, satisfying the invariance equations

π(A) = S π(dx)Q(x, A), A ∈ B(S). (2.29) 4
If the state space S is countable (this is not the case for (X(t), Y(t)) of proposition 2.2.), the convergence in total variation norm implies the more familiar pointwise convergence lim

n→∞ |Q n (x, y) -π(y)| = 0, ∀x, y ∈ S. (2.31) 5
That is, a function such that V(x) → ∞ as |x| → ∞. The condition V > 1 is of course arbitrary and 1 can be replaced by any positive constant. This is a strong form of ergodicity (note the geometric rate of convergence), and it can be characterized in terms of the drift operator D. Indeed, it is shown in [START_REF]Markov Chains and Stochastic Stability[END_REF][START_REF] Meyn | Stability of Markovian processes I: criteria for discrete-time chains[END_REF] that it is equivalent to the existence of a coercive function V (the "Lyapunov test function") such that

DV(x) ≤ -K 1 V(x) + K 2 I C (x) (Geometric drift condition.) (2.34)
for some positive constants K 1 and K 2 , and C ⊂ S a compact set. (Theorem 16.0.1 in [START_REF]Markov Chains and Stochastic Stability[END_REF].) Condition (2.34) is equivalent to

PV(x) ≤ θV(x) + K 3 I C (x) (2.35)
for some 0 < θ < 1. Intuitively, it says that the larger V(X n ) the stronger X is pulled back towards the center of the state space S.

Interestingly, it is possible to develop central limit theorems for functionals of V-uniformly ergodic Markov processes. This will be used to show that the price process in a stable Hawkes process-based order book model is asymptotically diffusive. Before that, we need the following auxiliary results.

Auxiliary Results

V-uniform Ergodicity of the Intensity of a Hawkes Process.

Let (N(t), λ(t)) t∈R + be a Hawkes process with parameters (µ, α, β), and (N n , λ n ) n∈N its embedded discrete-time process as constructed in section 2.3. Proposition 3.1. If α < β, then the process (λ n ) n∈N is V-uniformly ergodic, with V(λ) = e γλ , (3.1) and γ a suitably chosen positive number.

Proof. If τ n = T n+1 -T n be the waiting time between two successive jumps of (X(t)). There holds for t

′ ∈ [T n , T n+1 [, λ(t ′ ) = λ n + (λ n -µ)e -β(t ′ -T n ) . (3.2)
The hazard rate associated to τ n , conditional on

λ n = λ ∈ R + , is h(t) := µ + (λ -µ)e -βt , (3.3)
and the p.d.f. of τ n is

f (t) = h(t)e -t 0 h(s)ds = µ + (λ -µ)e -βt e -µt-λ-µ β (1-e -βt ) . (3.4) Let V(λ) := e γλ (3.5)
be a Lyapunov test function with γ > 0 an arbitrary parameter. Then

E [V (λ n+1 ) |λ n = λ] = ∞ 0 V(λ(t + )) f (t)dt = ∞ 0 V µ + (λ -µ)e -βt + α × µ + (λ -µ)e -βt e -µt-λ-µ β (1-e -βt ) dt = ∞ 0 e γ(µ+(λ-µ)e -βt +α) × µ + (λ -µ)e -βt e -µt-λ-µ β (1-e -βt ) dt (3.6) Hence PV(λ) V(λ) = e -γλ E [V (λ n+1 ) |λ n = λ] = ∞ 0 e -γ(λ-µ)(1-e -βt )+γα µ + (λ -µ)e -βt e -µt-λ-µ β (1-e -βt ) dt = e γα µ ∞ 0 e -γ(1+ 1 β )(λ-µ)(1-e -βt )-µt dt + e γα (λ -µ) ∞ 0 e -γ(1+ 1 β )(λ-µ)(1-e -βt )-(β+µ)t dt.
(3.7)

Using lemma 5.1, we get

PV(λ) V(λ) = e γα µ I (γ + 1 β )(λ -µ), β, µ + e γα (λ -µ) I (γ + 1 β )(λ -µ), β, β + µ , (3.8)
where

I(a, b, c) := ∞ 0 e -a(1-e -bt )-ct dt. (3.9) Then lim λ→∞ PV(λ) V(λ) = 0 + e γα β(γ + 1 β ) = e γα 1 + γβ . (3.10) And lim λ→∞ DV(λ) V(λ) = lim λ→∞ PV(λ) V(λ) -1 = e γα 1 + γβ -1 = e γα -1 -γβ 1 + γβ . (3.11)
A Taylor expansion in γ around 0 yields

lim λ→∞ DV(λ) V(λ) = γ(α -β) + o(γ), (3.12)
which has the sign of αβ. Finally, if α < β, one can choose γ > 0, λ 0 ∈ R * + and κ > 0 such that ∀λ > λ 0 DV(λ) ≤ -κV(λ), (3.13) and the V-uniform ergodicity of (λ n ) follows. Define also α max = max{α i j } 1≤i, j≤M ∈ R + , (3.17) and

V-uniform

β min = min{β i j } 1≤i, j≤M ∈ R * + , β max = max{β i j } 1≤i, j≤M ∈ R * + . (3.18)
We recall that Y(t) = (Y i j ) 1≤i, j≤M is defined by

Y i j (t) = α i j t 0 e -β i j s dN(s). (3.19)
As in the univariate case, let (N n ) and (Y n ) be the discrete time processes

N n = N(T + n ), and Y n = Y(T + n ), (3.20)
sampled at the jump times (T n ) of (X). We have the following stability result for (Y n ).

h i (t) = µ i + M j=1 e -β i j t y i j , 1 ≤ i ≤ M, (3.24)
and

h(t) = M i=1 h i (t). (3.25)
We first note that, conditional on τ n+1 = t, the probability that the next jump is on N i , i ∈ {1, . . . , M}, is

P E n+1 = i|Y n = y, τ n+1 = t = h i (t) h(t) . (3.26)
We have then

E V(Y n+1 )|Y n = y = ∞ 0 M i=1 e γ 1≤k,l≤M (e -β k,l t y kl +I(l=i)α k,l) h i (t) h(t) × h(t)e -t 0 h(s)ds dt = ∞ 0 M i=1 e γ M k=1 α ki +γ 1≤k,l≤M e -β k,l t y kl ×         µ i + M j=1 y i j e -β i j t         e -M k=1 µ k t-1≤k,l≤M (1-e -β kl t ) y kl
β kl dt.

(3.27)

Dividing by V(y) and rearranging the terms we get

PV(y) V(y) = ∞ 0 M i=1 e γ M k=1 α ki e 1≤k,l≤M (γ+ 1 β kl )(1-e 1-β k,l t )y kl         µ i + M j=1 y i j e -β i j t         dt = M i=1 e γ M k=1 α ki µ i I M 2        (γ + 1 β kl )y kl 1≤k,l≤M ; (β kl ) 1≤k,l≤M ; M k=1 µ k        + M i=1 e γ M k=1 α ki M j=1 y i j I M 2        (γ + 1 β kl )y kl 1≤k,l≤n ; β kl 1≤k,l≤M ; n k=1 µ k + β i j        , (3.28)
where

I p (a 1 , . . . , a p ; b 1 , . . . , b p ; c) := ∞ 0 e -a 1 (1-e -b 1 t )-•••-a p (1-e -bpt )-ct dt (3.29)
is defined in lemma 5.2. The first term in the r.h.s of (3.28) vanishes when |y| → ∞ by lemma 5.2. Again using lemma 5.2, as |y| → ∞, ∀1 ≤ i, j ≤ M,

I M 2        (γ + 1 β kl )y kl 1≤k,l≤M ; (β kl ) 1≤k,l≤M ; M k=1 µ k + β i j        ≤ 1 β min 1≤k,l≤M (γ + 1 β kl )y kl ≤ 1 β min (γ + 1 β max ) 1≤k,l≤M y kl .
(3.30)

Hence, the second term in the r.h.s of (3.28) is bounded by

M i=1 e γ M k=1 α ki |y| β min (γ + 1 β max )|y| ≤ e Mα max γ β min (γ + 1 β max )
.

And for large |y| we have

PV(y) V(y) ≤ e Mα max γ β min (γ + 1 β max ) . (3.31)
In order to conclude the proof, it is enough to show that there exists a suitably chosen γ > 0 such that

h(γ) = e Mα max γ β min (γ + 1 β max ) < 1. (3.32)
Minimizing h with respect to γ, the minimum is reached at

γ * = 1 Mα max - 1 β max > 0. (3.33)
and is equal to

h(γ * ) = Mα max β min e 1-Mαmax βmax .
(3.34)

Note that for γ * to be positive (and V to be coercive) we need

α max ≤ β max M , ( 3
then (Y n ) n∈N is V-uniformly ergodic. Remark 3.1. Note that for M = 1 the condition is α β e 1-α β < 1, (3.37) which is satisfied i.i.f. α β < 1. (3.38) (x → x(1 -e x
) is strictly increasing from 0 to 1 on [0, 1]). We get the result in the univariate case.

Remark 3.2. A sufficient condition is

α max < β min M .
(3.39) Remark 3.3. Stability condition (3.36) is not sharp: It is too stringent on the parameters α i, j and β i, j , and we suspect the stationarity condition of proposition 2.4 to be sufficient for V-uniform ergodicty.

V-uniform

Ergodicity of a "Birth-death" Hawkes Process. Let (N 1 (t), N 2 (t)) be a bivariate Hawkes process with intensities:

λ 1 (t) = µ 1 + α 11 t 0 e -β 11 s dN 1 (s) + α 12 t 0 e -β 12 s dN 2 (s), (3.40) λ 2 (t) = µ 2 + α 21 t 0 e -β 21 s dN 1 (s) + α 22 t 0 e -β 22 s dN 2 (s), (3.41)
and define the queue (X(t)) by

• X(t) → X(t) + 1 when N 1 (t) jumps. This happens with (infinitesimal) probability λ 1 (t)dt. • X(t) → X(t) -1 when N 2 (t) jumps and X(t) 0. This happens with probability λ 2 (t)dt.

• X(t) → X(t) -1 with probability λ 3 X(t)dt for a constant λ 3 > 0. This corresponds to a proportional death rate, or in the context of order book modelling, to a proportional cancellation rate.

We also denote by N 3 (t) a counting process with intensity λ 3 X(t) that jumps by 1 when X(t) jumps by -1 due to a "cancellation".

The queue X(t), albeit peculiar, is the building block of the order book model we present in the next section: N 1 represents the flow of limit orders, N 2 that of market orders and N 3 cancellations.

The following result is the key to the proof of the stability of the order book.

Proposition 3.3. Provided β min is large (specified precisely below), (X n , Y n ) is V-uniformly ergodic, where V(x, y) = e ωx+γ 1≤k,l≤2 y kl , (3.42)

and ω > 0 and γ > 0.

Proof. As usual we write This quantity can be made smaller than 1 if β min is large enough, hence the stated result.

PV(x, y) V(x, y) = e ω+γ(α
Remark 3.4. Intuitively, a large β corresponds to a short memory for the process (X n , Y n ).

Application in Order Book Modelling

4.1. Model Setup. We present a stylized order book model whose dynamics is governed by Hawkes processes. We have already discussed at length a similar Poissonian order book model in [START_REF] Abergel | A mathematical approach to order book modeling[END_REF], so the description provided here is brief.

We shall assume that each side of the order book is fully described by a finite number of limits K, ranging from 1 to K ticks away from the best available opposite quote. We use the notation

X(t) = (a(t); b(t)) = (a 1 (t), . . . , a K (t); b 1 (t), . . . , b K (t)) , (4.1)
where a = (a 1 , . . . , a K ) designates the ask side of the order book and a i the number of shares available i ticks away from the best opposite quote, and b = (b 1 , . . . , b K ) designates the bid side of the book.

3 types of events can happen:

• arrival of a new limit order;

• arrival of a new market order;

• cancellation of an already existing limit order.

Arrival of limit and market orders are described by 4 self and mutually exciting Hawkes processes:

• L ± (t): arrival of a limit order, with intensity λ L ± (t);

• M ± (t): arrival of new market order, with intensity λ M ± (t).

Cancellations are modelled by a (doubly stochastic) Poisson process whose intensity is proportional to the number of shares at each side of the order book, that is

λ C ± |x ± |. (4.2)
We denote by q the size of any new incoming order, and the superscript "+" (respectively "-") refers to the ask (respectively bid) side of the book. Buy limit orders L -(t) arrive below the ask price P A (t), and sell limit orders L + (t) arrive above the bid price P B (t).

Once a limit order arrives, its position is chosen randomly from 1 to K. Similarly once a cancellation order arrives, the order to be cancelled is chosen randomly among the outstanding orders. Furthermore, we impose constant boundary conditions outside the moving frame of size 2K: Every time the moving frame leaves a price level, the number of shares at that level is set to a ∞ (or b ∞ depending on the side of the book). a ∞ and b ∞ represent two "reservoirs of liquidity".

Our choice of a finite moving frame and constant boundary conditions has three motivations: Firstly, it assures that the order book does not empty and that P A , P B are always well defined. Secondly, it keeps the spread S = P A -P B and the increments of P A , P B and P = (P A + P B )/2 bounded-This will be important when addressing the scaling limit of the price. Thirdly, it makes the model Markovian as we do not keep track of the price levels that have been visited (then left) by the moving frame at some prior time.

Figure 1 is a schematic representation of the order book.

Stability of the Order Book.

We first specify the notations for the 4-variate Hawkes process. We set

λ i (t) = µ i + 4 j=1 α i j e -β i j s dN j (s), i ∈ 1, 4 , (4.3)
and by convention the index 1 corresponds to L + , 2 to M + , 3 to L -and 4 to M -.

Proposition 4.1. Provided β min is large (specified precisely below), the order book

(X n , Y n ) is V-uniformly ergodic, where V(x, y) = e ω K i=1 x ± i +γ 1≤k,l≤4 y kl , (4.4)
and ω > 0 and γ > 0.

Proof. We follow the same pattern as the proof of proposition 3.3, and only modify it to account for the fact that the order book is formed from

a ∞ = 4, b ∞ = -4.
The shape of the order book is such that a(t) = (0, 0, 0, 0, 1, 3, 5, 4, 2) and b(t) = (0, 0, 0, 0, -1, 0, -4, -5, -3). The spread S (t) = 5 ticks. Assume that at time t ′ > t a sell market order dM -(t ′ ) arrives, then a(t ′ ) = (0, 0, 0, 0, 0, 0, 1, 3, 5), b(t ′ ) = (0, 0, 0, 0, 0, 0, -4, -5, -3) and S (t ′ ) = 7. Assume instead that at t ′ > t a buy limit order dL - 1 (t ′ ) arrives one tick away from the best opposite quote, then a(t ′ ) = (1, 3, 5, 4, 2, 4, 4, 4, 4), b(t ′ ) = (-1, 0, 0, 0, -1, 0, -4, -5, -3) and S (t ′ ) = 1. multiple queues, and the role of the boundary conditions a ∞ and b ∞ :

PV(x, y) V(x, y) ≤ e ωq+ωa ∞ +γ 4 k=1 α k1 ∞ 0 e -1≤k,l≤4 (γ+ 1 β kl )(1-e -βklt )y kl -4 k=1 µ k t-λ C + K k=1 x + k t ×         µ 1 + 4 j=1 y 1 j e -β 1 j t         dt + e -ωq+γ 4 k=1 α k2 ∞ 0 e -1≤k,l≤4 (γ+ 1 β kl )(1-e -βklt )y kl -4 k=1 µ k t-λ C + K k=1 x + k t ×         µ 2 + 4 j=1 y 2 j e -β 2 j t         dt + e -ωq ∞ 0 e -1≤k,l≤4 (γ+ 1 β kl )(1-e -βklt )y kl -4 k=1 µ k t-λ C + K k=1 x + k t        λ C + K k=1 x + k        dt
+ similar terms for the bid side of the book (4.5)

PV(x, y) V(x, y) ≤        e ωq+ωa ∞ +γ 4 k=1 α k1 µ 1 + e -ωq+ωa ∞ +γ 4 k=1 α k2 µ 2 + e -ωq λ C + K k=1 x + k        × ∞ 0 e -1≤k,l≤4 (γ+ 1 β kl )(1-e -βklt )y kl -4 k=1 µ k t-λ C + K k=1 x + k t dt + e ωq+ωa ∞ +γ 4 k=1 α k1 4 j=1 y 1 j ∞ 0 e -1≤k,l≤4 (γ+ 1 β kl )(1-e -βklt )y kl -4 k=1 µ k t-λ C + K k=1 x + k t-β 1 j dt + e -ωq+γ 4 k=1 α k2 4 j=1 y 2 j ∞ 0 e -1≤k,l≤4 (γ+ 1 β kl )(1-e -βklt )y kl -4 k=1 µ k t-λ C + K k=1 x + k t-β 2 j dt
+ similar terms for the bid side of the book. (4.6)

Again, as |x| + |y| → ∞, PV(x, y) V(x, y) ≤ e -wq + 4 β min (γ + 1 β max ) × e ωq+ωa ∞ +γ 4 k=1 α k1 + e -ωq+γ 4 k=1 α k2 + e ωq+ωb ∞ +γ 4 k=1 α k3 + e -ωq+γ 4 k=1 α k4 . (4.7)
This quantity can be made smaller than 1 if β min is large enough, and this concludes the proof of the proposition.

4.3. Large-scale Limit of the Price Process. Given the state (X n-1 , Y n-1 ) of the order book at time n -1 and the event E n , the price increment at time n can be determined. We define the sequence of random variables

η n = Ψ(X n-1 , Y n-1 , E n ), = Φ(Z n , Z n-1 ), (4.8)
as the price increment at time n, where

Z n = (X n , Y n ). (4.9)
Ψ is a deterministic function giving the elementary "price-impact" of event E n on the order book at state X n-1 . Let µ be the stationary distribution of (Z n ), and M its transition probability function. We are interested in the random sums

P n := n k=1 η k = n k=1 Φ(Z k , Z k-1 ), (4.10)
where

η k := η k -E µ [η k ] = Φ k = Φ k -E µ [Φ k ], (4.11)
and the asymptotic behavior of the rescaled-centered price process

P (n) (t) := P ⌊nt⌋ √ n , (4.12)
as n goes to infinity.

Proposition 4.2. In event time, the large-scale limit of the price process is a Brownian motion. Formally, the series

σ 2 = E µ [η 2 0 ] + 2 ∞ n=1 E µ [η 0 η n ] (4.13)
converges absolutely, and

P (n) (t) n→∞ -→ σB(t), (4.14)
where (B(t)) is a standard Brownian motion.

Proof. This is an application of the functional central limit theorem for (stationary and ergodic) sequences of weakly dependent random variables with finite variance, and is identical to the proof of theorem 6.1 in [START_REF] Abergel | A mathematical approach to order book modeling[END_REF]. Firstly, we note that the variance of the price increments η n is finite since it is bounded by

K + 1. Secondly, the V-uniform ergodicity of (Z n ) is equivalent to ||M n (z, .) -µ(.)|| ≤ Rρ n V(z), n ∈ N, (4.15)
for some R < ∞ and ρ < 1. This implies thanks to theorem 16.1.5 in [START_REF]Markov Chains and Stochastic Stability[END_REF] 6 that for any g 2 , h 2 ≤ V, k, n ∈ N, and any initial condition z

|E z [g(Z k )h(Z n+k )] -E z [g(Z k )]E z [h(Z k )]| ≤ Rρ n [1 + ρ k V(z)], (4.16)
where

E z [.] means E[.|Z 0 = z]. This in turn implies |E z [h(Z k )g(Z k+n )]| ≤ R 1 ρ n [1 + ρ k V(z)] (4.17) for some R 1 < ∞, where h = h -E µ [h], g = g -E µ [g]
. By taking the expectation over µ on both sides of (4.17) and noting that E µ [V(Z 0 )] is finite by theorem 14.3.7 in [START_REF]Markov Chains and Stochastic Stability[END_REF], we get

|E µ [h(Z k )g(Z k+n )]| ≤ R 2 ρ n = ρ(n), k, n ∈ N. (4.18)
Hence the stationary version of (Z n ) satisfies a geometric mixing condition, and in particular

n ρ(n) < ∞. (4.19) 6
We refer to §16.1.2 "V-geometric mixing and V-uniform ergodicity" in [MT09] for more details.

Theorems 19.2 and 19.3 in [START_REF] Billingsley | Convergence of Probability Measures[END_REF] on functions of mixing processes let us conclude that

σ 2 := E µ [η 2 0 ] + 2 ∞ n=1 E µ [η 0 η n ] (4.20)
is well-defined-the series in (4.20) converges absolutely-and coincides with the asymptotic variance

lim n→∞ 1 n E µ        n k=1 (η k ) 2        = σ 2 . (4.21) Moreover P (n) (t) n→∞ -→ σB(t), (4.22)
where (B(t)) is a standard Brownian motion. The convergence in (4.22) happens in D[0, ∞), the space of R-valued càdlàg functions, equipped with the Skorohod topology.

Technical Lemmas

Lemma 5.1. Let a, b, c > 0 be three positive real numbers. Then

I(a, b, c) = ∞ 0 e -a(1-e -bt )-ct dt = (-a) -c/b b e -a Γ( c b ) -Γ( c b , -a) , (5.1)
(5.2)

where the Gamma function is defined for all complex numbers p such that R[p] > 0 as (5.6)

Γ(p) = ∞ 0 t p-1 e -t dt, ( 5.3
Proof. This representation and the limits can be obtained with a symbolic computation system such as Mathematica. Lemma 5.2. More generally, if

I p (a 1 , . . . , a p ; b 1 , . . . , b p ; c) = ∞ 0 e -a 1 (1-e -b 1 t )-•••-a p (1-e -bpt )-ct dt,
(5.7) (5.10)

with a i > 0, b i > 0 ∀i ∈ 1, p ,
Hence I p (a; b; c) ≤ 1 b min |a| , as |a| → ∞. (5.11)
CHAPTER 5

Numerical Results and Applications

This chapter is geared towards practical aspects of stochastic order book modelling, namely the identification of a realistic point process-based order book model and its use as a workbench for algorithmic trading strategies assessment.

1. Introduction 1.1. Notations. We note for the intensity of multivariate Hawkes processes,

λ i (t) = µ i + t 0 M j=1 ϕ i j (t -s)dN j (s), 1 ≤ i ≤ M. (1.1)
or in vector form,

λ(t) = µ + t 0 ϕ(t -s)dN(s), (1.2)
and for exponentially decaying kernels,

ϕ i j (t) = P p=1 α p i j e -β p i j t , P ≥ 1, (1.3)
λ i (t) = µ i + t 0 M j=1 P p=1 α p i j e -β p i j (t-s) dN j (s), 1 ≤ i ≤ M.
(1.4) 1.2. Classic Results. Firstly, the stationarity condition for Hawkes processes is

ρ ∞ 0 ϕ(t)dt < 1, (1.5)
where ρ designates the spectral radius of a matrix, that is, the maximum of the absolute values of its eigenvalues. A proof of this result can be found in [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF]. Stationarity was defined in chapter 4 (def. 2.6).

85

Secondly, taking the expectation of (1.2) and inverting, yields the following expression for the average intensity

λ = E [λ(t)] = I - ∞ 0 ϕ(t)dt -1 µ.
(1.6)

Maximum Likelihood Estimation of Multivariate Hawkes Processes

The standard method for estimating the parameters of Hawkes processes is via the maximization of the likelihood function [START_REF] Ozaki | Maximum likelihood estimation of Hawkes self-exciting point processes[END_REF]. We describe the method in this section and apply it to order book data.

First recall that the likelihood of a (regular univariate) point process, for an observation 0 < t 1 , • • • < t n ≤ T , reads (A proof can be found in [START_REF] Snyder | Random point processes in time and space[END_REF] or [START_REF] Rubin | Regular point processes and their detection[END_REF].)

L(t 1 , . . . , t n ) = e -T 0 λ(s)ds n i=1 λ(t i ).
(2.1)

L(t 1 , . . . , t n
) is the joint probability density of observing n jumps exactly at times t 1 , . . . , t n . The product n i=1 λ(t i ) represents the probability of the process jumping at instants t 1 , . . . , t n , while the term e -T 0 λ(s)ds corresponds to the probability of the process staying idle at any other instant in ]0, T ].

The log-likelihood is therefore

ln L(t 1 , . . . , t n ) = - T 0 λ(s)ds + n i=1 ln λ(t i ) = - T 0 λ(s)ds + T 0 ln λ(s)dN(s).
(2.2)

2.1. Likelihood of Univariate Hawkes Processes. Consider now the specific case of a univariate Hawkes process. The intensity reads

λ(t) = µ + t -∞ P p=1
α p e -β p (t-s) dN(s).

(2.3)

We are considering decay kernels of the form

ϕ(t) = P p=1
α p e -β p t .

(2.4)

Although this may seem restrictive at first, we shall show that this parametric form is rich enough to virtually approximate any kernel ϕ of interest (section 6.2).

The computations involved in equation (2.2) can be carried out explicitly. First set P = 1 (Thanks to the linearity of integration, the general case follows readily.)

λ(s) = µ + t k <s αe -β(s-t k ) .
(2.5)

We start with the first integral in (2.2). Assuming t n = T (This can always be made, as in practice we are only given an observation t 1 , . . . , t n and T is not specified [START_REF] Ozaki | Maximum likelihood estimation of Hawkes self-exciting point processes[END_REF].),

T 0 λ(s)ds = µt 1 + n i=2 t i t i-1
λ(s)ds.

(2.6)

t i t i-1 λ(s)ds = µ(t i -t i-1 ) + α t i t i-1 t k <s e -β(s-t k ) ds.
(2.7)

t i t i-1 t k <s e -β(s-t k ) ds = t i t i-1 t k ≤t i-1 e -β(s-t k ) ds = t k ≤t i-1 t i t i-1 e -β(s-t k ) ds = 1 β t k ≤t i-1 e -β(t i-1 -t k ) -e -β(t i -t k ) . (2.8) Then n i=2 t i t i-1 λ(s)ds = µ(t n -t 1 ) + n i=2 α β t k ≤t i-1 e -β(t i-1 -t k ) -e -β(t i -t k ) .
(2.9)

A careful look at the summand shows that the sum in the r.h.s. of (2.9) is telescopic. After simplification, it remains

n i=2 t i t i-1 λ(s)ds = µ(t n -t 1 ) + α β        n i=2 1 - n-1 k=1 e -β(t n -t k )        = µ(t n -t 1 ) + n-1 k=1 α β 1 -e -β(t n -t k ) , (2.10) whence T 0 λ(s)ds = µt n + n-1 k=1 α β 1 -e -β(t n -t k ) = µt n + n k=1 α β 1 -e -β(t n -t k ) . ( 2
λ m (t) = µ m + M j=1 t -∞ P p=1 α p m j e -β p m j (t-s) dN j (s) = µ m + M j=1 t j k <s P p=1 α p m j e -β p m j (t-t j k ) , m ∈ 1, M , (2.18)
where t j k is the k-th jump in the process (N j (t)), j ∈ 1, M .

We will denote the jump times of the multivariate process either as {t j k } 1≤k≤n j ,1≤ j≤M , t j k being the k-th jump in the process (N j (t)). Or equivalently as {t i } 1≤i≤n , in which case the marks {z i } 1≤i≤n ∈ 1, M , indicating which component has jumped at time t i , must also be given.

The log-likelihood of a multivariate Hawkes process is the sum of the log-likelihoods of each component

ln L ({t i , z i } 1≤i≤n ) = M m=1 ln L m ({t i , z i } 1≤i≤n ) , (2.19)
where

ln L m ({t i , z i } 1≤i≤n ) = - T 0 λ m (s)ds + T 0 ln λ m (s)dN m (s).
(2.20)

As in the univariate case, the difficult term to compute is the first integral. We proceed similarly. We write (P = 1 to simplify notations)

T 0 λ m (s)ds = µ m t 1 + n i=2 t i t i-1 λ m (s)ds.
(2.21)

t i t i-1 λ m (s)ds = µ m (t i -t i-1 ) + M j=1 t i t i-1 t j k <s α m j e -β m j (s-t j k ) ds = µ m (t i -t i-1 ) + M j=1 t j k ≤t i-1 α m j β m j e -β m j (t i-1 -t j k ) -e -β m j (t i -t j k ) .
(2.22)

Hence T 0 λ m (s)ds = µ m t n + M j=1 n i=2 t j k ≤t i-1 α m j β m j e -β m j (t i-1 -t j k ) -e -β m j (t i -t j k ) .
(2.23)

The sum

n i=2 t j k ≤t i-1 α m j β m j e -β m j (t i-1 -t j k ) -e -β m j (t i -t j k ) (2.24)
is telescopic and equals after simplification

t j k α m j β m j
1e -β m j (t n -t j k ) .

(2.25) Therefore,

T 0 λ m (s)ds = µ m t n + M j=1 t j k α m j β m j 1 -e -β m j (t n -t j k ) . ( 2
+ t m k ln         µ m + M j=1 P p=1 α p m j A p m j (k)         , (2.27
A p m j (k) = e -β p m j (t m k -t m k-1 ) A p m j (k -1) + t m k-1 ≤t j l <t m k e -β p m j (t m k -t j l ) , k ≥ 2, (2.30)
with the initial condition A p m j (1) = 0.

(2.31)

In a univariate setting, this simply reads A p (k) = e -β p (t k -t k-1 ) (A p (k -1) + 1) .

(2.32) Given a data set t 1 , . . . , t n and z 1 , . . . , z n , the maximum likelihood estimation method (MLE) consists in maximizing L({t i }, {z i }), or equivalently its logarithm ln L({t i }, {z i }), with respect to the model parameters µ, α, and β.

The rationale is that one maximizes the probability of observing the actual sample t 1 , . . . , t n , z 1 , . . . , z n . Formal foundations of the method for the identification of Hawkes processes (albeit in the univariate case only) can be found in [START_REF] Ogata | The asymptotic behavior of maximum likelihood estimators for stationary point processes[END_REF] (theorems 2 and 5). There, it is shown (under technical conditions) that the MLE estimator is consistent, that is tends to the actual parameters values, and asymptotically normal (around the parameters) as the sample size tends to infinity.

Separability of the log-likelihood of multivariate Hawkes processes.

We note that for each m ∈ 1, M , the partial log-likelihood ln L m (formula (2.27)), depends on µ m , (α p m j ) 1≤ j≤M,1≤p≤P and (β p m j ) 1≤ j≤M,1≤p≤P only. There is no dependence on (µ

m ′ ) m ′ m , (α p m ′ j ) m ′ m,1≤ j≤M,1≤p≤P nor (β p m ′ j ) m ′ m,1≤ j≤M,1≤p≤P
. Thus in order to maximize ln L (a function of M + 2M 2 P variables), one only needs to maximize the partial log-likelihood ln L m (a function of 1+2MP variables) separately for each m ∈ 1, M . This substantially reduces the complexity of the maximization.

Illustration

Goodness of Fit Assessment

Equally important to the identification of a point process is the assessment of the model's ability to mimic the data. We propose three approaches.

A Pragmatic Assessment.

The first approach is the easiest to implement. First estimate the parameters by MLE. Then simulate the process with the obtained parameters (details about the simulation can be found in section). Finally Q-Q-plot the inter-arrival times of the simulated sample against those of the data. If the obtained plot follows the 45 • line y = x, then the fit is satisfactory. Incidentally, this method takes into account the noise inherent to the simulation.

A Theoretically Sounder

Assessment: the Time-rescaling Theorem. We first state the underlying theorem.

Proposition 1 (Time-rescaling). Let N(t) = (N 1 (t), . . . , N M (t)) be a Mvariate Hawkes process. Fix m ∈ 1, M and define the sequence (τ m k ) 2≤k≤n m of random variables by

τ m k = t m k t m k-1 λ m (s)ds. (3.1)
Then (τ m k ) 2≤k≤n m are i.i.d. exponential random variables with parameter 1.

Proof. Can be found in [START_REF] Brown | The time-rescaling theorem and its application to neural spike train data analysis[END_REF] (univariate) or [START_REF] Bowsher | Modelling security market events in continuous time: intensity based, multivariate point process models[END_REF] (multivariate setting).

The principle of the method is to estimate the parameters by MLE, compute the quantities (τ m k ), then Q-Q-plot them against a standard exponential random variable. To compute the (τ m k ) in practice, we made use of equation (2.22)

t i t i-1 λ m (s)ds = µ m (t i -t i-1 ) + M j=1 t j k ≤t i-1 α m j β m j e -β m j (t i-1 -t j k ) -e -β m j (t i -t j k )
derived earlier for the integrated intensity, and the additivity of integration.

3.

3. An Even Sharper Assessment. If we transform the (τ m k ) to

ξ m k = 1 -e -τ m k , (3.2)
then (ξ m k ) 2≤k≤n m are independent uniform random variables on the interval ]0, 1[. It is known that the k-th order statistic ξ m (k) of a uniform random variable follows a Beta distribution with parameters k and n mk + 1

f (ξ | k, n m -k + 1) = n m ! (n m -k)!(k -1)! ξ k-1 (1 -ξ) n m -k , 0 < ξ < 1. (3.3)
In (3.3), we stress that the notation n m designates the total number of jumps in the porcess (N m (t))-the m is not for exponentiation.

We can set 95% confidence bounds on ξ m (k) using the 2.5% and 97.5% quantiles of the Beta distribution (These are available in most statistical software packages.) Alternatively, as noted in [START_REF] Brown | The time-rescaling theorem and its application to neural spike train data analysis[END_REF], for medium to large sized samples, we can use the Gaussian approximation of the Binomial distribution in equation (3.3). The 95% confidence intervals thus obtained are

ξ m (k) ± 1.96 ξ m (k) (1 -ξ m (k) ) n m . (3.4)
The idea is then, for each m ∈ 1, M , to compute (τ m k ) 2≤k≤n m as in section 3.2, transform them into (ξ m k ), order the latter into (ξ m Detailed examples of goodness of fit assessment are given in section 4. A fourth diagnosis method relying on the covariance function is thoroughly discussed in section 5.

Application to Tick By Tick Order Book Data

We refer to chapter 3, section 2 for the description of the tick by tick data we use and its processing. We estimate a 4-variate Hawkes process to the flow of:

(1) buy market orders M + (t), (2) sell market orders M -(t), (3) sell limit orders L + (t), (4) buy limit orders L -(t).

In figures 2 and 3 the resulting goodness of fit plots are displayed, and in figure 4 the estimated decay kernel.

The following practical remarks are in order:

• To achieve the fit accuracy of figures 2 and 3, we set the number of exponentials to P = 4. • Since the time resolution of the data is 1 millisecond, we add a random ∆t = 10 -3 U second, with U ∼ Uniform ]0, 1[ to points with the same time stamp (around 15% of the points for liquid CAC 40 stocks.) • We consider limit orders arriving 5 ticks away from the best opposite quote, lumping together limit orders at each side of the book in one order flow. More sophisticated models could have a point process attached to each limit.

• We illustrate with the particular stock TOTF.PA (TOTAL), but the quality of the fit is similar across all CAC 40 stocks. • We do not consider cancellations. Apart from reducing the dimensionality of the model from 6 to 4 order flows, this is motivated by the theoretical analysis of chapters 2 and 4, where it is shown that proportional cancellation rates are important for the stability of the order book. • We use warm-start optimization strategies to accelerate the estimation. We first optimize over a subset (30 min) of the data then use the result as a starting point for the optimization over the full dataset (2 hours).

Non-parametric Estimation of Multivariate Hawkes Processes

Although the results obtained by MLE are (very) satisfactory, we digress on a different estimation method recently proposed by Bacry and Muzy in [START_REF] Bacry | Hawkes model for price and trades high-frequency dynamics[END_REF]. We do this for the power of the method (when applicable) and the elegance of its underlying principle. While MLE exploits the structure of inter-event durations, the proposed method relies on the properties of counts of events, and is in a sense complementary to MLE. The idea dates back to the original papers [Haw71a, Haw71b] and is based on the covariance density matrix.

5.1. The Covariance Density Matrix. With usual notations, let N(t) = (N 1 (t), . . . , N M (t)) ⊤ be a M-variate regular 1 point process. Definition 5.1 (Covariance density matrix). The covariance density matrix of (N(t)) is the M × M matrix-valued function defined, for all x ∈ R * , by

C i j (x) = E dN i (t + x) dx dN j (t) dt -E dN i (t + x) dx E dN j (t) dt = E dN i (t + x) dx dN j (t) dt -λ i λ j , 1 ≤ i, j ≤ M, (5.1) 1
That is, the probability of multiple jumps occuring at the same time is null, or equivalently,

Prob [dN i (t) > 1] = o(dt), 1 ≤ i ≤ M.
C(x) = E dN(t + x) dx dN ⊤ (t) dt -E dN(t + x) dx E dN ⊤ (t) dt = E dN(t + x) dx dN ⊤ (t) dt -λ λ ⊤ , (5.2)
where

λ = E [λ(t)] = (λ 1 , . . . , λ M ) ⊤ , (5.3)
is the mean intensity of the point process.

In definition 5.1, dx and dt are better thought of as small (indeed, infinitesimal) time increments. The quantity dN i (t) is the number of jumps of (N i ) occurring in the interval [t, t + dt[. For small time intervals, the process being regular, this quantity is either 0 or 1.

Since

E (dN i (t)) 2 = [dN i (t)] = λ i dt,
(5.4) the diagonal terms C ii (x), 1 ≤ i ≤ M, are singular at x = 0. It is convenient to isolate the singularity, defining the complete covariance density for all x ∈ R as

C (c) (x) = Dδ(x) + C(x), (5.5)
whith D = diag(λ).

(5.6) δ(x) is Dirac's delta function, and C(x) continuous at the origin. The upper script (c) stands for "complete". Note that for all x ∈ R

C i j (-x) = C ji (x), 1 ≤ i, j ≤ M, ( 5.7)
or in more compact notation C(-x) = C ⊤ (x).

(5.8) Thus, knowing C on R + is enough to identify it.

Example 5.1 (Univariate Poisson process with constant intensity µ). Clearly in this case, for all x ∈ R,

C (c) (x) = µδ(x), C(x) = 0.
(5.9)

The Hawkes-Wiener-Hopf Integral Equation.

The following lemma will be used. (5.11)

Lemma 1. If H(t) = σ(N(s), s ≤ t)
E [λ(t)X] = E E dN(t) dt |H t X = E E dN(t) dt X|H t = E dN(t) dt X .
(5.12)

We now particularize to the case of (multivariate) Hawkes processes. We write for x > 0

C(x) = E dN(t + x) dx dN ⊤ (t) dt -λ λ ⊤ .
(5.13)

E dN(t + x) dx dN ⊤ (t) dt = E λ(t + x) dN ⊤ (t) dt = E µ + t+x -∞ ϕ(t + x -u)dN(u) dN ⊤ (t) dt = µλ ⊤ + E t+x -∞ ϕ(t + x -u)dN(u) dN ⊤ (t) dt .
(5.14)

E t+x -∞ ϕ(t + x -u)dN(u) dN ⊤ (t) dt = E x -∞ ϕ(x -v)dN(t + v) dN ⊤ (t) dt = x -∞ ϕ(x -v)E dN(t + v) dv dN ⊤ (t) dt dv = x -∞ ϕ(x -v) C (c) (v) -λ λ ⊤ dv = ϕ(x)D + x -∞ ϕ(x -v)C(v)dv - x -∞ ϕ(x -v)dv λ λ ⊤ .
(5.15)

C(x) = µ - x -∞ ϕ(x -v)dv λ -λ λ ⊤ + ϕ(x) D + x -∞ ϕ(x -v)C(v)dv. (5.16) Since µ - ∞ 0 ϕ(t)dt λ -λ = 0, (5.17) we get, x > 0, C(x) = ϕ(x) D + x -∞ ϕ(x -v)C(v)dv. (5.18)
This is an integral equation linking C and ϕ. We shall rewrite it in a classical integral equation form. Firstly

ϕ(x) = C(x)D -1 - x -∞ ϕ(x -v)C(v)D -1 dv. If we define G(x) = C(x) D -1 , (5.19) ϕ(x) = G(x) - x -∞ ϕ(x -v)G(v)dv.
(5.20)

Let the transposes of G and ϕ be

K(x) = G ⊤ (x).
(5.21)

φ(x) = ϕ ⊤ (x). (5.22) Then φ(x) = K(x) - x -∞ K(v)φ(x -v)dv.
(5.23)

A last change of variables in the integral yields the desired result.

Proposition 5.1 ([Haw71a, Haw71b]). With the above notations, it holds for all x > 0,

φ(x) = K(x) - ∞ 0 K(x -t)φ(t)dt.
(5.24)

We call this the Hawkes-Wiener-Hopf integral equation (HWH). It is a Fredholm integral equation of the second kind (actually a linear system thereof 2 ). As such, it lends itself to a variety of solution methods [START_REF] Delves | Computational Methods for Integral Equations[END_REF][START_REF] Press | Numerical Recipes, The Art of Scientific Computing[END_REF]. We choose the simplest: Nyström's.

In order to get a practical sense of the functions at play, we illustrate, before we proceed, with the following example. 2 Note that both K and φ are M by M matrices. In particular, the product inside the integral is not commutative.

Example 5.2. Consider the univariate exponential case K(x) = Ae -B|x| .

(5.25)

φ(x) = K(x) - ∞ 0 K(x -t)φ(t)dt.
(5.26)

With slightly more general notation

φ(x) = F(x) + κ ∞ 0 K(x -t)φ(t)dt. (5.27) κ = -1, F(x) = K(x).
(5.28)

This is close to the Lalesco-Picard integral equation. We rewrite

φ(x) = F(x) + κAe -Bx x 0 e Bt φ(t)dt + κAe Bx ∞ x e -Bt φ(t)dt.
(5.29)

Differentiation twice with respect to x yields (5.37)

φ ′ (x) = F ′ (x) -κABe -Bx x 0 e Bt φ(t)dt + κAe Bx ∞ x e -Bt φ(t)dt. (5.30) φ ′′ (x) = F ′′ (x) -B 2 F(x) + B(-2κA + B)φ(x). ( 5
We can therefore state.

Proposition 5.2. Suppose the covariance density of a univariate Hawkes process writes

C(x) = λ Ae -B|x| , x ∈ R, A ≥ 0, B > 0.
(5.38)

λ being the empirical mean intensity. Then, the process's kernel is

ϕ(x) = αe -βx , x ≥ 0, (5.39) ϕ(x) = 0, x < 0, (5.40) with β = √ B 2 + 2AB (5.41) α = A(B + √ B 2 + 2AB) A + B + √ B 2 + 2AB = A(B + β) A + B + β .
(5.42)

These are the inverse of equation ( 16) in [Haw71b]. When appropriate, that is, when the empirical covariance is well fit by one exponential 3 , formulae (5.39)-(5.42) make the identification of a univariate Hawkes process 4 immediate. As an aside, they could serve as a benchmark for numerical solutions of equation (5.24). 5.3. Nyström's Method. Consider the system of integral equations (5.24), x > 0,

φ(x) = F(x) + κ ∞ 0 K(x -t)φ(t)dt.
(5.48) 3 This is rarely the case in financial datasets. Typical covariances are power-law tailed.

β 2 + Bβ -AB ≥ β 2 + B 2 -AB = 2B 2 + AB > 0.
(5.47) φ is the unknown, a M by M matrix-valued function. Both K and F are given M by M matrix-valued functions. φ is defined on R + (By causality of the kernel, it is null outside.) F is defined on R + and K on R. κ is a constant.

Let {t m } 0≤m≤N and {w m } 0≤m≤N be the nodes and weights of a certain quadrature rule on [0, ∞[. For all 0 ≤ m ≤ N,

φ(t m ) = F(t m ) + κ ∞ 0 K(t m -t)φ(t)dt.
(5.49)

Discretize the integral,

φ(t m ) = F(t m ) + κ N n=0 w n K(t m -t n )φ(t n ). (5.50) Rearrange, N n=1 [Iδ mn -κw n K(t m -t n )] φ(t n ) = F(t m ), 0 ≤ m ≤ N. (5.51) Or, in matrix notation                   I -κw 0 K(0) -κw 1 K(t 0 -t 1 ) • • • -κw N K(t 0 -t N ) -κw 0 K(t 1 -t 0 ) I -κw 1 K(0) • • • -κw N K(t 1 -t N ) . . . . . . . . . . . . -κw 0 K(0) -κw 1 K(t N -t 1 ) • • • I -κw N K(0)                                     φ(t 0 ) φ(t 1 ) . . . φ(t N )                   =                   F(t 0 ) F(t 1 ) . . . F(t N )                   . (5.52) This is a M(N + 1) by M(N + 1) matrix equation in                   φ(t 0 ) φ(t 1 ) . . . φ(t N )                   ,
and can be readily solved by standard techniques.

In particular, using the trapezoïdal rule with step size h > 0 as the quadrature rule, and the original notations of the HWH equation (5.24)

F(x) = K(x),
(5.53)

κ = -1, (5.54) the linear system reads                   I + h 2 K(0) hK(t 0 -t 1 ) • • • h 2 K(t 0 -t N ) h 2 K(t 1 -t 0 ) I + hK(0) • • • h 2 K(t 1 -t N ) . . . . . . . . . . . . h 2 K(t N -t 0 ) hK(t N -t 1 ) • • • I + h 2 K(0)                                     φ(t 0 ) φ(t 1 ) . . . φ(t N )                   =                   K(t 0 ) K(t 1 ) . . . K(t N )                   . (5.55)
Note that, to evaluate K for negative t, we can use property (5.8)

K(-x) = K ⊤ (x).
The value of t N , which determines the domain of integration, should be chosen large enough for the approximation of the integral to be precise. In practise, we choose it such that the integral of the empirical density (or its rescaled transpose K) reaches a plateau before t N .

Because of the unbounded integration domain, and the fact that we are dealing with a system of integral equations, existence and uniqueness of solutions to the HWH equation are a priori not guaranteed. Nor are the stabiliy and convergence of the proposed numerical method, obvious. In practice, the method is successful on a wide range of datasets, but fails for some "difficult" cases (more on this below.) A systematic study of the properties of equation (5.24) and its numerical solution is therefore warranted 5 . These considerations are, however, beyond the scope of this thesis.

Empirical Covariance Density.

Estimating the empirical covariance density matrix is quite straightforward. One chooses a small sampling interval ∆ > 0, then sets C i j (x) = 1 ∆ 2 Cov N i (t + x + ∆) -N i (t + x), N j (t + ∆) -N j (t) , (5.56)

Cov being the sample's empirical covariance. Our experimentation with market data shows that this estimate is stable, and for reasonably small ∆, independent of the sampling interval. We note that the Fast Fourier Transform [START_REF] Press | Numerical Recipes, The Art of Scientific Computing[END_REF] can be used to speed up the evaluation of equations (5.56). to find µ.

5.6. Principle of the Non-parametric Kernel Estimation Method.

(1) Estimate the empirical mean intensity

λ i = # Jumps in N i in [0, T ] T , i ∈ 1, M .
(5.58)

(2) Estimate the empirical covariance density C. (5.59) 5.7. Illustration. Figure 6.

Simulation via Markovian Projection

In this section, a (new) technique for the fast simulation of multivariate Hawkes processes with arbitrary decay kernels is presented.

6.1. Ogata's Thinning Simulation Method. The standard method for the simulation of multivariate Hawkes processes is Ogata's modified thinning algorithm [START_REF]On Lewis' simulation method for point processes[END_REF] (See also Lewis and Schedler [LS79].) The idea is to simulate a non-homogeneous univariate Poisson process {t * n } whose intensity λ * (t) dominates the total intensity If the mark j 0 equals M + 1, the point is rejected. By proposition 1 in [START_REF]On Lewis' simulation method for point processes[END_REF], the remaining points form a point-process that has the desired intensity λ(t). For Hawkes processes, since the kernel ϕ(t) is decreasing (and therefore λ sum (t) after a jump), a natural choice for the dominant intensity λ * (t) is the value of λ sum immediately after a jump. One therefore gets a piecewise constant λ * between jumps (figure 8).

λ sum (t) = M i=1 λ i (t
We implement this idea in algorithm 2. Note that we heavily exploit the Markovian structure of the point process (chapter 4, proposition 2.2). We comment on this below.

Recall that α = α p i j 1≤i, j≤M,1≤p≤P , (6.4) Y jumps by 1 with each jump of the process. More precisely, if z n ∈ 1, M is the mark of the jump,

Y p iz n (t n+1 ) = Y p iz n (t n+1 ) + 1, 1 ≤ i ≤ M, 1 ≤ p ≤ P. (6.7)
Only a "slice" of Y jumps by 1, it corresponds the plan j = z n (figure 9).

Betwen t n and t n+1 , Y decays as 6 Y(t n+1 ) = e -β(t n+1 -t n ) Y(t n+1 ). (6.8)

These are the key recursive relationships we use in the simulation.

The exponential and multiplication are meant element by element. 6.2. General Kernel. Suppose now we are given a general decay kernel ϕ, a power-law, ϕ i j (t) = a i j (b i j + t) c i j , a i j , b i j , c i j > 0, (6.9) say. The problem at hand is not Markovian. But it can be very well approximated by one. The family of functions αe -βt (6.10) is rich enough to virtually approximate any kernel of interest (on a finite but arbitrarily large interval [0, T ] and with arbitrary accuracy). The idea is therefore to project the kernel ϕ on the "basis functions", αe -βt (6.11)

for instance by least squares fitting, then use algorithm 2 to simulate the process.

Algorithm 2 Multivariate Hawkes process simulation. Require: N, M, P, µ, α, β 1: n ← 1 2: Y ← 0 3: λ * ← M i=1 µ i 4: while n ≤ N, do Y ← e -βτ Y {eq. 6.8} 8:

for i = 2 to M do 9:

λ i = µ i + 1≤ j≤M,1≤p≤P
α p i j Y p i j 10:

end for 11:

λ sum = 1≤i≤M λ i 12:
Choose randomly the event type j 0 in 1, M + 1 according to the probability vector (λ 1 , λ 2 , . . . , λ M , λ *λ sum ) /λ * 13:

if j 0 ≤ M then 14:

{Keep the point} 15:

t n ← newT 16:

z n ← j 0 17:

Y(:, j 0 , :) ← Y(:, j 0 , :) + 1 {This is a shorthand notation of eq. (6.7).} end if 24: end while Interestingly, in the case of power-law decaying kernels, this has a nice financial interpretation. In econophysics, power-law behaviour is known to emerge from the aggregation of heterogeneous agents with different time α p e -β p t , (6.12) represents an agent (or group of agents) with a characteristic time of 1 β p . In a sense, Markovian projection identifies the most important groups. Remarkably, at least from a computing perspective, you only need about log T exponentials to approximate a power-law over the range [0, T ]. This is the key insight behind the algorithm.

Note that although we restricted the parameters to be positive, we can imagine them taking complex values (for fancier shaped kernels), exploit the Markovian structure to simulate, then return to the real line.

The speedup, from O(N 2 ) to O(N log N), is dramatic (figure 12).

Applications

In the previous sections, it was shown how to estimate and efficiently simulate a stochastic order book model based on Hawkes processes. We now turn to the measurement of two quantities of practical importance using the estimated model. The flexibility of the order book model allows the assessment of a wide range of algorithmic trading strategies, and this section is a proof-of-concept for the possibilities such models may offer. Note that we account for the effect of the agent's trading on the order flow (underlined terms). This would be impossible when merely replaying market data. The agent is submitting sell limit orders at a constant Poissonian rate, and affecting the order flow similarly to any other sell limit order. The result of the simulation is displayed in figure 13.

7.2. Execution Shortfall of a Limit Order. The second quantity is related to the cost of executing a limit order [START_REF] Nevmyvaka | Electronic trading in order-driven markets: efficient execution[END_REF]. Consider the socalled "submit and leave" strategy: A agent submits a limit order at price P with a certain volume v, waits T seconds for it to be (partially) executed, and finally submits a market order with the remaining shares. A question of interest is the determination of the optimal price and waiting time for such a simple strategy.

If we define the (algebraic) shortfall of this strategy (started at t 0 ) as 7 S = Execution Price -Prevailing mid-price at t 0 , (7.3) expressed in ticks. Then estimates of the average shortfall S and its variance are valuable indicators for the optimal choice. We run this strategy on simulated data and display the result in figure 14. For instance, it is clear, for figure 14's parameter set, that it is not optimal to place the order deep in the book. It is also possible to normalize by the mid-price S = Execution Price -Prevailing mid-price at t 0 Prevailing mid-price at t 0 , (7.2) to get the "return" of the strategy. Figure 13. Probability of fill of a sell limit order of size v = 50 shares. In panel (B), the legend indicates the price in ticks relative to the best offer. For instance, 0 means that the limit order is placed with its price set to the best ask, 1 at ask+1 tick etc. Model parameters are estimated for the stock TOTF.PA, June 2011.

Figure 1 .

 1 Figure 1. Comparative evolution of trading activity and computing power [Eco11]. y-axis is in logarithmic scale. Data sources: NYSE-Euronext / Wikipedia.

  Figure1. Order book schematic illustration: a buy market order arrives and removes liquidity from the ask side, then sell limit orders are submitted and liquidity is restored.

Figure 2 .

 2 Figure2. Order book dynamics: in this example, K = 9, q = 1,

  for a buyer initiated trade, i.e. a buy market order, (1.10) and, α(t) = -1, for a seller initiated trade, i.e. a sell market order.(1.11)

Figure 1 .

 1 Figure 1. Average depth profile. Simulation parameters are summarized in tables 3 and 4.

Figure 2 .

 2 Figure 2. Probability distribution of the spread. Note that the model (dark gray) predicts a tighter spread than the data.

1

  Figure 3. Autocorrelation of price increments. This figure shows the fast decay of the autocorrelation function, and the large negative autocorrelation of trades at the first lag.

  But

Figure 4 .

 4 Figure 4. Price sample path. At large time scales, the price process is a Brownian motion.

  Figure 5. Q-Q plot of mid-price increments. h is the time lag in seconds. This figure illustrates the aggregational normality of price increments.

FrequencyFigure 6 .

 6 Figure 6. Probability distribution of price increments. Time lag h = 1000 events.

Figure 7 .

 7 Figure7. Signature plot: σ 2 h := V [P(t + h) -P(t)]/h. y axis unit is tick 2 per trade for panel (a) and tick 2 .second -1 for panel (b). We used a 1,000,000 event simulation run for the model signature plots. Data signature plots are computed separately for each trading day [9 : 30-14 : 00] then averaged across 23 days. For calendar time signature plots, prices are sampled every second using the last tick rule. The inset is a zoom-in.

Figure 8 .

 8 Figure 8. Model parameters: arrival rates and average depth profile (parameters as in table 4). Error bars indicate variability across different trading days.

Figure 9 .

 9 Figure 9. Model parameters: volume distribution. Panels (a), (b) and (c) correspond respectively to market, limit and cancellation orders volumes. Dashed lines are lognormal fits (parameters as in table3).

Figure 10

 10 Figure 10. A cross-sectional comparison of liquidity and price diffusion characteristics between the model and data for CAC 40 stocks (March 2011).

CHAPTER 4

 4 On the Stability and Price Scaling Limit of a Hawkes Process-Based Order Book Model 1. Introduction Since their introduction in [Haw71b], Hawkes processes have been applied in a wide range of research areas from seismology (originally), to credit risk [EGG10], financial contagion [ASCDL13] and more recently market microstructure modelling [BDHM13a, BDHM13b, BM13, Bow07, Lar07, MT11b, MTP12, ZRA13].

1A

  more general definition would have

Figure 1 .

 1 Figure1. Order book dynamics: in this example, K = 9, q = 1,

)

  and the incomplete Gamma function is defined for all p ∈ C, R[p] > 0 and all z ∈ C asΓ(p, z) = ∞ z t p-1 e -t dt.(5.4)Inparticular, for all b > 0, c > 0 lim a→∞ I(a, b, c) = 0, (5.5) and lim a→∞ a I(a, b, c) = 1 b .

  and c > 0. Let b min = min{b i } 1≤i≤p . (5.8) Then I p (a 1 , . . . , a p ; b 1 , . . . , b p ; c) ≤ ∞ 0 e -p i=1 a i (1-e -b min t )-ct = I(|a|, b min , c),

  27) is the backbone of maximum likelihood estimation of linear Hawkes processes with exponential kernels.Remark 2.1. (Recursive expression of A(k) [MT11a, Oga81]) In order to accelerate the computation of the likelihood, it is custom to use the following recursive relanship:

  Figure 1.

Figure 1 .

 1 Figure 1. Negative of the partial log-likelihood,ln L 1 , for a 10000 points simulated sample as a function of (α 11 , α 12 ). The parameters values are M = 2, P = 1, µ = (0.1, 0.2) ⊤ , α = 1.0 2.0 3.0 4.0 , β = 10.0 20.0 30.0 40.0 . There is clearly a minimum at α * 11 , α * 12 = (1.0, 2.0).

  (k) ), and plot the result against the CDF of the uniform distribution b k = k-1 2 n m , with (local) confidence bounds defined by equation (3.4).

Figure 2 .

 2 Figure 2. Joint maximum likelihood estimation of 4 order flows (TOTF.PA, June 2011.) Continued in figure 3.

Figure 3 .

 3 Figure 3. Joint maximum likelihood estimation of 4 order flows (TOTF.PA, June 2011.) Continuation of figure 2.

Figure 4 .

 4 Figure 4. Maximum likelihood estimated kernel (TOTF.PA, June 2011.) This 4 × 4 matrix represents the decay kernels {ϕ i j (t)} 1≤i, j≤4 . x axis is time in seconds.

. 31 )

 31 With κ = -1 and F(x) = K(x) = Ae -Bx , x > 0 φ ′′ (x) = B(2A + B)φ(x). + β) A + B + β .

4

  It is easy to check that the stability condition α < β,(5.43) is always met for the Hawkes process thus defined. Ideed A(B + β) A + B + β <

5. 5 .Functional

 5 Estimating the Base Rate µ. Once ϕ has been estimated as the (numerical) solution to equation (5.24), we make use of the identity µ = Itransforms might be the right path to this end [Haw71a].

( 3 )

 3 Rescale it to K = C Diag( λ) -1 . (4) Solve the HWH equation for φ = ϕ ⊤ . (5) Set µ = I -∞ 0 ϕ(t)dt λ.

Figure 5 ..

 5 Figure 5. Non-parametric kernel estimation. Artificial data, 10000 points sample, is generated with the parameters M = 2, P = 1, µ = (0.1, 0.2) ⊤ , α = 5.0 10.0 1.0 2.0 , β = 20.0 15.0 3.0 10.0 . The covariance density matrix is then estimated, and the Hawkes-Wiener-Hopf equation solved by Nystöm's method. The number of integration nodes is N = 3000 and t N = 60 seconds. The estimated base intensity is µ = (0.0967, 0.1932) ⊤ .

Figure 6 .

 6 Figure 6. Non-parametric kernel estimation (Buy and sell market orders, TOTF.PA, June 2011.)

Figure 7 .Figure 8 .

 78 Figure 7. Data versus model covariance density (Buy and sell market orders, TOTF.PA, June 2011.)

Figure 9 .

 9 Figure 9. Illustration for equation (6.7).

  inter-event time according to an exponential r.v. with mean λ * τ ∼ Exp(1/λ * ) 6: newT ← newT + τ 7:

  1678 2.0049 3.1901 ϕ 12 1.1014 2.0511 3.1021 ϕ 21 1.1591 2.0395 3.1803 ϕ 22 1.0822 2.0690 3.0730 Table 1. Power-law kernel parameters. Each triplet (a, b, c) was randomly generated by (1, 2, 3) + 0.2(U 1 , U 2 , U 3 ), with U i ∼ Uniform]0, 1[. horizons. Each term in the exponential sum a (b + t) c ≈ P p=1

7. 1 .

 1 Probability of Fill of a Limit Order. The first quantity of interest is the probability of fill of a limit order. Its estimation from market data can be challenging (See e.g.[START_REF] Lo | Econometric models of limit-order executions[END_REF]) and we propose an indirect (and The first 5 order flows form a Hawkes process whose decay kernel matrix readsϕ =                       ϕ 11 ϕ 12 ϕ 13 ϕ 14 ϕ 13 ϕ 21 ϕ 22 ϕ 23 ϕ 24 ϕ 23 ϕ 31 ϕ 32 ϕ 33 ϕ 34 ϕ 33 ϕ 41 ϕ 42 ϕ 43 ϕ 44 ϕ 43

  7

Figure 14 .

 14 Figure 14. Execution shortfall (A) and its standard deviation (B) for a limit order of size v = 1000 shares. Model paramters are estimated for the stock TOTF.PA, June 2011.

Table 1 .

 1 Tick by tick data file sample. Note that the field "Level" does not necessarily correspond to the distance in ticks from the best opposite quote as there might be gaps in the book. Lines corresponding to the trades in table 2 are highlighted in italics.

	Timestamp Side Level Price Quantity
	33480.158 B	1	121.1	480
	33480.476 B	2	121.05	1636
	33481.517 B	5	120.9	1318
	33483.218	B	1	121.1	420
	33484.254 B	1	121.1	556
	33486.832 A	1	121.15	187
	33489.014 B	2	121.05	1397
	33490.473	B	1	121.1	342
	33490.473	B	1	121.1	304
	33490.473	B	1	121.1	256
	33490.473 A	1	121.15	237
	Timestamp Last Last quantity
	33483.097 121.1	60
	33490.380 121.1	214
	33490.380 121.1	38
	33490.380 121.1	48

Table 2 .

 2 Trades data file sample.

  ).

	b 1	b 2	R 2
	Log A (5) -0.42 (±0.11) 1.13 (±0.04) 0.99 Log S 0.20 (±0.06) 1.16 (±0.07) 0.97
	σ ∞ Table 5. CAC 40 stocks regression results. -0.012 (±0.05) 1.35 (±0.11) 0.94

  Ergodicity of the Intensity of a Multivariate Hawkes Process. Consider now a multivariate setting. Let X(t) = (N(t), Y(t)) be a M-variate Hawkes process with parameters

	µ = (µ 1 , . . . , µ M ) t ,	(3.14)
	α = (α i j ) 1≤i, j≤M ,	(3.15)
	and
	β = (β i j ) 1≤i, j≤M .	(3.16)

  11 +α 12 ) ω+γ(α 11 +α 12 ) µ 1 + e ω+γ(α 21 +α 22 ) µ 2 + e -ω λ 3 x)

	As |x| + |y| → ∞, PV(x, y) V(x, y) ≤ e -ω +	2 β min (γ + 1 β max )	(e ω+γ(α 11 +α 21 ) + e -ω+γ(α 12 +α 22 ) ). (3.44)
					∞	e	-1≤k,l≤2 (γ+ 1 β kl
					0
					∞	e -1≤k,l≤2 (γ+ 1 β kl
					0
	+ e -ω	∞	e -1≤k,l≤2 (γ+ 1 β kl
				0
	×	0	∞	e -1≤k,l≤2 (γ+ 1 β kl	)(1-e -βklt )y kl -(µ 1 +µ 2 )t-λ 3 xt dt
	+ e ω+γ(α 11 +α 21 ) y 11	∞	e	-1≤k,l≤2 (γ+ 1 β kl	)(1-e -βklt )y kl -(µ 1 +µ 2 )t-λ 3 xt-β 11 t dt
						0
	+ e ω+γ(α 11 +α 21 ) y 12	∞	e -1≤k,l≤2 (γ+ 1 β kl	)(1-e -βklt )y kl -(µ 1 +µ 2 )t-λ 3 xt-β 12 t dt
						0
	+ e -ω+γ(α 12 +α 22 ) y 21	∞	e -1≤k,l≤2 (γ+ 1 β kl	)(1-e -βklt )y kl -(µ 1 +µ 2 )t-λ 3 xt-β 21 t dt
						0
	+ e -ω+γ(α 12 +α 22 ) y 22	∞	e -1≤k,l≤2 (γ+ 1 β kl	(3.43)
						0

)(1-e -βklt )y kl -(µ 1 +µ 2 )t-λ 3 xt (µ 1 + y 11 e -β 11 t + y 12 e -β 12 t )dt

+ e -ω+γ(α 21 +α 22 ) )(1-e -βklt )y kl -(µ 1 +µ 2 )t-λ 3 xt (µ 2 + y 21 e -β 21 t + y 22 e -β 22 t )dt )(1-e -βklt )y kl -(µ 1 +µ 2 )t-λ 3 xt λ 3 x dt = (e

)(1-e -βklt )y kl -(µ 1 +µ 2 )t-λ 3 xt-β 22 t dt.

  Likelihood of Multivariate Hawkes Processes. Consider now the case of a M-variate Hawkes process, M ≥ 1. The intensity vector λ = (λ 1 , . . . , λ M ) ⊤ is specified by

							.11)
	Finally
	ln L(t 1 , . . . , t n ) = -µt n +	n k=1	α β	e -β(t n -t k ) -1 +	n i=1	ln (µ + αA(i)) . (2.12)
	with
	A(i) =	t k <t i	e -β(t i -t k ) , for i ≥ 2,	(2.13)
	and
				A(1) = 0.	(2.14)
	More generally for P ≥ 1, it holds,
	Proposition 2.1 (Log-likelihood of a univariate Hawkes process). ln L(t 1 , . . . , t n ) = -µt n + n k=1 P p=1 α p i=1 p=1 β p e -β p (t n -t k ) -1 + n ln         µ + P	α p A p (i)	        ,
							(2.15)
	where for all p ∈ 1, P , A p (i) =	t k <t i	e -β p (t i -t k ) , i ≥ 2,	(2.16)
	and
			A p (1) = 0.	(2.17)
	2.2.

  ({t i , z i } 1≤i≤n ) = -µ m t n +

						.26)
	More generally for P ≥ 1, we state,
	Proposition 2.2 (Log-likelihood of a multivariate Hawkes process).
	ln L m M j=1 t	j k	P p=1	α p m j β p m j	e -β p m j (t n -t	j k ) -1

MARKOVIAN ORDER BOOK MODELS I

Actually, taking for a ∞ and |b ∞ | independent positive random variables would not change much our analysis. We take constants for simplicity.

The results are qualitatively similar for all CAC 40 stocks.

NUMERICAL EXAMPLE

Thomson Reuters Tick History.

ϕ(u) = αe -βu .

Cf. Dynkyn's formula or its discrete-time formulation.

Remerciements

original) approach for its estimation. To do this, we simulate an order book model with 7 order flows:

(1) M + (t): buy market orders, (2) M -(t): sell market orders, (3) L + (t): sell limit orders, (4) L -(t): buy limit orders, (5) L + A (t): An additional agent's sell limit orders, (6) C + (t): cancellation of a sell limit order, (7) C -(t): cancellation of a buy limit order.