Xiaofei Lu
email: xiaofei.lu@centralesupelec.fr

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

Limit order book modelling with high dimensional Hawkes processes

come

Limit order book modelling with high dimensional

Hawkes processes

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.

Limit order books have been extensively studied, both from empirical and mathematical points of view, see e.g. [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 particularly popular class of order book models is that of Markovian models, starting 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.

More to the point, many empirical studies have identified some memory properties of financial markets. To name a few, [START_REF] Gopikrishnan | Statistical properties of share volume traded in financial markets[END_REF]][Bouchaud et al., 2009] 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. As a consequence, Hawkes processes have come up as a natural modelling choice, and triggered a lot of interest in recent studies on market microstructure and limit order books. [START_REF] Bacry | Hawkes model for price and trades highfrequency dynamics[END_REF], Bacry et al., 2013, Bacry et al., 2016] propose various models of price and order flow models. An extensive survey of the application of Hawkes process in finance can be found in [START_REF] Bacry | Hawkes processes in finance[END_REF]. In the specific context of limit order books, [Large, 2007] 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 [Abergel andJedidi, 2015][Zheng et al., 2014] 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.

This paper is a contribution to this latter strand of research.

In 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. They will be the main objects under scrutiny in the present work as well.

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: 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-processbased 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.

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, some 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, 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"2 . 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 relationships between events, and ratios greater than two are highlighted. 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:

Results of Table

• 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 M 1 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. 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.

M 1 , 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.

Table 1: Event types definitions

• 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 cancel-

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
Table 3: Conditional probability leverage lations 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 triggers 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 has now become 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 rich 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 }.
Recall that 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 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 modelled as 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 ticks3 . 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 ∈E down 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 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.

Here, we build on the results of Section 2 and study two classes models respectively based 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 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] Hawkes | A cluster process representation of a Self-Exciting process[END_REF]][Massoulié, 1998][Brémaud and Massoulié, 1996] for a more in-depth presentation of these processes, and to [Abergel andJedidi, 2015][Zheng et al., 2014] for their use in 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 [Hawkes andOakes, 1974][Massoulié, 1998] 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 φ(ts)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. [Massoulié, 1998]. 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 differ by 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.

Goodness of fit

It is well-known, see e.g. [Bowsher, 2007] that the transformed durations {τ } i of a Hawkes process

τ m i = T i+1
T i λ 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, 1-exponential 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

L 0 buy L 1 buy L 0 sell L 1 sell C 0 buy C 1 buy C 0 sell C 1 sell 0 5 M 0

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-processbased 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.

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 4 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:

Pair

• 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.

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 4 have norms equal to 0. 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.

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 C 1 buy 0.

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. [Brémaud andMassoulié, 1996][Zhu, 2015] for general results on nonlinear 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.

3.4 Performances of the nonlinear Hawkes models It appears clearly, simply by eyeballing the graphs, that the nonlinear Hawkes model leads to a 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.

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.

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)).

-6 -5 -4 -3 -2 -1 0 1 2 0.0 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 [Ozaki, 1979][Rubin, 1972], 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

+ T i α β e -β(T D -T i ) -1 + T i ln   µ + α T j <T i e -β(T i -T j )   ,
and, letting β tend to ∞, there holds

lim β→+∞ ln L(λ, {T }) = -µT + N (T ) ln µ,
which is finite. However, a strictly concave continuous function having a local maximum cannot tend to a finite limit at infinity. 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 Algorithm 1.

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:

for i ← 1, N do 7:

// Mutation 8:

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 r 1 ,g + F (x r 2 ,gx r 3 ,g ); 10:

// Crossover. Construct trial element u i,g+1

11:

I rand is a random integer from 1, D ;

12:

for j ← 1, D do 13:

rand j,i ∼ U(0, 1); 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 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

v i,g = x r 1 ,g + F (x r 2 ,g -x r 3 ,g )
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 r 2 ,gx r 3 ,g ) or use more points for deviation

