Bayesian semi-parametric estimation of compound inhomogeneous Poisson processes for ultra-high frequency financial transaction data

Abstract

Marked point processes provide a flexible framework for studying ultra-high frequency financial data that records the time and price for each transaction. This paper estimates compound, inhomogeneous Poisson processes (CIPP) where trades arrive according to an inhomogeneous Poisson process, and returns are drawn from a distribution after arrival. The intensity functions have Gaussian process priors to model time-of-day (diurnal) effects, which have bursts of activity at the open and close and a midday lull, and lagged returns from the previous trade, which are idiosyncratic and depend on the stock. The nonparametric density estimator for returns imposes a shape constraint to obtain a unimodal distribution that is asymmetric and leptokurtic. A mixture model separates trades into time sensitive and time insensitive trades with different intensity functions. The latter group tends to have larger orders and higher returns. The proposed model is compared to B-splines and geometric Brownian motion.

KEYWORDS:

• Gaussian process priors
• high frequency trading
• marked processes
• nonparametric density estimation
• shape constraints

1. Introduction

Diffusion processes for financial assets have had an enormous impact on financial research and practice since Bachelier (1900) and Black and Scholes (1973). The price series to estimate these models are usually sampled at regular intervals, such as the daily closing prices. Empirically, the movement in the sampled data often exceeds what would be expected from pure diffusion processes, which motivated numerous modifications and extensions (Durbin and Koopman 2000; Barndorff-Nielsen and Shephard 2001, 2006; Andersen et al. 2002; Andersen et al. 2003; Griffin and Steel 2006). Jump diffusion processes, which superimposes a jump process on the diffusion process (Press 1967, 1968; Merton 1976), model the outliers in the price series from the diffusion processes.

Despite the success of diffusion processes and their extensions, Engle (2000), Engle and Russell (1998), Russell and Engle (2010), and Fodra and Pham (2013), among others, maintain that ultra-high frequency data that records each transaction instead of sampling at fixed time intervals allows opportunities to discover stylised facts about the price process. Cont (2001), Kiem and Madhavan (1996), and Wood et al. (1985) provide empirical generalisations based on ultra-high frequency data such as the nature of the return distribution, auto-correlation, and stochastic volatility (non-constant variance) and the impact of time sampling on these statistics. Easley and O'Hara (1992) and Easley et al. (1997) test micro-economic models about the price discovery processes, to single out only a few representative instances in a rich and diverse literature (Campbell et al. 1997; Hautsch 2012). Additionally, a relatively new development in finance is high frequency trading (HFT), where firms with powerful computer systems, often located adjacent to clearinghouses’ servers, transact a large number of orders in fractions of seconds. Clark (2010) summarises risks posed by HFT. O’Hara (2015) summarises strategies for HFT, such as leveraging speed of trading, exploiting deterministic patterns from simple algorithms of other traders, and manipulating prices by triggering bidding wars. The study of HFT requires ultra-high frequency data because daily aggregates will mask the activity at the micro second level.

Marked point processes (Engle and Russell 1998; Engle 2000) provide a flexible and general framework for ultra-high frequency data and assume that the duration times between trades are stochastic and returns at trades are drawn from an appropriate model. Marked point processes provide a high resolution picture of the price discovery process that is not readily available from diffusion processes. Theoretically, as the sampling intervals for diffusion processes decrease, the continuity of the sample paths should persist. Empirically, realised prices are constant between transactions and can only change at transactions, if then. One bridge from continuous diffusion processes to ultra-high frequency data is to assume that there exists a latent price process that evolves between trades (Elerian et al. 2001) and that the empirical price process censors the latent one.

A special class of marked point process are compound, inhomogeneous Poisson process (CIPP) where transactions arrive according to an inhomogeneous Poisson process, and returns at arrivals are a random sample from a common distribution. Poisson processes are reasonable models for transactions of liquid or highly traded assets: they assume a large population of investors and a small probability that any particular investor makes a trade during a small time interval. CIPP can be consistent with the efficient market hypothesis, which implies that all information is encapsulated in stock prices (Fama 1970; Firth 1975; Basu 1977; Roll 1984). The efficient market hypothesis predicts that prices change in response to new, unexpected information that alters investors’ evaluation of the firm’s economic value. Though it may not hold in detail for all situations, it seems to be broadly true (Ying et al. 2019). A consequence of efficient markets is that past returns are not predictive of future returns, and CIPP assumes that returns are a random sample from a common distribution. The efficient market hypothesis is less proscriptive about the timing of trades. In a perfect world, economically rational investors would transact trades the moment they receive new information that changes their valuation of a firm’s value. In reality, information is noisy, various frictions delay updating and trading, and behavioural biases (Hirshleifer 2015) derail the best intentions to be ‘rational’.

We propose semi-parametric Bayes models for the intensity function and return distribution for CIPP. The intensity of the Poisson process has flexible functions of the time-of-day and lagged return from the last trade. These functions have Gaussian process priors. The parametric model for the intensity function includes day-of-week effects and the effect of the lagged duration between the last, two observed trades, which is similar to autoregressive conditional duration (ACD) model (Engle and Russell 1998; Ghysels et al. 2004; Russell and Engle 2006).

The model uses the spectral representation or Fourier series for Gaussian processes for the prior distribution of the return’s density (Lenk 1993, 2003) and arrival rates (Lenk 1999). The spectral representation projects the Gaussian process onto orthonormal basis functions to obtain a linear model for the Gaussian process. Smoothing priors (Whittaker 1923; Wahba 1978) for the spectral coefficients attenuate the effects of noisy, high-frequency basis functions by shrinking them towards zero. The smoothing prior determines both the smoothness of the functions and the tradeoff between the prior specification and the data. Our smoothing prior has a hierarchical Bayes specification, which allows the estimation of the smoothing parameter from the data. The nonparametric density estimator imposes a unimodal constraint by assuming that the density has an inverted U-shape (Lenk and Choi 2017), which is obtained by constraints on the derivatives of the nonparametric function (Ramsay 1998).

The Poisson intensity function has a bathtub shape for arrival rates during the trading day with a flurry of activity at the market open and close and a lull during mid-day. The relationship with the return from the previous trade seems idiosyncratic and depends on the security: overall, large price movements tend to trigger faster trading. Additionally, there are day-of-week effects and effects for the duration between the previous two trades. If the time between the last two trades is long, then the time until the next trade tends to be long too. The return distribution is slightly asymmetric and leptokurtic. We also compare the estimated model to geometric Brownian motion. Over longer periods of several days, the two are similar; however, they do deviate within a trading day.

We also propose a mixture model to illustrate segmenting the data into two components based on different functional forms of the Poisson intensity function. We call the two groups ‘time sensitive trades’ (TST) and ‘time insensitive trades’ (TIT) because the intensity function for TST includes diurnal effects and excludes information about the most recent trade, while TIT responds to the previous trade and excludes diurnal effects. Conceptually, TIT uses software and computer systems to trade without direct human intervention; it monitors the flow of trades and responds to it without regard to the time of day. Using their latency advantages, TIT are more responsive to the last trades than TST, which may be initiated by human decision making. Since TST is sensitive to diurnal factors and less responsive to the previous trade, which may have occurred in less than one second, we conceptually think of TST as being executed by humans instead of a trading algorithm. A logistic regression model for TIT indicates that they tend to have larger order sizes and returns than TST trades.

The paper has the following structure. Section 2 reviews compound, inhomogeneous Poisson processes for ultra-high frequency trading data. It also presents a simulation that highlights how time sampling a CIPP can lead to erroneous conclusions about the data-generating process. In particular, the time-sampled CIPP exhibits stochastic volatility (Engle 1982; Bollerslev 1986, 1987; Bollerslev et al. 1992) and leptokurtic returns, even though the CIPP simulation samples independent returns from a normal distribution with a constant variance. Section 3 develops the Bayesian nonparametric model for the Poisson intensity function and the return distribution. Section 4 presents the mixture model. Section 5 details the MCMC algorithm. Section 6 presents a simple model for off-hour trading. Off-hour trading could be ignored because it accounts for a very small proportion of all trades; however, it may be incumbent on researchers to include such details when investigating micro-structure with high-frequency data. Section 7 estimates the model with Trade and Quote data (TAQ) from Wharton Research Data Services and illustrates simulating from the model for option pricing. Section 8 concludes the paper.

2. Compound inhomogeneous Poisson process

The evolution of the stock price Y_t at time t can be written as a marked process that depends on the arrival of transactions and the rate of return at each transaction: (1) $Y_t=Y_s\prod _{i=1}^{N(s,t)}(1+r_i)\text{for}t>s$(1) where Y_s is the stock price at time s; N(s,t) is the number of transactions in the time interval s to t, and r_i is the return for the i^th transaction in s to t. Y_t is constant between trades and only jumps at trades. Equation (1)(1) $Y_t=Y_s\prod _{i=1}^{N(s,t)}(1+r_i)\text{for}t>s$(1) is a well-known identity from financial accounting and not a modelling assumption. Different assumptions about the number of transactions and return process result in different models.

