Frédéric Abergel

Rémy Tachet

Des Combes

Riadh Zaatour

Nonparametric model calibration for derivatives

published or not. The documents may come

Introduction

One of the most important challenge for real-life applications of a model to derivatives trading is the issue of calibration. Similar to common situations in many areas of physics and engineering, once a model has been suggested, its parameters have to be estimated using external data. In the case of derivative modelling, those data are the liquid (tradable) options, generally called vanilla products. It is well-known since the pioneering work of [START_REF] Breeden | State Contingent Prices Implicit in Option Prices[END_REF] and its celebrated extension by [START_REF] Dupire | Pricing and Hedging with Smiles[END_REF] that the knowledge of market data such as the prices of vanilla options across all strikes and maturities is equivalent to the knowledge of the risk-neutral marginals of the underlying stock distribution. Here, we are interested in applications of this result to three different cases. The first one is quite classic, and is the starting point of our work: local and stochastic volatility models. Such models are very useful in practice, since they offer both the flexibility and realistic dynamics of stochastic volatility, and the exact calibration properties of local volatility. The problem of calibrating local and stochastic volatility models has been dealt with for a while now, for instance by [START_REF] Lipton | The vol smile problem[END_REF]. However, practitioners seem to agree that the stability of its resolution becomes uncertain when the volatility's volatility is too large. We shall see that indeed some kind of instability appears, and offer explanations to the phenomenon. The second case we focus on is the correlation between assets. Empirical measures give a certain set of results. However, when modelling a basket on multiple underlyings, a problem occurs. If one uses local volatility models for each underlying and correlates their brownian motions using the empirical correlation, the basket obtained will not reproduce the vanillas quoted on the market. This raises significant issues when hedging products on multiple underlyings. One of the solution for this problem is the known 'local correlation' approach: the correlation matrix for the n underlyings is deformed using a parameter, function of the time and the basket level. Here, we use that approach to obtain a calibration equation for the basket, relatively similar to the one appearing in the local and stochastic volatility model, and then numerically solve said equation in a two-underlyings framework. The last topic we shall be interested in are interest rates, we study a hybrid model: local volatility with stochastic rates. Using a partial differential equation approach similar to the local correlation and the local and stochastic volatility, we write a calibration equation for the vanillas of this hybrid model, solve it and verify the accuracy of the fit. The general form of our calibration equations is nonlinear partial and integro-differential. For their resolution, we chose to adapt the alternating direction implicit scheme (very efficient to solve classic linear second order parabolic equations, [START_REF] Richtmyer | Difference methods for initial value problems[END_REF]). Being in a nonlinear non-local framework, many questions arise. Is it relevant to use ADI algorithms to solve the equations stemming from our calibration problems? How should we deal with the nonlocal term? Is the finite difference scheme we chose consistent, and what is the order of the truncation error? Can we detect an instability in certain cases? Is it possible to explain it? The aim of this work is to address and at least partially answer these questions. The paper is organized as follows. In Section 2, we quickly present the case of the local and stochastic volatility model. Section 3 is devoted to the local correlation, its calibration equation and the fit we obtain in the case of a basket on two underlyings. In Section 4, we do the same thing in the stochastic rates frame. Finally, in Section 5 and 6, we adress the questions raised previously concerning the ADI algorithm used for the resolution, Section 7 is a brief conclusion.

Local and Stochastic Volatility models 2.1 Partial Integro-Differential Equation for the calibration of LSV models

The diffusion model is assumed to be the following Proof. The exact assumptions and the existence proof can be found in [START_REF] Tachet Des Combes | Non-parametric model calibration in finance[END_REF]. Here, the main concern is the calibration result. Let us assume that the model fits exactly the surface C. Letting (S 0 , y 0 ) denote the initial state of the system, the joint density p(t, S, y) of the couple (S t , y t ) verifies Kolmogorov forward equation

dS t S t = r(t)dt + a(t, S t )b(y t )dW 1 t ( 1
)
dy t = µ(t
∂p ∂t - ∂ 2 ∂S 2 ( 1 2 a 2 b 2 S 2 p) - ∂ 2 ∂S∂y (ρabξSp) - ∂ 2 ∂y 2 ( 1 2 ξ 2 p)) + ∂ ∂S (rSp) + ∂ ∂y (µp) = 0 p(0, S, y) = δ S0,y0
Applying Fubini, let q(t, S) = R p(t, S, y)dy be the first marginal density of our couple. It is possible to integrate the previous equation and obtain

∂q ∂t - ∂ 2 ∂S 2 ( 1 2 a 2 S 2 ( R b 2 (y)p(t, S, y)dy)) + ∂ ∂S (rSq) = 0 q(0, S) = δ S0
In the case of a local volatility model (b = 1 and a = σ D ), the density q D of the Spot process solves the equation

∂q D ∂t - ∂ 2 ∂S 2 ( 1 2 σ 2 D S 2 q D ) + ∂ ∂S (rSq D ) = 0 q D (0, S) = δ S0 3
The vanillas of the LSV model being perfectly fitted, we have q = q D . Identifying the terms in the two last formulas gives

a 2 (t, S) = σ 2 D (t, S) q R b 2 pdy = σ 2 D (t, S) R p(t, S, y)dy R b 2 (y)p(t, S,

y)dy

Using this proposition, and reintroducing the value of a in Kolmogorov forward equation, the joint density p(t, S, y) is then solution of the nonlinear partial integro-differential equation

∂p ∂t - ∂ 2 ∂S 2 ( 1 2 σ 2 D b 2 S 2 R pdy R b 2 pdy p) - ∂ 2 ∂S∂y (ρσ D bξS( R pdy R b 2 pdy ) 1 2 p) (4)
- ∂ 2 ∂y 2 ( 1 2 ξ 2 p)) + ∂ ∂S (rSp) + ∂ ∂y (µp) = 0 p(0, S, y) = δ S0,y0 (5)
There is thus equivalence between the existence of a model (1-2) that calibrates the vanillas C and the existence of a solution p to the pide (4).

Remark 2.1. The quotient R b 2 pdy R pdy is nothing but the conditional expectation of the volatility squared, knowing the spot process. This result is not original in itself (by applying the theorem from [START_REF] Gyongy | Mimicking the One-Dimensional Marginal Distributions of Processes Having an Ito Differential[END_REF] for instance), the partial differential equation method however is unusual, and will be used on the other models as well.

The theoretical study of equation ( 4)-( 5) can be found in [START_REF] Abergel | A nonlinear partial integrodifferential equation from mathematical finance[END_REF]. Existence of solutions is proved under strong assumptions (especially on b, which must be sufficiently close to a constant). The general resolution remains an open problem.

Numerical results

It seems to be well-known among practitioners, that instabilities occur in their calibration when the volatility's volatility (in the notations, function ξ) is too large. This seems to confirm the theoretical limitations met trying to prove the global existence of a solution: when the function b oscillates too much (a change of scale in the factor y t clearly shows the equivalence between a b that moves a lot and a large ξ), the resolution of the equation is not guaranteed anymore. To assess these statements, we considered our problem from a practical viewpoint. In this section, the calibration that stems from solving the partial differential equation ( 4) is studied, for two stochastic volatility models: a lognormal one and a 'Cox-Ingersoll-Ross' process. The details of the algorithm used for the resolution and a study of the instabilities will be treated later, in sections 5 and 6.

Lognormal volatility

Starting with a simple mean reverting model for the volatility factor, the function b is chosen as an exponential

dS t S t = r(t)dt + a(t, S t )exp(y t )dW 1 t ( 6
)
dy t = κ (δ -y t ) dt + γ dW 2 t ( 7
)
with a 2 (t, S) = σ 2 D (t, S) R p(t, S, y)dy R exp(2y)p(t, S, y)dy Equation ( 4) is solved using the functions we just chose and the local volatility σ D associated to the Eu-roSTOXX 50 implied volatility surface of 2009/04/02. Once function p, density of the couple (S t , y t ), is found, we compute the vanilla prices for different strikes and maturities using this density. To have a point of comparison, the same prices are also computed with the local volatility σ D , both of them are then compared to the targeted prices (column TP). Let us then plot the error between the original vanillas and the ones obtained with the model. The calibration is quite efficient, the errors are equivalent to the ones of the local volatility model.

Strike

Cox-Ingersoll-Ross process

We also focus on the calibration of a model inspired from the interest rates framework: the volatility is assumed to follow a CIR process.

dS t S t = r(t)dt + a(t, S t )y t dW 1 t ( 8
)
dy t = κ (α -y t ) dt + γ √ y t dW 2 t ( 9
)
Detailled properties of this process are described by [START_REF] Brigo | Credit default swap calibration and derivatives pricing with the SSRD stochastic intensity model[END_REF]. In particular, as long as 2κα > γ 2 , y t is strictly positive a.e. Once again, equation ( 4) is solved with this stochastic volatility.

Application to the 'Local Correlation' model

In this chapter, we are interested in the calibration of a market with n stocks and a basket on those stocks. The purpose is to define a diffusion model for those underlyings that is able to reproduce their implied volatility surface as well as the one of the basket. The notations are the following, let (S i t ) 1≤i≤n denote the n stocks involved in our problem. The basket's value is given by [START_REF] Jourdain | Coupling Index and stocks[END_REF] where the set (w i ) 1≤i≤n stands for the weights of the different underlyings. They are assumed to be constant in the rest of our work. Let us also fix n + 1 surfaces of vanillas (C i (T, K)) 1≤i≤n and C B (T, K).

B t = n i=1 w i S i t

Inconsistencies between stock and basket options

The naive approach to solve this problem is simply to consider n local volatility models

dS i t S i t = r(t)dt + σ i (t, S i t )dW i t (11)
The functions σ i are easily determined to fit the surfaces (C i (T, K)) 1≤i≤n with this diffusion. The correlation matrix ρ = (ρ ij ) 1≤i,j≤n associated to the standard brownian motions W i t of each underlying can be estimated with historical data. The model is now entirely defined. By equation ( 10) of B t , the vanilla prices for the basket are completely determined and are equal to E[(

w i S i T -K) + ].
However, there is no particular reason for the surface computed in this framework to be equal to C B (T, K). In fact, the skew of the basket is more pronounced on the market than in a model with constant correlations between the underlyings, [START_REF] Qu | Pricing Basket Options With Skew[END_REF].

'Local Correlation' model

In the manner of B. Dupire who decided to let the volatility depend on the level of the spot process, a degree of freedom is added to our model by distorting the matrix of correlation with a function of B t . This method appears for instance in [START_REF] Tachet Des Combes | Non-parametric model calibration in finance[END_REF][START_REF] Jourdain | Coupling Index and stocks[END_REF] where the basket level induces some feedback on the values of the different underlyings. In our context, the new correlation matrix is taken as a linear combination of ρ and the constant matrix with only 1 as coefficients. This matrix ρ is equal to

ρij = λ + (1 -λ)ρ ij = ρ ij + λ(1 -ρ ij ) (12)
We shall see while writing the calibration equation that λ has to be chosen as a function of the time and of B t . The matrix ρij can be seen as an analogous of Dupire's local volatility, a 'Local Correlation' so to speak. Assuming that ρ is a proper correlation matrix (definite, positive...), and that the coefficients of the diffusion have a sufficient regularity, our diffusion model then possesses a density in the more general case of a matrix ρ function of the couple (t, B). It is also possible to write a condition for the vanillas of the model to be fitted Proposition 2. The diffusion model defined by [START_REF] Douglas | On the numerical solution of heat conduction problems in two and three space variables[END_REF]

w i w j ρ ij (t, B) ( βi βj p)(t, B, S 2 , .., S n )dS 2 ..dS n (13)
where

p(t, B, S 2 , .., S n ) = p(t, S 1 , S 2 , .., S n ) βi (t, B, S 2 , .., S n ) = S i σ i (t, S i ) with S 1 = 1 w 1 (B - n i=2 w i S i ) (14)
Proof. The existence of the transition density p(t, S 1 , .., S n ) stems from the assumptions made on the regularity of the coefficients, and on ρ, for more details see [START_REF] Tachet Des Combes | Non-parametric model calibration in finance[END_REF]. We can now write the calibration problem for the vanillas of the basket B t . The density just defined satisfies Kolmogorov forward equation

