Frédéric Abergel
email: frederic.abergel@ecp.fr

Côme Huré
email: hure@lpsm.paris

Huyên Pham
email: pham@lpsm.paris

Prof Frédéric Abergel

Prof Huyên Pham

Algorithmic trading in a microstructural limit order book model

Keywords: Limit order book, pure-jump controlled process, high-frequency trading, high-dimensional stochastic control, Markov Decision Process, quantization, local regression

Introduction

Most of the markets use a limit order book (order book) mechanism to facilitate trade. Any market participant can interact with the order book by sending either market orders or limit orders. In such type of markets, the market makers play a fundamental role by providing liquidity to other market participants, typically to impatient agents who are willing to cross the bid-ask spread. The profit made by a market-making strategy comes from the alternation of buy and sell orders.

From the mathematical modeling point of view, the market-making problem corresponds to the choice of an optimal strategy for the placement of orders in the order book. Such a strategy should maximize the expected utility function of the wealth of the market maker up to a penalization of their inventory. In the recent litterature, several works focused on the problem of market-making through stochastic control methods. The seminal paper by [START_REF] Avellaneda | High-frequency trading in a limit order book[END_REF] inspired by the work of [START_REF] Ho | Optimal dealer pricing under transactions and return uncertainty[END_REF] proposes a framework for trading in an order driven market. They modeled a reference price for the stock as a Wiener process, and the arrival of a buy or sell liquidity-consuming order at a distance δ from the reference price is described by a point process with an intensity in an exponential form decreasing with δ. They characterized the optimal market-making strategies that maximize an exponential utility function of terminal wealth. Since this paper, other authors have worked on related market-making problems. [START_REF] Guéant | Dealing with the Inventory Risk[END_REF] generalized the market-making problem of [START_REF] Avellaneda | High-frequency trading in a limit order book[END_REF] by dealing with the inventory risk. [START_REF] Cartea | Risk metrics and fine tuning of high-frequency trading strategies[END_REF] also designed algorithms that manage inventory risk. Fodra and Pham (2015b) and Fodra and Pham (2015a) considered a model designed to be a good compromise between accuracy and tractability, where the stock price is driven by a Markov Renewal Process, and solved the market-making problem. [START_REF] Guilbaud | Optimal High Frequency Trading with limit and market orders[END_REF] also considered a model for the mid-price, modeled the spread as a discrete Markov chain that jumps according to a stochastic clock, and studied the performance of the market-making strategy both theoretically and numerically. [START_REF] Cartea | Modeling Asset Prices for Algorithmic and High Frequency Trading[END_REF] employed a hidden Markov model to examine the intra-day changes of dynamics of the order book. Very recently, [START_REF] Cartea | Algorithmic and High-frequency trading[END_REF] and [START_REF] Guéant | The Financial Mathematics of Market Liquidity: From optimal execution to market making[END_REF] published monographs in which they developped models for algorithmic trading in different contexts. El Aoud and [START_REF] El Aoud | A stochastic control approach for options market making[END_REF] extended the framework of Avellaneda and Stoikov to the options market-making. A common feature of all these works is that a model for the price or/and the spread is considered, and the order book is then built from these quantities. This approach leads to models that predict well the long-term behavior of the order book. The reason for this choice is that it is generally easier to solve the market-making problem when the controlled process is lowdimensional. Yet, some recent works have introduced accurate and sophisticated micro-structural order book models. These models reproduce accurately the short-term behavior of the market data. The focus is on conditional probabilities of events, given the state of the order book and the positions of the market maker. [START_REF] Abergel | Limit Order Books[END_REF] proposed models of order book where the arrivals of orders in the order book are driven by Poisson processes or Hawkes processes. [START_REF] Cont | A stochastic model for order book dynamics[END_REF] also modeled the orders arrivals with Poisson processes. [START_REF] Huang | Simulating and analyzing order book data: The queue-reactive model[END_REF] proposed a queue-reactive model for the order book. In this model the arrivals of orders are driven by Cox point processes with intensities that only depend on the state of the order book (they are not time dependent). Other tractable dynamic models of order-driven market are available (see e.g. [START_REF] Cont | A stochastic model for order book dynamics[END_REF], [START_REF] Rosu | A dynamic model of the limit order book[END_REF], [START_REF] Cartea | Buy Low Sell High: a High Frequency Trading Perspective[END_REF]).

In this paper we adopt the micro-structural model of order book in [START_REF] Abergel | Limit Order Books[END_REF], and solve the associated trading problem. The problem is formulated in the general framework of Piecewise Deterministic Markov Decision Process (PDMDP), see [START_REF] Bäuerle | Markov Decision Processes with Applications to Finance[END_REF]. Given the model of order book, the PDMDP formulation is natural. Indeed, between two jumps, the order book remains constant, so one can see the modeled order book as a point process where the time becomes a component of the state space. As for the control, the market maker fixes their strategy as a deterministic function of the time right after each jump time. We prove that the value function of the market-making problem is equal to the value function of an associated non-finite horizon Markov decision process (MDP). This provides a characterization of the value function in terms of a fixed point dynamic programming equation. [START_REF] Jacquier | Optimal Liquidation in a level-I Limit order book for large tick stocks[END_REF] recently followed a similar idea to solve an optimal liquidation problem, while [START_REF] Baradel | Optimal inventory management and order book modeling[END_REF] and [START_REF] Lehalle | Optimal liquidity-based trading tactics[END_REF] also tackled this problem of reward functional maximization in a micro-structure model of order book framework.

The second part of the paper deals with the numerical simulation of the value functions. The computation is challenging because the micro-structural model used to model the order book leads to a high-dimensional pure jump controlled process, so evaluating the value function is computationally intensive. We rely on control randomization and Markovian quantization methods to compute the value functions. Markovian quantization has been proved to be very efficient for solving control problems associated with high-dimensional Markov processes. We first quantize the jump times and then quantize the state space of the order book. See [START_REF] Pagès | Optimal quantization methods and applications to numerical problems in finance[END_REF] for a general description of quantization applied to controlled processes. The projections are time-consuming in the algorithm, but Fast approximate nearest neighbors algorithms (see e.g. [START_REF] Muja | Fast Approximate Nearest Neighbors with Automatic Algorithm Configuration[END_REF]) can be implemented to alleviate the procedure. We borrow the values of intensities of the arrivals of orders for the order book simulations from [START_REF] Huang | Simulating and analyzing order book data: The queue-reactive model[END_REF] in order to test our optimal trading strategies.

The paper is organized as follows. The model setup is introduced in Section 2: we present the micro-structural model for the order book, and show how the market maker interacts with the market. In Section 3, we prove the existence and provide a characterization of the value function and optimal trading strategies. In Section 4, we introduce a quantization-based algorithm to numerically solve a general class of discrete-time control problem with finite horizon, and then apply it on our trading problem. We then present some results of numerical tests on simulated order book. Section 5 presents an extension of our model when order arrivals are driven by Hawkes processes, and finally the appendix collects some results used in the paper.

Model setup

Order book representation

We consider a model of the order book inspired by the one introduced in chapter 6 of [START_REF] Abergel | Limit Order Books[END_REF].

Let us fix K ≥ 0. An order book is supposed to be fully described by K limits 1 on the bid side and K limits on the ask side. Denote by pa t the best ask at time t, which is the cheapest price a participant in the market is willing to sell a stock at time t, and by pb t the best bid at time t, which is the highest price a participant in the market is willing to buy a stock at time t. We use the pair of vectors a t , b t = a 1 t , ..., a K t , b 1 t , . . . , b K t where

• a i t is the number of shares available i ticks away from pb t , • -b i t is the number of shares available i ticks away from pa t , to describe the order book. The vectors a t and b t describe respectively the ask and the bid sides at time t. The quantities a i t , 1 ≤ i ≤ K, live in the discrete space qN where q ∈ R * is the minimum order size on each specific market (lot size). The quantities b i t , 1 ≤ i ≤ K, live in the discrete space -qN. By convention, the a i are non-negative, and the b i are non-positive for 0 ≤ i ≤ K. The tick size represents the smallest intervall between different price levels. We assume in the sequel that the orders arrivals have the same size q = 1, and set the tick size to = 1 for simplicity.

Constant boundary conditions are imposed outside the moving frame of size 2K in order to guarantee that both sides of the LOB are never empty: we assume that all the limits up to the K-th ones are equal to a ∞ in the ask side, and equal to b ∞ in the bid side, with a ∞ , -b ∞ ∈ N. The order book can receive at any time three different kinds of orders from general market participants: market orders, limit orders and cancel orders. The orders arrivals are modeled by the following point processes:

• M + stands for the buy market orders flow, and we denote by λ M + its intensity, • M -stands for the sell market orders flow, and we denote by λ M -its intensity, • L + i , for i ∈ {1, ...K}, stands for the sell orders flow at the i th limit on the ask side, and we denote by λ L + i its intensity, 1 limit is also referred to as quote or price level in the literature.

• L - i , for i ∈ {1, ...K}, stands for the buy orders flow at the i th limit on the bid side, and we denote by λ L - i its intensity, • C + i , for i ∈ {1, ...K}, stands for the cancel orders flow at the i th limit on the ask side, and we denote by λ C + i its intensity, • C - i , for i ∈ {1, ...K}, stands for the cancel orders flow at the i th limit on the bid side, and we denote by λ C - i its intensity.

We assume in the sequel that (Harrivals) The orders arrivals from general market participants (market orders, limit orders and cancel orders) occur according to Markov jump processes which intensities only depends on the couple a, b . Moreover, we assume that the all the intensities are at most linear w.r.t. the couple a, b and are constant between two events.

Under (Harrivals), Let λ L , λ C , λ M be positive real constants such that Remark 2.1 The linear conditions on the intensities are required to prove that the control problem is well-posed.

K i=1 λ L + i (a, b) + K i=1 λ L - i (a, b) ≤ λ L |a| + |b| , K i=1 λ C + i (a, b) + K i=1 λ C - i (a, b) ≤ λ C |a| + |b| , λ M + (a, b) + λ M -(a, b) ≤ λ M |a| + |b| ,
Remark 2.2 We assume the intensity to be constant between jumps in (Harrivals) for simplicity.

All the results proposed in Section 3 can be extended to the case where the intensities of the jump processes are deterministic between jumps. Such an extension is considered in Section 5, where the arrivals are modeled using Hawkes processes with exponential kernel.

Remark 2.3 Some information can be integrated to the order book model by adding new processes. For example, some exogenous processes that send orders to the best-ask and best-bid limits can be added to model the predictions of the mid-price and its volatility that an agent may have. Doing so is critical to manage the risk-reward tradeoffs.

Market maker strategies

We assume that the order book matching is done on price/time priority, which means that each limit of the order book is a queue where the first order in the queue is the first one to be executed1 .

We consider a market maker who stands ready to send buy and sell limit orders on a regular and continuous basis at quoted prices. A usual assumption in stochastic control theory to characterize the value function as solution to a HJB equation is to constrain the control space to be compact. In this spirit, we shall make the following assumption:

(Hcontrol) Assume that at any time, the total number of limit orders placed by the marker maker does not exceed a fixed (possibly large) integer M .

2.2.1.

Control of the market maker. The market maker can choose at any time to keep, cancel or take positions in the order book (as long as she does not hold more than M positions in the order book). Her positions are fully described by the following M -dimensional vectors ra t , rb t , na t , nb t where ra (resp. rb) records the limits in which the market maker's sell (resp. buy) orders are located; and na (resp. nb) records the ranks in the queues of each market maker's sell (resp. buy) orders. In order to guarantee that the strategy of the market maker is predictable w.r.t. the natural filtration generated by the orders arrivals processes, we shall make the following assumption.

(Harrivals2) The intensities do not depend on the control. Moreover, the market maker does not cross the spread.

To simplify the theoritical analysis, we also make the following assumption: (Harrivals3) Assume that the market maker does not change their strategy between two orders arrivals of the order book. In other words, the market maker makes a decision right after one of the order arrivals processes L ± , C ± , M ± jumps, and keep it until the next the jump of an order arrival.