We assume that the number of trades N(s,t) between s and t, is an inhomogeneous Poisson process (IPP) with intensity function λ(t) at time t. The Poisson assumption is appropriate for liquid markets with a large number of investors who act independently and have a small probability of making a transaction at any time conditional on the intensity function. Over time, the probability of investing changes with new information or investor whimsy, so unconditionally it can appear that trades are correlated because trades cluster. The distribution and expected value of N(s,t) are (Çinlar 1975): (2) $\begin{array}{rl}P[N(s,t)=k]& =\exp [-\mathrm{\Lambda }(s,t)]\frac{\mathrm{\Lambda }{(s,t)}^k}{k!}\text{for}k=0,1,\dots \\ \text{E}[N(s,t)]& =\mathrm{\Lambda }(s,t)=\underset{s}{\overset{t}{\int }}\lambda (x)\mathrm{dx}\text{for}s<t\end{array}$(2) Additionally, the number of trades over disjoint time intervals are mutually independent. The distribution of the time between trades or the duration time can be derived from Equation (2)(2) $\begin{array}{rl}P[N(s,t)=k]& =\exp [-\mathrm{\Lambda }(s,t)]\frac{\mathrm{\Lambda }{(s,t)}^k}{k!}\text{for}k=0,1,\dots \\ \text{E}[N(s,t)]& =\mathrm{\Lambda }(s,t)=\underset{s}{\overset{t}{\int }}\lambda (x)\mathrm{dx}\text{for}s<t\end{array}$(2) . If a trade occurs at time s, then the time of the next trade is T. If T is larger than t, then there has not been any arrivals between s and t, and N(s,t) is zero. The probability that the time of the next trade exceeds t given that the last trade occured at s is: $P(T>t|\text{Last trade at}s)=P[N(s,t)=0]=\text{exp}[-\mathrm{\Lambda }(s,t)]\text{for}t>s$Thus, T has a generalised exponential distribution with cumulative distribution F_T and density f_T: (3) $\begin{array}{c}{F}_T(t|s)=1-\text{exp}[-\mathrm{\Lambda }(s,t)]\\ f_T(t|s)=\lambda (t)\text{exp}[-\mathrm{\Lambda }(s,t)]\end{array}\text{for}t>s$(3) The times between trades are mutually independent, generalised exponential distributions. The duration between trades is T-s. The log-likelihood function for the instantaneous rate function for n trades at time t₀ < t₁ < … < t_n where t₀ is the origin is: (4) $\mathcal{L}(\lambda )=\left\{\sum _{i=1}^n\text{log}[\lambda (t_i)]\right\}-\mathrm{\Lambda }(t_0,t_n)$(4)

With our data the time stamp for trades are recorded in seconds, and we assume that they are rounded down.¹ The likelihood in Equation (4)(4) $\mathcal{L}(\lambda )=\left\{\sum _{i=1}^n\text{log}[\lambda (t_i)]\right\}-\mathrm{\Lambda }(t_0,t_n)$(4) is easily modified for discrete time recording because numerical integration is used to compute Λ. The discrete version of Equation (3)(3) $\begin{array}{c}{F}_T(t|s)=1-\text{exp}[-\mathrm{\Lambda }(s,t)]\\ f_T(t|s)=\lambda (t)\text{exp}[-\mathrm{\Lambda }(s,t)]\end{array}\text{for}t>s$(3) is: (5) $\begin{array}{rl}P(T\le s+1|\text{Last trade at }s)& =1-\text{exp}[-\mathrm{\Lambda }(s,s+1)]\\ P(s+k<T\le s+k+1|\text{Last trade at }s)& =\text{exp}[-\mathrm{\Lambda }(s,s+k)]-\text{exp}[-\mathrm{\Lambda }(s,s+k+1)]\\ & \text{for}k=1,2,\dots \end{array}$(5) CIPP requires a distribution for the returns. The efficient market hypothesis implies that the distribution of stock returns does not depend on past information, so we assume that the returns r_i for each transaction are random draws from the density f_R.

Time sampling from the CIPP can lead to return processes that have non-constant variance or stochastic volatility even though returns for each trade are independent. The total return from time s to t for s < t from Equation (1)(1) $Y_t=Y_s\prod _{i=1}^{N(s,t)}(1+r_i)\text{for}t>s$(1) is: $R(s,t)=\frac{Y_t-Y_s}{Y_s}=\prod _{i=1}^{N(s,t)}(1+r_i)-1$If the returns r_i are mutually independent with mean µ and standard deviation σ, then the mean and variance of the total returns are: (6) $\begin{array}{rl}E[R(s,t)]& =\text{exp}[\mathrm{\Lambda }(s,t)\mu ]-1\\ V[R(s,t)]& =\text{exp}[\mathrm{\Lambda }(s,t)({\sigma }^2+{\mu }^2+2\mu )]-\text{exp}[2\mathrm{\Lambda }(s,t)\mu ]\end{array}$(6)

When time sampling prices in equal length intervals, the variance of the sampled process in Equation (6)(6) $\begin{array}{rl}E[R(s,t)]& =\text{exp}[\mathrm{\Lambda }(s,t)\mu ]-1\\ V[R(s,t)]& =\text{exp}[\mathrm{\Lambda }(s,t)({\sigma }^2+{\mu }^2+2\mu )]-\text{exp}[2\mathrm{\Lambda }(s,t)\mu ]\end{array}$(6) is not constant if the integral of the intensity function over each time period Λ(s,t) is not constant.

To illustrate non-constant variance for the time sampled process, we simulate data from a CIPP where the returns are independent random draws from a normal distribution with mean 0 and standard deviation 0.04, and the intensity function that has a 20 d period.² Sampling the price process at the daily close results in leptokurtic returns and stochastic volatility even though returns are normally distributed and have a constant variance. Figure 1 graphs some of the main features of the CIPP model. Panel A graphs a realisation of the CIPP for 100,000 simulated trades over 500 trading days. Panel B graphs the return process from sampling the daily close and exhibits stochastic volatility. A GARCH(0,3) model provides a good fit and has highly significant coefficients. Panel C is the histogram of daily returns with kurtosis of 2.78 even though the true kurtosis is zero. The autocorrelation function for the squared daily returns, which is a measure of daily volatility, in Panel D provides further evidence for stochastic volatility. The simulation illustrates that time sampling from a CIPP introduces statistical features that are not present in the actual data-generating process. The key reason is that the number daily trades are random, and the Poisson rate parameter is not constant from day to day.

Figure 1. Simulated CIPP and daily close. (A) Simulated price process; (B) Returns based on daily close; (C) Histogram of returns from daily close; (D) ACF of squared returns from daily close.

3. Bayesian semiparametric model for CIPP

We use the following conventions that t is the consecutive number of seconds since the beginning of the observational period and that t_d is the time of day since midnight of the current day in seconds. For the first day of trading, t and t_d are equal. If t is in trading day n, then t = nk + t_d where k is 86,400 s per day. Our model considers three time periods: pre-open (4.00am to 9.30am), regular trading (9.30am to 4.00pm), and post-close (4.00pm to 8.00pm). The next section presents the model for pre-open and post-close trading. The processes restart at the opening of the three periods.

During regular trading hours, the CIPP’s intensity function λ is a function of the time of day t_d: (7) $\lambda (t_d|x,r)=\text{exp}[x^{\mathrm{\prime }}\beta +A(r)+B(t_d)]$(7) where x is a vector of covariates with regression coefficients β; A(r) is a flexible function of the lagged return from the previous transaction; and B(t_d) is a flexible function of the time of day t_d. In the application, the covariates are indicator variables for the day of week and the lagged duration between the last two transactions. The prior distribution for β is a normal distribution with mean m_β0 and variance V_β0²: β ∼ N(m_β0, V_β0²).

Since A and B are unknown functions, we assume that they have Gaussian process priors. We use a Fourier series representation or spectral analysis (Lenk 1999) to linearise the Gaussian processes. If Z is a second-order Gaussian process on the set ℬ, the mean process is E[Z(s)] = µ(s), and covariance function is Cov[Z(s),Z(t)] = ν(s,t). Without loss of generality, we set the mean process to zero since the parametric model could include the mean. In our case, ℬ is either regular trading hours from 9.30am to 4.00pm for the NYSE or an interval for the returns. The Karhunen-Loève representation projects Z into a linear space: $Z(s)=\sum _{k=1}^{\mathrm{\infty }}{\theta }_k{\phi }_k(s)$where φ_k are orthonormal basis functions. The spectral coefficients θ_k are defined as: ${\theta }_k=\underset{\mathcal{B}}{\overset{}{\int }}Z(s){\overline{\phi }}_k(s)\mathrm{ds}$where ${\overline{\phi }}_k$ is the complex conjugate of φ_k. The distribution of the spectral coefficients is inherited from the prior distribution of Z. Because the basis functions are orthonormal and Z is a Gaussian process, the spectral coefficients are mutually independent and have normal distributions. Their means and variances are given by: $\begin{array}{rl}E({\theta }_k)& =\underset{\mathcal{B}}{\overset{}{\int }}\mu (s){\overline{\phi }}_k(s)\mathrm{ds}\equiv {\mu }_k\\ \text{Var}({\theta }_k)& =\underset{\mathcal{B}}{\overset{}{\int }}\underset{\mathcal{B}}{\overset{}{\int }}\upsilon (s,t){\overline{\phi }}_k(s){\phi }_k(t)\mathrm{ds}\mathrm{dt}\equiv {\nu }_k^2\end{array}$Conversely, given the means and variances of the spectral coefficients, the mean and covariance of the Gaussian process is: $\begin{array}{rl}\mu (s)& =\sum _{k=1}^{\mathrm{\infty }}{\mu }_k{\phi }_k(s)\\ \upsilon (s,t)& =\sum _{k=1}^{\mathrm{\infty }}{\nu }_k^2{\phi }_k(s){\overline{\phi }}_k(t)\end{array}$provided that the variances ν_k² have a finite sum.