v i,g = x r 1 ,g + F (x r 2 ,g -x r 3 ,g ) + F (x r 4 ,g -x r 5 ,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 non-hybrid 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 currenttopbest/1 strategy, where the new donner vectors are obtained by v i,g = x i,g + F i (x pbest,gx i,g ) + F i (x r 1 ,gx r 2 ,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,r i , 0.1), CR i = randn i (M CR,r i , 0.1)1 {M CR,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 r 2 ,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.

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-topbest/1". The parameter adaptation is also combined with "rand/1" to provide better performances.

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 µ mmin µ m < e a , σ(ρ mn ) < e r ρ mn or max ρ mnmin ρ 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 Markovian 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 estimation. 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 for ALLIANZ SE.

3. 4 . 1

 41 Goodness of fitSimilarly 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 Figure4.

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 :Figure 7 :

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

Figure 8 :

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

  |µ 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 2 :

 2 Conditional probabilities of occurrences per event type

	buy	27.94 9.73 26.24 20.99 1.09	0.50	4.98	1.78 0.08 3.60 2.65	0.44
	L 0 sell	9.63 28.43 20.36 26.50 0.58	1.01	1.83	4.94 3.55 0.08 0.43	2.68
	C 0 buy	21.64 24.72 29.33 7.51	0.66	0.57	2.92	4.28 4.97 1.32 0.97	1.11
	C 0 sell	24.04 21.97 7.45 29.82 0.64	0.58	4.30	2.85 1.32 5.00 1.08	0.95
	M 0 buy	20.69 8.02	6.96 11.18 9.08	0.59	9.64	0.86 1.12 6.65 24.27 0.94
	M 0 sell	7.19 20.72 10.18 6.90	0.64	9.63	0.72	9.39 6.67 1.10 0.93 25.92
	L 1 buy	32.48 10.83 1.17 26.57 0.92	0.94	4.38	1.68 8.89 4.93 2.43	4.77
	L 1 sell	10.24 33.61 26.27 1.14	0.94	0.90	1.71	4.32 4.88 8.97 4.54	2.47
	C 1 buy	14.46 12.40 51.27 4.59	0.26	0.10	8.42	4.06 2.83 0.96 0.51	0.15
	C 1 sell	11.85 14.60 4.32 52.25 0.10	0.22	3.96	8.41 0.91 2.77 0.14	0.48
	M 1 buy	12.23 6.18	4.56 30.18 1.04	0.64 24.94 4.39 1.16 8.70 3.35	2.64
	M 1 sell	5.93 12.67 29.80 4.63	0.71	1.09	4.36 24.68 8.88 1.13 2.59	3.52
	O	19.93 20.42 20.23 20.74 0.79	0.73	4.02	3.99 2.93 2.95 1.62	1.65

Table 4 :

 4 Conditional probability comparison between simulated order flows and real data

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

Table 5 :

 5 Kernel L 1 norm medians in the 2-exponential model Moreover, two other event pairs come out of the calibration with 0 kernel norms, M 0

	1563 0.2357 0.9392 0.0914	0	0	0.3845 0.1607	0	0	0.0013	0
	C 1 sell |C 1 buy .						buy and buy |C 1

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.

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 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

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.

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".

- -1980 -1965 -1950 lnL 6 ( 6 , {T i , X i } i 12 )

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.

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:

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

• 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


Ahmed Belhadjayed

Grégoire Loeper

Frédéric Abergel

Ahmed Bel Hadj Ayed

FORECASTING TRENDS WITH ASSET PRICES

come

Introduction

In theory, asset prices could be well described by random walks, according to the Efficient Market Hypothesis. If this were indeed the case, then future returns would not be predictable. Nevertheless, professionals in the finance industry tend to have divergent views on the subject, and trend following strategies are the principal sources of returns for Commodity Trading Advisors (see [START_REF] Lempérière | Two centuries of trend following[END_REF]). In fact, most quantitative strategies are based on the more or less explicit assumption that the trends of assets are known (see [START_REF] Markovitz | Portfolio Selection[END_REF], [START_REF] Merton | Continuouse-Time Finance[END_REF]) and can be extracted from the asset prices themselves. It is therefore natural to address the question of forecasting asset trends, and to seek to provide reliable statistical estimators. Unfortunately, the estimation of the trend of an asset is a statistically difficult problem, mainly because of a high measurement noise: consider for example a simple model with a constant trend dSt St = µdt + σ S dW S t . Then, the best estimate of the trend in the least square sense at time T is given by µ t = 1 T T 0 dSu Su . Student's t-test will reject the hypothesis

µ = 0 if | µ T | > 1.96σ S √ T
at a 5% significance level. Therefore, with σ S = 30%, the estimate µ T = 1% becomes statistically relevant for observation times T > 3457 years! The purpose of this work is to assess the feasibility of forecasting trends modelled by an unobserved mean-reverting diffusion. Using a Bayesian approach or maximum likelihood estimation (see e.g. [START_REF] Leroux | Maximum-likelihood estimation for hidden markov models[END_REF], [START_REF] Cappé | Inference in Hidden Markov Models[END_REF], [START_REF] Benmiloud | Parameter Estimation in Hidden Markov Chains and Segmentation of Images[END_REF], [START_REF] Casarin | Online data processing: comparison of Bayesian regularized particle filters[END_REF] or [START_REF] Dahia | Nouvelles méthodes en filtrage particulaire. Application au recalage de navigation inertielle par mesures altimétriques[END_REF]), inference methods for partially observed processes have been applied to financial time series, mostly in the framework of stochastic volatility models (see [START_REF] Jacquier | Bayesian analysis of stochastic volatility models[END_REF], [START_REF] Eraker | Mcmc analysis of diffusion models with application to finance[END_REF], [START_REF] Kim | Stochastic volatility: Likelihood inference and comparison with arch models[END_REF] or [START_REF] Chib | Markov chain monte carlo methods for stochastic volatility models[END_REF]). Closer in spirit to our motivation, several authors have considered the situation of an unobservable stochastic trend, and use filtering methods in this context (most of these methods are introduced in [START_REF] Genay | An Introduction to Wavelets and Other Filtering Methods in Finance and Economics[END_REF]). For example, in [START_REF] Whittle | Prediction and Regulation by Linear Least-Square Methods[END_REF], [START_REF] Maravall | Stochastic linear trends[END_REF] and [START_REF] Pollock | Wiener-kolmogorov filtering, frequency-selective filtering, and polynomial regression[END_REF] the Wiener-Kolmogorov filter (see [START_REF] Wiener | Extrapolation, Interpolation, and Smoothing of Stationary Time Series With engineering applications[END_REF] for details) is used; in [START_REF] John | Optimal properties of exponentially weighted forecasts[END_REF], [START_REF] Holt | Forecasting seasonals and trends by exponentially weighted moving averages[END_REF] and [START_REF] Harvey | Forecasting, structural time series models and the Kalman filter[END_REF] the trend is supposed to be a random walk and the Kalman filter is used, and in [START_REF] Harvey | General model-based filters for extracting cycles and trends in economic time series[END_REF] the Butterworth filter (see [START_REF] Butterworth | On the theory of filter amplifiers[END_REF] for details) is applied. As it turns out, most of these filters are based on a parametric stochastic model for the trend, and their usefulness in realistic trading strategies is therefore confronted to the problem of parameters estimation (see [START_REF] Lakner | Optimal trading strategy for an investor: the case of partial information[END_REF], [START_REF] Pham | Optimal Portfolio in Partially Observed Stochastic Volatility Models[END_REF], [START_REF] Laskry | Contrôle stochastique avec informations partielles et applications la finance[END_REF] or [START_REF] Brendle | Portfolio selection under incomplete information[END_REF] where filtering methods are used in the context of trading strategies). It is our aim to partly fill the existing gap in the quantitative finance literature, and shed a new light on the feasability of classical trading strategies based on the determination of asset trends. Such a problem is addressed e.g. in [START_REF] Harvey | Forecasting, structural time series models and the Kalman filter[END_REF] and [START_REF] Harvey | Estimation procedures for structural time series models[END_REF], where the emphasis is set on Kalman filtering and the asymptotic behaviour of the maximum likelihood estimator. In this paper, we focus on the question of parameters estimation with a unobserved mean-reverting trend, and measure the effect of parameters misspecification. The paper is organized as follows: in Section 1, we present the model and recall some results on Kalman filtering. Section 2 is devoted to parameters inference using discrete time observations. The performance of statistical estimators is evaluated by giving their asymptotic behaviours, and by providing, in closed form, the Cramer-Rao bound. Section 3 introduces the continuous time misspecified Kalman filter. We provide estimates for the impact of parameters misspecification on trend filtering, and compute the probability to have a positive trend, given a positive estimate. Finally, Section 4 contains numerical examples illustrating the relevance of parameters misspecification in trend filtering.

Framework

In this section, the model for the asset price and the mean-reverting dynamics of its trend is made precise. Then, the Kalman filtering method, in a time-discretized version, is recalled.

1.1. Model.

1.1.1. Continuous time model. Consider a financial market living on a stochastic basis (Ω, F, F, P), where F = {F t , t 0} is the natural filtration associated to a two-dimensional, uncorrelated Wiener process (W S , W µ ), and P is the objective probability measure. The dynamics of the risky asset S is given by

dS t S t = µ t dt + σ S dW S t , (1)
dµ t = -λ µ µ t dt + σ µ dW µ t , (2)
with µ 0 = 0. We also assume that (

λ µ , σ µ , σ S ) ∈ R * + × R * + × R * + . Denote by F S = F S t
be the natural filtration associated to the price process S. An important point is that only F S -adapted processes are observable, which implies that agents in this market do not observe the trend µ.

Remark 1.1. The linear and Gaussian framework has well-known shortcomings: in practice, financial asset returns are heavy-tailed (see [START_REF] Mandelbrot | The variation of certain speculative prices[END_REF] [START_REF] Chakraborti | Econophysics review: I. empirical facts[END_REF]) because of jumps and volatility fluctuations. Moreover, interactions between the trend and the volatility processes are possible. In this paper, we focus on a class of simple models so as to extract easily interpreted analytical expressions for various quantities of interest in the context of trend forecasting and investment strategies. 2) is now presented. Let then δ be a discrete time step, and denote by the subscript k the value of a process at time t k = kδ. The discrete time model is:

Discrete time model. A discrete-time version of (1)-(

y k+1 = S k+1 -S k δS k = µ k+1 + u k+1 , (3)
µ k+1 = e -λµδ µ k + v k , (4)
where

u k ∼ N 0, σ 2 S δ and v k ∼ N 0, σ 2 µ 2λµ 1 -e -2λµδ
. The system (3)-( 4) corresponds to an AR(1) model with noise.

1.2. Optimal trend estimator.

1.2.1. Discrete Kalman filter. In this subsection, the parameters θ = (λ µ , σ µ ) and σ S are supposed to be known. The discrete time system (3)-( 4) corresponds to a Linear Gaussian Space State model where the observation is y and the state, µ (see [START_REF] Brockwell | Introduction to time series and forecasting[END_REF] for details). In this case, the optimal estimator is the conditional expectation E [µ k |y 1 , ..., y k ], given by the Kalman filter (classical results on the discrete Kalman filter are recalled in Appendix A). Applications of this filter in Finance are numerous. Indeed, it can be used for trend filtering (see [START_REF] Harvey | Forecasting, structural time series models and the Kalman filter[END_REF]), for term structure models (see [START_REF] Manoliu | Energy futures prices: term structure models with kalman filter estimation[END_REF] or [START_REF] Babbs | Kalman filtering of generalised vasicek term structure models[END_REF]), for trading strategies (see [START_REF] Dilellio | A Kalman filter control technique in mean-variance portfolio management[END_REF] or [START_REF] Elliott | Pairs trading[END_REF]) and many other applications (see [START_REF] Wells | The Kalman Filter in Finance[END_REF], [START_REF] Bebendo | The dynamics of the volatility skew: A Kalman filter approach[END_REF], [START_REF] Das | Market risk beta estimation using adaptive kalman filter[END_REF] or [START_REF] Stasinakis | Forecasting us unemployment with radial basis neural networks, kalman filters and support vector regressions[END_REF]). For simplicity, we let X k/l denote E [X k |y 1 , ..., y l ]. The Kalman filter is decomposed in two distinct phases:

(1) An a priori estimate given μk+1/k and Γ k+1/k = E [(µ k+1 -μk+1/k )(µ k+1 -μk+1/k ) T . This estimate is done using the transition equation ( 4). (2) An a posteriori estimate. When the new observation is available, a correction of the first estimate is done to obtain μk+1/k+1 and Γ k+1/k+1 = E (µ k+1 -μk+1/k+1 )(µ k+1 -μk+1/k+1 ) T . The criterion for this correction is the least squares method. Thus, μk/k is the minimum variance linear unbiased estimate of the trend µ k . Formally, the iterative method is given by:

µ k+1/k+1 = e -λµδ µ k/k + K k+1 y k+1 -e -λµδ µ k/k , (5)
Γ k+1/k+1 = (1 -K k+1 ) Γ k+1/k , (6)
with

K k+1 = Γ k+1/k Γ k+1/k + σ 2 S δ , Γ k+1/k = e -2λµδ Γ k/k + σ 2 µ 2λ µ 1 -e -2λµδ .
1.2.2. Stationary limit and continuous-time representation. Solving the equation Γ k+1/k+1 = Γ k/k corresponding to the steady-state yields:

Γ ∞ = g (σ S , λ µ , σ µ ) -f (σ S , λ µ , σ µ ) 2e -2λµδ , where f (σ S , λ µ , σ µ ) = σ 2 S δ + σ 2 µ 2λ µ 1 -e -2λµδ ,
and

g (σ S , λ µ , σ µ ) = f (σ S , λ µ , σ µ ) 2 + 2σ 2 S σ 2 µ λ µ δ (e -2λµδ -e -4λµδ
).

Using the stationary covariance error Γ ∞ , a stationary gain K ∞ is defined and the estimate can be rewritten as a corrected exponential average:

µ n+1 = K ∞ ∞ i=0 e -λµδi (1 -K ∞ ) i y n+1-i . (7)
The steady-state Kalman filter has also a continuous-time limit that depends on the asset returns. This result is recalled in the Proposition 1. The steady-state Kalman filter µ solves the following stochastic differential equation

d µ t = -λ µ β (λ µ , σ µ , σ S ) µ t dt + λ µ (β (λ µ , σ µ , σ S ) -1) dS t S t , (8)
where

β (λ µ , σ µ , σ S ) = 1 + σ 2 µ λ 2 µ σ 2 S 1 2 . ( 9
)
Proof. Based on [START_REF] Lakner | Optimal trading strategy for an investor: the case of partial information[END_REF], the Kalman filter is given by:

E µ t |F S t = φ (t) µ 0 + 1 σ 2 S t 0 P (u) φ (u) dS u S u , φ (t) = e -λµt-1 σ 2 S t 0 P (u)du
, where the estimation error variance P is the solution of the following Riccati equation:

P (t) = -1 σ 2 S P (t) 2 -2λ µ P (t) + σ 2 µ .
In this steady-state regime, we have P (t) = 0. Then, the positive solution of this equation is given by

P ∞ = σ 2 S λ µ (β (λ µ , σ µ , σ S ) -1)
, and there holds:

µ t = φ ∞ (t) µ 0 + 1 σ 2 S t 0 P ∞ φ ∞ (u) dS u S u , φ ∞ (t) = e -λµβ(λµ,σµ,σ S )t .
Since:

dφ ∞ (t) φ ∞ (t) = -λ µ β (λ µ , σ µ , σ S ) dt, Equation (8) follows.
The Kalman filter is the optimal estimator for linear systems with Gaussian uncertainty, and such a continuous-time representation can be used for risk/return analysis of trend following strategies (see [START_REF] Bruder | Risk-Return Analysis of Dynamic Investment Strategies[END_REF] for details). Of course, in practice, the parameters θ = (λ µ , σ µ ) are unknown and must be estimated. This important question is addressed in the next section.

Inference of the trend parameters

In this section, the problem of parameters inference is tackled, based on the use of statistical estimators such as Maximum Likelihood or Bayesian estimators. We analyze the asymptotic behaviours of statistical estimators and provide the Cramer-Rao bound in closed form. Since these classical estimators are based on the likelihood, two on-line computations of this function are presented in Appendix B.

2.1. Asymptotic behaviour of statistical estimators. The discrete time model ( 3)-( 4) can be reformulated using the following proposition:

Proposition 2. Consider the model (3)-(4) with (λ µ , σ µ , σ S ) ∈ R * + × R * + × R * + .
In this case, the process (y i ) is ARMA(1, 1). The asymptotic behaviour of the classical estimators follows. Indeed, the identifiability property and the asymptotic normality of the maximum likelihood estimator are well known for stationary ARMA Gaussian processes (see [START_REF] Brockwell | Time series: Theory and Methods[END_REF], section 10.8). Moreover, the asymptotic behaviour of the Bayesian estimators are also guaranteed by the ARMA(1, 1) property of the process (y i ). If the prior density function is continuous and positive in an open neighbourhood of the real parameters, the Bayesian estimators are asymptotically normal (see [START_REF] Tamaki | The bernstein-von mises theorem for stationary processes[END_REF] in which a generalized Bernstein-Von Mises theorem for stationary "short memory" processes is given, or [START_REF] Monahan | Fully bayesian analysis of arma time series models[END_REF] for a discussion on the Bayesian analysis of ARMA processes).

2.2.

Cramer-Rao bound. This bound is the lowest variance of the unbiased estimators. We recall the following result, providing a formal description of the Cramer-Rao bound (CRB in short).

, • • • , y N ) T . Suppose that (λ µ , σ µ , σ S ) ∈ R * + × R * + × R * + .
If θN is an unbiased estimator of θ = (λ µ , σ µ ), we have:

Cov θ θN CRB (θ) .
This bound is given by CRB (θ) = I -1 N (θ), where I N (θ) is the Fisher Information matrix:

(I N (θ)) i,j = -E ∂ 2 log f (y 1 , ...y N |θ) ∂θ i ∂θ j ,
and I N (θ) = N I 1 (θ). Moreover, the maximum likelihood estimator θ M L N attains this bound:

√ N θ M L N -θ → N 0, I -1 1 (θ)
. This result is a consequence of Proposition 2 (see [START_REF] Brockwell | Time series: Theory and Methods[END_REF], section 10.8). The following result is an analytic representation of the Fisher information matrix:

Theorem 2.2. For the model (3)-(4), if (λ µ , σ µ , σ S ) ∈ R * + × R * + × R * + , we have: I 1 (θ) = 1 4Π Π -Π f -2 θ (ω) ∂f θ ∂θ i (ω) ∂f θ ∂θ j (ω) dω 1≤i,j≤2
, where f θ is the spectral density of the process (y i ):

f θ (ω) = σ 2 µ 2λµ 1 -e -2λµδ + σ 2 S δ 1 + e -2λµδ - 2e -λµδ σ 2 S δ cos (ω)
1 + e -2λµδ -2e -λµδ cos (ω) .

Proof. Whittle's formula (see [START_REF] Whittle | The analysis of multiple stationary time series[END_REF] for details) gives the integral representation of the Fisher information matrix. Since the process (y i ) is ARMA(1, 1), the expression of its spectral density follows (see [START_REF] Brockwell | Time series: Theory and Methods[END_REF], section 4.4).

Finally, the Cramer-Rao Bound of the trend parameters can be computed using Theorem 2.2.

Impact of parameters misspecification

In this section, we consider the continuous-time Kalman filter with a bad calibration in the steady-state regime. The law of the residuals between the filter (misspecified or not) and the hidden process is characterized, and the impact of parameters misspecification on the detection of a positive trend is studied.

3.1. Context. Suppose that the risky asset S is given by the model ( 1)-( 2) with θ * = σ * µ , λ * µ , and suppose that an agent thinks that the parameters are equal to θ = (σ µ , λ µ ). Assuming the steady-state regime and using these estimates and Proposition 1, the agent implements the continuous-time misspecified Kalman filter:

d µ t = -λ µ β µ t dt + λ µ (β -1) dS t S t , (10)
where β = β (λ µ , σ µ , σ S ) (see Equation ( 9)) and µ 0 = 0. The following lemma gives the law of the misspecified Kalman filter:

Lemma 3.1. Consider the model (1)-( 2) with θ * = σ * µ , λ * µ . The misspecified, continuous-time filter of Equation ( 10) is given by:

µ t = λ µ (β -1) e -λµβt t 0 e λµβs µ * s ds + σ S t 0 e λµβs dW S s . (11)
Moreover, µ is a centered Gaussian process and its variance is given by:

Var [ µ t ] = E µ 2 t = λ 2 µ (β -1) 2 σ * µ 2 λ * µ λ µ β -λ * µ 1 -e -(λµβ+λ * µ) t λ µ β + λ * µ + 2e -(λµβ+λ * µ) t -e -2λ * µ t -e -2λµβt λ µ β -λ * µ + e -2λµβt -1 2λ µ β + λ µ (β -1) 2 σ 2 S 2β 1 -e -2λµβt .
Proof. Applying Itō's lemma to the function f ( µ, t) = e λµβt µ t , and integrating from 0 to t yields Equation [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. Therefore, µ is also a Gaussian process. Its mean is zero (because µ * 0 = 0). Since the processes µ * and W S are supposed to be independent, the variance of µ is given by the sum of the variances of the terms in Equation [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. Moreover:

Var t 0 e λµβ dW S s = e 2λµβt -1 2λ µ β , Var t 0 e λµβs µ * s ds = t 0 t 0 e λµβ(s 1 +s 2 ) Cov µ * s 1 , µ * s 2 ds 1 ds 2 ,
and Cov µ * s 1 , µ * s 2 is given by Equation [START_REF] Harvey | Forecasting, structural time series models and the Kalman filter[END_REF]. The variance of the process µ t follows.

3.2.

Filtering with parameters misspecification. The impact of parameters misspecification on trend filtering can be measured using the difference between the filter and the hidden process. The following theorem gives the law of the residuals. Theorem 3.2. Consider the model (1)-( 2) with θ * = σ * µ , λ * µ and the trend estimate defined in Equation [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. The process µ -µ * is a centered Gaussian process and its variance has a stationary limit:

lim t→∞ Var [ µ t -µ * t ] = σ 2 S 2β λ µ (β -1) 2 + λ * µ (β * ) 2 -1 λ * µ β + λ µ λ µ β + λ * µ ,( 12
)
where β = β (λ µ , σ µ , σ S ) and β * = β λ * µ , σ * µ , σ S as given in Equation (9). Moreover, if (σ µ , λ µ ) = σ * µ , λ * µ , Equation (12) becomes: lim t→∞ Var [ µ * t -µ * t ] = λ * µ σ 2 S (β * -1) . ( 13
)
Proof. Equation [START_REF] Cappé | Inference in Hidden Markov Models[END_REF] implies that the process µ -µ * is a centered Gaussian process. Its variance can be computed in closed form:

Var [ µ t -µ * t ] = Var [ µ t ] + Var [µ * t ] -2 * Cov [ µ t , µ * t ] ,
where

Var [µ * t ] = (σ * µ) 2 2λ * µ
1 -e -2λ * µ t , and Var [ µ t ] is given by Lemma 3.1. Since the processes W S and µ * are supposed to be independent, there holds:

Cov [ µ t , µ * t ] = λ µ (β -1) σ * µ 2 2λ * µ 1 -e -(λµβ+λ * µ) t λ µ β + λ * µ - e -2λ * µ t -e -(λµβ+λ * µ) t λ µ β -λ * µ .
The asymptotic variance is obtained by letting t → ∞:

lim t→∞ Var [ µ t -µ * t ] = λ µ (β -1) 2β (β -1) σ 2 S - σ * µ 2 (β + 1) λ * µ λ µ β + λ * µ + σ * µ 2 2λ * µ ,
and Equation ( 12) follows. Finally, Equation ( 13) is obtained by letting θ → θ * .

Remark 3.3. Consider the well-specified case (σ µ , λ µ ) = σ * µ , λ * µ . Using Equation [START_REF] Chakraborti | Econophysics review: I. empirical facts[END_REF], it follows that:

lim t→∞ Var [ µ * t -µ * t ] Var [µ * t ] = 2 1 + 1 + (σ * µ) 2 (λ * µ) 2 σ 2 S . (14)
Then, the asymptotic relative variance of the well-specified residuals is an increasing function of λ * µ and a decreasing function of σ * µ . 3.3. Detection of a positive trend. In practice, the trend estimate (misspecified or not) will be used to make an investment decision. For example, a positive estimate leads to a long position. So, it is interesting to estimate the probability of a positive trend knowing a positive estimate. We derive this probability in closed form, based on the following proposition giving the asymptotic conditional law of the trend (µ * t |μ t = x): Proposition 3. Consider the model (1)-( 2) with θ * = σ * µ , λ * µ and the trend estimate defined in Equation [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. Then, there holds:

(µ * t |μ t = x) L → t→∞ N M ∞ µ * |μ , Var ∞ µ * |μ , (15)
with:

M ∞ µ * |μ = λ * µ β (β * ) 2 -1 (β -1) λ µ β + λ * µ (β * ) 2 x, (16)
Var ∞ µ * |μ = Var ∞ µ * 1 - λ * µ λ µ β (β * ) 2 -1 λ * µ + λ µ β λ µ β + λ * µ (β * ) 2 , ( 17
)
where [START_REF] Dahia | Nouvelles méthodes en filtrage particulaire. Application au recalage de navigation inertielle par mesures altimétriques[END_REF] becomes:

Var ∞ µ * = (σ * µ ) 2 2λ * µ . Moreover, if (σ µ , λ µ ) = σ * µ , λ * µ , Equation
(µ * t |μ * t = x) L → t→∞ N x, 2Var ∞ µ * β * + 1 , ( 18
)
where β * = β λ * µ , σ * µ , σ S (see Equation ( 9)). Proof. Since the estimate μ and the trend µ * are two centered and correlated Gaussian processes (see Lemma 3.1 and the proof of Theorem 3.2), the conditional law (µ * t |μ t = x) is Gaussian with a mean and a variance given by:

M µ * t |μt = Cov (μ t , µ * t ) Var [μ t ] x, Var µ * t |μt = Var [µ * t ] - Cov (μ t , µ * t ) 2 Var [μ t ]
. Using Lemma 3.1 and the expression of Cov (μ t , µ * t ) in the proof of Theorem 3.2 yields

lim t→∞ M µ * t |μt = M ∞ µ * |μ , lim t→∞ Var µ * t |μt = Var ∞ µ * |μ ,
and Equation ( 15) follows. Finally, Equation ( 18) is obtained by letting θ → θ * .

The following proposition is a consequence of the previous proposition. It gives the asymptotic probability to have a positive trend, knowing a positive estimate equal to x. Proposition 4. Consider the model (1)-( 2) with θ * = σ * µ , λ * µ and the trend estimate defined in Equation [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. In this case:

lim t→∞ P (µ * t > 0|μ t = x) = P ∞ (µ * > 0|μ = x) , (19)
where

P ∞ (µ * > 0|μ = x) = 1 -Φ   -M ∞ µ * |μ=x Var ∞ µ * |μ=x   , (20)
where M ∞ µ * |μ=x and Var ∞ µ * |μ=x are defined in Equations ( 16) and ( 17), and Φ is the cumulative distribution function of the standard normal law.

Moreover, if x > 0 and (σ µ , λ µ ) = σ * µ , λ * µ , this asymptotic probability becomes an increasing function of σ * µ and a decreasing function of λ * µ . Proof. Equations ( 19) and ( 20) follow from Proposition 3. Now, consider the well-specified case (σ µ , λ µ ) = σ * µ , λ * µ and x > 0. Using Equation [START_REF] Dilellio | A Kalman filter control technique in mean-variance portfolio management[END_REF], it follows that:

Var ∞ µ * |μ * =x = f σ * µ , λ * µ , σ S , where
f σ * µ , λ * µ , σ S = σ * µ 2 λ * µ 1 + 1 + (σ * µ) 2 σ 2 S (λ * µ) 2 . Since ∂f ∂λ * µ σ * µ , λ * µ , σ S = -σ * µ 2 λ * µ 2 1 + 1 + (σ * µ) 2 σ 2 S (λ * µ) 2 + (σ * µ) 2 σ 2 S ≤ 0, ∂f ∂σ * µ σ * µ , λ * µ , σ S = λ * µ σ * µ σ 2 S 1 + (σ * µ) 2 σ 2 S (λ * µ) 2 σ * µ 2 + σ 2 S λ * µ 2 ≥ 0,
the asymptotic well-specified probability to have a positive trend, knowing a positive estimate equal to x is an increasing function of σ * µ and a decreasing function of λ * µ . Remark 3.4. This probability is an increasing function of x. Indeed, it is easier to detect the sign of the real trend with a high estimate than with a low estimate. Moreover, this probability is always superior to 0.5. This is due to the non-zero correlation between the trend and the filter. As shown in the previous sections, trend filtering is easier with a small spot volatility σ S . Here, the probability to make a good detection is also a decreasing function of σ S .

Simulations

In this section, numerical simulations are performed, in order to make the reader aware of the trend filtering problem. First, the feasibility of trend forecasting with statistical estimators is illustrated on different trend regimes. Then, the effects of a bad forecast on trend filtering and on the detection of a positive trend are discussed. Finally, we study the numerical maximum likelihood inference of this model using Kalman filtering and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm.

4.1. Feasibility of trend forecasting. Suppose that only discretetime observations are available and that the time step is equal to δ = 1/252, corresponding to daily returns. Note that, if this time step is divided by M , the variance of the observation noise u k is multiplied by M (u k ∼ N 0, σ 2 S δ ), so that increasing the frequency of the observations will not be really helpful to calibrate the trend. We also assume that the agent uses an unbiased estimator. Given T years of observations, The Cramer-Rao bound is given by:

CRB T (θ) = I -1 1 (θ) T * 252 ,
where I 1 (θ) is given by Theorem 2.2. The smallest confidence region is obtained with this matrix. In practice, the real values of the parameters θ are unknown and asymptotic confidence regions are computed (replacing θ by the estimates θ in the Fisher information matrix I 1 θ ).

Since the goal of this subsection is to evaluated the feasibility of this estimation problem, we suppose that we know the real values of the parameters. In such a case, the exact Cramer-Rao bound can be computed. Suppose that a target standard deviation x i is fixed for the parameter θ i . In this case, to reach the precision x i , the length of the observations must be larger than

T x i = (I -1 1 (θ)) ii 252 * x 2 i
. We consider a fixed spot volatility σ S = 30%, two target precisions for each parameter θ i and we compute T x i for several configurations. Figures 1,2, 3 and 4 represent the results. It is well-known that for a high measurement noise, which means a high spot volatility, the problem is harder because of a low signal-to-noise ratio. The higher the volatility, the longer the observations must be. Here, we observe that with a higher drift volatility σ µ and a lower λ µ , the problem is easier. Indeed, the drift takes higher values and is more detectable. Moreover, the simulations show that the classical estimators are not adapted to such a weak signal-to-noise ratio: even after a long period of observations, the estimators exhibit high variances. Indeed, the shortest observation period is longer than 29 years. It corresponds to a target standard deviation equal to 0.5 for a parameter λ µ = 1, and a trend standard deviation equal to σ µ (2λ µ ) -1/2 ≈ 63%. Therefore, for this configuration, after 30 years of observations, the standard deviation is equal to 50% of the real parameter value λ µ . After 742 years, this standard deviation is equal to 10%. Even in such a regime, the trend forecast with a good precision is impossible. This subsection illustrates the impact of parameters misspecification on trend filtering. Using the results of Theorem 3.2, we represent, for different configurations, and for the well-and mis-specified case, the asymptotic standard deviation of the residuals between the trend and the filter. Figures 5 and6 represent the asymptotic standard deviation of the trend and of the residuals in the well-specified case (the agent uses the real values of the parameters) for different configurations. As seen in Equation ( 13), the asymptotic standard deviation of the wellspecified residuals is an increasing function of the drift volatility σ * µ and a decreasing function of the parameter λ * µ . For λ * µ = 1 and σ * µ = 90%, the standard deviation of the residuals ( 44%) is inferior to the standard deviation of the trend ( 64%). For a high λ * µ and a small drift volatility, the two quantities are approximately equal. This figure leads to the same conclusions than Equation [START_REF] Chib | Markov chain monte carlo methods for stochastic volatility models[END_REF]. Indeed, as for the calibration problem, the problem of trend filtering is easier with a small λ * µ and a high drift volatility σ * µ . Now consider the worst configuration σ S = 30%, λ * µ = 5 and σ * µ = 10%. Figure 7 represents the asymptotic standard deviation of the residuals for different estimates (λ µ , σ µ ). This regime corresponds to a standard deviation of the trend equal to σ * µ 2λ * µ -1/2 ≈ 3.2% and to a standard deviation of the residuals equal to 3.16% in the well-specified case. If the agent implements the Kalman filter with λ µ = 1 and σ µ = 90%, the standard deviation of the residuals becomes superior to 25%. Finally, consider the best configuration σ S = 30%, λ * µ = 1 and σ * µ = 90%. Figure 8 represents the asymptotic standard deviation of the residuals for different estimates (λ µ , σ µ ). This regime corresponds to a standard deviation of the trend equal to

σ * µ 2λ * µ -1/2 ≈ 63%
and to a standard deviation of the residuals equal to 44% in the wellspecified case. If the agent implements the Kalman filter with λ µ = 5 and σ µ = 10%, the standard deviation of the residuals becomes larger than 60%. Even in a favourable regime, the impact of parameters misspecification on trend filtering is not negligible.

-giving the asymptotic probability to have a positive trend, knowing a trend estimate equal to a threshold x -is illustrated.

In order to compare this probability for different trend regimes, we choose a threshold equal to the standard deviation of the filter μ: this quantity is tractable in practice and moreover, since the continuous time misspecified filter μ is a centered Gaussian process, the probability that μ becomes larger than its standard deviation is independent of the parameters σ * µ , λ * µ , σ µ , λ µ , σ S . Suppose first that the agent uses the real values of the parameters, and consider the asymptotic probability P (µ * > 0|μ * = V μ * to have a positive trend, knowing an estimate equal to its standard deviation. Figure 9 represents this probability for different regimes. As seen in Proposition 4, in the well-specified case, this probability is an increasing function of the trend volatility σ * µ and a decreasing function of λ * µ . Again, as in the calibration and filtering problems, the detection is easier with a small λ * µ and a high drift volatility. Now, suppose that the agent uses wrong estimates (σ µ , λ µ ). In this case, the agent implements the continuous time misspecified Kalman filter. Figures 10 and11 represent the asymptotic probability P µ * > 0|μ = V μ for the best and the worst configuration of Figure 9. As explained in Remark 3.4, this probability is always superior to 0.5, even with a bad calibration of the parameters. For each case, the probability to have a positive trend, knowing an estimate equal to its standard deviation, does not vary a lot with an error on the parameters. This quantity seems to be robust to parameters misspecifications. In this subsection, we focus on the numerical maximum likelihood inference of the model ( 1)-( 2) using Kalman filtering and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. We first illustrate the performance of a parametric bootstrap, and then give a typical application on real market data. and suppose that we want to use the maximum likelihood method. This model is Gaussian and the likelihood is thus characterized by the moments of first and second order. When the number of observations N is large, these moments can be difficult to compute numerically (to mitigate these problems, two methods are introduced in Appendix B).

Here, we focus on the computation using the Kalman filter. Since we are able to compute recursively this likelihood with a large number of observations, we choose the direct maximization of this function using the BFGS algorithm. An alternative way is to use the EM-algorithm (see [START_REF] Frydman | Maximum Likelihood Estimation of Hidden Markov Processes[END_REF] or [START_REF] Dembo | Parameter Estimation of Partially Observed Continious Time Stochastic Processes Via EM Algorithm[END_REF] for details). The aim of these tests is to compare the standard deviation of this numerical maximum likelihood estimator to the theoretical results of Theorem 2.2. So, we consider the following parametric bootstrap:

(1) Simulate M paths of length T years with the model ( 3)-( 4).

(2) Consider only one parameter (for example θ = σ µ ), the other one is supposed to be known. (3) On each path i ∈ {1..M }, we compute a numerical maximum likelihood estimate θML T,i with Kalman filtering and the BFGS algorithm (initialised with σ 0 µ = 0.1 or λ 0 µ = 0.1).

(4) The variance of θ M L T is estimated by:

V ar θ M L T = 1 M -1 M i=1 θML T,i - 1 M M k=1 θML T,k 2
. As a conclusion, one can see that for these particular examples, the asymptotic behaviour of the numerical maximum likelihood estimator is close to the asymptotic behaviour of the theoretical maximum likelihood estimator. 4.4.2. Example of application. This paragraph gives an example of application of this estimator on real market data. We consider the components {S j , j = 1..50} of the Eurostoxx 50 on 11/11/2014. The backtest period is from 6/19/2008 to 11/11/2014. In this test, we assume a risk free rate equal to zero and that each stock is tradable at its closing price and without transaction costs. We consider two self-financing portfolios {P i , i = 1, 2}:

Figures
dP i t P i t = 50 j=0 ω j,i t dS j t S j t , P i 0 = 100, based on two trend indicators {(α i t (S j )) , i = 1, 2}.
Let s be a threshold, for each portfolio {P i , i = 1, 2}, the allocation ω j,i t is given by the following rules:

(1) At each time t, we define:

N i,+ t = j ∈ {1..50} \ α i t S j > s , N i,- t = j ∈ {1..50} \ α i t S j < -s . ( 2
) If N i,+ t = ∅ N i,- t = ∅: ∀j ∈ {1..50} , ω j,i t = 0, else: ∀j ∈ {1..50} , ω j,i t = 1 #N i,+ t 1 j∈N i,+ t - 1 #N i,- t 1 j∈N i,- t Since 50
j=0 ω j,i t = 0, the two portfolios are market neutral. We consider the following trend indicators:

(1) (α 1 t (S j )) is the annualised sliding moving average estimator of daily returns:

α 1 t S j = 252 T T k=1 dS j t-k S j t-k
, with T = 252 business days. (2) (α 2 t (S j )) is the Kalman filter estimator of the model ( 1)-( 2) where the parameters λ t µ , σ t µ , σ t S are given by the numerical maximum likelihood (introduced above) over the last 252 business days. The initialisation of the BFGS algorithm at time t > 0 is given by the previous estimate λ t-1 µ , σ t-1 µ , σ t-1 S and the first one is equal to (0.1, 10%, 30%). Figure 14 illustrates the results with s = 10%. In this example, the portfolio with Kalman filtering outperforms the one with the sliding moving average. Moreover, Figure 14 shows that the annualised Sharpe ratio (see [START_REF] Sharpe | Mutual fund performance[END_REF]) of the Kalman filtering portfolio, computed with the daily returns, is higher and more robust to a threshold variation than the one with the sliding moving average. Satisfactory as the results may seem on this specific dataset, it is however important to note that the comparison performed in this section concerns only a very specific strategy. In a related work, a thorough study of single asset investment strategies based on various criteria is made, and the robustness of the investment strategy with respect to the model parameters is actually challenged. The interested reader is referred to [START_REF] Bel Hadj Ayed | Robustness of the optimal trading strategy[END_REF].

Conclusion

The present work tries to illustrate the difficulty of trend filtering in a model based on an unobserved mean-reverting diffusion.

This model belongs to the class of Linear Gaussian Space State models. The advantage of this kind of system is to have an on-line method of estimation: the Kalman filter.

In practice, the parameters of the model are unknown, and the calibration of filtering parameters becomes crucial. The linear and Gaussian case allows to compute, in closed form, the likelihood. The Kalman filter can also be used for this analysis. These methods can be generalized to a non-constant volatility, and classical estimators can be easily implemented in practice.

Although this framework is particularly convenient for forecasting, the results of the analysis show that the classical estimators are not adapted to such a weak signal-to-noise ratio. The horizons of observations needed for an acceptable precision are too long. Therefore, the convergence is not guaranteed and the impact of misspecification on trend filtering is not negligible.

Despite these difficulties, the non-zero correlation between the trend and its estimate can be used for trend detection. This fact can explain why trend following strategies work well in practice.

Appendix A: discrete Kalman filter

Framework. This section is based on [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF]. The discrete Kalman filter is a recursive method. Consider two objects: the observations {y k } and the states of the system {x k }. This filter is based on a Gauss-Markov first order model. Consider the following system:

x k+1 = F k x k + v k , y k = H k x k + u k .
The first equation is an a priori model, the transition equation of the system. The matrix F k is the transition matrix and v k is the transition noise. The second equation is the measurement equation. The matrix H k is named the measurement matrix and u k is the measurement noise.

The aim is to identify the underlying process {x k }. The two noises are supposed white, Gaussian, centered and decorrelated. In particular:

E u k v k u l v l T = R u k 0 0 R v k δ kl .
The two noises are also supposed independent of x k and the initial state is Gaussian. So, it can be proved with a recurrence that all states are Gaussian. Therefore, just the mean and the covariance matrix are needed for the characterization of the state. The estimation is given by two steps. The first one is an a priori estimation given xk+1/k and Γ k+1/k = E (x k+1 -xk+1/k )(x k+1 -xk+1/k ) T . When the new observation is available, a correction of the estimation is done to obtain xk+1/k+1 and Γ k+1/k+1 = E (x k+1 -xk+1/k+1 )(x k+1 -xk+1/k+1 ) T . This is the a posteriori estimation. The criterion considered for the a posteriori estimation is the least squares method, which corresponds to the minimization of the trace of Γ k+1/k+1 .

Filter. The prediction (a priori estimation) is given by

xk+1/k = F k xk/k , Γ k+1/k = F k Γ k/k F T k + R v k .
The a posteriori estimation is a correction of the a priori estimation. A gain is introduced to do this correction: xk+1/k+1 = xk+1/k + K k+1 y k+1 -H k+1 xk+1/k .

As explained above, the gain K k+1 is found by least squares method, which corresponds to ∂trace Γ k+1/k+1 ∂K k+1 = 0.

Using straightforward matrix differentiation, the gain is found: -1 e λµδ +e -λµδ -1

K k+1 = Γ k+1/k H T k+1 H k+1 Γ k+1/k H T k+1 + R u
0 ••• ••• 0 -1 e λµδ                               .
Therefore, the inverse of the matrix Σ µ 1:N +1 is given by:

Σ -1 µ 1:N +1 |θ =         
2λµ(e λµδ +e -λµδ ) Procedure. Finally, at time t, the inverse of the covariance matrix is given by the following protocol:

• Computation of the matrix A -1 t using A -1 t-1 . • Computation of the matrix Σ -1 µ 1:t |θ using Σ -1 µ 1:t-1 |θ . • Using Σ -1 y 1:t |θ = Σ -1 µ 1:t |θ -Σ -1 µ 1:t |θ A -1 t Σ -1 µ 1:t |θ , the matrix Σ -1 y 1:t |θ is obtained.

Determinant of the covariance matrix. The iterative computation of det Σ y 1:N |θ is based on the following lemma:

Corollary 2 . 1 .

 21 Consider the model (3)-(4) and N observations (y 1

Figure 1 .

 1 Figure 1. Time to reach a target standard deviation on σ µ equal to 0.05 (ln(years))

Figure 2 .

 2 Figure 2. Time to reach a target standard deviation on σ µ equal to 0.01 (ln(years))

Figure 3 .

 3 Figure 3. Time to reach a target standard deviation on λ µ equal to 0.5 (ln(years))

Figure 4 . 4 . 2 .

 442 Figure 4. Time to reach a target standard deviation on λ µ equal to 0.1 (ln(years))

Figure 5 .Figure 6 .

 56 Figure 5. Asymptotic standard deviation of the trend as a function of the trend parameters with σ S = 30%.

Figure 7 .Figure 8 . 4 . 3 .

 7843 Figure 7. Asymptotic standard deviation of the residuals of the misspecified Kalman filter as a function of the trend estimate parameters with σ S = 30%, λ * µ = 5 and σ * µ = 10%.

Figure 9 .

 9 Figure 9. Asymptotic probability to have a positive trend given a well-specified estimate equal to its standard deviation with σ S = 30%.

Figure 10 .Figure 11 . 4 . 4 .

 101144 Figure 10. Asymptotic probability to have a positive trend given a misspecified estimate equal to its standard deviation with σ S = 30%, λ * µ = 1 and σ * µ = 90%.

4. 4 . 1 .

 41 Parametric bootstrap. Consider the discrete time model (3)-(4)

Figure 12 .Figure 13 .

 1213 Figure 12. Standard deviation for the estimation of λ * µ = 1 with σ * µ = 90% and σ S = 30% on a double logarithmic scale.

  /2008 3/31/2009 1/4/2010 9/28/2010 7/1/2011 3/9/2012 12/4/2012 9/3/2013 6/2/2014 Portfolio with sliding moving average Portfolio with Kalman filtering

Figure 14 .Figure 15 .

 1415 Figure 14. Portfolio evolutions with a threshold s = 10%.

k+1 - 1 ,

 1 Γ k+1/k+1 = (I d -K k+1 H k+1 ) Γ k+1/k .

σ 2 µ

 2 ( e λµδ -e -λµδ ) -2λµ σ 2 µ( e λµδ -e -λµδ ) 0 • • • 0 -2λµ σ 2µ( e λµδ -e -λµδ ) 0

  [START_REF] Casarin | Online data processing: comparison of Bayesian regularized particle filters[END_REF] and 13 illustrate the results for the estimation of the trend mean reversion speed and of the trend volatility with M = 1000, σ *

								µ =
	90%, λ * µ = 1 and σ S = 30%.
			Standard deviation for the estimation of the trend mean reversion speed
		-2.55	+				-	with CRB
		-2.60	+				+	with bootstrap
		-2.65		+
	standard deviation	-2.75 -2.70		+	+	+
		-2.80				+
								+
		-2.85						+
		-2.90
		6.2	6.3	6.4	6.5	6.6	6.7	6.8
					years

Appendix B: Likelihood computation

The likelihood can be computed using two methods. The first one is based on a direct calculus while the second method uses the Kalman filter. 5.0.3. Direct computation of the likelihood. A first approach is to directly compute the likelihood. The vectorial representation of the discrete time model ( 3)-( 4) is:

, are two independent Gaussian processes. Therefore the vector (y 1 , • • • , y N ) T , knowing θ, is also a Gaussian process. The likelihood is then characterized by the mean M y 1:N |θ and the covariance Σ y 1:N |θ :

where Σ u 1:

Since the drift µ is an Ornstein Uhlenbeck process, then:

Finally, the likelihood is given by:

Remark 5.1. When the dimension N is large, it is difficult and numerically unstable to directly invert the covariance matrix Σ y 1:N |θ and compute its determinant. An iterative approach, the details of which are given in Appendix C, is preferred.

5.0.4.

Computation of the likelihood using the Kalman filter. The likelihood can also be evaluated via the prediction error decomposition (see [START_REF] Schweppe | Evaluation of likelihood functions for gaussian signals[END_REF] for details):

where the conditional laws are given in the following proposition:

Proposition 5. The process (y n |y 1 , ...y n-1 , θ) is Gaussian:

The a posteriori estimate of the trend μn-1/n-1 and the covariance error Γ n-1/n-1 are given by Kalman filtering (see Equations ( 5) and ( 6)).

Proof. Since the process y n is Gaussian, so is the process (y n |y 1 , ...y n-1 , θ). Moreover, using Equations ( 3)-( 4), we have:

and

Remark 5.2. In practice, the volatility is not constant. However, if the volatility σ S is F S -adapted, the two methods can be suitably modified and implemented. This assumption is satisfied if the volatility is a continuous process.

Appendix C: Iterative methods for the inverse and the determinant of the covariance matrix

We provide here, for the sake of completeness, some standard iterative methods to compute the inverse and the determinant of a covariance matrix.

Inverse of the covariance matrix. The use of the Matrix Inversion Lemma on Equation ( 21) gives:

. Then, we have to compute the inverse of the matrices A N and Σ µ 1:N |θ .

Inverse of the matrix

N is computed. The matrix A N +1 can be broken into four sub-matrices:

where

Therefore, the matrix A N +1 can be inverted blockwise.

Inverse of the matrix Σ µ 1:N |θ . The following lemma is used (see [START_REF] Akesson | Discrete eigenfunction expansion of multidimensional brownian motion and the ornstein-uhlenbeck process[END_REF] for details): Lemma 5.3. Let µ be an Ornstein Uhlenbeck process with parameters θ = (λ µ , σ µ ). The covariance matrix of µ 1 , .., µ N is Σ µ 1:N |θ . Then:

Lemma 5.4. The determinant of the matrix Σ y 1:N |θ is given by:

and for N ≥ 2, we have:

where

µ (e λµδ -e -λµδ )

.

Proof. The multiplication of Equation ( 21) by Σ -1 µ 1:N |θ gives: