Censored EM algorithm for Weibull mixtures: application to arrival times of market orders

Electronic Order Book (EOB) trading is now common to most stock exchanges. A set of EOB data from the FTSE for the period June to September 2010 was described and partly analysed in Ref.[1]. In this data set timestamps of changes in the EOB were given in milliseconds and because of the huge volume of EOB trading by (ultra-)high frequency trading algorithms, taking place on a microsecond scale, differences of timestamps appeared to be zero. This is why the timestamp time series appeared to be “zero-inflated”. Because of this rounding error, all differences in timestamps, $\Delta t$, are mapped to our actual data set as follows

	$\displaystyle\Delta t\in(0.0,0.5)=\mathcal{I}_{1}$	$\displaystyle\longrightarrow$	$\displaystyle 0=c_{1}$
	$\displaystyle\Delta t\in[0.5,1.5)=\mathcal{I}_{2}$	$\displaystyle\longrightarrow$	$\displaystyle 1=c_{2}$
	$\displaystyle\Delta t\in[1.5,2.5)=\mathcal{I}_{3}$	$\displaystyle\longrightarrow$	$\displaystyle 2=c_{3}$
		$\displaystyle etc.$

where the integer numbers $c_{1},c_{2},...$ are the rounded time stamp differences to the precision of milliseconds.

In Ref.[1] it was demonstrated that sets of non-negative time differences between subsequent market orders (MO) with $\Delta t>10$ milliseconds describe a random variable $X$ that could be fitted to parametric distributions, such as a left-truncated Weibull distribution, which passed a rigorous set of goodness-of-fit tests, such as Kolmogorof-Simirnov, Anderson-Darling, Cramer von Misses and Kuiper tests¹¹1The analysis in Ref.[1] also included limit orders (LO) but this will not be the topic of our letter..

Here we analyse the entire set of observations of the random variable $X$, including the ”zero-inflated” ones, by starting from a Weibull distribution for the whole range

$\displaystyle f_{i}(x|\alpha_{i},\beta_{i})=\frac{\beta_{i}}{\alpha_{i}}\left(\frac{x}{\alpha_{i}}\right)^{\beta_{i}-1}e^{-\left(\frac{x}{\alpha_{i}}\right)^{\beta_{i}}}\,,$

(1)

where $\alpha_{i}>0$ is the scale parameter and $\beta_{i}>0$ the shape parameter. Note that for $\beta_{i}=1$ we recover the exponential distribution. We denote the parameter vector by $\theta_{i}=(\alpha_{i},\beta_{i})$ and use the subscript index $i$ to denote various Weibull or exponential distributions with different parameter vectors $\theta_{i}$.

Now, the attempt to analyse the entire data set including the “zero-inflated” part (which sometimes was even more than half of the observed data set) requires a different procedure. One promising approach is to apply interval censoring to the observed data for the intervals of small time differences $\Delta t$: we use the above intervals $\mathcal{I}_{1},\mathcal{I}_{2},\mathcal{I}_{3},...$ together with the information of how many observations belong to them, i.e. $N_{1},N_{2},N_{3},...$. In this notation we can write an observed sample of a positive random variable $X$ as

	$\displaystyle x_{1},$	$\displaystyle x_{2},\cdots$	$\displaystyle,x_{n},x_{n+1},x_{n+2},...,x_{N}$		(2)
		$\displaystyle=$	$\displaystyle x_{1},x_{2},...,x_{n},\underbrace{c_{1},c_{1},...,c_{1}}_{N_{1}\text{times}},...,\underbrace{c_{L},c_{L},...,c_{L}}_{N_{L}\text{times}}\,,$

where, after grouping, all observations with index bigger than $n$ have been censored (with $L$ censoring intervals).

Despite censoring, the random variable $X$ itself will be assumed to come from a mixed distribution with a mixture of $M$ components. As a consequence of the findings in Ref.[1], we expect that a significant component of this mixed distribution should come from a Weibull distribution. Similar results where found in Ref.[3]: the authors found that the differences between subsequent MO’s for 30 DJIA stocks from the NYSE in October 1999 can be well described by Weibull distributions. Thus we want to investigate here mixtures of the following kind

$\displaystyle f(x|\theta)=\sum_{i=1}^{M}\Pi_{i}f_{i}(x|\theta_{i})$

(3)

where the weights, $\Pi_{i}$, satisfy $\Pi_{i}\geq 0$ with $\sum_{i=1}^{M}\Pi_{i}=1$ and the probability density functions $f_{i}(\cdot|\theta_{i})$ are given by Eq.(1), where the parameters $\theta_{i}$ need to be estimated by maximum likelihood estimation or other methods²²2In the following the index ”i” stands for the mixture components and the index ”j” for the number of observation.. Defining $\Theta=(\theta_{1},...,\theta_{M})$, the log-likelihood function $\mathcal{L}$ for $L$ censoring intervals $\mathcal{I}_{1},...,\mathcal{I}_{L}$ is now given as in Ref.[5] by

$\displaystyle\mathcal{L}=f(x_{1}|\Theta)\cdots\ f(x_{n}|\Theta)\cdot\left(\int_{\mathcal{I}_{1}}~{}dyf(y|\Theta)\right)^{N_{1}}\cdots\left(\int_{\mathcal{I}_{L}}~{}dyf(y|\Theta)\right)^{N_{L}}$