In the empirical section, the basis is the cosine basis: ${\phi }_k(s)=[2g_0(s){]}^{\frac{1}{2}}\text{cos}[\mathrm{k\pi }{G}_0(s)]\text{for}s\in \mathcal{B}\text{and}k>1$where g₀ and G₀ are the density and cumulative distribution functions for a base measure on ℬ. The cosine functions are an orthonormal basis on ℬ and are complete for piecewise continuous functions (Kreider et al. 1966). The cumulative distribution transforms G₀ are the Cauchy distribution for A(r) and the uniform distribution for B(t_d). The Fourier series do not include sine functions because the functions would be periodic on ℬ, and they do not include the constant function because it is part of the parametric model. Other choices of basis can be used.

An important issue with nonparametric estimation is overfitting the data, which leads to functions with spurious oscillations. The spectral representation provides a natural method to specify smoothness because the high frequency spectral coefficients will be small for smooth functions. If the spectral coefficients decay at rate o(k^m) for some m > 1, then the function is m times differentiable almost everywhere (Katznelson 2004). An exponential hierarchical Bayes (HB) smoothing prior shrinks the spectral coefficients to 0, and the amount of shrinkage increases with frequency: (8) ${\theta }_k\sim N[0,{\tau }^2\text{exp}(-\mathrm{k\gamma })]$(8) The variance parameter τ controls the tradeoff between the data and the smoothing prior. If τ is small (large), then the model puts more weight on the prior (data), and the function will be smoother (rougher). The variance parameter γ allows the amount of shrinkage towards zero to depend on frequency k: high-frequency spectral coefficients have more shrinkage to zero than low-frequency coefficients, thus dampening the high-frequency oscillations in the function that are due to noise and not signal.

Silverman (1985) observed that most nonparametric methods can be explicitly or implicitly viewed as kernel regression method. The covariance function for the Gaussian process priors plays the role of the kernel. Panel A of Figure 2 plots the spectral covariance function when τ = 2 and γ = 1 on the unit square for the exponential HB smoother. The covariance has a saddle shape and is maximum at the end points, which encodes the information that estimates of the regression functions at the endpoints tend to be more uncertain than those in the middle. The prior covariance is negative at the corners (0, 1) and (1, 0). One issue with flexible function estimation is end effects. This covariance attenuates end effects by increasing the prior variance at the end points of ℬ, which gives the observations near the end points more weight when estimating the function at the endpoints. The prior variance function is smallest at the middle, so observations at the middle have less weight in estimation. For comparison, Panel B of Figure 2 plots the implied covariance function³ for B-splines (Eilers and Marx 1996) where the prior standard deviation τ of the coefficients are 2, and there are equally spaced internal knots from 0 to 1 for a total of 13 knots from the function bs() of the R package splines. The covariance of the B-spline is concentrated along a diagonal band and is zero otherwise; only observations that are close to the abscissa are used in estimating the ordinate.

Figure 2. Covariance functions. (A) Spectral covariance function τ = 2 and γ = 1; (B) B-spline covariance function τ = 2 and 13 knots.

Instead of assuming fixed values for the smoothing parameters τ and γ or using cross-validation, we take a fully Bayesian approach with a hierarchical Bayes prior: (9) $\begin{array}{rl}{\tau }^2& \sim \mathrm{IG}\left(\frac{r_0}{2},\frac{s_0}{2}\right),\mathrm{the}\mathrm{inverse}\mathrm{Gamma}\mathrm{distribution}\\ \gamma & \sim G(1,w_{\gamma ,0}),\mathrm{the}\mathrm{Gamma}\mathrm{distribution}\text{.}\end{array}$(9) In Equation (9)(9) $\begin{array}{rl}{\tau }^2& \sim \mathrm{IG}\left(\frac{r_0}{2},\frac{s_0}{2}\right),\mathrm{the}\mathrm{inverse}\mathrm{Gamma}\mathrm{distribution}\\ \gamma & \sim G(1,w_{\gamma ,0}),\mathrm{the}\mathrm{Gamma}\mathrm{distribution}\text{.}\end{array}$(9) , the inverse Gamma distribution has shape parameter r₀/2, scale parameter s₀/2, and mean s₀/(r₀–2); and the Gamma distribution has shape parameter 1 (Exponential distribution), scale parameter w₀, and mean 1/w₀. The infinite series has to be truncated to implement it. The modelling strategy is to use a large number of cosine functions and let the smoothing prior shrink the insignificant ones towards zero instead of trying to pick the optimal number of cosine terms. Since smooth functions do not have high frequency oscillations, the estimates are insensitive to the number of cosines as long as the cutoff is sufficiently large.

(Lenk 1993, 2003) proposed a logistic-Gaussian prior for the nonparametric estimation of the density f_R: $f_R(r)=\frac{\text{exp}[Z(r)]}{{\int }_a^b\text{exp}[Z(u)]\mathrm{du}}\text{ for }a\le r\le b$where Z is a Gaussian process. However, this model is too flexible and will result in multiple modes due to the lumpiness in returns of ultra-high frequency data: empirically returns cluster because prices are discrete to the nearest $0.01. Unimodal shape constraints are imposed by assuming that the derivative of f_R before zero is positive and the derivative after zero is negative by introducing an auxiliary function h whose derivative changes signs at zero (Lenk and Choi 2017): (10) $\begin{array}{rl}{f}_R(r)& =\frac{g_0(r)\text{exp}\left\{{\int }_a^{r}h(s){[Z(s)]}^2\mathrm{ds}\right\}}{{\int }_a^b{g}_0(w)\text{exp}\left\{{\int }_a^{w}h(s){[Z(s)]}^2\mathrm{ds}\right\}\mathrm{dw}}\\ Z(s)& =\sum _{k=0}^K\sqrt{2}\text{cos}[\mathrm{k\pi }{G}_o(s)]{\theta }_{f,k}\text{for}a\le s,r\le b\\ h(s)& =\frac{1-\text{exp}(100s)}{1+\text{exp}(100s)}\end{array}$(10) where a and b are sufficiently large to include the domain of observed returns, and the base measure G₀ is a Cauchy distribution that has a mode at zero. The auxiliary logistic function h is between 1 and −1 and changes signs at 0. Its ‘slope’ is fixed to 100 so that it is essentially 1 before −0.1 and −1 after 0.1. The Gaussian process Z has a finite spectral representation that includes the constant function; however, the sign of the constant term θ_f,0 is not identified. We restrict it to be positive by assuming that it has a truncated normal distribution on the positive numbers. If the spectral coefficients θ_f,k are 0, then the density is the base density g₀. Z is squared in the exponent to make it a positive function. The first derivative of f_R: ${f_R}^{\mathrm{\prime }}(r)=\left\{h(r){[Z(r)]}^2+\frac{g_0^{\mathrm{\prime }}(r)}{g_0(r)}\right\}{f}_R(r)$is positive before zero and negative after zero because h and the derivative of g₀ are positive before 0 and negative after zero, and the other quantities are positive. Thus, the return density f_R is a unimodal with mode at zero. The spectral coefficients θ_f,k for k > 0 have the exponential HB smoothing priors given by Equations (8)(8) ${\theta }_k\sim N[0,{\tau }^2\text{exp}(-\mathrm{k\gamma })]$(8) and (9)(9) $\begin{array}{rl}{\tau }^2& \sim \mathrm{IG}\left(\frac{r_0}{2},\frac{s_0}{2}\right),\mathrm{the}\mathrm{inverse}\mathrm{Gamma}\mathrm{distribution}\\ \gamma & \sim G(1,w_{\gamma ,0}),\mathrm{the}\mathrm{Gamma}\mathrm{distribution}\text{.}\end{array}$(9) .

4. Mixture model

This section proposes a mixture model for the intensity functions where trades can be classified into two or more groups. These groups can be meaningful for economic theory and the simplest dichotomy is ‘rational’ versus ‘noise’ trades (Ramiah et al. 2015). Stocks are a contingent claim on future earnings, so the economic value of a stock is the expected present value of free cash flow. Rational trades occur when the market price, which is determined by supply and demand conditions, deviates from their economic value. Rational traders may disagree about economic value because they can have different opinions of future cash flow or use different discount rates, thus leading to different valuations; however, in the long run market prices should adhere to economic value for liquid assets in efficient markets with rational investors, and trading data reveals the price discovery process. In contrast, noise trades are triggered by various processes, including behavioural, psychological, and situational factors. These trades introduce noise into the price discovery process. Market anomalies, i.e. material deviations between economic value and stock price, may occur if enough investors are subject to the same biases, resulting in potential market failures.