Remark 2.4 Assumption (Harrivals3) is mild if the order book jumps frequently, since the market maker can change their decisions frequently in such a case. It can also easily be relaxed by considering piecewise constant controls between jumps (which seems well-adapted to most of the time-discretized control problems met in the industry) or any other parametric family of functions. The results and proofs can then be extended, relying mainly on the PDMDP1 .

We provide in Figure 1 a graphical representation of the controlled LOB. Notice that the market maker interacts with the order book by placing orders at some limits. The latter have ranks that evolve after each orders arrivals. .: Example of market maker's placements and decisions she might make. In this example: the ask-side of the order book is described by a = (7, 5, 4); and the bid-side by b = (-6, -4, -5). The market maker positions are described by ra = (0, 1, -1) and rb = (0, 2, -1...); and the associated ranks vectors are na = (2, 1, -1 . . . ) and nb = (4, 2, -1). We have i0 = 0. After each order arrival, she can send new limit orders (see action on the top right), cancel some positions (see dashed cross on the bottom left), or just keep their orders unchanged.

Price Volume & Rank in the queues

Denote by (T n ) n∈N the sequence of jump times of the order book. We denote by A the set of the admissible strategies, defined as the predictable processes ra t , rb t t≤T such that the control is constant between two consecutive arrivals of orders from the participants, and such that the order of the market maker do not cross the spread. These conditions reads:

• for all n ∈ N, ra t , rb t ∈ {1, ..., K} M × {1, ..., K} M are constant on T n , T n+1 • ra * , rb * ≥ i0 where, for every vector a: a * = min By convention, we set in the sequel ra i (t) = -1 if the ith market maker's order is not placed in the order book.

2.2.2.

Controlled order book. The order book, controlled by the market maker, is fully described by the following state process Z:

Z t := X t , Y t , a t , b t , na t , nb t , pa t , pb t , ra t , rb t ,
where, at time t:

• X t is the cash held by the market maker on a zero interest account.

• Y t is the inventory of the market maker, i.e. it is the (signed) number of shares held by the market maker. • pa t is the ask price, i.e. the cheapest price a general market participant is willing to sell stock.

• pb t is the bid price, i.e. the highest price a general market participant is willing to buy stock.

• a t = (a 1 (t), . . . , a K (t)) (resp. b t = (b 1 (t), . . . , b K (t))) describes the ask (resp. bid) side: i ∈ {1, . . . , K}, a i (t) is the sum of all the general market participants' sell orders which are i ticks away from the bid (resp. ask) price. • ra t (resp. ra t ) describes the market maker's orders in the ask (resp. bid) side: for i ∈ {1, ..., M }, ra t (i) is the number of ticks between the i-th market maker's sell (resp. bid) order and the bid (resp. ask) price. By convention, we set ra t (i) = -1 (resp. rb t (i) = -1) if the ith sell (resp. buy) order of the market maker is not placed in the order book. As a result ra t (i), rb t (i) ∈ {1, . . . , K} ∪ {-1}. • na t (resp. nb t ) describes the ranks of the market maker's orders in the ask (resp. bid) side.

For i ∈ {1, ..., M }, na t (i) ∈ -1, ..., |a| + M (resp. nb t (i) ∈ -1, ..., |b| + M ) is the rank of the i-th sell (resp. buy) orders of the market maker in the queue. By convention, we assume that na t (i) = -1 (resp. nb t (i) = -1) if the i-th sell (resp. buy) order of the market maker is not placed in the order book.

3. Presentation of the market-making problem. Theoretical resolution.

Definition of the market-making problem and well-posedness of the value function

We denote by V the value function for the following market-making problem:

V (t, z) = sup α∈A E α t,z T t f α s , Z s ds + g Z T , (t, z) ∈ [0, T ] × E, (3.1)
where:

• A is the set of the admissible strategies, defined in Section 2.2.1.

• f and g are respectively the instantaneous and terminal reward functions.

• E α t,z stands for the expectation conditioned by Z t = z and when strategy α = (α s ) t≤s<T is followed on [t, T ].

Example 3.1 The terminal reward g can be defined as the sum of the market maker's terminal wealth function and an inventory penalization term, i.e. g : z → x + L(y) -ηy 2 where L is the amount earned from the immediate liquidation of the inventory 1 . We remind that z stands for a state of order book; ε is the tick size of the LOB; x is the value of the risk-free account of the market maker; η is the penalization parameter of the market maker; and where we remind that y stands for the (signed) market maker's inventory.

The running reward f can stand for a penalization of inventory term: f (z) := -γy 2 , with γ > 0.

We shall assume the following conditions on the rewards to insure the well-posedness of the market-making problem.

(Hrewards) The expectation of the integrated running reward is uniformly upper-bounded w.r.t. the strategies in A, i.e.

sup α∈A E α t,z T t f + (Z s , α s ) ds < +∞ holds;
where for all state z and action a, we denote f + (z, a) := max(f (z, a), 0). Moreover, the terminal reward g(Z T ) is a.s. at most linear with respect to the number of events up to time T , denoted by N T in the sequel, i.e. there exists a constant c 1 > 0 such as g(Z T ) ≤ c 1 N T , a.s.. Remark 3.1 Under Assumption (Hcontrol), Assumption (Hrewards) holds when g is defined as the wealth of the market maker plus an inventory penalization. In particular, we have g(Z T ) ≤ N T M , where M is the maximal number of orders that can be sent by the market maker, which holds a.s. since the best profit the market maker can make is when their buy (resp. sell) limit orders are all executed, and then the price keeps going to the right (resp. left) direction. Hence the second condition of (Hrewards) holds with c 1 = M .

The following Lemma 3.1 tackles the well-posedness of the control problem.

Lemma 3.1 Under (Hrewards) and (Hcontrol), the value function is well-defined, i.e.

sup α∈A E α t,z g(Z T ) + T t f α s , Z s ds < +∞,
where, as defined previously, E α t,z [.] stands for the expectation conditioned by the event {Z t = z}, assuming that strategy α ∈ A is followed in [t, T ].

1 L is defined as follows:

L(z) =    -1 k=1 a k (pa + k ) + (y -a 0 -... -a -1 )(pa +  ) if y < 0 --1 k=1 b k (pb -k ) + (y + b 0 + ... + b -1 )(pb - ) if y > 0 0 if y = 0,
for all z = x, y, a, b, na, nb, pa, pb, ra, rb , and where we define:

 := min j j i=1 a i > -y if y < 0 min j j i=1 |b i | > y if y > 0.
Proof. Denote by (N t ) t the sum of all the arrivals of orders up to time t. Under (Hrewards), we can bound E α t,z T t f α s , Z s ds + g(Z T ) , the reward functional at time t associated to a strategy α ∈ A, as follows:

E α t,z T t f (α s , Z s ) ds + g(Z T ) ≤ sup α∈A E α t,z [g(Z T )] + sup α∈A E α T t f + (Z s , α s ) ds ≤ c 1 sup α∈A E α t,0 [N T ] + sup α∈A E α t,z T t f + (Z s , α s ) ds , (3.2)
where once again, for all general process M and all m ∈ E, E α t,m [M T ] stands for the expectation of M T conditioned by M t = m and assuming that the market maker follows strategy α ∈ A in [t, T ]. Let us show that the first term in the r.h.s. of (3.2) is bounded. On one hand, we have:

E α t,0 [N T ] ≤ λ ∞ T 0 E (|a| t + |b| t ) dt, (3.3)
where λ ∞ := λ L + λ C + λ M is a bound on the intensity rate of N t . On the other hand, there exists a constant c 2 > 0 such that d(|a| + |b|) t ≤ c 2 dL t so that:

E α t,|a|0+|b|0 [|a| t + |b| t ] ≤ |a| 0 + |b| 0 + c 3 t 0 E [|a| s + |b| s ]
ds. Applying Gronwall's inequality, we then get:

E α t,|a|0+|b|0 [|a| t + |b| t ] ≤ (|a| 0 + |b| 0 ) e c3t . (3.4)
Plugging (3.4) into (3.3) finally leads to:

E α t,0 [N T ] ≤ c 4 e c3T
wit c 3 and c 4 > 0 that do not depends on α, which proves that the first term in the r.h.s. of (3.2) is bounded. Also, its second term in the r.h.s. of (3.2) is bounded under (Hrewards). Hence, the reward functional is bounded uniformly in α, which proves that the value function of the considered market-making problem is well-defined.

Markov Decision Process formulation of the market-making problem

In this section, we first reformulate the market-making problem as a Markov Decision Process (MDP), and then characterize the value function as solution to a Bellman equation.

Let us denote (T n ) is the increasing sequence of the arrivals of market/limit/cancel order to the market; and let Z n := φ a(ZT n ) (Z Tn ), where φ a (z) ∈ E is the state of the order book at time t such that T n < t < T n+1 , given that Z Tn = z and given that the strategy a has been chosen by the market maker at time T n .

Let us consider the Markov Decision Process (T n , Z n ) n∈N , which is characterized by the following information

[0, T ] × E state space , A z market maker control , λ intensity of the jump , Q transitions kernel
, r reward where: y, a, b, na, nb, ra, rb, pa, pb where: x is the cash held by the market maker, y their inventory; a and b, introduced in Section 2.2.2. represent the orders in the ask and bid sides of the order book of all the participants except the market maker's; na (resp. nb) is the Mdimensional vector of the ranks of the market maker's sell (resp. buy) orders in the queues ; ra (resp. rb) is the M -dimensional vector of the number of ticks the M market maker's sell (resp. buy) orders are from the bid (resp. ask) price; pa (resp. pb) is the ask-price (resp. bid-price). • A z , for every state z ∈ E, is the set of the admissible actions (i.e.the actions the market maker can take) when the order book is at state z:

• [0, T ] × E is the state space of the time-continuous controlled process (T n , Z n ) n∈N ; and E := R ×N × N K × N K × N M × N M × N M × N M × R × R is the state space of (Z t ). For z ∈ E, z = x,
A z = ra, rb ∈ {1, ..., K} M × {1, ..., K} M rb * , ra * ≥ i0 ,
where we define c * = min 1≤i≤K {c i |c i = -1} and c0 = argmin

1≤i≤K {c i > 0} for c ∈ N M .
We recall that this condition means that the market maker is not allowed to cross the spread. • λ is the intensity of the controlled process (Z t ), and reads:

λ(z) := λ M + (z) + λ M -(z) + 1≤j≤K λ L + j (z) + 1≤j≤K λ L - j (z) + 1≤j≤K λ C + j (z) + 1≤j≤K λ C - j (z).
Observe that λ does not depend on the strategy α chosen by the market maker since we assumed that the general participants does not "see" the market maker's orders in the order book. Although we wrote z as argument for he intensity of the order book process, it cannot depend on any controlled component variable of the latter. To simplify, the reader can assume that the intensities only depend on the vectors a and b. • Q is the transition kernel of the MDP, which is defined as follows:

Q B ×C|t, z, α := λ(z) T -t 0 e -λ(z)s 1 B (t+s)Q C|φ α (z), α ds+e -λ(z)(T -t) 1 T ∈B,z∈C , (3.5)
for all Borelian sets B ⊂ R + and C ⊂ E, for all (t, z) ∈ [0, T ] × E, for all α ∈ A, and where Q is the transition kernel of (Z t ) defined for all state z as:

Q z |z, u =          λ M + (z) λ(z) if z = e M + (φ u (z)) . . . λ C + (z) λ(z) if z = e C + K (φ u (z)),
where φ u (z) is the new state of the controlled order book when decision u as been taken and when the order book was at state z before the decision; e M + (z) is the new state of the order book right after it received a buy market order, given that it was at state z before the jump; and e C ± i (z) is the new state of the order book right after it received a cancel order from a general market participant on its i th ask/bid limit, given that it was at state z.

• r : [0, T ] × E C → R is the running reward associated to the MDP with infinite horizon defined as follows:

r(t, z, a) := -c z, a e -λ(z)(T -t) (T -t)1 t>T + c z, a 1 λ(z) - e -λ(z)(T -t) λ(z) + e -λ(z)(T -t) g(z)1 t≤T , (3.6)
and its definition is motivated by Proposition 3.1 below.

The cumulated reward functional associated to the MDP (T n , Z n ) n∈N for an admissible policy (f n ) ∞ n=0 is defined as:

V ∞,(fn) (t, z) = E (fn) t,z ∞ n=0 r T n , Z n , f n (T n , Z n ) ,
and the associated value function is the supremum of the cumulated reward functional over all the admissible controls in A, i.e.

V ∞ (t, z) = sup (fn) ∞ n=0 ∈A V ∞,α (t, z), (t, z) ∈ [0, T ] × E, (3.7)
Notice that we used the same notation for admissible controls of the MDP and those of the continuous-time control problem.

Remark 3.2 Q is defined as in (3.5) because

P T n+1 -T n ≤ t, Z n+1 ∈ B|T 0 , Z 0 , . . . , T n , Z n = λ(Z n ) t 0 e -λ(Zn)s Q B|Z Tn , α Tn ds = λ(Z n ) t 0 e -λ(Zn)s Q B|Z Tn , f n (Z n ) ds,
holds for any admissible policy α = (f n ) ∞ n=0 ∈ A, for all Borelian B ⊂ E, and for all t ∈ [0, T ].

In the sequel, we denote [0, T ]×E

C := t, z, a ∈ E ×{1, . . . , K} 2 M t ∈ [0, T ], z ∈ E, a ∈ A z ,
and

E C := z, a ∈ E × {1, . . . , K} 2 M z ∈ E, a ∈ A z . Q is the stochastic kernel from E C to E
that describes the distribution of the jump goals, i.e., Q B|z, u is the probability that the order book jumps in the set B given that it was at state z ∈ E right before the jump, and the control action u ∈ A z has been chosen right after the jump time.

Remark 3.3 The MDP is defined in such a way that the control is feedback and constant between two consecutive arrivals of market/limit/cancel orders in the market, i.e. in the time-continuous setting: we restrict ourselves to the control α = (α t ) which are entirely characterized by the decision functions f n : [0, T ] × E → A, and such that

α t = f n (T n , Z n ) for t ∈ T n , T n+1
By abuse of notation, we denote in the sequel by α the sequence of controls

(f n ) ∞ n=0 .
The following Proposition 3.1 motivates the special choice of the running reward r as defined in (3.6):

Proposition 3.1 The value function of the MDP defined by (3.7) coincides with (3.1), i.e. we have for all (t, z) ∈ E C :

V ∞ (t, z) = V (t, z).
(3.8)

Proof. Let us show that for all α = (f n ) ∈ A and all (t, z) ∈ E C V α (t, z) = V (fn) ∞ (t, z). (3.9)
Let us first denote by H n := (T 0 , Z 0 , ..., T n , Z n ). Notice then that for all admissible strategy α:

V α (t, z) = E α t,z ∞ n=0 1 T >Tn+1 T n+1 -T n c Z n , α n + 1 [Tn≤T <Tn+1) g(Z T ) -ηY T 2 + (T -T n )c Z n , α n = ∞ n=0 E (fn) t,z r T n , Z n , f n (T n , Z n ) , (3.10)
where we conditioned by H n between the first and the second line. We recognize V (fn) ∞ in the r.h.s. of (3.10), so that the proof of (3.9) is completed.

It remains to take the supremum over all the admissible strategies A in (3.9) to get (3.8).

From Proposition 3.1, we deduce that the value function of the market-making problem is the same as the value function V ∞ of the discrete-time MDP with infinite horizon. We now aim at solving the MDP control problem. To proceed, we first define the maximal reward mapping for the infinite horizon MDP:

(T v)(t, z) := sup a∈Az r(t, z, a) + v(t , z )Q(t , z |t, φ a (z), a) = sup a∈Az r(t, z, a) + λ(z) T -t 0 e -λ(z)s v(t + s, z )Q dz |φ a (z), a ds , (3.11)
where we recall that:

• φ α (z) is the new state of the order book when the market maker follows the strategy α and the order book is at state z before the decision is taken. • λ(z) is the intensity of the order book process given that the order book is at state z.

We shall tighten assumption (Hrewards) in order to guarantee existence and uniqueness of a solution to (3.1), as well as characterizing the latter.

(HrewardsBis): The running and terminal rewards are at most quadratic w.r.t. the state variable, uniformly w.r.t. the control variable, i.e.

(i) The running reward f is such that |c| is uniformly bounded by a quadratic in z function, i.e. there exists c 5 > 0 such that:

∀(z, a) ∈ E × A, |f (z, a)| ≤ c 5 (1 + |z| 2 ).
(ii) The terminal reward g has no more than a quadratic growth, i.e. there exists c 6 > 0 such that:

∀z ∈ E, |g(z)| ≤ c 6 (1 + |z| 2 ).
Remark 3.4 Assumption (HrewardsBis) holds in the case where g is the terminal wealth of the market maker plus a penalization of their inventory, and where with no running reward, i.e. f = 0.

The main result of this section is the following theorem that gives existence and uniqueness of a solution to (3.1), and moreover characterizes the latter as fixed point of the maximal reward operator defined in (3.11).

Theorem 3.1 T admits a unique fixed point v which coincides with the value function of the MDP. Moreover we have:

v = V ∞ = V.
Denote by f * the maximizer of the operator T . Then f * , f * , ... is an optimal stationary (in the MDP sense) policy.

Remark 3.5 Theorem 3.1 states that the optimal strategy is stationary in the MDP formulation of the problem, but of course, it is not stationary for the original time-continuous trading problem with finite horizon (3.1), since the time component is not a state variable anymore in the original formulation. Actually, given n ∈ N and the state of order book z at that time, the optimal decision to take at time T n is given by f * T n , z .

We devote the next section to the proof of Theorem 3.1.

Proof of Theorem 3.1

Remind first that we defined in the previous section

E C := z, a ∈ E × {1, . . . , K} 2 M z ∈ E, a ∈ A z and [0, T ] × E C := t, z, a ∈ [0, T ] × E × {1, . . . , K} 2 M t ∈ [0, T ], z ∈ E, a ∈ A z .
Definition 3.1 A measurable function b : E → R + is called a bounding function for the controlled process (Z t ) if there exists positive constants c c , c g , c Q , c φ such that:

(i) |f (z, a)| ≤ c c b(z) for all (z, a) ∈ E C . (ii) |g(z)| ≤ c g b(z) for all z in E. (iii) b(z )Q (dz |z, a) ≤ c Q b(z) for all (z, a) ∈ E C . (iv) b(φ α t (z)) ≤ c φ b(z) for all (t, z, α) ∈ [0, T ] × E C .
Proposition 3.2 Let b be such that :

∀z ∈ E, b(z) := 1 + |z| 2 .
Then, b is a bounding function for the controlled process (Z t ), under Assumption (HrewardsBis).

Proof. Let us check that b defined in Proposition 3.2 satisfies the four assertions in Definition 3.1.

• Assertion 1 and 2 of Definition 3.1 holds under (HrewardsBis).

• First notice that ra, rb are bounded by √ M K (where we recall that K is the number of limits in each side of the order book, and M is the biggest number of limit orders that the market maker is allowed to send in the market). Secondly, pa ∈ B(pa, K), pb ∈ B(pb, K), where B(x, r) is the ball centered in x with radius r > 0, because of the limit conditions that we imposed in our LOB model. And last, we can see that |a | ≤ |a| + a ∞ K. These three bounds are linear w.r.t. z so that assertion 3 holds. Let us define

• φ α (z) = z α
Λ := (4K + 2)sup λ M ± |a| + |b| , λ L ± |a| + |b| , λ C ± |a(z)| + |b(z)| ,
which is well-defined under (Harrivals). t) , with γ(z) = γ 0 (4K + 2)Λ 1 + |a| + |b| and γ 0 > 0 is a bounding function for the MDP, i.e. for all t ∈ [0, T ], z ∈ E, a ∈ A z , we have:

Proposition 3.3 If b is a bounding function for (Z t ), then b(t, z) := b(z)e γ(z)(T -
|r(t, z, a)| ≤ c g b(t, z), b(s, z )Q(ds, dz |t, z, a) ≤ c φ c Q e C(T -t) 1 1 + γ 0 b(t, z), with C = γ 0 ΛK(4K + 2) |a| ∞ + |b| ∞ .
Proof. Let z = x , y , a , b , na , nb , ra , rb be the state of the order book after an exogenous jump occurs given that it was in state z before the jump. Since

|a | ≤ |a| + a ∞ K and |b | ≤ |b| + b ∞ ,
where a ∞ and b ∞ are defined as the border conditions of the order book, we have:

γ(z ) ≤ γ(z) + C, (3.12) with C = γ 0 ΛK(4K + 2)(a ∞ + b ∞ ). Then, we get: b(s, z )Q(ds, dz |t, φ α (z), α) = λ(z) T -t 0 e -λ(z)s b(t + s, z )Q dz |φ α s (z), α ds = λ(z) T -t 0 e -λ(z)s b(z )e γ(z )(T -(t+s)) Q dz |φ α s (z), α ds ≤ λ(z) T -t 0 e -λ(z)s b(z )e (γ(z)+C)(T -(t+s)) Q dz |φ α s (z), α ds ≤ λ(z) T -t 0 e -λ(z)s e (γ(z)+C)(T -(t+s)) b(z )Q dz |φ α s (z), α ds ≤ λ(z) T -t 0 e -λ(z)s e (γ(z)+C)(T -(t+s)) c Q c φ b(z)ds ≤ λ(z)c Q c φ λ(z) + γ(z) + C e (γ(z)+C)(T -t) 1 -e -(T -t)(λ(z)+γ(z)+C) b(z) ≤ c Q c φ λ(z) λ(z) + γ(z) + C e C(T -t) 1 -e -(T -t)(λ(z)+γ(z)+C) b(t, z),
where we applied (3.12) at the third line. It remains to notice that

λ(z) λ(z) + γ(z) + C = λ(z) λ(z) 1 + γ 0 + γ 0 Λ(|a| + |b|) -λ(z) ≥0 ≤ 1 1 + γ 0 ,
to complete the proof of the proposition.

Let us denote by . b the weighted supremum norm such that for all measurable function v :

E → R, v b := sup (t,z)∈E |v(t, z)| b(t, z) ,
and define the set:

B b := v : E → R|v is measurable and v b < ∞ .
Moreover let us define

α b := sup (t,z,α)∈E ×R b(s, z )Q(ds, dz |t, φ α (z), α) b(t, z) .
From the preceding estimations we can bound α b as follows:

α b ≤ c Q c φ 1 1 + γ 0 e CT ,
so that, by taking: γ 0 = c Q c φ e CT , we get: α b < 1. In the sequel, we then assume w.l.o.g. that α b < 1. Recall that the maximal reward mapping for the MDP has been defined as:

T v : (t, z) → sup a∈Az r(t, z, a) + λ(z) T -t 0 e -λ(z)s v(t + s, z )Q dz |φ a (z), a ds
It is straightforward to see that:

T v -T w b ≤ α b v -w b , (3.13)
which implies that T is contracting, since α b < 1.

Let M be the set of all the continuous function in B b . Since b is continuous, (M, . b ) is a Banach space.

T sends M to M. Indeed, for all continuous function v in B b , (t, z, a) → r(t, z, a) + λ(z)

T -t 0 e -λ(z)s v(t + s, z )Q dz |φ a (z), a ds is continuous on [0, T ] × E C . A z
is finite, so we get the continuity of the application:

T v : (t, z) → sup a∈Az r(t, z, a) + λ(z) T -t 0 e -λ(z)s v(t + s, z )Q dz |φ a (z), a ds .
Proposition 3.4 There exists a maximizer for T , i.e. let v ∈ M, then there exists a Borelian function f : [0, T ] × E → A such that for all (t, z) ∈ E :

T v t, z, f t, z = sup a∈A r(t, z, a) + λ(z) T -t 0 e -λ(z)s v(t + s, z )Q dz |φ a (z), a ds Proof. D * (t, z) = a ∈ A T a v(t, z) = T v(t, z) is finite, so it is compact. So (t, z) → D * (t, z) is a compact-valued mapping. Since the application (t, z, a) → T a (t, z) -T (t, z) is continuous, we get that D * = (t, z, a) ∈ E C T a v(t, z) = T v(t,
sup α∈A E α t,z ∞ k=n |r(T k , Z k )| ≤ α n b 1 -α b b(t, z),
and in particular, we have:

lim n→∞ sup α∈A E α t,z ∞ k=n r t k , Z k = 0.
Proof. By conditioning we get

E α t,z r(T k , Z k ) ≤ c g α b k b(t, z) for k ∈ N,
and for all α ∈ A. It remains to sum this inequality to complete the proof of Lemma 3.2.

We can now prove Theorem 3.1.

Proof. We divided the proof of Theorem 3.1 into four steps.

Step 1: Inequality (3.13) and Proposition 3.3 imply that T is a stable and contracting operator defined on the Banach space M. Banach's fixed point theorem states that T admits a fixed point, i.e. there exists a function v ∈ M such that v = T v, and moreover we have v = lim n→∞ T n 0. Notice that T N 0 coincides with v 0 defined recursively by the following Bellman equation:

v N = 0 v n = T v n+1 for n = N -1, ..., 0. (3.14)
The solution of the Bellman equation is always larger than the value function of the MDP associated (see e.g. Theorem 2.3.7 p.22 in [START_REF] Bäuerle | Markov Decision Processes with Applications to Finance[END_REF]). Then we have:

T n 0 ≥ sup (fk) E (fk) n n-1 k=0 r(t k , X k ) =: J n
, where J n is the value function of the MDP with finite horizon n and terminal reward 0, associated to (3.14). Moreover, by Lemma 7.1.4 p.197 in [START_REF] Bäuerle | Markov Decision Processes with Applications to Finance[END_REF], we know that J n n converges as n → ∞ to a limit that we denote by J. Passing at the limit in the previous inequality we get: lim n→∞ T n 0 ≥ J, i.e.

v ≥ J.

(3.15)

Step 2: Let us fix a strategy α ∈ A, and take n ∈ N. We denote J n (α) := E (αk) 0

n-1 k=0 r(t k , X k ) , the reward functional associated to the control α on the discrete finite time horizon {0, . . . , n}. By definition, we have J n (α) ≤ J n . We get by letting n → ∞: lim n→+∞ J n (α) =: J ∞ (α) ≤ J. Taking the supremum over all the admissible strategies α finally leads to:

V ∞ ≤ J.
(3.16)

Step 3: Let us denote by f a maximizer of T associated to v, which exists, as stated in Proposition 3.4. v is the fixed point of

T so that v = T n f (v), for n ∈ N. Moreover v ≤ δ where δ := sup α∈A E ∞ k=0 r + (Z k , α k ) , so that T n f (v) ≤ T n f 0 + T n o δ, where T n o δ = sup α E α n ∞ k=n r + (t k , Z k ) . Lemma 3.2 implies that T n o δ → 0 as n → ∞. Hence, we get: v ≤ J f .
(3.17)

Step 4: Conclusion. Since

J f ≤ V ∞ (3.18)
holds, we get by combining (3.15), (3.16), (3.17) and (3.18):

V ∞ ≤ J ≤ v ≤ J f ≤ V ∞ . (3.19)
All the inequalities in (3.19) are then equalities, which completes the proof of Theorem 3.1.

Numerical resolution of the market-making control problem

In this section, we first introduce an algorithm to numerically solve a general class of discrete-time control problem with finite horizon, and then apply it on the trading problem (3.1).

Framework

Let us consider a general discrete-time stochastic control problem over a finite horizon N ∈ N \ {0}.

The dynamics of the controlled state process

Z α = (Z α n ) n valued in R d is given by Z α n+1 = F (Z α n , α n , ε n+1 ), n = 0, . . . , N -1, Z α 0 = z ∈ R d , (4.1)
with (ε n ) n is a sequence of i.i.d. random variables valued in some Borel space (E, B(E)), and defined on some probability space (Ω, F, P) equipped with the filtration

F = (F n ) n generated by the noise (ε n ) n (F 0 is the trivial σ-algebra), the control α = (α n ) n is an F-adapted process valued in A ⊂ R q , and F is a measurable function from R d × R q ×E into R d .
Given a running cost function f defined on R d × R q , a terminal cost function g defined on R d , the cost functional associated to a control process α is

J(α) = E N -1 n=0 f (Z α n , α n ) + g(Z α N ) .
The set A of admissible controls is the set of control processes α satisfying some integrability conditions ensuring that the cost functional J(α) is well-defined and finite. The control problem, also called Markov decision process (MDP), is formulated as

V 0 (x 0 ) := sup α∈A J(α), (4.2)
and the goal is to find an optimal control α * ∈ A, i.e., attaining the optimal value: V 0 (z) = J(α * ).

Notice that problem (4.1)-(4.2) may also be viewed as the time discretization of a continuous time stochastic control problem, in which case, F is typically the Euler scheme for a controlled diffusion process.

Problem (4.2) is tackled by the dynamic programming approach. For n = N, . . . , 0, the value function V n at time n is characterized as solution of the following backward (Bellman) equation:

V N (z) = g(z) V n (z) = sup a∈A f (z, a) + E a n,z [V n+1 (Z n+1 )] , z ∈ R d , (4.3)
Moreover, when the supremum is attained in the DP formula at any time n by a * n (z), we get an optimal control in feedback form given by: α * = (a * n (Z * n )) n where Z * = Z α * is the Markov process defined by

Z * n+1 = F (Z * n , a * n (Z * n ), ε n+1 ), n = 0, . . . , N -1, Z * 0 = z.
There are two usual ways that have been studied in the literature, to solve numerically (4.3): one way is to use local regression methods, relying e.g. on quantization, k-nearest neighbors or kernel ideas to approximate the conditional expectations by cubature methods; another way is to rely on MC regress-now or later methods to regress the value functions V n+1 at time n for n = 0, . . . , N -1 on basis functions or neural networks. See e.g. [START_REF] Kharroubi | A numerical algorithm for fully nonlinear HJB equations: an approach by control randomization[END_REF] for the regress-now and [START_REF] Balata | Regress-Later Monte Carlo for optimal control of Markov processes[END_REF] for the regress-later methods for algorithms using basis functions, and e.g. [START_REF] Huré | Deep neural networks algorithms for stochastic control problems on finite horizon: convergence analysis[END_REF] for regression on neural networks based on regress-now or regress-later techniques.

Presentation and rate of convergence of the Qknn algorithm

In this section, we present an algorithm based on k-nn estimates for local non-parametric regression of the value function, and optimal quantization to quantize the exogenous noise, in order to numerically solve (4.3).

Let us first introduce some ingredients of the quantization approximation:

• We denote by ε a K-quantizer of the E-valued random variable ε n+1 ∼ ε 1 , that is a discrete random variable on a grid Γ = {e 1 , . . . , e K } ⊂ E K defined by ε = Proj Γ (ε 1 ) := K =1 e l 1 ε1∈Ci(Γ) ,
where C 1 (Γ), . . The grid points (e ) which minimize the L 2 -quantization error ε 1 -ε 2 lead to the so-called optimal L-quantizer, and can be obtained by a stochastic gradient descent method, known as Kohonen algorithm or competitive learning vector quantization (CLVQ) algorithm, which also provides as a byproduct an estimation of the associated weights (p ). We refer to [START_REF] Pagès | Optimal quantization methods and applications to numerical problems in finance[END_REF] for a description of the algorithm, and mention that for the normal distribution, the optimal grids and the weights of the Voronoi tesselations are precomputed and available on the website http://www.quantize.maths-fi.com • Recalling the dynamics (4.1), the conditional expectation operator is equal to

P âM n (z) W (x) = E W (Z âM n n+1 )|Z n = x = E W (F (z, âM n (z), ε 1 )) , z ∈,
that we shall approximate analytically by quantization via:

P âM n (z) W (z) := E W (F (z, âM n (z), ε)) = K =1 p W (F (z, âM n (z), e )).
Let us secondly introduce the notion of training distribution that will be used to build the estimators of value functions at time n, for n = 0, . . . , N -1. Let us consider a measure µ on the state space E. We refer to it in the sequel as the training measure. Let us take a large integer M , and for n = 0, . . . , N , introduce Γ n = Z The proposed algorithm reads as:

       V Q N (z) = g(z), for z ∈ Γ N , Qn (z, a) = K =1 p f (z, a) + V Q n+1 Proj n+1 F z, e , a , V Q n (z) = sup a∈A Qn (z, a), for z ∈ Γ n , n = 0, . . . , N -1. (4.4)
where, for n = 0, . . . , N , Proj n (z) stands for the closest neighbor of z ∈ E in the grid Γ n , i.e. the operator z → Proj n (z) is actually the euclidean projection on the grid Γ n .

In the sequel, we refer to (4.4) as the Qknn algorithm.

We shall make the following assumption on the transition probability of (Z n ) 0≤n≤N , to guarantee the convergence of the Qknn algorithm.

(Htrans) Assume that the transition probability P(Z n+1 ∈ A Z n = z, a) conditioned by Z n = z when control a is followed at time n admits a density r w.r.t. the training measure µ, which is uniformly bounded and lipschitz w.r.t. the state variable z, i.e. there exists r ∞ > 0 such that for all z ∈ E and control u taken at time n: Denote by Supp(µ) the support of µ. We shall assume smoothness conditions on µ and F to provide a bound on the projection error.

(Hµ) We assume Supp(µ) to be bounded, and denote by µ ∞ the smallest real such that Supp(µ) ⊂ B (0, µ ∞ ). Moreover, we assume x ∈ E → µ B(x, η) to be Lipschitz, uniformly w.r.t. η, and we denote by [µ] L its Lipschitz constant.

(HF) For x ∈ E and a ∈ A, assume F to be L 1 -Lipschitz w.r.t. the noise component ε, i.e., there exists [F ] L > 0 such that for all x ∈ E and a ∈ A, for all r.v. ε and ε , we have:

E F (x, a, ε) -F (x, a, ε ) ≤ [F ] L E ε -ε
We now state the main result of this section whose proof is postponed in Appendix A.

Theorem 4.1 Take K = M 2+d points for the optimal quantization of the exogenous noise ε n , n = 1, . . . , N . There exist constants [ V Q n ] L > 0, that only depends on the Lipschitz coefficients of f , g and F , such that, under (Htrans), it holds for n = 0, ..., N -1, as M → +∞:

V Q n (X n ) -V n (X n ) 2 ≤ N k=n+1 r N -k ∞ V Q k L ε proj k + [F ] L ε Q k + O 1 M 1/d , (4.5)
where ε Q k := εk -ε k 2 stands for the quantization error, and

ε proj n := sup a∈A Proj n+1 (F (X n , a, εn )) -F (X n , a, εn ) 2
stands for the projection error, when decision a is taken at time n.

Remark 4.1 The constants [ V Q n ] L > 0 are defined in (A8).
From Theorem 4.1, we can deduce consistency and provide a rate of convergence for the estimator V Q n , n = 0, . . . , N -1, under some rather tough yet usual compactness conditions on the state space.

Corollary 4.1 Under (Hµ) and (HF), the Qknn-estimator V Q n is consistent for n = 0, . . . , N -1, when taking M d+1 points for the quantization; and moreover, we have for n = 0, ..., N -1, as M → +∞:

V Q n (X n ) -V n (X n ) 2 = O 1 M 1/d .
Proof. We postpone the proof of Theorem 4.1 to Appendix A.

Qknn agorithm applied to the order book control problem (3.1)

We recall the expression of the time-continuous controlled order book process Z t = X t , Y t , a t , b t , na t , nb t , pa t , pb t , ra t , rb t admits a representation as a MDP as shown in Section 3.2. In Section 3.3, we proved that the value function associated to the MDP is characterized as limit as N → ∞ of the the Bellman equation (3.14). In this section, some implementation details on the Qknn algorithm are presented in order to numerically solve (3.14).

Training set design

Inspired by [START_REF] Fiorin | Product Markovian quantization of a diffusion process with applications to finance[END_REF], we use product-quantization method and randomization techniques to build the training set Γ n on which we project (T n , Z n ) that lies on [0, T ] × E, where T n and Z n stands for the n th jump of Z and the state of Z at time t n , i.e. Z n = Z Tn , for n ≥ 0. This basic idea of Control Randomization consists in replacing in the dynamics of Z the endogenous control by an exogenous control (I Tn ) n≥0 , as introduced in [START_REF] Kharroubi | A numerical algorithm for fully nonlinear HJB equations: an approach by control randomization[END_REF]. In order to alleviate the notations, we denote by I n the control taken at time T n , for n ≥ 0.

Initialization. Set: Γ E 0 = {z} and Γ T 0 = {0}. Randomize the control, using e.g. uniform distribution on A at each time step, and then simulate

D randomized processes to generate (T k n , Z k n ) N,D n=0,k=1 . For all n = 1, . . . , N , set Γ T n = {T k n , 1 ≤ k ≤ D},
which stands for the grid associated to the quantization of the n th jump time T n , and set Γ E n = {Z k n , 1 ≤ k ≤ D} which stands for the grid associated to the quantization of the state Z n of Z at time T n .

Remark 4.2 The way we chose our training sets is often referred to as an exploration strategy in the reinforcement learning literature. Of course, if one has ideas or good guess of where to optimally drive the controlled process, she should not follow an exploration-type strategy to build the training set, but should rather use the guess to build it, which is referred to as the exploitation strategy in the reinforcement learning and the stochastic bandits literature. We refer to [START_REF] Balata | A class of finite-dimensional numerically solvable McKean-Vlasov control problems[END_REF] for several other applications of the exploration strategy to build training sets. Note that this idea is the root of all the Q-learning based algorithms. See [START_REF] Sutton | Reinforcement Learning: An Introduction, Second[END_REF] for more details on Q-learning.

Let F and G be the Borelian functions such that Z n = F Z n-1 , d n , I n and T n = G T n-1 , n , I n , where n ∼ E(1) stands for the temporal noise, and d n is the state noise, for n ≥ 0.

Let us fix N ≥ 1 and consider T n , Z n N n=0 , the dimension-wise projection of T n , Z n N n=0 on the grids Γ T n × Γ E n , n = 0, . . . , N , i.e. T 0 = 0, Z 0 = z, and

   T n = Proj G T n-1 , n , I n , Γ T n , Z n = Proj F Z n-1 , d n , I n , Γ E n , for n = 1, . . . , N. T n , Z n , I n n∈{0,N } is a Markov chain. Define then T Q n , Z Q n N n=0 as temporal noise-quantized version of T n , Z n , I n N n=0
. Note that we do not need to quantize the spacial noise since this noise already takes a finite number of states. Let εn be the quantized process associated to n . The process

T Q n , Z Q n N n=0
is then defined as follows:

Z Q 0 = z, T Q 0 = 0 and ∀1 ≤ n ≤ N :    T Q n = Proj G tn-1 , εn , I n , Γ T n , Z Q n = Proj F Z n-1 , d n , I n , Γ E n . Denote by V Q,(N,D) n N n=0
the solution of the Bellman equation associated to

T Q n , Z Q n N n=0 : (B Q N,D ) :    V Q,(N,D) N = 0 V Q,(N,D) n (t, z) = r(t, z, a) + sup a∈A E a t,z V Q,(N,D) n+1 T Q n+1 , Z Q n+1
, for n = 0, . . . , N,

where E a t,z [.] stands for the expectation conditioned by the events

T Q n = t,Z Q n = z and when decision I n = a is taken at time t.
We wrote the pseudo-code of the Qknn algorithm to compute (B

Q N,D ) in Algorithm 1.
Algorithm 1 Generic Qknn Algorithm Inputs:

-N : number of time steps -z: state in E at time T 0 = 0 -Γ ε = {e 1 , . . . , e L } and (p ) L =1 : the grid and the weights for the optimal quantization of (ε n ) N n=1 . -Γ n and Γ E n the grids for the projection of respectively the time and the state components at time n, for n = 0, . . . , N .

1: for n = N -1, . . . , 0 do 2:

Compute the approximated Qknn-value at time n:

Qn (z, a) = r (T n , z, a) + L =1 p V Q n+1 Proj G(z, e , a), Γ T n+1 , Proj F (z, e , a), Γ E n+1 , for (z, a) ∈ Γ n × A z ; 3:
Compute the optimal control at time n

Ân (z) ∈ argmin a∈Az Qn (z, a), for z ∈ Γ n ,
where the argmin is easy to compute since A z is finite for all z ∈ E;

4:

Estimate analytically by quantization the value function:

V Q n (z) = Qn z, Ân (z) , ∀z ∈ Γ n ; 5: end for Output: -( V Q 0 ): Estimate of V (0, z);
We discuss in Remark 4.3 the reasons why we can apply Theorem 4.1.

Remark 4.3 When the number of jumps of the LOB N ≥ 1 is fixed, the set of all the states that can take the controlled order book by jumping less than N times, denoted by K in the sequel, is finite. Hence, the reward function r, defined in (3.6), is bounded and Lipschitz on K.

The following proposition states that V Q,(N,D) n

, built from the combination of time-discretization, k-nearest neighbors and optimal quantization methods, is a consistent estimator of the value function at time T n , for n = 0, . . . , N -1. It provides a rate of convergence for the Qknn-estimations of the value functions.

Proposition 4.1 The estimators of the value functions provided by Qknn algorithm are consistent. Moreover, it holds as M → +∞:

V Q,(N,M ) n T n , Z n -V n T n , Z n M,2 = O α N + 1 M 2/d , for n = 0, . . . , N -1,
where we denote by . M,2 the L 2 (µ) norm conditioned by the training sets that have been used to build the estimator V Q,(N,M ) n+1

.

Proof. Splitting the error of time cutting and quantization, we get:

V n T n , Z n -V (N,M ) n T n , Z n M,2 ≤ V n T n , Z n -V (N ) n T n , Z n M,2 + V (N ) n T n , Z n -V (N,M ) n T n , Z n M,2 . (4.8)
Step 1: Applying Lemma 3.2, we get the following bound on the first term in the r.h.s. of (4.8):

V n T n , Z n -V (N ) n T n , Z n M,2 ≤ α N 1 -α b ∞ , (4.9)
where b ∞ stands for the supremum of b over [0, T ] × E.

Step 2: Note that the assumptions of Theorem 4.1 are met as noticed in Remark 4.3, so that the latter provides the following bound for the second term in the r.h.s. of (4.8):

V (N ) n T n , Z n -V Q,(N,M ) n T n , Z n M,2 = M →∞ O 1 M 2/d . (4.10)
It remains to plug (4.9) and (4.10) into (4.8) to complete the proof of Proposition 4.1.

Numerical results

In this section, we propose several settings to test the efficiency of Qknn on simulated order books. We take no running reward, i.e. f = 0, and take the wealth of the market maker after liquidating their inventory as terminal reward, i.e. g(z) = x + L(y). The intensities are taken constant in the first tests, and state dependent in the second ones. The values of the state dependent intensities are similar to the ones in [START_REF] Huang | Simulating and analyzing order book data: The queue-reactive model[END_REF] (see ). Although the intensities are assumed to be uncontrolled in section 3 for predictability reasons, the latter are controlled processes in this section, i.e. the intensities of the order arrivals depends on the orders in the order book from all the participant plus the ones of the market maker. The optimal trading strategies have been computed 13). The intensities are estimated by the authors using their data from Alcatel Lucent.

among two different classes of strategies: in Section 4.4.1, we tested the algorithm to approximate the optimal strategy among those where the market maker is only allowed to place orders only at the best bid and the best ask; in Section 4.4.2, we computed the optimal trading strategy among the class of the strategies where the market maker allows herself to place orders on the two best limits on each side of the order book. Note that the second class of controls is more general than the first one. The code is available on https://github.com/comeh.

The search of the k nearest neighbors, that arises when estimating the conditional expectations using the Qknn algorithm, is very time-consuming; especially in the considered market-making problem which is of dimension more than 10. The efficiency of Qknn then highly depends on the algorithm used to find the k nearest neighbors in high-dimension. We implemented the method using the Fast Library for Approximate Nearest Neighbors algorithm (FLANN), introduced in [START_REF] Muja | Fast Approximate Nearest Neighbors with Automatic Algorithm Configuration[END_REF] and already available as a library of C++, Python, Julia and many other languages. This algorithm is based on tree methods. Note that recent algorithms based on graph also proved to perform well and can also be used.

Case 1:

The market maker only place orders at the best ask and best bid. Denote by A1lim the class of controls where the placements of orders are allowed on the best ask and best bid exclusively. We implement the Qknn algorithm to compute the optimal strategy among those in A1lim. We then compared the optimal strategy with a naive strategy which consists in always placing one order at the best bid and one order at the best ask. The naive strategy is called 11 in the plots, and can be seen as a benchmark. The naive strategy is a good benchmark when the model for the intensities of order arrivals is symmetrical, i.e. the intensities for the bid and the ask sides are the same. Indeed, in this case, the market maker can expect to earn the spread in average.

In Figure 3a, we take constant intensities to model the limit and market orders arrivals, and linear intensity to model the cancel orders. In this setting, as we can see in the figure, the strategy computed using Qknn algorithm performs as well as the naive strategy. Note that, obviously, the market maker has to take enough points for the state quantization in order for Qknn algorithm to perform well. In Figure 3b, we plotted the P&L of the market maker when the latter compute the optimal strategy using only 6000 points for the state space discretization, and for such a low number of points for the grid, Qknn algorithm performs poorly. In this setting, notice that the 11). We took symmetrical and constant intensities, and a short terminal time T = 1. Notice that the Qknn strategy looks to improve the P&L by reducing the losses when enough points are taken to build the grids (see Figure 3a), and that its performance is worse if less points are taken to build the grids (see Figure 3b). naive strategy seems to perform well.

In Figure 4, we plotted the empirical histogram of the P&L of the market maker using the Qknn-estimated optimal strategy, computed with grids of size N = 10 3 , 10 4 , 10 5 , 10 6 for the state space discretization; and the empirical histogram of the P&L of the market maker using the naive strategy. We took intensities that are state dependent. One can see that the larger the size of the grids are, the better the Qknn-estimation of the optimal strategy is.

In Figure 5, we plot the P&L of the market maker following the Qknn-estimated optimal strategy and the naive strategy. We took the same parameters as in Figure 4 to run the simulations except from the terminal time that we set to be equal to T = 10. In this setting, the Qknn-estimated optimal strategy performs much better than the naive strategy, which highlights the fact that the naive strategy is not optimal.

We plotted in Figure 6 the reaction of Qknn when a trend is added in the dynamic of the market. In this example, we took a higher intensity for the sell market order than the one for the buy market order, which creates an artificial positive trend in the dynamic of the price. Observe that Qknn understood correctly that it is better not to sell when the price goes up.

4.4.2.

Case 2: the market maker place orders on the first two limits of the Orders Book. We extend the class of admissible controls to the ones where the market maker places order on the first two limits on the bid and ask sides of the order book. Denote by A2lim the latter. We run simulations to test the Qknn algorithm on A2lim. In figure 7 and figure 8, we plot the empirical distributions of the P&L when the market maker follows the three different strategies:

• Qknn-estimated optimal strategy among those in A2lim (PLOpt2lim).

• Qknn-estimated optimal strategy among those in A1lim (PLOpt1lim).

• naive strategy, i.e. always place orders on the best bid and best ask queues (PL11).

Note that the P&L of the market maker is always better when the class of admissible controls is extended, see Figure 7; but in some numerical tests, the extended set of controls does not seem to improve the P&L. Indeed, we observed that the terminal P&L estimated using Qknn among A2lim and A1lim have the similar empirical distribution in the tests whose results are presented in Figure 8. 11). The intensities λ M , λ L i and λ C i depend on the state of the order book. The P&L of the market maker when following the Qknn-estimated optimal strategy is computed using 10 3 , 10 4 , 10 5 , and 10 6 points for the grids: see respectively Figures 4a, 4b, 4c and4d. The reader can see that the market maker increases their expected terminal wealth (after liquidation) by taking more and more points for the state space discretization. Also, the naive strategy is beaten when the intensities are state dependent.

Figure 5.: Distribution of the P&L when the market maker follows the optimal strategy (dark blue) and the naive strategy (light blue). We took symmetrical and state dependent intensities; and long terminal time: T = 10. Notice that the Qknn strategy does better than the naive strategy when the intensities are state dependent. Figure 6.: Reaction of Qknn to an artificial positive trend in the LOB. The y-axis represents the ratio of orders sent on the ask side by the market maker to their orders sent in the ask and bid sides in the order book. The x-axis is represents the ratio of the size of the best-ask limit to the sum of the sizes of the best-bid and best-ask limits. The blue points are those without trend in the market. The orange points are those with a positive trend in the market. We can see that Qknn took correctly the trend into account in its decisions: for two same order books, Qknn is less willing to sell when the price is expected to increase. Figure 7.: P&L of the market maker who follows optimal strategies and the naive strategy (PL11). Short Terminal Time. asymmetrical intensities for the market order arrivals: the intensity for the buying market order process is taken higher than the one for the selling market order process. The wealth of the market maker is greater when she places orders on the two first limits of each sides of the order book, rather than when she places orders only on the best limits at the bid and ask sides.

Model extension to Hawkes Processes

We consider in this section a market maker who aims at maximizing a function of their terminal wealth, penalizing their inventory at terminal time T in the case where the orders arrivals are driven by Hawkes processes. Let us first present the model with Hawkes processes for the LOB.

Model for the LOB:

We assume that the order book receives limit, cancel, and market orders. We denote by L + (resp. L -) the limit order arrivals process the ask (resp. bid) side; by C + (resp. C -) the cancel order on the ask (resp. bid) side; and by M + (resp. M -) the buy (resp. sell) market order arrivals processes. In this section, the limit orders arrivals are assumed to follow Hawkes processes dynamics, and moreover we assume the kernel to be exponential. The order arrivals are then modeled by a (4K+2)-variate Hawkes process (N t ) with a vector of exogenous intensities λ 0 and exponential kernel φ, i.e. φ ij (t) = α ij βe βt 1 t≥0 . Figure 8.: P&L when following the optimal strategy or the naive strategy (PL11). Long Terminal Time. Symmetrical intensities for the arrival of market orders. 4.10 4 points for the quantization. Notice that the Qknn strategy computed on the extended class of controls, i.e. order placements on the two first limits (StratOpt2lim), performs as well as the one computed on the original class of controls, i.e. order placements on the best-bid and best-ask (StratOpt1lim).

Let α=(α ij ) i,j . We assume: (α) i,j to have spectral radius strictly smaller than 1, which is a sufficient condition to guarantee stationarity of the model and convergence of E[λ u |F t ] as u → +∞, as shown e.g. in [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF].

Note that in the presented model, the following holds:

(Hλ) λ is assumed to be independent of the control.

Denoting by D = 4K + 2 the dimension of (N t ), the m th component of the intensity λ of N t writes, under (Hλ):

λ m t = λ m 0 + D j=1 α mj t 0 e -β(t-s) dN j s , for m = 1, . . . , D,
It is well-known that for this choice of intensity, the couple (N t , λ t ) t≥0 becomes Markovian, see e.g. Lemma 6 in [START_REF] Massoulié | Stability results for a general class of interacting point processes dynamics, and applications[END_REF] for a proof of this result, and moreover we have:

dλ m t = -βD λ m t -λ m 0 dt + D j=1 α mj dN j t , for m = 1, . . . , D,
with given initial conditions: λ m 0 ∈ R * + for m = 1, . . . , D. We can now rewrite the control problem (3.1) in the particular case where the order book is driven by Hawkes processes, there is no running reward, i.e. f = 0, and where the terminal reward G stands for the terminal wealth of the market maker penalized by their inventory. We then consider the following problem in this section:

V (t, λ, z) := sup α∈A E α t,z,λ G Z T , (5.1)
where G(z) denotes the wealth of the market maker when the controlled order book is at state z, plus a term of penalization of their inventory; and where A is the set of the admissible controls, i.e. the predictable decisions taken by the market maker until a terminal time T > 0.

We now present the main result of this section.

Theorem 5.1 V is characterized as the unique solution of the following HJB equation:

                 f (T, z, λ) = G(z), for z ∈ E 0 = ∂f ∂t (t, z, λ) -Dβ D m=1 λ m -λ m 0 ∂f ∂λ m (t, z, λ) +λ m sup a∈Az f t, e a m (z), λ + α m -f t, z, λ , for 0 ≤ t < T, and (t, z, λ) ∈ R + ×E × R * + .
(5.2)

where α m = (α 1m , . . . , α Dm ). Moreover, V admits the following representation:

V (t, z, λ) = sup α∈A ∞ n=0 E α t,z,λ 1 Tn≤T G Z α Tn exp -|λ 0 |(T -T n ) + D m=1 λ m Tn -λ m 0 Dβ e -D j=1 βmj(T -Tn) -1 , (5.3)
where, for n ≥ 0, T n stands for the n th jump time of Z after time t, and (Z α Tn , λ Tn ) ∞ n=0 is as a MDP controlled by α ∈ A; and where E α t,z,λ [.] stands for the expectation conditioned by Z t = z, λ t = λ when the control α is followed.

Remark 5.1 V is characterized in (5.3) as the value function associated to an MDP with infinite horizon, where the running reward reads:

r(t, z, λ) = 1 t≤T G (z) exp -|λ 0 | 1 (T -t) + D m=1 λ m -λ m 0 Dβ e -Dβ(T -t) -1 , where |.| 1 denotes the L 1 R D norm.
Proof: (of Theorem 5.1)

Step 1: Let us check that (5.3) holds, where V is defined as solution of (5.1).

First notice that (λ t , Z t ) t is a PDMDP 1 , i.e. (λ t , Z t ) t is deterministic between two jumping times. We then aim at rewriting the expression of the value function defined in (5.1) as the value function associated to a infinite horizon control problem of the PDMDP (λ t , Z t ) t . To do so, we first notice that by conditioning on the time jumps we get:

V (t, z, λ) = sup α∈A E α t,z,λ G Z α T = sup α∈A E α t,z,λ ∞ n=0 1 Tn≤T <Tn+1 G Z α Tn = sup α∈A ∞ n=0 E α t,z,λ 1 Tn≤T G Z α Tn P T -T n ≤ T n+1 -T n T n , (5.4)
where T n n is the sequence of jump times of N . This process is a jump process with intensity µ s = D m=1 λ m s . Since it holds, conditioned on F Tn :

µ s = D m=1 λ m 0 + λ m Tn -λ m 0 e -Dβ(s-Tn) , for s ∈ [T n , T n+1 ),
then, we have:

P T n+1 -T n ≥ T -T n T n = ∞ T -Tn µ s e -s 0 µudu ds = exp -|λ 0 |(T -T n ) + D m=1 λ m Tn -λ m 0 Dβ e -Dβ(T -Tn) -1 .
(5.5) Plugging (5.5) into (5.4), the value function rewrites:

V (t, z, λ) = sup α∈A ∞ n=0 E α t,z,λ 1 Tn≤T G Z α Tn exp -|λ 0 |(T -T n ) + D m=1 λ m Tn -λ m 0 Dβ e -Dβ(T -Tn) -1 , (5.6)
which completes the step 1. The r.h.s of (5.6) can be seen as the value function of an infinite horizon control problem associated to the PDMDP.

Step 2: Let us show that V is the unique solution to (5.2).

Notice first that the solutions to the following HJB equation

     G(z) = f (T, z, λ) 0 = ∂f ∂t -D m=1 Dβ λ m -λ m 0 ∂f ∂λ m + λ m sup a∈Az f t, e a m (z), λ + α m -f t, z, λ , for 0 ≤ t < T.
are the fixed points of the operator T = T 1 • T 2 where T 1 and T 2 are defined as follows: We now use the characteristic method to rewrite the image of T 1 . Let us take function F , and define f = T 1 (F ). Let us fix t ∈ [0, T ] and λ ∈ (R + ) D , and denote by g the function g(s, z) = f (s, z, λ 1 s , ..., λ D s ) where, for m = 1, . . . , D, s → λ m s is a differentiable function defined on [t, T ] as solution to the following ODE: where T 1 is the first jump time of N larger than t, we denote Z 1 = Z T1 . Equation (5.9) shows that the fixed point of T 1 • T 2 is characterized as the fixed point of the operator T defined for any smooth enough function f by:

T 1 : F → f solution of ∂f ∂t -Dβ D m=1 λ m -λ m 0 ∂f ∂λ m = F (t, z, λ) f (T, z, λ) = G(z),
T (f ) = G(z)e -D m=1 T t λ m s ds + sup a∈Az E a t,λ,z f T 1 , Z 1 , λ T1 + α m ,
where E a t,λ,z [.] stands for the expectation conditioned by the events λ t = λ and Z t = z, when decision a is taken at time t. We recognize here the maximal reward operator of the value function defined in (5.6). Basic theory on PDMDP shows that the maximal reward operator T admits V as unique fixed point, which completes step 2.

Conclusion

In this paper, we solved theoretically and numerically a general market-making problem with different microstructural models of order books, by rewriting the problem as a Markov decision process with infinite horizon. This new representation offers a nice and simple characterization of the optimal strategy of market-making, which is implementable relying for example on some quantization and control randomization ideas as proposed in this paper. Others algorithms, like those based on reinforcement learning (see e.g. Chapter 6.5 of [START_REF] Sutton | Reinforcement Learning: An Introduction, Second[END_REF] for an introduction to (deep) Q-learning, and/or [START_REF] Guéant | Deep reinforcement learning for market making in corporate bonds: beating the curse of dimensionality[END_REF] for its application to market-making), look particularly well-adapted to solve the market-making problem using its MDP representation, especially in the context of high-dimension. The proposed methodology can be adapted to solve theoretically and numerically control problems associated to a general class of controlled point processes.

where we remind that ε proj n := sup a∈A Proj n+1 (F (X n , a, εn )) -F (X n , a, εn ) 2 stands for the average projection error.

Proof. Let us take η > 0, and observe that P Proj n+1 F (X n , a, εn+1 ) -F (X n , a, εn+1 )

2 > η = E where for all x ∈ E and η > 0, B(x, η) denote the ball of center x and radius η. Since x → (1 -x) M is M -Lipschitz, we get by application of Zador's theorem:

P Proj n+1 F (X n , a, εn+1 ) -F (X n , a, εn+1 )

2 > η ≤ M [F ] L [µ] L εn+1 -ε n+1 2 + E 1 -µ B F (X n , a, ε n+1 ), √ η M = M [F ] L [µ] L K 1/d + E 1 -µ B F (X n , a, ε n+1 ), √ η M + O M K 1/d ,
as the number of points for the quantization of the exogenous noise K goes to +∞, and where M stands for the size of the grids Γ n .

Let us introduce A 1 , ..., A N (η) , a cubic partition of Supp(µ), which is bounded under (Hµ), such that for all j = 1, . . . , N (η), A j has diameter η. Also, Notice that there exists c > 0, which only depends on Supp(µ), such as

N (η) ≤ c η d . (A2)
If x ∈ A j , then A j ⊂ B(x, η), therefore: > 0 such that for

E (1 -µ (B(X n , η))) M = N (η) j=1
x, x ∈ E, it holds as M → ∞:

V Q n (x) -V Q n (x ) ≤ V Q n L x -x + O 1 M 1/d . (A7)
Moreover, following bounds holds on V Q n L

, for n = 0, . . . , N :

   V Q N L ≤ [g] L V Q n L ≤ [f ] L + [F ] L V Q n+1 L
, for n = 0, ..., N -1.

(A8)

Proof. Let us show that by induction that V Q N is Lipschitz. First, notice that (A7) holds at ter-

