{'fraction_non_alphanumeric': 0.05325252678774983, 'fraction_numerical': 0.02051827646477696, 'mean_word_length': 3.9023935922510944, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A stochastic partial differential equation model for limit order book dynamics

We propose an analytically tractable class of models for the dynamics of a limit order book, described through a stochastic partial differential equation (SPDE) with multiplicative noise for the order book centered at the mid-price, along with stochastic dynamics for the mid-price which is consistent with the order flow dynamics. We provide conditions under which the model admits a finite dimensional realization driven by a (low-dimensional) Markov process, leading to efficient estimation and computation methods. We study two examples of parsimonious models in this class: a two-factor model and a model with mean-reverting order book depth. For each model we analyze in detail the role of different parameters, the dynamics of the price, order book depth, volume and order imbalance, provide an intuitive financial interpretation of the variables involved and show how the model reproduces statistical properties of price changes, market depth and order flow in limit order markets.

Financial instruments such as stocks and futures are increasingly traded in electronic, order-driven markets, in which orders to buy and sell are centralized in a limit order book and market orders are executed against the best available offers in the limit order book. The dynamics of prices in such markets are not only interesting from the viewpoint of market participants -for trading and order executionbut also from a fundamental perspective, since they provide a detailed view of the dynamics of supply fand demand and their role in price formation.
The availability of a large amount of high frequency data on order flow, transactions and price dynamics on these markets has instigated a line of research which, in contrast to traditional market microstructure models which make assumptions on the behavior and preferences of various types of agents, focuses on the statistical modeling of aggregate order flow and its relation with price dynamics, in a quest to understand the interplay between price dynamics and order flow of various market participants (Cont, 2011).
A fruitful line of approach to these questions has been to model the stochastic dynamics of the limit order book, which centralizes all buy and sell orders, either as a queueing system (Luckock, 2003;Smith et al., 2003;Cont et al., 2010;Cont and De Larrard, 2012;Cont and de Larrard, 2013;Kelly and Yudovina, 2018) or, at a coarse-grained level, through a (stochastic) partial differential equation describing the evolution of the distribution of buy and sell orders (Lasry and Lions, 2007;Caffarelli et al., 2011;Burger et al., 2013;Carmona and Webster, 2013;Markowich et al., 2016;Hambly et al., 2020;Horst and Kreher, 2018). These PDE models may be viewed as scaling limits of discrete point process models (Cont and De Larrard, 2012;Hambly et al., 2020;Horst and Kreher, 2018).
Although joint modeling of order flow at all price levels in the limit order book is more appealling, (S)PDE models have lacked the analytical and computational tractability needed for applications; as a result, most analytical results have been derived using reduced-form models of the best bid-ask queues (Cont and De Larrard, 2012;Cont and de Larrard, 2013;Chavez-Casillas and Figueroa-Lopez, 2017;Huang et al., 2017).
We propose a class of stochastic models for the dynamics of the limit order book which represent the dynamics of the entire order book while retaining at the same time the analytical and computational tractability of low-dimensional Markovian models, and provides realistic dynamics for the joint dynamics of the market price and order book depth. Starting with a description of the dynamics of the limit order book via a stochastic partial differential equation (SPDE) with multiplicative noise, we show that in many cases, the solutions of this equation may be parameterized in terms of a low-dimensional diffusion process, which then makes the model computationally tractable. In particular, we are able to derive analytical relations between model parameters and various observable quantities. This feature may be used for calibrating model parameters to match statistical features of the order flow and leads to empirically testable predictions, which we proceed to test using high frequency time series of order flow in electronic equity markets.
Outline Section 1 introduces a description of the dynamics of a limit order book through a stochastic partial differential equation (SPDE). We describe the various terms in the equation, their interpretation and discuss the implications for price dynamics (Section 1.3). This class of models is part of a more general family of SPDEs driven by semimartingales, introduced in Sec. 1.5 and studied in Sec. 2.
We then focus on two analytically tractable examples: a two-factor model (Section 3) and a model with mean-reverting depth and imbalance (Section 4). For each model we perform a detailed analysis of the role of different parameters and study the dynamics of the price, order book depth, volume and order imbalance,  Figure 1. Snapshot of the NASDAQ limit order book for CISCO shares (Jan 30, 2018), displaying outstanding buy orders (green) and sell orders (red) awaiting execution at different prices. The highest buying price ($ 42.15 in this example) is the bid and the lowest selling price ($ 42.16) is the ask.
provide an intuitive financial interpretation of the variables involved and show how the model may be estimated from financial time series of price, volume and order flow.

A stochastic PDE model for limit order book dynamics
We consider a market for a financial asset (stock, futures contract, etc.) in which buyers and sellers may submit limit orders to buy or sell a certain quantity of the asset at a certain price, and market orders for immediate execution against the best available price. 1 Limit orders awaiting execution are collected in the limit order book, an example of which is shown in Figure 1: at any time t, the state of the limit order book is summarized by the volume V (t, p) of orders awaiting execution at price levels p on a grid with mesh size given by the minimum price increment or tick size δ. By convention we associate negative volumes with buy orders and positive volumes with sell orders, as shown in Figure 1. An admissible order book configuration is then represented by a function p → V (p) such that s b (V ) (resp. s a (V )) is called the bid (resp. ask) price and represents the price associated with the best buy (resp. sell) offer. The quantity is called the mid-price and the difference s a (V ) − s b (V ) is called the bid-ask spread.
In the example shown in Figure 1, s b (V ) = 42.15, s a (V ) = 42.16 and the bid-ask spread is equal in this case to the tick size, which is 1 cent. One modelling approach has been to represent the dynamics of V (t, p) as a spatial (marked) point process (Luckock, 2003;Cont et al., 2010;Cont and de Larrard, 2013;Kelly and Yudovina, 2018). These models preserve the discrete nature of the dynamics at high frequencies but can become computationally challenging as one tries to incorporate realistic dynamics. In particular, price dynamics, which is endogenous in such models, is difficult to study, even when the order flow is a Poisson point process.
When the bid-ask spread and tick size δ are much smaller than the price level, as is often the case, another modelling approach is to use a continuum approximation for the order book, describing it through its density v(t, p) representing the volume of orders per unit price: V (t, p) v(t, p)δ. The evolution of the density of buy and sell orders is then described through a partial differential equation (PDE). A deterministic description of the dynamics of order densities through a system of coupled partial differential equations was proposed by (Lasry and Lions, 2007) and studied in detail by (Chayes et al., 2009;Caffarelli et al., 2011;Burger et al., 2013). In the Lasry-Lions model, the evolution of the density of buy and sell orders is described by a pair of diffusion equations coupled through the dynamics of the price, which represents the free boundary between prices of buy and sell orders. This model is appealing in many respects, especially in terms of analytical tractability, but leads to a deterministic price process which decays to a constant price, so does not provide any insight into the relation between liquidity, depth, order flow and price volatility. (Markowich et al., 2016) explore some stochastic extensions of this model but essentially show that these extensions do not provide realistic price dynamics.
We adopt here this continuum approach for the description of the limit order book, but describe instead its dynamics through a stochastic partial differential equation, paying close attention to price dynamics and its relation with order flow.
The model we propose shares some features with (Lasry and Lions, 2007), but also has some essential differences. Unlike the Lasry-Lions model, which is a free boundary problem in which the dynamics of the price is implicitly determined, we formulate the model as a stochastic partial differential equation in relative price coordinates, which leads to a stochastic moving boundary problem in absolute price coordinates. This leads to a more realistic joint dynamics for the market price and order book depth which can be related to empirical observations. Our model also relates to the classes of models studied in (Horst and Kreher, 2018;Hambly et al., 2020) as scaling limits of discrete queueuing systems.
We now describe our model in some detail.
1.1. State variables and scaling transformations. We focus on the case where the tick size δ and the bid-ask spread are small compared to the typical price level and consider a limit order book described in terms of a mid-price S t and the density v(t, p) of orders at each price level p, representing buy orders for p < S t , and sell orders for p > S t . We use the convention, shown in Figure 1, of representing buy orders with a negative sign and sell orders with a positive sign, so v(t, p) ≤ 0 for p < S t and v(t, p) ≥ 0 for p > S t Limit orders are executed against market orders according to price priority and their position in the queue; execution of a limit order only occurs if they are located at the best (buy/sell) prices. This means that price dynamics is determined by the interaction of market orders with limit orders of opposite type at or near the interface defined by the best price (Cont et al., 2010). Due to this fact, most limit orders flow are submitted close to the best price levels: the frequency of limit order submissions is highly inhomogeneous as a function of distance to the best price and concentrated near the best price. As shown in previous empirical studies, order flow intensity at a given distance from the best price can be considered as a stationary variable in a first approximation (Bouchaud et al., 2009;Cont et al., 2010). For this reason, in a stochastic description it is more convenient to model the dynamics of order flow in the reference frame of the (mid-)price S t . We define where x represents a distance from the mid-price. We refer to u t as the centered order book density.
The simplest way of centering is to set x(p) = p−S t but other, nonlinear, scalings may be of interest. Although limit orders may be placed at any distance from the bid/ask prices, price dynamics is dominated by the behavior of the order book a few levels above and below the mid price (Cont and De Larrard, 2012). This region becomes infinitesimal if the tick size δ is naively scaled to zero, suggesting that the correct scaling limit is instead one in which we choose as coordinate a scaled version (p − S t ), as classically done in boundary layer analysis of PDEs (Schlichting and Gersten, 2017), in order to zoom into the relevant region: for bid and ask side, respectively. We will consider examples of such nonlinear scalings when discussing applications to high-frequency data in Sections 3 and 4. These arguments also justify limiting the range of the argument xto a bounded interval [−L, +L], setting u t (x) = 0 for x / ∈ (−L, L). This amounts to assuming that no orders are submitted at price levels at distances |x| ≥ L from the midprice and that orders previously submitted at some price p are cancelled as soon as |S t − p| ≥ L i.e. when the mid-price S t moves away from p by more than L. When L is a large multiple of daily volatility, this is a realistic assumption. In some market (for example futures contracts), limit orders can be in fact only submitted within a range ±L of the mid-price.
1.2. Dynamics of the centered limit order book. Empirical studies on intraday order flow in electronic markets reveal the coexistence of two, very different types of order flow operating at different frequencies (Lehalle and Laruelle, 2018).
On one hand, we observe the submission (and cancellation) of orders queueing at various price levels on both sides of the market price by regular market participants. Cancellation may occur in several ways: we distinguish outright cancellations, which we model as proportional to current queue size, from cancellations with replacement ('order modifications'), in which an order is cancelled and immediately replaced by another one of the same type, usually at a neighboring price limit. The former results in a net decrease in the volume of the order book whereas the latter is conservative and simply shifts orders across neighboring levels of the book. Further decomposing this conservative flow into a symmetric and antisymmetric part leads to two terms in the dynamics of u t : a diffusion term representing the cancellation of orders and their (symmetric) replacement by orders at neighboring price levels and a convection (or transport) term representing the cancellation of orders and their replacement by orders closer to the mid-price. Denoting by ∇ the gradient in the variable x, the net effect of this order flow on the order book may thus be described as a superposition of a term f b (x) (resp. f a (x) ) representing the rate of buy (resp. sell) order submissions at a distance x from the best price; a term α b u t (x) (resp. α a u t (x)) representing (outright) proportional cancellation of limit buy (resp. sell) orders at a distance x from the mid-price (where α a , α b ≤ 0). a convection term −β b ∇u t (x) (resp. +β a ∇u t (x)) with β a , β b > 0 which models the replacement of buy (resp. sell) orders by orders closer to the mid-price (i.e. closer to x = 0, hence the signs in these terms): in the reference frame where the origin is the mid-price, this translates into a flow of volume towards the origin; a diffusion term η b ∆u t (x) (resp. η a ∆u t (x)) which represents the cancellation and symmetric replacement of orders at a distance x from the midprice.
Another component of order flow is the one generated by high-frequency traders (HFT). These market participants buy and sell at very high frequency and under tight inventory constraints, submitting and cancelling large volumes of limit orders near the mid-price and resulting in an order flow whose net contribution to total order book volume is zero on average over longer time intervals but whose sign over small time intervals fluctuates at high frequency. At the coarse-grained time scale of the average (non-HF) market participants, these features may be modeled as a multiplicative noise term of the form σ b u t (x)dW b for buy orders (x < 0) and σ a u t (x)dW a for sell orders (x > 0) where (W a , W b ) is a two-dimensional Wiener process (with possibly correlated components). The multiplicative nature of the noise accounts for the high-frequency cancellations associated with HFT orders.
The impact of these different order flow components may be summarized by the following stochastic partial differential equation for the centered order book density u: Here η a , η b , β b , β a , σ a , σ b ∈ (0, ∞), α a , α b ≤ 0 and f a , f b : I → [0, ∞) although the equation may be equally considered without these sign restrictions.
Note that, unlike the Lasry-Lions model, there is no 'smooth pasting' condition at x = 0: in general ∇u t (0+) = ∇u t (0−): the difference ∇u t (0+) − ∇u t (0−) is in fact random and represents an imbalance in the flow of buy and sell orders, which drives price dynamics. This important feature is discussed in Section 1.3 below.
Remark 1.1. In simple price impact models used in the literature on optimal trade execution it is assumed that the relation between price impact and order size is deterministic. This corresponds to the case αu + f = β = σ = η = 0 which leads to a constant centered order book profile u t (.) = u 0 (.). These terms thus correspond to deformations of the centered order book profile due to new order book events and lead to a stochastic market impact of trades dependent on the current state of the order book.
The existence of a solution satisfying the boundary and sign constraints is not obvious but we will see in Section 2 that (1.2) is well-posed: it follows from (Da Prato and Zabczyk, 2014, Theorem 6.7) and (Milian, 2002, Theorem 3) that, when f a , f b ∈ L 2 (I), then for all u 0 ∈ L 2 (I) there exists a unique weak solution of (1.2) (see 6 Definition 2.2 below) and, when u 0 | (0,L) ≥ 0 and u 0 | (−L,0) ≤ 0 this solution satisfies We will study the mathematical properties of the solution in more detail below.
1.3. Price dynamics. The dynamics of the limit order book determines the dynamics of the bid and ask price, which corresponds to the location of the best (buy and sell) orders. The dynamics of the price should thus be related to the arrival and execution of orders in the order book.
To understand the relation between price dynamics and order flow, let us take a step back and consider an order book with discrete price levels, multiples of a tick size δ, D b orders per level on the bid side and D a orders per level on the ask side. Price changes during a time interval [t, t + ∆t] are triggered through the interaction of the net order flow, or order flow imbalance (OFI) and the outstanding limit orders at the top of the order book (Cont et al., 2014). As illustrated in Figure 2, an order flow imbalance of ∆D a t > 0 on the ask side over a short time interval [t, t + ∆t] represents an excess of buy orders, which will then be executed against limit sell orders sitting on the ask side and move the ask price by ∆D a t /D a ticks, resulting in a price move of δ ∆D a t /D a . Similarly, an order flow imbalance ∆D b t on the bid side will move the bid price up by ∆D b t /D b ticks. Using our sign conventions for buy/sell volumes, this leads to the following dynamics: so the dynamics of the mid price s t = (s b t + s a t )/2 is given by This relation is exact (up to rounding) in the case of a discrete order book with constant depth per level (and thus, no empty levels), as shown in Figure 2. However, in a dynamic setting where the order book may have an arbitrary profile which randomly shifts at each instant, one can only expect a 'homogenized' version of (1.4) to hold: where θ is an impact coefficient which relates order imbalance to price movements. This relation between order flow imbalance and price movements has been empirically verified in equity markets (Cont et al., 2014), and we shall use it as a basis for defining the relation between price dynamics and order flow in our model. Let us now see how the relation (1.5) translates in terms of the variables in our model. Denoting by D b t (resp. D a t ) the volume of buy (resp. sell) limit orders at the top of the book (i.e. the first or average of the first few levels). Given a mid-price S ∈ R + , we define a scaling transformation x : [S, S + L] → [0, ∞) as discussed in Section 1.1, with continuously differentiable inverse and such that x(S) = 0. The volume D a in the best ask queue is then given by D b may be similarly defined for the bid side. These quantities represent the depth at the top of the book; we will refer to them as 'market depth'. In the case of linear Similarly, for the bid side Substituting these expressions in (1.5), we obtain the following dynamics of the mid-price: We observe that price dynamics is entirely determined by the order flow at the top of the book and the depth of the limit order book around the mid-price. The tick size δ, used in the derivation, does not appear anymore in (1.9). The only trace of the microstructure is the impact coefficient θ which relates the order flow imbalance to the magnitude of the price change, and whose amplitude may vary across assets.
Remark 1.2. Equation (1.9) requires left and right-differentiability of u at the origin. This can be guaranteed whenever u t takes values in the Sobolev space H 2γ (I), for some γ > 3 /4 which will be the case in our model. Note however that, in contrast to (Lasry and Lions, 2007), in general ∇u(0+) = ∇u(0−): the difference between these two quantities is proportional to the order flow imbalance which drives price moves.
Remark 1.3. As noted in Remark 1.1, in the case αu + f = β = σ = η = 0 corresponds to a constant centered order book profile u t = u 0 . In this case, Equation (1.9) implies dS t = 0 i.e. the price is constant, which is consistent with a zero net order flow. This is a (desirable) consequence of the consistency between the price dynamics (1.9) and the order book dynamics (1.2).
1.4. Dynamics in absolute price coordinates. The model above describes dynamics of the order book in relative price coordinates, i.e. as a function of the (scaled) distance from the mid-price. The density of the limit order book parameterized by the (absolute) price level p ∈ R is given by where we extend u t to R by setting u t (y) = 0 for y ∈ R \ [−L, L]. Assume S t follows an (arbitrary) Itô process where θ > 0 and µ t is predictable and integrable and ξ a t and ξ b t are predictable and square-integrable processes. This includes the case of price dynamics (1.9), which can be used to express µ t , ξ a t , ξ b t in terms of u t and model parameters. We will not go into such detail here but will return to this in the examples in Section 3 and 4. Defineξ Using a (generalized) Itô-Wentzell formula (see Appendix A), we can show that v is the solution of a stochastic moving boundary problem (Mueller, 2018): We refer to (1.13) as a stochastic boundary condition at S t .
Here, we assumed for simplicity that f a , f b = 0. A more detailed discussion of this result is given in Appendix A.
1.5. Linear evolution models for order book dynamics. We will now describe a more general class of linear models for order book dynamics, rich enough to cover the examples we discussed so far, but also covering all level-1 models where the best bid and ask queue are modeled by positive semimartingales. Generally, the densities of orders in the bid and ask side will take values in some function spaces H b and H a , respectively. We assume that orders at relative price level x for |x| ≥ L ∈ (0, ∞] will be cancelled. The relative price levels are on the bid side I b := (−L, 0), and on the ask side I a := (0, L). Then, in order to preserve the interpretation of a density it will be reasonable to ask . From mathematical side, we will assume that H a and H b are real separable Hilbert spaces. For notational convenience we now also set I : The density of limit orders at relative price level x and time t is given by u : I × [0, ∞) × Ω → R, such that u := u| I is an H -valued adapted process. The initial state is described by h : I → R, such that h := h| I is an element in H . The (averaged) intra-book dynamics are modeled by linear operators A : dom(A ) ⊂ H → H , for ∈ {a, b}, which we assume to be densely defined and such that for ∈ {a, b} there exist weak solutions in H of the equations (1.14) The random order arrivals and cancellations are assumed to be proportional and are modelled by cadlag semimartingales X b and X a , which we assume to have jumps greater than −1 almost surely. We assume the initial order book state is denoted by h ∈ H and we write h a := h| I a , h b := h| I b .
Model 1.4 (Linear Homogeneous Evolution). The general form of the linear homogeneous model is for t ≥ 0, and u 0 = h. u can be alternatively expressed as where g b and g a are solutions of (1.14), see Theorem 2.5 below. If, in addition, t → ∇g b t (0−) and t → ∇g a t (0+) are of bounded variation, then we obtain the price dynamics (1.9).
Corollary 1.5. Assume the setting of Model 1.4 and, in addition, that h is an eigenfunction of −A with eigenvalue ν ∈ R, for = b and = a. Then, (1.14) can be solved explicitly and Remark 1.6. In case that X b and X a are (local) martingales, the eigenvalues −ν b and −ν a play the role of net order arrival rates on bid and ask side, respectively.
Model 1.7 (Linear models with source terms). A more realistic setting assumes in addition an influx/outflow of orders at a rate f a (x), f b (x) which depends on the distance x to the mid price (Cont et al., 2010). The equation then becomes: for t ≥ 0, with initial condition u 0 = h. As we will discuss in Section 4, an interesting case is when f b (resp. f a ) is an eigenfunction of −A b (resp. −A a ) associated with some eigenvalue ν b (resp. ν a ). Then by Theorem 2.10 we obtain where, for ∈ {a, b}, Z t is the solution of (1.20) Remark 1.8. If ν b , ν a > 0 the state of the order book is mean reverting to the state . We will give an example of such a mean-reverting order book model in Section 4.
Remark 1.9. Any model for the dynamics of the order book implies a model for price dynamics via (1.9). In particular this implies a relation between price volatility and parameters describing order flow, in the spirit of (Cont and de Larrard, 2013). We will derive this relation for the examples studied in the sequel and use it to construct a model-based intraday volatility estimator.
In the next section, we will study this class of models from a mathematical point of view. We will then continue with the analysis of the two examples mentioned above in Sections 3 and 4.

Linear stochastic PDE models with multiplicative noise
In order to further study the properties of the SPDE model (1.2), we require a more explicit characterization of the solution, in order to compute various quantities of interest and estimate model coefficients from observations. A useful approach is to look for a finite dimensional realization of the infinite-dimensional process u: Availability of a finite dimensional realization for the SPDE (1.2) makes simulation, computation and estimation problems more tractable, especially if the process Z is a low-dimensional Markov process. Existence of such finite-dimensional realizations for stochastic PDEs have been investigated for SPDEs arising in filtering (Lévine, 1991) and interest rate modelling (Filipovic and Teichmann, 2003;Gaspar, 2006).
We will now show that finite dimensional realizations may indeed be constructed for a class of SPDEs which includes (1.2), and use this representation to perform an analytical study of these models.
2.1. Homogeneous equations. We now consider a more general class of linear homogeneous evolution equations with multiplicative noise taking values in a real separable Hilbert space (H, ·, · H ). Typically, H will be a function space such as L 2 (I) for some interval I ⊂ R. We consider the following class of evolution equations: where X is a real càdlàg semimartingale whose jumps satisfy ∆X t > −1 a. s. and A : dom(A) ⊂ H → H a linear operator on H whose adjoint we denote by A * . We assume that dom(A) ⊂ H is dense, and A is closed. Since A is closed we have that also dom(A * ) ⊂ H is dense and that A * * = A (Yosida, 1995, Theorem VII.2.3).
Definition 2.2. An adapted H-valued stochastic process (u t ) is an (analytical) weak solution of (2.1) with initial condition h 0 if, for all ϕ ∈ dom(A * ), [0, ∞) t → u t , ϕ H ∈ R is càdlàg a. s. and for each t ≥ 0, a. s.
The case X ≡ 0 corresponds to a notion of weak solution for the PDE: That is, for all ϕ ∈ dom(A * ), where the integral on the right hand side is assumed to exist. 2 In particular, this yields that [0, ∞) t → g t , ϕ H ∈ R is continuous. Denote by Z t = E t (X) the stochastic exponential of X. We recall the following useful lemma (see e.g. (Karatzas and Kardaras, 2007, Lemma 3.4)): Then, E t (X)E t (Y ) = 1 almost surely, for all t ≥ 0. Moreover, Theorem 2.5. Let Z := E(X), h 0 ∈ H. Then every weak solution of (2.1) is of the form u t := Z t g t , t ≥ 0 where g is a weak solution of (2.2).
Remark 2.6. In particular, the SPDE (2.1) admits a two dimensional realization in the sense of Definition 2.1 with factor process (t, E t (X)) and φ(t, y) := yg t .
Proof. Set u t := g t Z t , t ≥ 0, and for ϕ ∈ D(A * ) write B ϕ t := g t , ϕ H , C ϕ t := B ϕ t Z t = u t , ϕ H . Since t → g t , ϕ H is continuous and Z is scalar and càdlàg, we get that t → u t , ϕ H is càdlàg. Note that B ϕ is of finite variation and Z is a semimartingale, so that also C ϕ is a semimartingale. Moreover, by Itô product rule and since B ϕ is of finite variation and continuous, is (2.1). Now, let u be a solution of (2.1) and set and Z t := E t (Y ), t ≥ 0. Recall that by Lemma 2.4 we have Z t E t (X) = 1 for all t ≥ 0. Set g t := Z t u t , and, as above, fix ϕ ∈ dom(A * ) and write B ϕ t := u t Z t , ϕ H = g t , ϕ H and C ϕ t := u t , ϕ H . By Itô's product rule and Lemma 2.4, Thus, g is a weak solution of (2.2).
Example 2.7. Let A be the generator of a strongly continuous semigroup (S t ) t≥0 . Then, for h 0 ∈ H define g t := S t h 0 , t ≥ 0, which is a weak solution of (2.2). By Theorem 2.5.
Remark 2.8. If h 0 is an eigenfunction of A with eigenvalue ν, then, g t = e νt h 0 is the unique locally H-integrable solution of (2.2), and the unique solution of (2.1) is given by u t := h 0 e νt E t (X).
2.2. Inhomogeneous equations. We keep the assumptions on A, h 0 and X from the previous section and let f ∈ H. We now consider the inhomogeneous linear evolution equations Definition 2.9. A weak solution of (2.7) is an adapted H-valued stochastic process u such that for all ϕ ∈ dom(A * ) the mapping [0, ∞) t → u t , ϕ H is càdlàg and We exclude the cases α = 0 or f ≡ 0 which correspond to the homogeneous case discussed above. Let us first consider the case where A admits at least one eigenfunction.
Theorem 2.10. Suppose that f ∈ dom(A) is an eigenfunction for A with eigenvalue λ ∈ R, and let z 0 > 0 and Z be the solution of Then: (i) The stochastic process defined by u t = Z t f , t ≥ 0, is a solution of (2.7) with initial condition h 0 := z 0 f . (ii) Let, in addition, h 0 ∈ H be such that there exists a weak solution g = (g t ) t≥0 of the deterministic equation Then, u t := g t E t (X) + f Z t is a solution of (2.7) with initial condition h 0 . (iii) Let h 0 ∈ H be such that there exists a weak solution u = (u t ) t≥0 of (2.7) with initial condition h 0 . Then, g := (u − f Z)E(X) −1 , is a weak solution of (2.9).
Remark 2.11. Let (Z 1 t ) t≥0 and (Z 2 t ) t≥0 be given by (2.8) with respective initial data Proof. Part (i) follows by direct a computation: Let ϕ ∈ H, then for t ≥ 0, Similarly, we obtain that any solution u of (2.7) with initial data h 0 ∈ H can be written as where u •,(h0−z0f ) is the solution of the homogeneous problem (2.1) with initial data h 0 − z 0 f and u (z0f ) is a solution of (2.7) with initial data z 0 f . Then, part (i) and Theorem 2.5 finish the proof of (ii) and (iii).
It is then readily verified using Itô's formula that the unique solution of (2.8) is given by We now focus on the case X = σW for a real Brownian motion W and a constant σ > 0. Then, we will consider regular two-dimensional realizations of the for measurable functions b, a : J → R, where J has non-empty interior, a(y) > 0 for all y ∈ J and 1/a is locally integrable on J.
Examples of such regular two-dimensional realizations are given by Theorem 2.10.(i).
Theorem 2.12. Let X t = σW t , t ≥ 0, for σ > 0 and a real Brownian motion W , and assume that (2.7) admits a regular finite-dimensional realization u t = Φ(t, Y t ), t ≥ 0. Then f is an eigenfunction of A for some eigenvalue λ ∈ R, and there exists an invertible transformation h : J → R + such that for t ≥ 0, almost surely where Z is given by (2.8).
Proof. Let ϕ ∈ dom(A * ), and An application of the Itô formula yields Comparing the martingale term with (2.7), we see that Φ ϕ satisfies the ODE dt ⊗ dP-a. e. and hence, Φ ϕ must be of the form for some g ϕ ∈ C 1 (R ≥0 ) and y 0 in the interior of J. The regularity property of the representation guarantees that h is well-defined and strictly monotone increasing. We stress that h is in fact independent of ϕ ∈ dom(A * ). Setting Z t = h(Y t ), we see that Z satisfies dZ t = m(Z t )dt + σZ t dW t for the drift function m = (bh ) • h −1 + 1 2 (a 2 h ) • h −1 .
Note that for each t ≥ 0, the mapping ϕ → g ϕ (t) is linear continuous from dom(A * ) ⊂ H into R. Since dom(A * ) ⊂ H is dense, by Riesz representation theorem for each t ≥ 0 there exists g(t) ∈ H such that (2.14) g(t), ϕ = g ϕ (t).
By density of dom(A * ) in H this yields that Af = λf , i.e. f must be an eigenfunction of A with eigenvalue λ.
Putting everything together, we have shown that u t = α c f Z t where dZ t = (λZ t + c) dt + σZ t dW t .
Rescaling Z by α c concludes the proof.

Linear SDEs & Pearson diffusions.
Let again X t = σW t for some σ > 0 and a real Brownian motion W . The factor processes Z appearing above are then special cases of the linear SDE studied e.g. in (Kloeden and Platen, 1992, Ch. 4) or (Kallenberg, 2002, Prop. 21.2). Well-known special cases are the geometric Brownian motion (c = d = 0) and the Ornstein-Uhlenbeck-process (b = 0). Relevant in our context is the less common case d = 0, on which we focus now. Using (2.11), the solution is given by Solutions of (2.18) have also been studied in the context of reciprocal gamma diffusions (see e.g. the 'Case 4' in (Forman and Sørensen, 2008)) or also Pearson diffusions. These are generalizations of (2.18) that allow for a square-root term in the diffusion coefficient.
Proposition 2.13. Assume that z 0 > 0, a < 0 and c > 0. Then, Z has unique invariant distribution , which is an Inverse Gamma distribution with shape parameter 1 − 2a b 2 and scale parameter b 2 2c and, for any bounded measurable function φ : (0, ∞) → R, Proof. First, note that define a scale density and speed measure for Z. Then, one can easily verify that Z is strictly positive and recurrent on (0, ∞), see e. g. (Karatzas and Shreve, 1987, Prop. 5.5.22). Moreover, m((0, ∞)) < ∞ and so the unique invariant distribution of Z is .
The remaining results then follow from e. g. (Borodin and Salminen, 2012, II.35) or (Revuz and Yor, 1999, X.3.12). Thus, Remark 2.14. Let a < 0, c > 0 and (Z t ) be the stationary solution of that is, Z 0 is chosen distributed according to inverse gamma distribution with shape parameter 1 − 2a b 2 and scale parameter b 2 2c . Then, as shown in (Bibby et al., 2005), the autocorrelation function of (Z t ) is given by To study price dynamics it is also useful to examine the reciprocal process Y = 1/Z. When d = 0, Y = 1/Z is the unique solution of In particular, with X given in (2.20), When a < b 2 , (2.24) is called the stochastic logistic equation.
2.4. Positivity, stationarity and martingale property. Let us first come back to the linear homogeneous situation. On average, market makers do not accumulate inventory, which suggests to consider the baseline case of balanced order flowfor which X is a (local) martingale. If X is a local martingale with ∆X > −1 a. s., then, from the properties of stochastic exponentials we obtain that: The weak solution u t of the homogeneous equation (2.1) is a local martingale, if and only if the initial condition h 0 is A−harmonic: h 0 ∈ dom(A) and Ah 0 = 0. If E(M ) is a martingale and Ah 0 = 0, then (u t ) t≥0 is a martingale.
In the Brownian motion case, from the discussion in the previous section we directly obtain: Corollary 2.15. Let X = σW where W is a standard Brownian motion and σ > 0, and u be the solution of the inhomogeneous equation (2.7), where f is an eigenfunction of A with eigenvalue −ν and h 0 = z 0 f , for some z 0 > 0. If ν > 0 and α > 0, then where Z ∞ has an Inverse Gamma distribution with shape parameter 1 + 2 ν σ 2 and scale parameter σ 2 2α .
Remark 2.16. The Inverse Gamma distribution has a Pareto (right) tail with tail index 1 + 2ν in this case: the k-th moment of E(Z k ∞ ) < ∞ if and only if k < 1 + 2ν.
So far, we have set aside the positivity constraint for u. By Theorem 2.5 this reduces to analysis of the deterministic equation. In the case of second-order elliptic operators, positivity results from the comparison principle, whenever the initial condition h 0 is positive: Assumption 2.17. Let I ⊂ R be an interval and suppose that A is a uniformly elliptic operator of the form with Dirichlet boundary conditions, and where η, β and α are smooth and bounded coefficients, and in particular η(x) ≥ η > 0 for all x ∈ I.
In addition, the principal eigenvalue of A, λ 1 has an eigenfunction f which is positive on I (Evans, 2010, Sec. 6.5). Note that the factor process Z t has state space (0, ∞) both in Theorem 2.5 and 2.5. We thus obtain the following corollary. This simple result thus guarantees the existence of a solution with the correct sign, thereby avoiding recourse to 'reflected' solutions as in (Hambly et al., 2020) and considerably simplifying the analysis of our model.

A two-factor model
We now study the simplest example of model satisfying Assumption 2.17, namely the case of constant coefficients η a , (3.1) together with the sign condition: In the following, we will write u b In particular the only positive eigenfunction is h 1 .
Proof. First we note that that φ is an eigenfunction of A with eigenvalue ν, if and only if x → e β 2η x φ(x) is an eigenfunction of A 0 := η∆ + α Id with zero Dirichlet boundary conditions, for eigenvalue ν + β 2 4η . Details of calculations are given in (Cont, 2005). The operator A 0 with domain dom(A 0 ) := dom(A) is self-adjoint, has compact resolvent (Cont, 2005) and eigenvalues Eigenfunctions of A 0 with eigenvalue ν ∈ R are solutions of the Sturm-Liouville problem with zero boundary conditions, which yields that g must be of the form The zero boundary conditions at 0 and ±L imply γ 2 = k L π for some k ∈ N so (3.8) ν = α − ηk 2 π 2 L 2 . Translating this from A 0 to A yields the result.
Define the following bilinear forms: which define equivalent inner products, respectively for L 2 (−L, 0) and L 2 (0, L).
For γ > 0, and k ∈ N, define Then (h b k ) k∈N is an orthonormal basis of L 2 (−L, 0), ·, · −γ b and (h a k ) k∈N is an orthonormal basis for L 2 (0, L), ·, · γa and solutions for the SPDE may be constructed using an expansion along these bases: is the unique continuous weak solution of (3.1) in the sense of Definition 2.2.
Proof. The unique continuous solutions of the respective deterministic equations are given by (S b t u b 0 ) t≥0 and (S a t u a 0 ) t≥0 , where (S b t ) t≥0 and (S a t ) t≥0 are the Dirichlet semigroups generated by on (−L, 0) and (0, L), respectively. Thus, from Theorem 2.5 we get (S a t ) and (S b t ) are linear continuous so that for each h a ∈ L 2 (0, L), h b ∈ L 2 (−L, 0), By Proposition 3.1 h a k (resp. h b k ) are eigenfunctions of A a (resp. A b ) and thus also of S a (resp. S b ). This yields the desired representation, where the series converge in L 2 . To obtain pointwise convergence, we note that for x ∈ [0, L] and t > 0, by Cauchy-Schwarz inequality, Parseval's identity and integral criterion for sequences, for ∈ {a, b}, L 2 y 2 dy.
When η > 0, then the weights e −ν k t of the spectral decomposition decay exponentially in k 2 for large k. This justifies approximating the solution by the first few terms. Note also that the only positive eigenfunctions are the principal eigenfunctions h a 1 and −h 1 b so the sign constraints (3.2) only if the projection of the solution along the principal eigenfunctions dominates the other terms in the expansion. This motivates us to focus on solutions which live in the first eigenspace. This occurs if the initial condition is a (positive) linear combination of h a 1 and h 1 b . We will later show that this assumption is supported by market data. This leads to a finite-dimensional realization which satisfies the sign constraints (3.2): The unique solution of (3.1)-(3.2) with initial condition The L 1 normalization (3.18) allows to interpret the variables in terms of order book volume and depth: represents the volume of sell (resp. buy) orders, while ∇u t (0+)θ = θπ L V a t (resp. ∇u t (0−).θ = θπ L V b t ) represents the depth at the top of the book. In this simple two-factor model, these two are proportional to each other: they may be decoupled by considering multifactor specifications involving higher-order eigenfunctions.
The drift parameter −ν a (resp. −ν b ) thus represents the net growth rate of decrease of the volume of sell (resp. buy) orders. As shown in (3.20), this net growth rate results from the superposition of several effects: submission/ cancellation of limit sell (resp. buy) orders by directional sellers (resp. buyers) at rate α a (resp. α b ); this may be interpreted as the 'low frequency' component of the order flow; replacement of limit orders by new ones closer to the mid-price, at rate β 2 a 4ηa (resp. β 2 b 4η b ); cancellation of limit orders as the mid-price moves away (i.e. at distance ±L from the mid-price), at rate ηaπ 2 L 2 (resp. η b π 2 L 2 ).

20
In the case of a balanced order flow for which there is no systematic accumulation or depletion of limit orders away from the mid-price, these terms compensate each other and the volume of limit orders in any interval [S t + x 1 , S t + x 2 ] is a (local) martingale. The following result follows from the remarks in Section 2.4: Corollary 3.4 (Balanced order flow). The order book density u is a local martingale (in L 2 ), if and only if Remark 3.5 (Balance between high-and low-frequency order flow). The balance condition (3.22) expresses a balance between the slow arrival of directional orders, represented by the terms α a and α b , and the fast replacement of orders inside the book, represented by the terms β 2 a 4η b and β 2 a 4η b , and finally the cancellation of limit orders deep inside the book, at rate η a π 2 /L 2 .
This balance between order flow at various frequencies may be seen as a mathematical counterpart of the observations made by (Kirilenko et al., 2017) on the nature of intraday order flow.
3.2. Shape of the order book. An implication of the above results is that the average profile of the order book is given, up to a constant, by the principal eigenfunctions H a 1 , H b 1 : H a 1 (x) Dropping the indices a, b, the normalized profile of the order book has the form: where c 1 is such that Lβ .
The position of the maximum moves closer to the origin as β /η is increased. For β = 0 we havex = L 2 , and, on the other handx 0 as β /η → ∞. Typically, the order book profile for liquid large-tick securities a few ticks from the mid price. Figure 4 shows the average order book profile for QQQ; similar results were found in (Bouchaud et al., 2009;Cont et al., 2010). This suggestsx is of the order of a few ticks, so we are interested in the parameter range for which β /η is large.
The value at the maximum is which grows linearly as β /2η → ∞, as shown in Figure 3, where we have plotted h, normalized by its L 1 -norm, for various values of β with L := 3π and η = 1. The above results are valuable for calibrating the model parameters β 2η , α and σ to reproduce the average profile (for each side) of the order book. The height of this maximum gives a further constraint on parameters, using (3.25). We will use this result for parameter estimation in Section 3.6.
3.3. Dynamics of order book volume. As noted in Corollary 3.3, V a t and V b t may be identified as the volume of sell (resp. buy) limit orders: they follow (correlated) geometric Brownian motions: for p ∈ (S t , S t + L), and in the following sense: (v t ) t≥0 is an continuous L 2 (R)-valued stochastic process and for all ϕ ∈ C ∞ 0 (R) and t ≥ 0, for x, y , y , y ∈ R.
Remark 3.8. Note that (3.36) is a stochastic boundary condition at S t .
3.6. Parameter estimation. We now describe a method for estimating model parameters. We use time series of order books for NASDAQ stocks and ETFs, from the LOBSTER database.
Given that we do not observe separately the various components of the order flow as in (3.1), we use the relations discussed in Sec. 3.2 to calibrate the parameters σ, ν and the shape parameter for each side of the order book. We set L to the largest value in our data set, (L := 1000). Parameters may be calibrated either through (a) a least squares fit of (3.23) to the average order book profile, or (b) calibrating parameters to reproduce the positionx and height of the maximum of the order book profile.
Remark 3.9. The estimator based on the maximum position of the peak is fast in computation but the fixed price level grid in the data restricts the set possible values for estimation of γ. In particular, the estimator is sensitive to the location of the maximum (i.e. the mode of the order book profile).  We show results for a set of NASDAQ stocks and ETFs. Figure 4 shows how the model reproduces the average book profile for QQQ at NASDAQ on 17th November 2017. In Figure 5 we see the coefficient γ estimated across various 30-min windows during the trading day. The one-factor model based on the principal eigenfunction yields a reasonable approximation for the average order book profile, which justifies our assumptions on the dynamics in Section 1.2. For low-price/large tick stocks, the average order book profiles may differ from the exponential-sine shape. For such stocks, we use the nonlinear scaling described in Section 1.1, leading to an average order book profile: where S t is the best price. Figure 6 shows such a nonlinear fit for the average order book profile of SIRI.
The solution of the SPDE may then be obtained as follows: Proposition 4.1.
(i) The unique L 2 -continuous solution of (1.2) -(4.1) for a general initial condition u 0 is given by (ii) For an initial condition of the form We obtain the general solution of the linear homogeneous equation from Proposition 3.2. The series representation of u results from the spectral decomposition, Proposition 3.1 and Theorem 2.10.

4.2.
Long time asymptotics and stationary solutions. In order to derive properties of the 'average' order book profile, we now examine whether the order book profile u t has an ergodic behavior and describe stationary solutions. The following result describes the long-term dynamics and shows that this dynamics is well approximated by projecting the initial condition on the principal eigenfunctions as done in (4.3): Proposition 4.2. Let u t be the unique solution of (1.2) -(4.1) for a general initial condition u 0 ∈ L 2 (−L, L) and define: If ν b 1 > 0 and ν a 1 > 0, then: (i) The long-term dynamics of the order book is well approximated by the dynamics (4.4) projected along the principal eigenfunctions: (4.5) P lim t→∞ u t −ǔ t ∞ = 0 = 1.
(ii) u t has a unique stationary distribution and where f a , f b are given by (4.1) and Z a (resp. Z b ) is an Inverse Gamma random variable with shape parameter 1 + 2 νa ) and scale parameter σ 2 a 2Va (resp.
Proof. For t 0 > 0, let This term is indeed finite by integral criterion for series, see e. g. proof of Proposition 3.2. Denote u • t (.; h) the unique solution of the linear homogeneous equation (3.1) for an initial condition h. Recall from Theorem 2.10 that . It suffices now to prove the results for the ask side and note that the calculations will be analogous for the bid side. Using the representation of u • t from Proposition 3.2 we get for all t > t 0 and all h ∈ L 2 (0, L), which, as t → ∞, converges to 0 provided that ν a 1 > 0. This proves (i). To show (iii), a similar calculation but using the orthogonality of the decomposition in Proposition 3.2 yields .
If σ 2 a < 2ν a 1 , then this converges to 0 as t → ∞. Since . βa 2ηa defines an equivalent norm on L 2 (0, L), this finishes the proof of (iii). Assertion (ii) follows from Proposition 2.13. Indeed, recall that V i , i ∈ {a, b} are ergodic processes whose unique invariant distribution is given by an Inverse Gamma distribution with shape parameter 1 + 2νi σ 2 i and scale parameters σ 2 i (Vi) 2 , i ∈ {a, b}. Denote by Z b and Z a random variables with these distribution For any x ∈ [−L, L], we have the convergence in distribution =⇒ Z a f a (.). Since almost sure convergence yields convergence in distribution, by part (i) this yields that (4.8) holds also for u t with arbitrary initial data u 0 ∈ L 2 (−L, L).

4.3.
Dynamics of order book volume. Consider now the 'projected' dynamics as in the setting of Proposition 4.1.(ii). The dynamics of the order book volume V t is then given by where V b and V a , defined in (4.2), represent the volume of buy (resp. sell) orders in the order book.
for some Brownian motion W , independent of W . Then, In particular, the quadratic variation ('realized variance') of the order book volume is given by dt For the symmetric and perfectly correlated case, V is itself a reciprocal gamma diffusion: Corollary 4.3. Assume the setting of Proposition 4.1.(ii) and, in addition, that ν a = ν b =: ν, σ a = σ b =: σ and a,b = 1. Then, V is the unique solution of (4.12) In all cases, we get from (2.22) that for i ∈ {a, b}, t ≥ 0, (4.13) and (4.14) Over a small time interval ∆t, where N 0,1 is a standard Gaussian variable. In particular the conditional probability of an upward mid-price move of size y is given by (4.17) Combining Proposition 2.13 and Remark 2.16 with (Leonenko andŠuvak, 2010, Theorem 6.3) we obtain that if D > 0 and c > 5, then V has finite 4th moment and the estimators are consistent and asymptotically normal. For the autocorrelation parameters ν a and ν b we use the martingale estimation function (Bibby and Sørensen, 1995 Convergence of this estimator is discussed in (Bibby and Sørensen, 1995, Theorem 3.2). We apply these estimators to high-frequency limit order book time series for NASDAQ stocks and ETFs, obtained from the LOBSTER database, arranged into equally spaced observations over time intervals of size ∆t = 10ms and dt = 50ms. For each observation we use as market depth the average volume of order in the first two price levels, respectively on bid and ask side. 3 Below we show sample results for ETFs (SPY and QQQ) and liquid stocks,(MSFT and INTC).
Figure Table 1 shows estimated parameter values across different days for INTC, MSFT, QQQ and SPY. We observe negative values of correlation a,b across bid and ask order flows which is consistent with observations in (Carmona and Webster, 2013). Figures 7 and 8 show intraday variation of estimators for ν a , ν b , σ a , σ b and a,b computed over 15-minute windows.
There are various estimators for intraday price volatility in this model, which allows to test the model. Recall that in (4.16) we expressed price volatility in terms of the parameters describing the order flow: (4.25)σ S := θ σ 2 b + σ 2 a − 2σ b σ a a,b . where θ is the impact coefficient. We call this the RV estimator.
Another estimator is obtained by first estimating σ b and σ a using the martingale estimation function (4.22) then computing the price volatility using Equation (4.25). We label this the RCG estimator.
Finally, one can compute the realized variance of the price over a 30 minute time window using price changes over 10 ms intervals. Comparing these different estimators is a qualitative test of the model. Figure 9 compares these estimators, computed over 30 minute time windows: we observe that the model-based estimators are of the same order and closely track the intraday realized price volatility, which shows that the model captures correctly the qualitative relation between order flow and volatility.
3 The source code for the implementation is available online (Cont and Mueller, 2018). . Comparison of various estimators for intraday price volatility σ S : standard deviation of price changes (blue), estimator based on realized variance/covariance of bid/ask depth (red), and estimator based on martingale estimation function (orange).