The distinctions between rational and noise trades are easier to make in theory than empirically. The ideal situation would be an algorithm that classifies trades as rational or noise in real time. An economic approach would require a valuation model to determine if a trade was consistent with economic theory; however, the algorithm’s valuation may be one of many possible realisations of the future economy. At the other extreme, a purely data-driven classification model would attempt to find clusters in the data. When applied to trading times, it tends to find outliers from the inhomogeneous Poisson process model.

This section proposes a mixture model between the two extremes of using an economic evaluation model and cluster analysis by including different variables in the Poisson intensity function for the different segments of trades. We do not claim to solve the rational versus noise question and refrain from using those labels. The motivation for the model is that economic valuation models typically do not include diurnal factors. The first group, called ‘time insensitive trades’ or TIT, do not systematically respond to diurnal factors. We vaguely conceptualise this group as trades that are executed by algorithms, which have the ability to respond to recent trades and do not become sleepy after lunch. The second group, called ‘time sensitive trades’ or TST, are subject to diurnal effects but lack the speed to respond to the most recent trades. We conceptualise TST as being executed by human traders. The intensity function for TST includes time of day and day of week and excludes lagged returns and lagged durations.

The mixture model modifies the intensity function in Equation (7)(7) $\lambda (t_d|x,r)=\text{exp}[x^{\mathrm{\prime }}\beta +A(r)+B(t_d)]$(7) for the time sensitive trades (TST) and time insensitive trades (TIT). The Poisson intensity functions for TST (Z = 0) and TIT (Z = 1) are: (11) $\begin{array}{rl}{\lambda }_0(t_d)& =\exp [{\beta }_1\mathrm{Mo}+{\beta }_2\mathrm{Tu}+{\beta }_3\mathrm{Wd}+{\beta }_4\mathrm{Th}+{\beta }_5\mathrm{Fr}+B(t_d)]\text{for TST}\\ {\lambda }_1(t_d)& =\exp [{\beta }_0+{\beta }_{6}D+A(r)]\text{ for TIT}\end{array}$(11) where Mo, Tu, Wd, Th, and Fr are day-of-week indicator variables; D is the duration between the previous two trades before time t; A(r) is the flexible function of the previous return, and B(t_d) is the flexible function for time of day. Additionally, we use a logistic regression model for time insensitive trades: (12) $\text{logit}(P[\text{TIT}])={\alpha }_0+{\alpha }_1\text{Return}+{\alpha }_2\text{log}(\text{Order Size})$(12) where Return is the return between the current trade and the next trade, and Order Size is the number of shares being transacted. Equation (12)(12) $\text{logit}(P[\text{TIT}])={\alpha }_0+{\alpha }_1\text{Return}+{\alpha }_2\text{log}(\text{Order Size})$(12) is used for classification purposes and not intended as a data generating mechanism since the return for the current trade is not known at the time of transaction.

5. MCMC algorithm

The MCMC algorithm uses random walk, Metropolis Hastings (MH) where the candidate values are drawn from normal distributions that are centred at the current draws of the MCMC algorithm. The variances of the generating distributions are adaptively determined from the past draws of the MCMC algorithm (Haario et al. 2001; Roberts and Rosenthal 2009). For a generic parameter δ of dimension d, the adaptive algorithm generates a candidate draw from a normal distribution at iteration t + 1: (13) ${\delta }_{t+1}^c={\delta }_t+c{\tilde{V}}_t^{1/2}{z}_t$(13) where δ_t is the current draw from iteration t, ${\tilde{V}}_t$ is an estimate of the variance of δ, c is a scale factor, and z_t are independent draws from a standard normal distribution. The suggested scale factor c is 2.38/d^1/2. The candidate draw is either accepted according to the standard MH probability or the current value δ_t is retained at iteration t + 1.

One method to estimate the variance is with exponentially weighted moving averages (EWMA): (14) $\begin{array}{rl}{\tilde{\delta }}_t& =w_t{\delta }_t+(1-w_t){\tilde{\delta }}_{t-1}\\ {\tilde{V}}_t& =w_t[({\delta }_t-{\tilde{\delta }}_t)({\delta }_t-{\tilde{\delta }}_t{)}^{\mathrm{\prime }}]+(1-w_t){\tilde{V}}_{t-1}\end{array}$(14) where w_t is a positive constant less than 1. The EWMA parameter w_t can depend on the iteration: for instance, w_t = w when t ≤ t₀; w_t = w/(t–t₀) when t > t₀; w = 0.2; and t₀ is an initial burn-in period. This specification reduces the adaptation after the chain reaches its stationarity distribution. The EWMAs in Equation (14)(14) $\begin{array}{rl}{\tilde{\delta }}_t& =w_t{\delta }_t+(1-w_t){\tilde{\delta }}_{t-1}\\ {\tilde{V}}_t& =w_t[({\delta }_t-{\tilde{\delta }}_t)({\delta }_t-{\tilde{\delta }}_t{)}^{\mathrm{\prime }}]+(1-w_t){\tilde{V}}_{t-1}\end{array}$(14) are updated when a candidate draw is accepted; otherwise, a long series of rejecting the candidates can result in variance estimates of zero. Roberts and Rosenthal (2009) propose a mixture distribution where 95% of the mass is the generating distribution in Equation (13)(13) ${\delta }_{t+1}^c={\delta }_t+c{\tilde{V}}_t^{1/2}{z}_t$(13) and 5% of the mass substitutes the identity matrix for the estimated variance to keep the estimated variances from becoming zero. We take a different approach where we replace eigenvalues of the estimated variance that are below a small, positive number with that small positive number, and only generate the candidate from Equation (13)(13) ${\delta }_{t+1}^c={\delta }_t+c{\tilde{V}}_t^{1/2}{z}_t$(13) .

We make several modifications to the standard algorithm. First, we use the current values of the HB standard deviations τexp(-kγ/2) when generating the spectral coefficients θ_k for the Gaussian process priors. Second, we found that the suggest scale factor in Equation (13)(13) ${\delta }_{t+1}^c={\delta }_t+c{\tilde{V}}_t^{1/2}{z}_t$(13) did not work with our data and model. Instead, we generate a scale parameter ξ_t from an inverse Gamma distribution where the shape parameter is fixed to a small integer (α = 3), and the scale parameter is estimated from an EWMA of the past ξ: (15) $\begin{array}{rl}{\tilde{\xi }}_t& =w_t({\xi }_t+{\xi }_{\min })+(1-w_t){\tilde{\xi }}_{t-1}\\ {\xi }_{t+1}& \sim \mathrm{IG}(\alpha ,(\alpha -1){\tilde{\xi }}_t)\end{array}$(15) where ξ_min is a small non-negative constant which keeps the EWMA positive, and the weight w_t is positive and less than one. The mean of the inverse Gamma distribution is ${\tilde{\xi }}_t$, and the standard deviation is ${\tilde{\xi }}_t$/√(α−2). The inverse Gamma distribution with a small shape parameter is right skew, so the step size will occasionally be relatively large, which can help with mixing. The EWMA ${\tilde{\xi }}_t$ is only updated when the candidate draw is accepted. We also monitor the acceptance rates with an EWMA where the datum is one if the candidate draw is accepted and zero if the current draw retained on each MCMC iteration. If the acceptance rates are too small (< 0.1), the step size of the generating distribution is too large, and ${\tilde{\xi }}_t$ in Equation (15)(15) $\begin{array}{rl}{\tilde{\xi }}_t& =w_t({\xi }_t+{\xi }_{\min })+(1-w_t){\tilde{\xi }}_{t-1}\\ {\xi }_{t+1}& \sim \mathrm{IG}(\alpha ,(\alpha -1){\tilde{\xi }}_t)\end{array}$(15) is reduced by a multiplicative factor that is randomly drawn between 0.1 and 0.5. If the acceptance rates to too big (> 0.5), the step size is too small, and ${\tilde{\xi }}_t$ is increased by a multiplicative factor that is randomly generated between 2 and 10.

