Anthony Giroud
email: giroud@univ-umlv.fr

Venceslas Biri
email: biri@univ-umlv.fr

Modeling and rendering heterogeneous fog in real-time using B-Spline wavelets

Keywords: Participating medium, Fog, Rendering, GPU

Heterogeneous fogs are often modeled with several layers of different density or using particle systems. However, layers are limited to vertical variations and using particles can involve long computation time with large outdoor scenes. In this article we present a simple method to render heterogeneous fog in real-time. The extinction function of our fog, related to its density, is first modeled in a B-Spline function basis. A wavelet transform is applied on this function to obtain a decomposition in both space and frequency domains. A grid traversal is used to render the fog in real time using GPU. Since no precomputation is required concerning the position of the camera or the fog, we can freely navigate or move the fog into the scene.

INTRODUCTION

Fog is massively used in rendering both for aesthetic purposes and to increase performances by providing an efficient way to cull surfaces that are far from the camera. Simple fog models, are straightforward to implement but, like OpenGL's fog model, only allow a basic representation of homogeneous fog as can be seen on figure 1. Most of the time, these models are barely convincing visually, as we know that natural fogs never reach such perfect homogeneity. Considering latest advances in GPU programming, design of heterogeneous fog should be simple, and its rendering reachable in real-time.

The fog phenomenon is due to small particles of water in suspension. Because it interacts with light rays, fog is considered as a participating medium in computer graphics. Fog effects take into account attenuation, caused by absorption and out-scattering, and also consider multiple scattering of light as isotropic and constant over the scene. If we consider an homogeneous fog in its simplest form, equations are simple enough to allow an analytical integration of its effects along a view ray. When rendering heterogeneous fog, the density of water particles is varying across the scene, thus dramatically complexifying the model, involving local changes in physical properties of the fog, such as its extinction coefficient. Therefore, in order to compute the light-fog interaction, we have no other choice than per-Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Considerable work has been achieved in the development of real time solutions to handle participating media. Physical simulations taken aside [START_REF] Rushmeier | The zonal method for calculating light intensities in the presence of participating medium[END_REF][START_REF] Sillion | A Unified Hierarchical Algorithm for Global Illumination with Scattering Volumes and Object Clusters[END_REF][START_REF] Jensen | Efficient Simulation of Light Transport in Scenes with Participating Media using Photon Maps[END_REF][START_REF] Pauly | Metropolis light transport for participating media[END_REF][START_REF] Jarosz | Irradiance Gradients in the Presence of Participating Media and Occlusions[END_REF], which do not reach realtime, researchers have been working on rendering complex exchanges of light within the medium, dealing, for example, with single scattering. They also considered simpler forms of fogs, with a density either varying along horizontal layers, defined by Perlin noise, or using particles. But few tried a direct and continuous mathematical representation of its density.

In this paper, we present a new method helping to shape and render complex heterogeneous fog in large outdoor scenes, lighted by a single light source (the sun). First, the fog is modeled in a B-Spline function basis, which allows a simple and efficient construction of its extinction function. As a preparation before rendering, Mallat's wavelet decomposition is applied on the extinction function in order to automatically generate different resolutions, enabling an optimized realtime rendering using the GPU. The use of wavelets offers several advantages :

• An easy modelization leading to a smooth and continuous fog density by opposition to particles approaches that are discrete. Analytical representation compresses data more efficiently and are, for example, easier to animate.

• Wavelet modelization is generic. It includes naturally, using Haar wavelets, discrete approaches like quad tree or octree representation.

• Wavelet decomposition leads to sparse data that can be used to improve rendering time.

Therefore the contribution of this paper is :

• Establishing a wavelet framework for the definition and modelization of an heterogeneous fog.

• Rendering the fog in real time using this representation without precomputation involving camera or fog position.

• Allowing a tradeoff between correctness and speed using the multiresolution offered by the wavelet decomposition.

In the next section, we review previous methods to render, in real time, the effects of participating media in a scene. Then, we briefly introduce to the wavelet theory along with the equation of transfer inside a participating medium. Section 4 presents our modelling scheme and our implementation for rendering. In section 5, we expose and discuss our results.

PREVIOUS WORK

Rendering participating media such as fog in real-time has been well studied. We will not consider global illumination algorithms concerning participating media. For more information on this subject, the readers should refer to the excellent survey of Cerezo et al. [START_REF] Cerezo | A survey on participating media rendering techniques[END_REF]. Algorithms dealing with single scattering, including volume based approaches [START_REF] Sun | A practical analytic single scattering model for real time rendering[END_REF] or direct representation [START_REF] Biri | Real Time Single Scattering Effects[END_REF], also handle fog naturally but due to complexity problems these techniques only consider homogeneous mediums (except [START_REF] Zhou | Fogshop: Real-time design and rendering of inhomogeneous, single-scattering media[END_REF] discussed bellow). Therefore, we limit our overview to other real time approaches for heterogeneous fog which can roughly be divided in, on the one hand, particle approaches and, on the other hand, layered or bounded approaches.

Particles provide a natural way to handle heterogeneous fog. They have been used efficiently in numerous works [START_REF] Foster | Modeling the motion of a Hot, Turbulent Gas[END_REF][START_REF] Stam | Stable Fluids[END_REF][START_REF] Adabala | Modeling and rendering of gaseous phenomena using particle maps[END_REF][START_REF] Fedkiw | Visual Simulation of Smoke[END_REF]. The idea is to consider particles as groups of water drops, allowing real time rendering of effect like smoke or physically based simulation. But is not well adapted to large scale fog recovering a whole scene. Moreover, animation of all particles in a large scene is computer time consuming. The same drawbacks hold for the hybrid approach of Zhou et al. [START_REF] Zhou | Fogshop: Real-time design and rendering of inhomogeneous, single-scattering media[END_REF] which handles single scattering in a heterogeneous participating medium combining particles and spherical harmonics. We can also cite the work of Zdrojewska [START_REF] Zdrojewska | Real time rendering of heterogeneous fog based on the graphics hardware acceleration[END_REF] which uses Perlin noise to alter the homogeneous density of the fog. Despite this good idea, the use of 3D random noise forbids any animation of this fog.

The idea behind layered or bounded approaches is to enclose fog density variations into layers [START_REF] Legakis | Fast multi-layer fog[END_REF][START_REF] Heidrich | Applications of pixel textures in visualization and realistic image synthesis[END_REF] or bounded volumes [START_REF] Nvidia | Fog polygon volumes -rendering objects as thick volumes[END_REF]. These works consider homogeneous fog enclosed in volume, inducing a discontinuous density function and creating artifacts on the border of these volumes. Moreover, intuitive or physically based animations of this kind of representation could be difficult to handle. Despite these limitations, it is often the kind of solution we can find in common graphic engines, along with particle rendering. Nevertheless, none of the previous methods offers a simple and efficient mathematical representation of heterogeneous fog adapted for both animation and rendering.

THEORETICAL BACKGROUND

Fog's illumination model

Our main goal is to render our fog in real-time, using conventional graphics cards. Although performances of GPUs have never been increasing so fast, we have to slightly simplify our fog model. Between points O and P, fog induces an attenuation (due to out-scattering and absorption) of the luminance L of P and an increase (in-scattering and emission) of light along the ray OP. We start directly with the integral transfer equation, see [START_REF] Siegel | Thermal Radiation Heat Transfert[END_REF] :

L(O) = τ(O, P)L(P) + P O τ(O, u)K t (u)J(u, ω)du (1)
L(O) being the radiance received by the observer, J(u, ω) being the incoming radiance along the ray, K t the extinction coefficient and τ(u, v) the transmittance of the fog along the ray going from u to v : First, when daylight passes through fog, it is immediately scattered such that light in-scattering can be simplified by a constant amount L fog . Moreover, if we con-sider that the light emitted by the fog itself can be neglected, the incoming radiance J(u, ω) equals L fog , and then equation (1) becomes :

τ(u, v) = e -v u K t (s)ds (2)
L(O) = τ(O, P)L(P) + P O τ(O, u)K t (u)L fog du (3)
The second part of equation ( 3) can be analytically integrated to obtain :

L(O) = τ(O, P)L(P) + L fog (1 -τ(O, P)) (4)

Wavelets

From equation ( 3), we can see that the density variation could be represented efficiently by the extinction function. Therefore, K t will be modeled using the wavelet framework whom principal used characteristics are detailed in this section. More details on the wavelet framework can be found in [START_REF] Stollnitz | Wavelets for computer graphics: a primer.1[END_REF].

The wavelet framework In a multiresolution analysis, data is represented using several approximation spaces. Different functions bases are used to represent a single signal, and each functions basis corresponds to a different resolution. Moreover, all basis functions are obtained by translating and scaling a single original pattern function f ∈ L 2 (R), in other words :

f j,k (x) = f (2 j x -k), with j ∈ N, k ∈ Z (5)
where f j,k represents the basis functions and j the resolution level. If we define F j as the closed subspace of L 2 using basis functions { f j,k } k∈N , the closure of j∈N F j is the space L 2 and represent all square integrable functions.

The wavelet framework uses, to build basis functions of spaces F j , a particular function called scaling function and often denoted by φ . It verifies equation ( 5) and generates a φ jk family, j ∈ N, k ∈ Z. This function φ also presents the property of being written as a linear combination of k/2 translated and 1/2 scaled versions of itself. It is the scaling relation of the scaling function, given by :

φ (x) = ∞ ∑ k=-∞ p k × φ (2x -k) ( 6
)
where {p k } are the coefficients of the scaling sequence of φ . Note that each subspace F j , j ∈ N will in fact use the same and unique function φ translated and scaled.