∂p ∂t - 1 2 1≤i,j≤n ∂ 2 ∂S i ∂S j (ρ ij S i σ i S j σ j p) + 1≤i≤n ∂ ∂S i (rS i p) + rp = 0
To ease the problem, it is useful to change the coordinates (S 1 , ..., S n ) into (B, S 2 , ..., S n ) with S 1 defined by [START_REF] Beam | Alternating direction implicit methods for parabolic equations with a mixed derivative[END_REF]. After computations, the equation becomes

∂ p ∂t - 1 2
1≤i,j≤n

w i w j ∂ 2 ∂B 2 (ρ ij βi βj p) - 1 2 2≤i,j≤n ( ∂ 2 ∂S i ∂S j + w i ∂ 2 ∂B∂S j + w j ∂ 2 ∂S i ∂B )(ρ ij βi βj p) - 2≤i≤n w 1 ∂ 2 ∂B∂S j (ρ 1i β1 βi p) + ∂ ∂B (rB p) + 2≤i≤n ∂ ∂S i (rS i p) + r p = 0
where p(B, S 2 , ..., S n ) and βi (B, S 2 , ..., S n ) are defined above. Integrating the equation against the variables (S 2 , ..., S n ), and writing q = 1 w1 pdS 2 ..dS n , the density of the marginal law of B satisfies

∂q ∂t - 1 2
1≤i,j≤n

w i w j ∂ 2 ∂B 2 (ρ ij βi βj pdS 2 ..dS n ) + ∂ ∂B (
∂ 2 ∂B 2 (σ 2 D B 2 q) + ∂ ∂B (rBq) + rq = 0
if our model reproduces the vanillas σ D , then the following equality must be verified

σ 2 D B 2 pdS 2 ..dS n = 1≤i,j≤n w i w j ρ ij βi βj pdS 2 ..dS n
Reciprocally, the condition we just wrote is clearly sufficient for the options to be calibrated Remark 3.1. Let us note that this condition is written as an equality between two functions of the time and of B. The other variables are no longer represented.

Assuming that condition [START_REF] Douglas | A General Formulation of Alternating Direction Methods[END_REF] is not verified, the model defined by [START_REF] Douglas | On the numerical solution of heat conduction problems in two and three space variables[END_REF] does not fit the vanillas of the basket B t . It has to be enriched to solve the calibration problem. Our choice is to distort the correlation matrix.

The new matrix ρ is described by [START_REF] Douglas | Alternating direction methods for three space variables[END_REF]. Hence, let Θ denote the matrix Θ ij = 1 for all 1 ≤ i, j ≤ n. We also notice that, the trace of ρ being equal to n, its smallest eigenvalue K ρ is smaller than 1.

Lemma 1. The matrix ρ = (1 -λ)ρ + λΘ is also a correlation matrix as long as λ is in ]r, 1[ with r = -min(max

i =j 1 + ρ ij 1 -ρ ij , K ρ n -K ρ ) < 0 ( 15
)
Proof. Clearly, for all (ξ i , ξ j ) ∈ R 2 1≤i,j≤n

ξ i ((1 -λ)ρ ij + λ)ξ j ≥ K ρ (1 -λ)|ξ| 2 + λ( 1≤i≤n ξ i ) 2
If λ is positive, since K ρ (1 -λ) is stricty positive, the matrix remains definite positive. Now, if λ < 0, Cauchy-Schwarz gives 1≤i,j≤n

ξ i ((1 -λ)ρ ij + λ)ξ j ≥ (K ρ (1 -λ) + λn)|ξ| 2 Since K ρ ≤ 1, λ > - Kρ n-Kρ is enough for K ρ (1 -λ)
+ λn to be stricly positive. The diagonal coefficients of ρ are still 1. As for the other terms, thanks to the first term in relation [START_REF] Alibaud | Existence, uniqueness and regularity for nonlinear parabolic equations with nonlocal terms equations with nonlocal terms[END_REF], they still belong to the interval ] -1, 1[. We introduce the new correlation matrix in condition [START_REF] Douglas | A General Formulation of Alternating Direction Methods[END_REF], this gives

λ(p) = σ 2 D B 2 pdS 2 ..dS n - 1≤i,j≤n w i w j ρ ij βi βj pdS 2 ..dS n 1≤i,j≤n w i w j (1 -ρ ij ) βi βj pdS 2 ..dS n (16)
λ(p) is a function of B and t. It is now possible to use this value to write a pide on the density p. Any solution of the following equation is a density that calibrates the vanillas of the basket

∂ p ∂t + L 1 p + L λ( p) 2 (p) = 0 (17)
where L 1 is linear and verifies

L 1 p = - 1 2 1≤i,j≤n w i w j ρ ij ∂ 2 ∂B 2 ( βi βj p) - 1 2 2≤i,j≤n ρ ij ( ∂ 2 ∂S i ∂S j + w i ∂ 2 ∂B∂S j + w j ∂ 2 ∂S i ∂B )( βi βj p) - 2≤i≤n w 1 ρ 1i ∂ 2 ∂B∂S j ( β1 βi p) + ∂ ∂B (rB p) + 2≤i≤n ∂ ∂S i (rS i p) + r p and L λ( p) 2
is the nonlinear part of the equation

L λ( p) 2 (p) = - 1 2 1≤i,j≤n w i w j (1 -ρ ij ) ∂ 2 ∂B 2 (λ(p) βi βj p) - 1 2 2≤i,j≤n (1 -ρ ij )( ∂ 2 ∂S i ∂S j + w i ∂ 2 ∂B∂S j +w j ∂ 2 ∂S i ∂B )(λ(p) βi βj p) - 2≤i≤n w 1 (1 -ρ 1i ) ∂ 2 ∂B∂S j (λ(p) β1 βi p)
Remark 3.2. The operator L 1 + L 2 stems from a change of coordinates on a uniformly elliptic operator. It is also elliptic, uniformly on any domain where the βi are bounded away from 0 by a strictly positive constant.

Furthermore, the initial condition is

p(0, B, S 2 , ..., S n ) = δ( w i S i 0 , S 2 0 , ..., S n 0 )
where S i 0 is the market value at instant 0 of the i-th stock. Applying this initial condition to (16), the initial value of λ is equal to

λ(p)(0, B) = σ 2 D (0, B)B 2 - 1≤i,j≤n w i w j ρ ij S i 0 σ i (0, S i 0 )S j 0 σ j (0, S j 0 )
1≤i,j≤n

w i w j (1 -ρ ij )S i 0 σ i (0, S i 0 )S j 0 σ j (0, S j 0 ) (18)
For a theoretical study of the calibration equation, we refer the reader to [START_REF] Tachet Des Combes | Non-parametric model calibration in finance[END_REF].

Resolution of the equation for the calibration of a basket

This subsection focuses on the results of the calibration for a two-underlyings basket. Let us consider two assets, both of them are assumed to generate the following implied volatility surface KT 0,088 0,167 0,25 0,50 0,75 Here are some results for other tests. We calibrate a model with the following parameters ρ 12 = 0, w 1 = 0.3 and w 2 = 0.7, the targeted surface is the theoretical one distorted with two factors: first 0.95 and second 1.05. At last, a different surface for the second underlying is chosen, mutliplying the first one (described in 3.3) by 0.9, the correlation is this time taken as 0. The results are rather satisfactory, especially at the money. To keep the computations to a reasonnable duration, a sparse initial surface was used. This explains why the calibration is not better far from the money, the fitting method is nontheless valid. Now follows an outlook of the values taken by the new correlation ρ at different maturities when the theoretical surface is distorted by factors 0.95 and 1.05. The parameters are: ρ 12 = 0.5, w 1 = 0.7 and w 2 = 0.3. As expected, the Local Correlation and the distorsion factor evolve in the same direction. The underlyings must be more correlated when the implied volatility of the basket is higher, and reciprocally. Furthermore, it appears that in the case of the 0.95 distorsion, the correlation has to violently decrease for high values of B: the two underlyings must be anti-correlated when they are both large. As for the influence of the maturity, let us first state that in the computation of λ, when the denominator is smaller than 10 -6 , we chose not to change the correlation, to avoid numerical errors. It appears that, as long as B is in a zone where λ was actually computed, the framework chosen to test the calibration actually generates a local correlation constant in time.

Application to stochastic interest rates

This section is dedicated to hybrid local volatility models with stochastic rates. The interest rate is assumed to be stochastic and to follow a diffusion equation. The volatility depends on the level of the spot process exactly as in a local volatility model. The idea is to compute its exact value for the vanillas in this model to be calibrated.

Calibration of the hybrid local volatility model

The risk-neutral diffusion of the model is written as The two brownian motions are correlated with a constant correlation denoted by ρ. Classic regularity and ellipticity assumptions are made on the coefficients of the diffusion (described in [START_REF] Tachet Des Combes | Non-parametric model calibration in finance[END_REF]) to get the Proposition 3. The diffusion model defined above has a transition density with respect to Lebesgue's measure. The value of σ that fits its vanillas is given by This equation needs to be matched with

dS t S t =
σ 2 (t, S) = σ 2 D (t, S) + 2 r(t) R +∞ S p(t,
∂q D ∂t - ∂ 2 ∂S 2 ( 1 2 σ 2 D S 2 q D ) + ∂ ∂S (rSq D ) + rq D = 0 q D (0, S) = δ S0
Both of them can be written as In order to calibrate the vanillas, all that remains to be done is match the marginal density q with q D , giving the necessary condition

1 2 σ 2 S 2 = +∞ 0 (s -S) + ( ∂q ∂t + ∂ ∂s (s rpdy) + rpdy)ds q 1 2 σ 2 D S 2 = +∞ 0 (s -S) + ( ∂q D ∂t + ∂ ∂s (srq D ) + rq D )ds q D Computing +∞ 0 (s -S) + ∂ ∂s ( s
1 2 σ 2 S 2 + S +∞ S rpdyds q = 1 2 σ 2 D S 2 + rS +∞ S

qds q

Replacing q by R p(t, S, y)dy completes the proof. The calibration equation for the vanillas of our hybrid model is thus

∂p ∂t - ∂ 2 ∂S 2 ( 1 2 σ 2 S 2 p) - ∂ 2 ∂S∂y (ρσλSp) - ∂ 2 ∂y 2 ( 1 2 λ 2 p)) + ∂ ∂S (rSp) + ∂ ∂y (µp) + rp = 0 (21)
with σ given by formula (20). Using similar technics to the other cases, an existence result can be obtained under certain assumptions, but this is not the scope of this paper. It is however noteworthy that one of the necessary hypothesis is the small variation of function r with respect to the deterministic curve r.

Numerical calibration

In this section, the theoretical results above are applied to calibrate a given diffusion model. Assuming the instantaneous rate to obey a Vasicek model (or in other words an Ornstein-Uhlenbeck process), the diffusion equations become

dS t S t = r t dt + σ(t, S t )dW 1 t dr t = a(b -r t )dt + γdW 2 t
The ADI algorithm described in section 5 is applied to equation ( 21) with the coefficients associated to this diffusion. The initial condition is p 0 (S, r) = δ S0,r0 .

As in the two previous sections, this partial integro-differential equation is solved with a variable change for the spot process x = ln(S). The grid chosen is [-10σ √ t, 10σ √ t] × [-0.1, 0.2] with σ = 0.2 and r 0 = 0.04. We discretize it with 300 points in both the spot and the rate direction. The initial condition (Dirac mass at the point (ln(S 0 ), r 0 )) is approximated by a bivariate Gaussian centred at that point with a very small variance.

The following numerical values are taken for the diffusion To assess the quality of the calibration, call and put options on the spot process are computed by integration on the grid and compared to the targeted prices (target columns). Convergence is quite satisfactory. For instance, for 6 months and 1 year maturity vanillas, one finds

Algorithm for the resolution of the calibration equation

In the previous sections, we wrote the calibration equations and gave graphs for their efficiency. In this one, we describe the algorithm used to solve them: a classic alternating direction implicit scheme. The nonlinear term is handled using a forward induction at first, and then a predictor-corrector method.

The strong feature of an ADI scheme is its convergence rate in time and space: O(δx 2 ) + O(δt 2 ). The nonlinearity of the equation challenges this assertion. It is however possible to prove that it remains true in this case too.

Alternating Direction Implicit scheme

The calibration equation is a second-order parabolic equation. One of the most efficient method to solve such equations is a finite-difference approximation with alternating direction methods. For more informations on the subject, the reader can look into [START_REF] Richtmyer | Difference methods for initial value problems[END_REF], numerous articles have also been published, in particular [START_REF] Douglas | On the numerical solution of heat conduction problems in two and three space variables[END_REF][START_REF] Douglas | Alternating direction methods for three space variables[END_REF][START_REF] Douglas | A General Formulation of Alternating Direction Methods[END_REF][START_REF] Beam | Alternating direction implicit methods for parabolic equations with a mixed derivative[END_REF]. Let us now consider the following equation

∂p ∂t - ∂ 2 ∂x 2 ( 1 2 f 2 I 2 (p)p) - ∂ 2 ∂x∂y (ρf gI(p)p) - ∂ 2 ∂y 2 ( 1 2 g 2 p)) + ∂ ∂x (αp) + ∂ ∂y (βp) + γp = 0 (22) p(0, S, y) = δ S0,y0 (23)
where f, g, α, β and γ are functions of t, x and y, ρ a constant and I(p) the quotient of integrals

I 2 (p)(t, x) = h 2 (y)p(t,
x, y)dy p(t, x, y)dy (24)

Remark 5.1. We restricted ourselves to two-dimensional equations since they cover all the concrete examples studied previously. But the computations that follow are true in the general case. The cost in time however becomes an issue in higher dimensions.

The domain for the numerical resolution is ]0, T [×]x * , x * [×]y * , y * [. The first step is to take care of the initial condition. Instead of the Dirac, the initial condition p 0 is chosen as a gaussian distribution with very small variance. It obviously approximates our initial condition. It also verifies (on any bounded domain) the properties of regularity and strict positivity required in the theoretical study of the equation (though in the present section, we are only interested in its numerical resolution). Let p(t, x, y) = p 0 (x, y) if x = x * , x * or y = y * , y * be the boundary conditions.

The algorithm is based upon a predictor-corrector approach. Let ∆x, ∆y and ∆t be increments of the variables x, y and t, where ∆x = x * -x * I , ∆y = y * -y * J and ∆t = T N with I, J and N integers. The sets of points in the x,y,t-plane is given by x n = x * + i ∆x

x * -x * , y j = y * + j ∆y y * -y * and t n = n ∆t T , for 0 ≤ i ≤ I, 0 ≤ j ≤ J and 0 ≤ n ≤ N . We construct four sequences p n , p * n , q n and q * n of space-dependent functions with n between 0 and N.

A classic alternating direction implicit scheme functions as follows: let us define the initial functions p * 0 (i, j) = q 0 (i, j) = q * 0 (i, j) = p 0 (x i , y j ) and then by induction

q * n+1 -p n ∆t = δ 2 x (I 2 (p n ) f 2 n+1 q * n+1 + f 2 n p n 4 ) + δ 2 xy (ρf n g n I(p n )p n ) + δ 2 y ( 1 2 g 2 n p n ) -δ x ( α n+1 q * n+1 + α n p n 2 ) -δ y (β n p n ) - γ n+1 q * n+1 + γ n p n 2 (25)
q n+1 -p n ∆t = δ 2 x (I 2 (p n ) f 2 n+1 q * n+1 + f 2 n p n 4 ) + δ 2 xy (ρf n g n I(p n )p n ) + δ 2 y ( g 2 n+1 q n+1 + g 2 n p n 4 ) -δ x ( α n+1 q * n+1 + α n p n 2 ) -δ y ( β n+1 q n+1 + β n p n 2 ) - γ n+1 q n+1 + γ n p n 2 ( 26
)
where f n (i, j) designates f (n∆t, x i , y j ) (the same thing being true for the other coefficients of the equation).

δ is a difference operator for the space derivatives. For instance, with 1 ≤ i ≤ I -1 and 1 ≤ j ≤ J -1

δ 2 x f n = f n (i -1, j) -2f n (i, j) + f n (i + 1, j) ∆x 2 δ 2 xy f n = f n (i -1, j -1) -f n (i + 1, j) -f n (i, j + 1) + f n (i + 1, j + 1) 4∆x∆y
δ x f n = f n (i + 1, j) -f n (i -1, j) 2∆x . . .
The equations ( 25) and (26) form two tridiagonal systems that can be solved very efficiently. A recursion formula can be computed on the functions q n defined above

q n+1 -q * n+1 ∆t = δ 2 y ( g 2 n+1 q n+1 -g 2 n p n 4 ) -δ y ( β n+1 q n+1 -β n p n 2 ) -γ n+1 q n+1 -q * n+1 2 Thus q * n+1 = q n+1 - 2∆t 2 + γ n+1 ∆t (δ 2 y ( g 2 n+1 q n+1 -g 2 n p n 4 ) -δ y ( β n+1 q n+1 -β n p n 2 ))
And eventually

q n+1 -p n ∆t = δ 2 x (I 2 (p n ) f 2 n+1 q n+1 + f 2 n p n 4 ) + δ 2 xy (ρf n g n I(p n )p n ) + δ 2 y ( g 2 n+1 q n+1 + g 2 n p n 4 ) -δ x ( α n+1 q n+1 + α n p 2 ) -δ y ( β n+1 q n+1 + β n p n 2 ) - γ n+1 q n+1 + γ n p n 2 - ∆t 4 + 2γ n+1 ∆t δ 2 x (I 2 (p n )f 2 n+1 (δ 2 y ( g 2 n+1 q n+1 -g 2 n p n 4 ) -δ y ( β n+1 q n+1 -β n p n 2 ))) + ∆t 2 + γ n+1 ∆t δ x (α n+1 (δ 2 y ( g 2 n+1 q n+1 -g 2 n p n 4 ) -δ y ( β n+1 q n+1 -β n p n 2 ))) (27)
In the litterature concerning alternating direction implicit schemes, the functions q n are often called the predicted value of the solution p, a second 'corrector' step usually follows. In our case, the equation being nonlinear, we start by chosing p n+1 = q n+1 and study the finite-difference approximation that results.

Proposition 4. This last finite-difference equation is consistent with the partial differential equation ( 22) on a bounded domain with smooth initial condition p 0 , the truncation error is

O(∆t) + O(∆x 2 ) + O(∆y 2 ) + O( ∆x 3 ∆y ).
Proof. Let p be a classic solution of equation ( 22) on a bounded domain, with p |t=0 = p 0 . We assume that p is strictly positive and sufficiently differentiable for all the quantities in the sequel to be properly defined. All the derivatives that appear are bounded as a consequence of the regularity assumptions on p. Writing

p n i,j = p n (x i , y j ) = p(t n , x i , y j ),
∂ 3 p n ∂x 3 (x i + θ 4 , y j-1 ) + ∂ 3 p n ∂x 3 (x i -θ 5 , y j+1 ) - ∂ 3 p n ∂x 3 (x i -θ 6 , y j+1 )] δ x p n = ∂p n ∂x (x i , y j ) + ∆x 2 12 ( ∂ 3 p n ∂x 3 (x i + θ 1 , y j ) + ∂ 3 p n ∂x 3 (x i -θ 2 , y j ))
where the different constants θ are between 0 and ∆t, ∆x or ∆y depending on the context, they may change from one formula to another. Let E denote the truncation error for scheme (27), we have

E = p n+1 -p n ∆t -δ 2 x (I 2 (p n ) f 2 n+1 p n+1 + f 2 n p n 4 ) -δ 2 xy (ρf n g n I(p n )p n ) 2 y ( g 2 n+1 p n+1 + g 2 n p n 4 ) + δ x ( α n+1 p n+1 + α n p n 2 ) + δ y ( β n+1 p n+1 + β n p n 2 ) + γ n+1 p n+1 + γ n p n 2 + ∆t 4 + 2γ n+1 ∆t δ 2 x (I 2 (p n )f 2 n+1 (δ 2 y ( g 2 n+1 p n+1 -g 2 n p n 4 ) -δ y ( β n+1 p n+1 -β n p n 2 ))) - ∆t 2 + γ n+1 ∆t δ x (α n+1 (δ 2 y ( g 2 n+1 p n+1 -g 2 n p n 4 ) -δ y ( β n+1 p n+1 -β n p n 2 )))
Applying the previous expansions and using the fact that p verifies equation ( 22

(p n )f 2 n p n )(x i + θ 1 , y j ) + ∂ 4 ∂x 4 (I 2 (p n )f 2 n p n )(x i -θ 2 , y j ) + ∂ 4 ∂x 4 (I 2 (p n )f 2 n+1 p n+1 )(x i + θ 3 , y j ) + ∂ 4 ∂x 4 (I 2 (p n )f 2 n+1 p n+1 )(x i -θ 4 , y j )] e 22 = ∆y 2 12 ( ∂ 4 ∂x∂y 3 (ρf n g n I(p n )p n )(x i , y j + θ 1 ) + ∂ 4 ∂x∂y 3 (ρf n g n I(p n )p n )(x i , y j -θ 2 )) + ∆x 2 24∆y ( ∂ 3 ∂x 3 (ρf n g n I(p n )p n )(x i + θ 3 , y j-1 ) - ∂ 3 ∂x 3 (ρf n g n I(p n )p n )(x i + θ 4 , y j-1 ) + ∂ 3 ∂x 3 (ρf n g n I(p n )p n )(x i -θ 5 , y j+1 ) - ∂ 3 ∂x 3 (ρf n g n I(p n )p n )(x i -θ 6 , y j+1 )) e 23 = ∆y 2 96 [ ∂ 4 ∂y 4 (g 2 n p n )(x i , y j + θ 1 ) + ∂ 4 ∂y 4 (g 2 n p n )(x i , y j -θ 2 ) + ∂ 4 ∂y 4 (g 2 n+1 p n+1 )(x i , y j + θ 3 ) + ∂ 4 ∂y 4 (g 2 n+1 p n+1 )(x i , y j -θ 4 )]
e 3 enables us to compensate for both the nondiagonal term and the nonlocal term I(p n ) that cannot be computed implicitely

e 3 = ∂ 2 ∂x 2 [(I 2 (p n ) -I 2 (p n+1 )) f 2 n+1 p n+1 4 ] + ∂ 2 ∂x∂y [ ρ 2 (f n g n I(p n )p n -f n+1 g n+1 I(p n+1 )p n+1 )]
e 41 and e 42 are the terms corresponding to the first order space-derivatives

e 41 = ∆x 2 24 [ ∂ 3 ∂x 3 (α n p n )(x i + θ 1 , y j ) + ∂ 3 ∂x 3 (α n p n )(x i -θ 2 , y j ) + ∂ 3 ∂x 3 (α n+1 p n+1 )(x i + θ 3 , y j ) + ∂ 3 ∂x 3 (α n+1 p n+1 )(x i -θ 4 , y j )] e 42 = ∆y 2 24 [ ∂ 3 ∂y 3 (β n p n )(x i , y j + θ 1 ) + ∂ 3 ∂y 3 (β n p n )(x i , y j -θ 2 ) + ∂ 3 ∂y 3 (β n+1 p n+1 , y j + θ 3 ) + ∂ 3 ∂y 3 (β n+1 p n+1 )(x i , y j -θ 4 )]
At last, e 5 is the correction term stemming from the alternating direction implicit scheme

e 5 = ∆t 4 + 2γ n+1 ∆t δ 2 x (I 2 (p n )f 2 n+1 (δ 2 y ( g 2 n+1 p n+1 -g 2 n p n 4 ) -δ y ( β n+1 p n+1 -β n p n 2 ))) - ∆t 2 + γ n+1 ∆t δ x (α n+1 (δ 2 y ( g 2 n+1 p n+1 -g 2 n p n 4 ) -δ y ( β n+1 p n+1 -β n p n 2 )))
Thanks to the regularity of p and of the coefficients of ( 22), the upper bounds e 1 ≤ K∆t 2 , e 21 +e 41 ≤ K∆x 2 , e 23 + e 42 ≤ K∆y 2 and e 22 ≤ K(∆y 2 + ∆x 3 ∆y ) are obtained. It is also easily proven that e 5 ≤ K∆t 2 . The term that prevents us from getting an error in O(∆t 2 ) is e 3 . All we have is e 3 ≤ K∆t. This concludes the proof.

Remark 5.2. In a case with no I(p) term, the equation is a classic linear and parabolic one. In that case, when the off-diagonal 1 term is absent (ρ = 0 for instance), the previous scheme has an error in O(∆t 2 ). To obtain such an error in the general case, a second 'corrector' step is generally used: the predicted value q n+1 is introduced as an approximation of p n+1 in the cross-derivatives. Here, we use it in the nonlocal term too.

The correction step is the following

p * n+1 -p n ∆t = δ 2 x ( I 2 (q n+1 )f 2 n+1 p * n+1 + I 2 (p n )f 2 n p n 4 ) + δ 2 xy (ρ (f g) n+1 I(q n+1 )q n+1 + (f g) n I(p n )p n 2 ) + δ 2 y ( 1 2 g 2 n p n )) -δ x ( α n+1 p * n+1 + α n p n 2 ) -δ y (β n p n ) - γ n+1 p * n+1 + γ n p n 2 p n+1 -p n ∆t = δ 2 x ( I 2 (q n+1 )f 2 n+1 p * n+1 + I 2 (p n )f 2 n p n 4 ) + δ 2 xy (ρ (f g) n+1 I(q n+1 )q n+1 + (f g) n I(p n )p n 2 ) + δ 2 y ( g 2 n+1 p n+1 + g 2 n p n 4 ) -δ x ( α n+1 p * n+1 + α n p n 2 ) -δ y ( β n+1 p n+1 + β n p n 2 ) - γ n+1 p n+1 + γ n p n 2
Here too, one can compute p

* n+1 p * n+1 = p n+1 - 2∆t 2 + γ n+1 ∆t (δ 2 y ( g 2 n+1 p n+1 -g 2 n p n 4 ) -δ y ( β n+1 p n+1 -β n p n 2
((f g) n+1 I(q n+1 )q n+1 )(x i , y j + θ 1 ) + ∂ 4 ∂x∂y 3 ((f g) n+1 I(q n+1 )q n+1 )(x i , y j -θ 2 )) + ∆x 2 24∆y ( ∂ 3 ∂x 3 ((f g) n+1 I(q n+1 )q n+1 )(x i + θ 3 , y j-1 ) - ∂ 3 ∂x 3 ((f g) n+1 I(q n+1 )q n+1 )(x i + θ 4 , y j-1 ) + ∂ 3 ∂x 3 ((f g) n+1 I(q n+1 )q n+1 )(x i -θ 5 , y j+1 ) - ∂ 3 ∂x 3 ((f g) n+1 I(q n+1 )q n+1 )(x i -θ 6 , y j+1 )) e * 5 = ∆t 4 + 2γ n+1 ∆t δ 2 x (I 2 (q n+1 )f 2 n+1 (δ 2 y ( g 2 n+1 p n+1 -g 2 n p n 4 ) -δ y ( β n+1 p n+1 -β n p n 2 ))) - ∆t 2 + γ n+1 ∆t δ x (α (δ 2 y ( g 2 n+1 p n+1 -g 2 n p n 4 ) -δ y ( β n+1 p n+1 -β n p n 2 )))
which also verify e * 22 ≤ K(∆y 2 + ∆x 3 ∆y ) and e * 5 ≤ K∆t 2 . The real difference can be seen in

e * 3 = ∂ 2 ∂x 2 [(I 2 (q n+1 ) -I 2 (p n+1 )) f 2 n+1 p n+1 4 ] + ∂ 2 ∂x∂y [ ρ 2 ((f g) n+1 I(q n+1 )q n+1 -(f g) n+1 I(p n+1 )p n+1 )]
The important feature of the predictor is that the difference q n+1 -p n+1 is O(∆t 2 ). This gives e * 3 ≤ K(∆t 2 ) and concludes the proof.

Time Convergence Rate of the modified ADI algorithm

In this brief section, we compare the convergence of the algorithm with the theoretical rates computed in the previous part. To do so, the calibrated value of 1-year at-the-money vanillas is computed for different number N of time steps. We then plot the error between this price and the targeted value against N. The next graph is obtained with the one-step predictor algorithm

The error is clearly in O(∆t) as was proved in Proposition 4. We conduct the same experiment with this time both the predictor and the corrector steps. This time too, numerical experiments seem to agree with theory. The error appears to be in O(∆t 2 ). The predictor/corrector scheme serves its purpose.

Instabilities of the solutions: numerical explosion for 'oscillating' volatilities

In the previous sections, we described the algorithm used to solve the different calibration equations concerned by our work, and also mentioned the theoretical research performed on the subject. Though a partial existence result for equation ( 4) was found, it was still impossible to prove it in the general case: a strongly variable function b. This last section is devoted to the local and stochastic volatility model, from that point of view. The numerical resolution of the calibration is performed for strongly variable functions b.

At first, let us go back to the equation for the calibration of a local and stochastic volatility model. For the sake of simplicity, the interest rate is assumed to be zero. The equation is the following The solution is not smooth anymore. On the contrary, some kind of instability seems to occur.

∂p ∂t - ∂ 2 ∂S 2 (
Remark 6.1. To check that no other numerical effects are involved in the instability, the adjoint equation of ( 4) was also studied. From a theoretical viewpoint, it admits a solution without any restrictive assumption on b, see [START_REF] Alibaud | Existence, uniqueness and regularity for nonlinear parabolic equations with nonlocal terms equations with nonlocal terms[END_REF][START_REF] Tachet Des Combes | Non-parametric model calibration in finance[END_REF]. On both these graphics (Figure 2), no sign of instability can be seen. With higher values of C, for instance C = 20, the function p computed numerically takes meaningless values (10 80 ), but at no time does it start to oscillate.

[] []

Conclusion

Using methods inspired from local and stochastic volatility models, we were able to write calibration equations for two other cases: the so-called local correlation model and a hybrid local volatility and stochastic rates model. Their numerical resolution, based upon an alternating direction implicit scheme, produces a satisfactory fit under certain assumptions (confirmed by the theoretical difficulties met when studying them). When those hypothesis are not verified, an instability occurs. Explaining it brought us to consider Hadamard stability and a certain class of integro-differential equations. Unfortunately, the criterion we found can not be applied in the case of the local and stochastic volatility model, it remains an open problem.

As far as the ADI algorithm is concerned, we managed to adapt it to deal with the nonlinearity of our equation. Its consistency was also proved, with a convergence rate in time in O(∆t 2 ). A result confirmed numerically.

  rBq) + rq = 0 Comparing this equation to Dupire's equation for the local volatility σ B ∂q ∂t -1 2

  rpdy)ds = +∞ 0 (s -S) + (s ∂ ∂s ( rpdy) + rpdy)ds =where the second line stems from a simple integration by parts. Reintroducing this into the previous equationsS) + ∂q D ∂t ds -rS +∞ S q D ds q D

  S, y) = δ S0,y0The existence of a solution to this equation (on a bounded domain, with regularized boundary conditions) was previously stated under certain assumptions on b, the essential one being: there exists a constant b * such that |b -b(y 0 )| ≤ b * for a given y 0 . Using the algorithm described in the previous section, we are now interested in the behavior of the numerical solution of (28) when the function b violates the assumption. The model chosen for this study is the mean-reverting volatility already expounded b(y) = exp(y) α(t, y) = γ β(t, y) = κ (δ -y) with γ, κ and δ three strictly positive constants. The values chosen in section 2 for the different parameters of the model and of the algorithm led us to a satisfactory calibration. Using them once again, we plot on the left (Figure 1 (a)) the density p(T, x, y) for T = 1 year and x = ln(S/S 0 ) close to 0. As expected, the solution p is perfectly smooth. Now, bouncing on the idea of strongly variable functions b, another density is plotted on the right (Figure 1 (b)) with this time a function b equal to b(y) = exp(10y).

Figure 1 :Figure 2 :

 12 Figure 1: Solution p of the equation for b(y) = exp(y) and b(y) = exp(10y).

  with a correlation function of the couple (t, B), has a density p(t, S 1 , .., S n ) with respect to Lebesgue's measure. Furthermore, this model calibrates the surface C B (T, K) of B t 's vanillas (represented by its local volatility surface σ B ) if and only if σ 2 B (t, B)B 2 p(t, B, S 2 , .., S n )dS 2 ..dS n =

	1≤i,j≤n

  Carlo simulation and the local volatilities stemming from those surfaces, the theoretical prices for the basket B t are computed, with weights w 1 = w 2 = 0.5 and a correlation ρ 12 = -0.5 Distorting this theoretical surface by a factor of 0.9, and thus making the prices of the basket inconsistent with the prices of the underlyings, the calibration algorithm is applied. Solving the partial integro-differential equation (17) gives the following vanillas (quoted in implied volatility)

	KT	0,25			0,50			0,75	1,0
	Theor.	LC	Theor.	LC	Theor.	LC	Theor.	LC
	80% 19,62% 19,77% 18,65% 19,03% 18,00% 18,55% 17,65% 18,14%
	90% 18,66% 18,83% 17,80% 18,21% 17,40% 17,91% 17,13% 17,70%
	100% 17,85% 17,66% 17,24% 17,18% 17,02% 17,08% 16,81% 16,99%
	110% 16,89% 16,66% 16,72% 16,54% 16,63% 16,34% 16,44% 16,35%
	120% 16,23% 16,02% 16,29% 16,24% 16,34% 16,11% 16,24% 16,26%
									1,0	1,1
		80% 48,0% 47,9% 46,5% 45,0% 43,8% 43,0% 43,0%
		90% 45,0% 44,9% 44,1% 43,0% 42,0% 41,4% 41,4%
	100% 43,0% 42,9% 42,2% 41,5% 40,8% 40,4% 40,4%
	110% 41,5% 41,4% 40,8% 40,3% 39,9% 39,6% 39,6%
	120% 41,0% 40,9% 40,6% 40,2% 39,7% 39,4% 39,4%
	Using a Monte-KT	0,088		0,167	0,25		0,50	0,75	1,0	1,1
	80% 22,90% 22,54% 21,80% 20,72% 20,00% 19,61% 19,54%
	90% 21,52% 21,21% 20,73% 19,78% 19,34% 19,03% 19,04%
	100% 20,01% 19,94% 19,83% 19,15% 18,91% 18,68% 18,71%
	110% 18,73% 18,43% 18,77% 18,58% 18,47% 18,26% 18,31%
	120% 17,91% 17,67% 18,04% 18,11% 18,15% 18,04% 18,13%

  5.

	KT	0,25		0,50		0,75		1,0
		Theor.	LC	Theor.	LC	Theor.	LC	Theor.	LC
	80%	40.78% 40.35% 39.16% 38.77% 37.87% 37.49% 37.16% 36.84%
	90%	39.73% 39.51% 38.39% 38.10% 37.26% 36.97% 36.74% 36.46%
	100% 38.44% 38.37% 37.52% 37.27% 36.71% 36.40% 36.32% 36.01%
	110% 37.54% 37.40% 37.06% 36.78% 36.20% 35.85% 36.05% 35.71%
	120% 37.31% 37.24% 36.89% 36.68% 36.10% 35.81% 35.98% 35.68%

  a simple Taylor expansion with remainder gives (x i + θ 1 , y j ) + ∂ 4 p n ∂x 4 (x i -θ 2 , y j )) (x i , y j + θ 1 ) + ∂ 4 p n ∂x∂y 3 (x i , y j -θ 2 ))

	As for the space derivatives, clearly
		δ 2 x p n = ∂x 4 δ 2 ∂ 2 p n ∂x 2 (x i , y j ) + ∆x 2 24 ∂ 4 p n ( xy p n = ∂ 2 p n ∂x∂y (x i , y j ) + ∆y 2 12 ∂ 4 p n ( ∂x∂y 3 + ∆x 2 24∆y [ ∂ 3 p n ∂x 3 (x i + θ 3 , y j-1 ) -
	p n+1 i,j -p n i,j ∆t	=	∂p ∂t (t n , x i , y j ) + ∂p ∂t (t n+1 , x i , y j ) 2	+	∆t 2	(	∂ 2 p ∂t 2 (t n + θ, x i , y j ) -	∂ 2 p ∂t 2 (t n + θ * , x i , y j ))

  ) at times t n and t n+1 gives us E = e 1 -e 21 -e 22 -e 23 -e 3 + e 41 + e 42 + e 5 with e 1 the error coming from the time-derivative (t n + θ * , x i , y j )) e 21 , e 22 and e 23 come from the second order space-derivatives

	e 1 = ∂t 2 e 21 = ∆t 2 ( ∂ 2 p ∂ 2 p ∂t 2 (t n + θ, x i , y j ) -∆x 2 96 [ ∂ 4 ∂x 4 (I 2

Eventually, a Crank-Nicholson like formula appears

Let us study the consistency of this new scheme Proposition 5. The algorithm with a corrector step is also consistent. The truncation error is

Proof. To prove the consistency, let q n+1 be defined as

))

The computations are almost identical to the previous proposition. This time the error is equal to

Using the same decomposition, the errors e


Xiaofei Lu
email: xiaofei.lu@centralesupelec.fr

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

High dimensional Hawkes processes for limit order books Modelling, empirical analysis and numerical calibration

Keywords: Hawkes processes, limit order books, high frequency data, non-convex optimization JEL Classification: C58, C61, G14

High-dimensional Hawkes processes with exponential kernels are used to describe limit order books in order-driven financial markets. The dependencies between orders of various types are carefully studied and modelled, based on a thorough empirical analysis. The observation of inhibition effects is particularly interesting, and leads us to the use of non-linear Hawkes processes. A specific attention is devoted to the calibration problem, in order to account for the high dimensionality of the problem and the very poor convexity properties of the MLE. Our analyses show a good agreement between the statistical properties of order book data and those of the model.

Introduction

Limit order books have attracted a considerable amount of attention since the electronification of financial markets in the early '90s. The historical quote-driven markets, where designated market makers used to provide liquidity to all participants, have largely evolved into order-driven markets, where buy and sell orders are matched continuously in a double auction queueing system.

In an order-driven market, participants can submit orders of three basic types: limit order, market order and cancellation:

• Limit order: An order that specifies an upper/lower price limit (also called "quote") at which one (commonly called "liquidity provider") is willing to buy/sell a certain number of shares.

The advantage of the limit order is that the transaction price is better than the instantaneous mid-price. However, there is no certainty that the limit order will be executed. Currently most markets adopt the "first in first out" rule, i.e. the priorities of limit orders are decided first according to price, and then to arrival time. A limit order can be entirely, partly or not executed. • Market order: An order that triggers an immediate buy/sell transaction for a certain number of shares at the best available opposite quote(s). The advantage is to offer an immediate execution, however the price is worse than the mid-price. A market order can be executed with different limit orders as counterparties. The price is not necessarily the best limit price, if the quantity is big enough that the order eats up completely the first limit and hits the second or higher limits. • Cancellation : An order that removes an existing limit order.

In addition to these three main types of orders, there exist various order services provided by the exchanges such as "stop loss", "good til' canceled"... Also note that some markets allow orders, such as "iceberg" orders, to provide hidden liquidity, making their presence difficult to infer from the order flow. Nevertheless, it is commonly agreed upon -and verified in practice -that the basic orders carry enough information for market microstructure studies.

An example is given in Figure 1.

V 4 b V 3 b V 2 b V 1 b V 1 a V 2 a V 3 a V 4 a V 5
a Figure 1. Illustrative order book example. Blue bars on the left half of the figure represent the available buy orders with prices P • b and total quantities V • b . These correspond to the buyer side, also called the bid side. Participants in the bid side are providing liquidity with prices at which they are ready to buy some quantities of the stock. The right hand bars represent the sell side, commonly called the ask side or offer side, where participants willing to sell post their orders with the prices they are ready to sell the stock. The line in the middle corresponds to the mid-price level and is computed as the average between the best (highest) bid price and the best (lowest) ask price. A transaction occurs when a sell order and a buy order are at least partially matched. A queue of limit orders with the same price is called a limit. Different colors in the same limit represent orders with different priority with darker bars having higher execution priority.

Due to their obvious relevance to the understanding of financial markets, limit order books have been extensively studied in the past decade, both from empirical and mathematical points of view, see [START_REF] Abergel | Limit order books[END_REF] for a survey of their properties. In particular, the mathematical modelling of limit order books is itself an active research area that has many useful and practical applications, and this paper is a contribution to the field.

A popular class of order book models is that of Markovian models, originating with the so-called zero-intelligence models as in [START_REF] Smith | Statistical theory of the continuous double auction[END_REF], then enriched with more complex and realistic contributions such as [START_REF] Cont | A stochastic model for order book dynamics[END_REF] or [START_REF] Huang | Simulating and analyzing order book data: The queue-reactive model[END_REF]. In Markovian models the order flows are described by point processes with state-dependent conditional intensities.

In parallel with the development of Markovian models, and serving as an inspiration for the present work, many empirical studies have identified some time dependency properties of financial markets. To name a few, [START_REF] Gopikrishnan | Statistical properties of share volume traded in financial markets[END_REF], [START_REF] Bouchaud | How Markets Slowly Digest Changes in supply and Demand[END_REF] underline a significant positive autocorrelation and slow decay of the trade flow. [START_REF] Chakraborti | Econophysics review: I. Empirical facts[END_REF] confirms that the Poisson hypothesis for the arrival of orders is not empirically satisfied, whereas [START_REF] Eisler | The price impact of order book events: market orders, limit orders and cancellations[END_REF] is an in-depth study of the correlation between, and price impact of, orders of all types. These findings advocate for a direct modelling of the temporal dependencies between order arrivals and in such a context, Hawkes processes have come up as a very adequate choice.

For this reason, Hawkes processes have triggered a lot of interest in recent studies on market microstructure and limit order books. To name a few, [START_REF] Bacry | Hawkes model for price and trades high-frequency dynamics[END_REF], [START_REF] Bacry | Modelling microstructure noise with mutually exciting point processes[END_REF][START_REF] Bacry | Estimation of slowly decreasing Hawkes kernels: application to highfrequency order book dynamics[END_REF] propose various models of price and order flow models, whereas [START_REF] Bacry | Hawkes processes in finance[END_REF] is an extensive survey of the application of Hawkes process in finance. In the specific context of limit order books, [START_REF] Large | Measuring the resiliency of an electronic limit order book[END_REF] is an early study of Hawkes processes applied to order book modelling, Hawkes-process-based limit order book models are introduced and mathematically investigated in [START_REF] Abergel | Long-Time Behavior of a Hawkes Process-Based Limit Order Book[END_REF], [START_REF] Zheng | Modelling Bid and Ask Prices Using Constrained Hawkes Processes: Ergodicity and Scaling Limit[END_REF] and, in a sligthly different direction, [START_REF] Rambaldi | The role of volume in order book dynamics: a multivariate Hawkes process analysis[END_REF] models the order volumes -in addition to their types -based on a multivariate Hawkes process.

In this paper, we contribute to this very fruitful and original strand of research by designing, analyzing and calibrating high-dimensional Hawkes processes describing limit order books. The effects that are most important to model are the mutual excitations between orders of various types, including some inhibition effects that are hightlighted in our empirical analysis.

The quality of various Hawkes-process-based order book models will be assessed using some objective criteria: a model will be deemed satisfactory if it can reproduce as many as possible of the stylized facts of financial data. Our approach starts with a precise empirical analysis of the dependencies between order arrivals of various types. Then, models built from multivariate, possibly nonlinear, Hawkes processes with multiple exponential kernels are introduced. Once a model is designed, it is evaluated. Following most papers involving Hawkes processes for order book modelling, the natural quantities of interest are the inter-event durations -or : inter-event forward recurrence times, which will be the main objects under scrutiny in the present work. The distribution of forward recurrence times, as well as the signature plot1 , are used as selection criteria. With this approach, we are able to discriminate between various Hawkes-process-based models, and provide a financial interpretation of the more successful ones in terms of their behaviour at various time scales, and the presence of inhibition as well as excitation effects.

In a slightly different, yet very important perspective, we also address carefully the estimation problem: it is well-known that the Maximum Likelihood Estimator of a Hawkes process is not the solution to a convex optimization problem, and therefore, great care must be taken when using it to obtain the values of the parameters. Actually, we found this part of our work so relevant and influential in practice that we decided it was worth a full section of the paper, hopefully sharing some useful knowledge with other potential users of high dimensional Hawkes processes.

The paper is organized as follows: Section 2 presents our most relevant empirical findings, laying the ground for the modelling based on linear and nonlinear Hawkes processes discussed in Section 3. In Section 4, the numerical aspects of model calibration are discussed in detail. Finally, some concluding remarks are presented in Section 5.

Empirical findings: the interdependencies of order book events

In this section, we present our main empirical findings on the dependencies between order arrivals. These findings pave the way for the modelling avenues followed in the next sections.

Data and Framework

This paper focuses on the DAX listed 30 stocks trading in XETRA -the electronic trading venue of the Frankfurt Stock Exchange. Three months (February to April 2016) of tick-by-tick data are used in this study. The data consist in the list of all trades and order book states any time a modification or a transaction occurs -with a resolution of 1µs (10 -6 s). As is classical for high frequency financial data, see e.g. Muni Toke (2017) for a recent survey on order book reconstruction,

Notation Definition

M , L, C, O market order, limit order, cancellation, any order.

M 0 , L 0 , C 0 , O 0 market order, limit order, cancellation, any order, that does not change the mid-price.

M1 , L 1 , C 1 , O 1 market order, limit order, cancellation, any order, that changes the mid-price.

M buy , M sell buy/sell market order.

M 0 buy , M 0 sell buy/sell market order that does not change the mid-price: i.e. order quantity < best ask/bid available quantity.

M 1 buy , M 1 sell buy/sell market order that changes the mid-price: ie. order quantity ≥ best ask/bid available quantity.

L buy , L sell buy/sell limit order.

L 0 buy , L 0 sell buy/sell limit order that does not change the mid-price: i.e. order price ≤ / ≥ best bid/ask price.

L 1 buy , L 1 sell buy/sell limit order that changes the mid-price: ie. order price > / < best bid/ask price.

C buy , C sell buy/sell cancellation.

C 0 buy , C 0 sell buy/sell cancellation that does not change the mid-price: i.e. partial cancellation at best bid/ask limit or cancellation at another limit.

C 1 buy , C 1 sell buy/sell cancellation that changes the mid-price: ie. total cancellation of best bid/ask limit order. some data cleaning was involved in order to identify limit orders, market orders and cancellations given the states of the order book and the list of trades. Due to the large quantity of data, problems such as mismatches of quantities and lack of synchronization were expected. However, such anomalies represent less than 3% of the data, and our results are thus reliable.

Event definitions

In this study, any change that modifies the best limits of the order book is called an "event" 1 . More precisely, an event can be a limit order, a market order, or a cancellation, and can affect the best bid or best ask. Moreover, events will be tagged according to whether they change the mid-price or not. Table 1 summarizes the definitions and notations for the various event types considered in this paper.

Statistical dependencies between order book events

Table 2 represents the empirical probabilities of occurrence of an event of type j (in column), conditioned on the fact that the last observed event is of type i (in row). The last row represents the unconditional probabilities of each type of events.

To simplify the interpretation of the results, Table 3 represents the ratio of conditional probabilities to unconditional probabilities, rounded to one decimal. It aims at revealing the mutual

L 0 buy L 0 sell C 0 buy C 0 sell M 0 buy M 0 sell L 1 buy L 1 sell C 1 buy C 1 sell M 1 buy M 1 sell L 0
L 0 buy L 0 sell C 0 buy C 0 sell M 0 buy M 0 sell L 1 buy L 1 sell C 1 buy C 1 sell M 1 buy M 1 sell L 0 buy 1.4 0.5 1.3 1.0 1.4 0.7 1.2 0.4 0.0 1.2 1.6 0.3 L 0 sell 0.5 1.4 1.0 1.3 0.7 1.4 0.5 1.2 1.2 0.0 0.3 1.6 C 0 buy 1.1 1.2 1.4 0.4 0.8 0.8 0.7 1.1 1.7 0.4 0.6 0.7 C 0 sell 1.2 1.1 0.4 1.4 0.8 0.8 1.1 0.7 0.5 1.7 0.7 0.6 M 0 buy 1.0 0.4 0.3 0.5 11.5 0.8 2.4 0.2 0.4 2.3 15.0 0.6 M 0 sell 0.4 1.0 0.5 0.3 0.8 13.2 0.2 2.4 2.3 0.4 0.6 15.7 L 1 buy 1.6 0.5 0.1 1.3 1.2 1.3 1.1 0.4 3.0 1.7 1.5 2.9 L 1 sell 0.5 1.6 1.3 0.1 1.2 1.2 0.4 1.1 1.7 3.0 2.8 1.5 C 1 buy 0.7 0.6 2.5 0.2 0.3 0.1 2.1 1.0 1.0 0.3 0.3 0.1 C 1 sell 0.6 0.7 0.2 2.5 0.1 0.3 1.0 2.1 0.3 0.9 0.1 0.3 M 1 buy 0.6 0.3 0.2 1.5 1.3 0.9 6.2 1.1 0.4 2.9 2.1 1.6 M 1 sell 0.3 0.6 1.5 0.2 0.9 1.5 1.1 6.2 3.0 0.4 1.6 2.1
relationships between events, and ratios greater than two are highlighted.

Results of Table 3 are quite symmetric and no significant differences are observed between the buy and the sell side. Therefore, only the buy side is interpreted in detail below:

• L 0 buy : adds liquidity to the first limit, signalling an increase of market demand at the current price level. This stimulates L 1 buy and M 1 buy events, based on the new consensus for a higher price. The corresponding probabilities for orders of type 0 are also increased, based on a similar reasoning but in a less aggressive way. On the other hand, the selling activity decreases in general except for C 1 sell , because some newly added limit orders may be cancelled shortly after. One notable thing is the sharp decrease of C 1 buy , as the newly added limit order probably comes from another trader, making it very unlikely that the first limit should be cancelled.

• C 0 buy : decreases liquidity on the buy side. It triggers successive cancellations C 0 buy and C 1 buy : cancellations tend to follow themselves. M 0 buy , M 0 sell , M 1 buy and M 1 sell become less likely, revealing the influence of low liquidity on the participants' willingness to generate executions.

• M 0 buy : largely increases the probability of M 0 buy and M1 buy . This is commonly attributed to order splitting and the momentum effect (other participants following the move). L 1 buy and C 1 sell are also stimulated as a new price consensus emerges. • L 1 buy : improves the offered price to buy. The first effect is a strong increase in the probability of M 1 sell , i.e., participants entirely consume the new liquidity as the offered price has become higher. The second effect is an increase in the probability of C 1 buy , i.e., the new liquidity is rapidly cancelled. This is consistent with a similar observation made for L 0 buy orders, and might reflect some sort of market manipulation where agents are posting fake orders. Not surprisingly, the conditional probability of C 0 buy is almost zero, because after one limit order is submitted, it cannot be partially cancelled. The fact that probability is not exactly 0 may be due to poor data synchronization, or the existence of hidden liquidity.

• C 1 buy : a total cancellation of the best buy limit increases the probability of C 0 buy -order cancellations come in succession as market makers lose interest to provide liquidity even at the new best limit -and that of L 1 buy events, as traders may re-offer at the previous best price to gain priority. Events of other types become less frequent.

• M 1 buy : consumes all the offered liquidity at the best ask. It stimulates L 1 buy and C 1 sell as a higher price consensus emerges among market participants. The probability of M 1 buy increases, indicating a short term momentum effect, and order splitting.

As a conclusion to this empirical section, let us just say that strong temporal dependencies between events are identified. Some orders actually trigger other events, a fact that can be seen as self-or cross-excitation phenomena. There are also some inhibition effects, when incoming orders prevent other events to occur. These two important features will be the target of the modelling approach presented in the next section.

Modelling dependencies using Hawkes processes

It is now widely accepted in the high frequency and market microstructure community that limit order books are worth modelling, and that the price dynamics can easily be extracted from that of the order book. In fact, the complexity of inter-event dependencies is so high that most significant features of the price dynamics : co-existence of time scales, leverage effect, signature plot, long term diffusivity... can be derived from advanced order book models.

In this section, point-process-based order book models are studied, building on the 12 event types previously introduced:

E = {L 0 buy , L 0 sell , C 0 buy , C 0 sell , M 0 buy , M 0 sell , L 1 buy , L 1 sell , C 1 buy , C 1 sell , M 1
buy , M 1 sell }, where events with superscript 1 have an instantaneous price impact. In particular, it is clear that events in E up = {L 1 buy , C 1 sell , M 1 buy } lead to a price increase, whereas those in

E down = {L 1 sell , C 1 buy , M 1 sell } result in price decrease.
The arrival of order book events is therefore modelled by a 12-variate simple point process

N (t) = (N L 0 buy (t), . . . , N M 1 sell (t)).
Of interest is the associated intensity process (λ L 0 buy (t), . . . , λ M 1 sell (t)). Assuming that the process is simple means that two events cannot occur at the same time, a fairly realistic assumption due to the high time resolution of modern stock exchanges.

Since the focus in this paper is on temporal interdependencies, N is actually a 12-variate counting process, and the marks determining the price jump when an event of type 1 occurs are not modelled. Rather, a simplifying assumption is made, namely, that the jump of the best bid or ask price following an event of type 1 is always one tick. This approximation reduces the dimensionality of the point process, while being consistent with the real behaviour of the chosen data set, for which the average jump size of the best bid and ask prices is 1.08 ticks 1 . Under this assumption, the reconstructed mid-price dynamics easily obtains as a by-product of event arrivals:

S(t) = S(0) +   e∈Eup N e (t) - e ∈Edown N e (t)   × η 2 , t > 0
where η > 0 is the tick size. This simplification will be taken into account when comparing the performances of the model with the behaviour of real data.

Clearly, events of type 0 do not directly influence the price, rather, their impact will come from their influence on the intensities of the type 1 event arrival process.

As already said in the introduction, it has long been recognized that the class of Hawkes processes is particularly well suited to the modelling of point processes interacting via their conditional intensities. In the context of limit order book and high-frequency financial data, the paper Bacry et al. ( 2016) is a very interesting and relevant application of Hawkes processes to level-one limit order book modelling.

Let us now build on the results of Section 2 and study two classes of models, based respectively on linear and nonlinear Hawkes processes, that capture well the main characteristics of market dynamics. The performances of the models are presented in this section, while some more technical aspects pertaining to their numerical calibration are deferred until Section 4.

Linear Hawkes process models

In this short paragraph, we recall some essential definitions and results on the particularly interesting class of point processes introduced in [START_REF] Hawkes | A Cluster Process Representation of a Self-Exciting Process[END_REF]. We refer the interested readers to [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF], [START_REF] Massoulié | Stability results for a general class of interacting point processes dynamics, and applications[END_REF] for a more in-depth presentation of these processes, and to [START_REF] Zheng | Modelling Bid and Ask Prices Using Constrained Hawkes Processes: Ergodicity and Scaling Limit[END_REF], [START_REF] Abergel | Long-Time Behavior of a Hawkes Process-Based Limit Order Book[END_REF] [START_REF] Bacry | Estimation of slowly decreasing Hawkes kernels: application to highfrequency order book dynamics[END_REF] for recent applications to order book modelling.

A multivariate point process ((T i , X i )) i∈N * , associated to a counting process (N (t)) t∈R+ = (N 1 (t), . . . , N M (t)) t∈R+ with conditional intensity process (λ(t)) t∈R+ = (λ 1 (t), . . . , λ M (t)) t∈R+ , is called a (linear, multivariate) Hawkes process [START_REF] Hawkes | A Cluster Process Representation of a Self-Exciting Process[END_REF], [START_REF] Massoulié | Stability results for a general class of interacting point processes dynamics, and applications[END_REF] if there holds for m ∈ {1, . . . , M } :

λ m (t) = µ m + M n=1 t 0 φ mn (t -s)dN n (s)
where µ m are positive real numbers and φ mn are nonnegative functions.

The µ m are the base intensities and can be viewed as background intensities. Whenever an event occurs, the intensities increase, making subsequent events arrive at a higher frequency. Such effects are controlled by φ mn . The functions φ mn , the kernel functions, control the instantaneous increases and the relaxation speeds of the intensities in response to excitations.

For a multivariate Hawkes process, φ mm describe the self-excitations, while φ mn for m = n measure the cross-(or: mutual) excitations, that is, the impact of an event of type n on the arrival of an event of type m.

A convenient, alternate way to express the intensity process is provided by the following equation:

λ(t) = µ + Φ dN (1)
where Φ(t) is the M × M matrix whose entries are φ mn (t), " " denotes the "matrix convolution"

Φ dN = R φ(t -s)dN (s)
and φ(t -s)dN (s) stands for the standard matrix-vector product.

It is clear that Hawkes processes are fully determined by their baseline intensity µ and the matrix Φ of kernel functions. In the following, we will concentrate on exponential kernels. This particular choice is classical, one of its main advantages being the Markovianity of the joint process (N, λ), see e.g. [START_REF] Massoulié | Stability results for a general class of interacting point processes dynamics, and applications[END_REF]. For the models considered in this work, the intensities follow Equation (1), where Φ is a 12 × 12 kernel function matrix describing the excitation between events of various types:

Φ = (φ ij ) i,j∈E .
What we call 1-exponential and 2-exponential Hawkes models simply differ in the number of exponential functions used to define each kernel, namely:

• For the 1-exponential Hawkes model, φ ij (t) = α ij exp(-β ij t) • For the 2-exponential Hawkes model, φ ij (t) = 2 p=1 α ijp exp(-β ijp t)

Performances of the linear Hawkes models

The adequacy of a linear Hawkes-process-based order book model is now evaluated, according to two criteria: a goodness-of-fit criterion for the distribution of forward recurrence times, and a criterion based on the signature plot generated by the model. As a matter of fact, it is generally agreed upon that such statistical properties of the price process as the unconditional distribution of returns or the diffusive behaviour at large time scale, can easily be reproduced even with simpler models, whereas the signature plot and the inter-event durations offer a better challenge to discriminate among order book models.

3.2.1.

Goodness of fit. It is well-known, see e.g. [START_REF] Bowsher | Modelling security market events in continuous time: Intensity based, multivariate point process models[END_REF] that the transformed durations {τ } i of a Hawkes process

τ m i = Ti+1
Ti λ m (s)ds are i.i.d. exponential random variables with parameter 1. This property is used to test the goodnessof-fit of the model to the data, by drawing Q-Q plots of the empirical quantiles with respect to the theoretical exponential distribution quantiles. Though a global test can be conducted by concatenating all the transformed durations, plotting each dimension separately provides more information. This can be viewed as a marginal distribution fit test, i.e.: Given the law of other types of orders, how well can we fit the order under scrutiny ?

The procedure is as follows: first, the parameters for several order book models (Poisson, 1exponential linear Hawkes, 2-exponential linear Hawkes) are calibrated, for each day in the study period. Then the transformed durations in the model are computed, and a Q-Q plot test is then performed. The results are shown in Figure 2.

As a first conclusion, one can easily see that a Poisson-process-based model globally fails to capture the distributional properties of recurrence times. The performances of the 1-and 2-exponential Hawkes models are similar, except for orders of type 0: the 2-exponential model significantly outperforms the 1-exponential model for L 0 and, to a lesser extent, for C 0 events. However, what is annoying is the behaviour for C 1 events: the distributions of the transformed durations in 1-and 2-exponential models are extremely close to one another, but neither is close to the theoretical exponential distribution. This is an important, negative feature of the linear Hawkes models that will be revisited in the upcoming Subsection 3.3

3.2.2.

Signature plots. The signature plot reveals some of the most important stylized facts about high frequency financial data. It is a plot of the realized variance as a function of the sampling frequencies. The realized variance for a stochastic process X t over a time period [0, T ] at a sampling frequency h is simply

RV (h) = 1 T T /h n=0 (X((n + 1)h) -X(nh)) 2 . ( 2
)
An important stylized fact of financial markets is that the quantity RV generally increases when h becomes small. This phenomenon is associated to the mean reverting behaviour of the price at short time scales. It has long been observed and was already reported in Andersen et al. [START_REF] Andersen | Great realizations[END_REF]. It is noteworthy that the signature plot becomes even steeper when computed on transaction prices rather than mid-prices because of the bid-ask bounce, and we will focus on mid-prices to avoid this spurrious effect.

Once the model parameters are calibrated, the mid-price is easily simulated using Equation (3). Realized variances are calculated with sampling periods from 1 to 50 seconds, with a step of 1 second.

The results for the models and the real data are shown in Figure 3. Not surprisingly, the signature plot of the Poisson model is flat -this is expected, as the price dynamics in this model is that of a mid-price model with Poisson jumps, due to the mapping of orders that increase (resp. decrease) the price into upward (resp.downward) jumps.

The 1-exponential and 2-exponential Hawkes process models behave similarly: the realized volatility decreases when the sampling interval increases, but the long-term volatility level is too high compared to the data. Though reproducing the overall shape of the signature plot, the linear Hawkes-process-based order book models are not satisfactory.

Nonlinear Hawkes process model

This subsection addresses the shortcomings of linear Hawkes models in reproducing some characteristics of forward recurrence times and signature plots. Nonlinear Hawkes processes are introduced to overcome these difficulties, and their performances are studied.

3.3.1.

Order dependencies: inconsistencies between real data and linear Hawkes models. The results presented in Paragraph 3.2.1 are now revisited in event time, temporarily ignoring the durations. When comparing the average conditional probability matrix of the 2-exponential Hawkes model with that of real data, one can check that most of the conditional probabilities are pretty close. However, for several pairs, there exist huge differences between the model and the real data, in particular for C 1 buy |L 0 buy , M 0 buy |C 1 buy and M 0 sell |C 1 buy . Table 3.3.1 below gives the list of all pairs (X, Y ) for which the probability of an event of type X, conditioned on following an event of type Y , in the simulated order flow is either smaller than 50% or greater than 5 times the real conditional probability (only the buy side is shown, the sell side behaves similarly).

From a financial point of view, these discrepancies can easily be accounted for:

• A C 1 buy |L 0 buy sequence almost never happens, because L 0 buy is a limit order added to the current first limit and it is highly unlikely that two orders should be cancelled at the same microsecond.

• The low probabilities of L 1 buy |L 1 buy and L 1 sell |L 1 buy comes from the constraint of the bid-ask spread: an aggressive limit order decreases the spread, and when the spread becomes one tick wide, other price-changing limit orders are no longer possible. • The remaining cases correspond to orders following a C 1 buy order that increases the bid-ask spread. There is no physical constraint preventing the spread from being wide, but participants in the market are not seemingly ready to sell when a cancellation order has already decreased the best bid price.

L 0 buy L 0 sell C 0 buy C 0 sell M 0 buy M 0 sell L 1 buy L 1 sell C 1 buy C 1 sell M 1 buy M 1 sell 0.
From a mathematical point of view, this poor fit comes from an inherent shortcoming of the linear Hawkes process model: the intensity for the arrival of an order of type e is written as

λ e (t) = µ e + e ∈E T e <t φ ee (t -T e )
where φ e ≥ 0 and, by construction, λ e (t) < µ e cannot happen ! Consequently, inhibition effects, leading to a temporary decrease of certain short term conditional probabilities, are not modeled.

Note that, when calibrating the linear model (see Section 4 for details), the kernels corresponding to inhibitory behaviours are indeed forced to 0.

Below are the median values of the L 1 norms of the kernels stemming from the calibration results for C 1 buy stimulations in Table 5: clearly, kernels corresponding to the event pairs listed in Table 3.3.1 have norms equal to 0.

Moreover, two other event pairs come out of the calibration with 0 kernel norms, M 0 buy |C 1 buy and C 1 sell |C 1 buy . Although less obvious from the conditional probability matrix, this phenomenon is easy to interpret: a defensive cancellation on the bid side indicates a consensus of a fair price decrease in the market, therefore traders are less willing to buy at the previous ask price or cancel an existing ask order as it has already gained some queue priority with a profitable price.

Model definition.

In order to incorporate inhibitory behaviours in the model, negative kernels are introduced in the Hawkes process. Then, a truncation is applied to avoid meaningless negative process intensities.

In the new model, the intensities satisfy the equation

λ(t) = (µ + Φ dN ) + , (3)
where the entries of the matrix Φ are no longer constrained to take on positive values, and () + denotes the elementwise positive part function.

When enriched with the nonlinearity, the 2-exponential Hawkes process model retains its Markovian nature, see e.g. [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF] [START_REF] Zhu | Large deviations for Markovian nonlinear Hawkes processes[END_REF] for general results on nonlinear October 29, 2017 Quantitative Finance

LuAbergelHawkes

Hawkes processes. The negative kernels are chosen under the following form

φ mn = 2 p=1 -α mnp exp(-β mnp t)
where the α's and β's are nonnegative real numbers. Note that we fix the same sign for the two exponentials, in order to avoid overfitting -it is actually unexpected for interdependencies to have different time regimes, for example an inhibitory effect in the short term that would become an excitation in the long run. .1593 -0.0541 -0.1439 -0.1908 Table 8. Median of optimal likelihood functions for each type of order in 2-exponential and non-linear models. statistically more satisfactory fit than the linear 2-exponential Hawkes model previously studied. This better performance will be confirmed by the analysis of the signature plots and forward recurrence times.

M 0 buy M 0 sell C 1 buy C 1 sell M 1 sell -0.0319 -0
Model L 0 buy L 0 sell C 0 buy C 0 sell M 0 buy M 0 sell 2-
In Tables 6 and7, the calibration results of the nonlinear model are compared to those of the 2-exponential linear model. Clearly, the conditional probabilities are closer to real data, showing that a definite improvement is achieved by the nonlinear model. Also note that the kernels that were formerly set to 0 now take on quite significant negative norms, a fact which confirms that the inhibition effect plays an important role in the order dynamics.

Another interesting insight is provided in Table 8, by comparing the optimal values of the loglikelihood functions for both models. One can actually see where inhibition effects become more pregnant: orders of type 1 are clearly more influenced than orders of type 0, confirming the improvement already observed for the Q-Q plots in the goodness-of-fit test.

3.4.2.

Signature plots. The signature plots of linear and nonlinear 2-exponential Hawkes models are shown in Figure 5, and compared to that of real data. The asymptotic volatility level significantly improves with the nonlinear model, and the resulting signature plot is overall a very good fit.

3.4.3.

Analysis of self-and cross-excitation recurrence times. The rationale behind the introduction of nonlinear Hawkes models was the empirically observed presence of inhibitory effects among events. As a consequence, one should hope that the inter-event recurrence times would behave in a more realistic way with these models.

Figure 6 and 7 show the cumulative distribution function (CDF) and the probability density of the (logarithm of) the forward recurrence times for all events of type 1 -that is, the forward recurrence times of (or: duration between) price jumps.

Specifically, define the inter-jump duration as

∆T i = T i+1 -T i
where T i are the timestamps of the event arrivals.

According to the type of event causing the jump, these durations are furthermore separated into two subgroups: self-excitation durations ∆T a ∈ {∆T i |X i = X i+1 } and cross-excitation durations

∆T c ∈ {∆T i |X i = X i+1 }.
Inter-jump durations predicted by the model are then computed, and compared to data: although the linear Hawkes model already performs well in reproducing the inter-jump duration distributions both for self-and cross-excitations, one can see that the nonlinear Hawkes process further improves the fit in the range between milliseconds and seconds (log 10 (∆T ) ∈ (-3, 1)).

As a conclusion, one can say that the nonlinear Hawkes model provides a very satisfactory enhancement to the classical one, whether one uses Q-Q plots, signature plots or inter-jump recurrence times as benchmarks. This improvement is in fact quite natural, and is related to the empirical evidence presented in 3.3.1 on inhibition effects between events.

Some numerical aspects of model calibration

This section is devoted to an analysis of the numerical algorithms used to calibrate the various models introduced in Section 3. Although rather technical, we think it is relevant -actually, very useful -for readers interested in calibrating high-dimensional Hawkes-processes to high frequency financial data (or other types of data).

Several optimization procedures are discussed and compared, and the best performer among those we have tested is thorougly investigated.

Calibration with maximum likelihood estimation

Let ((T i , X i )) i∈N * be a multivariate point process with associated counting process (N 1 (t), . . . , N M (t)), whose intensities are to be estimated.

The log-likelihood, see [START_REF] Ozaki | Maximum likelihood estimation of Hawkes' self-exciting point processes[END_REF][START_REF] Rubin | Regular point processes and their detection[END_REF], of given intensities (λ 1 (t), . . . , λ M (t)), and a sample of observation {T i , X i } i∈{1,...,M } , is defined by the sum of the log-likelihood of each component:

ln L(λ, {T i , X i } i∈{1,...,D} ) = m ln L m (λ m , {T i , X i } i≤D ) = M m=1 T 0 ln λ m (s)dN m (s) + T 0 (-λ m (s))ds .
In the case of a Hawkes process with exponential kernels, a straightforward computation gives:

T 0 ln λ m (s)dN m (s) = T m i ln µ m + M n=1 α mn A mn (i) and T 0 λ m (s)ds = µ m T - M n=1 T n k α mn β mn e -βmn(T -T n k ) -1 , where A mn (i) = T n k <T m i e -βmn(T m i -T n k
) can be computed iteratively as

A mn (i) = A mn (i -1)e -βmn(T m i -T m i-1 ) + T m i-1 ≤T n k <T m i e -βmn(T m i -T n k ) so that ln L m (λ m , {T i , X i } i≤D ) = -µ m T + M n=1 T n k α mn β mn e -βmn(T -T n k ) -1 + T m i ln µ m + M n=1 α mn A mn (i) .
It is however clear, and quite unfortunate, that the likelihood function is not strictly concave. For example, in the 1-dimensional case, its expression simplifies to ln L(λ, {T }) = -µT + In fact, not only is the likelihood function not concave, but it actually has several local maxima. An illustrative example is given in Figure 8 where we draw the contour plot of the partial likelihood function ln L 2 of a simulated 2-dimensional Hawkes process. The kernels are exponential functions with parameters specified in Equation (4). While µ 2 , α 21 and α 22 are kept fixed, the likelihood values are plotted as functions of β 21 and β 22 . The two axes are presented in logarithmic scale.

It is clear that there are at least two local minima in this example. The existence of several local maxima make gradient-type algorithms less relevant for the maximum likelihood procedure and a global optimization algorithm appears necessary. The Nelder-Mead simplex algorithm (NM ) has been widely used in previous works on the calibration of Hawkes processes; however we find it not stable enough when a good a priori guess is not available. For these reasons, the Differential Evolution algorithm (DE ) [START_REF] Storn | Differential evolution-a simple and efficient adaptive scheme for global optimization over continuous spaces[END_REF] has been chosen to perform the optimization. DE is an efficient genetic evolutionary algorithm that has been adopted in various engineering domains such as electrical power systems, artificial neural networks, operation research, image processing... Starting from a population of randomly generated points, the algorithm performs a mutation-crossover-selection procedure, where the population is updated to have better objective function values and a large tentative space is scanned.

A pseudocode is given in Appendix A.

Benchmarking the DE algorithm

Simulation-based numerical experiments are performed in order to compare the efficiency of the NM and DE algorithms. More specifically, we consider a 2-dimensional Hawkes process or use more points for deviation v i,g = x r1,g + F (x r2,g -x r3,g ) + F (x r4,g -x r5,g )

Combinations of these ideas are of course possible, which create vast candidate strategies. • Crossover. Apart from the idea of the binomial/uniform crossover, another method called exponential crossover is also considered. The trial vector u takes the value of the donner vector v for adjacent coordinates. The benefit is limited to special structures of problems where neighboring variables are linked but relatively independent of other variables. As a result the binomial crossover is more frequently used. • Adaptation of control parameter (F and CR) and strategy. It aims at adding learning performances to the offspring generation. Either the strategies are randomly chosen from fixed ensemble of strategies and parameters, which are designed to aid the algorithm to converge or explore larger space so that the combination can balance the two effects; or the mutation strategy is fixed, but the parameters can adapt to the evolution. • Population control. The most natural idea is the reduction of population as they approach to each other and concentrate in a small region. Such reduction can be pre-scheduled or dynamically controlled based on the computational budget. On the other hand, varied population (instead of monotonically decreasing) is also introduced as a choice to adapt to the evolution of the algorithm.

Other extensions actually go beyond the classical framework, for example using new initialization techniques, adding clustering technique for the sub-population topology, and so on. Hybridization opens another branch of research: on the one hand DE is combined with other heuristic methods to explore the advantages of exploration strategies, and on the other hand, local search methods are injected into the DE algorithm to boost convergence and precision.

In the interest of tractability, we choose to concentrate on the non-hybrid extensions. In [START_REF] Das | Recent advances in differential evolution-An updated survey[END_REF], the algorithm L-SHADE is reported to have the "best competitive performance among nonhybrid algorithms at the CEC 2014 competition on real parameter single-objective optimization". Compared to the classical algorithm, L-SHADE combines adaptation in every respect -mutation, parameter control and population control:

• Mutation use the current -to -pbest/1 strategy, where the new donner vectors are obtained by

v i,g = x i,g + F i (x pbest,g -x i,g ) + F i (x r1,g -x r2,g ),
where x pbest,g is randomly selected from the best pN members in generation g, where (p ∈ [0, 1]) . This strategy exhibits some greediness towards the current best points, but the existence of p leaves the flexibility for tradeoff between exploitation and exploration. • Parameter control In order to dynamically adapt the parameters F and CR, a record of past candidates is maintained. Two lists of size H, M CR and M F , are kept. For each generation, F i and CR i are drawn randomly with certain distributions depending on randomly chosen means from the lists:

F i = randc i (M F,ri , 0.1), CR i = randn i (M CR,ri , 0.1)1 {MCR,r i =N ull} ,
where randn follows a normal distribution and randc, a Cauchy distribution. For each generation, the kth element (k = g mod H) of the list is updated, according to CR i and F i that succeed to find ameliorated points. Such mechanism introduces learning characteristics for the F and CR selection, in order to overcome the stagnation problem. • External archive introduction To maintain diversity, a external archive is used so that parent vectors that are worse than the trial vectors are preserved in A. When generating donner vectors, x r2,g can be selected from P ∪ A. • Linear population size reduction The whole population N g decreases according to the allowed total number of generations.

N g+1 = round N min -N init G * g + N init ,
where N init is the classical initial population size, and N min is the smallest possible population size for a mutation strategy.

The L-SHADE is a combination of interesting ideas. Roughly stated, the current-to-pbest/1 mutation helps approach the best candidates in the population, accelerating the convergence of the algorithm; the parameter control aims at learning the trade-off between exploration and exploitation; the external archive is to help keep diversity of the population so that exploration is partly internalized by the exploitation of the abandoned history; and the population size reduction saves computational cost to allow larger initial populations.

4.3.2.

Calibrating high dimensional Hawkes order book models. Let us now turn towards the actual application of L-shade to the task at hand.

Starting from 100 different initial populations for each strategy with the same number of points and maximum generations, we plot the histograms of the final log-likelihood function values for one dimension of the 12-dimensional Hawkes model with real data in Figure 9, for different modifications of DE. The classical strategy, noted as "rand/1", serves as a reference for the suggested "current-to-pbest/1". The parameter adaptation is also combined with "rand/1" to provide better performances. We finally introduce a version with a refinement of the initial parameter intervals, noted as "better guess". The right subplot is a zoom of the one on the left, to further show the improvement due to "better guess".

Some comments are in order :

• In the "current-to-pbest/1" strategy, the closer p is to 0, the greedier the algorithm is, and the more probable it is that the optimization gets trapped at a local maximum. The closer p is to 1, the more the algorithm favors exploration. • The learning mechanism for F and CR in the adaptative version leads to some improvements.

Parameters are initialized according to the following distributions:

µ m ∼ U(0, 0.2N m T ), ρ mn ∼ (0, min( 0.2N m N n , 0.5)), β mn = u 1 1 {u0=0} + u 2 1 {u0=1} for u 0 ∼ B(1, 0.5), u 1 ∼ U(0, 1) u 2 ∼ U(0, 100)
derived from the physical interpretation of µ as the baseline intensity, of ρ as the integrated intensity of the influence from event arrival, and based on the relation

E[λ m,∞ ]T = µ m T + E[ n ρ mn N n ].
• Although different runs starting from different initial populations do not converge to the global maximum, some improvement may be gained from a "better guess" of the initial intervals.

Clearly, an increase of the population size plays a major role in boosting the convergence: a larger population prevents points from getting trapped around the same local maximum. On the other hand, it is useless to keep all the population as the algorithm approaches the end of its iterations, since points tend to form clusters. As a consequence, it makes sense to consider effective population reduction techniques and use the saved computational budget to cover a larger search space.

Building on the linear population reduction method inspired by the combination of DE with clustering algorithms in [START_REF] Li | A differential evolution algorithm with dynamic population partition and local restart[END_REF] and the use of pairwise Euclidean distance for dynamic population control in [START_REF] Yang | An improved adaptive differential evolution algorithm with population adaptation[END_REF], we propose an additional reduction mechanism which allows, not only to decrease the function evaluation times, but also to avoid convergence to local maxima.

The algorithm is said to have converged if each coordinate of all the points in the population has converged. The convergence conditions of the coordinates are σ(µ m ) < e r µ m or max µ m -min µ m < e a , σ(ρ mn ) < e r ρ mn or max ρ mn -min ρ mn < e a , σ(β mn ) < e r β mn or ρ mn < e a , where σ(.) and • are the standard deviation and the mean value respectively, and e r and e a are the relative and absolute error tolerance. At each generation, we eliminate points that are close to the current b best ones, using a criterion similar to the termination conditions for the population: suppose the points are sorted according to their objective function values by descending order. For a given point x i , if ∃j ∈ 1, b /{i} such that all the following conditions are satisfied: The decrease of population size saves some computational budget for the algorithm, which is very beneficial as the computation of the likelihood function is costly. The combination of these population reduction techniques allows to increase the initial population by a factor of 5 to 10 with no significant impact on the total computation time, and the convergence is largely improved.

As a conclusion, one can say that the improved version L-SHADE of the DE algorithm drastically enhances the performances of the calibration, but despite all these efforts, we are still left with an average failure rate of approximately 5%.

Conclusion

This paper is a study of Hawkes processes applied to high frequency limit order book data. Suitably designed nonlinear Hawkes processes that include inhibitory effects and a co-existence of time scales are shown to successfully model the dependencies between the arrival of order book events. Thanks to the particularly well-suited distinction between events that trigger, or do not trigger, an immediate change in the current price, the dynamics of the model fully reflect that of the price. Such a description helps cope with some shortcomings of order book models that were previsouly observed, particularly concerning the realized spot price volatility.

The paper also gives a detailed analysis on a very important, albeit technical, topic: the choice of the optimization algorithm for the Maximum Likelihood Estimator. The L-SHADE algorithm is a significant improvement over the classical Differential Evolution algorithm, thanks to better initializations and population control.

As a conclusion, one can say that nonlinear Hawkes processes capture well such fundamental features of market dynamics as conditional probabilities, forward recurrence times, or the signature plot. They provide an accurate description of the order book in the high frequency realm, as well as a realistic behaviour of more macroscopic quantities. While leading to a better understanding of the mechanisms driving the markets, their use in the simulation of order driven markets can also lead to a host of potential applications.

Figure 2 .

 2 Figure 2. Q-Q plot goodness of fit tests of order book models.

Figure 3 .

 3 Figure 3. Mean signature plots of simulated price compared with real data.

3. 4 .

 4 Performances of the nonlinear Hawkes models 3.4.1. Goodness of fit. Similarly to the analysis presented in Paragraph 3.2.1, the Q-Q plots of the empirical quantiles with respect to those of the theoretical exponential distribution are shown on Figure 4.

Figure 4 .

 4 Figure 4. Q-Q plot goodness of fit tests of nonlinear Hawkes model.

Figure 5 .

 5 Figure 5. Mean signature plots of linear and nonlinear 2-exponential Hawkes process compared with real data.

Figure 6 .

 6 Figure 6. Cumulative distribution functions of log inter-jump durations for simulated price processes compared with real data.

Figure 7 .

 7 Figure 7. Probability density of log inter-jump durations.

Figure 8 .

 8 Figure 8. Example of local maxima in 2-d Hawkes process likelihood function ln L2(λ2)

Figure 9 .

 9 Figure 9. Distribution of optimal objective likelihood functions in different optimization strategy tests. The right one is zoomed at the optimal zone for further illustration.

  |µ mi -µ mj | < e r µ mj or |µ mi -µ mj | < e a |ρ mni -ρ mnj | < e r ρ mnj or |ρ mni -ρ mnj | < e a |β mni -β mnj | < e r β mnj or ρ mnj < e a ,then x i is eliminated from the population. In practice, it is convenient to select a small value for b.

Table 1 .

 1 Event types definitions.

Table 2 .

 2 Conditional probabilities of occurrences per event type.

Table 3 .

 3 Conditional probability leverage.

Table 4 .

 4 Conditional probability comparison between simulated order flows and real data.

	Pair	P simu	P real
	C 1 buy |L 0 buy L 1 buy |L 1 buy L 1 sell |L 1 buy M 0 sell |C 1 buy C 1 buy |C 1 buy M 1 sell |C 1 buy	0.402 1.628 1.288 0.545 0.548 0.854	0.048 0.141 0.171 0.068 0.072 0.037

Table 5 .

 5 Medians of L1 norm of kernels φ .C 1

	buy	in the 2-exponential model.

Table 6 .

 6 Improvement of conditional probabilities with the nonlinear Hawkes model.

	Pair	P simu	P N L simu	a	P real
	C 1 buy |L 0 buy L 1 buy |L 1 buy L 1 sell |L 1 buy M 0 sell |C 1 buy C 1 buy |C 1 buy M 1 sell |C 1 buy	0.402 1.628 1.288 0.545 0.548 0.854	0.105 0.254 0.235 0.135 0.163 0.108		0.048 0.141 0.171 0.068 0.072 0.037

a Conditional probabilities in nonlinear model.

Table 7 .

 7 Medians of L1 norm of kernels φ .C 1

	buy	in nonlinear model for
	those that were 0 in 2-exponential model.

A characterization of the realized price volatility at various frequencies.

This simplifying choice essentially means that a level-1 order book is considered.

Actually, for some large tick stocks, the average is even smaller than 1.01.

October 29, 2017 Quantitative Finance

where the parameters are specified in (4). 100 process paths are simulated for each T ∈ {100, 250, 500, 1000, 2500, 5000, 10000, 25000}, and the parameters are calibrated from each simulated path with various algorithms.

NM is used with different initialization methods. For NM random, the initial reference point is drawn from uniform distributions. Denoting by ρ the L 1 norm of the kernel (ρ = α β ), we choose

and optimize with respect to ρ instead of α.

The algorithm NM perfect refers to NM where the true input parameters are used as reference point.

The empirical probability of error for each optimization algorithm is show in Table 9 for T ∈ {250, 2500, 25000}:

Clearly, with the possible exception of short time horizon, DE almost always finds the optimal point, getting very close to the NM perfect algorithm.

Improvement in high dimensions

The local maximum problem is more severe when dealing with higher dimension and real data instead of simulated data. In this section, we present some treatments designed to mitigate the numerical issues and boost the convergence towards a global maximum.

4.3.1.

Some evolutions of the DE algorithm: a quick guided tour. Thanks to its wide variety of applications, DE has attracted a lot of interest, and the recent survey paper [START_REF] Das | Recent advances in differential evolution-An updated survey[END_REF] documents a host of novel ideas to improve its classical form. Below is a brief summary of some of the proposed improvements (notations are those used in Algorithm 1):

• Mutation strategy. The donner vector v i,g in mutation can be generated with different strategies. The classical algorithm adopts a so-called "DE/rand/1" strategy

where r g 1 , r g 2 and r g 3 are mutually exclusive integers randomly chosen in 1, N {i}. It could be preferable to approach the current best value v i,g = x best,g + F (x r2,g -x r3,g ) Appendix A: Pseudocode of basic Differential Evolution Algorithm 1 Differential Evolution algorithm 1: Input. Maximum total generation G, population size N ≥ 4, mutation factor F ∈ (0, 2), crossover rate CR ∈ (0, 1), parameter domain Ω, termination criteria. 2: Output. optimal point (optimal function value, termination generation etc.). 3: // Initialization phase 4: g=1; Initialize the initial population (x 1,1 , . . . , x N,1 ) randomly such that x i,1 ∈ Ω; 5: while g ≤ G and termination criteria not met do 6:

Choose randomly r 1 , r 2 and r 3 in 1, N such that i, r 1 , r 2 and r 3 are distinct;

9:

Construct donner v i,g+1 ← x r1,g + F (x r2,g -x r3,g ); if f (u i,g+1 ) ≤ f (x i,g ) then 23:

x i,g+1 ← u i,g+1 ; 24: else 25:

x i,g+1 ← x i,g ;

26:

end if 27: end for 28: g ← g + 1; 29: end while