The duration model has three sets of parameters: the regression parameters β, the spectral coefficients θ_A,k for A(r), and the spectral coefficients θ_B,k for B(t_d). Their generating distributions at iteration t + 1 are: (16) $\begin{array}{rl}{\beta }_{t+1}^c& ={\beta }_t+{\xi }_{\beta ,t}^{1/2}{\tilde{V}}_{t,\beta }^{1/2}{z}_{t,\beta }\\ {\theta }_{A,k,t+1}^c& ={\theta }_{A,k,t}+{\xi }_{A,t}^{1/2}{\tau }_{A,t}\text{exp}\left(-\frac{k{\gamma }_{A,t}}{2}\right)z_{t,A,k}\text{ for}k=1,\dots ,K_A\\ {\theta }_{B,k,t+1}^c& ={\theta }_{B,k,t}+{\xi }_{B,t}^{1/2}{\tau }_{B,t}\text{exp}\left(-\frac{k{\gamma }_{B,t}}{2}\right)z_{t,B,k}\text{ for}k=1,\dots ,K_B.\end{array}$(16) There are three separate scale parameters ξ_t for the three blocks of parameters in Equation (16)(16) $\begin{array}{rl}{\beta }_{t+1}^c& ={\beta }_t+{\xi }_{\beta ,t}^{1/2}{\tilde{V}}_{t,\beta }^{1/2}{z}_{t,\beta }\\ {\theta }_{A,k,t+1}^c& ={\theta }_{A,k,t}+{\xi }_{A,t}^{1/2}{\tau }_{A,t}\text{exp}\left(-\frac{k{\gamma }_{A,t}}{2}\right)z_{t,A,k}\text{ for}k=1,\dots ,K_A\\ {\theta }_{B,k,t+1}^c& ={\theta }_{B,k,t}+{\xi }_{B,t}^{1/2}{\tau }_{B,t}\text{exp}\left(-\frac{k{\gamma }_{B,t}}{2}\right)z_{t,B,k}\text{ for}k=1,\dots ,K_B.\end{array}$(16) . Each scale factor follows the updating and generating distributions in Equation (15)(15) $\begin{array}{rl}{\tilde{\xi }}_t& =w_t({\xi }_t+{\xi }_{\min })+(1-w_t){\tilde{\xi }}_{t-1}\\ {\xi }_{t+1}& \sim \mathrm{IG}(\alpha ,(\alpha -1){\tilde{\xi }}_t)\end{array}$(15) . Another modification sets the variance of β to a diagonal matrix with a separate ξ for each parameter. We tried both methods, and believe the latter is more robust than the former. When the chain is near the posterior mode, the estimated variance of β is accurate, and the MH sampler works at least as well as the diagonal variance. However, when the chain is in a low probability region, the estimated variance can cause the chain to become stuck.

The log-likelihood uses the discretised form of the nonhomogeneous distribution in Equation (5)(5) $\begin{array}{rl}P(T\le s+1|\text{Last trade at }s)& =1-\text{exp}[-\mathrm{\Lambda }(s,s+1)]\\ P(s+k<T\le s+k+1|\text{Last trade at }s)& =\text{exp}[-\mathrm{\Lambda }(s,s+k)]-\text{exp}[-\mathrm{\Lambda }(s,s+k+1)]\\ & \text{for}k=1,2,\dots \end{array}$(5) since the trades are rounded down to the nearest second. The integral of the rate parameters for regular trading hours are computed with numerical methods where B(t_d) is computed on a grid with gaps of one second. We experimented with MH by generating β, θ_A,k, and θ_B,k simultaneously and in three blocks. The advantage of three blocks is that the scale parameters ξ can be turned separately. For a given number of iterations, the single block method is faster since the log-likelihood is computed once per iteration, while it is computed three times per iteration for the blocked MH. However, we did not consistently find that one method dominates the other.

The nonparametric model for the return distribution in Equation (10)(10) $\begin{array}{rl}{f}_R(r)& =\frac{g_0(r)\text{exp}\left\{{\int }_a^{r}h(s){[Z(s)]}^2\mathrm{ds}\right\}}{{\int }_a^b{g}_0(w)\text{exp}\left\{{\int }_a^{w}h(s){[Z(s)]}^2\mathrm{ds}\right\}\mathrm{dw}}\\ Z(s)& =\sum _{k=0}^K\sqrt{2}\text{cos}[\mathrm{k\pi }{G}_o(s)]{\theta }_{f,k}\text{for}a\le s,r\le b\\ h(s)& =\frac{1-\text{exp}(100s)}{1+\text{exp}(100s)}\end{array}$(10) also has a Gaussian process with spectral coefficients θ_R,k. The candidate distribution for MCMC iteration t + 1 is similar to those in Equation (16)(16) $\begin{array}{rl}{\beta }_{t+1}^c& ={\beta }_t+{\xi }_{\beta ,t}^{1/2}{\tilde{V}}_{t,\beta }^{1/2}{z}_{t,\beta }\\ {\theta }_{A,k,t+1}^c& ={\theta }_{A,k,t}+{\xi }_{A,t}^{1/2}{\tau }_{A,t}\text{exp}\left(-\frac{k{\gamma }_{A,t}}{2}\right)z_{t,A,k}\text{ for}k=1,\dots ,K_A\\ {\theta }_{B,k,t+1}^c& ={\theta }_{B,k,t}+{\xi }_{B,t}^{1/2}{\tau }_{B,t}\text{exp}\left(-\frac{k{\gamma }_{B,t}}{2}\right)z_{t,B,k}\text{ for}k=1,\dots ,K_B.\end{array}$(16) except the Gaussian process includes the constant (k = 0), which is assumed to be positive to identify the model: $\begin{array}{rl}{\theta }_{R,0,t+1}^c& ={\theta }_{R,0,t}+{\xi }_{R,0,t}^{1/2}{z}_{t,R,0}\text{where}{\theta }_{R,0,t+1}^c>0\\ {\theta }_{R,k,t+1}^c& ={\theta }_{R,k,t}+{\xi }_{R,t}^{1/2}{\tau }_{R,t}\text{exp}\left(-\frac{k{\gamma }_{R,t}}{2}\right)z_{t,R,k}\text{for}k=1,\dots ,K_R\end{array}$while ξ_R,0,t and ξ_R,t are generated as in Equation (15)(15) $\begin{array}{rl}{\tilde{\xi }}_t& =w_t({\xi }_t+{\xi }_{\min })+(1-w_t){\tilde{\xi }}_{t-1}\\ {\xi }_{t+1}& \sim \mathrm{IG}(\alpha ,(\alpha -1){\tilde{\xi }}_t)\end{array}$(15) . The constant is drawn from a truncated normal distribution by using the inverse cdf method. The Gaussian process and density function are computed on a fine grid that extend slightly beyond the smallest and largest observed returns.

The HB smoothing parameters (τ²,γ) for each Gaussian process are updated from their full conditional distributions given the spectral coefficients θ_k. In the following, we drop the subscripts for the different Gaussian processes. Their prior distributions are in Equation (9)(9) $\begin{array}{rl}{\tau }^2& \sim \mathrm{IG}\left(\frac{r_0}{2},\frac{s_0}{2}\right),\mathrm{the}\mathrm{inverse}\mathrm{Gamma}\mathrm{distribution}\\ \gamma & \sim G(1,w_{\gamma ,0}),\mathrm{the}\mathrm{Gamma}\mathrm{distribution}\text{.}\end{array}$(9) . The full conditional distribution of τ² is conjugate: $\begin{array}{rl}{\tau }^2& \sim \mathrm{IG}\left(\frac{r_K}{2},\frac{s_K}{2}\right)\\ r_K& =r_0+K\\ s_K& =s_0+\sum _{k=1}^K{\theta }_k^2\text{exp}(\mathrm{k\gamma })\end{array}$The full conditional distribution of γ is the product of K normal distributions and an Exponential distribution: (17) $f(\gamma )\propto \text{exp}\left[\left(\frac{K(K+1)}{2}-w_{\gamma ,0}\right)\gamma -\frac{1}{2{\tau }^2}\sum _{k=1}^K{\theta }_k^2\text{exp}(\mathrm{k\gamma })\right].$(17) We use slice sampling (Damien et al. 1999), which decomposes complex distributions with auxiliary, uniformly distributed random variables. The full conditional distribution in Equation (17)(17) $f(\gamma )\propto \text{exp}\left[\left(\frac{K(K+1)}{2}-w_{\gamma ,0}\right)\gamma -\frac{1}{2{\tau }^2}\sum _{k=1}^K{\theta }_k^2\text{exp}(\mathrm{k\gamma })\right].$(17) with the auxiliary random variables can be written as $\begin{array}{rl}f(\gamma ,V_1,\dots ,V_k)& \propto \prod _{k=1}^{K}I[0\le V_k\le \text{exp}\{-a_k\text{exp}(\mathrm{k\gamma })\}]\text{exp}(w_{\gamma ,K}\gamma )\\ a_k& =\frac{{\theta }_k^2}{2{\tau }^2}\\ w_{\gamma ,K}& =\frac{K(K+1)}{2}-w_{\gamma ,0}\end{array}$At iteration t + 1, the full conditional of distributions of V_k are independent, uniform distributions and are generate as: v_k,t + 1= u_kexp[-a_kexp(kγ_t)] where u_k ∼ U(0,1). Given v_k,t + 1, the full conditional density of γ is proportional to exp(w_Kγ) given the bounds: $u_k\text{exp}[-a_k\text{exp}(k{\gamma }_t)]\le \text{ exp}[-a_k\exp (\mathrm{k\gamma })]\text{for }k=1,\dots ,K$

Inverting the inequalities as a function of γ gives: (18) $\gamma \le b_k={\gamma }_t+\frac{1}{k}[\text{log}\{a_k-\log (u_k)\exp (-k{\gamma }_t)\}-\log (a_k)]\text{for}k=1,\dots ,K.$(18) Combining the bounds restricts γ to be less than b = min(b_k). The full conditional distribution of γ is proportional to exp(w_k γ) given γ ≤ b. The draw γ_t + 1 is generated with the inverse cdf method for bounded variables: (19) ${\gamma }_{t+1}=b+w_{\gamma ,K}^{-1}\text{log}[u+(1-u)\text{exp}(-b{w}_{\gamma ,K})]$(19) The forms of Equations (18)(18) $\gamma \le b_k={\gamma }_t+\frac{1}{k}[\text{log}\{a_k-\log (u_k)\exp (-k{\gamma }_t)\}-\log (a_k)]\text{for}k=1,\dots ,K.$(18) and (19)(19) ${\gamma }_{t+1}=b+w_{\gamma ,K}^{-1}\text{log}[u+(1-u)\text{exp}(-b{w}_{\gamma ,K})]$(19) avoid overflow because the index k and w_γ,K can be large positive numbers.

The MCMC algorithm for the mixture model that classifies trades into time sensitive trades (TST) and time insensitive trades (TIT) adapts standard MCMC procedures with the modification that the two groups have structurally different distribution functions in Equation (11)(11) $\begin{array}{rl}{\lambda }_0(t_d)& =\exp [{\beta }_1\mathrm{Mo}+{\beta }_2\mathrm{Tu}+{\beta }_3\mathrm{Wd}+{\beta }_4\mathrm{Th}+{\beta }_5\mathrm{Fr}+B(t_d)]\text{for TST}\\ {\lambda }_1(t_d)& =\exp [{\beta }_0+{\beta }_{6}D+A(r)]\text{ for TIT}\end{array}$(11) and the ‘prior’ classification probabilities depend on covariates. The MCMC algorithm introduces indicator variables z_i which are one if trade i is classified as TIT and 0 if trade i is classified as TST. The ‘prior’ probability that trade i belongs to TIT (z_i = 1) is: $P(z_i=1|u_i)=\frac{\text{exp}(u_i^{\mathrm{\prime }}\alpha )}{1+\text{exp}(u_i^{\mathrm{\prime }}\alpha )}$where u_i is a vector of covariates for trade i. Given the parameter draws at time t, the MCMC algorithm generates draws at time t + 1 according to the following steps.

1. Given the group indicators z_i of the trades, generate β₀ (parameters for days of week) and θ_B,k (spectral coefficients for the Gaussian process for time of day) using the observations in the TST group, and generate β₁ (parameters for constant and lagged durations) and θ_A,k (spectral coefficients for the Gaussian process for lagged returns) using the observations in the TIT group. The algorithm modifies the generation distribution in Equation (16)(16) $\begin{array}{rl}{\beta }_{t+1}^c& ={\beta }_t+{\xi }_{\beta ,t}^{1/2}{\tilde{V}}_{t,\beta }^{1/2}{z}_{t,\beta }\\ {\theta }_{A,k,t+1}^c& ={\theta }_{A,k,t}+{\xi }_{A,t}^{1/2}{\tau }_{A,t}\text{exp}\left(-\frac{k{\gamma }_{A,t}}{2}\right)z_{t,A,k}\text{ for}k=1,\dots ,K_A\\ {\theta }_{B,k,t+1}^c& ={\theta }_{B,k,t}+{\xi }_{B,t}^{1/2}{\tau }_{B,t}\text{exp}\left(-\frac{k{\gamma }_{B,t}}{2}\right)z_{t,B,k}\text{ for}k=1,\dots ,K_B.\end{array}$(16) for the homogenous model.

2. Given the group indictors, generate α by using standard, MH procedures for logistic regression.

3. Reclassify the trades according to the posterior probabilities of group memberships. The duration between trade i-1 at time t_i-1 and trade i at time t_i is y_i = t_i – t_i–1. The posterior classification probability for trade i is: $P_{t+1}(z_i=1|y_i,u_i)=\frac{P(z_i=1|u_i)f_1(y_i|t_{i-1})}{P(z_i=0|u_i)f_0(y_i|t_{i-1})+P(z_i=1|u_i)f_1(y_i|t_{i-1})}$ where f₀ and f₁ are the probabilities (Equation 5(5) $\begin{array}{rl}P(T\le s+1|\text{Last trade at }s)& =1-\text{exp}[-\mathrm{\Lambda }(s,s+1)]\\ P(s+k<T\le s+k+1|\text{Last trade at }s)& =\text{exp}[-\mathrm{\Lambda }(s,s+k)]-\text{exp}[-\mathrm{\Lambda }(s,s+k+1)]\\ & \text{for}k=1,2,\dots \end{array}$(5) ) for the durations under the TIT and TST duration models.

Because the models for TIT and TST differ in their specifications, the MCMC does not have label switching (Titterington et al. 1987).

In the empirical section, the MCMC had a total of 100,000 iterations and the last 10,000 iterations were used for estimation. Execution times averaged from two to three minutes per 10,000 iteration per 10,000 trades for 30 cosine functions using R with a Xeon E5-1660 CPU running at 3.20 GHz with 64 GB of RAM. Convergence diagnostics (Geweke and Heidel) from the CODA package in R indicate that the MCMC converged.⁴

6. Off-hours trading

For the NYSE pre-open hours are 4am to 9.30am, and post-close hours are 4pm to 8pm. Off-hour trading incurs additional transaction costs, which limits trading. There may be societal factors that limit off-hour trading: professional traders who are located in the United States are not particularly keen to extend their work-day from 4am to 8pm. Off-hours trading tends to be too sporadic to use nonparametric methods for inhomogeneous Poisson processes. If trading were 24 h a day, then it would be possible to extend B to the entire day by including both sine and cosine terms since B would be a periodic function. In our empirical example, only 97 of 117,276 trades over 80 days were off-hours for UGI, and only 0.15% of PG trades where off-hours. Additionally, most of the off-hour trades are clustered around the open or close, so the data are very thin at other times. Off-hour trades are not central to this paper’s focus, but we include models for pre-open and post-close trading for completeness. We assume that the distribution of returns for off-hour trades is the same as the distribution during regular hours. Because off-hour trading is light, the Poisson process intensity function for regular hours does not apply to off-hours.

The models for off-hour trading has multiple stages. Not all days have pre-open trading. A given day has pre-open trading with probability p, and days are mutually independent. If a day has pre-open trading, there usually is a lag between the 4am start of pre-open trading and the first trade, which usually occurs between 7am and 9am. This lag has a normal distribution.⁵ After the first trade, subsequent trades are from a homogeneous Poisson process with rate parameter λ₁ until the market opens. Post-close trading tends to cluster around the market close at 4.30pm with very few, if any, trades after 6pm, even though post-close trading does not end until 8pm. The number of trades in post-close trading has a Poisson distribution with parameter or expected number λ_2a. Given the total number of trades, we assume that the trades arrive according to a homogeneous Poisson process with rate parameter λ_2b. We use conjugate priors: p has a Beta distribution, the rate Poisson parameters have Gamma distributions, the mean of the lag parameter given the variance has a normal distribution, and the variance has an inverse Gamma distribution.

7. Data and results

TAQ (Trade and Quote) data from Wharton Research Data Service for four stocks: Citigroup (CITI), Procter and Gamble (PG), UGI Utilities (UGI), and Walmart (WMT) are used to estimate the models. The price process was observed from 12.00am on 3 January 2011 to 8.00pm on 23 April 2011. These four stocks were selected for different industries and different sizes. The first 60 trading days are used for estimation. The top of Table 1 gives total numbers of trades in the 60 days, the numbers for regular, pre-open, and post-close trading for the four stocks. Table 1 also presents the estimated posterior means for the parametric model for the intensity function for the four securities. The top third of the table is for the homogeneous model, while the bottom two-thirds of the table are for the mixture model. There are significant differences among the daily effects, and the coefficients for lagged durations are negative. If the duration of the last trade was long, then the current duration will also tend to be long. In the mixture model, time sensitive trades (TST) only have the day-of-week indicators, while time insensitive trades (TIT) have a constant and lagged duration. The probability of TIT increases with contemporaneous return and order size, except the coefficients of order size for WMT and return for UGI are negative. TIT or algorithmic-generated trades tend to have larger order sizes and reap larger returns to the detriment of TST or human-generated trades. The last row of Table 1 reports the proportion of TIT trades. CITI has the highest proportion of TIT, and WMT has the highest proportion of TST.

Table 1. Estimated parametric components for duration models. All of the coefficients are significant at 0.05 according to the Bayesian p-value^a except Tuesday and Wednesday for CITI in the mixture model. TST are time sensitive trades. TIT are time insensitive trades.

