Online learning of order flow and market impact with Bayesian change-point detection methods

Abstract

Financial order flow exhibits a remarkable level of persistence, wherein buy (sell) trades are often followed by subsequent buy (sell) trades over extended periods. This persistence can be attributed to the division and gradual execution of large orders. Consequently, distinct order flow regimes might emerge, which can be identified through suitable time series models applied to market data. In this paper, we propose the use of Bayesian online change-point detection (BOCPD) methods to identify regime shifts in real-time and enable online predictions of order flow and market impact. To enhance the effectiveness of our approach, we have developed a novel BOCPD method using a score-driven approach. This method accommodates temporal correlations and time-varying parameters within each regime. Through empirical application to NASDAQ data, we have found that: (i) Our newly proposed model demonstrates superior out-of-sample predictive performance compared to existing models that assume i.i.d. behavior within each regime; (ii) When examining the residuals, our model demonstrates good specification in terms of both distributional assumptions and temporal correlations; (iii) Within a given regime, the price dynamics exhibit a concave relationship with respect to time and volume, mirroring the characteristics of actual large orders; (iv) By incorporating regime information, our model produces more accurate online predictions of order flow and market impact compared to models that do not consider regimes.

Keywords:

• Regime switching models
• Bayesian methods
• Price impact
• Metaorders
• Online learning algorithm

1. Introduction

The study and modeling of order flow and market impact in financial markets hold paramount importance for understanding the incorporation of private information into prices and designing effective trading algorithms that consider transaction costs. A substantial body of literature (see for example Bouchaud et al. (2009) and Lillo (2023) and references therein) has revealed that the joint modeling of impact and order flow is more intricate than initially presumed.

The persistence and autocorrelation of signed trade order flow¹ have been extensively documented since the works of Lillo and Farmer (2004) and Bouchaud et al. (2004). This persistence aligns with a long memory process, suggesting that a realistic market impact model should combine statistically efficient prices with correlated order flow. The introduction of transient impact models, also known as propagator models (Bouchaud et al. 2004), successfully accomplishes this goal. Additionally, empirical evidence has attributed the temporal persistence of order flow primarily to order splitting, as discussed in Tóth et al. (2015). Order splitting refers to the common practice of large investors incrementally executing their orders, termed ‘metaorders,’ through several smaller trades known as ‘child orders.’ The model proposed by Lillo et al. (2005) quantitatively establishes a relationship between the autocorrelation of order flow and the distribution of metaorder sizes. In essence, this model postulates that the strong serial dependence arises from the optimal execution strategies employed by institutional investors, which leads to persistent order flow.

From an econometric perspective, this notion might be connected to the well-known fact (see Granger and Hyung 1999, Mikosch and Starica 1999, Diebold and Inoue 2001) that long memory time series can be (approximately) generated by regime-shift models, where each regime exhibits short memory and heterogeneous lengths. Regime shift models gained popularity around two decades ago when they were proposed to explain the long-range memory of volatility, as seen in Diebold and Inoue (2001).

This paper proposes to use regime shift models to describe order flow time series, with the objectives of: (i) econometrically explaining the long memory of order flow; (ii) enhancing the prediction of order flow and price dynamics through detected regimes; and (iii) suggesting a connection between regimes and the execution of metaorders. Unlike many conventional regime shift approaches that require a predetermined number of regimes (e.g. Hidden Markov Models), we focus on a model that allows online detection of change-points (CPs) to identify the occurrence of new regimes. Existing algorithms typically operate offline, recursively segmenting the time series ex-post into increasingly smaller regimes.

In this paper, we specifically concentrate on online CP detection, employing Bayesian approaches. The Bayesian framework is well-suited for quantifying our uncertainty regarding CPs using the posterior distribution, as illustrated in Ghahramani (2015). We primarily consider the class of algorithms known as Bayesian Change-Point Detection Methods (BOCPD), pioneered by Adams and MacKay (2007) as an improvement on the ideas developed by Fearnhead and Liu (2007). Since 2007, BOCPD and its extensions have found applications in various financial settings. Most applications have focused on stock returns, as demonstrated in Fan and Mackey (2017), Lleo et al. (2020) and Zhao et al. (2022), and more recently, Lleo et al. (2022, 2023) utilized BOCPD as an exit-entry model for long-short prediction in the stock market. The work by Fan and Mackey (2017) extended the BOCPD approach to the multi-sequence setting to analyze changes in 401 U.S. stocks within the S&P 500 index. BOCPD utilizes a message-passing algorithm to recursively compute the posterior distribution of the time since the last CP, termed the ‘run length’. This elapsed time is continuously updated upon receiving new data points. To perform online inference, the underlying predictive model (UPM) is computed, representing the distribution of data given the current run length. For instance, the UPM may assume a Gaussian model with different means across regimes.

In the BOCPD model, the data is assumed to be independently and identically distributed (i.i.d.) within each regime. However, this assumption is unrealistic for most financial time series. In this work, we extend BOCPD to accommodate Markovian data within each regime. While Xuan and Murphy (2007) consider the correlation structure of multivariate time series, their CP detection algorithm is offline. To the best of our knowledge, this is the first work that combines an online learning algorithm for CPs with a Markovian data structure. Furthermore, we propose a second extension of the BOCPD algorithm that relaxes the assumption of constant parameters within a regime, allowing for time-varying autocorrelation. To achieve this, we employ the class of Score Driven models introduced by Creal et al. (2013) and Harvey (2013), which provide an observation-driven framework for real-time learning of time-varying parameters. Thus, our newly proposed method combines the online CP detection approach of BOCPD with the online learning of time-varying autocorrelation parameters within each regime.

We present an empirical application of the proposed methods using order flow samples from stocks traded on NASDAQ. In an out-of-sample forecasting exercise, we find that the Score-Driven-based method outperforms other models, including autocorrelated time series models without regimes. Our analysis demonstrates that the model is correctly specified, and the residuals within each regime exhibit no correlation. By investigating the price dynamics during identified regimes, we discover that they follow concave functions of time, with the total price change in a regime also exhibiting a concave relationship with volume. These findings resemble those observed for real metaorders, consistent with the square root impact law. Finally, we demonstrate that knowledge of order flow regimes can be effectively utilized to improve predictions of order flow and price dynamics. We accomplish this by exploiting the well-known correlation between order flow and simultaneous/future price changes through market impact.

The paper is organized as follows: In section 2, we present the dataset, the variable of interest, and provide the motivation for applying a regime shift model to order flow time series. Section 3 covers the methodological aspects of the paper, outlining the main properties of BOCPD and introducing the two proposed extensions. In section 4, we describe the estimation results of the models on NASDAQ data and analyze the obtained findings. We examine the average price dynamics during an order flow regime and quantify the relationship between the total price change and the net order flow exchanged within a regime. Moreover, it is presented the forecasting analysis of order flow and the correlation with price dynamics through market impact. Finally, in section 5, we draw conclusions and offer suggestions for further research.

2. Data set and motivation

In this paper, we consider the order flow of trades and in particular the aggregated signed volume. Let M be the number of trades in a given day and let $v_i$ ($i=1,\dots ,M$) the signed volume (positive for buyer initiated and negative for seller initiated trades) of the i-th trade and let us indicate with N the number of trades we aggregate, so that our time series is composed by $T=\lfloor M/N\rfloor$ observations. The time series of interest is the aggregated order flow $x_t$ on the interval $\mathbb{N}\cap \left[N\right(t-1)+1,N(t-1)+N]$ given by (1) $x_t=\sum _{j=1}^N{v}_{N(t-1)+j},t=1,\dots ,T.$(1) Our data set consists of executed orders during March 2020 for Microsoft Corp. (MSFT) and of December 2021 for Tesla Inc. (TSLA). In order to investigate the role of aggregation time scale, we choose N = 240 and N = 730 for TSLA and N = 400 and N = 1200 for MSFT, which, in both cases correspond to an average time interval of 1 and 3 minutes, respectively. The length of the two time series is 8686 data points for 1 minute and 2856 data points for 3 minutes of TSLA and 8723 data points for 1 minute and 2908 data points for 3 minutes of MSFT. The choice of 1 and 3 minutes is of course arbitrary. However we decided to choose these interval lengths in order to avoid any microstructure noise but being still at high frequency.

To motivate our analysis, figure 1 shows the autocorrelation function of the order flow for TSLA and MSFT at the 3 minutes aggregation scale. Consistently with the literature (see, e.g. Lillo and Farmer 2004, Bouchaud et al. 2018), we observe that the autocorrelation function in all the cases has a slow decay. i.e. a buy (sell) trade is more likely followed by a buy (sell) trade. More quantitatively, it has been documented in many papers (see Bouchaud et al. 2004, Lillo and Farmer 2004, Bouchaud et al. 2006, 2009, Yamamoto and Lebaron 2010, Tóth et al. 2015, Taranto et al. 2018a, 2018b) that the autocorrelation function of trade signs $\rho \left(\tau \right)$ decays asymptotically as a power law with an exponent smaller than one $\rho \left(\tau \right)\sim \frac{1}{{\tau }^{\gamma }}$ where $\gamma <1$ which implies that the time series is long memory with a Hurst exponent $H=1-\gamma /2>1/2$.