The particularity of the wavelet framework is its ability to decompose a function of F j+1 using several functions of F j and of its orthogonal complement G j . Therefore, if J is the maximum resolution level, the F J space can be written :

F J = F 0 ∪ J-1 j=0 G j (7)
This equation means that a function (up to a resolution J) can be described using only one scaling function (space F 0 ) and several functions of spaces G j . The basis functions of spaces G j are called wavelet function and verify equation [START_REF] Heidrich | Applications of pixel textures in visualization and realistic image synthesis[END_REF]. They can also be built using the scaling relation for wavelets, which we will call the wavelet relation :

ψ(x) = ∞ ∑ k=-∞ q k × φ (2x -k) ( 8
)
where {q k } are the coefficients of the wavelet sequence of ψ. Note that, similarly to F j , each subspace G j uses the same and unique function ψ translated and scaled.

Decomposition and multiresolution using wavelet

The advantage of the wavelet framework is that it provides an efficient way to decompose a function into multiresolution spaces. The fast decomposition can be assured by the Mallat's wavelet transform which uses, as entry data, coefficients of the function modeled directly in the maximum resolution level. Therefore, our fog extinction function will be modeled using scaling function. Mallat's algorithm takes advantage of equation ( 7) and consists, for each step, in extracting from the approximation at level n (represented in a scaling functions basis) first the approximation at level n -1 (F n-1 which is twice less precise), and then the corresponding layer of details (G n-1 represented by a wavelet basis). We simply repeat this process until we obtain the approximation at level 0. Mallat's transform is lossless, therefore when we simply sum up the coarsest approximation with all layers of details, we recover the original signal untouched. Wavelets in two dimensions Now that we know how to build scaling functions and wavelets in one dimension, going 2D will actually be quite straightforward. In a nutshell, it simply consists in assigning the corresponding 1D function to each axis, and the result is given by the product of these two 1D functions. Basically :

φ φ (x, y) = φ (x)φ (y) (9)
where φ φ is a 2D scaling function and φ is the corresponding 1D scaling function. Things go exactly the same way with wavelet functions. Obtaining a 2D wavelet transform is slightly harder and requires to process rows and columns separately. There are two different decomposition methods : the standard decomposition and the nonstandard decomposition. These two types of decomposition output exactly the same kind of result :

• A single coarse approximation at level 0, modeled with 2D scaling functions φ φ (x, y) = φ (x)φ (y).

• J -1 layers of vertical details, modeled with hybrid functions φ ψ(x, y) = φ (x)ψ(y).

• J -1 layers of horizontal details, modeled with hybrid functions ψφ (x, y) = ψ(x)φ (y).

• J -1 layers of 2D details, modeled with 2D wavelets ψψ(x, y) = ψ(x)ψ(y).

For example, our fog's extinction function can be written as :

K t = ∑ i j α i j φ φ i j + (10) J-1 ∑ n=1 ∑ i j β n i j φ ψ n i j + δ n i j ψφ n i j + γ n i j ψψ n i j 4 OUR METHOD 4.

Modeling the fog

The two-dimensional framework Unlike other types of participating media from the same family, fog almost always appears in large outdoor scenes as a horizontal layer of varying thickness. This is quite different from smoke, which can evolve indifferently in all directions in terms of shape and movement, and thus really need to be defined with the same precision along all three dimensions. For this reason, and in order to ease the shape definition as much as possible and, later, the rendering step, we have chosen to restrict our main framework to two dimensions. The optical properties of our fog, similarly to most other participating media rendering techniques, are proportional to its density, which depends itself on its extinction function. Therefore, the fog's main shape will actually be modeled as horizontal layers containing horizontal extinction function projected in a two-dimensional function basis. Further parameters, starting with a vertical extinction coefficient, will then thicken the fog vertically and give its final appearance.

Designing the fog's shape The horizontal variations of our fog's density are modeled by specifying the value of each coefficient in the extinction function basis. These coefficients can be adjusted by hand, or be, for example, the result of a simulation, which was exported as a fogmap (see figure 3), i.e. a greyscale image, and then loaded back in our implementation.

Compared with other techniques such as RBF or particle-based methods, shaping our fog using a grayscale image is straightforward. The fogmap represents, in some extent, a direct preview of its aspect, what can be interesting for some applications where great intuition is needed. To ease the manual setting of the coefficients, we also developed a small application where the values of the density can be directly adjusted using a drag-and-drop interface. Choosing the basis functions The appearance of the fog's density is a key criteria to choose our basis function. It is clear that abrupt changes in density would not look natural, so we would ideally like continuous functions to design smooth fogs using as few coefficients as possible. In order to avoid border effects, the scaling function must tend to zero on both sides of its support, which eliminates, for example, Legendre scaling functions.

For design and optimisation purposes, our rendering algorithm also needs the scaling function never to oscillate under zero. Whatever the trajectory of the ray within the function in 2D, and more generally within the fog, we would like to be sure that the sum of the density it intersects can only increase as it traverses the fog from the observer to the nearest object. Daubechies wavelets, which, by the way, are not symmetrical, might not be the way to go.

Finally, we have to consider the fact that, as will be discussed in the next section, the cost of using a particular type of wavelet is quadratically proportional to the support of the scaling function in one dimension.

According to their shape, the most adapted candidates seem the linear or quadratic B-Splines, which are shaped like a hill (see figure 4), and have a relatively compact support.

Although we are limited to wavelet scaling functions for the fog's representation, our method is not reduced to a particular type of wavelet. Our implementation specifically handles all degrees of B-Spline wavelets, but can be extended to other families, as long as they are compatible with Mallat's decomposition.

Preparing data for rendering

Generating multiple resolutions One of our main goals is to take profit of multiresolution. Indeed, multiresolution helps to omit details that could be expensive to render, while being of limited visual importance. Therefore, perform a wavelet decomposition on our fog, which generate multiple level of details (i.e. multiple resolutions) from the original extinction function, and use them at the rendering phase. The most adapted solution seems Mallat's fast wavelet transform, which is lossless, but requires data to be modeled in a scaling functions basis of the same type as the wavelets used for the decomposition. Therefore, each pixel of the fogmap will represent the coefficient of a scaling wavelet function.

Computing textures

From the fogmap we generate four multiple-level function bases : the approximation on a single level (i.e. a single 2D grid of values), and three different kinds of details for each level which was decomposed. All details bases have the same depth, which corresponds to the number of decomposition steps that were executed, value which must be decided by the user, depending on how much details can be omitted. Coefficients from the approximation and details basis will be stored in packed textures, and transmitted to the GPU under this form.

Rendering the fog

Overview The purpose of our algorithm is to alter the original color of each pixel of the image using equation (3), blending L(P), the color of the object behind the fog, and the fog color to obtain L(O) the new color to compute.

For each pixel, we perform a ray-marching from the camera to the nearest surface, in which we integrate over the fog's extinction function to obtain the transmittance τ(O, P) along the view ray OP. The grid As a result from the wavelet decomposition, the fog's density is scattered in several multiple-level function bases, having their own vector space and definition domain in 2D. Each single level can be assimilated to a rectangular grid, each cell being associated to both a coefficient and a basis function. Since all bases have the same definition domain, grids from different bases match at a given level.

Since our fog is only modeled in two dimensions, we do not take into account vertical variations and consider the fog as homogeneous on that direction. However, a vertical extinction coefficient taken as parameter allows to fade the fog out while its vertical distance from the viewer increases. But note that this is only a quick approximation over the exact equations. Integration along the ray The algorithm is iterative, but instead of advancing regularly along the ray, we move cell by cell. Each step corresponds to a new intersection between the ray and the grid, thus we always integrate between two intersections, i.e. between two positions on the perimeter of a square cell. This is a brute-force method, and some optimisations will be discussed in the next section.

We start by transposing both positions of the camera and the object from the scene to the fog's vector space. Our algorithm performs the entire integration level by level, and then, for each single function basis level, cell by cell.

To initiate the integration on a given level, we first determine both entry and exit points of our integration on the grid. The entry point corresponds to either the nearest intersection between the ray and the current level's bounding box, or the viewer's position in case he stands within the fog. Similarly, the exit point corresponds to the intersection with either the farthest plane of the bounding box, or with the nearest object if situated within the fog.

When integrating a given level, the contribution of each single cell can be obtained by the product of both the function basis coefficient and the integral of the basis function associated to that cell along the view ray.

Mathematically, considering each cell c intersected by OP and using the extinction function decomposition of (10), we have :

τ(O, P) = ∑ cell:c c∩OP K t = ∑ c c∩OP α c φ φ c + (11) J-1 ∑ n=0 c∩OP β n c φ ψ n c + δ n c ψφ n c + γ n c ψψ n c
J being the maximum decomposition level of our fog. Thanks to multiresolution analysis, each function in-dexed by cell c and level n is indeed a translated and scaled version of φ φ , φ ψ, ψφ or ψψ. Therefore, we can precompute on the CPU a bunch of integrals for a set of sampled paths (complete or partial) within 1 × 1 squares on each function's definition domain, so that these values are directly available at runtime, transmitted on the GPU in packed textures. Integration on partial paths allow handling particular cases when the ray either starts and/or ends at the center of a cell within the fog's bounding box.

Figure 5 shows ray OP traversing a single level's grid from entry point S to exit point E. Integration steps (i.e. intersections with the grid) are shown in red. The basis function (in this example : linear B-Spline scaling function) associated to the orange cell's coefficient c is shown in blue.

When using basis functions which are supported on an 1 × 1 square (e.g. Haar scaling functions and wavelets), their contribution area matches exactly that of the cell it is attached to, therefore we know that the cells which contribute to the pixel being rendered are exactly those traversed by the ray.

