State-dependent Hawkes processes and their application to limit order book modelling

Abstract

We study statistical aspects of state-dependent Hawkes processes, which are an extension of Hawkes processes where a self- and cross-exciting counting process and a state process are fully coupled, interacting with each other. The excitation kernel of the counting process depends on the state process that, reciprocally, switches state when there is an event in the counting process. We first establish the existence and uniqueness of state-dependent Hawkes processes and explain how they can be simulated. Then we develop maximum likelihood estimation methodology for parametric specifications of the process. We apply state-dependent Hawkes processes to high-frequency limit order book data, allowing us to build a novel model that captures the feedback loop between the order flow and the shape of the limit order book. We estimate two specifications of the model, using the bid–ask spread and the queue imbalance as state variables, and find that excitation effects in the order flow are strongly state-dependent. Additionally, we find that the endogeneity of the order flow, measured by the magnitude of excitation, is also state-dependent, being more pronounced in disequilibrium states of the limit order book.

Keywords:

• Hawkes process
• High-frequency financial data
• Market microstructure
• Limit order book
• Maximum likelihood estimation
• Endogeneity

JEL Classification:

• C58
• C51
• G17

1. Introduction

Hawkes processes are a class of self- and cross-exciting point processes where events of different types may increase the rate of new events of the same or other type (Hawkes 1971, Laub et al. 2015), dispensing with the independence-of-increments property of Poisson processes. Given their ability to capture clustering and contagion effects, Hawkes processes have found numerous applications in finance in the last decade, most notably in the modelling of high-frequency financial data (Embrechts et al. 2011, Bacry et al. 2015). Concurrently, the last two decades have witnessed a major transformation of financial markets through the proliferation of electronic trading in order-driven markets, where traders submit buy and sell orders to an electronic trading platform that attempts to match the arriving orders (Gould et al. 2013). A significant proportion of the order flow is driven by high-frequency trading algorithms, whereby successive orders are sometimes separated in time only by a few microseconds. As billions of orders are submitted each day, this has given rise to a profusion of new market data to study, giving an unprecedented opportunity to statistically analyse and model price formation and market microstructure at the shortest possible timescales.

From the point of view of statistical modelling, a key objective here is to find an accurate, yet parsimonious, statistical dynamic description of the limit order book (LOB), that is, the record of orders for a given asset that currently remain unfilled (Gould et al. 2013). As one of the main approaches to LOB modelling, Hawkes processes have been used to describe the order flow, which inherently drives the evolution of the LOB (Large 2007, Bacry et al. 2016, Rambaldi et al. 2017, Cartea et al. 2018b, Lu and Abergel 2018). This approach proceeds by specifying a set of order types (i.e. orders classified in terms of their impact on the LOB) and by fitting a multivariate Hawkes process to their timestamps. Hawkes processes have essentially provided us with a prism through which we can see the market reaction to the different order types and their interaction. More generally, the Hawkes-process approach fits into the longer trend of using point processes to model sequences of irregularly-spaced market events in continuous time (Engle and Russell 1998, Hautsch 2004, Bauwens and Hautsch 2009, Bacry et al. 2013, Bacry and Muzy 2014).

Whilst LOB models based on Hawkes processes have been rather successful in describing the dynamics of the order flow, they are unable to incorporate any endogenous state variables describing the LOB, such as prices, volumes or the bid–ask spread, nor their influence on the arrival rate of orders. We note that Hawkes processes have been used to build full-fledged LOB models (Muni Toke 2011, Abergel and Jedidi 2015), but the arrival rate of orders in these models is not influenced by the state of the LOB. In fact, Hawkes processes are complemented by the other main approach to LOB modelling, based on continuous-time Markov chains (Huang et al. 2015, Huang and Rosenbaum 2017, Cartea et al. 2018a). This alternative approach extends some of the earlier zero-intelligence models (Smith et al. 2003, Cont et al. 2010, Cont and De Larrard 2011) and postulates that the arrival rate of orders is driven by the state of the LOB alone, making the state part of the model. This has the merit of introducing a feedback loop between the order flow and the state of the LOB, but it omits the self- and cross-excitation effects evidenced by the empirical work based on Hawkes processes.

In short, each of these two approaches to LOB modelling has desirable qualities that the other lacks, as also noticed by Bacry et al. (2016), Gonzalez and Schervish (2017), Taranto et al. (2018), and Morariu-Patrichi and Pakkanen (2018). In fact, Gonzalez and Schervish (2017) present empirical results, demonstrating that the type of the next order depends on both the type of the previous order and the state of the LOB, which highlights the need for a more general modelling framework that can amalgamate the two approaches.

The purpose of this paper is to introduce a novel, state-dependent extension of a Hawkes process, in the context of statistical modelling. The new model, formulated in Section 2, pairs a multivariate point process $\mathit{N}$, governed by a past-dependent stochastic intensity, with an observable state process X. The intensity of the point process $\mathit{N}$ is Hawkes process-like, characteristically depending on the past of $\mathit{N}$, but also influenced by the past of X. The process X switches state when there is an event in $\mathit{N}$, according to a Markov transition matrix that depends on the type of the event. This two-way interaction between $\mathit{N}$ and X makes them fully coupled, just like the order flow and the state of the LOB in an order-driven market. It also distinguishes the model from the existing regime switching Hawkes processes (Wang et al. 2012, Cohen and Elliott 2013, Vinkovskaya 2014, Swishchuk and Huffman 2020), where the state process evolves exogenously, receiving no feedback from the point process. Mathematically, we lift the pair $(\mathit{N},X)$ to a higher-dimensional ordinary point process, which allows us to invoke the theory of hybrid marked point processes (Morariu-Patrichi and Pakkanen 2018) to establish the existence and uniqueness of non-explosive solutions and derive a simulation algorithm for the process.

Further, in Section 3 we develop a maximum likelihood (ML) estimation method for parametric specifications of a state-dependent Hawkes process, extending the ML methodology for ordinary Hawkes processes pioneered by Ozaki (1979). We find that the likelihood function of a state-dependent Hawkes process has very convenient separable form that makes it possible to estimate the transition probabilities of X independently of the parameters of the point process $\mathit{N}$, which is remarkable given that $\mathit{N}$ and X are fully coupled. The upshot of this useful property is that, as far as ML estimation is concerned, state-dependent Hawkes processes are no harder to estimate than ordinary multivariate Hawkes processes.

In Section 4 we finally apply the new model and methodology to high-frequency financial data, estimating for the first time an LOB model that accommodates an explicit feedback loop between the state of the LOB and the order flow with self- and cross-excitation effects. More concretely, we estimate two specifications of a state-dependent Hawkes process on level-I LOB data on the stock of Intel (INTC) from Nasdaq, from June 2017 until May 2018, using the bid–ask spread and queue imbalance, respectively, as the state variable. The latter variable can be interpreted as an indicator of the shape of the LOB and is a popular price-predictive signal (Cartea et al. 2018a). Our estimation results reveal that the magnitude and duration of excitation effects in the order flow depend significantly on both state variables, being effective at timescales ranging from 100 microseconds to 100 milliseconds. Since human reaction times have been measured to be of the order of 200 milliseconds, and longer for tasks involving choice (Kosinski 2006), this is a clear sign of the highly algorithmic nature of the order flow in modern electronic markets. Moreover, we find that the high level of endogeneity of the order flow, dubbed critical reflexivity by Filimonov and Sornette (2012) and Hardiman et al. (2013), and quantifiable here by the spectral norm of the Hawkes excitation kernel, is state-dependent and, intriguingly, more pronounced in what can be seen as disequilibrium states of the LOB.

In the Appendix of this paper, we prove the theoretical results of Sections 2 and 3 and gather some auxiliary details on ML estimation. A Python library called mpoints that implements the models and estimation methodology of this paper is available from https://mpoints.readthedocs.io.

2. State-dependent Hawkes processes

2.1. Preliminaries and definition

To set the stage for state-dependent Hawkes processes, we first review some central concepts of point process theory (Brémaud 1981, Daley and Vere-Jones 2003). A $d_e$-dimensional multivariate point process consists of an increasing sequence of positive random times $(T_n{)}_{n\in \mathbb{N}}$ and a matching sequence of random marks $(E_n{)}_{n\in \mathbb{N}}$ in $\mathcal{E}:=\{1,\dots ,d_e\}$. We interpret, for any $n\in \mathbb{N}$, the pair $(T_n,E_n)$ as an event of type $E_n$ occurring at time $T_n$. A point process can be identified by its counting process representation $\mathit{N}:=(N_1,\dots ,N_{d_e})$, where $N_e\left(t\right):=\sum _{n\in \mathbb{N}}{\text{1}}_{\{T_n\le t,E_n=e\}},t\ge 0,e\in \mathcal{E},$ counts how many events of type e have occurred by time t. The counting process $\mathit{N}$ is said to be non-explosive if $\underset{n\to \mathrm{\infty }}{lim}{T}_n=\mathrm{\infty }$ with probability one. Given a filtration $\mathbb{F}=({\mathcal{F}}_t{)}_{t\ge 0}$ to which $\mathit{N}$ is adapted, we say that a non-negative $\mathbb{F}$-predictable process $\mathit{\lambda }=({\lambda }_1,\dots ,{\lambda }_{d_e})$ is the $\mathbb{F}$-intensity of $\mathit{N}$ if $\begin{array}{rl}& \mathbb{E}\left[N_e\right(t)-N_e(s\left)|{\mathcal{F}}_s\right]=\mathbb{E}\left[{\int }_s^t{\lambda }_e\right(u\left)\mathrm{d}u|{\mathcal{F}}_s\right],t\ge s\ge 0,\\ & e\in \mathcal{E}.\end{array}$ Intuitively, ${\lambda }_e\left(t\right)$ is the infinitesimal rate of new events of type e at time t.

To construct a state-dependent Hawkes process, we couple the point process $\mathit{N}$ with a càdlàg state process $\left(X\right(t){)}_{t\ge 0}$ that takes values in a finite state space $\mathcal{X}:=\{1,\dots ,d_x\}$ and denote by ${\mathbb{F}}^{\mathit{N},X}$ the natural filtration of the pair $(\mathit{N},X)$.

Definition 2.1

Let $\mathit{\varphi }=({\mathit{\varphi }}_e{)}_{e\in \mathcal{E}}$ be a collection of $d_x\times d_x$ transition probability matrices, $\mathit{\nu }=({\nu }_1,\dots ,{\nu }_{d_e})\in {\mathbb{R}}_{>0}^{d_e}$ and $\mathit{k}=(k_{e^{\prime }e}{)}_{e^{\prime },e\in \mathcal{E}}$, where $k_{e^{\prime }e}:{\mathbb{R}}_{>0}\times \mathcal{X}\to {\mathbb{R}}_{\ge 0}$. The pair $(\mathit{N},X)$ is a state-dependent Hawkes process with transition distribution $\mathit{\varphi }$, base rate vector $\mathit{\nu }$ and excitation kernel $\mathit{k}$, abbreviated as $\text{sdHawkes}(\mathit{\varphi },\mathit{\nu },\mathit{k})$, if

1. $\mathit{N}$ has ${\mathbb{F}}^{\mathit{N},X}$-intensity $\mathit{\lambda }$ that satisfies (1) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(s\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E};\end{array}$(1)

2. X is piecewise constant and jumps only at the event times $(T_n{)}_{n\in \mathbb{N}}$, so that (2) $\begin{array}{rl}& \mathbb{P}\left(X\right(T_n)=x|E_n,{\mathcal{F}}_{T_n-}^{\mathit{N},X})={\varphi }_{E_n}\left(X\right(T_n-),x),\\ & n\in \mathbb{N},x\in \mathcal{X},\end{array}$(2)

where ${\mathcal{F}}_{T_n-}^{\mathit{N},X}:=\underset{ϵ>0}{\bigvee }{\mathcal{F}}_{T_n-ϵ}^{\mathit{N},X}$ and $X(T_n-):=\underset{t\uparrow T_n}{lim}X\left(t\right)$.

In applications to LOB modelling, the counting process $\mathit{N}$ and state process X will represent the order flow and the state of the LOB, respectively. An event of type $e^{\prime }\in \mathcal{E}$ at time $T=T_n$, for some $n\in \mathbb{N}$, increases the infinitesimal rate of new events of type $e\in \mathcal{E}$ at time t>T by $k_{e^{\prime }e}(t-T,X(T\left)\right)$. The novelty here is that the self- and cross-excitation effects now depend on the state process X. Reciprocally, at each event in $\mathit{N}$, the state process X may switch to a new state according to the probability (2(2) $\begin{array}{rl}& \mathbb{P}\left(X\right(T_n)=x|E_n,{\mathcal{F}}_{T_n-}^{\mathit{N},X})={\varphi }_{E_n}\left(X\right(T_n-),x),\\ & n\in \mathbb{N},x\in \mathcal{X},\end{array}$(2) ) that depends on the event type. Thereby, mimicking the mechanics of the LOB, the processes $\mathit{N}$ and X are fully coupled. In the case $d_x=1$, Definition 2.1 reduces to that of an ordinary linear Hawkes process. A simulated sample path of a state-dependent Hawkes process with $d_e=1$, $d_x=2$ and exponential kernel is shown in Figure 1. In this example, the process exhibits self-excitation only in the second state.

