Fabien Guilbaud Exqim

Huyên Pham
email: pham@math.jussieu.fr

Optimal high frequency trading with limit and market orders

Keywords: Market making, limit order book, inventory risk, point process, stochastic control

We propose a framework for studying optimal market making policies in a limit order book (LOB). The bid-ask spread of the LOB is modelled by a Markov chain with finite values, multiple of the tick size, and subordinated by the Poisson process of the tick-time clock. We consider a small agent who continuously submits limit buy/sell orders at best bid/ask quotes, and may also set limit orders at best bid (resp. ask) plus (resp. minus) a tick for getting the execution order priority, which is a crucial issue in high frequency trading. By trading with limit orders, the agent faces an execution risk since her orders are executed only when they meet counterpart market orders, which are modelled by Cox processes with intensities depending on the spread and on her limit prices. By holding non-zero positions on the risky asset, the agent is also subject to the inventory risk related to price volatility. Then the agent can also choose to trade with market orders, and therefore get immediate execution, but at a least favorable price because she has to cross the bid-ask spread.

The objective of the market maker is to maximize her expected utility from revenue over a short term horizon by a tradeoff between limit and market orders, while controlling her inventory position. This is formulated as a mixed regime switching regular/impulse control problem that we characterize in terms of quasi-variational system by dynamic programming methods. In the case of a mean-variance criterion with martingale reference price or when the asset price follows a Levy process and with exponential utility criterion, the dynamic programming system can be reduced to a system of simple equations involving only the inventory and spread variables.

Calibration procedures are derived for estimating the transition matrix and intensity parameters for the spread and for Cox processes modelling the execution of limit orders. Several computational tests are performed both on simulated and real data, and illustrate the impact and profit when considering execution priority in limit orders and market orders.

Introduction

Most of modern equity exchanges are organized as order driven markets. In such type of markets, the price formation exclusively results from operating a limit order book (LOB), an order crossing mechanism where limit orders are accumulated while waiting to be matched with incoming market orders. Any market participant is able to interact with the LOB by posting either market orders or limit orders 1 .

In this context, market making is a class of strategies that consists in simultaneously posting limit orders to buy and sell during the continuous trading session. By doing so, market makers provide counterpart to any incoming market orders: suppose that an investor A wants to sell one share of a given security at time t and that an investor B wants to buy one share of this security at time t ′ > t; if both use market orders, the economic role of the market maker C is to buy the stock as the counterpart of A at time t, and carry until date t ′ when she will sell the stock as a counterpart of B. The revenue that C obtains for providing this service to final investors is the difference between the two quoted prices at ask (limit order to sell) and bid (limit order to buy), also called the market maker's spread. This role was traditionally fulfilled by specialist firms, but, due to widespread adoption of electronic trading systems, any market participant is now able to compete for providing liquidity. Moreover, as pointed out by empirical studies (e.g. [START_REF] Menkveld | High frequency trading and the new-market makers[END_REF], [START_REF] Hendershott | Does algorithmic trading improve liquidity?[END_REF]) and in a recent review [START_REF] Grillet-Aubert | Négociation d'actions: une revue de la littérature à l'usage des régulateurs de marché[END_REF] from AMF, the French regulator, this renewed competition among liquidity providers causes reduced effective market spreads, and therefore reduced indirect costs for final investors.

Empirical studies (e.g. [START_REF] Menkveld | High frequency trading and the new-market makers[END_REF]) also described stylized features of market making strategies. First, market making is typically not directional, in the sense that it does not profit from security price going up or down. Second, market makers keep almost no overnight position, and are unwilling to hold any risky asset at the end of the trading day. Finally, they manage to maintain their inventory, i.e. their position on the risky asset close to zero during the trading day, and often equilibrate their position on several distinct marketplaces, thanks to the use of high-frequency order sending algorithms. Estimations of total annual profit for this class of strategy over all U.S. equity market were around 10 G$ in 2009 [START_REF] Grillet-Aubert | Négociation d'actions: une revue de la littérature à l'usage des régulateurs de marché[END_REF].

Popular models of market making strategies were set up using a risk-reward approach. Three distinct sources of risk are usually identified: the inventory risk, the adverse selection risk and the execution risk. The inventory risk [START_REF] Avellaneda | High frequency trading in a limit order book[END_REF] is comparable to the market risk, i.e. the risk of holding a long or short position on a risky asset. Moreover, due to the uncertain execution of limit orders, market makers only have partial control on their inventory, and therefore the inventory has a stochastic behavior. The adverse selection risk, popular in economic and econometric litterature, is the risk that market price unfavourably deviates, from the market maker point of view, after their quote was taken. This type of risk 1 A market order of size m is an order to buy (sell) m units of the asset being traded at the lowest (highest) available price in the market, its execution is immediate; a limit order of size ℓ at price q is an order to buy (sell) ℓ units of the asset being traded at the specified price q, its execution is uncertain and achieved only when it meets a counterpart market order. Given a security, the best bid (resp. ask ) price is the highest (resp. lowest) price among limit orders to buy (resp. to sell) that are active in the LOB. The spread is the difference, expressed in numéraire per share, of the best ask price and the best bid price, positive during the continuous trading session (see [START_REF] Gould | The limit order book: a survey[END_REF]).

appears naturally in models where the market orders flow contains information about the fundamental asset value (e.g. [START_REF] Frey | Liquidity supply and adverse selection in a pure limit order book market[END_REF]). Finally, the execution risk is the risk that limit orders may not be executed, or be partially executed [START_REF] Kühn | Optimal portfolios of a small investor in a limit order market: a shadow price approach[END_REF]. Indeed, given an incoming market order, the matching algorithm of LOB determines which limit orders are to be executed according to a price/time priority2 , and this structure fundamentally impacts the dynamics of executions.

Some of these risks were studied in previous works. The seminal work [START_REF] Avellaneda | High frequency trading in a limit order book[END_REF] provided a framework to manage inventory risk in a stylized LOB. The market maker objective is to maximize the expected utility of her terminal profit, in the context of limit orders executions occurring at jump times of Poisson processes. This model shows its efficiency to reduce inventory risk, measured via the variance of terminal wealth, against the symmetric strategy. Several extensions and refinement of this setup can be found in recent litterature: [START_REF] Guéant | Dealing with inventory risk[END_REF] provides simplified solution to the backward optimization, an in-depth discussion of its characteristics and an application to the liquidation problem. In [START_REF] Bayraktar | Liquidation in limit order books with controlled intensity[END_REF], the authors develop a closely related model to solve a liquidation problem, and study continuous limit case. The paper [START_REF] Cartea | Modeling Asset Prices for Algorithmic and High Frequency Trading[END_REF] provides a way to include more precise empirical features to this framework by embedding a hidden Markov model for high frequency dynamics of LOB. Some aspects of the execution risk were also studied previously, mainly by considering the trade-off between passive and aggressive execution strategies. In [START_REF] Kühn | Optimal portfolios of a small investor in a limit order market: a shadow price approach[END_REF], the authors solve the Merton's portfolio optimization problem in the case where the investor can choose between market orders or limit orders; in [START_REF] Veraart | Optimal Investment in the Foreign Exchange Market with Proportional Transaction Costs[END_REF], the possibility to use market orders in addition to limit orders is also taken into account, in the context of market making in the foreign exchange market. Yet the relation between execution risk and the microstructure of the LOB, and especially the price/time priority is, so far, poorly investigated.

In this paper we develop a new model to address these three sources of risk. The stock mid-price is driven by a general Markov process, and we model the market spread as a discrete Markov chain that jumps according to a stochastic clock. Therefore, the spread takes discrete values in the price grid, multiple of the tick size. We allow the market maker to trade both via limit orders, which execution is uncertain, and via market orders, which execution is immediate but costly. The market maker can post limit orders at best quote or improve this quote by one tick. In this last case, she hopes to capture market order flow of agents who are not yet ready to trade at the best bid/ask quote. Therefore, she faces a trade off between waiting to be executed at the current best price, or improve this best price, and then be more rapidly executed but at a less favorable price. We model the limit orders strategy as continuous controls, due to the fact that these orders can be updated at high frequency at no cost. On the contrary, we model the market orders strategy as impulse controls that can only occur at discrete dates. We also include fixed, per share or proportional fees or rebates coming with each execution. Execution processes, counting the number of executed limit orders, are modelled as Cox processes with intensity depending both on the market maker's controls and on the market spread. In this context, we optimize the expected utility from profit over a finite time horizon, by choosing optimally between limit and market orders, while controlling the inventory position. We study in detail classical frameworks including mean-variance criterion and exponential utility criterion.

The outline of this paper is as follows. In section 2, we detail the model, and comment its features. We also provide direct calibration methods for all quantities involved in our model. We formulate in Section 3 the optimal market making control problem and derive the associated Hamilton-Jacobi-Bellman quasi variational inequality (HJBQVI) from dynamic programming principle. We show how one can reduce the number of state variables to the inventory and spread for the resolution to this QVI in two standard cases. In the last section 4, we provide some numerical results and empirical performance analysis.

A market-making model 2.1 Mid price and spread process

Let us fix a probability space (Ω, F, P) equipped with a filtration F = (F t ) t≥0 satisfying the usual conditions. It is assumed that all random variables and stochastic processes are defined on the stochastic basis (Ω, F, F, P).

The mid-price of the stock is an exogenous Markov process P , with infinitesimal generator P and state space P. For example, P is a Lévy process (e.g. an arithmetic Brownian motion), or an exponential of Lévy process (e.g. geometric Brownian motion). In the limit order book (LOB) for this stock, we consider a stochastic bid-ask spread resulting from the behaviour of market participants, taking discrete values, which are finite multiple of the tick size δ > 0, and jumping at random times. This is modelled as follows: we first consider the tick time clock associated to a Poisson process (N t ) t with deterministic intensity λ(t), and representing the random times where the buy and sell orders of participants in the market affect the bid-ask spread. We next define a discrete-time stationary Markov chain ( Ŝn ) n∈N , valued in the finite state space S = δI m , I m := {1, . . . , m}, m ∈ N \ {0}, with probability transition matrix (ρ ij ) 1≤i,j≤M , i.e. P[ Ŝn+1 = jδ| Ŝn = iδ] = ρ ij , s.t. ρ ii = 0, independent of N , and representing the random spread in tick time. The spread process (S t ) t in calendar time is then defined as the time-change of Ŝ by N , i.e.

S t = ŜNt , t ≥ 0. (2.1)
Hence, (S t ) t is a continuous time (inhomogeneous) Markov chain with intensity matrix R(t) = (r ij (t)) 1≤i,j≤m , where r ij (t) = λ(t)ρ ij for i = j, and r ii (t) =j =i r ij (t). We assume that S and P are independent. The best-bid and best-ask prices are defined by:

P b t = P t -St 2 , P a t = P t + St 2 .

Trading strategies in the limit order book

We consider an agent (market maker), who trades the stock using either limit orders or market orders. She may submit limit buy (resp. sell) orders specifying the quantity and the price she is willing to pay (resp. receive) per share, but will be executed only when an incoming sell (resp. buy) market order is matching her limit order. Otherwise, she can post market buy (resp. sell) orders for an immediate execution, but, in this case obtain the opposite best quote, i.e. trades at the best-ask (resp. best bid) price, which is less favorable.

Limit orders strategies. The agent may submit at any time limit buy/sell orders at the current best bid/ask prices (and then has to wait an incoming counterpart market order matching her limit), but also control her own bid and ask price quotes by placing buy (resp. sell) orders at a marginal higher (resp. lower) price than the current best bid (resp. ask), i.e. at P b + t := P b t + δ (resp. P a - t := P a t -δ). Such an alternative choice is used in practice by a market maker to capture market orders flow of undecided traders at the best quotes, hence to get priority in the order execution w.r.t. limit order at current best/ask quotes, and can be taken into account in our modelling with discrete spread of tick size δ.