Number of Trades	CITI	PG	UGI	WMT
Total for Estimation	1,242,560	1,003,464	91,193	1,177,086
Regular Hours Trading	1,165,358	1,001,970	91,122	1,175,578
Pre-Open Trading	55,019	1157	4	916
Post-Close Trading	22,183	337	67	592
Homogeneous	CITI	PG	UGI	WMT
Monday	−0.2496	−0.5574	−1.2451	−0.4462
Tuesday	−0.1786	−0.5122	−1.3153	−0.3318
Wednesday	−0.1800	−0.5101	−1.3582	−0.3153
Thursday	−0.2451	−0.4599	−1.3544	−0.3365
Friday	−0.2167	−0.5159	−1.5052	−0.3752
Lag Duration	−0.1508	−0.0741	−0.0074	−0.0960
Mixture: TST	CITI	PG	UGI	WMT
Monday	0.0234	0.6310	−0.4511	1.4692
Tuesday	0.0011	0.8134	−0.4049	−0.2448
Wednesday	−0.0483	0.7832	0.4971	0.5749
Thursday	−0.2181	0.5139	0.5140	1.2445
Friday	−0.9907	1.2802	0.8205	−0.0675
Mixture: TIT	CITI	PG	UGI	WMT
Constant	0.9755	−1.5892	−1.6550	−1.8171
Lag Duration	−0.4544	−0.0365	−0.0052	−0.0437
Logit P(TIT)	CITI	PG	UGI	WMT
Intercept	−0.0701	−1.5936	0.9257	−1.7713
Return	0.3241	1.1007	0.7962	−2.3495
Log Order Size	0.1820	0.0133	−0.2140	0.0090
Proportion of TIT Trades	0.7331	0.1784	0.4724	0.1512

^a The coefficients β is significant at 0.05 if at least 95% of the posterior distribution is shifted away from zero: P(β >0|Data) > 0.025 or P(β <0|Data)<0.025.

Figure 3 illustrates some of the major features of ultra-high frequency data for PG. Panel A plots the stock price process for 13,854 trades on 1 April 2011. The stock price can only move at trades, and it does not change on 85.5% of the trades. There are 13,829 trades during regular trading hours and 25 off-hours trades. The average duration between trades during regular hours is 1.69 s; the standard deviation is 5.11 s, and the maximum duration is 104 s. Seventy per cent of the durations are less than one second. Panel B graphs the average number of arrivals during 15 min trading intervals during regular trading hours from 3 January to 31 March. Trading is more intense at the opening and closing bell than the midday lull from noon to 2pm. Panel C plots the histogram of returns. 85% of the trades did not have a price change, which results in the dominate spike at zero in the histogram in Panel C. The much smaller side lobes are an artefact of price being measured to the nearest $0.01. Returns occur between the lobes and are not visible due to the dominate bar at 0%. Engle (1982) noted that sampling from a marked process can result in inducing features that are not apparent in the non-sampled data. Panel D shows the histogram from sampling at 30 min intervals, and looks more like a t-distribution than Panel C. Panels E and F plot the autocorrelation function of the squared returns, which is a measure of stochastic volatility. Because there are one million trades, there is statistical significant autocorrelations for the raw price process in Panel E, but they are small. The sampled process illustrates more autocorrelation, in Panel F, due to the diurnal trading. The lag 12 for the sampled process in Panel F is 6 h, while lag 12 for the original process in Panel E is less than one minute.

Figure 3. Price process for PG. (A) Prices for 2011-04-01; (B) Average number of arrivals per 15 min; (C) Return histogram all trades; (D) Return histogram 30 min sampling; (E) ACF of return squared all trades; (F) ACF of returned squared 30 min sampling.

Figure 4 compares the exponential HB smoother that attenuates high-frequency spectral coefficients with two other models: (1) an exchangeable prior with equal prior variances for the spectral coefficients (γ = 0) and (2) cubic B-splines (De Boor 1978; Eubank 1988; Eilers and Marx 1996) with an exchangeable prior on the coefficients in place of the cosine series. The three models are estimated to 91,122 trades of UGI over 60 trading days. Panels A and C plot A(r) for the lagged percent returns, and Panels B and D plot B(t_d) for the time of day. Panels A and B plot the cosine series models, and Panels C and D plot the B-splines with two choices of knots. The spectral expansions use 30 cosine functions and the splines models use (1) knots at 1 percentiles for both functions and (2) customised knots for A(r) and knots at 12 min intervals for B(t_d). The B-splines were computed from the bs() function in the R package splines.

Figure 4. Cosine series and B-spline for duration model of UGI. (A) 30 term cosine series of lagged percent return. (B) 30 term cosine series of time of day. (C) B-splines with knots at 1 percentile points (35 internal knots) and customised (12 internal knots) for lagged per cent return. (D) B-splines with knots at 1% points (99 internal knots) and 12 min intervals (32 internal knots) for time of day.

Panels A and C of Figure 4 show that the (log) arrival rate is minimum when the previous return was zero, and the arrival rates increase for larger absolute returns, but this effect attenuates beyond ±1%. The exponential HB smoothing prior (solid blue line) is smoother than the exchangeable prior (dashed red line) in Panel A. Estimating A(r) with B-splines is challenging because the return data are long-tailed and clumpy: 34% of the returns are within ±0.01% and only seven trades had absolute returns larger than 1%, while the minimum return was −3.8% and the maximum return was 4.5%. The 5 percentiles from 5% to 95% have seven unique values: −0.0318, −0.0312, −0.0305, 0, 0.0305, 0.0312, .0318, and the 1 percentiles from 1% to 99% result in 35 unique values due to the clusters. The unique 1% percentiles are the internal knots for the solid blue line in Panel C, while the red dashed line used 12 customised internal knots at ±0.01, ± 0.03, ± 0.04, ± 0.06, ± 0.07, ± 1, which were chosen by trial and error to better match Panel A. The cubic splines are sensitive to the number and placement of knots, especially outside of ±1%, which highlights a qualitative difference in the cosine basis and B-splines. B-splines are a local regression method where only observations that are close to the abscissa are used to estimate the ordinate because the basis functions are zero outside of subintervals. In contrast, spectral analysis is a global regression method: the cosine functions are non-zero except for integer multiples of π. The coefficients for the B-splines that are outside ±1% in Panel B are effectively estimated with very few observation, which makes the B-spline sensitive to the knots. In contrast, all of the data are used to estimate the spectral coefficients in Panel A. Because the observed durations after the seven trades outside of ±1% are zero, the plateaus in Panel A better represent the data than the cubic spline in Panel C. An alternative approach is to delete these outliers and only estimate the function between ±1%. Even though large price movements occur very infrequently relative to the number of trades (7/91,122), they are not rare relative to the number of trading days (7/60). We decided to keep the full data to model these ‘black swan’ events (Taleb 2010) and to exploit fully the information in high frequency trading.

Panels B and D of Figure 4 show the bath-tub shape for inter-day arrivals B(t_d) where there is a lull in trading around noon and a surge later in the afternoon. The two priors for the cosine series result in similar shapes with the exchangeable prior having more high-frequency oscillations than the exponential HB smoother. The B-splines are more robust to the choice of knots in Panel D than Panel C because the arrival data are dense during regular trading hours and do not have clusters. Unsurprisingly, the B-spline is smoother with fewer knots (dashed red line) than with more knots (solid blue line). B-splines are sensitive to the choice and locations of knots, and cross-validation methods can provide guidance on their choice. The exponential HB smoothing prior automatically determines the smoothness in spectral analysis.

Figure 5 displays the estimated functions A(r) for lagged return in the trading intensity of Equations (7)(7) $\lambda (t_d|x,r)=\text{exp}[x^{\mathrm{\prime }}\beta +A(r)+B(t_d)]$(7) and (11)(11) $\begin{array}{rl}{\lambda }_0(t_d)& =\exp [{\beta }_1\mathrm{Mo}+{\beta }_2\mathrm{Tu}+{\beta }_3\mathrm{Wd}+{\beta }_4\mathrm{Th}+{\beta }_5\mathrm{Fr}+B(t_d)]\text{for TST}\\ {\lambda }_1(t_d)& =\exp [{\beta }_0+{\beta }_{6}D+A(r)]\text{ for TIT}\end{array}$(11) . The panels are for the four stocks, and each panel plots the functions in the homogeneous and mixture model: only the intensity function for TIT includes A(r). The centre lines are the posterior mean and the bounds, which are very narrow, are 95% highest posterior distribution (HPD) intervals. The plots indicate that the functions vary across assets. Further studies are needed to determine the extent that these differences are systematic or random. For instance, they may depend on the industry and the stage of the investment/economic cycle. The estimated functions for the homogeneous and mixture models are similar, with the mixture model tending to have taller peaks and lower valleys – except the valley for WMT.

Figure 5. Nonparametric estimates of lagged return in duration intensity function. Posterior means and 2.5 and 97.5 posterior percentiles.

Figure 6 plots the functions B(t_d) for time of day in the intensity function for arrivals in Equations (7)(7) $\lambda (t_d|x,r)=\text{exp}[x^{\mathrm{\prime }}\beta +A(r)+B(t_d)]$(7) and (11)(11) $\begin{array}{rl}{\lambda }_0(t_d)& =\exp [{\beta }_1\mathrm{Mo}+{\beta }_2\mathrm{Tu}+{\beta }_3\mathrm{Wd}+{\beta }_4\mathrm{Th}+{\beta }_5\mathrm{Fr}+B(t_d)]\text{for TST}\\ {\lambda }_1(t_d)& =\exp [{\beta }_0+{\beta }_{6}D+A(r)]\text{ for TIT}\end{array}$(11) during trading hours. In the mixture model, only TST includes B(t_d). These are more consistent with bursts of activity after the opening bell an before closing with a minimum around the lunch hour.