When the functions are supported on a domain larger than 1 × 1, part of the contribution of each cell gets superimposed on that of its neighbouring cells, thus also contributing to rays which do not necessarily pass through those cells themselves. Actually, a ray passing through a cell must take into account the contribution of that cell, plus the contributions of the dx -1 previous cells on the X axis, times the dy -1 previous cells on the Y axis, where dx and dy are the dimensions of the basis function's definition domain.

When the ray encounters a new function, we only integrate the density on the portion of that function which overlays the current cell, and then resume the integration for another 1 × 1 square of the same function when the ray traverses the next cell. If we directly integrate on the whole function's support at once, we omit the contributions of the functions attached to cells which are not encountered by the ray.

When the ending point has been reached, the whole process must be repeated with each level of each wavelet basis that was generated by Mallat's wavelet transform.

Optimizations & multiresolution

Our idea consists in omitting an increasing quantity of details from layers whom resolution is above a threshold which decreases as the observer moves away from the fog. When integrating the fog's density from the observer O to point P, the maximum integration distance d max l on level l ∈ N is given by :

d max l ( OP) = OP × µ l ( 12
)
where µ ∈ [0, 1] is the optimization coefficient. When µ = 1, the integration is performed entirely on all levels ; on the contrary, when µ = 0, only the upper level of the basis is rendered.

As seen previously, when using scaling functions that are defined on more than an 1×1 square, the integration cost is no longer proportional to the fog's size, since more than each single particular cell traversed by the ray brings a contribution on these cell's area. That's why although the total number of coefficients modeling the fog stays almost unchanged, the rendering cost increases dramatically after the wavelet decomposition, since B-Spline wavelets always have a larger support than their scaling function.

When using such basis functions, for example linear or quadratic B-Splines, it can be interesting to use the two-scale relation for wavelets 8 to deconstruct the three wavelet bases. This turns them back into scaling function bases, which can then be merged (i.e. added) together. When using scaling functions with a large support, this operation, performed on the CPU just after the decomposition, can reduce the rendering cost by up to 2, while keeping the multi-resolution aspect brought by the decomposition. Moreover, if we perform a deconstruction, we can stop the integration as soon as the sum reaches a particular threshold, close to a great opacity. Deconstruction is important since it assures than each new cell will only add opacity.

RESULTS AND DISCUSSION

This algorithm has been implemented using GLSL, an Intel Core 2 Quad 2.8Ghz processor and a NVidia GeForce GTX 280 graphics card. Screen resolution is 800x600.

Performance

Table 1 show FPS results obtained when using our raymarching alone to directly render raw Haar, linear or quadratic fogmaps, without any decomposition. Each type of basis functions is defined on an area which size is increasing linearly in 1D, which involves a quadratically increasing number of neighbouring cells contributing to the density on each 1 × 1 square on the grid.

Table 2 show FPS results obtained when rendering a 64 × 64 linear B-Spline fog using our details dropping optimization, for different values of the tolerance parameter µ. With µ = 1, no details are dropped, and we are performing a simple ray-marching. If, in addition, we do not apply any decomposition step, we are 3: FPS improvement when turning back into scaling function bases the four b-spline/wavelet generated by the decomposition (before → after). directly rendering the fogmap, like in table 1, therefore this value stands for the threshold above which we have a substential acceleration.

❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳

Visual quality

The higher the degree of the B-Spline wavelet is, the smoother each basis function looks. With Haar wavelets, we can see, in figure 6.B, that the visual result is a bit unsatisfactory, with abrupt changes in density which betray the discontinuity of Haar functions. With linear B-Spline wavelets (figure 6.C) the framerate decreases but the visual result is a lot smoother and artifacts and peaks are now practically imperceptible. Finally, with quadratic B-Spline wavelets (figure 6.D), we loose in performance but this time, the quality gain is relatively low compared to linear B-Spline wavelets.

Discussion

Linear B-Spline seems a good trade-off between speed and quality but Haar could be used if rendering time is an issue. The advantage of using wavelets, beside their property of good data compression, is to have a mathematical representation of heterogeneous fog from physical simulation to rendering. Indeed, animating such fogs is easy, since wavelet decomposition can be performed in real-time. Moreover, unless previous approaches, we perform a precise numerical integration of density along the view ray, without any approximation. In comparison to particle approaches like [START_REF] Zhou | Fogshop: Real-time design and rendering of inhomogeneous, single-scattering media[END_REF], our method is more adapted to large outdoor scenes when camera is moving in the fog, and the modelling is far more intuitive than using particles.

CONCLUSION AND FUTURE WORK

In this paper, we presented a new method for modeling heterogeneous fog using wavelet scaling functions.

Rendering is performed through a simple decomposition scheme of the fog density function represented in a scaling function basis leading to sparse data. Wavelets and scaling functions allow and ease a certain number of precomputations, such as the integrals of the wavelets along each ray. A brute force rendering algorithm using the GPU has been presented allowing real-time rendering for moderated complex fog along with an optimized version taking profit of the sparcity of data induced by the wavelet decomposition. We have shown that our method outperforms brute force integration and allows exact computation of the effects of fog, without exotic approximations. Moreover, our method do not depends on the position of either the light or the fog, allowing simple transformations of the fog. The use of wavelets opens the door to other major optimisations for our method. Mainly, the rendering algorithm can be improved by focusing only on the grid's cells which actually contain a non-negligible value, in order to be able to directly jump to the interesting zones of the fog when performing the integration along the ray. For this purpose, we aim at designing a simple GPU traversal of the graph generated by the wavelet decomposition. Since wavelets can be used to solve fluids equations, we also plan to link our rendering algorithm to a physical simulation involving wavelets, allowing a real-time physical animation and rendering of heterogeneous fog. Finally, we plan to add single scattering and volumetric shadows in our model.

Figure 1 :

 1 Figure 1: Example of heterogeneous fog

Figure 2 :

 2 Figure 2: Ray of incoming light from P to O through a participating medium.

Figure 3 :

 3 Figure 3: Left : snapshot of our modeling tool. Right : greyscale image representing highest resolution coefficients

Figure 4 :

 4 Figure 4: Shape of Haar, Linear and Quadratic B-Spline

Figure 5 :

 5 Figure 5: Ray-marching through a single level, designed with linear B-Splines scaling functions (sup-port=2).

Algorithm 1

 1 Pseudo code of the shader for each pixel do sum = 0 for l = 0 to nb_levels do compute 2D entry point on grid compute 2D exit point on grid while pos = exit do inter = compute next intersection with grid if (l = 0) then coef = get cell coef on approx basis approx = integrate on φ φ between pos & inter sum += coef*approx end if coef = get cell coef on details1 basis det1 = integrate on φ ψ between pos & inter sum += coef*det1 coef = get cell coef on details2 basis det2 = integrate on ψφ between pos & inter sum += coef*det2 coef = get cell coef on details3 basis det3 = integrate on ψψ between pos & inter sum += coef*det3 pos = inter end while end for pixel color=sum*obj color + (1-sum)*fog color end for

Figure 6 :

 6 Figure 6: Quality difference with a large 30x30 fog. Fog taken from above (A), and the associated fogmap (E). Zoom on the red part when using Haar (B), Linear (C) and Quadratic (D) wavelets.

Figure 7 :

 7 Figure 7: Quality difference when removing a layer of details. Above : Haar fog with a resolution of 32 × 32. Below : Same fog, minus the lower layer of details.

Table 1 :

 1 FPS results with our ray-marching without optimizations.

	Dimensions	Type Haar Linear Quadratic
	16 × 16			199	142	71
	32 × 32			124	83		31
	64 × 64			90	45		15
	P P P P P P P P Nb levels µ	1 0.8 0.6 0.4 0.2	0
	0		45	-	-	-	-	-
	1		35 39 47 55 66	83
	2		31 45 55 71 90 124
	3		27 39 66 76 99 166

Table 2 :

 2 FPS results with our optimization, using a 64× 64 linear B-Spline fogmap.

	Nb levels	1	2	3
	Linear 16x16	58 → 111 47 → 99 35 → 83
	Linear 32x32	20 → 62 15 → 55 13 → 49
	Quadratic 32x32	7 → 23	5 → 20	5 → 18
	Table


Rama Cont

Adrien De Larrard

Price dynamics in a Markovian limit order market

Keywords: limit order book, market microstructure, queueing, diffusion limit, high-frequency data, liquidity, duration analysis, point process

We propose and study a simple stochastic model for the dynamics of a limit order book, in which arrivals of market order, limit orders and order cancelations are described in terms of a Markovian queueing system. Through its analytical tractability, the model allows to obtain analytical expressions for various quantities of interest such as the distribution of the duration between price changes, the distribution and autocorrelation of price changes, and the probability of an upward move in the price, conditional on the state of the order book. We study the diffusion limit of the price process and express the volatility of price changes in terms of parameters describing the arrival rates of buy and sell orders and cancelations. These analytical results provide some insight into the relation between order flow and price dynamics in order-driven markets.

Introduction

An increasing number of stocks are traded in electronic, order-driven markets, in which orders to buy and sell are centralized in a limit order book available to market participants 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 execution [START_REF] Alfonsi | Optimal execution strategies in limit order books with general shape function[END_REF], [START_REF] Predoiu | Optimal execution of a general one-sided limit-order book[END_REF])-but also from a fundamental perspective, since they provide a rare glimpse into the dynamics of supply and demand and their role in price formation.

Equilibrium models of price formation in limit order markets [START_REF] Parlour | Price dynamics in limit order markets[END_REF], [START_REF] Rosu | A dynamic model of the limit order book[END_REF]) have shown that the evolution of the price in such markets is rather complex and depends on the state of the order book. On the other hand, empirical studies on limit order books [START_REF] Bouchaud | How markets slowly digest changes in supply and demand[END_REF], [START_REF] Farmer | What really causes large price changes?[END_REF], [START_REF] Gourieroux | Intra-day market activity[END_REF], [START_REF] Hollifield | Empirical analysis of limit order markets[END_REF], [START_REF] Smith | Statistical theory of the continuous double auction[END_REF]) provide an extensive list of statistical features of order book dynamics that are challenging to incorporate in a single model. While most of these studies have focused on unconditional/steady-state distributions of various features of the order book, empirical studies [START_REF] Harris | The information content of the limit order book: evidence from NYSE specialist trading decisions[END_REF], Cont et al. (2010a)) show that the state of the order book contains information on short-term price movements so it is of interest to provide forecasts of various quantities conditional on the state of the order book. Providing analytically tractable models which enable to compute and/or reproduce conditional quantities which are relevant for trading and intraday risk management has proven to be challenging, given the complex relation between order book dynamics and price behavior.

The search for tractable models of limit order markets has led to the development of stochastic models which aim to retain the main statistical features of limit order books while remaining computationally manageable. Stochastic models also serve to illustrate how far one can go in reproducing the dynamic properties of a limit order book without resorting to detailed behavioral assumptions about market participants or introducing unobservable parameters describing agent preferences, as in more detailed market microstructure models.

Starting from a description of order arrivals and cancelations as point processes, the dynamics of a limit order book is naturally described in the language of queueing theory. [START_REF] Engle | Trades and Quotes: A Bivariate Point Process[END_REF] formulate a bivariate point process to jointly analyze trade and quote arrivals. Cont et al. (2010b) model the dynamics of a limit order book as a tractable multiclass queueing system and compute various transition probabilities of the price conditional on the state of the order book, using Laplace transform methods.

Summary

We propose a Markovian model of a limit order market, which captures some salient features of the dynamics of market orders and limit orders and their influence on price dynamics, yet is even simpler than the model of Cont et al. (2010b) and enables a wide range of properties of the price process to be computed analytically.

Our approach is motivated by the observation that, if one is primarily interested in the dynamics of the price, it is sufficient to focus on the dynamics of the (best) bid and ask queues. Indeed, empirical evidence shows that most of the order flow is directed at the best bid and ask prices [START_REF] Biais | An empirical analysis of the order flow and order book in the Paris Bourse[END_REF]) and the imbalance between the order flow at the bid and at the ask appears to be the main driver of price changes (Cont et al. (2010a)).

Motivated by this remark, we propose a parsimonious model in which the limit order book is represented by the number of limit orders (q b t , q a t ) sitting at the bid and the ask, represented as a system of two interacting queues. The remaining levels of the order book are treated as a 'reservoir' of limit orders represented by the distribution of the size of the queues at the 'next-to-best' price levels.

Through its analytical tractability, the Markovian version of our model allows to obtain analytical expressions for various quantities of interest such as the distribution of the duration until the next price change, the distribution and autocorrelation of price changes, and the probability of an upward move in the price, conditional on the state of the order book. Compared with econometric models of high frequency data [START_REF] Engle | Autoregressive conditional duration: a new model for irregularly-spaced transaction data[END_REF], [START_REF] Engle | Trades and Quotes: A Bivariate Point Process[END_REF] where the link between durations and price changes is specified exogenously, our model links these quantities in an endogenous manner, and provides a first step towards joint 'structural' modeling of high frequency dynamics of prices and order flow.

A second important observation is that order arrivals and cancelations are very frequent and occur at millisecond time scale, whereas, in many applications such as order execution, the metric of success is the volume-weighted average price (VWAP) so one is interested in the dynamics of order flow over a large time scale, typically tens of seconds or minutes. As shown in Table 1.1, thousands of order book events may occur over such time scales. This observation enables us to use asymptotic methods to study the link between price volatility and order flow in this model by studying the diffusion limit of the price process. In particular, we prove a functional central limit theorem for the price process and express the volatility of price changes in terms of parameters describing the arrival rates of buy and sell orders and cancelations. For example, we show (Theorem 1) that the variance of intraday price changes in a 'balanced' limit order market is given by the following simple relation:

σ 2 = πδ 2 λ D(f )
where δ is the 'tick size', λ is the intensity of order arrivals and D(f ) is a measure of market depth. These analytical results provide insights into the relation between order flow and price dynamics in order-driven markets. Comparison of these results with empirical data shows the main insights of the model to be correct.

Outline

The chapter is organized as follows. Section 2 introduces a reduced-form representation of a limit order book and presents a Markovian model in which limit orders, market orders and cancelations occur according to Poisson processes. Section 3 presents various analytical results for this model: we compute the distribution of the duration until the next price change (section 3.1), the probability of upward move in the price (section 3.2) and the dynamics of the price (section 3.4). In Section 4, we show that the price exhibits diffusive behavior at longer time scales and express the variance of price changes in terms of the parameters describing the order flow, thus establishing a link between volatility and order flow statistics.

2. A Markov model of limit order book dynamics 2.1. A stylized representation of a limit order book Empirical studies of limit order markets suggest that the major component of the order flow occurs at the (best) bid and ask price levels (see e.g. [START_REF] Biais | An empirical analysis of the order flow and order book in the Paris Bourse[END_REF]). All electronic trading venues also allow to place limit orders pegged to the best available price (National Best Bid Offer, or NBBO); market makers used these pegged orders to liquidate their inventories. Furthermore, studies on the price impact of order book events show that the net effect of orders on the bid and ask queue sizes is the main factor driving price variations (Cont et al. (2010a)). These observations, together with the fact that queue sizes at the best bid and ask of the consolidated order book are more easily obtainable (from records on trades and quotes) than information on deeper levels of the order book, motivate a reduced-form modeling approach in which we represent the state of the limit order book by • the bid price s b t and the ask price s a t

• the size of the bid queue q b t representing the outstanding limit buy orders at the bid, and • the size of the ask queue q a t representing the outstanding limit sell orders at the ask Figure 1 summarizes this representation.

If the stock is traded in several venues, the quantities q b and q a represent the best bids and offers in the consolidated order book, obtained by aggregating over all (visible) trading venues. At every time t, q b t (resp. q a t ) corresponds to all visible orders available at the bid price s b t (resp. s a t ) across all exchanges.

q a q b δ Quantities s b Price s a Figure 1
Simplified representation of a limit order book.

The state of the order book is modified by order book events: limit orders (at the bid or ask), market orders and cancelations (see Cont et al. (2010b,a), [START_REF] Smith | Statistical theory of the continuous double auction[END_REF]). A limit buy (resp. sell) order of size x increases the size of the bid (resp. ask) queue by x, while a market buy (resp. sell) order decreases the corresponding queue size by x. Cancelation of x orders in a given queue reduces the queue size by x. Given that we are interested in the queue sizes at the best bid/ask levels, market orders and cancelations have the same effect on the queue sizes (q b t , q a t ). The bid and ask prices are multiples of the tick size δ. When either the bid or ask queue is depleted by market orders and cancelations, the price moves up or down to the next level of the order book. The price processes s b t , s a t is thus a piecewise constant process whose transitions correspond to hitting times of the axes {(0, y), y ∈ N} ∪ {(x, 0), x ∈ N} by the process q t = (q b t , q a t ). If the order book contains no 'gaps' (empty levels), these price increments are equal to one tick:

• when the bid queue is depleted, the (bid) price decreases by one tick.

• when the ask queue is depleted, the (ask) price increases by one tick. If there are gaps in the order book, this results in 'jumps' (i.e. variations of more than one tick) in the price dynamics. We will ignore this feature in what follows but it is not hard to generalize our results to include it.

The quantity s a t -s b t is the bid-ask spread, which may be one or several ticks. As shown in Table 2, for liquid stocks the bid-ask spread is equal to one tick for more than 98% of observations. When either the bid or ask queue is depleted, the bid-ask spread widens immediately to more than one tick. Once the spread has increased, a flow of limit sell (resp. buy) orders quickly fills the gap and the spread reduces again to one tick. When a limit order is placed inside the spread, all the limit orders pegged to the NBBO price move in less than a millisecond to the price level corresponding to this new order. Once this happens, both the bid price and the ask price have increased (resp. decreased) by one tick.

Bid

The histograms in Figure 2 show that this 'closing' of the spread takes place very quickly: as shown in Figure 2 (left) the lifetime of a spread larger than one tick is of the order of a couple of milliseconds, which is negligible compared to the lifetime of a spread equal to one tick (Figure 2, right). In our model we assume that the second step occurs infinitely fast: once the bid-ask spread widens, it returns immediately to one tick. For the example of Dow Jones stocks (Figure 2), this is a reasonable assumption since the widening of the spread lasts only a few milliseconds. This simply means that we are not trying to describe/model how the orders flow inside the bid-ask spread at the millisecond time scale and, when we describe the state of the order book after a price change we have in mind the state of the order book after the bid-ask spread has returned to one tick. Under this assumption, each time one of the queues is depleted, both the bid queue and the ask queues move to a new position and the bid-ask spread remains equal to one tick after the price change. Thus, under our assumptions the bid-ask spread is equal to one tick, i.e. s a t = s b t + δ, resulting in a further reduction of dimension in the model.

Once either the bid or the ask queue are depleted, the bid and ask queues assume new values. Instead of keeping track of arrival, cancelation and execution of orders at all price levels (as in Cont et al. (2010b), [START_REF] Smith | Statistical theory of the continuous double auction[END_REF]), we treat the queue sizes after a price change as stationary variables drawn from a certain distribution f on N 2 which represents, in a statistical sense, the depth of the order book after a price change: f (x, y) represents the probability of observing (q b t , q a t ) = (x, y) right after a price increase. Similarly, we denote f (x, y) the probability of observing (q b t , q a t ) = (x, y) right after a price decrease. More precisely, denoting by F t the history of prices and order book events on [0, t],

• if q a t-= 0 then (q b t , q a t ) is a random variable with distribution f , independent from F t-. • if q b t-= 0 then (q b t , q a t ) is a random variable with distribution f , independent from F t-. The distributions f and f summarize the interaction of the queues at the best bid/ask levels with the rest of the order book, viewed here as a 'reservoir' of limit orders. Figure 3 shows the (joint) empirical distribution of bid and ask queue sizes after a price move for Citigroup stock on June 26th 2008.

In summary, state of the limit order book is thus described by a continuous-time process X t = (s b t , q b t , q a t ) which takes values in the discrete state space δZ × N 2 , with piecewise constant sample paths whose transitions correspond to the order book events. Denoting by (T a i , i ≥ 1) (resp. T b i ) the durations between two consecutive orders arriving at the ask (resp. the bid) and V a i (resp. V b i ) the size of the associated change in queue size, the above assumptions translate into the following dynamics for X t = (s b t , q b t , q a t ): • If an order or cancelation arrives on the ask side at time T :

(s b T , q b T , q a T ) = (s b T -, q b T -, q a T -+ V a i )1 q a T ->-V a i + (s b T -+ δ, R b i , R a i )1 q a T -≤-V a i ,
• If an order or cancelation arrives on the bid side i.e.

T ∈ {T b i , i ≥ 1}: (s b T , q b T , q a T ) = (s b T -, q b T -+ V b i , q a T -)1 q b T ->-V b i + (s b T --δ, Rb i , Ra i )1 q b T -≤-V b i , and (R i ) i≥1 = (R b i , R a i ) i≥1 is a sequence of IID variables with (joint) distribution f , and ( Ri ) i≥1 = ( Rb i , Ra i ) i≥1 is a sequence of IID variables with (joint) distribution f .

A Markov model for order book dynamics

To give a complete statistical description of the dynamics of the limit order book, we need to describe the distributional properties of the sequences T a i , T b i , V a i , V b i describing the timing and size of order book events.

We assume that these events occur according to independent Poisson processes:

• Market buy (resp. sell) orders arrive at independent, exponential times with rate µ,

• Limit buy (resp. sell) orders at the (best) bid (resp. ask) arrive at independent, exponential times with rate λ,

• Cancelations occur at independent, exponential times with rate θ.

• These events are mutually independent.

• All orders sizes are equal (assumed to be 1 without loss of generality). Denoting by (T a i , i ≥ 1) (resp. T b i ) the durations between two consecutive queue changes at the ask (resp. the bid) and V a i (resp. V b i ) the size of the associated change in queue size, the above assumptions translate into the following properties for the sequences

T a i , T b i , V a i , V b i : (i) (T a i ) i≥0
is a sequence of independent random variables with exponential distribution with parameter λ + θ + µ, (ii) (T b i ) i≥0 is a sequence of independent random variables with exponential distribution with parameter λ + θ + µ, (iii) (V a i ) i≥0 is a sequence of independent random variables with

P[V a i = 1] = λ λ + µ + θ and P[V a i = -1] = µ + θ λ + µ + θ , (1)
(iv) (V b i ) i≥0 is a sequence of independent random variables with

P[V b i = 1] = λ λ + µ + θ and P[V b i = -1] = µ + θ λ + µ + θ (2)
• These sequences are independent. Under these assumptions q t = (q b t , q a t ) is thus a Markov process, taking values in N 2 , whose transitions correspond to the order book events {T a i , i ≥ 1} ∪ {T b i , i ≥ 1}: • At the arrival of a new limit buy (resp. sell) order the bid (resp. ask) queue increases by one unit. This occurs at rate λ.

• At each cancelation or market order, which occurs at rate θ + µ, either: (a) the corresponding queue decreases by one unit if it is > 1, or (b) if the ask queue is depleted then q t is a random variable with distribution f . (c) if the bid queue is depleted then q t is a random variable with distribution f . The values of λ and µ + θ are readily estimated from high-frequency records of trades and quotes (see Cont et al. (2010b) for a description of the estimation procedure). Table 3 gives examples of such parameter estimates for the stocks mentioned above. We note that in all cases λ < µ + θ but that the difference is small: |(µ + θ) -λ| λ.

Remark 1 (Independence assumptions). The IID assumption for the sequences (R n ), ( Rn ) is only used in Section 4. The results of Section 3 do not depend on this assumption. Estimates for the intensity of limit orders and market orders+cancellations, in number of batches per second (each batch representing 100 shares) on June 26th, 2008).

Quantities of interest

In applications, one is interested in computing various quantities that intervene in high frequency trading such as:

• the conditional distribution of the duration between price moves, given the state of the order book (Section 3.1),

• the probability of a price increase, given the state of the order book (Section 3.2),

• the dynamics of the price autocorrelations and distribution and autocorrelations of price changes (section 3.4), and

• the volatility of the price (section 4). We will show that all these quantities may be characterized analytically in this model, in terms of order flow statistics.

Analytical results

The high-frequency dynamics of the price may be described in terms of durations between successive price changes and the magnitude of these price changes. It is of interest to examine what information the current state of the (consolidated) order book gives about the dynamics of the price. We now proceed to show how the model presented above may be used to compute the conditional distributions of durations and price changes, given the current state of the order book, in terms of the arrival rates of market orders, limit orders and cancelations. The results of this section do not depend on the assumptions on the sequences (R n ), ( Rn ).

Duration until the next price change

We consider first the distribution of the duration until the next price change, starting from a given configuration (x, y) of the order book. We define • σ a the first time when the ask queue (q a t , t ≥ 0) is depleted, • σ b the first time when the bid queue (q b t , t ≥ 0) is depleted Since the queue sizes are constant between events, one can express these stopping times as:

σ a = inf{T a 1 + ... + T a i , q a T a 1 +...+T a i -+ V a i = 0} σ b = inf{T b 1 + ... + T b i , q b T b 1 +...+T b i -+ V b i = 0}
The price (s t , t ≥ 0) moves when the queue q t = (q b t , q a t ) hits one of the axes: the duration until the next price move is thus τ = σ a ∧ σ b . The following theorem gives the distribution of the duration τ , conditional on the initial queue sizes:

Proposition 1 (Distribution of duration until next price move). The distribution of τ conditioned on the state of the order book is given by:

P[τ > t|q b 0 = x, q a 0 = y] = ( µ + θ λ ) x+y ψ x,λ,θ+µ (t)ψ y,λ,θ+µ (t) (3)
where

ψ n,λ,θ+µ (t) = ∞ t n u I n (2 λ(θ + µ)u)e -u(λ+θ+µ) du (4)
and I n is the modified Bessel function of the first kind. The conditional law of τ has a regularly varying tail

• with tail exponent 2 if λ < µ + θ • with tail exponent 1 if λ = µ + θ. In particular, if λ = µ + θ, E[τ |q b 0 = x, q a 0 = y] = ∞ whenever x > 0, y > 0.
Proof. Since (q a t , t ≥ 0) follows a birth and death process with birth rate λ and death rate µ + θ, L(s, x) := E[e -sσa |q a 0 = x] satisfies:

L(s, x) = λL(s, x + 1) + (µ + θ)L(s, x -1) λ + µ + θ + s .
We can find the roots of the polynomial:

λX 2 -(λ + µ + θ + s)X + µ + θ; one root is > 1, the other is < 1; since L(s, 0) = 1 and lim x→∞ L(s, x) = 0, L(s, x) = ( (λ + µ + θ + s) -((λ + µ + θ + s)) 2 -4λ(µ + θ) 2λ ) x .
Moreover if we use the relation

P[τ > t|q b 0 = x, q a 0 = y] = P[σ b > t|q b 0 = x]P[σ a > t|q a 0 = y], P[τ > t|q b 0 = x, q a 0 = y] = ∞ t L(u, x)du ∞ t L(u, y)du.
This Laplace transform may be inverted (see [START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF], XIV.7)) and the inversion yields λ+θ+µ) , which gives us the expected result.