Figure 1. Autocorrelation function of the order flow aggregated at the 3 minute time scale corresponding to N = 730 executions for TSLA (a), and N = 1200 executions for MSFT (b).

The origin of this large persistence has been investigated both empirically and theoretically. Making use of labeled data allowing to identify the market member initiating each trade, Tóth et al. (2015) empirically showed that the long-range persistence observed at the London Stock Exchange is strongly driven by order splitting, i.e. the same trader sequentially placing trades with the same sign, very likely as part of an optimal execution program. On the contrary, herding, i.e. groups of investors trading in the same direction in the same period, plays a much minor role.

From a theoretical point of view, the connection between order splitting and long memory of order flow has been elucidated by Lillo et al. (2005) (LMF). They proposed a simple model that postulates that market participants who intend to execute large orders split them into smaller orders and trade them incrementally. The large orders are called metaorders and the small trades in which they are split and sequentially traded are termed child orders. Under the assumption that metaorders are randomly sampled from a size distribution $p_L$ ($L\in \mathbb{N})$, the LMF model predicts the form of the autocorrelation function of trade signs. In particular, when the distribution $p_L$ is a power law (2) $p_L=\frac{\alpha }{L^{1+\alpha }}$(2) then the autocorrelation function of trade signs decays asymptotically as (3) $\rho \left(\tau \right)\sim \frac{1}{{\tau }^{\alpha -1}}$(3) i.e. the model predicts that $\gamma =1-\alpha$. Interestingly, very recently Sato and Kanazawa (2023) theoretically showed that predictions of the LMF model remain valid also when there is heterogeneity in trading frequency and size distribution across market participants.

The empirical validation of the LMF model poses some challenges because of the need for complete information on metaorders traded in the market. Lillo et al. (2005) used off-market trades as a proxy of metaorders. An alternative to such proxy is suggested by Vaglica et al. (2008), Moro et al. (2009) and Vaglica et al. (2010), who proposed segmentation algorithms and Hidden Markov Models to identify metaorders from brokerage data. Without relying on noisy proxies, Bershova and Rakhlin (2013) and Zarinelli et al. (2015) used private data of real metaorders by financial companies to test the power law hypothesis. In such a case, the results lacks generality since information on metaorders is company-specified. Recently, Sato and Kanazawa (2023) used account-level data of the whole Tokyo Stock Exchange to directly test the LMF model. The predicted relationship between the exponent γ and α has been very accurately verified, both at the market and at the single stock level.

Summarizing, the LMF model proposes that most of the autocorrelation of order flow comes from the execution of metaorders but, due to the anonymous nature of financial markets, their presence cannot be easily and directly inferred from public market data. However, the start of the execution of a new metaorder should lead to a regime change in the order flow time series, which could be detected with suitable statistical methods. Thus the identification of a CP in the order flow might signal the arrival of a new metaorder execution. Moreover, from the market (or econometric) point of view, the identification of a CP modifies the forecasting of future order flow and price dynamics, since only past data after the last CP are useful for predictions. The practical use of such methods for CP detection requires that they work online, i.e. that the identification is done in real time and not ex-post (such as in the segmentation algorithms used, for example, in Vaglica et al. (2008)). Finally, the method should allow to perform one-step ahead prediction in an online mode.

3. Methods: Bayesian online change-point detection algorithms

In the following, we briefly review the BOCPD algorithm, introduced by Adams and MacKay (2007), and present two novel extensions that we propose here, namely the Markovian BOCPD (MBO) and the Markovian BOCPD for Correlated data (MBOC) algorithm. The original BOCPD algorithm relies on the assumption that the data are independent and identically distributed within each regime. To relax such a strong assumption, we propose a new MBO algorithm that considers Markovian dynamics within each regime, thus allowing for serial correlation. At this stage, both BOCPD and MBO assume that parameters (i.e. mean, variance, autocorrelation) are constant within each regime. We further relax this assumption in a generalization of the MBO, named MBOC, that accounts for time-varying correlations. Such a generalization is based on the Score Driven approach introduced by Creal et al. (2013). In a companion paper Tsaknaki et al. (2024), we introduce in full generality the new class of regime-shift score-driven models, where other parameters can change over time within each regime.

3.1. The BOCPD algorithm

The BOCPD algorithm for i.i.d. data has been introduced by Adams and MacKay (2007). Let $x_{1:T}=\{x_1,\dots ,x_T\}$ be a sample time series. The model assumes that data are non-stationary (because of regimes) and satisfy the product partition model (PPM), see Barry and Hartigan (1992), meaning that data can be partitioned into regimes. Moreover, the parameters ${\theta }_R$ within each regime R are i.i.d. random variables drawn from some given distribution. Such a distribution needs to belong to the exponential family. Throughout the paper, we consider normal distributions. This assumption is tested using a Jarque-Bera (JB) statistic in the empirical application as shown in section 4.1.

Adams and MacKay (2007), assume a time series as the realization of i.i.d. random variables from a normal distribution with unknown mean ${\theta }_R$ and known variance ${\sigma }^2$, (4) $x_i\sim \mathcal{N}({\theta }_R,{\sigma }^2).$(4) Regimes and CPs separating them are not directly observable but must be inferred from data. To this end, the goal is to infer the elapsed time since the last CP, a quantity named run length and defined as follows.

Definition 3.1

The run length $r_t$ is a non-negative discrete variable defined as: (5) $r_t=\{\begin{array}{ll}0,& \text{if a CP occurs at time }t\\ r_{t-1}+1,& \mathrm{else}.\end{array}$(5)

In the BOCPD algorithm, the arrival of a CP is modeled as a Bernoulli process² with hazard rate 1/h: (6) $p\left(r_t|r_{t-1}\right)=\{\begin{array}{ll}1/h,& \mathrm{if}{r}_t=0\\ 1-1/h,& \mathrm{if}{r}_t=r_{t-1}+1\\ 0,& \mathrm{otherwise}.\end{array}$(6) The primary quantity of interest is the computation of the run length posterior $p\left(r_t|x_{1:t}\right)$ which characterizes probabilistically the number of time steps since the last CP given the data observed so far, (7) $p\left(r_t|x_{1:t}\right)=\frac{p(r_t,x_{1:t})}{p\left(x_{1:t}\right)}.$(7) The joint distribution over both the run length and the observed data can be written recursively, (8) $\begin{array}{rl}p(r_t,x_{1:t})& =\sum _{r_{t-1}}p(r_t,r_{t-1},x_t,x_{1:t-1})\end{array}$(8) (9) $\begin{array}{rl}& =\sum _{r_{t-1}}\underset{\mathrm{UPM}}{\underbrace{p(x_t|r_{t-1},x_{1:t-1})}}\underset{\mathrm{Hazard}}{\underbrace{p\left(r_t|r_{t-1}\right)}}\underset{\mathrm{Message}}{\underbrace{p(r_{t-1},x_{1:t-1})}}.\end{array}$(9) An important assumption that simplifies the computation is about the changepoint prior, namely $r_t$ is conditionally dependent on $r_{t-1}$ only. The quantity $p\left(x_{1:t}\right)$ is named evidence and is computed as (10) $p\left(x_{1:t}\right)=\sum _{r_t}p(r_t,x_{1:t}).$(10) The Underlying Predictive Model (UPM) is defined as the predictive posterior distribution given the current run length. Because of the assumption on PPM, such a distribution depends only on the last $r_{t-1}$ observations and can be stated in a more compact form as (11) $p(x_t|r_{t-1},x_{1:t-1})=p\left(x_t|x_{t-1}^{\left(r_{t-1}\right)}\right)$(11) where (12) $x_{t-1}^{\left(r_{t-1}\right)}=x_{t-r_{t-1}:t-1}\mathrm{and}{x}_{t:t-1}=\mathrm{\varnothing }.$(12) By using the conjugacy property of the exponential family when data are i.i.d., one obtains closed-form solutions for the UPM term, see Diaconis and Ylvisaker (1979) and Murphy (2007). For normal distributions, the UPM term is (13) $p\left(x_t|x_{t-1}^{\left(r_{t-1}\right)}\right)=\mathcal{N}({\mu }_{r_{t-1}},{\sigma }^2+{\sigma }_{r_{t-1}}^2),$(13) with posterior parameters given by (14) $\begin{array}{rl}{\mu }_{r_{t-1}}& =\frac{\frac{\sum _{i=t-r_{t-1}}^{t-1}{x}_i}{{\sigma }^2}+\frac{{\mu }_0}{{\sigma }_0^2}}{\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_{t-1}}^2={(\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_{t-1}\in \{1,\dots ,t-1\}.\end{array}$(14) Let us stress that the run length in equation (14(14) $\begin{array}{rl}{\mu }_{r_{t-1}}& =\frac{\frac{\sum _{i=t-r_{t-1}}^{t-1}{x}_i}{{\sigma }^2}+\frac{{\mu }_0}{{\sigma }_0^2}}{\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_{t-1}}^2={(\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_{t-1}\in \{1,\dots ,t-1\}.\end{array}$(14) ) is a latent variable we must infer. As such, the value of the posterior parameters varies depending on $r_{t-1}$, i.e. where we put the last change point, whose probability is in equation (7(7) $p\left(r_t|x_{1:t}\right)=\frac{p(r_t,x_{1:t})}{p\left(x_{1:t}\right)}.$(7) ).

The BOCPD algorithm works as follows. At time t = 0, we initialize the prior values ${\mu }_0,{\sigma }_0^2$ and the known variance ${\sigma }^2$ (see below for details). At a generic time t>0, a new data point $x_t$ becomes available, and the UPM in equation (13(13) $p\left(x_t|x_{t-1}^{\left(r_{t-1}\right)}\right)=\mathcal{N}({\mu }_{r_{t-1}},{\sigma }^2+{\sigma }_{r_{t-1}}^2),$(13) ) is computed for any possible ${\mu }_{r_{t-1}}$ and ${\sigma }_{r_{t-1}}^2$, as a function of the run length $r_{t-1}$ that takes value from 0 to t−1. Then, the joint distribution over both the run length and the observed data point, see equation (8(8) $\begin{array}{rl}p(r_t,x_{1:t})& =\sum _{r_{t-1}}p(r_t,r_{t-1},x_t,x_{1:t-1})\end{array}$(8) ), is computed for all the possible values of the run length. Thus we obtain:

1. the growth probabilities, (15) $p(r_t=l,x_{1:t}),\mathrm{for}l=1,\dots ,t;$(15)

2. the CP probability (16) $p(r_t=0,x_{1:t}).$(16)

After computing the evidence in equation (10(10) $p\left(x_{1:t}\right)=\sum _{r_t}p(r_t,x_{1:t}).$(10) ), the run length posterior is obtained by equation (7(7) $p\left(r_t|x_{1:t}\right)=\frac{p(r_t,x_{1:t})}{p\left(x_{1:t}\right)}.$(7) ). Finally, ${\mu }_{r_t}$ and ${\sigma }_{r_t}^2$ are updated as in equation (14(14) $\begin{array}{rl}{\mu }_{r_{t-1}}& =\frac{\frac{\sum _{i=t-r_{t-1}}^{t-1}{x}_i}{{\sigma }^2}+\frac{{\mu }_0}{{\sigma }_0^2}}{\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_{t-1}}^2={(\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_{t-1}\in \{1,\dots ,t-1\}.\end{array}$(14) ) in order to be used at the next time t + 1.

3.2. The MBO algorithm

The assumption about the independence of data is clearly restrictive. Here we introduce the MBO algorithm as an extension of the BOCPD algorithm to the case of Markovian dependence. Similarly to before, we consider normally distributed data. As such, within a regime R, a time series is a realization of an AR(1) process with normal innovations, (17) $\begin{array}{rl}& x_t\sim \mathcal{N}({\theta }_R,{\sigma }^2),\end{array}$(17) (18) $\begin{array}{rl}& x_t|x_{t-1}\sim \mathcal{N}({\theta }_R+\rho (x_{t-1}-{\theta }_R),{\sigma }^2(1-{\rho }^2\left)\right).\end{array}$(18) As in the previous model, in each regime, the unconditional distribution is normal with unknown mean ${\theta }_R$ and known variance ${\sigma }^2$. Moreover the conditional distribution is normal with constant correlation $\rho =\frac{\mathrm{Cov}(x_t,x_{t-1})}{{\sigma }^2}$. In the next Section, we introduce a further generalization, allowing ρ to be time-varying within each regime³.

The key observation is that the conjugacy property still holds when the data are Markovian since the conditional distribution of any member in the exponential family is still in the family, see Wainwright and Jordan (2008). After some computations, one can obtain a closed form for the UPM term as (19) $\begin{array}{rl}p\left(x_{t+1}|x_t^{\left(r_t\right)}\right)& =\mathcal{N}({\mu }_{r_t}+\rho (x_t-{\mu }_{r_t}),\\ & {\sigma }^2(1-{\rho }^2)+{\sigma }_{r_t}^2(1-\rho {)}^2),\end{array}$(19) where the posterior parameters are (20) $\begin{array}{rl}{\mu }_{r_t}& =\frac{b_{r_t}+\frac{{\mu }_0}{{\sigma }_0^2}}{a_{r_t}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_t}^2={(a_{r_t}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_t\in \{1,\dots ,t\}\end{array}$(20) and (21) $\begin{array}{rl}{a}_{r_t}& =\frac{1}{{\sigma }^2}+\frac{(r_t-1)(1-\rho {)}^2}{{\sigma }^2(1-{\rho }^2)}\end{array}$(21)

(22) $\begin{array}{rl}{b}_{r_t}& =\{\begin{array}{ll}{\displaystyle \frac{x_t}{{\sigma }^2}},& \mathrm{for}{r}_t=1\\ {\displaystyle \frac{x_{t-1}}{{\sigma }^2}}+{\displaystyle \frac{(1-\rho )(x_t-\rho x_{t-1})}{{\sigma }^2(1-{\rho }^2)}},& \mathrm{for}{r}_t=2\\ {\displaystyle \frac{x_{t+1-r_t}}{{\sigma }^2}}+{\displaystyle \frac{(1-\rho {)}^2\sum _{i=t+2-r_t}^{t-1}{x}_i+(1-\rho \left)\right(x_t-\rho x_{t+1-r_t})}{{\sigma }^2(1-{\rho }^2)}},& \mathrm{for}{r}_t\in \{3,\dots ,t\}\end{array}.\end{array}$(22)

Let us notice that we recover the previous case when $\rho =0$.

3.3. The MBOC algorithm

Both the baseline model and its Markovian generalization assume that the parameters within each regime are constant. This might be unrealistic in many empirical cases. For example, heteroscedasticity, i.e. time-varying variance, is ubiquitous in financial time series. Since our interest here mostly focuses on the correlation of the order flow and its temporal dependencies, we consider a model where temporal correlation (i.e. ρ) is time-varying within the regime. This might capture the variability and temporal persistence of trading volume, which in turn depends on the available liquidity of the market.

Time-varying parameters models display typically some difficulties for estimation. Following Cox (1981), we consider the class of observation-driven models where the parameters are unconditionally random variables, but evolve in time based on some nonlinear deterministic function of past observations.

In particular, we consider the class of Score-Driven models introduced by Creal et al. (2013) and Harvey (2013), which assume that the dynamics of the time-varying parameter(s) is autoregressive with an innovation term depending on the so-called score⁴. The score is then re-scaled by the inverse of the Fisher matrix⁵ , which is used to modulate the importance of the innovation according to the concavity of the log-likelihood. The intuition is simple: the scaled score adjusts the value of the parameter(s) in order to maximize the likelihood of the observed data. It is worth noticing that many standard models in financial econometrics, such as the GARCH, ACD, MEM, etc., are special cases of score-driven models (see www.gasmodel.com for more details).

We extend the MBO model by promoting the correlation coefficient ρ to a time-varying parameter ${\rho }_t$ described by the Score-Driven version of the AR(1) process (see Blasques et al. (2014)). We name such an extension as MBOC. We then introduce an online method to estimate both the time-varying parameter ${\rho }_t$ and the regime characteristics, namely the mean ${\theta }_R$ characterizing the regime and the run length $r_t$.

More specifically, within a regime R, the data generating process is assumed to be (23) $x_t={\rho }_t(x_{t-1}-{\theta }_R)+{\theta }_R+u_t,u_t\sim \mathcal{N}(0,{\sigma }^2),$(23) where ${\theta }_R$ and ${\sigma }^2$ are unknown. According to the Score-Driven AR(1) process, the time-varying correlation ${\rho }_t$ is described by the recursive relation⁶ (24) ${\rho }_t=\omega +\alpha s_{t-1}+\text{β}{\rho }_{t-1}$(24) where $s_t$ is the scaled score defined as (25) $\begin{array}{rl}& s_t={\mathcal{I}}_{t|t-1}^{-d}\cdot {\mathrm{\nabla }}_t,d\in [0,1]\end{array}$(25) (26) $\begin{array}{rl}& {\mathrm{\nabla }}_t=\frac{\mathrm{\partial }\log p_u\left(u_t\right)}{\mathrm{\partial }{\rho }_t}\end{array}$(26) (27) $\begin{array}{rl}& {\mathcal{I}}_{t|t-1}={\mathbb{E}}_{t|t-1}\left[{\mathrm{\nabla }}_t^T{\mathrm{\nabla }}_t\right]\end{array}$(27) and $u_t$ is the prediction error associated with the observation $x_t$. The parameter d controls the role of the Fisher matrix in modulating the role of the score in the update of the parameters. Standard choices for d are $d=0,1/2,1$. Each of them lead to different dynamic scale models. The vector of parameters $\overrightarrow{\lambda }=[\omega ,\alpha ,\text{β},{\sigma }^2{]}^{\prime }$ is estimated through a Maximum Likelihood Estimation method. In the analysis below, we set d = 0, i.e. we consider a not re-scaled score. It is (28) $s_t={\mathrm{\nabla }}_t=\frac{u_t}{{\sigma }^2}(x_{t-1}-{\theta }_R).$(28) Then the UPM term becomes (29) $p\left(x_{t+1}|x_t^{\left(r_t\right)}\right)=\mathcal{N}({\mu }_{r_t}+{\rho }_t(x_t-{\mu }_{r_t}),{\sigma }^2+{\sigma }_{r_t}^2),$(29) where the posterior parameters are (30) $\begin{array}{rl}{\mu }_{r_t}& =\frac{b_{r_t}+\frac{{\mu }_0}{{\sigma }_0^2}}{a_{r_t}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_t}^2={(a_{r_t}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_t\in \{1,\dots ,t\},\end{array}$(30) and (31) $\begin{array}{rl}{a}_{r_t}& =\frac{1}{{\sigma }^2}+\frac{(r_t-1)(1-{\rho }_t{)}^2}{{\sigma }^2(1-{\rho }_t^2)}\end{array}$(31)

(32) $\begin{array}{rl}{b}_{r_t}& =\{\begin{array}{ll}{\displaystyle \frac{x_t}{{\sigma }^2}},& \mathrm{for}{r}_t=1\\ {\displaystyle \frac{x_{t-1}}{{\sigma }^2}}+{\displaystyle \frac{(1-{\rho }_t)(x_t-{\rho }_t{x}_{t-1})}{{\sigma }^2(1-{\rho }_t^2)}},& \mathrm{for}{r}_t=2\\ {\displaystyle \frac{x_{t+1-r_t}}{{\sigma }^2}}+{\displaystyle \frac{(1-{\rho }_t{)}^2\sum _{i=t+2-r_t}^{t-1}{x}_i+(1-{\rho }_t\left)\right(x_t-{\rho }_t{x}_{t+1-r_t})}{{\sigma }^2(1-{\rho }_t^2)}},& \mathrm{for}{r}_t\in \{3,\dots ,t\}\end{array}.\end{array}$(32)

The vector of parameters $\overrightarrow{\lambda }$ is estimated at each time-step within the time window associated with the most likely regime after we demean the data with the posterior mean see equation (30(30) $\begin{array}{rl}{\mu }_{r_t}& =\frac{b_{r_t}+\frac{{\mu }_0}{{\sigma }_0^2}}{a_{r_t}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_t}^2={(a_{r_t}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_t\in \{1,\dots ,t\},\end{array}$(30) ). In particular at each time step t>1 we find $i={\mathrm{argmax}}_{i\in \{1,\dots ,t\}}p(r_t=i|x_{1:t})$ and we consider the demeaned data set $x_t^{\left(i\right)}-{\mu }_i=\{x_{t+1-i}-{\mu }_i,\dots ,x_t-{\mu }_i\}$ in which we infer $\overrightarrow{\lambda }$ and we filter ${\rho }_t$ with the use of the Score-Driven model. In order to robustify the algorithm and accomplish better Mean Squared Error (MSE) (see section 4.1 for more details) we define a threshold value η, then we filter the time-varying correlation and we infer the variance (see algorithm 1 for more details on the inference procedure of the MBOC model) whenever $i>\eta$. In that way, the demeaned data set contains at least η data points. The threshold value is a hyperparameter that is tuned in a preliminary phase, see below for the implementation details.

4. Results

4.1. Model estimation and empirical analysis

We estimate the three models, BOCPD, MBO, and MBOC, on the time series of aggregated order flow $x_t$ of TSLA and MSFT, both for the 1 and the 3 min aggregation time scale. The application of the models requires a careful choice of several hyperparameters. For both TSLA and MSFT, the prior value of the mean for all models is set to ${\mu }_0=0$ shares due to the symmetry between buy and sell orders. The tuning of the other hyperparameters is obtained by minimizing the MSE in the first day of each month and each stock. For TSLA, for all the models, the prior value of the variance of the mean is set to ${\sigma }_0^2={10}^7$ shares${}^2$, while for MSFT ${\sigma }_0^2=15\cdot {10}^7$ shares${}^2$. The hazard rate is set to $h=\frac{1}{80}$. For both BOCPD and MBO the known variance is set to ${\sigma }^2={10}^8$ shares${}^2$ for TSLA while for MSFT ${\sigma }^2=15\cdot {10}^8$ shares${}^2$. The same values for TSLA and MSFT are being used for the initial variance ${\sigma }_i^2$ of the MBOC model (see algorithm 1). Moreover, the initial correlation ${\rho }_1$ is set to 0.2 for TSLA and 0.3 for MSFT while the initial parameters of the Score-Driven dynamics are set to $\overrightarrow{\lambda }=[0.08,0.02,0.05,{10}^8{]}^{\prime }$ and the η value is set equal to 20 shares for the 1 minute and 10 shares for the 3 minute data set. The hyperparameter η is being tuned in such a way in order to obtain the best MSE. Finally, for the constant correlation coefficient ρ of MBO we have tested different specifications. In the table below we will report the results for three of them, showing that the predictive capacity of the model slightly depend on it.

4.1.1. Model comparison

We present here the results of an online prediction study for order flow data by using BOCPD, MBO, and MBOC models introduced in section 3, and we compare their performances by computing the Mean Squared Error for the predictive mean of each model. The three models are then compared with the ARMA(1,1) model, estimated on the whole time period by assuming the absence of regimes. As such, the ARMA(1,1) model represents a natural benchmark to test whether including regime-switching dynamics does improve the forecasting of order flow.

The predictive mean ${\hat{\mu }}_t$ at time t as one step ahead forecast is defined by using observations up to time t and to predict out-of-sample the realization at time t + 1. That is (33) $\begin{array}{rl}{\hat{\mu }}_t& =\sum _{r_t}p(x_{t+1}|x_{1:t},r_t,)p\left(r_t|x_{1:t}\right)\\ & =\sum _{r_t}{\mu }_{r_t}p\left(r_t|x_{1:t}\right)\end{array}$(33) where ${\mu }_{r_t}$ is defined in equations (14(14) $\begin{array}{rl}{\mu }_{r_{t-1}}& =\frac{\frac{\sum _{i=t-r_{t-1}}^{t-1}{x}_i}{{\sigma }^2}+\frac{{\mu }_0}{{\sigma }_0^2}}{\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_{t-1}}^2={(\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_{t-1}\in \{1,\dots ,t-1\}.\end{array}$(14) ), (20(20) $\begin{array}{rl}{\mu }_{r_t}& =\frac{b_{r_t}+\frac{{\mu }_0}{{\sigma }_0^2}}{a_{r_t}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_t}^2={(a_{r_t}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_t\in \{1,\dots ,t\}\end{array}$(20) ), and (30(30) $\begin{array}{rl}{\mu }_{r_t}& =\frac{b_{r_t}+\frac{{\mu }_0}{{\sigma }_0^2}}{a_{r_t}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_t}^2={(a_{r_t}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_t\in \{1,\dots ,t\},\end{array}$(30) ) for BOCPD, MBO and MBOC model respectively. Then, the MSE can be computed as (34) $\mathrm{MSE}=\frac{1}{T}\sum _{t=1}^T({\hat{\mu }}_{t-1}-x_t{)}^2.$(34) Table 1 shows the MSE of the three aforementioned models along with the ARMA(1,1) for both TSLA and MSFT stocks. As mentioned above, we consider three different values of the correlation coefficient ρ in the MBO model. We observe that the MBOC model outperforms all competitors. In particular, the proposed models (MBO and MBOC) outperform systematically the ARMA(1,1) benchmark, while the baseline BOCPD model displays comparable performances. Finally, notice that the role of the hyperparameter ρ for the MBO model is relatively marginal.

Table 1. Comparison of out of sample one-step-ahead MSE of ARMA(1,1), BOCPD, MBO and MBOC.

	TSLA			MSFT
ρ	0.1	0.2	0.3	0.2	0.3	0.4
ARMA	0.908	0.908	0.908	0.860	0.860	0.860
BOCPD	0.907	0.907	0.907	0.878	0.878	0.878
MBO	0.895	0.896	0.911	0.835	0.834	0.844
MBOC	0.890	0.890	0.890	0.831	0.831	0.831

Notes: The correlation ρ is the one used in MBO. Data refer to TSLA and MSFT at the 3 minute resolution.

In conclusion, the online prediction study with order flow data suggests that regime-switching models accounting for a Markovian correlation structure outperform both the baseline BOCPD model and the benchmark. The MBOC model displays the best forecasting performance and high flexibility in data description. In the following Sections, we exploit such flexibility in modeling regime-switching dynamics in the presence of time-varying correlations to empirically show a clear connection between regimes for aggregated order flows and the market impact of associated trades (likely including metaorders).

4.1.2. Empirical analysis of identified regimes

Here we investigate the statistical properties of the identified regimes for the aggregated order flows. We consider the MBOC model because of the best performances. Let us first introduce the adopted definition of an identified regime.

Definition 4.1

Let $x_{1:T}$ be a time series and $t,s\in \mathbb{N}\cap [1,T]$ times with t<s such that (35) $\begin{array}{rl}& {\text{argmax}}_{i\in \{0,1,\dots ,t\}}p(r_t=i|x_{1:t})\\ & ={\text{argmax}}_{i\in \{0,1,\dots ,s\}}p(r_s=i|x_{1:s})=0\end{array}$(35) and for any $u\in \mathbb{N}\cap (t,s)$, ${\text{argmax}}_{i\in \{0,1,\dots ,u\}}p(r_u=i|x_{1:u})\ne 0$. Then the subset $x_{t:s-1}$ of the time series $x_{1:T}$ is defined as a regime.

The top panel of figures 2 and 3 shows $x_t$ for TSLA (MSFT) aggregated every 730 (1200) executions corresponding to an average time interval of 3 minutes. The vertical red dashed lines indicate the CPs identified by MBOC, according to the definition above. Interestingly, many CPs are observed at the end of a trading day for both stocks. On one side this is expected since overnight is a natural separation between regimes, but on the other side, this is an indication that the proposed method is able to identify regime changes.

Figure 2. Top: Order flow of Tesla at the 3 minutes aggregation time scale during December 2021. The black solid line is the predictive mean and the dashed line, the predictive standard deviation. The red dashed lines indicate the CPs the MBOC finds and as a result the regimes. The tick labels on the x-axis indicate the end of each trading day. Bottom: Run length posterior of the MBOC model. The darkest the color, the highest the probability of the run length. The red line highlights the most likely path, i.e. value of $r_t$ with the largest run length posterior $p\left(r_t|x_{1:t}\right)$ for each t.

Figure 3. Top: Order flow of Microsoft at the 3 minutes aggregation time scale during March 2020. The black solid line is the predictive mean and the dashed line, the predictive standard deviation. The vertical red dashed lines indicate the CPs the MBOC finds and as a result the regimes. The tick labels on the x-axis indicate the end of each trading day. Bottom: Run length posterior of the MBOC model. The darkest the color, the highest the probability of the run length. The red line highlights the most likely path, i.e. value of $r_t$ with the largest run length posterior $p\left(r_t|x_{1:t}\right)$ for each t.

The bottom panels of figures 2 and 3 show the run length posterior of the MBOC model for the two assets. For each time (on the abscissa) the vertical axis displays in grayscale the probability that the run length has a given value (on the ordinate). Darker gray regions correspond to higher probabilities. The red line highlights the most likely path, i.e. the value of $r_t$ with the largest run length posterior $p\left(r_t|x_{1:t}\right)$ for each t. Finally, we also show the (one step ahead) predictive standard deviation defined as (36) ${\hat{\sigma }}_t=\sqrt{\sum _{r_t}{\sigma }_{r_t}^{2}p\left(r_t|x_{1:t}\right)},$(36) where ${\sigma }_{r_t}^2$ is as in equations (14(14) $\begin{array}{rl}{\mu }_{r_{t-1}}& =\frac{\frac{\sum _{i=t-r_{t-1}}^{t-1}{x}_i}{{\sigma }^2}+\frac{{\mu }_0}{{\sigma }_0^2}}{\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_{t-1}}^2={(\frac{r_{t-1}}{{\sigma }^2}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_{t-1}\in \{1,\dots ,t-1\}.\end{array}$(14) ), (20(20) $\begin{array}{rl}{\mu }_{r_t}& =\frac{b_{r_t}+\frac{{\mu }_0}{{\sigma }_0^2}}{a_{r_t}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_t}^2={(a_{r_t}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_t\in \{1,\dots ,t\}\end{array}$(20) ), and (30(30) $\begin{array}{rl}{\mu }_{r_t}& =\frac{b_{r_t}+\frac{{\mu }_0}{{\sigma }_0^2}}{a_{r_t}+\frac{1}{{\sigma }_0^2}}\mathrm{and}{\sigma }_{r_t}^2={(a_{r_t}+\frac{1}{{\sigma }_0^2})}^{-1}\\ & \mathrm{for}{r}_t\in \{1,\dots ,t\},\end{array}$(30) ) for BOCPD, MBO, and MBOC models, respectively.

Regime length distribution : For the 1 (3) minute(s) data set, we find 911 (546) regimes for TSLA and 1394 (690) regimes for MSFT. Figure 4 shows the histograms of the length of the detected regimes. Consistently with the constant hazard function, we find that the regime length is approximately exponentially distributed with a mean regime length of 10 for the 1 minute data set and 5 for the 3 minutes data set intervals for TSLA and 7 for the 1 minute data set and 4 intervals for the 3 minutes data set of MSFT corresponding to 10 and 15 trading minutes for TSLA and to 7 and 12 trading minutes for MSFT respectively. The length of the regimes ranges in the interval $[1,92]$ for 1 minute and in $[1,30]$ for 3 minutes for TSLA while for MSFT in $[1,65]$ for 1 minute and in $[1,30]$ for 3 minutes for MSFT.

Figure 4. Histograms of regimes length for Tesla (a) with $\mathrm{\Delta t}=1\min$, (b) with $\mathrm{\Delta t}=3\min$ and for Microsoft (c) with $\mathrm{\Delta t}=1\min$, (d) with $\mathrm{\Delta t}=3\min$.

Gaussianity inside regimes : The main assumption of both BOCPD and MBO is that the variable $x_t$ is Gaussian within each regime, with constant parameters. For MBOC, we expect that $x_t$ is only conditionally Gaussian, because of the time-varying autocorrelation, but not unconditionally over the whole period and only approximately within a regime. The Jarque-Bera (JB) test rejects the Gaussianity hypothesis of unconditional aggregated order flow $\{x_t{\}}_{t=1,\dots ,T}$ at a $1\mathrm{\%}$ significance level, for both stocks and time scales. We then perform the JB test within each regime detected by the MBOC model. For the 3 minute time scale, we cannot reject the null hypothesis at $5\mathrm{\%}$ confidence level for $94\mathrm{\%}$ ($95\mathrm{\%}$) of the regimes for TSLA (MSFT). When we consider the 1 minute time scale, the frequency of rejection is $86\mathrm{\%}$ ($87\mathrm{\%}$). These findings support the choice of MBOC to identify regimes with order flows as approximately Gaussian within each regime.

Autocorrelation of residuals inside regimes : As a final model checking we test for the lack of serial correlation in the residuals of our model within each regime. We have seen above that, coherently with the literature, $v_t$ is strongly autocorrelated. Following Lillo et al. (2005), our assumption is that this correlation is driven by the presence of regimes, which in turn are likely associated with metaorders. We thus apply the Ljung-Box test to the residuals in each regime. For the 3 minute time scale, we cannot reject the null hypothesis of uncorrelated residuals for $98\mathrm{\%}$ ($99\mathrm{\%}$) of the regimes of TSLA (MSFT), with $5\mathrm{\%}$ confidence level. For the 1 minute time scale, the frequency of rejection is $97\mathrm{\%}$ ($93\mathrm{\%}$). It is possible to conclude that the MBOC model captures most of the serial correlation of aggregated order flow. Notice that, according to the model, the unconditional slow decay of the autocorrelation of order flow observed in the literature (see also figure 1) is due largely to regime-switching dynamics and, only partially, to Markovian temporal dependencies.

4.2. Price impact during order flow regimes

In this Section, we empirically study the average price dynamics inside a detected order flow regime and we measure the relation between the total price change and the net volume exchanged in the same regime.

4.2.1. Price as a function of time inside an order flow regime

As said, we first study the average price dynamics during an order flow regime. This type of analysis mirrors the one performed, for example, by Bacry et al. (2015) and Zarinelli et al. (2015) which studied the average price dynamics during the execution of a metaorder. Using labeled data allowing to identify when an institutional investor executed a metaorder, these papers find that (i) the average price dynamics is correlated with the conditioning metaorder sign, i.e. the price increases (decreases) when a buy (sell) metaorder is executed; (ii) the price dynamics is concave in time, i.e. the price increases faster at the beginning of a buy metaorder and slowly toward the end.

Here we take a step forward by asking what is the average price dynamics during a regime of aggregated order flow detected with the MBOC model. To this end, for each detected regime R (see definition 4.1), characterized by an initial time $t_R$ and a final time $s_R>t_R$, we denote with ${ϵ}_R=\mathrm{sign}\left(\sum _{t_R\le t<s_R}{x}_t\right)$ the sign of the order flow in the regime, being equal to $+1$ ($-1$) when the regime is dominated by the volume of buyer (seller) initiated trades. Since during a metaorder execution we expect a significant net imbalance of buy or sell volume, we will consider subsets of regimes for which (37) $Z_R:=\left|\frac{\sum _{t_R\le t<s_R}{x}_t}{\sum _{t_R\le t<s_R}|x_t|}\right|>\mathrm{\Theta },0\le \mathrm{\Theta }<1$(37) i.e. when the difference between buy and sell volume divided by their sum is larger than Θ. Notice that when $\mathrm{\Theta }=0$ the subset coincides with the entire set of regimes identified by the MBOC. Appendix A.1 reports the number of regimes in the different subsets.

Indicating with $p_t$ the log-price of the last transaction in the interval labeled by time t, we compute the log-price change between the beginning of the regime and $t_R+k$, where $0\le k<s_R-t_R$ and we take the average (38) ${\mathcal{I}}^{\mathrm{\Theta }}\left(k\right)={\mathbb{E}}_R\left[{ϵ}_R\right(p_{t_R+k}-p_{t_R-1})|t_R+k<s_R,Z_R>\mathrm{\Theta }].$(38) With ${\mathbb{E}}_R[\cdot ]$, we denote the sample average over the regimes, i.e. that $t_R$ is the first interval of a regime, and the conditioning restricts it to those regimes for which the observation at $t_R+k$ is in the same regime as the one at $t_R$ as well as to those regimes that satisfy condition in equation (37(37) $Z_R:=\left|\frac{\sum _{t_R\le t<s_R}{x}_t}{\sum _{t_R\le t<s_R}|x_t|}\right|>\mathrm{\Theta },0\le \mathrm{\Theta }<1$(37) ).

Figure 5 shows the impact function ${\mathcal{I}}^{\mathrm{\Theta }}\left(k\right)$ as a function of k for the two stocks and the two time scales when $\mathrm{\Theta }=0,0.5$, and 0.9. Error bars are standard errors in each bin. We notice that in all cases impact is positive and increasing. This is somewhat expected since we are implicitly conditioning on the sign of order flow in the whole regime, thus the observed behavior is coherent with the known correlation between aggregated order flow and contemporaneous price change, see Bouchaud et al. (2009), Patzelt and Bouchaud (2018). Interestingly the price dynamics is a concave function of time, similarly to what is observed when conditioning on metaorders execution instead of on order flow regimes. Moreover the degree of concavity increases with Θ. Clearly the concavity is not expected by the mere fact that the regime is characterized by a net order flow sign, while it could instead be explained by the Transient Impact Model of Bouchaud et al. (2004) or by the LLOB model of Donier et al. (2015), which predicts a concave average price temporal profile when the order flow has a non-zero average as during a metaorder execution.

Figure 5. Market Impact ${\mathcal{I}}^{\mathrm{\Theta }}\left(k\right)$ defined in equation (38(38) ${\mathcal{I}}^{\mathrm{\Theta }}\left(k\right)={\mathbb{E}}_R\left[{ϵ}_R\right(p_{t_R+k}-p_{t_R-1})|t_R+k<s_R,Z_R>\mathrm{\Theta }].$(38) ) for Tesla (a) with $\mathrm{\Delta t}=1$ min, (b) with $\mathrm{\Delta t}=3$ min and for Microsoft (c) with $\mathrm{\Delta t}=1$ min, (d) with $\mathrm{\Delta t}=3$ min. The plotted quantities are in basis points. Error bars are standard errors. The black dashed line indicates the y = 0 axis.

4.2.2. Price impact as a function of volume

Finally, we study the relation between the total price change in a regime and the total net volume in the same time span. A large body of empirical literature Torre (1997), Almgren et al. (2005), Moro et al. (2009), Tóth et al. (2011), Bershova and Rakhlin (2013), Zarinelli et al. (2015) and Tóth et al. (2016) have shown that on average the total price impact during a metaorder execution scales with a sublinear power law of metaorder volume, a relation well fit by a power law with an exponent ranging in $[0.4,0.7]$. This is the celebrated square root impact law. It is therefore natural to investigate empirically the relation between the same two quantities within a regime identified with our method.

To this end, defining $\mathrm{\Delta }{p}_R=p_{s_R}-p_{t_R}$ we consider the non-linear regression: (39) ${ϵ}_R\mathrm{\Delta }{p}_R=A{\left({ϵ}_R\sum _{t_R\le t<s_R}{x}_t\right)}^{\gamma }+\mathrm{noise}.$(39) Since the measurement of market impact is notoriously very noisy, we have performed the estimation both on the original dataset and on a dataset where potential outliers are removed. In the latter approach, we used the standard procedure of considering outliers datapoints corresponding to regimes whose price change is smaller (larger) than the first (third) quartile minus (plus) 1.5 times the interquartile range (see appendix A.2 for details).

Table 2 reports the estimated parameters when outliers are removed, while appendix A.2 reports the results for the entire dataset and presents the scatter plots of the data and the fitted curve. Both Tables indicate that the exponent γ is smaller than one and for the data without outliers it is remarkably close to 0.5, as postulated by the square root law.

Table 2. Estimated parameters and their standard errors (SE) of the regression of equation (39(39) ${ϵ}_R\mathrm{\Delta }{p}_R=A{\left({ϵ}_R\sum _{t_R\le t<s_R}{x}_t\right)}^{\gamma }+\mathrm{noise}.$(39) ).

	TSLA
	$\mathrm{\Delta t}=1\min$				$\mathrm{\Delta t}=3\min$
Θ	A	SE of A	γ	SE of γ	A	SE of A	γ	SE of γ
0	0.121	0.05	0.592	0.041	0.298	0.139	0.52	0.045
0.5	0.159	0.066	0.567	0.041	0.284	0.133	0.528	0.045
0.9	0.172	0.084	0.552	0.049	0.387	0.179	0.504	0.044
	MSFT
	$\mathrm{\Delta t}=1\min$				$\mathrm{\Delta t}=3\min$
Θ	A	SE of A	γ	SE of γ	A	SE of A	γ	SE of γ
0	0.498	0.209	0.4	0.038	0.365	0.203	0.458	0.048
0.5	0.469	0.209	0.408	0.04	0.314	0.185	0.47	0.05
0.9	0.334	0.144	0.444	0.038	0.332	0.199	0.47	0.051

Note: Results refer to data with outliers removed.

Clearly these results are preliminary and should be validated on larger panels of stocks, also pooling them together with the usual rescaling by daily volatility and volume. However we find these results very encouraging and suggestive of a relation between the identified regimes and the execution of metaorders.

4.3. Online prediction of order flow and market impact

The possibility of performing an online detection of regimes and regime changes in the order flow opens the question of how to use this information to predict subsequent order flow and price changes. In the PPM, regimes are independent, hence in forecasting future values only data from the current regime are useful, while older data add noise to the prediction. This idea will be used to build online predictions of order flow. Additionally, through market impact, price dynamics is correlated with order flow. Thus a proper modeling of order flow is useful to forecast future price.

Since we have seen in the last Section that order flow sign of a regime correlates with contemporaneous price change, we can ask the question of whether the knowledge that a new regime in order flow has just started allows to predict the future order flow and, more importantly, the future price dynamics. Consistently with the results of section 4.1 showing that the MBOC model outperforms the competitors (BOCPD, MBO, ARMA) in one step ahead prediction of order flow, here we focus our analyses on the regimes identified by MBOC. It is important however to stress that qualitatively similar results are obtained with the other two simpler models, BOCPD and MBO. In other words the relation between price dynamics and aggregated order flow is importantly understood by using regimes, while the choice of the specific regime shift model improves the short-term prediction of aggregated order flow. However, as shown in in the appendix 2, the MBOC method achieves higher predictability wrt the other CP detection methods.

To better understand the role of regimes in prediction, let consider the following argument. If the data generating process of order flow is truly consistent with a product partition model (i.e. independent regimes), the knowledge of the data of the previous regimes is not useful for prediction. Thus in this case, it is better to use only the data in the current regime and its learned statistical properties. However, even if we are relatively sure that a new regime has just started, its parameters could be quite uncertain at the beginning. Thus, in order to form a forecast, it is better to wait for few observations into a new regime. Moreover, if many regimes are very short (e.g. composed by one or two intervals), as observed empirically in figure 4, it might be better to build predictions after the observation of a few intervals in a regime.

Online order flow prediction : Following the above argument, we adopt the following procedure. Whenever we detect a CP in an online fashion for the time series of order flow, we measure the correlation (40) ${\mathcal{I}}_{ϵ}^{\left(1\right)}\left(k\right)={\mathbb{E}}_R\left[\mathrm{sign}\right(x_{t_R})\cdot \mathrm{sign}(x_{t_R+k}\left)\right]k=1,2,\dots$(40) where ${\mathbb{E}}_R$ indicates that we are conditioning on the fact that $t_R$ is the first observation of a new regime. Notice that (i) the correlation is extended to values of k possibly beyond the end of the detected regime (at time $s_R$); (ii) we do not consider the case k = 0 since in this case the correlation is trivially equal to 1. The superscript ${}^{\left(1\right)}$ in the above expression means that we take the sign of the aggregated order flow in the first interval (see below for an extension).

The continuous yellow line in figure 6 shows that ${\mathcal{I}}_{ϵ}^{\left(1\right)}\left(k\right)$ is a poor predictor of order flow. To better quantify this statement, the dashed yellow line in the figure is the unconditional correlation (41) ${\stackrel{\mathcal{~}}{\mathcal{I}}}_{ϵ}^{\left(1\right)}\left(k\right)=\mathbb{E}\left[\mathrm{sign}\right(x_t)\cdot \mathrm{sign}(x_{t+k}\left)\right]k=1,2,\dots$(41) which makes no use of regime detection (and for this reason the expectation does not have the subscript ${}_R$; the tilde refers to the expectation without considering regimes). Clearly, the unconditional correlation is larger than the conditional one.

As said above, one of the reasons of the comparable performance of ${\mathcal{I}}_{ϵ}^{\left(1\right)}$ with respect to ${\stackrel{\mathcal{~}}{\mathcal{I}}}_{ϵ}^{\left(1\right)}$ is the fact that there are many regimes of length one and also that the sign of the new regime, ${ϵ}_R$, might be poorly measured by the sign of the first interval $\mathrm{sign}\left(x_{t_R}\right)$. A better option is to wait few more intervals within the regime before building the predictor. Thus, defining $m=1,2,\dots$ the number of intervals in a regime we wait before forming the prediction, we introduce the correlation (42) $\begin{array}{rl}{\mathcal{I}}_{ϵ}^{\left(m\right)}\left(k\right)& ={\mathbb{E}}_R[\mathrm{sign}\left(\sum _{t=t_R}^{t_R+m-1}{x}_t\right)\cdot \mathrm{sign}(x_{t_R+m-1+k}\left)\right]\\ & k=1,2,\dots \end{array}$(42) which is the correlation between the sign of the order flow in the first m intervals of a regime and the sign of the order flow in an interval k steps after these m intervals. Notice that we are not conditioning on the fact that $t_R+m-1+k$ is in the same regime as $t_R$, so the two observation could belong to different regimes. However, we condition on the fact that $t_R+m-1$ is in the same regime as $t_R$. Similarly we use as a benchmark case the predictor (43) $\begin{array}{rl}{\stackrel{\mathcal{~}}{\mathcal{I}}}_{ϵ}^{\left(m\right)}\left(k\right)& =\mathbb{E}[\mathrm{sign}\left(\sum _{s=t}^{t+m-1}{x}_s\right)\cdot \mathrm{sign}(x_{t+m-1+k}\left)\right]\\ & k=1,2,\dots .\end{array}$(43) The orange, red, and dark red continuous lines in figure 6 show ${\mathcal{I}}_{ϵ}^{\left(m\right)}\left(k\right)$ for m = 2, 3, 4 respectively, while the corresponding dashed lines refer to ${\stackrel{\mathcal{~}}{\mathcal{I}}}_{ϵ}^{\left(m\right)}\left(k\right)$. We observe that the correlations based on regimes are larger than the corresponding ones without regimes, especially for large m. This empirical evidence indicates that the knowledge of the order flow regimes improves the short-term predictability of order flow.

Figure 6. Online order flow prediction for Tesla (a) with $\mathrm{\Delta t}=1\min$, (b) with $\mathrm{\Delta t}=3\min$ and for Microsoft (c) with $\mathrm{\Delta t}=1\min$, (d) with $\mathrm{\Delta t}=3\min$. The plotted quantities are defined in equations (42(42) $\begin{array}{rl}{\mathcal{I}}_{ϵ}^{\left(m\right)}\left(k\right)& ={\mathbb{E}}_R[\mathrm{sign}\left(\sum _{t=t_R}^{t_R+m-1}{x}_t\right)\cdot \mathrm{sign}(x_{t_R+m-1+k}\left)\right]\\ & k=1,2,\dots \end{array}$(42) ) and (43(43) $\begin{array}{rl}{\stackrel{\mathcal{~}}{\mathcal{I}}}_{ϵ}^{\left(m\right)}\left(k\right)& =\mathbb{E}[\mathrm{sign}\left(\sum _{s=t}^{t+m-1}{x}_s\right)\cdot \mathrm{sign}(x_{t+m-1+k}\left)\right]\\ & k=1,2,\dots .\end{array}$(43) ) and error bars are standard errors.

Online market impact prediction : We now consider the prediction of price change based on the knowledge of being in a regime of order flow. To this end we introduce the online impact (44) $\begin{array}{rl}{\mathcal{I}}_{\mathrm{\Delta }p}^{\left(m\right)}\left(k\right)& ={\mathbb{E}}_R[\mathrm{sign}\left(\sum _{t=t_R}^{t_R+m-1}{x}_t\right)\cdot (p_{t_R+m-1+k}-p_{t_R+m-1}\left)\right]\\ & k=1,2,\dots \end{array}$(44) which is the correlation between the sign of the total order flow in the m initial intervals of a regime and the subsequent price change over k intervals. Compared to equation (38(38) ${\mathcal{I}}^{\mathrm{\Theta }}\left(k\right)={\mathbb{E}}_R\left[{ϵ}_R\right(p_{t_R+k}-p_{t_R-1})|t_R+k<s_R,Z_R>\mathrm{\Theta }].$(38) ), two important differences are worth to be highlighted. First, the sign inside the expectation is taken only on the aggregated order flow of the m intervals used to build the prediction, while in equation (38(38) ${\mathcal{I}}^{\mathrm{\Theta }}\left(k\right)={\mathbb{E}}_R\left[{ϵ}_R\right(p_{t_R+k}-p_{t_R-1})|t_R+k<s_R,Z_R>\mathrm{\Theta }].$(38) ) ${ϵ}_R$ considers the sign of the whole regime and therefore is non-causal. Second, in ${\mathcal{I}}_{\mathrm{\Delta }p}^{\left(m\right)}\left(k\right)$ we do not condition on $s_R-t_R>k$ as in equation (38(38) ${\mathcal{I}}^{\mathrm{\Theta }}\left(k\right)={\mathbb{E}}_R\left[{ϵ}_R\right(p_{t_R+k}-p_{t_R-1})|t_R+k<s_R,Z_R>\mathrm{\Theta }].$(38) ) since after having observed m intervals in a regime we do not know when the regime is going to end. In other words, for a given k, we take the average both on cases when $t_R$ and $t_R+k$ belong to the same regime and when they do not. Finally, as before we use as a benchmark an impact predictor that is based on the sign of the order flow, (45) $\begin{array}{rl}{\stackrel{\mathcal{~}}{\mathcal{I}}}_{\mathrm{\Delta }p}^{\left(m\right)}\left(k\right)& =\mathbb{E}[\mathrm{sign}\left(\sum _{s=t}^{t+m-1}{x}_s\right)\cdot (p_{t+m-1+k}-p_{t+m-1}\left)\right]\\ & k=1,2,\dots .\end{array}$(45) When m = 1 the quantity ${\stackrel{\mathcal{~}}{\mathcal{I}}}_{\mathrm{\Delta }p}^{\left(m\right)}\left(k\right)$ becomes the response function widely investigated in the market impact literature, mainly in transaction time, see Bouchaud et al. (2004) and Bouchaud et al. (2009). We choose this more general definition in order to make a fairer comparison between impact predictors using the same number of past order flow observations.

Figure 7 shows these different quantities for online market impact prediction, considering both stocks and both time intervals. It is evident that as soon as m>1, ${\mathcal{I}}_{\mathrm{\Delta }p}^{\left(m\right)}\left(k\right)$ (orange to dark red lines) is much larger than the corresponding response function ${\stackrel{\mathcal{~}}{\mathcal{I}}}_{\mathrm{\Delta }p}^{\left(m\right)}\left(k\right)$ (dashed orange to dark red lines) which does not make use of regimes. Moreover, the larger m, the larger the correlation between the order flow sign and the future price change, in all four investigated cases. Thus the (online) knowledge that a regime has started provides a significant additional forecasting power to future price change with respect to the response function, which is an unconditional cross-correlation between current order flow and future price change.

Figure 7. Online market impact prediction for Tesla (a) with $\mathrm{\Delta t}=1$ min, (b) with $\mathrm{\Delta t}=3$ min and for Microsoft (c) with $\mathrm{\Delta t}=1$ min, (d) with $\mathrm{\Delta t}=3$ min. The plotted quantities (in basis points) are defined in equations (44(44) $\begin{array}{rl}{\mathcal{I}}_{\mathrm{\Delta }p}^{\left(m\right)}\left(k\right)& ={\mathbb{E}}_R[\mathrm{sign}\left(\sum _{t=t_R}^{t_R+m-1}{x}_t\right)\cdot (p_{t_R+m-1+k}-p_{t_R+m-1}\left)\right]\\ & k=1,2,\dots \end{array}$(44) ) and (45(45) $\begin{array}{rl}{\stackrel{\mathcal{~}}{\mathcal{I}}}_{\mathrm{\Delta }p}^{\left(m\right)}\left(k\right)& =\mathbb{E}[\mathrm{sign}\left(\sum _{s=t}^{t+m-1}{x}_s\right)\cdot (p_{t+m-1+k}-p_{t+m-1}\left)\right]\\ & k=1,2,\dots .\end{array}$(45) ) and error bars are standard errors.

5. Conclusion

In this work, we proposed the use of Bayesian Online Change Point Detection Methods to identify (in a real-time setting) regimes in time series of aggregated order flow of financial assets. Since the existing methods make very strong assumptions on the data generating process, in particular for what concerns the serial correlation of data within each regime, we proposed here two extensions of the regime detection algorithm: the first one assumes a Markovian dynamics inside each regime, while the second one makes use of an observation driven dynamics based on the score-driven mechanism. As shown by the recent econometric literature, the score driven approach is extremely flexible also as a filter of a misspecified dynamics (tantamount to GARCH). The companion paper (Tsaknaki et al. 2024) provides more methodological details of this new class of models by discussing different specifications where other parameters (e.g. the variance) are time-varying within each regime. The analysis of two liquid stocks traded in the NASDAQ market shows that the new algorithms presented here, particularly the latter, outperform the baseline model in out-of-sample forecasting. In general, we find that the regime-switching methods outperform standard econometric time series models like ARMA(1,1). Moreover, a careful model checking shows that the algorithm outputs well specified regimes both in terms of Gaussianity of data and of lack of serial correlation of residuals, within each regime.

From the financial point of view, the identification of weakly autocorrelated regimes in the order flow time series suggests that the observed unconditional long memory might be explained by regime switching. This is in line with the mechanism proposed by Lillo et al. (2005) who connected the long memory to order splitting by heterogeneous institutional investors. It is natural at this point to try to identify the detected regimes with time periods when one or a few institutional investors are trading a large order. Of course, we do not have any empirical evidence in support of this idea which, at this point, can be considered as a conjecture to be tested with suitable data (for example those used by Zarinelli et al. (2015) or Sato and Kanazawa (2023)).

The paper shows how the online identification of regimes can be used to significantly improve the forecasting of order flow and of price dynamics. Using the knowledge of the order flow during the current regime provides better predictions when compared with methods using unconditionally the past history of order flow. We foresee that such improvement could be fruitfully used in several financial applications, such as optimal trading, market making, and alpha signal detection. Similarly, if our interpretation above is correct, the online regime detection method could be used to statistically identify the execution of a large institutional execution from anonymous market data.

Disclosure statement

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

Data availability statement

The data that support the findings of this study are not publicly available but can be purchased at https://lobsterdata.com.

Code availability statement

The code that reproduces this study for section 3 can be found at https://sobigdata.d4science.org/catalogue-sobigdata?path=/dataset/score-driven_bayesian_online_change_point_detection_sd-bocpd_ or https://gitlab.com/YvTsak/ScoreDrivenBOCPDhttps://gitlab.com/YvTsak/ScoreDrivenBOCPD and for section 4 at https://sobigdata.d4science.org/catalogue-sobigdata?path=/dataset/online_learning_of_order_flow_and_market_impact_olofmi_ or https://gitlab.com/YvTsak/olofmihttps://gitlab.com/YvTsak/olofmi.

Additional information

Funding

This paper is funded by the European Union Next Generation EU with grant number PNRR IR0000013 ‘SoBigData.it’.

Notes

1 I.e. the sequence of signed trade volume, positive (negative) when buyer (seller) initiated.

2 Other assumptions on the distribution of regimes can be implemented, for example, those leading to a non-exponential distribution of regime length. This more realistic extension is left for future research.

3 We have also explored a further generalization where the variance inside each regime is time-varying, similarly to a GARCH model. The results for the order flow are qualitatively similar and we do not present them here. The reader interested in this model can consult the companion paper Tsaknaki et al. (2024).

4 We remind that the score is the derivative of the log-likelihood with respect to the parameter(s).

5 The Fisher matrix is defined as ${\mathcal{I}}_{t|t-1}={\mathbb{E}}_{t|t-1}\left[{\mathrm{\nabla }}_t^T{\mathrm{\nabla }}_t\right]$ where ${\mathrm{\nabla }}_t$ is defined in equation (26(26) $\begin{array}{rl}& {\mathrm{\nabla }}_t=\frac{\mathrm{\partial }\log p_u\left(u_t\right)}{\mathrm{\partial }{\rho }_t}\end{array}$(26) ).

6 This specification does not guarantee that $|{\rho }_t|\le 1$, thus sometimes one uses a link function (e.g. an inverse logistic) which maps $[-1,1]$ in $\mathbb{R}$, see Blasques et al. (2014). In our empirical analysis, we observe that the filtered $|{\rho }_t|$ is larger than 1 in less than one per thousand observations, thus we simply set a threshold $|{\rho }_t|\le 1$.

Appendices

Appendix 1.

Price dynamics inside an order flow regime

A.1. Concavity of order flow regimes

In table A1 we report the number of regimes that satisfy the condition: (A1) $Z_R:=\left|\frac{\sum _{t_R\le t<s_R}{x}_t}{\sum _{t_R\le t<s_R}|x_t|}\right|>\mathrm{\Theta }$(A1) for various values of Θ.

Table B1. Number of regimes satisfying the condition in equation (A1(A1) $Z_R:=\left|\frac{\sum _{t_R\le t<s_R}{x}_t}{\sum _{t_R\le t<s_R}|x_t|}\right|>\mathrm{\Theta }$(A1) ) for $\mathrm{\Theta }=0,0.5$ and 0.9, for both TSLA and MSFT with $\mathrm{\Delta t}=1$ min and $\mathrm{\Delta t}=3$ min.

	TSLA		MSFT
Θ	1 min	3 min	1 min	3 min
0	911	546	1394	690
0.5	674	466	1195	651
0.9	344	321	827	550

A.2. Testing the square root impact law

Table A2 presents the estimation of the parameters of the regression in equation (39(39) ${ϵ}_R\mathrm{\Delta }{p}_R=A{\left({ϵ}_R\sum _{t_R\le t<s_R}{x}_t\right)}^{\gamma }+\mathrm{noise}.$(39) ) when we consider the entire data sets without removing any outliers. We observe that the exponent γ is a bit larger than 1/2 being typically close to 0.7.

Table B2. Estimated parameters and their standard errors (SE) of the regression of equation (39(39) ${ϵ}_R\mathrm{\Delta }{p}_R=A{\left({ϵ}_R\sum _{t_R\le t<s_R}{x}_t\right)}^{\gamma }+\mathrm{noise}.$(39) ).

	TSLA
	1 min				3 min
Θ	A	SE of A	γ	SE of γ	A	SE of A	γ	SE of γ
0	0.032	0.013	0.732	0.04	0.08	0.04	0.654	0.047
0.5	0.045	0.018	0.7	0.039	0.096	0.047	0.639	0.046
0.9	0.05	0.023	0.691	0.044	0.221	0.101	0.564	0.043
	MSFT
	1 min				3 min
Θ	A	SE of A	γ	SE of γ	A	SE of A	γ	SE of γ
0	0.015	0.006	0.737	0.037	0.005	0.003	0.835	0.053
0.5	0.015	0.006	0.736	0.038	0.004	0.003	0.845	0.055
0.9	0.022	0.01	0.706	0.042	0.006	0.004	0.827	0.061

Note: The whole dataset is considered.

The outlier removal is obtained by the standard interquartile approach. Namely, we compute the first and third quartile $Q1$ and $Q3$, respectively, of log-returns. Then the data points outside the range $[Q1-1.5\mathrm{IQR},Q3+1.5\mathrm{IQR}]$, where IQR = Q3−Q1,are considered as outliers.

Figure A1 illustrates the data and the fits for the two stocks and the two timescales. The red points are those which are identified as outliers.

Figure A1. Data points and best fitting curves according to the regression of equation (39(39) ${ϵ}_R\mathrm{\Delta }{p}_R=A{\left({ϵ}_R\sum _{t_R\le t<s_R}{x}_t\right)}^{\gamma }+\mathrm{noise}.$(39) ) with $\mathrm{\Theta }=0$ for Tesla (a) with $\mathrm{\Delta t}=1$ min, (b) with $\mathrm{\Delta t}=3$ min and for Microsoft (c) with $\mathrm{\Delta t}=1$ min, (d) with $\mathrm{\Delta t}=3$ min. The plotted quantities are in basis points. The red points are the outliers. The green line is the fitting curve of the entire data set, while the orange line is the fitting curve when the outliers are excluded.

Appendix 2.

Comparison of order flow and market impact predictions under different regime shift models

Figure A2 compares the correlation function ${\mathcal{I}}_{ϵ}^{\left(m\right)}\left(k\right)$ of order flow for the MBOC and the BOCPD model. Figure A3 compares the predictor of market impact ${\mathcal{I}}_{\mathrm{\Delta }p}^{\left(m\right)}\left(k\right)$ for the same models.

Figure A2. Online order flow prediction according to MBOC (solid lines) and BOCPD (dotted lines) models with $\mathrm{\Delta t}=3$ min (a) for Tesla and (b) for Microsoft.

Figure A3. Online market impact prediction according to MBOC (solid lines) and BOCPD (dotted lines) models with $\mathrm{\Delta t}=3$ min (a) for Tesla and (b) for Microsoft.

From both figures, it is evident that the MBOC model outperforms the BOCPD when m>1. This justifies why in the main text we present the results obtained with the MBOC model.