Figure 1. Simulation of a state-dependent Hawkes process with $d_e=1$, $d_x=2$. The upper plot shows the evolution of the state process. The blue dots indicate the event times and the lower plot represents the intensity. The process is specified so that ${\nu }_1=1$ and $k_{11}(t,x)=\exp (-4t){\text{1}}_{\{x=2\}}$, that is, in state 2 the process exhibits exponential self-excitation whereas no self-excitation occurs in state 1.

Remark 2.2

1. The state space $\mathcal{X}$ is chosen to be finite here on the grounds of practical estimation. Theoretically, a state-dependent Hawkes process can however be defined in more general (infinite) state spaces (Morariu-Patrichi and Pakkanen 2018, Example 2.15).

2. The relationship (2(2) $\begin{array}{rl}& \mathbb{P}\left(X\right(T_n)=x|E_n,{\mathcal{F}}_{T_n-}^{\mathit{N},X})={\varphi }_{E_n}\left(X\right(T_n-),x),\\ & n\in \mathbb{N},x\in \mathcal{X},\end{array}$(2) ) does not imply that X is a pure jump Markov process. In particular, the time to next state transition need not be conditionally exponentially distributed, since it is determined by the counting process $\mathit{N}$, which nests a multitude of non-Markovian ordinary Hawkes processes, for instance.

3. If for any $e\in \mathcal{E}$ the transition distribution $\left({\varphi }_e\right(x^{\prime },x){)}_{x\in \mathcal{X}}$ does not depend on the previous state $x^{\prime }\in \mathcal{X}$, the process $\mathit{N}$ reduces to a marked Hawkes process, as described in Bacry et al. (2015, Section 2.2.1).

2.2. Existence and uniqueness of non-explosive solutions

Since Definition 2.1 is implicit in the sense that the counting process $\mathit{N}$ is defined by intensity $\mathit{\lambda }$ that depends on the past of $\mathit{N}$ and X, care is needed to establish the existence and uniqueness of a pair $(\mathit{N},X)$ that solves (1(1) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(s\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E};\end{array}$(1) ) and (2(2) $\begin{array}{rl}& \mathbb{P}\left(X\right(T_n)=x|E_n,{\mathcal{F}}_{T_n-}^{\mathit{N},X})={\varphi }_{E_n}\left(X\right(T_n-),x),\\ & n\in \mathbb{N},x\in \mathcal{X},\end{array}$(2) ) so that $\mathit{N}$ is non-explosive. To this end, we lift $(\mathit{N},X)$ to a $d_e{d}_x$-dimensional multivariate point process $\stackrel{\mathbf{~}}{\mathit{N}}=({\tilde{N}}_{ex}{)}_{e\in \mathcal{E},x\in \mathcal{X}}$, where ${\tilde{N}}_{ex}$ counts the number of events of type e after which the state is x, formally, ${\tilde{N}}_{ex}\left(t\right)=\sum _{n\in \mathbb{N}}{\text{1}}_{\{T_n\le t,E_n=e,X(T_n)=x\}},t\ge 0,e\in \mathcal{E},x\in \mathcal{X}.$ The marks corresponding to $\stackrel{\mathbf{~}}{\mathit{N}}$ are given by $(E_n,X_n)$, $n\in \mathbb{N}$, where $E_n$ is as above and $X_n:=X\left(T_n\right)$ is the value of the state process following the nth event and $X_0$ is the initial state. Thus, given $\stackrel{\mathbf{~}}{\mathit{N}}$, the state $X\left(t\right)$ can be recovered from the most recent mark at time t, which makes the relationship between $\stackrel{\mathbf{~}}{\mathit{N}}$ and $(\mathit{N},X)$ bijective. In fact, state-dependent Hawkes processes were first introduced using this representation in Morariu-Patrichi and Pakkanen (2018).

By applying the general characterisation result in Morariu-Patrichi and Pakkanen (2018), the dynamics of $\stackrel{\mathbf{~}}{\mathit{N}}$ can be expressed in terms of the dynamics of $(\mathit{N},X)$, and vice versa. The natural filtration of $\stackrel{\mathbf{~}}{\mathit{N}}$ is denoted by ${\mathbb{F}}^{\stackrel{\mathbf{~}}{\mathit{N}}}$.

Theorem 2.3

The pair $(\mathit{N},X)$ is a non-explosive $\text{sdHawkes}(\mathit{\varphi },\mathit{\nu },\mathit{k})$ process if and only if $\stackrel{\mathbf{~}}{\mathit{N}}$ is non-explosive, admitting ${\mathbb{F}}^{\stackrel{\mathbf{~}}{\mathit{N}}}$-intensity $\stackrel{\mathbf{~}}{\mathit{\lambda }}$ that satisfies (3) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\varphi }_e\left(X\right(t),x)\\ & \times ({\nu }_e+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,x^{\prime }\left)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\right(s\left)\right),\\ & t\ge 0,e\in \mathcal{E},x\in \mathcal{X}.\end{array}$(3)

The above theorem in fact shows that state-dependent Hawkes processes belong to the class of hybrid marked point processes studied in Morariu-Patrichi and Pakkanen (2018). The general existence and uniqueness results therein apply to the present class of processes as follows.

Theorem 2.4

A unique, non-explosive $\text{sdHawkes}(\mathit{\varphi },\mathit{\nu },\mathit{k})$ process exists if one of the following two conditions is satisfied:

1. the components of $\mathit{k}$ are bounded functions;

2. $\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_0^{\mathrm{\infty }}{k}_{e^{\prime }e}(t,x^{\prime })dt<\left(\underset{x^{\prime }\in \mathcal{X}}{max}{\varphi }_e\right(x^{\prime },x){)}^{-1}$ for all $e\in \mathcal{E}$ and $x\in \mathcal{X}$.

Condition (i) above suffices for the purposes of the present paper, where we apply bounded, exponential kernels. However, condition (ii) is included here for completeness as it reduces to the classical stability condition of Massoulié (1998) in the case $d_x=1$.

2.3. Simulation

Another implication of Theorem 2.3 is that the simulation of a state-dependent Hawkes process can be reduced to the simulation of a multivariate point process with an intensity given by (3(3) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\varphi }_e\left(X\right(t),x)\\ & \times ({\nu }_e+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,x^{\prime }\left)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\right(s\left)\right),\\ & t\ge 0,e\in \mathcal{E},x\in \mathcal{X}.\end{array}$(3) ) and, thus, many simulation techniques from point process theory can be reused (Lewis and Shedler 1976, Daley and Vere-Jones 2003).

In fact, the sample path in Figure 1 has been generated using Ogata's thinning algorithm (Ogata 1981), which is an exact simulation algorithm and is adaptable for state-dependent Hawkes processes as follows. We write $R\left(t\right):=\sum _{e\in \mathcal{E},x\in \mathcal{X}}{\tilde{\lambda }}_{ex}\left(t\right)=\sum _{e\in \mathcal{E}}{\lambda }_e\left(t\right)$, which is a function of all $(T_n,E_n,X_n)$ such that $T_n<t$.

Remark 2.6

1. In Algorithm 2.5, we are implicitly assuming that the components of the kernel $\mathit{k}$ are non-increasing, which guarantees that $R\left(t^{\prime }\right)\le R\left(t\right)$ for all $t^{\prime }\in [t,T_n]$. In the general case, one needs to define $R\left(t\right)$ so that it bounds the total intensity $\sum _{e\in \mathcal{E}}{\lambda }_e\left(t^{\prime }\right)$ for all $t^{\prime }\in [t,T_n]$.

2. Lines 9–10 of Algorithm 2.5 use the product form (3(3) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\varphi }_e\left(X\right(t),x)\\ & \times ({\nu }_e+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,x^{\prime }\left)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\right(s\left)\right),\\ & t\ge 0,e\in \mathcal{E},x\in \mathcal{X}.\end{array}$(3) ) to simulate the marks, which avoids the computation of $d_e{d}_x$ products at the cost of generating an additional random number.

2.4. Comparison with related models of limit order books

We will now briefly compare state-dependent Hawkes processes, in the sense of Definition 2.1, to some recent, closely related models, which similarly aim to couple a point process to a state process in the context of LOB modelling.

The regime-switching Hawkes process of Vinkovskaya (2014) can be seen as a special case of Definition 2.1, with the exception that the dynamics of the state process X are not modelled. (Effectively, this means that the counting process $\mathit{N}$ is then specified conditional on a realisation of X, precluding two-way interaction between $\mathit{N}$ and X.) In her empirical application, Vinkovskaya (2014) estimates a four-dimensional version of the model for arriving orders of four different types in the New York Stock Exchange Trades and Quotes (TAQ) data. She employs the bid–ask spread as a state variable, like we do in one of our models in Section 4.

Swishchuk and Huffman (2020) construct a compound Hawkes process using a one-dimensional ordinary Hawkes process N (possibly non-linear) and a Markov chain $(Y_n{)}_{n\in \mathbb{N}}$ in a finite state space $\mathcal{Y}$, independent of N. The actual process $\left(S\right(t){)}_{t\ge 0}$, which is used as a model of prices, is given by $S\left(t\right):=S\left(0\right)+\sum _{n=1}^{N\left(t\right)}a\left(Y_n\right),t\ge 0,$ where a is a function from $\mathcal{Y}$ to $\mathbb{R}$. Whilst having similar ingredients, the compound Hawkes process is otherwise unrelated to state-dependent Hawkes processes. In particular, the Markov chain $(Y_n{)}_{n\in \mathbb{N}}$ in the compound Hawkes process does not influence the rate of events, that is, price changes. Swishchuk and Huffman (2020) establish limit theorems for the process, describing its long-term behaviour, and also empirically estimate it using Nasdaq Stock Market LOB data.

Note that in Definition 2.1 the excitation effect of each event is predicated upon the state prevailing at the time the event occurs. Alternatively, we could make the level of excitation track the current state. We could also make the base rates track the current state. Applying these modifications to equation (1(1) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(s\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E};\end{array}$(1) ) in Definition 2.1 yields an analogous intensity satisfying (4) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e\left(X\right(t-\left)\right)+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(t-\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E}.\end{array}$(4) Whether it is ultimately more natural to use $X\left(s\right)$ or $X(t-)$ in the excitation kernel is open to debate, but in the current context of LOB modelling, where the most significant excitation effects tend to be ephemeral, the difference is unlikely to be large in practice.

Now, several recent LOB models conform to (4(4) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e\left(X\right(t-\left)\right)+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(t-\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E}.\end{array}$(4) ). Cohen and Elliott (2013) introduce a one-dimensional Markov-modulated Hawkes process following (4(4) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e\left(X\right(t-\left)\right)+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(t-\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E}.\end{array}$(4) ), where X is an exogenous Markov process in a finite state space (exogenous in the sense that the counting process does not influence X). The key feature of their work is that they assume X to be unobservable, leading them to derive a filtering procedure for the estimation of the current state of X. Cohen and Elliott (2013) illustrate their methodology by estimating a regime switch in TAQ trade data during the US equity market flash crash on 6 May 2010.

Coinciding with the first preprint version of the present paper, Wu et al. (2019) develop a queue-reactive Hawkes process based on (4(4) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e\left(X\right(t-\left)\right)+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(t-\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E}.\end{array}$(4) ). In their model, X is endogenous and carries information about queue lengths in the LOB, while the multi-dimensional counting process driven by the intensity (4(4) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e\left(X\right(t-\left)\right)+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(t-\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E}.\end{array}$(4) ) models events pertaining to these queues. Wu et al. (2019) estimate their model on German bond (Bund) and index (DAX) futures LOB data. Subsequently, Mounjid et al. (2019) generalise the queue-reactive Hawkes process to a more general point process framework that allows for non-linearity and quadratic Hawkes structure. Mounjid et al. (2019) additionally establish ergodicity for the model and also derive functional limit theorems for its long-term behaviour. They apply the model to evaluate and rank equities market makers on Euronext Paris.

Finally, Fosset et al. (2020) develop a Hawkes process model of endogenous liquidity crises. They model the increase and decrease, respectively, in the bid–ask spread as a two-dimensional point process that conforms to (4(4) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e\left(X\right(t-\left)\right)+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(t-\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E}.\end{array}$(4) ). The first component of the process follows an ordinary Hawkes process (possibly non-linear), whilst the second component has a state-dependent base rate, depending on the the current bid–ask spread. Fosset et al. (2020) then analyse the long-term stability of the model.

3. Parametric estimation via maximum likelihood

Estimating the base rate vector $\mathit{\nu }$ and kernel $\mathit{k}$ of an ordinary Hawkes process has become a vibrant research topic in statistics and statistical finance literature. Whilst the recent focus has mostly been on non-parametric methodology (Bacry and Muzy 2016, Kirchner 2017, Eichler et al. 2017, Sancetta 2018, Achab et al. 2018), our aim in this paper is to extend the classical parametric framework (Ozaki 1979), comprehensively summarised in Bowsher (2007), to state-dependent Hawkes processes, leaving non-parametric methodology for future work. From now on, we work with a kernel $\mathit{k}={\mathit{k}}^{\mathit{\theta }}$ parametrised by a vector $\mathit{\theta }\in {\mathbb{R}}^p$.

3.1. Likelihood function

We know from Subsections 2.2 and 2.3 that a state-dependent Hawkes process $(\mathit{N},X)$ can be lifted to a $d_e{d}_x$-dimensional point process $\stackrel{\mathbf{~}}{\mathit{N}}$. Given a realisation $(t_n,e_n,x_n{)}_{n=1,\dots ,N}$ of $\stackrel{\mathbf{~}}{\mathit{N}}$ over a time horizon $[0,T]$, the likelihood function $\mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })$ can be informally understood as the probability that $\sum _{e\in \mathcal{E},x\in \mathcal{X}}{\tilde{N}}_{ex}\left(T\right)=N$ and $(T_n,E_n,X_n{)}_{n=1,\dots ,N}$ lies in a small neighbourhood of $(t_n,e_n,x_n{)}_{n=1,\dots ,N}$, under the assumption that $(T_n,E_n,X_n{)}_{n\in \mathbb{N}}$ is generated by a state-dependent Hawkes processes with parameters $(\mathit{\varphi },\mathit{\nu },\mathit{\theta })$. More rigorously, the likelihood function is the density of the Janossy measure with respect to the Lebesgue measure on ${\mathbb{R}}^N$ (Daley and Vere-Jones 2003, p. 125, 213).

For ordinary Hawkes processes, the likelihood function $\mathcal{L}(\mathit{\nu },\mathit{\theta })$ can be expressed directly in terms of $\mathit{\lambda }$ (Daley and Vere-Jones 2003) and the maximum likelihood (ML) estimator $(\hat{\mathit{\nu }},\hat{\mathit{\theta }})$ is obtained by maximising $\mathcal{L}(\hat{\mathit{\nu }},\hat{\mathit{\theta }})$, in practice numerically. For state-dependent Hawkes processes, we are able to express $\mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })$ in terms of $\mathit{\varphi }$ and $\mathit{\lambda }$, and find that maximising the likelihood is conveniently achieved by solving two independent optimisation problems.

Theorem 3.1

The log likelihood function of an $\text{sdHawkes}(\mathit{\varphi },\mathit{\nu },{\mathit{k}}^{\mathit{\theta }})$ process is given by (5) $\begin{array}{rl}\ln \mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })& =\sum _{n=1}^N\ln {\varphi }_{e_n}(x_{n-1},x_n)+\sum _{n=1}^n\ln {\lambda }_{e_n}\left(t_n\right)\\ & -{\int }_0^T\sum _{e\in \mathcal{E}}{\lambda }_e\left(t\right)\mathrm{d}t.\end{array}$(5) Furthermore, $(\hat{\mathit{\varphi }},\hat{\mathit{\nu }},\hat{\mathit{\theta }})\in {\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_{\mathit{\varphi },\mathit{\nu },\mathit{\theta }}\mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })$ if and only if $\begin{array}{r}\{\begin{array}{l}{\hat{\varphi }}_e(x,x^{\prime })={\displaystyle \frac{\sum _{n=1}^N\text{1}(x_{n-1}=x,e_n=e,x_n=x^{\prime })}{\sum _{n=1}^N\text{1}(x_{n-1}=x,e_n=e)}},e\in \mathcal{E},\\ x,x^{\prime }\in \mathcal{X},\\ {\displaystyle (\hat{\mathit{\nu }},\hat{\mathit{\theta }})\in {\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_{\mathit{\nu },\mathit{\theta }}\sum _{n=1}^N\ln {\lambda }_{e_n}\left(t_n\right)-{\int }_0^T\sum _{e\in \mathcal{E}}{\lambda }_e\left(t\right)\mathrm{d}t.}\end{array}\end{array}$

The upshot of Theorem 3.1 is that the ML estimation of state-dependent Hawkes processes is no harder than that of ordinary Hawkes processes. Namely, $\mathit{\varphi }$ is estimated in a straightforward manner by the empirical transition probabilities, whilst $(\mathit{\nu },\mathit{\theta })$ is estimated by maximising the log quasi-likelihood of $\mathit{N}$, that is, the Radon-Nikodym derivative of a change of measure that transforms a standard Poisson process into $\mathit{N}$ (Brémaud 1981, Theorem 10, p. 241), which is similar to the log likelihood of a multivariate ordinary Hawkes process. It is remarkable that, in spite of the strong coupling between the events and the state process, the estimation of $\mathit{\varphi }$ and $(\mathit{\nu },\mathit{\theta })$ is decoupled due to the separable form (5(5) $\begin{array}{rl}\ln \mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })& =\sum _{n=1}^N\ln {\varphi }_{e_n}(x_{n-1},x_n)+\sum _{n=1}^n\ln {\lambda }_{e_n}\left(t_n\right)\\ & -{\int }_0^T\sum _{e\in \mathcal{E}}{\lambda }_e\left(t\right)\mathrm{d}t.\end{array}$(5) ) of the log likelihood function.

Remark 3.2

It should be stressed that the separable form (5(5) $\begin{array}{rl}\ln \mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })& =\sum _{n=1}^N\ln {\varphi }_{e_n}(x_{n-1},x_n)+\sum _{n=1}^n\ln {\lambda }_{e_n}\left(t_n\right)\\ & -{\int }_0^T\sum _{e\in \mathcal{E}}{\lambda }_e\left(t\right)\mathrm{d}t.\end{array}$(5) ) of the log likelihood function $\mathcal{L}$ is not a trivial consequence of Definition 2.1. It is once again the lift $\stackrel{\mathbf{~}}{\mathit{N}}$ that allow us to transpose the problem to the setting of point process theory, where classical results apply (see the proof of Theorem 3.1 in Appendix A.1).

We note that, in the case of ordinary Hawkes processes ($d_x=1$), consistency and asymptotic normality results for the ML estimator $(\hat{\mathit{\nu }},\hat{\mathit{\theta }})$ are available in the literature, see for example (Ogata 1978, Clinet and Yoshida 2016). However, these results rely on the stationarity and ergodicity of the underlying process and, unfortunately, these properties need not extend to state-dependent Hawkes processes in general. If the excitation kernel $\mathit{k}$ is dominated in all states by a non-state-dependent kernel that satisfies a classical stability condition (Massoulié 1998), we conjecture that then stationarity and ergodicity hold. This can potentially be shown by adapting the arguments used in Mounjid et al. (2019) to prove ergodicity for closely related point processes, or in the special case of exponential kernels (see Subsection 3.3) using Markov process techniques, as outlined in Section 4.7 in Morariu-Patrichi (2019). The case where $\mathit{k}$ fails to be dominated by a stable kernel in some of the states appears however to be empirically more relevant given our results in Subsection 4.8. In this case, stationarity and ergodicity become more elusive, whilst they may still plausibly hold as long as the sojourns of the process in the ‘unstable states’ do not dominate its time evolution. All in all, the analysis of the asymptotic properties of the ML estimator is unfortunately not straightforward and remains beyond the scope of the present paper.

In Appendix A.2.3, we present Monte Carlo results that exemplify the favourable finite-sample performance of the ML estimator. In particular, the results provide evidence that $(\hat{\mathit{\varphi }},\hat{\mathit{\nu }},\hat{\mathit{\theta }})$ is consistent as the length of the estimation window increases.

3.2. Goodness-of-fit diagnostics using residuals

The goodness of fit of an estimated $\text{sdHawkes}(\hat{\mathit{\varphi }},\hat{\mathit{\nu }},{\mathit{k}}^{\hat{\mathit{\theta }}})$ process can be assessed as follows. Denote by $(T_n^e{)}_{n=1,\dots ,n_e}$ the sequence of times at which an event of type e occurred and by $(T_n^{ex}{)}_{n=1,\dots ,n_{ex}}$ the sequence of times at which an event type e occurred and after which the state was x, and set $T_0^e=T_0^{ex}=0$ for all $e\in \mathcal{E}$ and $x\in \mathcal{X}$. Introduce additionally the event residuals $r_n^e$ and total residuals ${\tilde{r}}_n^{ex}$ $\begin{array}{rl}{r}_n^e& :={\int }_{T_{n-1}^e}^{T_n^e}{\hat{\lambda }}_e\left(t\right)\mathrm{d}t,n=1,\dots ,n_e,e\in \mathcal{E},\\ {\tilde{r}}_n^{ex}& :={\int }_{T_{n-1}^{ex}}^{T_n^{ex}}{\widehat{\stackrel{~}{\lambda }}}_{ex}\left(t\right)\mathrm{d}t,n=1,\dots ,n_{ex},e\in \mathcal{E},x\in \mathcal{X},\end{array}$ where $\hat{\mathit{\lambda }}$ and $\widehat{\stackrel{~}{\mathit{\lambda }}}$ are computed by plugging the estimated values $(\hat{\mathit{\varphi }},\hat{\mathit{\nu }},{\mathit{k}}^{\hat{\mathit{\theta }}})$ into (1(1) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(s\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E};\end{array}$(1) ) and (3(3) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\varphi }_e\left(X\right(t),x)\\ & \times ({\nu }_e+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,x^{\prime }\left)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\right(s\left)\right),\\ & t\ge 0,e\in \mathcal{E},x\in \mathcal{X}.\end{array}$(3) ), respectively.

It is a classical result that under the right time change, $\mathit{N}$ and $\stackrel{\mathbf{~}}{\mathit{N}}$ become unit-rate Poisson processes (Meyer 1971), with the event residuals and total residuals as their respective time increments, which applies to the present setting as follows.

Theorem 3.3

1. Suppose that $(T_n^e{)}_{n=1,\dots ,n_e}$ are generated by a $d_e$-dimensional multivariate point process $\mathit{N}$ with an ${\mathbb{F}}^{\mathit{N},X}$-intensity $\mathit{\lambda }$ satisfying (1(1) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(s\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E};\end{array}$(1) ), where X is a given state process, $\mathit{\nu }=\hat{\mathit{\nu }}$ and $\mathit{k}={\mathit{k}}^{\hat{\mathit{\theta }}}$. Then the event residuals $(r_n^e{)}_{n=1,\dots ,n_e}$ for any $e\in \mathcal{E}$ are i.i.d. and follow the $\text{Exp}\left(1\right)$ distribution.

2. Suppose that $(T_n^{ex}{)}_{n=1,\dots ,n_{ex}}$ are generated by an $\text{sdHawkes}(\hat{\mathit{\varphi }},\hat{\mathit{\nu }},{\mathit{k}}^{\hat{\mathit{\theta }}})$ process $(\mathit{N},X)$. If moreover ${\tilde{N}}_{ex}\left(t\right)\to \mathrm{\infty }$ as $t\to \mathrm{\infty }$ with probability one for any $e\in \mathcal{E}$ and $x\in \mathcal{X}$, then the total residuals $({\tilde{r}}_n^{ex}{)}_{n=1,\dots ,n_{ex}}$ for any $e\in \mathcal{E}$ and $x\in \mathcal{X}$ are i.i.d. and follow the $\text{Exp}\left(1\right)$ distribution.

Consequently, the goodness-of-fit of the process can be assessed by comparing the empirical distribution of each sequence among $(r_n^e{)}_{n=1,\dots ,n_e}$, $e\in \mathcal{E}$ and $({\tilde{r}}_n^{ex}{)}_{n=1,\dots ,n_{ex}}$, $e\in \mathcal{E}$, $x\in \mathcal{X}$ to the $\text{Exp}\left(1\right)$ distribution, which can by achieved by inspecting a Q–Q plot (see Subsection 4.7). Goodness-of-fit diagnostics can also check whether the residuals are mutually independent, which can be assessed by examining their correlogram. One can of course go beyond these visual assessments by applying formal statistical tests (Bowsher 2007).

3.3. The special case of exponential kernels

The direct computation of the term $\sum _{n=1}^N\ln {\lambda }_{e_n}\left(t_n\right)$ in the log likelihood function (5(5) $\begin{array}{rl}\ln \mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })& =\sum _{n=1}^N\ln {\varphi }_{e_n}(x_{n-1},x_n)+\sum _{n=1}^n\ln {\lambda }_{e_n}\left(t_n\right)\\ & -{\int }_0^T\sum _{e\in \mathcal{E}}{\lambda }_e\left(t\right)\mathrm{d}t.\end{array}$(5) ) involves a double sum requiring $\mathcal{O}\left(N^2\right)$ operations to evaluate, which may render ML estimation numerically infeasible when the sample is large. However, for ordinary Hawkes processes ($d_x=1$) with exponential kernels, it is known that this computation can be achieved in $\mathcal{O}\left(N{d}_e^2\right)$ operations (Ozaki 1979, Ogata 1981). Fortunately, this property carries over to state-dependent Hawkes processes with kernel $\mathit{k}$ having exponential form (6) $k_{e^{\prime }e}(t,x)={\alpha }_{e^{\prime }xe}\exp (-{\text{β}}_{e^{\prime }xe}t),t>0,e^{\prime },e\in \mathcal{E},x\in \mathcal{X},$(6) where the impact coefficients $\mathit{\alpha }:=({\alpha }_{e^{\prime }xe}{)}_{}$ and decay coefficients $\text{β}:=({\text{β}}_{e^{\prime }xe}{)}_{}$ are non-negative. The reduction in computational cost becomes apparent from the derivation in Appendix A.2.1. In particular, the term $\sum _{n=1}^N\ln {\lambda }_{e_n}\left(t_n\right)$ can be expressed in terms of sums that satisfy a convenient recursive relationship, reducing the computational cost from $\mathcal{O}\left(N^2\right)$ to $\mathcal{O}\left(N\right)$ operations. A similar remark holds for the computation of the gradient and the residuals. Moreover, the $\mathcal{O}\left(N\right)$ complexity extends also to any kernel whose components are linear combinations of functions of the form (6(6) $k_{e^{\prime }e}(t,x)={\alpha }_{e^{\prime }xe}\exp (-{\text{β}}_{e^{\prime }xe}t),t>0,e^{\prime },e\in \mathcal{E},x\in \mathcal{X},$(6) ).

4. Application to high-frequency limit order book data

4.1. Limit order book mechanism

We first briefly review the mechanics of a limit order book and recall the definitions of some key market quantities, following Gould et al. (2013).

In order-driven markets, market participants submit orders to buy or sell an asset (e.g. a stock) at the price of their choice. Formally, an order is defined by its submission time, direction (buy or sell), price p and size q. The quantities p and q must typically be multiples of the tick size and lot size, respectively, which are fixed by the exchange or the market regulator. When a buy order is submitted, if there are unfilled sell orders with prices $p^{\prime }<p$, the buy order is matched at the smallest price $p^{\prime }$ with the oldest order(s), following the price–time priority rule. Otherwise the order becomes active, that is, it enters the queue of unfilled orders with price p. Such orders are called limit orders. The collection of current limit orders is called the limit order book (LOB) and can thus be understood as a snapshot of the expressed supply and demand or visible liquidity. Orders that result in an instant match with pre-existing limit orders are called market orders. Orders may also be cancelled, that is, an unfilled or partially filled order is withdrawn from the LOB.

The highest (respectively lowest) price among buy (respectively sell) active limit orders is called the bid (respectively ask) price. For instance, the bid price is the best price at which one can instantly sell by sending a market order. The difference between the ask price and bid price is called the bid–ask spread. The orders and cancellations with prices that are equal to or more aggressive than the bid and ask prices at their submission time form the level-I order flow.

4.2. Data

We analyse tick-by-tick level-I LOB data on the stock of Intel Corporation (INTC), traded on the Nasdaq Stock Market, from 1 June 2017 to 31 May 2018. We discard the data on 3 June 2017 (the day before Independence Day) and 24 November 2017 (Black Friday, i.e. the day after Thanksgiving) because of early market close on those days at 13:00, leaving 250 full trading days in the data set. Besides INTC, we have also studied in the same manner the stocks of Advanced Micro Devices Inc. (AMD), Micron Technology Inc. (MU), Snap Inc. (SNAP) and Twitter Inc. (TWTR) from January to April 2018, but due to space constraints we only report our results on INTC as they are largely representative of the findings on the other stocks. (Full results are available from the authors upon request.)

Whilst the primary listing of INTC is the Nasdaq Stock Market, it can be traded on several stock exchanges and alternative trading systems (ATS) across the increasingly fragmented US equity market system (O'Hara and Ye 2011). In particular, there is no single LOB that would in real time aggregate the available liquidity on all trading venues—the Nasdaq LOB of INTC represents only a part of the visible liquidity. The LOBs of different exchanges cannot, however, diverge significantly, at least in terms of prices, due to the arbitrage opportunities that would ensue and the National Best Bid and Offer (NBBO) rule that mandates brokers to execute trades at the best available prices in the market system. Moreover, we deem the Nasdaq LOB of INTC to be representative of the state of the market since it had the largest market share of INTC among all exchanges during the observation period—according to Nasdaq trading volume statistics (https://www.nasdaqtrader.com/trader.aspx?ID=marketsharedailyhttps://www.nasdaqtrader.com/trader.aspx?ID=marketsharedaily), around 30% of the total INTC market volume was traded on Nasdaq, which is a typical figure for a stock with Nasdaq as primary listing.

In the US equity markets, the tick size is fixed to $0.01 by Rule 612 of Regulation National Market System (Reg NMS), with the exception of stocks priced below $1.00 per share. This means that for a stock with a low price per share, the uniform tick size is relatively large compared to the price. Such stocks are dubbed large-tick stocks. Whilst the exact criterion is subjective, the price of INTC during the sample period was until January 2018 below the threshold of $50 per share used by Bonart and Gould (2017) to distinguish large-tick stocks and below $60 during the entire period. A characteristic feature of large-tick stocks is that their bid–ask spread is most of the time equal to one tick, which is confirmed for INTC in Figure 3(a). Another important feature of large-tick stocks is that most of their liquidity and trading activity is concentrated on the bid and ask levels, which is also our rationale for focusing on level-I data on INTC and eschewing deeper levels of the LOB.

In the modelling part of our analysis, we focus exclusively on trading activity between 12:00 and 14:30 for the following reason. It is well-known, and also exemplified in Figure 2(a), that the intensity of trading activity is not constant throughout the day but follows on average a U-shaped curve. Because of this diurnal pattern, imposing a constant base rate vector $\mathit{\nu }$ over the entire trading day might result in overestimation of the self- and cross-excitation effects (Rambaldi et al. 2015, Omi et al. 2017). The intensity of trading activity is relatively constant over the chosen intraday period, whilst the choice still leaves an ample amount of observations for model estimation (at least 50,000 level-I orders on each day, see Figure 2(b)).

Figure 2. Descriptive statistics of level-I order flow of INTC. Except for Figure 2(a), only the data between 12:00 and 14:30 are used. In Figure 2(a), the sample mean of the arrival rate of level-I orders is computed over 10-minute bins. The translucent area represents the range of the arrival rate across the 250 trading days, excluding the bottom and top 5% values. (a) Arrival rate of level-I orders. (b) Number of level-I orders. (c) Fraction of level-I orders with non-unique timestamp. (d) Distribution of level-I order types.

Figure 3. Joint distribution of events and states for INTC, depicting the empirical distribution of the marks $(E_n,X_n)$ for the two considered state processes. (a) ${\text{Model}}_{\text{S}}$ (state process: bid–ask spread). (b) ${\text{Model}}_{\text{QI}}$ (state process: queue imbalance).

The data was supplied by LOBSTER (https://lobsterdata.comhttps://lobsterdata.com) in a form where the LOB at the time of each event has already been reconstructed. The timestamp of each event is recorded with nanosecond precision. As a result, more than 99% of level-I orders have a unique timestamp (Figure 2(c)). This contrasts with some earlier studies (Bowsher 2007, Large 2007) where lower timestamp resolutions (e.g. one second) were used, leading to a significant amount of tied timestamps, shared by multiple events. As Hawkes processes capture the Granger causality between different event types (Eichler et al. 2017, Embrechts and Kirchner 2018), being able to accurately establish the order of events, even at the shortest timescales, is essential for accurate estimation of Hawkes processes. It is worth pointing out that in LOBSTER data, a market order that is matched with, say, n multiple limit orders is recorded as n individual market orders sharing the same timestamp. However, in our analysis, including Figure 2(c), we address this artefact by aggregating market orders with a tied timestamp and counting them as a single level-I event.

4.3. Model specification

We work with models, based on state-dependent Hawkes processes, that distinguish two event types, denoted ask and bid (that is, $\mathcal{E}=\left\{ask,bid\right\}$, $d_e=2$). In our analysis, bid events consist of buy market orders, level-I buy limit orders and level-I sell cancellations, whilst ask events inversely consist of sell market orders, level-I sell limit orders and level-I buy cancellations. Consequently, $N_{\mathsf{\text{{ bid}}}}-N_{\mathsf{\text{{ ask}}}}$ can be interpreted as a proxy of the order flow imbalance, which was shown by Cont et al. (2013) to be the main driver of price changes. More concretely, bid events tend to push the price up and ask events down. From Figure 2(d), we see that the distribution between these two event types is very balanced, with market orders accounting for less than 5% of the level-I activity. As discussed above, only the level-I events are modelled and events occurring at deeper levels of the LOB discarded. One could of course increase $d_e$ to have a more granular classification of event types. Still, the focus of this paper is on the modelling of state dependence and the choice $\mathcal{E}=\left\{ask,bid\right\}$ already leads to interesting results whilst keeping the dimensionality low, which makes the results easier to visualise.

As the state variable we consider the bid–ask spread and the queue imbalance, giving rise to two models dubbed ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$, respectively. In ${\text{Model}}_{\text{S}}$, we set $\mathcal{X}=\{1,\mathbf{\text{{ 2+}}}\}$ ($d_x=2$), where the states correspond to the bid–ask spread being one tick ($X\left(t\right)=\mathbf{\text{{ 1}}}$) and two ticks or more ($X\left(t\right)=\mathbf{\text{{ 2+}}}$). Increasing the number of states beyond $d_x=2$ in this setting would not be practically relevant since the bid–ask spread is very rarely strictly wider than two ticks. The queue imbalance, used in ${\text{Model}}_{\text{QI}}$, is nowadays recognised as a popular trading signal with predictive power on the direction of the next price move (Cartea et al. 2018a). Denoting the total size of limit orders sitting at the ask price by $Q_{\text{ask}}\left(t\right)$, and defining $Q_{\text{bid}}\left(t\right)$ analogously, the queue imbalance can be expressed as $\text{QI}\left(t\right):=\frac{Q_{\text{bid}}\left(t\right)-Q_{\text{ask}}\left(t\right)}{Q_{\text{bid}}\left(t\right)+Q_{\text{ask}}\left(t\right)}\in [-1,1].$ (The denominator $Q_{\text{bid}}\left(t\right)+Q_{\text{ask}}\left(t\right)$ equals zero if and only if the LOB is empty at time t. In this case, it would be natural to define $\text{QI}\left(t\right):=0$, but this never occurs in our data set.) For example, the condition $\text{QI}\left(t\right)>0$ signals buy pressure and tends to be followed by an upwards price move. As in Cartea et al. (2018a), we split the interval $[-1,1]$ into $d_x=5$ bins of equal width which we label as follows: (7) $\mathcal{X}=\{\underset{[-1,-0.6)}{\underbrace{sell++}},\underset{[-0.6,-0.2)}{\underbrace{sell+}},\underset{[-0.2,0.2)}{\underbrace{\phantom{\text{+}}\mathbf{\text{{ neutral}}}}},\underset{[0.2,0.6)}{\underbrace{buy+}},\underset{[0.6,1]}{\underbrace{buy++}}\},$(7) whereby in ${\text{Model}}_{\text{QI}}$ the state variable $X\left(t\right)$ indicates bin where $\text{QI}\left(t\right)$ is located.

Finally, given the large number of observations we are dealing with (see Figure 2(b)), we use the exponential specification (6(6) $k_{e^{\prime }e}(t,x)={\alpha }_{e^{\prime }xe}\exp (-{\text{β}}_{e^{\prime }xe}t),t>0,e^{\prime },e\in \mathcal{E},x\in \mathcal{X},$(6) ) of the kernel $\mathit{k}$ in both models, as it leads to a significant reduction of computational cost in estimation, as discussed in Subsection 3.3. The full specifications of the models are summarised in Table 1.

Table 1. Summary of ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$.

	${\text{Model}}_{\text{S}}$	${\text{Model}}_{\text{QI}}$
Event types	ask and bid events, $d_e=2$, $\mathcal{E}=\left\{ask,bid\right\}$
State process	Bid–ask spread, $d_x=2$  $\mathcal{X}=\{1,\mathbf{\text{{ 2+}}}\}$	Queue imbalance, $d_x=5$  $\mathcal{X}=\{\mathbf{\text{{ sell++}}},\mathbf{\text{{ sell+}}},\mathbf{\text{{ neutral}}},\mathbf{\text{{ buy+}}},\mathbf{\text{{ buy++}}}\}$
Kernel $\mathit{k}$	Exponential, of the form (6(6) $k_{e^{\prime }e}(t,x)={\alpha }_{e^{\prime }xe}\exp (-{\text{β}}_{e^{\prime }xe}t),t>0,e^{\prime },e\in \mathcal{E},x\in \mathcal{X},$(6) )
No. of parameters	26	92

4.4. Visualising estimated excitation effects

We use the following approach to present our estimation results on self- and cross-excitation effects. For each triple $(e,e^{\prime },x)\in {\mathcal{E}}^2\times \mathcal{X}$, ML estimation produces an estimated excitation profile $t\mapsto k_{e^{\prime }e}(t,x)$, parameterised by the estimated impact coefficient ${\hat{\alpha }}_{e^{\prime }xe}$ and decay coefficient ${\widehat{\text{β}}}_{e^{\prime }xe}$. However, instead of reporting ${\hat{\alpha }}_{e^{\prime }xe}$ and ${\widehat{\text{β}}}_{e^{\prime }xe}$, we visualise the excitation profile by plotting the truncated $L^1$-norm $t\mapsto \Vert {\hat{k}}_{e^{\prime }e}(\cdot ,x){\Vert }_{1,t}:={\int }_0^t{\hat{k}}_{e^{\prime }e}(s,x)\mathrm{d}s,$ from which the magnitude and the effective timescale of the excitation effect is easier to gauge than from the numerical values of ${\hat{\alpha }}_{e^{\prime }xe}$ and ${\widehat{\text{β}}}_{e^{\prime }xe}$. In fact, taking guidance from the cluster representation of ordinary Hawkes processes (Hawkes and Oakes 1974), we can conveniently interpret $\Vert {\hat{k}}_{e^{\prime }e}(\cdot ,x){\Vert }_{1,t}$ as the average number of events of type e that have been directly triggered by an event of type $e^{\prime }$ in state x within t seconds of its occurrence. Further, we note that the full $L^1$-norm is given by $\Vert {\hat{k}}_{e^{\prime }e}(\cdot ,x){\Vert }_{1,\mathrm{\infty }}=\underset{t\to \mathrm{\infty }}{lim}\Vert {\hat{k}}_{e^{\prime }e}(\cdot ,x){\Vert }_{1,t}=\frac{{\hat{\alpha }}_{e^{\prime }xe}}{{\widehat{\text{β}}}_{e^{\prime }xe}}.$

4.5. Estimation results for ${\text{Model}}_{\text{S}}$

We estimate the model parameters $({\hat{\mathit{\varphi }}}^{\left(i\right)},{\hat{\mathit{\nu }}}^{\left(i\right)},{\hat{\mathit{\alpha }}}^{\left(i\right)},{\widehat{\text{β}}}^{\left(i\right)})$ for each trading day i in the sample by ML estimation, as explained in Section 3. Practical details on the numerical solution of the underlying optimisation problem can be found in the Appendix A.2.2.

The estimated transition distribution $\hat{\mathit{\varphi }}$ of ${\text{Model}}_{\text{S}}$, obtained by averaging over the daily estimates ${\hat{\mathit{\varphi }}}^{\left(i\right)}$, is presented in Figure 4(a). The state process describing the bid–ask spread exhibits persistent behaviour, in the sense that the probability of remaining in the current state is very high. We also observe higher likelihood of moving from state 2+ to 1 than vice versa, which is consistent with one-tick bid–ask spread being the equilibrium state for a large-tick stock like INTC. We also find that the transition probabilities are not sensitive to the event type, which is natural since the bid–ask spread not is expected to be influenced by the direction of orders, ceteris paribus.

Figure 4. Estimated transition distributions $\hat{\mathit{\varphi }}$ of ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$. We report the average of ${\hat{\mathit{\varphi }}}^{\left(i\right)}$ across the 250 trading days. (Daily estimates vary little from these averaged values.). (a) Transition probabilities of the bid–ask spread (${\text{Model}}_{\text{S}}$). (b) Transition probabilities of the queue imbalance (${\text{Model}}_{\text{QI}}$).

The estimation results on the excitation kernel, given in Figure 5(a), indicate that self-excitation effects surpass cross-excitation effects in both states. Their magnitude and timescales, however, vary between the two states. The effective timescales of these effects range from 0.1 to 100 milliseconds, which is in agreement with the predominantly algorithmic origin and multiscale nature of trading in modern electronic markets.

Figure 5. The estimated kernel $\hat{\mathit{k}}$ under ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$. Each panel describes self- or cross-excitation as indicated by its title, whilst each colour corresponds to a different state. For example, in Figure 5(a), the red curves in the second panel represent the estimates ${\hat{k}}_{e^{\prime }e}^{\left(i\right)}(\cdot ,x)$ where $e^{\prime }=\mathbf{\text{{ ask}}}$, $e=\mathbf{\text{{ bid}}}$ and $x=\mathbf{\text{{ 1}}}$. All daily estimates are superposed with one translucent curve for each day. An ‘aggregate’ kernel is represented by a solid line, computed using the median of ${\hat{\mathit{\alpha }}}^{\left(i\right)}$ and ${\widehat{\text{β}}}^{\left(i\right)}$ across the 250 trading days. (a) ${\text{Model}}_{\text{S}}$ (state variable: bid–ask spread). (b) ${\text{Model}}_{\text{QI}}$ (state variable: queue imbalance).

When the bid–ask spread increases to 2+, the magnitude of the self-excitation effects doubles whilst their timescale remains roughly the same. The timescale of cross-excitation, however, lengthens drastically, whilst their magnitude increases slightly. A plausible microstructural explanation for this pattern goes as follows. When the bid–ask spread is in state 2+, a trader can submit an aggressive limit order inside the spread, gaining queue priority at the cost of a less favourable price. A bid event can then be seen as a signal for an upwards price move, which may prompt limit orders from buyers seeking a favourable position in the new queue and cancellations of limit orders from sellers trying to avoid adverse selection. Traders who submit sell orders are, per contra, incentivised to wait and see if the expected price increase actually materialises.

The evolution of the base rate vector ${\hat{\mathit{\nu }}}^{\left(i\right)}$ throughout the 250 days of data is displayed in Figure 6. We find a remarkable balance between buyers and sellers (i.e. ${\nu }_{bid}\approx {\nu }_{ask}$). We also notice that the evolution of ${\hat{\mathit{\nu }}}^{\left(i\right)}$ mimics that of the total number of orders (Figure 2(b)), which suggests that the day-to-day variation in market activity is mainly due to exogenous factors (cf. the analysis of endogeneity in Subsection 4.8).

Figure 6. Estimated base rate vector ${\hat{\mathit{\nu }}}^{\left(i\right)}$ for ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$ over time (in number of events per second). (a) ${\text{Model}}_{\text{S}}$. (b) ${\text{Model}}_{\text{QI}}$.

4.6. Estimation results for ${\text{Model}}_{\text{QI}}$

The estimated transition probabilities of ${\text{Model}}_{\text{QI}}$, presented in Figure 4(b), convey a tendency to stick to the current state, similar to what is seen in ${\text{Model}}_{\text{S}}$. Here, however, this behaviour is more of an artefact—each ask and bid event, by definition, changes the queue imbalance but not necessarily the state variable that is confined to the bins (6(6) $k_{e^{\prime }e}(t,x)={\alpha }_{e^{\prime }xe}\exp (-{\text{β}}_{e^{\prime }xe}t),t>0,e^{\prime },e\in \mathcal{E},x\in \mathcal{X},$(6) ).

In contrast to ${\text{Model}}_{\text{S}}$, the estimated transition probabilities now depend on the event type and we observe remarkable mirror symmetry, whereby ${\hat{\varphi }}_{\mathsf{\text{{ bid}}}}$ equals, up to 1 percentage point, ${\hat{\varphi }}_{\mathsf{\text{{ ask}}}}$ with the order of states reversed. This symmetry is natural, given the definition $\text{QI}\left(t\right)$—a sell order always decreases the queue imbalance unless it is submitted inside the bid–ask spread or it depletes the current bid queue, whilst an analogous statement is true for buy orders. As in ${\text{Model}}_{\text{S}}$, we find again that the probability of a state transition is higher when it is towards the equilibrium state, here neutral, consistent with ideas about the resilience of the LOB (Large 2007).

Looking at the estimation results for the excitation kernel in Figure 5(b), we observe that, like in ${\text{Model}}_{\text{S}}$, self-excitation surpass cross-excitation, with the magnitude of the former and the timescale of the latter being manifestly sensitive to the current state. The mirror symmetry seen above in the context of the transition probabilities holds here as well, whereby it suffices to only speak about the results for ask events, whilst analogous conclusions can be drawn on bid events.

Even though the day-to-day variation of the estimates is more pronounced in this model compared to ${\text{Model}}_{\text{S}}$, some clear patterns emerge again. The self-excitation of ask events consistently increases under heavy sell pressure (state sell++). This can plausibly be explained by a combination of a flight to liquidity, through the submission of sell orders, and fear of adverse selection, leading to cancellations of buy orders, by traders expecting a downwards price move. Under mild buy pressure (state buy+), self-excitation decreases so that the corresponding kernel norm nearly halves.

In the state buy++ (heavy buy pressure), a closer look at the daily estimates of the self-excitation of ask events plotted in purple in Figure 5(b) reveals two distinct groups of curves above and below of the median curve. The daily estimates of the (full) kernel norms plotted in time in Figure 7 suggest that the self-excitation effect has undergone a structural break in early February 2018—the moment that suddenly marked the end of a year-long period of unusually low volatility in the US equity markets, dubbed ‘the return of volatility’ by some financial journalists. It is also worth mentioning that at the same time the price of INTC went above $50, at which point the large-tick character of the stock starts to weaken. The behaviour characterised by the lower group of curves (pre-February 2018) can be interpreted as traders expecting a price increase and thus deferring the submission of sell orders, whereas the upper group of curves (post-February 2018) hints at a tendency of sellers to seek an advantageous position in the ask queue. Indeed, a queue imbalance close to one implies that that only a very small amount of liquidity is available at the ask price. Thus, placing a sell limit order following a succession of sell limit orders from other traders allows one to acquire a good position in the queue, should it be replenished (which brings the queue imbalance back to equilibrium). If, however, the ask queue progresses towards depletion, one has still time to cancel the order whilst the sell limit orders at the front of the queue are matched with incoming buy market orders.

Figure 7. Estimated kernel norms $\Vert \hat{\mathit{k}}{\Vert }_{1,\mathrm{\infty }}$ for ${\text{Model}}_{\text{QI}}$ over time. The dashed line marks 5 February 2018 (‘the return of volatility’), a day when the CBOE Volatility Index (VIX) jumped by 116% to 38 points, a level not seen since August 2015. This day seems to have introduced a systematic change in the magnitude of cross- and self-excitation. The spike in cross-excitation that occurs on 21 February 2018 is linked to a sudden change in market behaviour around 14:00 on that day, when Intel Corporation rolled out patches for its most recent generation of processors.

Under buy pressure (states buy++ and buy+), the timescale of the cross-excitation from bid to ask events becomes almost as short as that of the self-excitation of bid events. This reflects the resilience of the LOB—in response to a bid event, ask events compete neck and neck with bid events precisely when they push the queue imbalance back towards the equilibrium state. Recall that the resilience of the LOB is also reinforced by the estimated state transition probabilities, as discussed above.

To exemplify the estimated intensity processes and their state dependence, a very brief extract from the estimated dynamics of ${\text{Model}}_{\text{QI}}$ is presented in Figure 8. In particular, we observe how pronounced the self-excitation of bid events becomes when the queue imbalance drops below $-0.6$, that is, to state sell++.

Figure 8. The upper panel depicts the evolution of the queue imbalance and level-I order flow of INTC on 13 February 2018 and the ask (red dots) and bid (blue dots) events. The lower panel displays the estimated intensities of ${\text{Model}}_{\text{QI}}$. (The self- and cross-excitation kernel norms of ${\text{Model}}_{\text{QI}}$ on 13 February 2018 are visualised in Figure A2.)

4.7. Goodness-of-fit diagnostics

To assess the goodness of fit of the estimated models, we examine the event residuals $r_n^e$ in sample (daily, 12:00–14:30) and out of sample (daily, 14:30–15:00) for both ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$. For comparison, we also estimated an ordinary Hawkes process ($d_x=1$) with exponential kernel for the same event types $\mathcal{E}=\left\{ask,bid\right\}$ and computed its event residuals. Since the state-dependent Hawkes process nests the ordinary Hawkes process, the former will by construction provide a better fit in sample than the latter. The Q–Q plots in Figure 9 show that this improvement in the goodness of fit extends out of sample, albeit the improvement is smaller than in sample. This is a confirmation that that ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$, and their state-dependent features in particular, are not overfitted.

Figure 9. In-sample (12:00–14:30) and out-of-sample (14:30–15:00) Q–Q plots of event residuals under ${\text{Model}}_{\text{S}}$, ${\text{Model}}_{\text{QI}}$ (state-dependent) and an ordinary Hawkes process (simple). The residuals of the ith day are computed using the ML estimates $({\hat{\mathit{\nu }}}^{\left(i\right)},{\hat{\mathit{\alpha }}}^{\left(i\right)},{\widehat{\text{β}}}^{\left(i\right)})$ obtained from the 12:00–14:30 period. The empirical quantiles are obtained by pooling the residuals of all 250 trading days. The two panels in each sub-figure correspond to the sequences of residuals $\left(r_n^e\right)$ for $e\in \mathcal{E}=\{\mathrm{a}\mathrm{s}\mathrm{k},\mathrm{b}\mathrm{i}\mathrm{d}\}$. (a) ${\text{Model}}_{\text{S}}$: in sample. (b) ${\text{Model}}_{\text{S}}$: out of sample. (c) ${\text{Model}}_{\text{QI}}$: in sample. (d) ${\text{Model}}_{\text{QI}}$: out of sample.

It should not come as a surprise that the improvement in goodness of fit provided by the state-dependent model looks meagre when one examines the Q–Q plots. Indeed, the behaviour of ${\text{Model}}_{\text{S}}$, ${\text{Model}}_{\text{QI}}$ and their ordinary Hawkes process alternative is quite similar when the bid–ask spread and queue imbalance are in their most likely states. It is only in the less likely states, 2+ in ${\text{Model}}_{\text{S}}$ and sell++ and buy++ in ${\text{Model}}_{\text{QI}}$, where the difference between the state-dependent and ordinary Hawkes process models becomes more pronounced. Thereby, the unconditional distribution of residuals does not vary much between these three models.

4.8. Endogeneity is state-dependent

The recent popularity of Hawkes processes in the modelling of high-frequency financial data partly stems from their ability to quantify the endogeneity of market activity. Indeed, based on the cluster representation of Hawkes processes (Hawkes and Oakes 1974) and the theory of branching processes (Harris 1963), the expected number of new events triggered by each event through self-excitation in a univariate ordinary Hawkes process equals the $L^1$-norm of its self-excitation kernel, provided it is less than one. The threshold one is the critical boundary for the stability of the process and validity of the cluster representation. Hawkes processes fitted to high-frequency financial data often exhibit kernel norms slightly below one, whilst the base rates tend to be relatively low in comparison. This phenomenon has been interpreted by Filimonov and Sornette (2012) and Hardiman et al. (2013) as evidence that most market events are endogenous, mere responses to earlier events, dwarfing the flow of less frequent exogenous events that are driven by new information. They dub the phenomenon critical reflexivity, which is a nod to George Soros and his reflexivity theory on the endogeneity of financial markets (Soros 1989).

In the context of state-dependent Hawkes processes, the kernel $\left(k_{e^{\prime }e}\right(\cdot ,x){)}_{e^{\prime },e\in \mathcal{E}}$, for any state $x\in \mathcal{X}$, defines a multivariate ordinary Hawkes process. Whilst in the multivariate case there is no direct analogue of the kernel norm—at least none with an equally clear-cut interpretation—the spectral radius $\rho \left(x\right)$ of the $d_e\times d_e$ matrix $\left(m_{ij}\right)$, where $m_{ij}:=\Vert k_{ji}(\cdot ,x){\Vert }_{1,\mathrm{\infty }}$, can be understood as a measure of endogeneity in this ordinary Hawkes process for each $x\in \mathcal{X}$. It also characterises the stability of the process, whereby $\rho \left(x\right)<1$ is a sufficient condition for the existence of and convergence to a stationary version (Brémaud and Massoulié 1996). Figure 10 displays the daily estimates of $\rho \left(x\right)$ for both ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$ as a function of $x\in \mathcal{X}$. We observe a remarkably clear pattern of $\rho \left(x\right)$ being almost uniformly higher in the disequilibrium states (2+, sell++, buy++) than in the equilibrium states (1, sell+, neutral, buy+). In particular, in ${\text{Model}}_{\text{S}}$, the spectral radius $\rho \left(\mathbf{\text{{ 2+}}}\right)$ is systematically above the critical value 1, whilst in ${\text{Model}}_{\text{QI}}$, the values of $\rho \left(sell++\right)$ and $\rho \left(buy++\right)$ are above 0.9 half of the time, exceeding 1 occasionally. These results thus open a new perspective on critical reflexivity, showing that it is in fact a largely state-dependent phenomenon, observed only in particular circumstances. They also lend credence to Soros's remark that ‘[e]ven in the financial markets demonstrably reflexive processes occur only intermittently’ (Soros 2008, p. 29).

Figure 10. The estimated spectral radius $\hat{\rho }\left(x\right)$ as a function of $x\in \mathcal{X}$ under ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$. The daily profiles $x\mapsto \hat{\rho }(x{)}^{\left(i\right)}$, $i=1,\dots ,250$, are represented by the red translucent curves. (a) ${\text{Model}}_{\text{S}}$. (b) ${\text{Model}}_{\text{QI}}$.

The increase in endogeneity in the disequilibrium states seems attributable to strategies employed by high-frequency traders (HFTs), who become active in these states, in anticipation of a price move. In a recent study, Lehalle and Neuman (2019) analyse a unique data set from Nasdaq Stockholm, where the identities of the buyer and seller in each transaction were disclosed by the exchange until 2014. In particular, they show that when the queue imbalance increases, the trading activity of market participants they classify as proprietary HFTs is amplified, in the direction of the imbalance. The pronounced sub-millisecond self-excitation effects seen in Figure 5(b), which are the key driver behind the high spectral radii $\rho \left(x\right)$, concur with the trading patterns observed by Lehalle and Neuman (2019). Besides its use as a trading signal, Lehalle and Neuman (2019) find the queue imbalance to be mean-reverting, which is similarly compatible with our results. A higher spectral radius in disequilibrium states corresponds here to an increase in market activity, which, reinforced by the structure of the estimated transition probabilities (Figure 4), is more likely to push the queue imbalance towards equilibrium than vice versa.

4.9. Event–state structure of limit order books

We could alternatively build a state-dependent variant of a Hawkes process in the following, conceptually simpler way. Using the representation $\stackrel{\mathbf{~}}{\mathit{N}}$ of the counting and state processes $(\mathit{N},X)$, one could specify an intensity of the form (8) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\nu }_{ex}+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }{x}^{\prime }ex}(t-s)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\left(s\right),t\ge 0,\\ & e\in \mathcal{E},x\in \mathcal{X},\end{array}$(8) instead of (3(3) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\varphi }_e\left(X\right(t),x)\\ & \times ({\nu }_e+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,x^{\prime }\left)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\right(s\left)\right),\\ & t\ge 0,e\in \mathcal{E},x\in \mathcal{X}.\end{array}$(3) ). This approach would in fact be tantamount to simply using a $d_e{d}_x$-dimensional ordinary Hawkes process.

The intensity (8(8) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\nu }_{ex}+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }{x}^{\prime }ex}(t-s)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\left(s\right),t\ge 0,\\ & e\in \mathcal{E},x\in \mathcal{X},\end{array}$(8) ) makes self- and cross-excitation state-dependent, but the simple structure of the state process X in Definition 2.1 is lost. Namely, under (8(8) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\nu }_{ex}+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }{x}^{\prime }ex}(t-s)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\left(s\right),t\ge 0,\\ & e\in \mathcal{E},x\in \mathcal{X},\end{array}$(8) ), the transition probabilities of the state process depend not only on the current state but on the entire history. However, LOBs enjoy a certain event–state structure—knowing the current state of the LOB and the characteristics of the next order suffices to (approximately) determine the next state. State-dependent Hawkes processes are by construction able to reproduce such an event-state structure and, therefore, compared to the alternative (8(8) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\nu }_{ex}+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }{x}^{\prime }ex}(t-s)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\left(s\right),t\ge 0,\\ & e\in \mathcal{E},x\in \mathcal{X},\end{array}$(8) ), we expect them to provide in general a better statistical description of the LOB. Moreover, the model given by (8(8) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\nu }_{ex}+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }{x}^{\prime }ex}(t-s)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\left(s\right),t\ge 0,\\ & e\in \mathcal{E},x\in \mathcal{X},\end{array}$(8) ) requires a kernel with $d_e^2{d}_x^2$ components whereas a state-dependent Hawkes process can be specified more parsimoniously, using only $d_e^2{d}_x$ components.

To compare this alternative model to state-dependent Hawkes processes, we estimate an ordinary Hawkes process with intensity (8(8) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\nu }_{ex}+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }{x}^{\prime }ex}(t-s)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\left(s\right),t\ge 0,\\ & e\in \mathcal{E},x\in \mathcal{X},\end{array}$(8) ), choosing an exponential form of the kernel $\mathit{k}$ and specifying the event types in $\mathcal{E}$ and the state space $\mathcal{X}$ as in ${\text{Model}}_{\text{QI}}$. Because of the higher dimensionality of this alternative model, estimation becomes computationally more expensive, which is why we use INTC data for May 2018 only. In Figure 11, this alternative model is compared to ${\text{Model}}_{\text{QI}}$ via Q–Q plots of the total residuals. We choose to display the results for a day when the alternative model provides one of its best fits, although the goodness of fit does not vary markedly across the 22 trading days in May 2018. Even though it has fewer parameters (92 as opposed to 210), we observe that ${\text{Model}}_{\text{QI}}$ provides a better in-sample fit. These empirical results thus underline the significance of the event–state structure of LOBs.

Figure 11. In-sample (12:00–14:30) Q–Q plots of total residuals of ${\text{Model}}_{\text{QI}}$ (state-dependent) and the alternative model given by (8(8) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\nu }_{ex}+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }{x}^{\prime }ex}(t-s)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\left(s\right),t\ge 0,\\ & e\in \mathcal{E},x\in \mathcal{X},\end{array}$(8) ) (complex) on 11 May 2018. Each panel corresponds to a sequence of residuals $\left({\tilde{r}}_n^{ex}\right)$ for all $e\in \mathcal{E}$ and $x\in \mathcal{X}$.

5. Discussion

State-dependent Hawkes processes enable us to model two-way interaction between a self- and cross-exciting point process, governing the temporal flow of events, and a state process, describing a system. In the context of LOB modelling, they provide a probabilistic foundation for a novel class of continuous-time models that encapsulate the feedback loop between the order flow and the shape of the LOB. Our estimation results, using these models and one year's worth of high-frequency LOB data on the stock of Intel Corporation, reveal that state dependence is indeed significant, as we uncover several robust patterns that persist throughout the daily estimation results. In particular, we find that market endogeneity, measured through the magnitude of self- and cross-excitation is state-dependent, being most pronounced in disequilibrium states of the LOB.

Our results also validate the event–state structure of LOBs that is embedded in the definition of a state-dependent Hawkes process. However, we do not claim that ${\text{Model}}_{\text{S}}$ and ${\text{Model}}_{\text{QI}}$ would be the best possible representations of the aforementioned feedback loop—we recognise that they could be refined as follows:

1. The exponential excitation kernels could be replaced by power laws, motivated by the non-parametric estimation results on ordinary Hawkes processes (Hardiman et al. 2013, Bacry et al. 2016). A numerically more convenient alternative would be to use a linear combination of exponentials within the parametric framework to mimic slow power-law decay (Rambaldi et al. 2015, Lu and Abergel 2018).

2. The base rates could be made state-dependent by replacing ${\nu }_e$ with ${\nu }_e\left(X\right(t\left)\right)$ in (1(1) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(s\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E};\end{array}$(1) ). Note that in this case, the model would contain both a continuous-time Markov chain ($\mathit{k}\equiv 0$) and an ordinary Hawkes process as special cases.

3. One might argue for excitation kernels that allow for negative values, to capture inhibition effects that are known to exist in LOB data (Lu and Abergel 2018). To ensure the non-negativity of the intensity, this would require transforming the right-hand side of (1(1) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(s\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E};\end{array}$(1) ) by a non-linear function.

4. More granular event types and states would provide more nuanced understanding of the LOB dynamics. Besides, notice that the present framework accommodates multiple state variables. For example, X could in fact jointly represent both the bid–ask spread and the queue imbalance using the state space $\begin{array}{rl}\mathcal{X}& =\{(1,sell++),\dots ,(1,buy++),\\ & (2+,sell++),\dots ,\text{({{ 2+}, { buy++})}\}.\end{array}$

However, no matter how one modifies the intensity (1(1) $\begin{array}{rl}{\lambda }_e\left(t\right)& ={\nu }_e+\sum _{e^{\prime }\in \mathcal{E}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,X(s\left)\right)\mathrm{d}{N}_{e^{\prime }}\left(s\right),\\ & t\ge 0,e\in \mathcal{E};\end{array}$(1) ) by implementing any of (i)–(iv) to incorporate one's views on the feedback loop, the event–state structure of the model remains intact. The process still falls within the class of hybrid marked point processes (Morariu-Patrichi and Pakkanen 2018), implying that the theoretical results of this paper (existence, uniqueness, separability of the likelihood function) still apply.

Acknowledgments

Maxime Morariu-Patrichi gratefully acknowledges the Mini-DTC scholarship awarded by the Mathematics Department of Imperial College London. Mikko S. Pakkanen acknowledges partial support from CREATES (DNRF78), funded by the Danish National Research Foundation, from EPSRC through the Platform Grant EP/I019111/1 and from the CFM–Imperial Institute of Quantitative Finance.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Appendix

A.1. Proofs

Proof of Theorem 2.3.

The statement follows by applying Theorem 2.13 in Morariu-Patrichi and Pakkanen (2018). We only need to check that if $(\mathit{N},X)$ is a state-dependent Hawkes process, then $\stackrel{\mathbf{~}}{\mathit{N}}$ admits an ${\mathbb{F}}^{\stackrel{\mathbf{~}}{\mathit{N}}}$-intensity. By Proposition 3.1 in Jacod (1975), $\stackrel{\mathbf{~}}{\mathit{N}}$ admits an ${\mathbb{F}}^{\stackrel{\mathbf{~}}{\mathit{N}}}$-compensator Λ given by $\mathrm{\Lambda }(\mathrm{d}t,\mathrm{d}e,\mathrm{d}x)=\sum _{n\in \mathbb{N}}\frac{G_n(\mathrm{d}t-T_n,\mathrm{d}e,\mathrm{d}x)}{G_n\left(\right[t-T_n,\mathrm{\infty }],\mathcal{E},\mathcal{X})}{\text{1}}_{\{T_n<t\le T_{n+1}\}},$ and thus this holds if the conditional distributions $G_n$ of $(T_{n+1}-T_n,E_{n+1},X_{n+1})$ with respect to ${\mathcal{F}}_{T_n}^{\stackrel{\mathbf{~}}{\mathit{N}}}$ are absolutely continuous with respect to $\mathrm{d}t{\mu }_e\left(\mathrm{d}e\right){\mu }_x\left(\mathrm{d}x\right)$ on $(0,T_{n+1}-T_n]\times \mathcal{E}\times \mathcal{X}$, where ${\mu }_e$ and ${\mu }_x$ are the counting measures on $\mathcal{E}$ and $\mathcal{X}$, respectively. Moreover, by the definition of $\stackrel{\mathbf{~}}{\mathit{N}}$, $\mathrm{\Lambda }(dt,de,\mathcal{X})$ is an ${\mathbb{F}}^{\stackrel{\mathbf{~}}{\mathit{N}}}$-compensator of $\mathit{N}$. But since $\mathit{N}$ admits an ${\mathbb{F}}^{\stackrel{\mathbf{~}}{\mathit{N}}}$-intensity, by uniqueness of the compensator (Kallenberg 2017, Theorem 1.25, p. 39), we necessarily have with probability one that $\begin{array}{rl}& \sum _{x\in \mathcal{X}}{G}_n(dt-T_n,\{e\},\{x\left\}\right){\text{1}}_{\{T_n<t\le T_{n+1}\}}=f(t,e)\mathrm{d}t{\text{1}}_{\{T_n<t\le T_{n+1}\}},\\ & e\in \mathcal{E},n\in \mathbb{N},\end{array}$ for some ${\mathcal{F}}_{T_n}^{\stackrel{\mathbf{~}}{\mathit{N}}}\otimes \mathcal{B}\left(\mathbb{R}\right)\otimes \mathcal{P}\left(\mathcal{E}\right)$-measurable function f, which concludes the proof.

Proof of Theorem 2.4.

By Theorem 2.3, it is sufficient to show that each of the two conditions ensures the existence and uniqueness of a non-explosive point process $\stackrel{\mathbf{~}}{\mathit{N}}$ with an ${\mathbb{F}}^{\stackrel{\mathbf{~}}{\mathit{N}}}$-intensity given by (3(3) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\varphi }_e\left(X\right(t),x)\\ & \times ({\nu }_e+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,x^{\prime }\left)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\right(s\left)\right),\\ & t\ge 0,e\in \mathcal{E},x\in \mathcal{X}.\end{array}$(3) ), which is achieved by applying Theorems 2.17 and 2.21 in Morariu-Patrichi and Pakkanen (2018).

Proof of Theorem 3.1.

By Theorem 2.3, we know that $\stackrel{\mathbf{~}}{\mathit{N}}$ has intensity $\stackrel{\mathbf{~}}{\mathit{\lambda }}$, which is given by (3(3) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\varphi }_e\left(X\right(t),x)\\ & \times ({\nu }_e+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,x^{\prime }\left)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\right(s\left)\right),\\ & t\ge 0,e\in \mathcal{E},x\in \mathcal{X}.\end{array}$(3) ). Hence, by applying Proposition 7.3.III in Daley and Vere-Jones (2003, p. 251), we can express the log likelihood function as (A1) $\ln \mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })=\sum _{n=1}^N\ln {\tilde{\lambda }}_{e_n{x}_n}\left(t_n\right)-{\int }_0^T\sum _{e\in \mathcal{E},x\in \mathcal{X}}{\tilde{\lambda }}_{ex}\left(t\right)\mathrm{d}t.$(A1) Plugging (3(3) $\begin{array}{rl}{\tilde{\lambda }}_{ex}\left(t\right)& ={\varphi }_e\left(X\right(t),x)\\ & \times ({\nu }_e+\sum _{e^{\prime }\in \mathcal{E},x^{\prime }\in \mathcal{X}}{\int }_{[0,t)}{k}_{e^{\prime }e}(t-s,x^{\prime }\left)\mathrm{d}{\tilde{N}}_{e^{\prime }{x}^{\prime }}\right(s\left)\right),\\ & t\ge 0,e\in \mathcal{E},x\in \mathcal{X}.\end{array}$(3) ) in (A1(A2) ${ϵ}_{rel}:=\frac{{\hat{\theta }}_{i_{j^{\star }}}-{\theta }_{i_{j^{\star }}}}{{\theta }_{i_{j^{\star }}}},\text{where}{j}^{\star }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_j\frac{|{\hat{\theta }}_{i_j}-{\theta }_{i_j}|}{{\theta }_{i_j}}.$(A2) ) and using that $\sum _{x^{\prime }\in \mathcal{X}}{\varphi }_e(x,x^{\prime })=1$, $e\in \mathcal{E}$, $x\in \mathcal{X}$, yields (5(5) $\begin{array}{rl}\ln \mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })& =\sum _{n=1}^N\ln {\varphi }_{e_n}(x_{n-1},x_n)+\sum _{n=1}^n\ln {\lambda }_{e_n}\left(t_n\right)\\ & -{\int }_0^T\sum _{e\in \mathcal{E}}{\lambda }_e\left(t\right)\mathrm{d}t.\end{array}$(5) ), from which it is immediate that $(\hat{\mathit{\varphi }},\hat{\mathit{\nu }},\hat{\mathit{\theta }})\in \text{arg}{\text{max}}_{\mathit{\varphi },\mathit{\nu },\mathit{\theta }}\mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })$ if and only if $\{\begin{array}{ll}\hat{\mathit{\varphi }}& \in {\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_{\mathit{\varphi }}{\displaystyle \sum _{n=1}^N\ln {\varphi }_{e_n}(x_{n-1},x_n),}\\ (\hat{\mathit{\nu }},\hat{\mathit{\theta }})& \in {\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_{\mathit{\nu },\mathit{\theta }}{\displaystyle \sum _{n=1}^N\ln {\lambda }_{e_n}\left(t_n\right)-{\int }_0^T\sum _{e\in \mathcal{E}}{\lambda }_e\left(t\right)\mathrm{d}t,}\end{array}$ where the first optimisation problem is performed under the constraint that ${\mathit{\varphi }}_e$ is a transition probability matrix, $e\in \mathcal{E}$. By solving this optimisation problem with the method of Lagrange multipliers, we obtain the claimed expression for ${\hat{\varphi }}_e(x,x^{\prime })$.

Proof of Theorem 3.3.

The statement follows directly from Theorem 1 in Brown and Nair (1988). Note that this theorem requires that, with probability one, ${\int }_0^t{\lambda }_e\left(s\right)\mathrm{d}s\to \mathrm{\infty }$ and ${\int }_0^t{\tilde{\lambda }}_{ex}\left(s\right)\mathrm{d}s\to \mathrm{\infty }$ as $t\to \mathrm{\infty }$, $e\in \mathcal{E}$, $x\in \mathcal{X}$. The first condition is satisfied because, in Definition 2.1, we assume that all the base rates are strictly positive ($\mathit{\nu }\in {\mathbb{R}}_{>0}^{d_e}$). By Lemma 17 in Brémaud (1981, p. 41), the second condition is equivalent to ${\tilde{N}}_{ex}\left(t\right)\to \mathrm{\infty }$, $t\to \mathrm{\infty }$, with probability one, $e\in \mathcal{E}$, $x\in \mathcal{X}$, which is assumed here.

A.2. Details on maximum likelihood estimation

A.2.1. Formulae for state-dependent Hawkes processes with exponential kernels

In the case of exponential kernel $\mathit{k}$ given by (6(6) $k_{e^{\prime }e}(t,x)={\alpha }_{e^{\prime }xe}\exp (-{\text{β}}_{e^{\prime }xe}t),t>0,e^{\prime },e\in \mathcal{E},x\in \mathcal{X},$(6) ), the following formulae can be derived for the second and third terms of the log likelihood function in (5(5) $\begin{array}{rl}\ln \mathcal{L}(\mathit{\varphi },\mathit{\nu },\mathit{\theta })& =\sum _{n=1}^N\ln {\varphi }_{e_n}(x_{n-1},x_n)+\sum _{n=1}^n\ln {\lambda }_{e_n}\left(t_n\right)\\ & -{\int }_0^T\sum _{e\in \mathcal{E}}{\lambda }_e\left(t\right)\mathrm{d}t.\end{array}$(5) ), denoted by $l_{+}$ and $l_{-}$, respectively. Here, we consider a general time horizon $(t_0,T]$, meaning that the origin of time is not necessarily $t_0=0$ and that the times $t_n^{ex}\le t_0$ are treated like an initial condition, and we have $\begin{array}{rl}{l}_{+}& =\sum _e\sum _{t_n^e>t_0}\ln {\lambda }_e\left(t_n^e\right)\\ & =\sum _e\sum _{t_n^e>t_0}\ln ({\nu }_e+\sum _{e^{\prime }{x}^{\prime }}{\alpha }_{e^{\prime }{x}^{\prime }e}{S}_{e^{\prime }{x}^{\prime }e}(t_n^e\left)\right),\\ S_{e^{\prime }{x}^{\prime }e}(s,t)& :=\sum _{s<t_i^{e^{\prime }{x}^{\prime }}<t}\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t-t_i^{e^{\prime }{x}^{\prime }}\left)\right),\\ S_{e^{\prime }{x}^{\prime }e}\left(t\right)& :=S_{e^{\prime }{x}^{\prime }e}(-\mathrm{\infty },t),\\ S_{e^{\prime }{x}^{\prime }e}\left(t\right)& =\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t-s\left)\right)S_{e^{\prime }{x}^{\prime }e}(-\mathrm{\infty },s]+S_{e^{\prime }{x}^{\prime }e}(s,t),\\ l_{-}& ={\int }_{t_0}^T\sum _e{\lambda }_e\left(t\right)\mathrm{d}t\\ & ={\int }_{t_0}^T\sum _e({\nu }_e+\sum _{e^{\prime }{x}^{\prime }}\sum _{t_i^{e^{\prime }{x}^{\prime }}<t}{\alpha }_{e^{\prime }{x}^{\prime }e}\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t-t_i^{e^{\prime }{x}^{\prime }}))\mathrm{d}t\\ & =\sum _e{\nu }_e(T-t_0)\\ & +\sum _{e^{\prime }{x}^{\prime }}\sum _{t_i^{e^{\prime }{x}^{\prime }}\le t_0}\frac{{\alpha }_{e^{\prime }{x}^{\prime }e}}{{\text{β}}_{e^{\prime }{x}^{\prime }e}}(\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t_0-t_i^{e^{\prime }{x}^{\prime }}))\\ & -\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(T-t_i^{e^{\prime }{x}^{\prime }}\left)\right)\\ & +\sum _{t_i^{e^{\prime }{x}^{\prime }}>t_0}\frac{{\alpha }_{e^{\prime }{x}^{\prime }e}}{{\text{β}}_{e^{\prime }{x}^{\prime }e}}(1-\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(T-t_i^{e^{\prime }{x}^{\prime }})).\end{array}$ The gradients can then also be computed via $\begin{array}{rl}\frac{\mathrm{\partial }{l}_{+}}{\mathrm{\partial }{\nu }_k}& =\sum _{t_n^k>t_0}{({\nu }_k+\sum _{e^{\prime }{x}^{\prime }}{\alpha }_{e^{\prime }{x}^{\prime }k}{S}_{e^{\prime }{x}^{\prime }k}(t_n^k\left)\right)}^{-1},\\ \frac{\mathrm{\partial }{l}_{+}}{\mathrm{\partial }{\alpha }_{ijk}}& =\sum _{t_n^k>t_0}\frac{S_{ijk}\left(t_n^e\right)}{{\nu }_k+\sum _{e^{\prime }{x}^{\prime }}{\alpha }_{e^{\prime }{x}^{\prime }k}{S}_{e^{\prime }{x}^{\prime }k}\left(t_n^k\right)},\\ \frac{\mathrm{\partial }{l}_{+}}{\mathrm{\partial }{\text{β}}_{ijk}}& =-\sum _{t_n^k>t_0}\frac{{\alpha }_{ijk}{S}_{ijk}^{\left(1\right)}\left(t_n^e\right)}{{\nu }_k+\sum _{e^{\prime }{x}^{\prime }}{\alpha }_{e^{\prime }{x}^{\prime }k}{S}_{e^{\prime }{x}^{\prime }k}\left(t_n^k\right)},\\ S_{e^{\prime }{x}^{\prime }e}^{\left(1\right)}(s,t)& :=\sum _{s<t_i^{e^{\prime }{x}^{\prime }}<t}(t-t_i^{e^{\prime }{x}^{\prime }})\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t-t_i^{e^{\prime }{x}^{\prime }}\left)\right),\\ S_{e^{\prime }{x}^{\prime }e}^{\left(1\right)}\left(t\right)& =S_{e^{\prime }{x}^{\prime }e}^{\left(1\right)}(-\mathrm{\infty },t),\\ S_{e^{\prime }{x}^{\prime }e}^{\left(1\right)}\left(t\right)& =\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t-s\left)\right)\\ & \times \left(S_{e^{\prime }{x}^{\prime }e}^{\left(1\right)}\right(-\mathrm{\infty },s]+(t-s\left)S_{e^{\prime }{x}^{\prime }e}\right(-\mathrm{\infty },s\left]\right)+S_{e^{\prime }{x}^{\prime }e}^{\left(1\right)}(s,t),\\ \frac{\mathrm{\partial }{l}_{-}}{\mathrm{\partial }{\nu }_k}& =T-t_0,\\ \frac{\mathrm{\partial }{l}_{-}}{\mathrm{\partial }{\alpha }_{e^{\prime }{x}^{\prime }e}}& =\sum _{t_i^{e^{\prime }{x}^{\prime }}\le t_0}\frac{1}{{\text{β}}_{e^{\prime }{x}^{\prime }e}}(\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t_0-t_i^{e^{\prime }{x}^{\prime }}))\\ & -\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(T-t_i^{e^{\prime }{x}^{\prime }}\left)\right)\\ & +\sum _{t_i^{e^{\prime }{x}^{\prime }}>t_0}\frac{1}{{\text{β}}_{e^{\prime }{x}^{\prime }e}}(1-\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(T-t_i^{e^{\prime }{x}^{\prime }})),\\ \frac{\mathrm{\partial }{l}_{-}}{\mathrm{\partial }{\text{β}}_{e^{\prime }{x}^{\prime }e}}& =\sum _{t_i^{e^{\prime }{x}^{\prime }}\le t_0}\frac{{\alpha }_{e^{\prime }{x}^{\prime }e}}{{\text{β}}_{e^{\prime }{x}^{\prime }e}}\left(\right(T-t_i^{e^{\prime }{x}^{\prime }})\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(T-t_i^{e^{\prime }{x}^{\prime }}))\\ & -(t_0-t_i^{e^{\prime }{x}^{\prime }})\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t_0-t_i^{e^{\prime }{x}^{\prime }}\left)\right))\\ & -\frac{{\alpha }_{e^{\prime }{x}^{\prime }e}}{{\text{β}}_{e^{\prime }{x}^{\prime }e}^2}(\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t_0-t_i^{e^{\prime }{x}^{\prime }}))\\ & -\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(T-t_i^{e^{\prime }{x}^{\prime }}\left)\right)\\ & +\sum _{t_i^{e^{\prime }{x}^{\prime }}>t_0}\frac{{\alpha }_{e^{\prime }{x}^{\prime }e}}{{\text{β}}_{e^{\prime }{x}^{\prime }e}}(T-t_i^{e^{\prime }{x}^{\prime }})\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(T-t_i^{e^{\prime }{x}^{\prime }}\left)\right)\\ & -\frac{{\alpha }_{e^{\prime }{x}^{\prime }e}}{{\text{β}}_{e^{\prime }{x}^{\prime }e}^2}(1-\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(T-t_i^{e^{\prime }{x}^{\prime }})\left)\right).\end{array}$ Besides, the efficient computation of the residuals is based on the identity $\begin{array}{rl}{\int }_s^t{\lambda }_e\left(u\right)\mathrm{d}u& =(t-s){\nu }_e+\sum _{e^{\prime }{x}^{\prime }}{C}_{e^{\prime }{x}^{\prime }e}\left(t\right)-C_{e^{\prime }{x}^{\prime }e}\left(s\right)\\ & -\sum _{e^{\prime }{x}^{\prime }}\frac{{\alpha }_{e^{\prime }{x}^{\prime }e}}{{\text{β}}_{e^{\prime }{x}^{\prime }e}}\left(S_{e^{\prime }{x}^{\prime }e}\right(t)-S_{e^{\prime }{x}^{\prime }e}(s\left)\right),\\ C_{e^{\prime }{x}^{\prime }e}\left(t\right)& :=\sum _{t_i^{e^{\prime }{x}^{\prime }}<t}\frac{{\alpha }_{e^{\prime }{x}^{\prime }e}}{{\text{β}}_{e^{\prime }{x}^{\prime }e}}\exp (-{\text{β}}_{e^{\prime }{x}^{\prime }e}(t_0\vee t_i^{e^{\prime }{x}^{\prime }}-t_i^{e^{\prime }{x}^{\prime }}\left)\right).\end{array}$ The accompanying Python library mpoints implements the above formulae in C via the Cython extension of Python, allowing us to drastically reduce the computation time (up to 300 times faster computations compared to a plain Python implementation using NumPy). This played a crucial role in making the present study computationally feasible.

Table A1. Parameter values for Specification 1.

(a) Impact coefficients.
${\alpha }_{i1j}$		${\alpha }_{i2j}$		${\alpha }_{i3j}$		${\alpha }_{i4j}$		${\alpha }_{i5j}$
11,314	2,311	8,098	3,556	17,016	270	3,198	535	34,747	220
139	33,460	724	4,497	178	15,497	3,778	7,943	2,353	12,452
(b) Decay coefficients.
${\text{β}}_{i1j}$		${\text{β}}_{i2j}$		${\text{β}}_{i3j}$		${\text{β}}_{i4j}$		${\text{β}}_{i5j}$
17,480	9,563	15,187	19,637	29,890	1,456	6,493	3,071	45,546	1,097
613	44,994	4,614	8,384	1,026	27,522	20,537	14,272	10,443	18,988
(c) Transition distribution.
${\varphi }_1(i,j)$					${\varphi }_2(i,j)$
.89	.02	.05	.03	.01	.78	.19	.02	<.01	<.01
.13	.83	.01	<.01	.02	<.01	.75	.2	.02	<.01
.02	.17	.79	<.01	.02	.02	<.01	.79	.17	.02
<.01	.02	.21	.75	.01	.02	<.01	<.01	.82	.14
<.01	<.01	.02	.19	.78	<.01	.03	.06	.02	.89
(d) Base rate vector.
${\nu }_i$
2.2
2.2

Table A2. Parameter values for Specification 2.

(a) Impact coefficients.
${\alpha }_{i1j}$		${\alpha }_{i2j}$		${\alpha }_{i3j}$		${\alpha }_{i4j}$		${\alpha }_{i5j}$
2	10	1	3	1,000	30	2,000	40	100	1,000
10	2	3	1	20	3,000	2,000	50	60	2,000
(b) Decay coefficients.
${\text{β}}_{i1j}$		${\text{β}}_{i2j}$		${\text{β}}_{i3j}$		${\text{β}}_{i4j}$		${\text{β}}_{i5j}$
10	15	8	4	3,000	500	6,000	160	500	8,000
15	10	4	8	1000	5,000	10,000	300	120	5,000
(c) Transition distribution.
${\varphi }_1(i,j)$					${\varphi }_2(i,j)$
.7	.3	.0	.0	.0	.0	.1	.2	.3	.4
.1	.8	.1	.0	.0	.2	.1	.4	.2	.1
.0	.1	.6	.2	.1	.1	.3	.1	.3	.2
.2	.2	.3	.1	.2	.0	.0	.1	.8	.1
.1	.3	.3	.1	.2	.1	.0	.1	.1	.7
(d) Base rate vector.
${\nu }_i$
5
1

A.2.2. Numerical optimisation

Similarly to ordinary Hawkes processes (Lu and Abergel 2018), the maximum likelihood (ML) search (6(6) $k_{e^{\prime }e}(t,x)={\alpha }_{e^{\prime }xe}\exp (-{\text{β}}_{e^{\prime }xe}t),t>0,e^{\prime },e\in \mathcal{E},x\in \mathcal{X},$(6) ) is broken down into $d_e$ separate optimisation problems. For every $e\in \mathcal{E}$, we have an independent optimisation problem that involves only the parameters ${\nu }_e$, $({\alpha }_{e^{\prime }xe}{)}_{e^{\prime }\in \mathcal{E},x\in \mathcal{X}}$, $({\text{β}}_{e^{\prime }xe}{)}_{e^{\prime }\in \mathcal{E},x\in \mathcal{X}}$, and which we solve using a conjugate gradient method.

More precisely, we call the minimize function in the scipy.optimize Python package and use the method TNC. Three random sets of parameters are generated and used as alternative initial guesses. An ordinary Hawkes process ($d_x=1$) is estimated before the state-dependent one and its parameters are also used as an initial guess. Moreover, the estimate of the previous day is used as another initial guess and, thus, a total of five different initial guesses are employed. For ${\text{Model}}_{\text{QI}}$, around 50 iterations and 400 function evaluations are required for each day, event type $e\in \mathcal{E}$ and initial guess.

As this optimisation problem is non-convex, the above procedure might converge to a mere local maximum instead of the global one (Lu and Abergel 2018). Nevertheless, in a Monte Carlo experiment, reported in the following subsection, this conjugate gradient method returns estimates that are consistently concentrated near the true parameter values (see also Figure A2).

A.2.3. Finite-sample performance

We assess the finite-sample performance of the ML estimator in a small Monte Carlo experiment using a state-dependent Hawkes process with exponential kernel of the form (6(6) $k_{e^{\prime }e}(t,x)={\alpha }_{e^{\prime }xe}\exp (-{\text{β}}_{e^{\prime }xe}t),t>0,e^{\prime },e\in \mathcal{E},x\in \mathcal{X},$(6) ) when $d_e=2$ and $d_x=5$. The parameters are naturally split into four groups: the transition probabilities in $\mathit{\varphi }$, the base rates in $\mathit{\nu }$, the impact coefficients in $\mathit{\alpha }$ and the decay coefficients in $\text{β}$. For each group of parameters $({\theta }_{i_j}{)}_{j=1,\dots ,p_i}$ and their estimator $({\hat{\theta }}_{i_j}{)}_{j=1,\dots ,p_i}$, we define the worst relative error as (A2) ${ϵ}_{rel}:=\frac{{\hat{\theta }}_{i_{j^{\star }}}-{\theta }_{i_{j^{\star }}}}{{\theta }_{i_{j^{\star }}}},\text{where}{j}^{\star }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_j\frac{|{\hat{\theta }}_{i_j}-{\theta }_{i_j}|}{{\theta }_{i_j}}.$(A2)

For two different sets of parameter values and four different sample sizes, we estimate the distribution of ${ϵ}_{rel}$ for each of the four groups of parameters using Monte Carlo via Algorithm 2.5. However, for the transition probabilities, we measure instead the worst absolute error, defined by replacing the denominator in (A2(A2) ${ϵ}_{rel}:=\frac{{\hat{\theta }}_{i_{j^{\star }}}-{\theta }_{i_{j^{\star }}}}{{\theta }_{i_{j^{\star }}}},\text{where}{j}^{\star }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_j\frac{|{\hat{\theta }}_{i_j}-{\theta }_{i_j}|}{{\theta }_{i_j}}.$(A2) ) by one. The first set of parameter values (Specification 1, Table A1) is constructed simply by averaging the daily estimates of ${\text{Model}}_{\text{QI}}$ using INTC data over the 19 trading days of February 2018. The second set of parameter values (Specification 2, Table A2) is artificial, chosen to produce more drastic changes in behaviour from one state to another. The results are displayed in Figure A1 and indeed support the conjectured consistency of the ML estimators of the state-dependent Hawkes process. Note that the observed bimodality is due to the fact that the worst relative error inherently alternates between positive and negative values.

Figure A1. Violin plots of the worst estimation errors (A2(A2) ${ϵ}_{rel}:=\frac{{\hat{\theta }}_{i_{j^{\star }}}-{\theta }_{i_{j^{\star }}}}{{\theta }_{i_{j^{\star }}}},\text{where}{j}^{\star }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_j\frac{|{\hat{\theta }}_{i_j}-{\theta }_{i_j}|}{{\theta }_{i_j}}.$(A2) ) under two different sets of parameter values (specifications 1 and 2). For every specification and sample size N, we simulate 100 paths with sample size N and perform ML estimation for each of them. The true parameters are used as the initial guess in the optimisation procedure to speed up estimation and reduce the computational cost. (a) Specification 1. (b) Specification 2.

Figure A2. Uncertainty quantification of the estimated excitation profiles for INTC on 13 February 2018 under ${\text{Model}}_{\text{QI}}$. The parametric bootstrap procedure involves simulating 100 paths covering 2.5 hours of trading using the ML estimates of the parameters on the considered day and applying ML estimation again to each of the simulated paths. Three random sets of parameters are used as the initial guess in the optimisation procedure. We use the 100 estimates to compute a 99%-confidence interval for the truncated kernel norm (translucent area). The solid line corresponds to the ML estimates using the original INTC data.

A.2.4. Uncertainty quantification via parametric bootstrap

The uncertainty of ML estimates of the parameters of a state-dependent Hawkes process can be quantified using a parametric bootstrap procedure. The procedure entails simulating realisations of the state-dependent Hawkes process under the estimated parameter values using Algorithm 2.5 and then applying ML estimation again to each simulated realisation, producing a sample of estimates that approximates the distribution of the ML estimator. To exemplify the method, we apply it here to ${\text{Model}}_{\text{QI}}$ estimates for INTC on 13 February 2018. The results are presented in Figure A2, and we observe that the uncertainty of the estimates is negligible compared to the estimated excitation patterns.