There is then a tradeoff between a larger performance for a quote at the current best bid (resp. ask) price, and a smaller performance for a quote at a higher bid price, but with faster execution. The submission and cancellation of limit orders are for free, as they provide liquidity to the market, and are thus stimulated. Actually, market makers receive some fixed rebate once their limit orders are executed. The agent is assumed to be small in the sense that she does not influence the bid-ask spread. The limit order strategies are then modelled by a continuous time predictable control process:

α make t = (Q b t , Q a t , L b t , L a t ), t ≥ 0,
where L = (L b , L a ) valued in [0, l] 2 , l > 0, represents the size of the limit buy/sell order, and Q = (Q b , Q a ) represent the possible choices of the bid/ask quotes either at best or at marginally improved prices, and valued in

Q = Q b × Q a , with Q b = {Bb, Bb + }, Q a = {Ba, Ba -}:
• Bb: best-bid quote, and Bb + : bid quote at best price plus the tick

• Ba: best-ask quote, and Ba -: ask quote at best price minus the tick Notice that when the spread is equal to one tick δ, a bid quote at best price plus the tick is actually equal to the best ask, and will then be considered as a buy market order. Similarly, an ask quote at best price minus the tick becomes a best bid, and is then viewed as a sell market order. In other words, the limit orders

Q t = (Q b t , Q a t ) take values in Q(S t -), where Q(s) = Q b × Q a when s > δ, Q(s) = {Bb} × {Ba} when s = δ. We shall denote by Q b i = Q b for i > 1, and Q b i = {Bb} for i = 1
, and similarly for Q a i for i ∈ I m . We denote at any time t by π b (Q b t , P t , S t ) and π a (Q a t , P t , S t ) the bid and ask prices of the market maker, where the functions π b (resp. π a ) are defined on Q b × P × S (resp. Q a × P × S) by:

π b (q b , p, s) = p -s 2 , for q b = Bb p -s 2 + δ for q b = Bb + . π a (q a , p, s) = p + s 2 , for q a = Ba p + s 2 -δ for q a = Ba -.
We shall denote by π b i (q b , p) = π b (q b , p, s), π a i (q a , p) = π a (q a , p, s) for s = iδ, i ∈ I m .

Remark 2.1 One can take into account proportional rebates received by the market makers, by considering; [START_REF] Avellaneda | High frequency trading in a limit order book[END_REF], or per share rebates with: π b (q b , p, s) = p -s 2 + δ1 q b =Bb + -ρ, π a (q a , p, s) = p + s 2 -δ1 q a =Ba -+ ρ, for some ρ > 0.