minal time n = N , if one define V Q N L as V Q N L = [g] L . Let us take x, x ∈ E. Assume V Q n+1 (x) -V Q n+1 (x ) ≤ V Q n+1 L |x -x | + O 1 M 1/d
holds for some n = 0, . . . , N -1. Let us show that

V Q n (x) -V Q n (x ) ≤ V Q n L x -x + O 1 M 1/d , where V Q n L
is defined in (A8). Notice that, by the dynamic programming principle and the triangular inequality, it holds:

| V Q n (x) -V Q n (x )| ≤ [f ] L x -x + sup a E a n V Q n+1 Proj n+1 (F (x, a, εn+1 )) -V Q n+1 Proj n+1 F (x , a, εn+1 ) ≤ [f ] L x -x + V Q n+1 L
sup a E Proj n+1 (F (x, a, εn+1 )) -F (x, a, εn+1 )

+ O 1 M 1/d ≤ [f ] L + V Q n+1 L [F ] L x -x + O 1 M 1/d ≤ V Q n L x -x + O 1 M 1/d ,
which completes the proof of (A7).

We now proceed to the proof of Theorem 4.

  for all state (a, b) of the LOB, where |a| := K k=1 -a k and |b| := K k=1 b k .

  Figure1.: Example of market maker's placements and decisions she might make. In this example: the ask-side of the order book is described by a = (7, 5, 4); and the bid-side by b = (-6, -4, -5). The market maker positions are described by ra = (0, 1, -1) and rb = (0, 2, -1...); and the associated ranks vectors are na = (2, 1, -1 . . . ) and nb = (4, 2, -1). We have i0 = 0. After each order arrival, she can send new limit orders (see action on the top right), cancel some positions (see dashed cross on the bottom left), or just keep their orders unchanged.

  1≤i≤K {a i s.t. a i = -1}; and: i0 = argmin 1≤i≤K a i s.t. a i > 0 . The control is the double vector of the positions of the M market maker's orders in the order book.

  only differs from z by its na, nb, and ra, rb components. But |na| ≤ √ M |a| + M and |nb| ≤ √ M |b| + M are bounded by a linear function of (a, b), also |ra| and |rb| are bounded by the universal constant √ M K, so assertion 4 in Definition 3.1 holds.

  ., C K (Γ) are Voronoi tesselations of Γ, i.e., Borel partitions of the Euclidian space (E, |.|) satisfying C (Γ) ⊂ e ∈ E : |e -e | = min j=1,...,K |e -e j | . The discrete law of ε is then characterized by p := P[ε = e ] = P[ε 1 ∈ C (Γ)],= 1, . . . , K.

  .i.d. sequence of r.v. following law µ. Γ n should be seen as a training sampling to estimate the value function V n at time n.

  |r(y; n, x, a)| ≤ r ∞ and |r(y; n, x, a) -r(y; n, x , a)| ≤ [r] L |x -x | and r is defined as follows: P(Z n+1 ∈ O Z n = z, u) = O r(y; n, x, a)dµ(y).and where we denoted by [r] L the Lipschitz constant of r w.r.t. x.

  Figure 2.: Values of the intensities w.r.t. the queue size used in the numerical tests of Section 4.4. The plot is taken from Huang et al. (2015) (see its Figure13). The intensities are estimated by the authors using their data from Alcatel Lucent.

  (a) 10 5 points for the grids. (b) 6000 points for the grids.

Figure 3 .

 3 Figure 3.: Histogram of the P&L when following the Qknn estimated optimal strategy (Opt) and the naive strategy (11). We took symmetrical and constant intensities, and a short terminal time T = 1. Notice that the Qknn strategy looks to improve the P&L by reducing the losses when enough points are taken to build the grids (see Figure3a), and that its performance is worse if less points are taken to build the grids (see Figure3b).

  Figure 4.: Histogram of the P&L when the market maker follow the Qknn estimated optimal strategy (Opt) and the naive strategy (11). The intensities λ M , λ Li and λ C i depend on the state of the order book. The P&L of the market maker when following the Qknn-estimated optimal strategy is computed using 10 3 , 10 4 , 10 5 , and 10 6 points for the grids: see respectively Figures4a, 4b, 4c and 4d. The reader can see that the market maker increases their expected terminal wealth (after liquidation) by taking more and more points for the state space discretization. Also, the naive strategy is beaten when the intensities are state dependent.

  , λ + α m -f t, z, λ .

  m = 1, . . . , D, basic theory on ODE provides existence and uniqueness of a solution to (5.7), which is given by:λ m s = λ m 0 + λ m -λ m 0 e -Dβ(s-t) , for s ∈ [t, T ], and m = 1, . . . , D. then g(t, z) = G(z) -T t F (s, z, λ s ) ds, which finally leads to the following expression of T 1 (F ):T 1 (F ) = f (t, z, λ) = G(z) -T 2 (f ) in (5.8), we get that f is fixed point of T 1 • T 2 if and only if: , λ s + α m ds. Notice ∂f (s, λ s , z)e -e a m (z), λ s + α m , so that: f (t, λ, z) = G(z)e -,z f T 1 , Z 1 , λ T1 + α m ,(5.9)

E 1

 1 |X t,(m) n+1 -F (Xn,a,εn+1)|> √ η X n , εn+1 = E 1 -µ B F (X n , a, εn+1 ), √ η M ,

  ) and (A4) leads toE (1 -µ (B(X n , η))) M ≤ e -1 N (η) M . (A5)Lemma A.3 Under (HF), for n = 0, . . . , N there exists constant V Q n L

  1.Proof. (of Theorem 4.1) Combining inequality|u 1 + u 2 + u 3 | 2 ≤ 3 |u 1 | 2 + |u 2 | 2 + |u 3 | 2 that holds for all u 1 , u 2 , u 3 ∈ R with inequality sup i∈I a i -sup i∈I b i ≤ sup i∈I |a i -b i | thatholds for all families (a i ) i∈I

  z) is borelian. Applying the measurable selection theorem yields to the existence of the maximizer. (seeBäuerle and Rieder (2011) p.352)

	Lemma 3.2 The following holds:

such an order book is sometimes referred to in the literature as an order book governed by a FIFO (First In First Out) rule.

PDMDP stands for Piecewise Deterministic Markov Decision Process, and refers to the the control processes which have deterministic dynamics between (random) jumps.

PDMDP stands for Piecewise Deterministic Markov Decision Process, which is a MDP whose dynamic is deterministic between jumps.

Acknowledgment

We would like to thank the anonymous referees for their valuable comments on the first version of the paper.

Appendix A: Proof of Theorem 4.1 and Corollary 4.1

We divided the proofs of Theorem 4.1 and Corollary 4.1 into the following Lemmas.

Lemma A.1 aims at bounding the projection error. It relies on [START_REF] Gyorfi | A Distribution-Free Theory of Nonparametric Regression[END_REF], see p.93, as well as Zador's theorem, stated in Section B for the sake of completeness.

Lemma A.1 Assume d ≥ 3, and take K = M d+2 points for the optimal quantization of ε n , then it holds under (Hµ) and (HF), as M → +∞,

Let L = 2 µ ∞ stands for the diameter of the support of µ. We then get, as M → +∞,

where c is defined as c :=

, and where we used (A5) and (A2) to go from the second to the third line. It remains to take K = M d+1 points for the optimal quantization of the exogenous noise, and then take square root of equality (A6), in order to derive (A1).

Lemma A.2 Assume d ≥ 3, take K = M d+2 points for the optimal quantization of ε n , and let x ∈ E. Then it holds under (Hµ) and (HF), as M → +∞:

where ε proj n (x), defined as ε proj n (x) := sup a∈A Proj n+1 (F (x, a, εn )) -F (x, a, εn ) 2 , stands for the later-projection error at state x.

Proof. Following the same steps as those used to prove Lemma A.1, we show that:

as K → +∞, and moreover,

holds, which is enough to complete the proof of Lemma A.2. and (b i ) i∈I of reals, and all set I, we have:

2 where E n,Xn stands for the expectation conditioned by the state X n at time n. It holds as M → +∞, using Lemma A.3:

Under (HF), (A9) can then be rewritten as: Assume that E |ε n | 2+η < +∞ for some η > 0. Then, there exists a universal constant C > 0 such that:

Proof. We refer to [START_REF] Graf | Foundations of quantization for probability distributions[END_REF] for a proof of Theorem B.1.


Pascal Lu
email: pascal.lu@outlook.com

Olivier Colliot
email: olivier.colliot@upmc.fr

Multilevel Survival Analysis with Structured Penalties for Imaging Genetics data

Keywords: Survival analysis, Cox model, Imaging genetics, Alzheimer's Disease

Predicting the future occurrence of Alzheimer's disease (AD) in patients with mild cognitive impairment (MCI) is a topic of active research. Many papers have formulated this question as a classification problem: one considers a fixed time of conversion and aims to discriminate between the patients who have converted to AD at that time and those who have not. However, a clinically more relevant question is to predict the date at which a patient to AD. Survival analysis is an adequate statistical framework for such a task.

Multimodal data (imaging and genetic) provide complementary information for the prediction. While imaging data provides an estimate of the current patient's state, genetic variants can be associated to the speed of progression to AD. Although they do not provide the same type of information, most papers in classification or regression put imaging and genetic variables on the same level in order to predict the current or future patient's state.

In this work, we propose a survival model using multimodal data to estimate the conversion date to AD, by considering joint effects between the imaging and genetic modalities. We introduce an adapted penalty in the survival model, the group lasso penalty, over joint groups of genes and brain regions.

The model is evaluated on genetic (single nucleotide polymorphisms) and imaging (anatomical MRI measures) data from the ADNI database, and compared to a standard Cox model.

INTRODUCTION

Early diagnosis of Alzheimer's disease (AD) is an active issue in medical imaging. In Alzheimer's disease, the group of patients with mild cognitive impairment is heterogeneous, some are more likely than others to convert to AD; and therefore giving an accurate diagnosis can be difficult. [START_REF] Razvan | TADPOLE Challenge: Prediction of Longitudinal Evolution in Alzheimer's Disease[END_REF] In this paper, we focus on predicting accurately the conversion to AD given only one time-point.

Most existing approaches to predict the conversion define the problem as a classification task at fixed time. [START_REF] Rathore | A review on neuroimaging-based classification studies and associated feature extraction methods for Alzheimer's disease and its prodromal stages[END_REF][START_REF] Cheng | Domain transfer learning for MCI conversion prediction[END_REF][START_REF] Davatzikos | Prediction of MCI to AD conversion, via MRI, CSF biomarkers, and pattern classification[END_REF] The problem of using classification at a fixed time is that we have to arbitrarily choose a date at which the conversion to AD is observed, we create two groups of homogeneous patients (converting before or after the fixed date), and we do not ensure that conversion in time is monotonous.

Survival models provide a regression framework which directly estimates the conversion date using only one time-point. However, by using survival models, we make several hypotheses: the conversion event is sometimes not observed (because the study ended before conversion), sometimes never occurs (some patients will never convert to AD), and only happens once. For all patients in the study, the starting point is the entry in the study, and the conversion is the time of progression of the disease to AD.

In a single visit, several kind of data can be acquired (clinical data, neuroimages, fluid biomarkers, genetic data. . .). Some data, such as clinical data or neuroimaging data provide a picture of the patient's state at the time they were acquired. [START_REF] Rathore | A review on neuroimaging-based classification studies and associated feature extraction methods for Alzheimer's disease and its prodromal stages[END_REF] On the contrary, others data, such as genetic data, help to identify whether or not a patient could develop AD in the future (for instance, some alleles of APOE increase the risk of developing AD). An issue raised by collecting data from different sources concerns their combination. In the area of imaging genetics, most papers focus either on the association between neuroimaging covariates and genetic data, [START_REF] Liu | A review of multivariate analyses in imaging genetics[END_REF][START_REF] Batmanghelich | Probabilistic modeling of imaging, genetics and diagnosis[END_REF] or on building machine learning predictors for a disease at fixed time using classification (logistic regression, SVM [START_REF] Peng | Structured sparse kernel learning for imaging genetics based Alzheimers Disease Diagnosis[END_REF] ). All these models combine neuroimaging and genetic data by using an additive framework. However, adding the effect of both modalities and putting them on the same level is not optimal, as these modalities do not provide the same type of information. We proposed a multilevel framework [START_REF] Lu | Multilevel Modeling with Structured Penalties for Classification from Imaging Genetics Data[END_REF] for combining imaging and genetic data for classification.

In this paper, we propose a survival model, based on the Cox Proportional Hazard model and using a multilevel framework. Survival models (Cox Proportional Hazard model) have been applied for combining multimodal data [START_REF] Bovelstad | Survival prediction from clinico-genomic models -a comparative study[END_REF] and for predicting the conversion to AD; 15 in both case using an additive framework. Learning the conversion date to AD with modalities taken separately shows that the genetic modality has weaker predictive value than the neuroimaging modality. Instead of summing both genetic and neuroimaging contributions (which could lead to a weaker contribution of the genetic modality in the model), we propose that the parameters, combined with the neuroimaging covariates, could be modulated by the patient's genetic data. This hypothesis leads to a multilevel model where genetic data express themselves through interactions with neuroimaging covariates.

Adding interactions leads to high-dimensional models, and adapted penalties for each modality is essential to avoid overfitting. For instance, SNPs can be grouped by genes [START_REF] Silver | ADNI: Identification of gene pathways implicated in Alzheimer's disease using longitudinal imaging phenotypes with sparse regression[END_REF][START_REF] Lu | Multilevel Modeling with Structured Penalties for Classification from Imaging Genetics Data[END_REF] and a group lasso penalty can be applied on the groups formed by genes. In this paper, we will use the group lasso penalty on the interactions for the parameters coupling (genes, brain region). We use a proximal gradient descent algorithm to learn all the parameters. This model is evaluated on genetic single-nucleotide polymorphisms (SNPs) and neuroimaging data (MRI modality) from the ADNI database (http://adni.loni.usc.edu), and is compared to standard Cox models.

METHODS

Model set-up

We aim to model the time to conversion to Alzheimer's Disease (AD) for MCI patients. For all patients, the starting point is the entry in the study, and the conversion is the time of progression to AD.

Notations For each patient i, we denote T * i his real conversion date from MCI to AD, C i the date of his final visit and T i = min(T * i , C i ) the duration observed in the study. We introduce

δ i = I{T * i ≤ C i } indicating if the conversion has occurred.
We denote x i G the vector of single-polymorphism nucleotides (SNP) counted by number of minor variants, x i I the vector of imaging variables (brain regions), |G| the number of SNPs and |I| the number of imaging variables.

The conversion date T is a continuous random variable with cumulative distribution function

F : t → P{T < t|x i G , x i I } = 1 -exp - t 0 h(u|x i G , x i I )
du where h is the hazard function, representing the instantaneous rate of occurrence of the event.

Multilevel framework We propose the multilevel framework, based on the Cox proportional hazard assumption, defined by h(t|x i G ,

x i I ) = h 0 (t)e β(x i G ) x i I , where β(x i G )

is the parameter vector depending on genetic data x i

G and h 0 is the baseline hazard function describing the risks for individuals whose covariates are null. We make the assumption that β is an affine function depending on genetic data:

β(x i G ) = W x i G + β I , where W ∈ M |G|,|I| (R). Then, h(t|x i G , x i I ) = h 0 (t)e (x i G ) Wx i I +β I x i I
The survival function is given by:

S(t|x i G , x i I ) = (S 0 (t)) (x i G ) Wx i I +β I x i I where S 0 (t) = exp - t 0 h 0 (u)du
The effect of the covariate (x i G , x i I ) on the survival function is to raise it to a power given by the prognostic index PI(x i G ,

x i I ) = e (x i G ) Wx i I +β I x i I .
The genetic modality has a much weaker predictive power compared to imaging or clinical features. And separately taken, the genetic modality provides poor results (see table 1). By combining SNPs and imaging features in that way, we ensure that SNPs will add a significant contribution to the model.

Optimization

Given the dataset {(x i G , x i I , T i , δ i ), i = 1, . . . , N } where the covariates x G , x I and x G x I are centered and normalized, the negative partial log-likelihood is given by:

(W, β I ) = - 1 N N i=1 δ i   (x i G ) Wx i I + β I x i I -log j∈R(Ti) e (x i G ) Wx i I +β I x i I  
where R(T i ) is the set of patients j such that T j ≥ T i .

Penalties As the number of parameters to estimate is much larger than the number of patients, we need to add penalties on W and β I . For imaging parameters β I , we considered the ridge penalty, as Alzheimer's Disease has a diffuse anatomical pattern of alteration. For the matrix W, we start by mapping SNPs to genes G ( ≤ L, where L is the number of genes), and we use a group lasso with overlap penalty, where groups are (genes, imaging covariate). This penalty enforces sparsity between groups and regularity inside the same group. Finally, we add the following penalty to the negative partial log-likelihood:

Ω(W, β I ) = λ |I| i=1 L =1 |G | W G ,i 2 + λ I β I 2 2
where λ > 0, λ I > 0 are the hyperparameters.

The parameters W, β I are obtained by minimizing the quantity (W, β I ) + Ω(W, β I ). The usual approach for dealing with the penalty Ω is to use a proximal gradient descent on the convex set defined by Ω. [START_REF] Hastie | Statistical Learning with Sparsity -The Lasso and Generalizations[END_REF][START_REF] Beck | Gradient-based algorithms with applications to signal recovery problems[END_REF] Algorithm 1: Optimization procedure 1 Input: {(x i G , x i I , T i , δ i ), i = 1, . . . , N }, δ = 0.5, ε 0 = 0.01, η = 10 -5 ; 2 Initialization: W = 0, β I = 0, converged = False 3 while not(converged) do 4 γ = (flatten(W), β I ) and ω = γ -ε∇ ;

5 W G ,i = max   0, 1 - ελ G θ G ω (t) G +i|G| 2   ω G +i|G| for (i, ) ∈ 1, |I| × 1, L ; 6 β I = ω I+|G||I| 1 + 2ελ I (imaging modality) ; 7 if ( + Ω)( W, β I ) > ( + Ω)(W, β I ) then 8 ε = δε 9 else 10 W = W, β I = β I , ε = ε 0 ; 11 end 12 end 13 converged = ( + Ω)(W, β I ) -( + Ω)( W, β I ) ? < η|( + Ω)(W, β I )|
Implementation We flatten the cross-product covariates x G x I and the matrix W and transform them into a vector. We create a vector γ = (flatten(W), β I ) containing the coefficients of W, β I . To recreate W, we just need to unflatten γ. The vector γ is updated using a proximal gradient descent described in Algorithm 1. The stopping criterion for this algorithm is

( + Ω)(W, β I ) -( + Ω)( W, β I ) < η|( + Ω)(W, β I )|

EXPERIMENTS AND RESULTS

Dataset

The ADNI1 GWAS dataset contains 326 MCI patients with 172 MCI patients at baseline who progressed to AD and 154 MCI patients who remained stable (censored data).

Covariates In this dataset, 620,901 SNPs have been genotyped. Based on the 44 first top genes related to AD (from AlzGene, http://www.alzgene.org) and on the Illumina annotation using the Genome build 36.2, we select 1,107 SNPs. For cross-validation purposes, SNPs for whose the variance among the dataset is smaller than 0.01 are removed, leading to 679 SNPs. Regarding the MRI modality, we use the segmentation of FreeSurfer which gives the volume of subcortical regions (44 features) and the average thickness in cortical regions (68 features).

Baseline survival function S 0 The baseline survival function S 0 is computed using the Kaplan-Meier estimate. On figure 1, is shown on the right the Kaplan-Meier estimated baseline survival function S 0 using the distribution of T * and C displayed on the left. The follow-up date is τ hor = 100 months; patients who convert after this date are truncated. The median survival time (the smallest survival time for which the survivor function is less than or equal to 0.5) is 36 months.

Evaluation

Baseline models We compare the multilevel framework to the Cox Proportional Hazard model using one modality or using an additive framework. In this later case, the hazard function for patient i is given by h(t|x i G ,

x i I ) = h 0 (t)e β G x i G +β I x i I .
Metrics We define the three following measures to asses the quality of the prediction:

• the concordance index (or C-index) [START_REF] Steck | On Ranking in Survival Analysis: Bounds on the Concordance Index[END_REF] which checks if the model orders the conversion dates in the same order as the ground truth. As a generalization of AUC, the range of the C-index is [0, 1], but typical values are between 0.55 and 0.7.

• the integrated Brier score, [START_REF] Graf | Assessment and Comparison of Prognostic classification schemes for survival data[END_REF] for uncensored data:

Brier = 1 τ hor τ hor 0 1 N N i=1 δ i (1(T * i > t) -S(t|x G , x I )] 2 dt
The Brier score measures the accuracy of probabilistic predictions. The model performs better when the Brier score is lower.

• the integrated Area Under Curve, defined as iAUC = τ hor 0 AUC(t)f (t)dt where τ hor is the follow-up date, f is the probability density function of T and AUC(t) is the cumulative AUC. [START_REF] Chambless | Estimation of time-dependent area under the ROC curve for long-term risk prediction[END_REF] The cumulative/dynamic AUC plays the same role as the classical Area Under Curve in classification. As for the AUC, the range of the iAUC is [0, 1], and the higher the iAUC is, the more predictive the model is.

Cross validation To determine the hyperparameters, we use a nested cross validation. We perform a 5-fold cross validation, and within each fold, we find the optimal hyperparameters using a 5-fold cross validation on the training set and taking the hyperparameters that maximize the C-index over the inner test set.

The hyperparameters are optimized between {10 -4 , 10 1 show that genetic modality, taken alone, have a much weaker predictive value than the imaging modality. The imaging modality already provides good performances, and adding the genetic modality improves the model performance, but not significantly (in both additive and multilevel frameworks). Adding a penalty on the imaging parameter also increases the performances. Finally, the multilevel model provides slightly better results than the additive model.

Effect of cross-product covariates on conversion

For interpretation purposes, we compute a new reduced matrix W ∈ M |I|,L (R), shown on figure 2, where for each brain region j and gene , W j, = max s∈G |W s,j |.

The matrix on figure 2 shows that some rows and columns have more non-null coefficients than others.

The strongest effects are found for the ventricles, which are enlarged in AD but also in aging and other degenerative diseases, and several medial temporal lobes (MTL) structures (entorhinal cortex, hippocampus, amygdala) which are altered early in AD. It is interesting to note that the ventricles have a strong interaction with all genes except APOE, while the MTL structures have interactions with a more restricted but quite consistent set of genes. Regarding the genes, the strongest effects are found for the single-nucleotide polymorphism rs6503018 (TNK1), rs429358 (APOE) and rs3093662 (TNF).

CONCLUSION

This paper proposes a novel approach to estimate the conversion date to AD for MCI patients, using the Cox Proportional Hazard model, from genetic and neuroimaging data. On the contrary of additive models, the multilevel model captures interactions between genes and brain regions. The use of adapted penalties avoids overfitting by providing a sparse matrix and highlighting brain regions and genes both related to the progression to AD.

Figure 1 :

 1 Figure 1: ADNI1 Dataset: baseline survival function and hasard function

Figure 2 :

 2 Figure 2: Reduced matrix W

Table 1 :

 1 -3 , . . . , 10, 10 2 }. Results for different modalities and methods (mean value across the test folds ± standard deviation) Results on table

	Modality	Method	C-index	Brier score	iAUC
	SNPs only	Cox PH model ( 1 penalty)	0.521 ±0.040	0.166 ±0.009	0.515 ±0.031
	MRI only	Cox PH model (no penalty)	0.636 ±0.034	0.190 ± 0.017	0.636 ±0.050
	MRI only	Cox PH model ( 2 penalty)	0.671 ± 0.022	0.149 ±0.008	0.663 ±0.044
	All	Additive Cox PH model ( 1 penalty)	0.677 ± 0.020	0.148 ± 0.006	0.680 ±0.030
	All	Multilevel model (ours)	0.681 ±0.018	0.147 ± 0.006	0.686 ±0.031