(4)

which is usually a difficult expression to handle. However, given our data sample Eq.(2), the estimation problem would be a standard maximum likelihood estimation (MLE) problem if we knew by some indicator $z_{ij}$ (taking values 0 or 1) denoting which observation $x_{j}$ belongs to which probability density function $f_{i}(\cdot|\theta_{i})$ and likewise some indicator $\tilde{z}_{i\ell}$ (also taking values 0 or 1) which censored observation $c_{\ell}$ belongs to which distribution. In this case, we would just group the observations and thus factorize the likelihood $\mathcal{L}=\mathcal{L}_{1}\cdots\mathcal{L}_{M}$ and separately maximize each group likelihood $\mathcal{L}_{i}$ corresponding to one mixture component $f_{i}(\cdot|\theta_{i})$.

The expectation-maximization (EM) algorithm Ref.[6, 7, 8] will be used to iteratively generate estimates for the indicators $z_{ij}$, for uncensored observations, $x_{j}$ and $\tilde{z}_{i\ell}$ for censored observations, $c_{\ell}$, and solve the estimation problem (provided it converges). The problem for censored mixtures has been discussed in some detail in Ref.[9]. In our study we apply this analysis for mixtures of exponential and Weibull distributions and study market order (MO) inter-arrival times. In section 2 we give a very brief review of the censored EM algorithm as given by Ref.[9] and in section 3 the relevant equations for Weibull distributions are given. Section 4 summarizes our results and conclusions.

The EM algorithm is an iterative procedure where an expectation-step (E-step) tries to estimate the unobservable indicators $z_{ij}$ and $\tilde{z}_{i\ell}$, respectively, and then uses the result in the maximization step (M-step) to estimate the parameters of the mixtures by an MLE. After the M-step there is another E-step and so on. Here we follow closely Ref.[9], where a hint is given that this iteration converges for certain function families (including exponential and Weibull distribution functions) in a certain limited parameter range. In the following sections the integer index $k$ denotes the number of the iteration.

For simplicity we now consider a 2-component mixture consisting of an exponential distribution, denoted by $i=1$, and a general Weibull distribution, denoted by $i=2$, as given in Eq.(1). Note that here the $\alpha$ and $\beta$ have a subscript $i$ for the mixture and a superscript $(k)$ for the iteration step in the EM algorithm. The expression corresponding to Eq.(8) is given in Appendix A. Our algorithmic results extend the results as given by Ref.[10] for exponential mixtures. It is clear how to modify our computations for arbitrary mixtures, e.g. $(p+r)$ mixtures consisting of $p$ exponentials and $r$ Weibulls,with $p,r=0,1,2,...$

The censored EM-algorithm in combination with a bootstrapping argument applied to the Baysean Information Criterion (BIC) allows us to choose as a model for the MO arrival times a two-component mixture distribution consisting of “1 exponential + 1 Weibull”. Of course, this conclusion is not applicable in the exceptional case of extraordinary trading activity in a stock when critical information is disclosed to the market participants. The first component of this mixed distribution is an exponential distribution with a relative weight of approximately 20% and a rather short scale parameter, in the range of 10 to 20 milliseconds. This result is independent of the stock under consideration and almost constant during the trading day. The second component, with a relative weight of approximately 80%, is a Weibull distribution for which the scale parameter varies intra-day with trading activity and lies typically between 1000 and 25000 milliseconds, albeit with universal shape parameter $\beta\approx 0.57$.

This result can be understood with a view to a typical stock exchange computer architecture (e.g. Ref.[20]): market orders are captured via various order gateways which forward the electronic orders to the so called accumulator, where timestamps are given. Low latency order gateways typically cater for high-frequency traders (e.g. algorithmic hedge funds) whose buy- and sell-orders are generated by computers located in the stock exchange’s immediate neighbourhood to minimize the transmission time. Here an exponential waiting time between these signals (with a short time scale parameter) is to be expected resulting in a Poisson process for high-frequency buy- and sell-orders. On the other hand, VSAT gateways cater for stock brockers whose customers include both institutional and private clients and whose orders are usually transmitted via an internet connection. It has been shown that internet traffic sends information as TCP-“parcels” with arrival times that can be described by Weibull distributions with a shape parameter less than 1 (see Ref.[21], Ref.[22]). In Ref.[21] it was shown that a shape parameter $\beta=0.569$ provided an excellent fit to the observed data. This is close to the value that we find for the Weibull distribution. The natural conclusion would be that for the time period we have studied trading on the LSE approximately 10% of the trading was done by market participants whose orders were generated electronically by a computer in the immediate neighborhood, whereas the other 90% of market participants were using the internet to generate their buy- and sell-orders.

Finally we want to emphasize that by Ref.[23] any Weibull distribution function with scale parameter $\beta$ smaller than 1 can be expressed as a superposition of exponential waiting time distributions. Thus our favorite model “1 Weibull + 1 exponential” is equivalent to a suitable mixture of exponential waiting time distributions. Consequently, the proposed censored EM algorithm for finite Weibull mixtures maybe compared to the problem discussed firstly in Ref.[17] and later expanded in Ref.[24] of how to fit a waiting time distribution consisting of a finite number of exponential waiting time distributions to the observed tick-by-tick data. Although exponential mixtures with more than 5 components seem to describe our actual data sets very well, they are excluded by the BIC due to their too high dof-value.