π b (q b , p, s) = (p -s 2 + δ1 q b =Bb + )(1 -ρ), π a (q a , p, s) = (p + s 2 - δ1 q a =Ba -)(1 + ρ), for some ρ ∈ [0,
The limit orders of the agent are executed when they meet incoming counterpart market orders. Let us then consider the arrivals of market buy and market sell orders, which completely match the limit sell and limit buy orders of the small agent, and modelled by independent Cox processes N a and N b . The intensity rate of N a t is given by λ a (Q a t , S t ) where λ a is a deterministic function of the limit quote sell order, and of the spread, satisfying λ a (Ba, s) < λ a (Ba -, s). This natural condition conveys the price/priority in the order execution in the sense that an agent quoting a limit sell order at ask price P a -will be executed before traders at the higher ask price P a , and hence receive more often market buy orders. Typically, one would also expect that λ a is nonincreasing w.r.t. the spread, which means that the larger is the spread, the less often the market buy orders arrive. Likewise, we assume that the intensity rate of N b t is given by λ b (Q b t , S t ) where λ b is a deterministic function of the spread, and λ b (Bb, s) < λ b (Bb + , s). We shall denote by λ a i (q a ) = λ a (q a , s),

λ b i (q b ) = λ b (q b , s) for s = iδ, i ∈ I m .
For a limit order strategy α make = (Q b , Q a , L b , L a ), the cash holdings X and the number of shares Y hold by the agent (also called inventory) follow the dynamics

dY t = L b t dN b t -L a t dN a t , (2.2)
dX t = -π b (Q b t , P t -, S t -)L b t dN b t + π a (Q a t , P t -, S t -)L a t dN a t . (2.3)
Market order strategies. In addition to market making strategies, the investor may place market orders for an immediate execution reducing her inventory. The submissions of market orders, in contrast to limit orders, take liquidity in the market, and are usually subject to fees. We model market order strategies by an impulse control:

α take = (τ n , ζ n ) n≥0 ,
where (τ n ) is an increasing sequence of stopping times representing the market order decision times of the investor, and ζ n , n ≥ 1, are F τn -measurable random variables valued in [-ē, ē], ē > 0, and giving the number of stocks purchased at the best-ask price if ζ n ≥ 0, or selled at the best-bid price if ζ n < 0 at these times. Again, we assumed that the agent is small so that her total market order will be executed immediately at the best bid or best ask price. In other words, we only consider a linear market impact, which does not depend on the order size. When posting a market order strategy, the cash holdings and the inventory jump at times τ n by: In the sequel, we shall denote by A the set of all limit/market order trading strategies α = (α make , α take ).

Y τn = Y τ - n + ζ n , (2.4)
X τn = X τ - n -c(ζ n , P τn , S τn ) (2.

Parameters estimation

In most order-driven markets, available data are made up of Level 1 data that contain transaction prices and quantities at best quotes, and of Level 2 data containing the volume updates for the liquidity offered at the L first order book slices (L usually ranges from 5 to 10). In this section, we propose some direct methods for estimating the intensity of the spread Markov chain, and of the execution point processes, based only on the observation of Level 1 data. This has the advantage of low computational cost, since we do not have to deal with the whole volume of Level 2 data. Yet, we mention some recent work on parameters estimation from the whole order book data [START_REF] Cont | A stochastic model for order book dynamics[END_REF], but involving heavier computations based on integral transforms.

Estimation of spread parameters. Assuming that the spread S is observable, let us define the jump times of the spread process:

θ 0 = 0, θ n+1 = inf {t > θ n : S t = S t-} , ∀n ≥ 1.
From these observable quantities, one can reconstruct the processes:

N t = # {θ j > 0 : θ j ≤ t} , t ≥ 0, Ŝn = S θn , n ≥ 0.
Then, our goal is to estimate the deterministic intensity of the Poisson process (N t ) t , and the transition matrix of the Markov chain ( Ŝn ) n from a path realization with high frequency data of the tick-time clock and spread in tick time over a finite trading time horizon T , typically of one day. From the observations of K samples of Ŝn , n = 1, . . . , K, and since the Markov chain ( Ŝn ) is stationary, we have a consistent estimator (when K goes to infinity) for the transition probability

ρ ij := P[ Ŝn+1 = jδ| Ŝn = iδ] = P[( Ŝn+1 , Ŝn ) = (jδ, iδ)]/P[ Ŝn = iδ]
given by:

ρij = K n=1 1 {( Ŝn+1 , Ŝn)=(jδ,iδ)} K n=1 1 { Ŝn=iδ} (2.6)
For the estimation of the deterministic intensity function λ(t) of the (non)homogeneous Poisson process (N t ), we shall assume in a first approximation a simple natural parametric form. For example, we may assume that λ is constant over a trading day, and more realistically for taking into account intra-day seasonality effects, we consider that the tick time clock intensity jumps e.g. every hour of a trading day. We then assume that λ is in the form:

λ(t) = λ k 1 {t k ≤t<t k+1 }
where (t k ) k is a fixed and known increasing finite sequence of R + with t 0 = 0, and (λ k ) k is an unknown finite sequence of (0, ∞). In other words, the intensity is constant equal to λ k over each period [t k , t k+1 ], and by assuming that the interval length t k+1 -t k is large w.r.t. the intensity λ k (which is the case for high frequency data), we have a consistent estimator of λ k , for all k, and then of λ(t) given by:

λk = N t k+1 -N t k t k+1 -t k . (2.7)
We performed these two estimation procedures (2.6) and (2.7) on SOGN.PA stock on April 18, 2011 between 9:30 and 16:30 in Paris local time. We used tick-by-tick level 1 data provided by Quanthouse, and fed via a OneTick Timeseries database.

In table 1, we display the estimated transition matrix: first row and column indicate the spread value s = iδ and the the cell ij shows ρij . For this stock and at this date, the tick size was δ = 0.005 euros, and we restricted our analysis to the first 6 values of the spread (m = 6) due to the small number of data outside this range. In table 2, we display the estimate for the tick time clock intensity λ(t) by assuming that it is piecewise constant and jumps only at 10:30, 11:30, ..., 16:30 in Paris local time.

In figure 1, we plot the tick time clock intensity by using an affine interpolation, and observed a typical U-pattern. Hence, a further step for the estimation of the intensity could be to specify a parametric form for the intensity function fitting U pattern, e.g. parabolic functions in time, and then use a maximum likelihood method for estimating the parameters.

Estimation of execution parameters. When performing a limit order strategy α make , we suppose that the market maker permanently monitors her execution point processes N a and N b , representing respectively the number of arrivals of market buy and sell orders matching the limit orders for quote ask Q a and quote bid Q b . We also assume that there is no latency so that the observation of the execution processes is not noisy. Therefore, observable variables include the quintuplet:

(N a t , N b t , Q a t , Q b t , S t ) ∈ R + × R + × Q a × Q b × S , t ∈ [0, T ]
Moreover, since N a and N b are assumed to be independent, and both sides of the order book can be estimated using the same procedure, we shall focus on the estimation for the intensity function λ b (q b , s),

q b ∈ Q b = {Bb, Bb + }, s ∈ S = δI m , of the Cox process N b .
The estimation procedure for λ b (q b , s) basically matchs the intuition that one must count the number of executions at bid when the system was in the state (q b , s) and normalize this quantity by the time spent in the state (q b , s). This is justified mathematically as follows. For any (q b , s = iδ) ∈ Q b × S, let us define the point process

N b,q b ,i t = t 0 1 {Q b u =q,S u-=iδ} dN b u , t ≥ 0,
which counts the number of jumps of N b when (Q b , S) was in state (q b , s = iδ). Then, for any nonnegative predictable process (H t ), we have

E ∞ 0 H t dN b,q b ,i t = E ∞ 0 H t 1 {Q b t =q b ,S t-=iδ} dN b t = E ∞ 0 H t 1 {Q b t =q b ,S t-=iδ} λ b (Q b t , S t )dt = E ∞ 0 H t 1 {Q b t =q b ,S t-=iδ} λ b i (q b )dt , (2.8)
where we used in the second equality the fact that λ b (Q b t , S t ) is the intensity of N b . The relation (2.8) means that the point process N b,q b ,i admits for intensity

λ b i (q b )1 {Q b t =q b ,S t-=iδ} . By defining T b,q b ,i t = t 0 1 {Q b u =q,S u-=iδ} du
as the time that (Q b , S) spent in the state (q b , s = iδ), this means equivalently that the process M b,q b ,i t = N b,q b ,i A b,q b ,i t , where A b,q b ,i t = inf{u ≥ 0 : T b,q b ,i u ≥ t} is the càd-làg inverse of T b,q b ,i , is a Poisson process with intensity λ i (q ). By assuming that T b,q b ,i T is large w.r.t. λ i (q b ), which is the case when ( Ŝn ) is irreducible (hence recurrent), and for high-frequency data over [0, T ], we have a consistent estimator of λ b i (q b ) given by:

λb i (q b ) = N b,q b ,i T T b,q b ,i T .
(2.9)

Similarly, we have a consistent estimator for λ a i (q a ) given by:

λa i (q a ) = N a,q a ,i T T a,q a ,i T , (2.10)
where N a,q a ,i T counts the number of executions at ask quote q a and for a spread iδ, and T a,q a ,i T is the time that (Q a , S) spent in the state (q a , s = iδ) over [0, T ]. Let us now illustrate this estimation procedure on real data, with the same market data as above, i.e. tick-by-tick level 1 for SOGN.PA on April 18, 2011, provided by Quanthouse via OneTick timeseries database. Actually, since we did not perform the strategy on this real-world order book, we could not observe the real execution processes N b and N a . We built thus simple proxies Ñ b,q b ,i and Ñ a,q a ,i , for q b = Bb, Bb + , q a = Ba, Ba -, i = 1, . . . , m, based on the following rules. Let us also assume that in addition to (S θn ) n , we observe at jump times θ n of the spread, the volumes (V a θn , V b θn ) offered at the best ask and best bid price in the LOB together with the cumulated market order quantities ϑ BU Y θ n+1 and ϑ SELL θ n+1

arriving between two consecutive jump times θ n and θ n+1 of the spread, respectively at best ask price and best bid price. We finally fix an arbitrarily typical volume V 0 , e.g. V 0 = 100 of our limit orders, and define the proxys Ñ b,q b ,i and Ñ a,q a ,i at times θ n by:

Ñ b,Bb + ,i θ n+1 = Ñ b,Bb + ,i θn + 1 V 0 <ϑ SELL θ n+1 ,S θn =iδ , Ñ b,Bb + ,i 0 = 0 Ñ b,Bb,i θ n+1 = Ñ b,Bb + ,i θn + 1 V 0 +V b θn <ϑ SELL θ n+1
,S θn =iδ , Ñ b,Bb,i

0 = 0 Ñ a,Aa -,i θ n+1 = Ñ a,Aa -,i θn + 1 V 0 <ϑ BU Y θ n+1
,S θn =iδ , Ñ a,Aa -,i

0 = 0 Ñ a,Aa,i θ n+1 = Ñ b,Aa,i θn + 1 V 0 +V a θn <ϑ BU Y θ n+1
,S θn =iδ , Ñ a,Aa,i 0 = 0, together with a proxy for the time spent in spread iδ:

T i θ n+1 = T i θn + (θ n+1 -θ n )1 {S θn =iδ} , T i 0 = 0.
The interpretation of these proxies is the following: we consider the case where the (small) market maker instantaneously updates her quote Q b (resp. Q a ) and volume L b ≤ V 0 (resp. L a ≤ V 0 ) only when the spread changes exogenously, i.e. at dates (θ n ), so that the spread remains constant between her updates, not considering her own quotes. If she chooses to improve best price i.e Q b θn = Bb + (resp. Q a θn = Ba -) she will be in top priority in the LOB and therefore captures all incoming market order flow to sell (resp. buy). Therefore, an unfavourable way for (under)-estimating her number of executions is to increment Ñ b (resp. Ñ a ) only when total traded volume at bid ξ SELL θ n+1 (resp. total volume traded at ask ξ BU Y θ n+1 ) was greater than V 0 . If the market maker chooses to add liquidity to the best prices i.e. Q b θn = Bb (resp. Q a θn = Ba), she will be ranked behind V b θn (resp. V a θn ) in LOB priority queue. Therefore, we increment Ñ b (resp. Ñ a ) only when the total traded volume at bid ϑ SELL θ n+1 (resp. total volume traded at ask ϑ BU Y θ n+1 ) was greater than V 0 + V b θn (resp. V 0 + V a θn ). We then provide a proxy estimate for λ b i (q b ), λ a i (q a ) by:

λb i (q b ) = Ñ b,q b ,i θn T i θn , λa i (q a ) = Ñ a,q a ,i θn T i θn . (2.11)
We performed the estimation procedure (2.11). In table 3 we computed λa i (q a ) and λb i (q b ) for i = 1, . . . , 6, and limit order quotes q b = Bb + , Bb, q a = Ba, Ba -. Due to the lack of data, estimate for large values of the spread are less robust. In this table, each row corresponds to a choice of the spread and each column to a choice of the quotes.

In figure 2 we plotted this estimated intensity as a function of the spread, i.e. s = iδ → λb i (q b ), λa i (q a ) for q b ∈ Q b , and q a ∈ Q a . As one would expect, ( λa i (.), λb i (.)) are decreasing functions of i for the small values of i which matches the intuition that the higher are the (indirect) costs, the smaller is market order flow. Surprisingly, for large values of i this function becomes increasing, which can be due either to an estimation error, caused by the lack of data for this spread range, or a "gaming" effect, in other word liquidity providers increasing their spread when large or autocorrelated market orders come in. i (q a ) and λb i (q b ), expressed in s -1 for SOGN.PA on April 18, 2011, between 9:30 and 16:30 Paris local time. Each row corresponds to a spread value as a multiple of δ = 0.005. Each column corresponds to a limit order quote. 3 Optimal limit/market order strategies

Control problem formulation

Our market model in the previous section is completely determined by the state variables (X, Y, P, S) controlled by the limit/marker order trading strategies α ∈ A.

The objective of the market maker is the following. She wants to maximize over a finite horizon T the profit from her transactions in the LOB, while keeping under control her inventory (usually starting from zero), and getting rid of her inventory at the terminal date:

maximize E U (X T ) -γ T 0 g(Y t )dt (3.1)
over all limit/market order trading strategies α = (α make , α take ) in A such that Y T = 0.

Here U is an increasing reward function, γ is a nonnegative constant, and g is a nonnegative convex function, so that the last integral term T 0 g(Y t )dt penalizes the variations of the inventory. Typical frameworks include the two following cases:

• Mean-quadratic criterion: U (x) = x, γ > 0, g(y) = y 2 .
• Exponential utility maximization:

U (x) = -exp(-ηx), γ = 0.
Let us first remove mathematically the terminal constraint on the inventory: Y T = 0, by introducing the liquidation function L(x, y, p, s) defined on R 2 × P × S by:

L(x, y, p, s) = x -c(-y, p, s) = x + yp -|y| s 2 -ε.
This represents the value that an investor would obtained by liquidating immediately by a market order her inventory position y in stock, given a cash holdings x, a mid-price p and a spread s. Then, problem (3.1) is formulated equivalently as

maximize E U (L(X T , Y T , P T , S T )) -γ T 0 g(Y t )dt (3.2)
over all limit/market order trading strategies α = (α make , α take ) in A. Indeed, the maximal value of problem (3.1) is clearly smaller than the one of problem (3.2) since for any α ∈ A s.t. Y T = 0, we have L(X T , Y T , P T , S T ) = X T . Conversely, given an arbitrary α ∈ A, let us consider the control α ∈ A, coinciding with α up to time T , and to which one add at the terminal date T the market order consisting in liquidating all the inventory Y T . The associated state process ( X, Ỹ , P, S) satisfies: Xt = X t , Ỹt = Y t for t < T , and XT = L(X T , Y T , P T , S T ), ỸT = 0. This shows that the maximal value of problem (3.2) is smaller and then equal to the maximal value of problem (3.1).

We then define the value function for problem (3.2) (or (3.1)):

v(t, z, s) = sup α∈A E t,z,s U (L(Z T , S T )) -γ T t g(Y u )du (3.3) for t ∈ [0, T ], z = (x, y, p) ∈ R 2 × P, s ∈ S.
Here, given α ∈ A, E t,z,s denotes the expectation operator under which the process (Z, S) = (X, Y, P, S) solution to (2.

-valued function v = (v i ) i∈Im defined on [0, T ] × R 2 × P.

Dynamic programming equation

For any q = (q b , q a ) ∈ Q, ℓ = (ℓ b , ℓ a ) ∈ [0, l] 2 , we consider the second-order nonlocal operator:

L q,ℓ ϕ(t, x, y, p, s) = Pϕ(t, x, y, p, s) + R(t)ϕ(t, x, y, p, s)

+ λ b (q b , s) ϕ(t, Γ b (x, y, p, s, q b , ℓ b ), p, s) -ϕ(t, x, y, p, s)
+ λ a (q a , s) ϕ(t, Γ a (x, y, p, s, q a , ℓ a ), p, s) -ϕ(t, x, y, p, s) , (

for (t, x, y, p, s)

∈ [0, T ] × R 2 × P × S, where
R(t)ϕ(t, x, y, p, s) = m j=1 r ij (t) ϕ(t, x, y, p, jδ) -ϕ(t, x, y, p, iδ) , for s = iδ, i ∈ I m , and Γ b (resp. Γ a ) is defined from R 2 × P × S × Q b × R + (resp. R 2 × P × S × Q a × R + into R 2 ) by Γ b (x, y, p, s, q b , ℓ b ) = (x -π b (q b , p, s)ℓ b , y + ℓ b ) Γ a (x, y, p, s, q a , ℓ a ) = (x + π a (q a , p, s)ℓ a , y -ℓ a ).
The first term of L q,ℓ in (3.4) correspond to the infinitesimal generator of the diffusion mid-price process P , the second one is the generator of the continuous-time spread Markov chain S, and the two last terms correspond to the nonlocal operator induced by the jumps of the cash process X and inventory process Y when applying an instantaneous limit order control (Q t , L t ) = (q, ℓ).

Let us also consider the impulse operator associated to market order control, and defined by

Mϕ(t, x, y, p, s) = sup e∈[-ē,ē]
ϕ(t, Γ take (x, y, p, s, e), p, s), where Γ take is the impulse transaction function defined from R 2 × P × S × R into R 2 by: Γ take (x, y, p, s, e) = x -c(e, p, s), y + e ,

The dynamic programming equation (DPE) associated to the control problem (3.3) is the quasi-variational inequality (QVI):

min - ∂v ∂t - sup (q,ℓ)∈Q(s)×[0, l] 2 L q,ℓ v + γg , v -Mv = 0, on [0, T ) × R 2 × P × S(3.5)
together with the terminal condition:

v(T, x, y, p, s) = U (L(x, y, p, s)), ∀(x, y, p) ∈ R 2 × P × S. (3.6)
This is also written explicitly in terms of system of QVIs for the functions

v i , i ∈ I m : min - ∂v i ∂t -Pv i - m j=1 r ij (t)[v j (t, x, y, p) -v i (t, x, y, p)] - sup (q b ,ℓ b )∈Q b i ×[0, l] λ b i (q b )[v i (t, x -π b i (q b , p)ℓ b , y + ℓ b , p) -v i (t, x, y, p)] - sup (q a ,ℓ a )∈Q a i ×[0, l] λ a i (q a )[v i (t, x + π a i (q a , p)ℓ a , y -ℓ a , p) -v i (t, x, y, p)] + γg(y) ; v i (t, x, y, p) -sup e∈[-ē,ē] v i (t,
x -c i (e, p), y + e, p) = 0, for (t, x, y, p) ∈ [0, T ) × R 2 × P, together with the terminal condition:

v i (T, x, y, p) = U (L i (x, y, p)), ∀(x, y, p) ∈ R 2 × P,
where we set L i (x, y, p) = L(x, y, p, iδ). By methods of dynamic programming, one can show by standard arguments that the value function v is the unique viscosity solution to the QVI (3.5)-(3.6) under suitable growth conditions depending on the utility function U and penalty function g. We next focus on some particular cases of interest for reducing remarkably the number of states variables in the dynamic programming equation DPE.

Two special cases

Mean criterion with penalty on inventory

We first consider the case as in [START_REF] Stoikov | Option market making under inventory risk[END_REF] where: U (x) = x, x ∈ R, and (P t ) t is a martingale.

(3.7)

The martingale assumption of the stock price under the historical measure under which the market maker performs her criterion, reflects the idea that she has no information on the future direction of the stock price. Moreover, by starting typically from zero endowment in stock, and by introducing a penalty function on inventory, the market maker wants to keep an inventory that fluctuates around zero. In this case, and exploiting similarly as in [START_REF] Bayraktar | Liquidation in limit order books with controlled intensity[END_REF] the martingale property that Pp = 0, we see that the solution v = (v i ) i∈Im to the dynamic programming system (3.5)-(3.6) is reduced into the form:

v i (t, x, y, p) = x + yp + φ i (t, y), (3.8)
where φ = (φ i ) i∈Im is solution the system of integro-differential equations (IDE):

min - ∂φ i ∂t - m j=1 r ij (t)[φ j (t, y) -φ i (t, y)] - sup (q b ,ℓ b )∈Q b i ×[0, l] λ b i (q b )[φ i (t, y + ℓ b ) -φ i (t, y) + iδ 2 -δ1 q b =Bb + ℓ b ] - sup (q a ,ℓ a )∈Q a i ×[0, l] λ a i (q a )[φ i (t, y -ℓ a ) -φ i (t, y) + iδ 2 -δ1 q a =Ba -ℓ a ] + γg(y) ; φ i (t, y) -sup e∈[-ē,ē] [φ i (t, y + e) - iδ 2 |e| -ε] = 0,
together with the terminal condition:

φ i (T, y) = -|y| iδ 2 -ε,
These one-dimensional IDEs can be solved numerically by a standard finite-difference scheme by discretizing the time derivative of φ, and the grid space in y. The reduced form (3.8) shows that the optimal market making strategies are price independent, and depend only on the level of inventory and of the spread, which is consistent with stylized features in the market. Actually, the IDEs for (φ i ) even show that optimal policies do not depend on the martingale modeling of the stock price.

Exponential utility criterion

We next consider as in [START_REF] Avellaneda | High frequency trading in a limit order book[END_REF] a risk averse market marker:

U (x) = -exp(-ηx), x ∈ R, η > 0, γ = 0, (3.9)
and assume that P follows a Bachelier model:

dP t = bdt + σdW t .
Such price process may take negative values in theory, but at the short-time horizon where high-frequency trading take place, the evolution of an arithmetic Brownian motion looks very similar to a geometric Brownian motion as in the Black-Scholes model. In this case, we see, similarly as in [START_REF] Guéant | Dealing with inventory risk[END_REF], that the solution v = (v i ) i∈Im to the dynamic programming system (3.5)-(3.6) is reduced into the form:

v i (t, x, y, p) = U (x + yp)ϕ i (t, y), (3.10)
where ϕ = (ϕ i ) i∈Im is solution the system of one-dimensional integro-differential equations (IDE):

max - ∂ϕ i ∂t + (bηy - 1 2 σ 2 (ηy) 2 )ϕ i - m j=1 r ij (t)[ϕ j (t, y) -ϕ i (t, y)] - inf (q b ,ℓ b )∈Q b i ×[0, l] λ b i (q b )[exp -η iδ 2 -δ1 q b =Bb + ℓ b ϕ i (t, y + ℓ b ) -ϕ i (t, y)] - inf (q a ,ℓ a )∈Q a i ×[0, l] λ a i (q a )[exp -η iδ 2 -δ1 q a =Ba -ℓ a ϕ i (t, y -ℓ a ) -ϕ i (t, y)] ϕ i (t, y) -inf e∈[-ē,ē]
[exp η|e| iδ 2 + ηε ϕ i (t, y + e)] = 0, together with the terminal condition:

ϕ i (T, y) = exp(η|y| iδ 2 ).
Actually, we notice that the reduced form (3.10) holds more generally when P is a Lévy process, by using the property that in the case: PU (x + yp) = ψ(y)U (x + yp) for some function ψ depending on the characteristics of the generator P of the Lévy process. As in the above mean-variance criterion case, the reduced form (3.10) shows that the optimal market making strategies are price independent, and depend only on the level of inventory and of the spread. However, it depends on the model (typically the volatility) for the stock price.

Computational results

In this section, we provide numerical results in the case of a mean criterion with penalty on inventory, that we will denote within this section by α ⋆ . We used parameters shown in table [START_REF] Cont | A stochastic model for order book dynamics[END_REF] Shape of the optimal policy. The reduced form (3.8) shows that the optimal policy α ⋆ does only depend on time t, inventory y and spread level s. One can represent α ⋆ as a mapping α ⋆ : R + × R × S → A with α ⋆ = (α ⋆,make , α ⋆,take ) thus it divides the space R + × R × S in two zones M and T so that α ⋆ |M = (α ⋆,make , 0) and α ⋆ |T = (0, α ⋆,take ). Therefore we plot the optimal policy in one plane, distinguishing the two zones by a color scale. For the zone M, due to the complex nature of the control, which is made of four scalars, we only represent the prices regimes. Moreover, when using constant tick time intensity λ(t) ≡ λ and in the case where T ≫ 1 λ we can observe on numerical results that the optimal policy is mainly time invariant near date 0; on the contrary, close to the terminal date T the optimal policy has a transitory regime, in the sense that it critically depends on the time variable t. This matches the intuition that to ensure the terminal constraint Y T = 0, the optimal policy tends to get rid of the inventory more aggressively when close to maturity. In figure 3, we plotted a stylized view of the optimal policy, in the plane (y, s), to illustrate this phenomenon.

Benchmarked empirical performance analysis. We made a backtest of the optimal strategy α ⋆ , on simulated data, and benchmarked the results with the three following strategies:

• Optimal strategy without market orders (WoMO), that we denote by α w : this strategy is computed using the same IDEs, but in the case where the investor is not allowed to use market orders, which is equivalent to setting ē = 0.

• Constant strategy, that we denote by α c : this strategy is the symmetric best bid, best ask strategy with constant quantity l0 on both sides, or more precisely α c := (α c,make , 0) with α c,make t ≡ (Bb, Ba, l0 , l0 ).

• Random strategy, that we denote by α r : this strategy consists in choosing randomly the price of the limit orders and using constant quantities on both sides, or more precisely α r := (α r,make , 0) with α r,make t

= (ς b t , ς a t , l0 , l0 ) where (ς b . , ς a . ) is s.t. ∀t ∈ [0; T ] , P(ς b t = Bb) = P(ς b t = Bb + ) = P(ς a t = Ba) = P(ς a t = Ba -) = 1 2
. Our backtest procedure is described as follows. For each strategy α ∈ {α ⋆ , α w , α c , α r }, we simulated N M C paths of the tuple (X α , Y α , P, S, N a,α , N b,α ) on [0, T ], according to eq. (2.1)-(2.2)-(2.3)-(2.4)-(2.5), using a standard Euler scheme with time-step ∆t. Therefore we can compute the empirical mean (resp. empirical standard deviation), that we denote by m(.) (resp. σ(.)), for several quantities shown in table [START_REF] Frey | Liquidity supply and adverse selection in a pure limit order book market[END_REF].

Optimal strategy α ⋆ demonstrates significant improvement of the information ratio IR(X T ) := m(X T )/σ(X T ) compared to the benchmark, which is confirmed by the plot of Even if absolute values of m(X T ) are not representative of what would be the real-world performance of such strategies, these results prove that the different layers of optimization are relevant. Indeed, one can compute the ratios m(X α ⋆ T ) -m(X α c T ) /σ(X α ⋆ T ) = 0.194 and m(X α ⋆ T ) -m(X α w T ) /σ(X α ⋆ T ) = 0.124 that can be interpreted as the performance gain, measured in number of standard deviations, of the optimal strategy α ⋆ compared respectively to the constant strategy α c and the WoMO strategy α w . Another interesting statistics is the surplus profit per trade m(X

α ⋆ T ) -m(X α c T ) / m(N b,α ⋆ T ) + m(N a,α ⋆ T ) + m(N market,α ⋆ T ) =
0.056 euros per trade, recalling that the maximum volume we trade is l = ē = 100. Note that for this last statistics, the profitable effects of the per share rebates ρ are partially neutralized because the number of executions is comparable between α ⋆ and α c ; therefore the surplus profit per trade is mainly due to the revenue obtained from making the spread.

To give a comparison point, typical clearing fee per execution is 0.03 euros on multilateral trading facilities, therefore, in this backtest, the surplus profit per trade was roughly twice the clearing fees. We observe in the synthesis table that the number of executions at bid and ask are symmetric, which is also confirmed by the plots of their empirical distributions in figure [START_REF] Frey | Liquidity supply and adverse selection in a pure limit order book market[END_REF]. This is due to the symmetry in the execution intensities λ b and λ a , which is reflected by the symmetry around y = 0 in the optimal policy. Moreover, notice that the maximum absolute inventory is efficiently kept close to zero in α ⋆ and α w , whereas in α c and α r it can reach much higher values. The maximum absolute inventory is higher in the case of α ⋆ than in the case α w due to the fact that α ⋆ can unwind any position immediately by using market orders, and therefore one may post higher volume for limit orders between two trading at market, profiting from reduced execution risk.

Efficient frontier. An important feature of our algorithm is that the market maker can choose the inventory penalization parameter γ. To illustrate its influence, we varied the inventory penalization γ from 50 to 6.10 -2 , and then build the efficient frontier for both the optimal strategy α ⋆ and for the WoMO strategy α w . Numerical results are provided in table [START_REF] Gould | The limit order book: a survey[END_REF] and a plot of this data is in figure [START_REF] Gould | The limit order book: a survey[END_REF].

We display both the "gross" information ratio IR(X α ⋆ T ) := m(X α ⋆ T )/σ(X α ⋆ T ) and the "net" information ratio NIR(X α ⋆ T ) := m(X α ⋆ T ) -m(X α c T ) /σ(X α ⋆ T ) to have more precise interpretation of the results. Indeed, m(X α T ) seems largely overestimated in this simulated data backtest compared to what would be real-world performance, for all α ∈ {α ⋆ , α w , α c , α r }. Then, to ease interpretation, we assume that α c has zero mean performance in real-world conditions, and therefore offset the mean performance m(X α ⋆ T ) by γ σ(X α ⋆ T ) m(X α ⋆ T ) σ(X α w T ) σ(X α w T ) IR(X α ⋆ T ) NIR(X α ⋆ T ) 50. 000 Observe that highest (net) information ratio is reached for γ ≃ 0.8 for this set of parameters. At this point γ ≃ 0.8, the annualized value of the NIR (obtained by simple extrapolation) is 47, but this simulated data backtest must be completed by a backtest on real data. Qualitatively speaking, the effect of increasing the inventory penalization parameter γ is to increase the zone T where we trade at market. This induces smaller inventory risk, due to the fact that we unwind our position when reaching relatively small values for |y|. This feature can be used to enforce a soft maximum inventory constraint directly by choosing γ.

  5) where c(e, p, s) = ep + |e| s 2 + ε represents the (algebraic) cost function indicating the amount to be paid immediately when passing a market order of size e, given the mid price p, a spread s, and a fixed fee ε > 0. We shall denote by c i (e, p) = c(e, p, s) for s = iδ, i ∈ I m . Remark 2.2 One can also include proportional fees ρ ∈ [0, 1) paid at each market order trading by considering a cost function in the form: c(e, p, s) = (e + ε|e|)p + (|e| + ρe) s 2 + ε, or fixed fees per share with c(e, p, s) = ep + |e|( s 2 + ρ) + ε.

Figure 1 :

 1 Figure 1: Plot of tick time clock intensity estimate for the stock SOGN.PA on April 18, 2011 expressed in second -1 (affine interpolation).

  0.0539 0.1485 0.0718 0.1763 0.010 0.0465 0.0979 0.0520 0.1144 0.015 0.0401 0.0846 0.0419 0.0915 0.020 0.0360 0.0856 0.0409 0.0896 0.025 0.0435 0.1009 0.0452 0.0930 0.030 0.0554 0.1202 0.0614 0.1255

Figure 2 :

 2 Figure 2: Plot of execution intensities estimate as a function of the spread for the stock SOGN.PA on the 18/04/2011, expressed in s -1 (affine interpolation).

Figure 3 :

 3 Figure 3: Stylized shape of the optimal policy sliced in YS.

Figure 4 :

 4 Figure 4: Empirical distribution of terminal wealth X T (spline interpolation).

  (a) N Bid empirical distribution (b) N Ask empirical distribution

Figure 5 :

 5 Figure 5: Empirical distribution of the number of executions on both sides.

Figure 6 :

 6 Figure 6: Efficient frontier plot

Table 1 :

 1 Estimation of the transition matrix (ρ ij ) for the underlying spread of the stock SOGN.PA on April 18, 2011.

	spread 0.005 0.01	0.015 0.02	0.025 0.03
	0.005	0	0.410 0.220 0.160	0.142 0.065
	0.01	0.201 0	0.435 0.192	0.103 0.067
	0.015	0.113 0.221 0	0.4582 0.147 0.059
	0.02	0.070 0.085 0.275 0	0.465 0.102
	0.025	0.068 0.049 0.073 0.363	0	0.446
	0.03	0.077 0.057 0.059 0.112	0.692 0

Table 2 :

 2 Time 10:30 11:30 12:30 13:30 14:30 15:30 16:30 Clock Intensity λ(t) 1.654 0.799 0.516 0.377 0.632 1.305 2.113 Estimation of the tick time clock (hourly basis) for the stock SOGN.PA on April 18, 2011. Tick time clock intensity λ(t) is expressed in second -1 .

	2.5
	2
	1.5
	1
	0.5
	0
	10:30:00 AM	11:30:00 AM	12:30:00 PM	1:30:00 PM	2:30:00 PM	3:30:00 PM	4:30:00 PM

Table 3 :

 3 Estimate of execution intensities λa

  Since the spread takes finite values in S = δI m , it will be convenient to denote for i ∈ I m , by v i (t, z) = v(t, z, iδ). By misuse of notation, we shall often identify the value function with the R m

	1)-(2.2)-
	(2.3)-(2.4)-(2.5), with initial state (Z t -, S t -) = (z, s), is taken. Problem (3.3) is a mixed
	regular/impulse control problem in a regime switching jump-diffusion model, that we shall
	study by dynamic programming methods.

  together with transition probabilities (ρ ij ) 1≤i,j≤M calibrated in table[START_REF] Avellaneda | High frequency trading in a limit order book[END_REF] and execution intensities calibrated in table[START_REF] Cartea | Modeling Asset Prices for Algorithmic and High Frequency Trading[END_REF], slightly modified to make the bid and ask sides symmetric.

	Parameter Signification δ Tick size ρ Per share rebate ǫ Per share fee ǫ0 Fixed fee λ(t) Tick time intensity ≡ 1s -1 Value 0.005 0.0008 0.0012 10 -6	Parameter Signification U (x) Utility function γ Inventory penalization 5 Value x l Max. volume make 100 ē Max. volume take 100
	(a) Market parameters		(b) Optimization parameters
	Parameter Signification T Length in seconds ymin Lower bound shares ymax Upper bound shares n Number of time steps 100 Value 300 s -1000 1000 m Number of spreads 6	Parameter Signification N M C Number of paths for MC simul. 10 5 Value ∆t Euler scheme time step 0.3 s l0 B/A qty for bench. strat. 100 x0 Initial cash 0 y0 Initial shares 0 p0 Initial price 45
	(c) Discretization/localization parameters		(d) Backtest parameters

Table 4 :

 4 Parameters

  optimal α ⋆ WoMO α w constant α c random α r

	Terminal wealth	m(XT )/σ(XT )	2.117	1.999	0.472	0.376
		m(XT )			26.759	25.19	24.314	24.022
		σ(XT )			12.634	12.599	51.482	63.849
	Num. of exec. at bid	m(N b T )			18.770	18.766	13.758	21.545
		σ(N b T )			3.660	3.581	3.682	4.591
	Num. of exec. at ask	m(N a T )			18.770	18.769	13.76	21.543
		σ(N a T )			3.666	3.573	3.692	4.602
	Num. of exec. at market m(N market T	)	6.336	0	0	0
		σ(N market T	)	2.457	0	0	0
	Maximum Inventory	m(sup s∈[0;T ] |Ys|) 241.019	176.204	607.913	772.361
		σ(sup s∈[0;T ] |Ys|)	53.452	23.675	272.631	337.403

Table 5 :

 5 Performance analysis: synthesis of benchmarked backtest (10 5 simulations).

the whole empirical distribution of X T (see figure (

4

)).

Table 6 :

 6 Efficient frontier data the constant -m(X α c T ) when computing the NIR. This has simple visual interpretation as shown in figure[START_REF] Gould | The limit order book: a survey[END_REF].

		5.283	12.448	4.064	9.165	2.356	-2.246
	25.000 7.562	18.421	7.210	16.466	2.436	-0.779
	12.500 9.812	22.984	9.531	20.971	2.343	-0.135
	6.250	11.852	25.932	11.749	24.232	2.188	0.136
	3.125	14.546	28.153	14.485	26.752	1.935	0.263
	1.563	15.819	28.901	16.830	28.234	1.827	0.289
	0.781	19.088	29.952	19.593	29.145	1.569	0.295
	0.391	20.898	30.372	20.927	29.728	1.453	0.289
	0.195	23.342	30.811	23.247	30.076	1.320	0.278
	0.098	25.232	30.901	24.075	30.236	1.225	0.261
	0.049	26.495	31.020	24.668	30.434	1.171	0.253
	0.024	27.124	30.901	25.060	30.393	1.139	0.242
	0.012	27.697	31.053	25.246	30.498	1.121	0.243
	0.006	28.065	30.998	25.457	30.434	1.105	0.238

A different type of LOB operates under pro-rata priority, e.g. for some futures on interest rates. In this paper, we do not consider this case and focus on the main mechanism used in equity market.


CL Marc Arnaudon

Ément Dombry

AND Anthony Phan

L E Yang

STOCHASTIC ALGORITHMS FOR COMPUTING MEANS OF PROBABILITY MEASURES

Consider a probability measure µ supported by a regular geodesic ball in a manifold. For any p ≥ 1 we define a stochastic algorithm which converges almost surely to the p-mean ep of µ. Assuming furthermore that the functional to minimize is regular around ep, we prove that a natural renormalization of the inhomogeneous Markov chain converges in law into an inhomogeneous diffusion process. We give an explicit expression of this process, as well as its local characteristic.

Introduction

The geometric barycenter of a set of points is the point which minimizes the sum of the distances at the power 2 to these points. It is the most common estimator is statistics, however it is sensitive to outliers, and it is natural to replace power 2 by p for some p ∈ [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF][START_REF] Arnaudon | Barycenters of measures transported by stochastic flows[END_REF], which leads to the definition of p-mean. When p = 1, the minimizer is the median of the set of points, very often used in robust statistics. In many applications, p-means with some p ∈ (1, 2) give the best compromise.

The Fermat-Weber problem concerns finding the median e 1 of a set of points in an Euclidean space. Numerous authors worked out algorithms for computing e 1 . The first algorithm was proposed by Weiszfeld in [START_REF] Weiszfeld | Sur le point pour lequel la somme des distances de n points donnés est minimum[END_REF]. It has been extended to sufficiently small domains in Riemannian manifolds with nonnegative curvature by Fletcher and al in [START_REF] Fletcher | The geometric median on Riemannian manifolds with application to robust atlas estimation[END_REF]. A complete generalization to manifolds with positive or negative curvature, including existence and uniqueness results (under some convexity conditions in positive curvature), has been given by one of the authors in [START_REF] Yang | Riemannian median and its estimation[END_REF].

The Riemannian barycenter or Karcher mean of a set of points in a manifold or more generally of a probability measure has been extensively studied, see e.g. [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF], [START_REF] Kendall | Probability, convexity and harmonic maps with small image I: uniqueness and fine existence[END_REF], [START_REF] Kendall | Convexity and the hemisphere[END_REF], [START_REF] Emery | Sur le barycentre d'une probabilité dans une variété[END_REF], [START_REF] Picard | Barycentres et martingales dans les variétés[END_REF], [START_REF] Arnaudon | Barycenters of measures transported by stochastic flows[END_REF], where questions of existence, uniqueness, stability, relation with martingales in manifolds, behaviour when measures are pushed by stochastic flows have been considered. The Riemannian barycenter corresponds to p = 2 in the above description. Computation of Riemannian barycenters by gradient descent has been performed by Le in [START_REF] Le | Estimation of Riemannian barycentres[END_REF].

In [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF] Afsari proved existence and uniqueness of p-means, p ≥ 1 on geodesic balls with radius r < 1 2 min inj(M ), π 2α if p ∈ [1, 2), and r < 1 2 min inj(M ), π α if p ≥ 2. Here inj(M ) is the injectivity radius of M and α > 0 is such that the sectional curvatures in M are bounded above by α 2 . The point is that in the case p ≥ 2, the functional to minimize is not convex any more, which makes the situation much more difficult to handle.

In this paper, under the assumptions of [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF] we provide in Theorem 2.3 stochastic algorithms which converge almost surely to p-means in manifolds, which are easier to implement than gradient descent algorithm since computing the gradient of the function to minimize is not needed. The idea is at each step to go in the direction of a point of the support of µ. The point is chosen at random according to µ and the size of the step is a well chosen function of the distance to the point, p and the number of the step. For general convergence results on recursive stochastic algorithms, see [START_REF] Ljung | Analysis of recursive stochastic algorithms[END_REF] Theorem 1. However they do not cover the manifold case and nonlinearity of geodesics. Here we give a proof using martingale convergence theorem, and the main point consists in determining and estimating all the geometric quantities, checking that under our curvature conditions all the convergence assumptions are fulfilled, since our processes live in manifolds. See also [START_REF] Benveniste | Analysis of stochastic approximation schemes with discontinuous and dependent forcing terms with applications to data communication algorithm[END_REF] for convergence in probability of recursive algorithms.

The speed of convergence is studied, and in theorem 2.6 we prove that the renormalized inhomogeneous Markov chain of Theorem 2.3 converges in law to an inhomogeneous diffusion process. This is an invariance principle type result, see e.g. [START_REF] Khas | On stochastic processes defined by differential equations with a small parameter[END_REF], [START_REF] Mcleish | Dependent central limit theorems and invariance principles[END_REF], [START_REF] Berger | An almost sure invariance principle for stochastic approximation procedures in linear filtering theory[END_REF], [START_REF] Gouëzel | Almost sure invariance principle for dynamical systems by spectral methods The annals of Probability[END_REF] for related works. Here again the main point is to obtain the characteristics of the limiting process from the curvature conditions, the conditions on the support of the mesure and estimates on Jacobi fields. Moreover we consider convergence in law for the Skorohod topology, and the limit depends in a crucial way on the decreasing steps of the algorithms.

Results

p-means in regular geodesic balls.

Let M be a Riemannian manifold with pinched sectional curvatures. Let α, β > 0 such that α 2 is a positive upper bound for sectional curvatures on M , and -β 2 is a negative lower bound for sectional curvatures on M . Denote by ρ the Riemannian distance on M .

In M consider a geodesic ball B(a, r) with a ∈ M . Let µ be a probability measure with support included in a compact convex subset K µ of B(a, r). Fix p ∈ [1, ∞). We will always make the following assumptions on (r, p, µ):

Assumption 2.1. The support of µ is not reduced to one point. Either p > 1 or the support of µ is not contained in a line, and the radius r satisfies

(2.1) r < r α,p with r α,p = 1 2 min inj(M ), π 2α if p ∈ [1, 2) r α,p = 1 2 min inj(M ), π α if p ∈ [2, ∞) Note that B(a, r) is convex if r < 1 2 min inj(M ), π α .
Under assumption 2.1, it has been proved in [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF] (Theorem 2.1) that the function

H p : M → R + x → M ρ p (x, y)µ(dy) (2.2)
has a unique minimizer e p in M , the p-mean of µ, and moreover e p ∈ B(a, r). If p = 1, e 1 is the median of µ.

It is easily checked that if p ∈ [1, 2), then H p is strictly convex on B(a, r). On the other hand, if p ≥ 2 then H p is of class C 2 on B(a, r). Proposition 2.2. Let K be a convex subset of B(a, r) containing the support of µ. Then there exists C p,µ,K > 0 such that for all x ∈ K,

(2.3) H p (x) -H p (e p ) ≥ C p,µ,K 2 ρ(x, e p ) 2 .
Moreover if p ≥ 2 then we can choose C p,µ,K so that for all x ∈ K,

(2.4) grad x H p 2 ≥ C p,µ,K (H p (x) -H p (e p )) .
In the sequel, we fix

(2.5) K = B(a, r -ε) with ε = ρ(K µ , B(a, r) c ) 2 .
We now state our main result: we define a stochastic gradient algorithm (X k ) k≥0 to approximate the p-mean e p and prove its convergence.

Theorem 2.3. Let (P k ) k≥1 be a sequence of independent B(a, r)-valued random variables, with law µ. Let (t k ) k≥1 be a sequence of positive numbers satisfying

(2.6) ∀k ≥ 1, t k ≤ min 1 C p,µ,K , ρ(K µ , B(a, r) c ) 2p(2r) p-1 , (2.7)
∞ k=1 t k = +∞ and ∞ k=1 t 2 k < ∞.
Letting x 0 ∈ K, define inductively the random walk (X k ) k≥0 by (2.8)

X 0 = x 0 and for k ≥ 0 X k+1 = exp X k -t k+1 grad X k F p (•, P k+1 )
where F p (x, y) = ρ p (x, y), with the convention grad x F p (•, x) = 0. The random walk (X k ) k≥1 converges in L 2 and almost surely to e p .

In the following example, we focus on the case M = R d and p = 2 where drastic simplifications occur.

Example 2.4. In the case when M = R d and µ is a compactly supported probability measure on R d , the stochastic gradient algorithm (2.8) simplifies into

X 0 = x 0 and for k ≥ 0 X k+1 = X k -t k+1 grad X k F p (•, P k+1 ).
If furthermore p = 2, clearly e 2 = E[P 1 ] and grad x F p (•, y) = 2(x -y), so that the linear relation X k+1 = (1 -2t k+1 )X k + 2t k+1 P k+1 , k ≥ 0 holds true and an easy induction proves that (2.9)

X k = x 0 k-1 j=0 (1 -2t k-j ) + 2 k-1 j=0 P k-j t k-j j-1 ℓ=0 (1 -2t k-ℓ ), k ≥ 1. Now, taking t k = 1 2k , we have k-1 j=0 (1 -2t k-j ) = 0 and j-1 ℓ=0 (1 -2t k-ℓ ) = k -j k so that X k = k-1 j=0 P k-j 1 k = 1 k k j=1 P j .
The stochastic gradient algorithm estimating the mean e 2 of µ is given by the empirical mean of a growing sample of independent random variables with distribution µ. In this simple case, the result of Theorem 2.3 is nothing but the strong law of large numbers. Moreover, fluctuations around the mean are given by the central limit theorem and Donsker's theorem.

2.2.

Fluctuations of the stochastic gradient algorithm. The notations are the same as in the beginning of section 2.1. We still make assumption 2.1. Let us define K and ε as in (2.5) and let (2.10)

δ 1 = min 1 C p,µ,K , ρ(K µ , B(a, r) c ) 2p(2r) p-1 .
We consider the time inhomogeneous M -valued Markov chain (2.8) in the particular case when (2.11)

t k = min δ k , δ 1 , k ≥ 1
for some δ > 0. The particular sequence (t k ) k≥1 defined by (2.11) satisfies (2.6) and (2.7), so Theorem 2.3 holds true and the stochastic gradient algorithm (X k ) k≥0 converges a.s. and in L 2 to the p-mean e p .

In order to study the fluctuations around the p-mean e p , we define for n ≥ 1 the rescaled

T ep M -valued Markov chain (Y n k ) k≥0 by (2.12) Y n k = k √ n exp -1 ep X k .
We will prove convergence of the sequence of process (Y n [nt] ) t≥0 to a non-homogeneous diffusion process. The limit process is defined in the following proposition:

Proposition 2.5. Assume that H p is C 2 in a neighborhood of e p , and that δ > C -1 p,µ,K . Define Γ = E grad ep F p (•, P 1 ) ⊗ grad ep F p (•, P 1 )
and G δ (t) the generator

(2.13) G δ (t)f (y) := d y f, t -1 (y -δ∇dH p (y, •) ♯ ) + δ 2 2 Hess y f (Γ)
where ∇dH p (y, •) ♯ denotes the dual vector of the linear form ∇dH p (y, •).

There exists a unique inhomogeneous diffusion process (y δ (t)) t>0 on T ep M with generator G δ (t) and converging in probability to 0 as t → 0 + . The process y δ is continuous, converges a.s. to 0 as t → 0 + and has the following integral representation:

(2.14) y δ (t) = d i=1 t 1-δλi t 0 s δλi-1 δσ dB s , e i e i , t ≥ 0,
where B t is a standard Brownian motion on T ep M , σ ∈ End(T ep M ) satisfies σσ * = Γ, (e i ) 1≤i≤d is an orthonormal basis diagonalizing the symmetric bilinear form ∇dH p (e p ) and (λ i ) 1≤i≤d are the associated eigenvalues.

Note that the integral representation (2.14) implies that y δ is the centered Gaussian process with covariance

(2.15) E y i δ (t 1 )y j δ (t 2 ) = δ 2 Γ(e * i ⊗ e * j ) δ(λ i + λ j ) -1 t 1-δλi 1 t 1-δλj 2 (t 1 ∧ t 2 ) δ(λi+λj )-1 ,
where y i δ (t) = y δ (t), e i , 1 ≤ i, j ≤ d and t 1 , t 2 ≥ 0. Our main result on the fluctuations of the stochastic gradient algorithm is the following: Theorem 2.6. Assume that either e p does not belong to the support of µ or p ≥ 2.

Assume furthermore that δ > C -1 p,µ,K . The sequence of processes Y n [nt] t≥0 weakly converges in D((0, ∞), T ep M ) to y δ . Remark 2.7.
The assumption on e p implies that H p is of class C 2 in a neighbourhood of e p . In the case p > 1, in the "generic" situation for applications, µ is a discrete measure and e p does not belong to its support. For p = 1 one has to be more careful since if µ is equidistributed in a random set of points, then with positive probability e 1 belongs to the support of µ.

Remark 2.8. From section 2.1 we know that, when p ∈ (1, 2], the constant

C p,µ,K = p(2r) p-2 (min (p -1, 2αr cot (2αr)))
is explicit. The constraint δ > C -1 p,µ,K can easily be checked in this case.

Remark 2.9. In the case

M = R d , Y n k = k √ n (X k -e p )
and the tangent space T ep M is identified to R d . Theorem 2.6 holds and, in particular, when t = 1, we obtain a central limit Theorem:

√ n(X n -e p ) converges as n → ∞ to a centered Gaussian d-variate distribution (with covariance structure given by (2.15) with t 1 = t 2 = 1). This is a central limit theorem: the fluctuations of the stochastic gradient algorithm are of scale n -1/2 and asymptotically Gaussian.

Proofs

For simplicity, let us write shortly e = e p in the proofs.

Proof of Proposition 2.2.

For p = 1 this is a direct consequence of [START_REF] Yang | Riemannian median and its estimation[END_REF] Theorem 3.7.

Next we consider the case p ∈ (1, 2). Let K ⊂ B(a, r) be a compact convex set containing the support of µ. Let x ∈ K\{e}, t = ρ(e, x), u ∈ T e M the unit vector such that exp e (ρ(e, x)u) = x, and γ u the geodesic with initial speed u : γu (0) = u. For y ∈ K, letting h y (s) = ρ(γ u (s), y) p , s ∈ [0, t], we have since p > 1

h y (t) = h y (0) + th ′ y (0) + t 0 (t -s)h ′′ y (s) ds
with the convention h ′′ y (s) = 0 when γ u (s) = y. Indeed, if y ∈ γ([0, t]) then h y is smooth, and if y ∈ γ([0, t]), say y = γ(s 0 ) then h y (s) = |s -s 0 | p and the formula can easily be checked.

By standard calculation,

h ′′ y (s) ≥ pρ(γ u (s), y) p-2 × (p -1) γu (s) T (y) 2 + γu (s) N (y) 2 αρ(γ u (s), y) cot (αρ(γ u (s), y)) (3.1)
with γu (s) T (y) (resp. γu (s) N (y) ) the tangential (resp. the normal) part of γu (s)

with respect to n(γ u (s), y) = 1 ρ(γ u (s), y) exp -1 γu(s) (y): γu (s) T (y) = γu (s), n(γ u (s), y) n(γ u (s), y), γu (s) N (y) = γu (s) -γu (s) T (y) .
From this we get

h ′′ y (s) ≥ pρ(γ u (s), y) p-2 (min (p -1, 2αr cot (2αr))). (3.2) Now H p (γ u (t ′ )) = B(a,r) h y (γ u (t ′ )) µ(dy) = B(a,r) h y (0) µ(dy) + t ′ B(a,r) h ′ y (0) µ(dy) + t ′ 0 (t ′ -s) B(a,r) h y (s) ′′ µ(dy) ds
and H p (γ u (t ′ )) attains its minimum at t ′ = 0, so

B(a,r)
h ′ y (0) µ(dy) = 0 and we have

H p (x) = H p (γ u (t)) = H p (e) + t 0 (t -s) B(a,r) h y (s) ′′ µ(dy) ds.
Using Equation (3.2) we get

H p (x) ≥ H p (e) + t 0 (t -s) B(a,r) pρ(γ u (s), y) p-2 (min (p -1, 2αr cot (2αr))) µ(dy) ds. (3.3) Since p ≤ 2 we have ρ(γ u (s), y) p-2 ≥ (2r) p-2 and H p (x) ≥ H p (e) + t 2 2 p(2r) p-2 (min (p -1, 2αr cot (2αr))) . (3.4) So letting C p,µ,K = p(2r) p-2 (min (p -1, 2αr cot (2αr)))
we obtain

(3.5) H p (x) ≥ H p (e) + C p,µ,K ρ(e, x) 2 2 .
To finish let us consider the case p ≥ 2.

In the proof of [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF] Theorem 2.1, it is shown that e is the only zero of the maps x → grad x H p and x → H p (x) -H p (e), and that ∇dH p (e) is strictly positive. This implies that (2.3) and (2.4) hold on some neighbourhood B(e, ε) of e. By compactness and the fact that H p -H p (e) and grad H p do not vanish on K\B(e, ε) and H p -H p (e) is bounded, possibly modifying the constant C p,µ,K , (2.3) and (2.4) also holds on K\B(e, ε).

Proof of Theorem 2.3.

Note that, for x = y,

grad x F (•, y) = pρ p-1 (x, y) -exp -1 x (y) ρ(x, y) = -pρ p-1 (x, y)n(x, y),
whith n(x, y) := exp -1 x (y) ρ(x, y) a unit vector. So, with the condition (2.6) on t k , the random walk (X k ) k≥0 cannot exit K: if X k ∈ K then there are two possibilities for X k+1 :

• either X k+1 is in the geodesic between X k and P k+1 and belongs to K by convexity of K; • or X k+1 is after P k+1 , but since

t k+1 grad X k F p (•, P k+1 ) = t k+1 pρ p-1 (X k , P k+1 ) ≤ ρ(K µ , B(a, r) c ) 2p(2r) p-1 pρ p-1 (X k , P k+1 ) ≤ ρ(K µ , B(a, r) c ) 2 ,
we have in this case

ρ(P k+1 , X k+1 ) ≤ ρ(K µ , B(a, r) c ) 2 which implies that X k+1 ∈ K. First consider the case p ∈ [1, 2). For k ≥ 0 let t → E(t) := 1 2 ρ 2 (e, γ(t)) , γ(t) t∈[0,t k+1 ] the geodesic satisfying γ(0) = -grad X k F p (•, P k+1 ). We have for all t ∈ [0, t k+1 ] (3.6) E ′′ (t) ≤ C(β, r, p) := p 2 (2r) 2p-1 β cotanh(2βr)
(see e.g. [START_REF] Yang | Riemannian median and its estimation[END_REF]). By Taylor formula,

ρ(X k+1 , e) 2 = 2E(t k+1 ) = 2E(0) + 2t k+1 E ′ (0) + t 2 k+1 E ′′ (t) for some t ∈ [0, t k+1 ] ≤ ρ(X k , e) 2 + 2t k+1 grad X k F p (•, P k+1 ), exp -1 X k (e) + t 2 k+1 C(β, r, p).
Now from the convexity of x → F p (x, y) we have for all x, y ∈ B(a, r)

(3.7) F p (e, y) -F p (x, y) ≥ grad x F p (•, y), exp -1
x (e) . This applied with x = X k , y = P k+1 yields (3.8) We want to prove that a = 0. We already proved that

ρ(X k+1 , e) 2 ≤ ρ(X k , e) 2 -2t k+1 (F p (X k , P k+1 ) -F p (e, P k+1 )) + C(β, r, p)t 2 k+1 . Letting for k ≥ 0 F k = σ(X ℓ , 0 ≤ ℓ ≤ k), we get E ρ(X k+1 , e) 2 |F k ≤ ρ(X k , e) 2 -2t k+1 B(a,r) (F p (X k , y) -F p (e, y)) µ(dy) + C(β, r, p)t 2 k+1 = ρ(X k , e) 2 -2t k+1 (H p (X k ) -H p (e)) + C(β, r, p)t 2 k+1 ≤ ρ(X k , e) 2 + C(β, r, p)t 2
(3.11) E ρ(X k+1 , e) 2 |F k ≤ ρ(X k , e) 2 -2t k+1 (H p (X k ) -H p (e)) + C(β, r, p)t 2 k+1 .
Taking the expectation and using Proposition 2.2, we obtain

(3.12) E ρ(X k+1 , e) 2 ≤ E ρ(X k , e) 2 -t k+1 C p,µ,K E ρ(X k , e) 2 + C(β, r, p)t 2
k+1 . An easy induction proves that for ℓ ≥ 1,

(3.13) E ρ(X k+ℓ , e) 2 ≤ ℓ j=1 (1 -C p,µ,K t k+j )E ρ(X k , e) 2 + C(β, r, p) ℓ j=1 t 2 k+j .
Letting ℓ → ∞ and using the fact that

∞ j=1 t k+j = ∞ which implies ∞ j=1 (1 -C p,µ,K t k+j ) = 0, we get (3.14) a ≤ C(β, r, p) ∞ j=1 t 2 k+j .
Finally using ∞ j=1 t 2 j < ∞ we obtain that lim k→∞ ∞ j=1 t 2 k+j = 0, so a = 0. This proves L 2 and almost sure convergence.

Next (see e.g. [START_REF] Yang | Riemannian median and its estimation[END_REF]). By Taylor formula,

H p (X k+1 ) = E p (t k+1 ) = E p (0) + t k+1 E ′ p (0) + t 2 k+1 2 E ′′ p (t) for some t ∈ [0, t k+1 ] ≤ H p (X k ) + t k+1 d X k H p , grad X k F p (•, P k+1 ) + t 2 k+1 C(β, r, p).
We get

E [H p (X k+1 )|F k ] ≤ H p (X k ) -t k+1 d X k H p , B(a,r) grad X k F p (•, y)µ(dy) + C(β, r, p)t 2 k+1 = H p (X k ) -t k+1 d X k H p , grad X k H p (•) + C(β, r, p)t 2 k+1 = H p (X k ) -t k+1 grad X k H p (•) 2 + C(β, r, p)t 2 k+1 ≤ H p (X k ) -C p,µ,K t k+1 (H p (X k ) -H p (e)) + C(β, r, p)t 2 k+1
(by Proposition 2.2) so that the process (Y k ) k≥0 defined by (3.16)

Y 0 = H p (X 0 ) -H p (e) and for k ≥ 1 Y k = H p (X k ) -H p (e) -C(β,
E [H p (X k ) -H p (e)] .
We want to prove that a = 0. We already proved that

E [H p (X k+1 ) -H p (e)|F k ] ≤ H p (X k ) -H p (e) -C p,µ,K t k+1 (H p (X k ) -H p (e)) + C(β, r, p)t 2 k+1 .
(3.18)

Taking the expectation we obtain

(3.19) E [H p (X k+1 ) -H p (e)] ≤ (1 -t k+1 C p,µ,K )E [H p (X k ) -H p (e)] + C(β, r, p)t 2 k+1
so that proving that a = 0 is similar to the previous case. Finally (2.3) proves that ρ(X k , e) 2 converges in L 1 and almost surely to 0.

y t = d i=1 y i ε ε δλi-1 + t ε s δλi-1 δσ dB s , e i t 1-δλi e i , t ≥ ε.
Now by definition of C p,µ,K we clearly have

(3.22) C p,µ,K ≤ min 1≤i≤d λ i .
So the condition δC p,µ,K > 1 implies that for all i, δλ i -1 > 0, and as ε → 0, Assume that a continuous solution y t converging in probability to 0 as t → 0 + exists. Since y i ε ε δλi-1 → 0 in probability as ε → 0, we necessarily have using (3.23)

(3.24) y t = d i=1 t 1-δλi t 0 s δλi-1 δσ dB s , e i e i , t ≥ 0.
Note y i δ is Gaussian with variance

tδ 2 Γ(e * i ⊗ e * i ) 2δλ i -1
, so it converges in L 2 to 0 as t → 0.

Conversely, it is easy to check that equation (3.24) defines a solution to (3.20).

To prove the a.s. convergence to 0 we use the representation

t 0 s δλi-1 δσ dB s , e i = B i ϕi(t)
where B i s is a Brownian motion and ϕ i (t) =

δ 2 Γ(e * i ⊗ e * i ) 2δλ i -1 t 2δλi-1 .
Then by the law of iterated logarithm lim sup

t↓0 t 1-δλi B i ϕi(t) ≤ lim sup t↓0 t 1-δλi 2ϕ i (t) ln ln ϕ -1 i (t)
But for t small we have 2ϕ i (t) ln ln ϕ -1 i (t) ≤ t δλi-3/4 so lim sup

t↓0 t 1-δλi B i ϕi(t) ≤ lim t↓0 t 1/4 = 0.
This proves a.s. convergence to 0. Continuity is easily checked using the integral representation (3.24).

3.4. Proof of Theorem 2.6. Consider the time homogeneous Markov chain (Z n k ) k≥0 with state space [0, ∞) × T e M defined by

(3.25) Z n k = k n , Y n k .
The first component has a deterministic evolution and will be denoted by t n k ; it satisfies

(3.26) t n k+1 = t n k + 1 n , k ≥ 0.
Let k 0 be such that

(3.27) δ k 0 < δ 1 .
Using equations (2.8), (2.12) and (2.11), we have for k ≥ k 0 , (3.28)

Y n k+1 = nt n k + 1 √ n exp -1 e exp exp e 1 √ nt n k Y n k - δ nt n k + 1 grad 1 √ nt n k Y n k F p (•, P k+1 ) .
Consider the transition kernel P n (z, dz ′ ) on (0, ∞) × T e M defined for z = (t, y) by

P n (z, A) = P t + 1 n , nt + 1 √ n exp -1 e exp exp e 1 √ nt y - δ nt + 1 grad exp e 1 √ nt y F p (•, P 1 ) ∈ A (3.29)
where A ∈ B((0, ∞) × T e M ). Clearly this transition kernel drives the evolution of the Markov chain (Z n k ) k≥k0 . For the sake of clarity, we divide the proof of Theorem 2.6 into four lemmas. Lemma 3.1. Assume that either p ≥ 2 or e does not belong to the support supp(µ) of µ (note this implies that for all x ∈ supp(µ) the function F p (•, x) is of class C 2 in a neighbourhood of e). Fix δ > 0. Let B be a bounded set in T e M and let 0 < ε < T . We have for all

C 2 function f on T e M n f nt + 1 √ n exp -1 e exp exp e 1 √ nt y - δ nt + 1 grad exp e 1 √ nt y F p (•, x) -f (y) = d y f, y t - √ n d y f, δ grad e F p (•, x) -δ∇dF p (•, x) grad y f, y t + δ 2 2 Hess y f (grad e F p (•, x) ⊗ grad e F p (•, x)) + O 1 √ n (3.30) uniformly in y ∈ B, x ∈ supp(µ), t ∈ [ε, T ].
Proof. Let x ∈ supp(µ), y ∈ T e M , u, v ∈ R sufficiently close to 0, and q = exp e uy t . For s ∈ [0, 1] denote by a → c(a, s, u, v) the geodesic with endpoints c(0, s, u, v) = e and c(1, s, u, v) = exp exp e ( uy t ) -vs grad exp e ( uy t ) F p (•, x) : c(a, s, u, v) = exp e a exp -1 e exp exp e ( uy t ) -sv grad exp e ( uy t ) F p (•, x) .

This is a C 2 function of (a, s, u, v) ∈ [0, 1] 2 × (-η, η) 2 , η sufficiently small. It also depends in a C 2 way of x and y. Letting c(a, s) = c a, s,

1 √ n , δ nt + 1 , we have exp -1 e exp exp e 1 √ nt y - δ nt + 1 grad exp e 1 √ nt y F p (•, x) = ∂ a c(0, 1).
So we need a Taylor expansion up to order n -1 of nt + 1 √ n ∂ a c(0, 1).

We have c(a, s, 0, 1) = exp e (-as grad e F p (•, x)) and this implies

∂ 2 s ∂ a c(0, s, 0, 1) = 0, so ∂ 2 s ∂ a c(0, s, u, 1) = O(u). On the other hand the identities c(a, s, u, v) = c(a, sv, u, 1) yields ∂ 2 s ∂ a c(a, s, u, v) = v 2 ∂ 2 s ∂ a c(a, s, u, 1), so we obtain ∂ 2 s ∂ a c(0, s, u, v) = O(uv 2 ) and this yields ∂ 2 s ∂ a c(0, s) = O(n -5/2 ), uniformly in s, x, y, t. But since ∂ a c(0, 1) -∂ a c(0, 0) -∂ s ∂ a c(0, 0) ≤ 1 2 sup s∈[0,1] ∂ 2 s ∂ a c(0, s)
we only need to estimate ∂ a c(0, 0) and ∂ s ∂ a c(0, 0). Denoting by J(a) the Jacobi field ∂ s c(a, 0) we have

nt + 1 √ n ∂ a c(0, 1) = nt + 1 √ n ∂ a c(0, 0) + nt + 1 √ n J(0) + O 1 n 2 .
On the other hand

nt + 1 √ n ∂ a c(0, 0) = nt + 1 √ n y √ nt = y + y nt
so it remains to estimate J(0). The Jacobi field a → J(a, u, v) with endpoints J(0, u, v) = 0 e and

J(1, u, v) = -v grad exp e ( uy t ) F p (•, x) satisfies ∇ 2 a J(a, u, v) = -R(J(a, u, v), ∂ a c(a, 0, u, v))∂ a c(a, 0, u, v) = O(u 2 v). This implies that ∇ 2 a J(a) = O(n -2
). Consequently, denoting by P x1,x2 : T x1 M → T x2 M the parallel transport along the minimal geodesic from x 1 to x 2 (whenever it is unique) we have (3.31) P c(1,0),e J(1) = J(0) + J(0) + O(n -2 ) = J(0) + O(n -2 ).

But we also have

P c(1,0,u,v),e J(1, u, v) = P c(1,0,u,v),e -v grad c(1,0,u,v) F p (•, x) = -v grad e F p (•, x) -v∇ ∂ac(0,0,u,v) grad • F p (•, x) + O(vu 2 ) = -v grad e F p (•, x) -v∇dF p (•, x) uy t , • ♯ + O(vu 2 )
where we used ∂ a c(0, 0, u, v) = uy t and for vector fields A, B on T M and a

C 2 function f 1 on M ∇ Ae grad f 1 , B e = A e grad f 1 , B e -grad f 1 , ∇ Ae B = A e df 1 , B e -df 1 , ∇ Ae B = ∇df 1 (A e , B e ) which implies ∇ Ae grad f 1 = ∇df 1 (A e , •) ♯ .
We obtain

P c(1,0),e J(1) = - δ nt + 1 grad e F p (•, x) - δ √ n(nt + 1) ∇dF p (•, x) y t , • ♯ + O(n -2 ).
Combining with (3.31) this gives

J(0) = - δ nt + 1 grad e F p (•, x) - δ nt + 1 ∇dF p (•, x) y √ nt , • ♯ + O 1 n 2 . So finally (3.32) nt + 1 √ n ∂ a c(0, 1) = y + y nt - δ √ n grad e F p (•, x) -δ∇dF p (•, x) y nt , • ♯ + O n -3/2 .
To get the final result we are left to make a Taylor expansion of f up to order 2.

Define the following quantities:

(3.33) b n (z) = n {|z ′ -z|≤1} (z ′ -z)P n (z, dz ′ ) and (3.34) a n (z) = n {|z ′ -z|≤1} (z ′ -z) ⊗ (z ′ -z)P n (z, dz ′ ).
The following property holds: Proof.

(1) We use the notation z = (t, y) and z ′ = (t ′ , y ′ ). We have

1 {|z ′ -z|>1} P n (z, dz ′ ) = 1 {max(|t ′ -t|,|y ′ -y|)>1} P n (z, dz ′ ) = 1 {max( 1 n ,|y ′ -y|)>1} P n (z, dz ′ ) = P nt + 1 √ n exp -1 e exp exp e 1 √ nt y - δ nt + 1 grad exp e 1 √ nt y F p (•, P 1 ) -y > 1 .
On the other hand, since

F p (•, x) is of class C 2 in
b n (z) = n (z ′ -z) P n (z, dz ′ ) = n 1 n , E nt + 1 √ n exp -1 e exp exp e y √ nt - δ nt + 1 grad exp e y √ nt F p (•, P 1 ) -y .
We have by lemma 3.1 Proof. First consider the case p ∈ [1, 2). We know by (3.12) that there exists some constant C(β, r, p) such that

n nt + 1 √ n exp -1 e exp exp e 1 √ nt y - δ nt + 1 grad exp e 1 √ nt y F p (•, P 1 ) -y = 1 t y -δ √ n grad e F p (•, P 1 ) -δ∇dF p (•, P 1 ) 1 t y, • ♯ + O 1 n 1/2
(3.43) E ρ 2 (e, X k+1 ) ≤ E ρ 2 (e, X k ) exp (-C p,µ,K t k+1 ) + C(β, r, p)t 2 k+1 .
From this (3.42) is a consequence of Lemma 0.0.1 (case α > 1) in [START_REF] Nedic | Convergence rate of incremental subgradient algorithms, Stochastic optimization : algorithms and applications[END_REF]. We give the proof for completeness. We deduce easily by induction that for all k ≥ k 0 ,

E ρ 2 (e, X k ) ≤ E ρ 2 (e, X k0 ) exp   -C p,µ,K k j=k0+1 t j   + C(β, r, p) k i=k0+1 t 2 i exp   -C p,µ,K k j=i+1 t j   , (3.44)
where the convention k j=k+1 t j = 0 is used. With t n = δ n , the following inequality holds for all i ≥ k 0 and k ≥ i:

(3.45) k j=i+1 t j = δ k j=i+1 1 j ≥ δ k+1 i+1 dt t ≥ δ ln k + 1 i + 1 .
Hence, This implies that the sequence kE ρ , the sequence of processes. We prove that from any subsequence Ỹ φ(n) ε n≥1

E ρ 2 (e, X k ) ≤ E ρ 2 (e, X k0 ) k 0 + 1 k + 1 δCp,µ,K + δ 2 C(β, r, p) (k + 1) δCp,µ,K k i=k0+1 (i + 1) δCp,µ,K i 2 . ( 3
, we can extract a further subsequence

Ỹ ψ(n) ε n≥1
that weakly converges in D([ε, 1], R d ).

Let us first prove that Ỹ φ(n) weakly converges to a diffusion (y t ) ε≤t≤T with generator G δ (t) given by (2.13) and such that y ε has law ν ε . This achieves the proof of lemma 3. This subsequence converges in D((0, T ], R d ) to (y ′ (t)) t∈(0,T ] . On the other hand, as in the proof of lemma 3.4, we have Ỹ η(n) (t) 2 2 ≤ Ct for some C > 0. So Ỹ η(n) (t) 2 2 → 0 as t → 0, which in turn implies y ′ (t) 2 2 → 0 as t → 0. The unicity statement in Proposition 2.5 implies that (y ′ (t)) t∈(0,T ] and (y δ (t)) t∈(0,T ] are equal in law. This achieves the proof.

ε (ε) n≥1 is bounded in L 2 . Ỹ φ(n) ε (ε) 2 2 = [φ(n)ε] 2 φ(n) E ρ 2 (e, X [φ(n)ε] ) ≤ ε sup

2 j

 2 k+1 so that the process (Y k ) k≥0 defined by (3.9) Y 0 = ρ(X 0 , e) 2 and for k ≥ 1 Y k = ρ(X k , e) 2 -C(β, r, p) is a bounded supermartingale. So it converges in L 1 and almost surely. Consequently ρ(X k , e) 2 also converges in L 1 and almost surely. Let (3.10) a = lim k→∞ E ρ(X k , e) 2 .

s 0 s

 0 δλi-1 δσ dB s , e i → t δλi-1 δσ dB s , e i in probability.

Lemma 3 . 2 . 1 )( 2 )( 3 )

 32123 Assume that either p ≥ 2 or e does not belong to the support supp(µ). (For all R > 0 and ε > 0, there exists n 0 such that for all n ≥ n 0 and z ∈ [ε, T ] × B(0 e , R), where B(0 e , R) is the open ball in T e M centered at the origin with radius R, (3.35) 1 {|z ′ -z|>1} P n (z, dz ′ ) = 0. For all R > 0 and ε > 0and L δ (y) = y -δ∇dH(y, •) ♯ . For all R > 0 and ε > 0, (3.38) lim n→∞ sup z∈[ε,T ]×B(0e,R) |a n (z) -a(z)| = 0 with (3.39) a(z) = δ 2 diag(0, Γ) and Γ = E [grad e F p (•, P 1 ) ⊗ grad e F p (•, P 1 )] .

2 ∼ δ 2 C 1 andE ρ 2

 2212 .46) For δC p,µ,K > 1 we have as k → ∞ (3.47)δ 2 C(β, r, p) (k + 1) δCp,µ,K k i=k0+1 (i + 1) δCp,µ,K i (β, r, p) (k + 1) δCp,µ,K k δCp,µ,K -1 δC p,µ,K -1 ∼ δ 2 C(β, r, p) δC p,µ,K -1 k -(e, X k0 ) k 0 + 1 k + 1 δCp,µ,K = o(k -1 ).

n≥1 nE ρ 2

 2 (e, X n ) and the last term is bounded by lemma 3.3.Consequently Ỹ φ(n) ε (ε) n≥1 is tight. So there is a subsequence Ỹ ψ(n) ε (ε) n≥1that weakly converges in T e M to the distribution ν ε . Thanks to Skorohod theorem which allows to realize it as an a.s. convergence and to lemma 3.2 we can apply Theorem 11.2.3 of[START_REF] Stroock | Multidimensional diffusion processes[END_REF], and we obtain that the sequence of processes Ỹ ψ(n) ε n≥1

4 .

 4 Proof of Theorem 2.6.Let Ỹ n = Y n [nt] 0≤t≤T. It is sufficient to prove that any subsequence of Ỹ n n≥1 has a further subsequence which converges in law to (y δ (t)) 0≤t≤T . So let Ỹ φ(n) n≥1 a subsequence. By lemma 3.4 with ε = 1/m there exists a subsequence which converges in law on [1/m, T ]. Then we extract a sequence indexed by m of subsequence and take the diagonal subsequence Ỹ η(n) .

  (y t ) dt + δσ dB t where L δ (y) = y -δ∇dH p (y, •) ♯ and B t and σ are as in Proposition 2.5. This sde can be solved explicitely on [ε, ∞). The symmetric endomorphism y → ∇dH p (y, •) ♯ is diagonalisable in the orthonormal basis (e i ) 1≤i≤d with eigenvalues (λ i ) 1≤i≤d . The endomorphism L δ = id -δ∇dH p (e)(id, •) ♯ is also diagonalisable in this basis with eigenvalues (1 -δλ i ) 1≤i≤d . The solution y t =

	3.3. Proof of Proposition 2.5. Fix ε > 0. Any diffusion process on [ε, ∞) with generator G δ (t) is solution of a sde of the type
	(3.20)	dy t =	1 t	L δ d
				y i t e i of (3.20) started at
				i=1
	d
	y ε =	y i ε e i is given by
	i=1
	(3.21)

  = δ 2 E [grad e F p (•, P 1 ) ⊗ grad e F p (•, P 1 )] + o(1)uniformly in z ∈ [ε, T ] × B(0 e , R), so this yields (3.38).

		converges to
	(3.41)		1 t	y -E δ∇dF p (•, P 1 )	1 t	y, •	♯	=	1 t	y -δ∇dH p	1 t	y, •	♯	.
		Moreover the convergence is uniform in z ∈ [ε, T ] × B(0 e , R), so this yields (3.36).
	(3) In the same way, using lemma 3.1,
		n (y ′ -y) ⊗ (y ′ -y) P n (z, dz ′ )
		=	1 n	E -	√ nδ grad e F p (•, P 1 ) ⊗ -	√ nδ grad e F p (•, P 1 ) + o(1)
	Lemma 3.3. Suppose that t n =	δ n	for some δ > 0. For all δ > C -1 p,µ,K ,
	(3.42)						sup
		a.s. uniformly in n, and since
							E δ	√ n grad e F p (•, P 1 ) = 0,
		this implies that
	n E	nt + 1 √ n	exp -1 e	exp exp e 1 √ nt y -	δ nt + 1	grad exp e 1 √ nt y F p (•, P 1 )	-y

n≥1 nE ρ 2 (e, X n ) < ∞.

  2 (e, X k ) is bounded. Next consider the case p ≥ 2. Now we have by(3.19) that (3.48) E [H p (X k+1 ) -H p (e)] ≤ E [H p (X k ) -H p (e)] exp (-C p,µ,K t k+1 ) + C(β, r, p)t 2k+1 . From this, arguing similarly, we obtain that the sequence kE [H p (X k ) -H p (e)] is bounded. We conclude with (2.3). Lemma 3.4. Assume δ > C -1 p,µ,K and that H p is C 2 in a neighbourhood of e. For all 0 < ε < T , the sequence of processes Y n [nt] ε≤t≤T is tight in D([ε, T ], R d ).

	Proof. Denote by Ỹ n ε = Y n [nt] ε≤t≤T n≥1