L(t, x) = x t ( µ + θ λ ) x I x (2 λ(θ + µ)t)e -t(
The above result allows in particular to study the tail behavior of the conditional distribution of τ :

• If λ < µ + θ:

L(s, x) = α(s) x ∼ s→0 1 - x(λ + µ + θ) 2λ(µ + θ -λ) s,
so Karamata's Tauberian theorem (Feller 1971, XIII.5) yields

P[σ a > t|q a 0 = x] ∼ t→∞ x(λ + µ + θ) 2λ(µ + θ -λ) 1 t ;
therefore the conditional law of the duration τ is a regularly varying with tail index 2

P[τ > t|q b 0 = x, q a 0 = y] ∼ t→∞ xy(λ + µ + θ) 2 λ 2 (µ + θ -λ) 2 1 4t 2 .
(5)

• If the order flow is balanced i.e. λ = µ + θ then

L(s, x) = α(s) x ∼ s→0 1 - x √ λ √ s,
the law of σ a is regularly-varying with tail index 1/2 and

P[σ a > t|q a 0 = x] ∼ t→∞ x √ πλ 1 √ t .
The duration then follows a heavy-tailed distribution with infinite first moment: The expression given in (3) is easily computed by discretizing the integral in (4). Plotting (3) for a fine grid of values of t typically takes less than a second on a laptop. Figure 3 gives a numerical example, with λ = 12 sec -1 , µ + θ = 13 sec -1 , q a 0 = 4, q b 0 = 5 (queue sizes are given in multiples of average batch size).

P[τ > t|q b 0 = x, q a 0 = y] ∼ t→∞ xy πλ 1 t ; (6)

Probability of price increase: balanced limit order book

Starting from a given configuration of the limit order book, the probability that the next price move is an increase is given by the probability that the process (q b t , q a t ) hits the x-axis before the y-axis When λ = µ + θ, i.e. when the flow of limit orders is balanced by the flow of market orders and cancelations, this probability can be computed analytically in terms of hitting time distributions of a random walk in the orthant: Proposition 2. For (n, p) ∈ N 2 , the probability p up 1 (n, p) that the next price move is an increase, conditioned on having the n orders on the bid side and p orders on the ask side is:

p up 1 (n, p) = 1 π π 0 (2 -cos(t) -(2 -cos(t)) 2 -1) p sin(nt) cos( t 2 ) sin( t 2 )
dt.

(7)

Proof. One can notice that ∀t ≤ τ, q t = M N 2λt where (M n , n ≥ 0) is a symmetric random walk in the positive orthant Z 2 + killed at the boundary and (N 2λt , t ≥ 0) is a Poisson process with parameter 2λ. Hence the probability of an upward move in the price starting from a configuration q b t = n, q a t = p for the order book is equal to the probability that the random walk M starting from (n, p) hits the x-axis before the y-axis. The generator of the bivariate random walk (M n , n ≥ 1) is the discrete Laplacian so

p up 1 (n, p) = P[σ a < σ b |q b 0-= n, q a 0-= p] satisfies, for all n ≥ 1 and p ≥ 1, 4p up 1 (n, p) = p up 1 (n + 1, p) + p up 1 (n -1, p) + p up 1 (n, p + 1) + p up 1 (n, p -1), (8)
with the boundary conditions: p up 1 (0, p) = 0 for all p ≥ 1 and p up 1 (n, 0) = 1 for all n ≥ 1. This problem is known as the discrete Dirichlet problem; solutions of ( 8) are called discrete harmonic functions. (Lawler and Limic 2010, Ch. 8) show that for all t ≥ 0, the functions f t (x, y) = e xr(t) sin(yt), and ft (x, y) = e -xr(t) sin(yt) with r(t) = cosh -1 (2 -cos t)

are solutions of (8). In (Lawler and Limic 2010, Corollary 8.1.8) it is shown that the probability that a simple random walk (M k , k ≥ 1) starting at (n, p) ∈ Z + × Z + reaches the axes at (x, 0) is 2 π π 0 e -r(t)p sin(nt) sin(tx)dt, therefore

p up 1 (n, p) = ∞ k=1 2 π π 0 e -r(t)p sin(tn) sin(tk)dt. Since m k=1 sin(kt) = sin( mt 2 ) sin( (m+1)t 2 ) sin(t/2) = cos( t 2 ) -cos((m + 1 2 )t) 2 sin(t/2) ,
using integration by parts we see that the second term leads to the integral:

π 0 e -r(t)p sin(nt) sin(t/2)

g(t) cos((m + 1/2)t)dt = - 1 m + 1 2 π 0 g (t) sin((m + 1 2 )t)dt → m→∞ 0.
since g is bounded. So finally:

p up 1 (n, p) = 1 π π 0 e -r(t)p sin(tn) cos( t 2 ) sin( t 2 )
dt.

Noting that e -r(t) = (2 -cos(t) -(2 -cos(t)) 2 -1) we obtain the result. Note that the conditional probabilities (7) are, in the case of a balanced order book, independent of the parameters describing the order flow.

The expression ( 7) is easily computed numerically: Figure 3.2 displays the shape of the function p up 1 . Comparison with empirical data for CitiGroup stock (June 2008) shows good agreement between the theoretical value (7) and the empirical transition frequencies of the price conditional on the state of the consolidated order book.

Probability of price increase: asymmetric order flow

In this subsection we relax the symmetry assumptions above and allow the intensity of limit and market orders at the bid and the ask to be different; more precisely we assume that:

• Limit orders at the ask arrive at independent, exponential time with parameter λ a • Market orders and cancelations at the ask arrive at independent, exponential time with parameter µ a + θ a

• Limit orders at the bid arrive at independent, exponential time with parameter λ b • Market orders and cancelations at the bid arrive at independent, exponential time with parameter µ b + θ b and that these Poisson processes are independent. The dynamics of bid and ask queues may be then represented as

q t = M N Λt for Λ = λ a + µ a + θ a + λ b + µ b + θ b ,
where N Λt is a Poisson process with intensity Λ and (M n , n ≥ 0) is a random walk on Z 2 killed when it hits either the x-axis or the y-axis whose the transition probabilities are:

p 0,1 = λ a Λ p 1,0 = λ b Λ p 0,-1 = µ a + θ a Λ p -1,0 = µ b + θ b Λ . (9)
The following result generalizes Proposition 2 for an asymmetric order flow.

Proposition 3. Given (q b , q a ) = (n, p), the probability p up 1 (n, p) that the next price move in an increase in:

p up 1 (n, p) = 1 - 1 π µ a + θ a λ a p 2[λ a (µ a + θ a )] 1/2 µ a + θ a + λ a π 0 Z(t) n sin(pt) sin(t)× × 2λ b Z t -(Λ -2[λ a (µ a + θ a )] 1/2 cos(t)) Λ 2[(µ a +θ a )λ a ] 1/2 λ b +µ b +θ b cos(t) -Λ Λdt (Λ -2[(µ a + θ a )λ a ] 1/2 cos(t)) 2 -4(µ b + θ b )λ b ,
where (Z t , t ≥ 0) is the function defined by:

∀t ≥ 0, Z t = Λ -2[(µ a + θ a )λ a ] 1/2 cos(t) -(Λ -2[(µ a + θ a )λ a ] 1/2 cos(t)) 2 -4λ b (θ b + µ b ) 2λ a .
Proof. Using results from [START_REF] Kurkova | Random walk in Z 2 + with non-zero drift absorbed at the axes[END_REF], it is shown in [START_REF] Raschel | Random walks in the quarter plane: the probability to hit the x-axis before the y-axis[END_REF] that the probability that M starting from (n, p) hits the x-axis before the y-axis is given by

1 - 1 π p 0,-1 p 0,1 p 2[p 0,1 p 0,-1 ] 1/2 p 0,1 + p 0,-1 π 0 Z(t) n sin(pt) sin(t)× × 2p 1,0 Z t -(1 -2[p 0,1 p 0,-1 ] 1/2 cos(t)) 2[p 0,1 p 0,-1 ] 1/2 p 0,1 +p 0,-1 cos(t) -1 dt (1 -2[p 0,1 p 0,-1 ] 1/2 cos(t)) 2 -4p 1,0 p -1,0 ,
where (Z t , t ≥ 0) is the function defined by:

∀t ≥ 0, Z t = 1 -2[p 0,1 p 0,-1 ] 1/2 cos(t) -(1 -2[p 0,1 p 0,-1 ] 1/2 cos(t)) 2 -4p 1,0 p -1,0 2p 1,0 .
Applying this formula to the random walk whose transition probabilities are given in Equation ( 9) we obtain that the probability of a price increase at the next price move, starting from a configuration q b t = n, q a t = p of the queues is equal to the probability that the random walk M starting from (n, p) hits the x-axis before the y-axis; thus, it is given by equation (3).

Dynamics of the price

The high-frequency dynamics of the price is described by a piecewise constant, right continuous process (s t , t ≥ 0) whose jumps times correspond to times when the order book process (q t , t ≥ 0) hits one of the axes. Denote by (τ 1 , τ 2 , ...) the successive durations between price changes. The number of price changes that occur during [0, t] is given by

N t := max{ n ≥ 0, τ 1 + ... + τ n ≤ t } At t = τ i , s τ i = s τ i-+ 1 if q τ a i-= 0 and s τ i = s τ i--1 if q τ b i- = 0. (X 1 , X 2 , X 3 , .
.., X n , ...) are the successive moves in the price which, in general, are not independent random variables. We have

s t = Z N t where Z n = n i=1 X i
is the price after n price changes have occurred. Hence, for all t ≥ 0, s t = Z N t . We are interested in the n-step ahead distribution of the price change:

p up n (x, y) = P[X n = +δ | (q b 0 , q a 0 ) = (x, y)] (10)
For n = 1 this corresponds to the probability p up 1 (x, y) = p 1 (x, y) of an upward price move, computed in Theorem 2. To simplify the analysis we use, in this Section and the next one, the following symmetry assumption:

Assumption 1 (Bid-ask symmetry). f (x, y) = f (y, x).

A key quantity for studying the dynamics of the price is the probability of two successive price changes in the same direction,

p cont = P[X k+1 = δ|X k = δ] = P[X k+1 = -δ|X k = -δ] ( 11
)
which may be expressed in terms of the distribution of queue sizes f after a price change:

Proposition 4. Let p cont = P[X 2 = δ|X 1 = δ] = P[X 2 = -δ|X 1 = -δ]
be the probability of two successive price moves in the same direction.

• p cont = ∞ i=1 ∞ j=1 f (i, j)p up 1 (i, j). • ∀k ≥ 1, E[X k |q b 0 = x, q a 0 = y] = (2p up 1 (x, y) -1)(2p cont -1) k-1 . • cov(X 1 , X 2 |q b 0 = x, q a 0 = y) = δ 2 (2p cont -1)(1 -(2p up 1 (x, y) -1) 2 ).
• Conditional on the current state of the limit order book, the n-step ahead distribution of the price change is given by:

p up n (x, y) := P[X n = δ|q b 0 = x, q a 0 = y] = 1 + (2p cont -1) n-1 (2p up 1 (x, y) -1) 2 . ( 12
)
where f is the joint distribution of queue sizes after a price increase. This condition is satisfied for most high-frequency data sets of Dow Jones stocks we have examined. For example, for CitiGroup stock we find

p cont = ∞ i=1 ∞ j=1 f (i, j)p up 1 (i, j) 0.35
This asymmetry of f corresponds to the fact that, after an upward price move, the new bid queue is generally smaller than the ask queue since the ask queue corresponds to the limit order previously sitting at second best ask level, while the bid queue results from the accumulation of orders over the very short period since the last price move. Under this condition, high frequency increments of the price are negatively correlated: an increase in the price is more likely to be followed by a decrease in the price.

Remark 3. The sequence of price increments (X 1 , X 2 , ...) is uncorrelated if and only if

p cont = i,j≥1
f (i, j)p up 1 (i, j) = 1/2 3.5. "Efficient" price Various authors (see e.g. Robert and Rosenbaum ( 2011)) have considered models in which the evolution of transaction prices is based on a non-observed (semi)martingale ŝ, sometimes called the "efficient" price: the observed price is then either a noisy version of ŝ or the value of ŝ rounded to the nearest tick. Given the probability p up 1 (q b t , q a t ) that the next price move is an "uptick" (Equation ( 7)), we can construct an auxiliary process ŝ whose value ŝt represents the expected value of the price after its next move: ∀t ≥ 0, ŝt = (s t + δ) p up 1 (q b t , q a t ) + (s t -δ) 1 -p up 1 (q b t , q a t ) , ŝt = s t + δ(2p up 1 (q b t , q a t ) -1), (ŝ t , t ≥ 0) is a continuous-time stochastic process with values between s t -δ and s t + δ:

∀t ≥ 0, s t -δ ≤ ŝt ≤ s t + δ.
The process ŝ incorporates the information on the price s t and the state of the order book (q b , q a ) insofar as it affects the next price move. The following result shows

Proposition 5. If p cont = 1/2 then (ŝ t , t ≥ 0) is a martingale. Remark 4. This condition is verified in particular if ∀(i, j) ∈ N 2 , f (i, j) = f (j, i) but more gener- ally if i,j≥1 p up 1 (i, j)f (i, j) = 1 2
Proof. Let (τ 1 , τ 2 , ..., τ k ) the sequence of times when the price s moves and (X 1 , ..., X n ) the sequence of consecutive price moves. Since p cont = 1/2, (X 1 , ..., X k , ...) is a sequence of I.I.D bernoulli random variables with parameter 1/2. Therefore we have the following property:

∀(i, j) ∈ N 2 , i < j, E[s τ j |F τ i ] = s τ i .
The function p up 1 , from equation ( 7), satisfies the equation Lp up 1 = 0, where L is the generator of the process (q b , q a ). Hence p up 1 is an harmonic function for the process (q b , q a ), and the process (p up 1 (q b t , q a t ), t ≥ 0) is a martingale. We proved that

∀ s ≤ t < τ 1 , E[ŝ t |F s ] = ŝs .
By recurrence on k, one can easily notice that

∀ k ≥ 1, ∀τ k ≤ t < τ k+1 , E[ŝ t |F τ 1 ] = ŝτ 1 . Assuming s ≤ τ 1 ≤ t, E[ŝ t |F s ] = E[ŝ t |F s , X 1 = 1]P[X 1 = 1|F s ] + E[ŝ t |F s , X 1 = -1]P[X 1 = -1|F s ], = E[ŝ t |F s , X 1 = 1]p up 1 (q b t , q a t ) + E[ŝ t |F s , X 1 = -1](1 -p up 1 (q b t , q a t )), = (s s + δ)p up 1 (q b s , q a s ) + s s (1 -p up 1 (q b s , q a s )), = s s + δ(2p up
1 (q b s , q a s ) -1) = ŝs , which completes the proof.

Contrarily to the 'latent price' models alluded to above, here ŝ is a function of the state variables (s t , q b t , q a t ) and thus is observable, provided one observes trades and quotes. Remark 5. When p cont = 1/2, the process (ŝ t , t ≥ 0) fails to possess the martingale property. When p cont < 1/2, the jump of ŝ is negative after a price increase and positive after a price decrease.

Diffusion limit of the price process

As discussed in Section 3.4, the high frequency dynamics of the price is described by a piecewise constant stochastic process s t = Z N t where

Z n = X 1 + ... + X n and N t = sup{k; τ 1 + ... + τ k ≤ t} is the number of price changes during [0, t].
However, over time scales much larger than the interval between individual order book events, prices are observed to have diffusive dynamics and modeled as such. To establish the link between the high frequency dynamics and the diffusive behavior at longer time scales, we shall consider a time scale t n over which the average number of order book events is of order n and exhibit conditions under which the scaled price process (s n t :=

s tn √ n , t ≥ 0) n≥1
satisfies a functional central limit theorem i.e. converges in distribution to a non-degenerate process (p t , t ≥ 0) as n → ∞. The choice of the time scale t n is such that

τ 1 + ... + τ n t n
has a well-defined limit: it is imposed by the distributional properties of the durations which, as observed in Section 3.1, are heavy tailed. In this section, we show that, under a symmetry condition, this limit can be identified as a diffusion process whose diffusion coefficient may be computed from the statistics of the order flow driving the limit order book. Assume λ + θ ≤ µ and that the joint distribution f of the queue sizes after a price move satisfies:

D(f ) = ∞ i=1 ∞ j=1 ijf (i, j) < ∞ (14)
The quantity D(f ) represents a measure of market depth: more precisely, D(f ) is the geometric average of the size of the bid queue and the size of the ask queue after a price change.

In this section we assume that f is symmetric: ∀i, j ≥ 0, f (i, j) = f (j, i). Under this assumption, price increments form a sequence (X i , i ≥ 0) of independent random variables with distribution:

P[X 1 = δ] = P[X 1 = -δ] = 1 2 .
We will show that the limit p is then a diffusion process which describes the dynamics of the price at lower frequencies. In particular, we will compute the volatility of this diffusion limit and relate it to the properties of the order flow.

In the following D denotes the space of right continuous paths ω : [0, ∞) → R 2 with left limits, equipped with the Skorokhod topology J 1 , and ⇒ will designate weak convergence on (D, J 1 ) (see [START_REF] Billingsley | Convergence of Probability Measures[END_REF], [START_REF] Whitt | Stochastic Process Limits[END_REF] for a discussion).

Balanced order book

We first consider the case of a balanced order flow for which the intensity of market orders and cancelations is equal to the intensity of limit orders. The study of high-frequency quote data indicates that this is an empirically relevant case for many liquid stocks: as shown in Table 3, the imbalance between arrival of limit orders on one hand and market orders/ cancelations on the other hand is around 5% or less for these stocks.

For balanced order flow, we proved in Section 3.1 that the distribution of price duration τ conditioned on observing i shares at the bid and j shares at the ask at t = 0 has a tail index 1:

P[τ > t|q b 0 = i, q a 0 = j] ∼ ij πλt ; (15)
The unconditional distribution of price durations keeps a tail index of one:

P[τ > t] t→∞ ∼ ∞ i=1 ∞ j=1 ijf (i, j) πλt = D(f ) πλt . ( 16
)
The sequence of durations between consecutive move of the price consecutive (τ 1 , τ 2 , τ 3 , ...) is a sequence of IID random variables with tail index 1. The following lemma 1 gives a central limit theorem for this sequence of durations.

Lemma 1. The sequence of durations (τ 1 , τ 2 , τ 3 , ...) satisfies

τ 1 + τ 2 + ... + τ k k log k k→∞ ⇒ D(f ) πλ .
Proof. The sequence (τ 1 , τ 2 , τ 3 , ...) is a sequence of i.i.d, regularly varying random variables, with tail index equal to 1. Let L(s) be the Laplace transform of the distribution of τ 2 :

L(s) = 1 -s ∞ 0 exp(-st)P[τ > t]dt. Since 1 log(n) n 0 exp(-st/n log(n))P[τ > t]dt → D(f )
πλ using an integration by parts we obtain

1 log(n) ∞ n exp(-st/n log(n))P[τ > t]dt → 0. Therefore s n log(n) ∞ 0 exp(- st n log n )P[τ > t]dt = D(f ) πλ s n + o( 1
n ) which implies (L(s/p up 1 (n))) n = (1 + D(f ) πλ + o( 1 n )) n → exp (- D(f ) πλ ).
So one can conclude that

τ 1 + τ 2 + ... + τ n n log n ⇒ D D(f ) πλ .
Using the rescaling t → t n = tn log n, the process N tn counting the number of price change in the interval [0, n log n] satisfies the implicit equation:

t n N tn log N tn n→∞ ∼ D(f ) πλ . ( 17
)
Hence if one defines ρ : (1, ∞) → (1, ∞) the inverse function of t → t log t:

∀t > 1, ρ(t) log ρ(t) = t, (18)
and knowing that ρ(t) ∼ t→∞ t log(t) , one can notice that the asymptotic number of price change occurring during [0, t n ] is approximately:

N tn ∼ ρ πλt n D(f ) ∼ πλt n D(f ) log(πλt n /D(f )) as n → ∞, (19)
N tn ∼ πλtn log n D(f ) log(n log n) ∼ πλt D(f ) n as n → ∞. (20)
Equation ( 20) shows that the number of price move occurring during [0, tn log n] is an order of n for large n, and proportional to the ratio of the order intensity λ to the order book depth D(f ).

Since each price change is ±δ, this factor πλ D(f ) also shows up in the expression of price volatility:

Theorem 1. If λ = µ + θ, s tn log n √ n , t ≥ 0 n→∞ ⇒ δ πλ D(f ) W t , t ≥ 0
where δ is the tick size, D(f ) is given by (14) and W is a standard Brownian motion.

Proof. Let t n = tn log n. One can decompose the process ( s tn log n √ n , t ≥ 0) as

s tn log n √ n = Z([tnπλ/D(f )])δ √ n + Z(N tn )δ √ n - Z([tnπλ/D(f )])δ √ n (21)
Since (X 1 , X 2 , ...) is a sequence of IID random variables with mean zero, one can apply the Donsker's invariance principle to the sequence of processes

( Z([tπλn/D(f )])δ √ n , t ≥ 0) n≥1 ⇒ δ πλ D(f ) W t , t ≥ 0 as n → ∞, (22)
where the convergence is in distribution on the Skorokhod space. As shown in (20),

N tn log n n→∞ ∼ ntπλ D(f ) , therefore for t ≥ 0, Z(N tn )δ √ n - Z([tnπλ/D(f )])δ √ n n→∞ ⇒ 0. ( 23
)
Hence the finite dimensional distributions of the sequence of processes

Z(N tn )δ √ n - Z(tπλ/D(f ))δ √ n , t ≥ 0 n≥1 (24)
converge to a point mass at zero. Since this sequence of processes is tight on (D, J 1 ), this sequence of processes converges weakly to zero on (D, J 1 ) (see [START_REF] Whitt | Stochastic Process Limits[END_REF]). So finally

s tn log n √ n , t ≥ 0 n→∞ ⇒ δ πλ D(f ) W. (25)
4.2. Empirical test using high-frequency data Theorem 1 relates the 'coarse-grained' volatility of intraday returns at lower frequencies to the high-frequency arrival rates of orders. Denote by τ 0 = 1/λ the typical time scale separating order book events. Typically τ 0 is of the order of milliseconds. In plain terms, Theorem 1 states that, observed over a time scale τ 2 >> τ 0 (say, 10 minutes), the price has a diffusive behavior with a diffusion coefficient given by

σ n = δ nπλ D(f ) ( 26
)
where δ is the tick size, n is an integer verifying n ln n τ 0 = τ 2 which represents the average number of orders during an interval τ 2 and D(f ), the geometric mean of the size of the bid queue and the size of the ask queue after a price change, is a measure of market depth. Formula (26) links properties of the price to the properties of the order flow. the left hand side represents the variance of price changes, whereas the right hand side only involves the tick size and quantities: it yields an estimator for price volatility which may be computed without observing the price!

The relation ( 26) has an intuitive interpretation. It shows that, in two 'balanced' limit order markets with the same tick size and same rate of arrival of orders at the next bid/ask, the market with higher depth of the next-to-best queues will lead to lower price volatility.

More precisely, this formula shows that the microstructure of order flow affects price volatility through the ratio λ/D(f ) where λ is the rate of execution/cancelation of limit orders and D(f ), given by ( 14), is a measure of market depth: in fact, our model predicts a proportionality between the variance of price increments and this ratio. This is an empirically testable prediction.

Figure 4.2 compares, for stocks in the Dow Jones index, the standard deviation of 10-minute price increments with λ/D(f ). We observe that, indeed, stocks with a higher value of the ratio λ/D(f ) have a higher variance, and standard deviation of price increments increases roughly proportionally to λ/D(f ). Remark 6. When the intensity of all orders coming in the limit order book is multiplied by the same factor x,

• The intensity of limit orders becomes λx • The intensity of market orders and cancelations becomes (µ + θ)x • The depth of the limit order book increases by a factor x, so D(f ) increases by a factor x 2 . Substituting in the above formula, we then see that price volatility is decreased by a factor 1

x . Interestingly, [START_REF] Rosu | A dynamic model of the limit order book[END_REF] shows the same dependence in 1/ √ x of price volatility using an equilibrium approach.

Case when market orders and cancelations dominate

We now consider the case in which the flow of market orders and cancelations dominates that of limit orders: λ < θ + µ. In this case, price changes are more frequent since the order queues are depleted at a faster rate than they are replenished by limit orders. We also obtain a diffusion limit though with a different scaling: where ψ x,λ,µ+θ is given by (4). Then

s nt √ n , t ≥ 0 n→∞ ⇒ 1 m(λ, θ + µ, f ) δW t , t ≥ 0
where W is a standard Brownian motion.

Proof. The sequence (τ 2 , τ 3 , ...) is a sequence of i.i.d random variables with finite mean equal to m(λ, θ + µ, f ). We apply the law of large numbers:

τ 1 + τ 2 + ... + τ n n → m(λ, θ + µ, f ) as n → ∞.
Therefore, ∀t ≥ 0, N tn ∼[ tn m(λ, θ + µ, f ) ] as n → ∞.

The rest of the proof follows the lines of the proof of theorem 1. We start to decompose for all n ≥ 1 the process ( s[nt] √ n , t ≥ 0) in two terms:

s [nt] √ n = δZ [nt/m(λ
The second term of the decomposition converges to zero:

( δZ N nt √ n - δZ [nt/m(λ,θ+µ,f )] √ n , t ≥ 0) n→∞ ⇒ 0, (29)
which concludes the proof.

Variance of price change at intermediate frequency Theorem 2 leads to an expression of the variance of the price at a time scale τ τ 0 , where τ 0 (∼ ms) is the average interval between order book events:

σ 2 = τ τ 0 1 m(λ, θ + µ, f ) δ 2 (30)
Here, m(λ, θ + µ, f ) represents the expected hitting time of the axes by the Markovian queuing system q with parameters (λ, θ + µ) and random initial condition with distribution f . This equation relates the variance of price changes (over a time scale τ 2 ) to the tick size and the statistical properties of the order flow.

Conclusion

We have exhibited a simple model of a limit order market in which order book events are described in terms of a Markovian queueing system. The analytical tractability of our model allow us to compute various quantities of interest such as

• the distribution of the duration until the next price change,

• the distribution of price changes, and • the diffusion limit of the price process and its volatility. in terms of parameters describing the order flow. These results provide analytical insights into the relation between price dynamics and order flow, in particular the relation between liquidity and volatility, in a limit order market.

We view this stylized model as a first step in the analytical study of realistic stochastic models of order book dynamics. Yet, comparison with empirical data shows that even our simple modeling setup is capable of yielding useful analytical insights into the relation between volatility and order flow, worthy of being further pursued. Moreover, the connection with two-dimensional queueing systems allow us to use the rich analytical theory developed for these systems (see [START_REF] Cohen | Boundary value problems in queueing system analysis[END_REF]) to compute many other quantities. We hope to pursue further some of these ramifications in future work.

A relevant question is to examine which of the above results are robust to departures from the model assumptions and whether the intuitions conveyed by our model remain valid in a more general context where one or more of these assumptions are dropped. This issue is further studied in [START_REF] Cont | Order book dynamics in liquid markets: limit theorems and diffusion approximations[END_REF] where we explore a more general queueing model relaxing some of the assumptions above.

Figure 2

 2 Figure 2Left: Average lifetime, in milliseconds of a spread larger than one tick for Dow Jones stocks. Right: Average lifetime, in milliseconds of a spread equal to one tick.

Figure 3

 3 Figure 3Joint density of bid and ask queue sizes after a price move (Citigroup, June 26th 2008).

Figure 4

 4 Figure 4Above: P (τ > t|q b 0 = 4, q a 0 = 5) as a function of t for λ = 12, µ + θ = 13. Below: same figure in log-log coordinates. Note the Pareto tail which decays as t -2 .

Figure 5

 5 Figure 5 Above: Conditional probability of a price increase, as a function of the bid and ask queue size. Below: comparison with empirical transition frequencies for CitiGroup stock price tick-by-tick data on June 26, 2008.

Figure 6

 6 Figure 6 λ/D(f ), estimated from tick-by-tick order flow (vertical axis) vs standard deviation of 10-minute price increments (horizontal axis) for stocks in the Dow Jones Index, estimated from high frequency data on June 26, 2008. Each point represents one stock. Red line indicates the best linear approximation.

Theorem 2 .

 2 Let λ < θ + µ and assume that the probability distribution f satisfiesm(λ, θ + µ, f ) = ∞ i=1 ∞ j=1 m(λ, θ + µ, i, j)f (i, j) < ∞,where for all (x, y) ∈ (N * ) 2 , m(λ, θ + µ, x, y) =

Table 1

 1 Average number of orders in 10 seconds and number of price changes(June 26th, 2008).

		Average no. of Price changes orders in 10s in 1 day
	Citigroup	4469	12499
	General Electric	2356	7862
	General Motors	1275	9016

Proof. First, let us prove that P[X 2 = δ|X 1 = δ] = P[X 2 = -δ|X 1 = -δ]:

where p up 1 is given in 7, by symmetry of the bid and the ask, for all (n, p)

By assumption 1, for all (i, j) ∈ N 2 , f (i, j) = f (j, i). Therefore,

x, y) defined by ( 12) is then characterized by the following relation:

The eigenvalues of this matrix are 1 and 2p cont -1:

Moreover for all n ≥ 2,

. and the correlation between two consecutive price moves is given by: cov

Remark 2 (Negative autocorrelation of price changes at first lag). It is empirically observed that high frequency price movements have a negative autocorrelation at the first lag; this observation is often attributed to the 'bid-ask' bounce of transaction prices, but in fact it also holds for the time series of bid or ask prices [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF]. Our model links the value of this autocorrelation at first lag to the properties of the distribution f of order book depth. As observed from (13), the sign of the cov(X 1 , X 2 |q b 0 , q a 0 ) does not depend on the initial configuration (q b 0 , q a 0 ) of the bid/ask queues, so cov(X k , X k+1 ) = cov(X 1 , X 2 ) < 0 if and only if

f (i, j)p up 1 (i, j) < 1/2