Figure 6. Nonparametric estimates of time of day in duration intensity function. Posterior means and 2.5 and 97.5 posterior percentiles.

Figure 7 plots the densities for the return distributions. The ‘prior’ density g₀ in Equation (10)(10) $\begin{array}{rl}{f}_R(r)& =\frac{g_0(r)\text{exp}\left\{{\int }_a^{r}h(s){[Z(s)]}^2\mathrm{ds}\right\}}{{\int }_a^b{g}_0(w)\text{exp}\left\{{\int }_a^{w}h(s){[Z(s)]}^2\mathrm{ds}\right\}\mathrm{dw}}\\ Z(s)& =\sum _{k=0}^K\sqrt{2}\text{cos}[\mathrm{k\pi }{G}_o(s)]{\theta }_{f,k}\text{for}a\le s,r\le b\\ h(s)& =\frac{1-\text{exp}(100s)}{1+\text{exp}(100s)}\end{array}$(10) is Cauchy with mode zero, so the tails for the estimated return densities are similar to those of the Cauchy density. The x-axes are truncated at ±0.15% in the plots. The densities are nearly symmetric and leptokurtic.

Figure 7. Return distributions. Posterior means and 2.5% and 97.5% percentiles.

Table 2 displays the estimated model for off-hours trading. UGI has the least amount of off-hour trading, and CITI has the most. The probability that a day has pre-open trading ranges from 0.05 for UGI to 0.99 for CITI. The lag between the 4am open and the first trade ranges from 3 to 5 h. In post-close trading, the expected number of trades (Poisson parameter) ranges from 0.2 for UGI to 369 for CITI. The Poisson rate of trading ranges from 0.0002 for UGI to 0.1035 for CITI in pre-open trading and 0.0011 for WMT and 0.0257 for CITI in post-close trading.

Table 2. Estimated off-hour trading models.

	Posterior Mean				Posterior SD
Pre-Open Trading	CITI	PG	UGI	WMT	CITI	PG	UGI	WMT
P(Pre-Open Trades)^a	0.9878	0.9390	0.0488	0.8780	0.0122	0.0265	0.0239	0.0364
Lag: Mean (Hours)^b	3.0580	4.2428	4.7493	4.4616	0.0016	0.0079	0.1098	0.0080
Lag: SD (Hours)^b	0.1281	0.6050	0.4390	0.5772	0.0097	0.0467	0.0946	0.0460
Poisson Arrival Rate^c	0.1035	0.0043	0.0002	0.0041	0.0004	0.0001	0.0002	0.0001
Post-Close Trading	CITI	PG	UGI	WMT	CITI	PG	UGI	WMT
Poisson Mean Trades^d	369.12	5.4077	0.2000	9.7671	2.4783	0.3000	0.0476	0.4031
Poisson Arrival Rate^e	0.0257	0.0009	0.0023	0.0011	0.2620	0.0004	0.0009	0.0004

^a P(Pre-Open Trades) is the probability that a day has pre-open trades.

^b The lag between 4am and the first trade has a normal distribution.

^c Trades after the first trade arrival according to a homogeneous Poisson process.

^d The daily total of post-close trades have a Poisson distribution.

^e Given the total number of daily post-close trades, trades arrive according to a homogeneous Poisson process.

Next, we illustrate simulating the price process by using the first 60 days of data for UGI to estimate the model and the next 20 days as a hold-out sample. The simulations use a different draw of the parameters from the posterior distribution as generated from the MCMC algorithm for each simulated price process. We simulated 5000 CIPP over 20 days after the training period. Figure 8 summarises the simulation. Panel A plots the distribution of the simulated price history at the daily close for 20 days and the observed series, and Panel B summarises the price distribution in 15 min increments for the first day in the holdout sample or day 61. The coloured bands are percentiles in 5% increments while the minimum and maximum areas represent 95% prediction intervals. The long dashed line is the mean for geometric Brownian motion, and the dotted lines are the 95% prediction intervals for geometric Brownian motion. The simulated price process from CIPP and Brownian motion are nearly the same over 20 days in Panel A, but they diverge in Panel B because CIPP better reflects the micro-structure of high-frequency data. For UGI it appears that daily sampling of CIPP at the daily close can be closely approximated by geometric Brownian motion, while sampling 15 min time intervals shows greater deviation.

Figure 8. Simulated CIPP in validation sample and option prices for UGI. The colour bands are intervals that contain 5% of the simulated values from 5% to 95%. The outmost limits are 2.5% and 97.5%. (A) Simulated price process in validation sample; (B) Simulated price process in day 61; (C) Option prices for simulated CCIP; (D) Ratio of CCIP to black-sholes option prices.

Panels C and D in Figure 8 use the simulated CIPP to price call options and compare them to Black–Scholes for a European call option (Hull 2009). We assume that the discount rate is 5% and the strike price K_D with duration of D days is: $K_D=S_{60,C}\text{exp}\left(\mu +\frac{{\sigma }^{2}D}{2}+0.1\sigma \sqrt{D}\right)$where µ and σ are estimated from the returns in 15 min gaps, and S_60,C is the closing price in the last day of the calibration period: strike prices increase with durations. Panel C gives the distribution of option prices based on the simulated CIPP (solid line) along with the Black-Sholes pricing (dashed line). The mean option price from CIPP is slightly less than Black-Sholes. Panel D plots the ratio of the CIPP option distribution to the Blacks-Sholes price, where the mean of the ratio is the black line. Overall, option pricing using the CIPP and geometric Brownian motion are similar when durations are 1 to 20 days when computed at the daily close, and they have greater differences within a day.

8. Conclusion

We propose semi-parametric Bayesian analysis for compound, inhomogeneous Poisson processes and apply it to ultra-high frequency trading data. Trades arrive according to an inhomogeneous Poisson process. The intensity function is a flexible function of time of day during regular trading hours and the lagged return from the previous trade, which are modelled with Gaussian process priors. Time-of-day effects have a bathtub shape with frequent arrivals at the open and close of regular trading and a lull around midday. Arrival rates tend to increase after large changes in prices and tend to decrease after small changes; however, the behaviour depends on the stock and a more detailed study is required. The parametric model for the intensity function includes day-of-week effects and the lagged duration times, which have carryover effects: long durations tend to follow long durations, and short durations follow short ones. We use a logistic-Gaussian process that is unimodal for the return distribution. The estimated density is asymmetric and leptokurtic with a mode at 0.

Spectral representations of the Gaussian processes linearise the analysis, simplify the computations, and facilitate shape constraints for the return density. The spectral coefficients have exponential HB smoothing priors with hyper-parameters that control the tradeoff between the prior means and the data and that control the smoothness of the functions by attenuating high-frequency oscillations that are due to noise. These smoothing parameters are estimated in a fully Bayesian analysis that avoids cross-validation. We use a cosine basis, but researchers could also use other basis functions such as B-splines or Hermite polynomials. We also propose a mixture model where the specification of the intensity function for the two components helps to classify observations into time sensitive trades and time insensitive trades. We conceptualise the first group as trades that are generated by human investors (time sensitive trades), and the second group as trades that are software generated (time insensitive trades). The probability of time insensitive trades is a logistic function of the contemporaneous return and order size of the trade. In the empirical example, time insensitive trades, tend to have larger order sizes and returns than time sensitive trades, though further research is required to verify the result. We also compare simulations of the estimated process to geometric Brownian motion. Over long time period (20 days), the two are very close, while they deviate over shorter periods within the day. The model could be extended to have different diurnal functions for the days of the week or a hierarchical Bayes model for separate functions for each day.

The semi-parametric CIPP methodology can be used for other research purposes. For example, the intensity function could include indicator variables for event analysis. For instance, what is the impact on the arrival of trades before and after unexpected earnings announcements, changes in interest rates, or major news events? Increased activity before the event may indicate insider trading. Trades from high-frequency trading firms could be flagged, and their impact on the price discovery process ascertained. The proposed mixture model could be modified to include more components or to investigate other structural classification issues. The purpose of the paper was to develop the methodology for semi-parametric, CIPP. A more complete empirical investigation of the price discovery processes would require a larger number of stocks and time periods.

Disclosure statement

No potential conflict of interest was reported by the authors.

Notes

1 The rounding assumption does not seem to make material difference on the results.

2 λ(t) = exp(–6 + 2cos[2πt/s]) where t is cumulative number of seconds and s = 20*10*3600 for 10 h in regular trading per day, and 3600 s in an hour.

3 We assume that the prior distributions for the coefficients of the B-spline are independent with common standard deviation τ. Technically, B-splines are not an orthonormal basis, so the Karhunen-Loève representation theorem does not hold; however, the resulting B-spline function has the covariance kernel in Panel B of Figure 2.

4 Some of the spectral coefficients did not pass CODA’s diagnostics, which may be due to applying single parameter tests simultaneously to a large number of parameters. Spectral coefficients near zero were flagged more often; however, they had little impact on the estimated functions.

5 Theoretically, the normal distribution is truncated between the open of pre-open trading and the open of regular trading. In practice the first trades tend to cluster around 8 am, which is sufficiently far from the endpoints that the truncation is not an issue. If it is an issue, then one could use a truncated normal distribution.