Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
92,359 characters · 22 sections · 53 citation commands
A Flexible Stochastic Conditional Duration Model
Duration models for financial transactions describe the irregular timing of trades, or other events such as price changes. They are useful because trading intensity is one measure of market liquidity. And unlike count models, duration models use all the information in trading times. They also shed light on market microstructure phenomena.
Modelers have to confront the fact that some trades are related to others and that related trades are nearly simultaneous. We will argue that widely accepted rules for aggregating seemingly related trades into clusters, together with unsuitable parametric duration distributions, mislead inference in important ways. The most commonly used rule treats two trades as related if they are executed, to within recording precision, simultaneously. In this paper, as in most studies, transaction times are truncated to the second, and in this context we will call it the {\em same-second} aggregation rule.\footnote{Many datasets, especially more recent ones, feature millisecond or finer recording precision. The issues raised in this paper apply widely, but since the unfortunate consequences of deterministic aggregation rules are qualitatively different when recording precision is finer, we leave this case for future research.} The rule is usually correct, case by case, but in a reasonably sized sample from a liquid market, it is likely that many pairs of unrelated trades occur within the same clock second by happenstance.
Our paper makes three main contributions. First, we put into our model the uncertainty about whether each duration is a {\em cluster} duration---one between related trades---or a {\em regular} duration---one between unrelated trades. Not surprisingly, the posterior probability that a duration recorded as $0s$ (zero seconds) is due to coincidence, and therefore regular, is always low. But it is never zero, and the posterior mean of the number of these coincidences is much larger than its posterior standard deviation. We find that the conditional hazard function for regular durations varies smoothly near zero, in contrast to the abrupt changes found in most studies, which we argue are artifacts arising from classifying {\em all} durations recorded as $0s$ as cluster durations.
Our second contribution is a new flexible model for regular durations. We identify some undesirable features of existing parametric conditional duration densities and propose a new family of conditional distributions that is flexible, but also amenable to shrinkage towards the exponential distribution, which we will argue is a theoretically appealing first-order model.
Our third contribution is computational, and promotes highly numerically efficient posterior simulation. We use the HESSIAN method, introduced by McCa12, to draw, in a single Gibbs block, the full sequence of state variables describing trade intensity, together with various parameters. To date, the HESSIAN method has only been applied to non-Gaussian state-space models with parametric distributions for observed variables and homogenous state transitions. Here we use it in a model with a flexible distribution for observables and heterogenous state transitions. The numerical efficiency we achieve is in fact considerably higher than that achieved using auxiliary mixture model methods, which rely on special features of parametric distributions.
Two other features of our model are appealing and uncommon, but not original. First, we estimate the regular diurnal (time of day) pattern of trading intensity jointly with other features of the model.\footnote{ VereRodrEspa02 and BrowVann13 jointly estimated the diurnal patterns and the parameters for ACD models using respectively a semi-parametric approach and MCMC methods within a Bayesian framework. } Second, the latent state process is an irregularly sampled Ornstein-Uhlenbeck (OU) process, rather than a homogenous autoregressive process; autocorrelations depend on the elapsed time, rather than the number of intervening trades, between two durations.
We now illustrate in detail the problem with treating all durations recorded as $0s$ as cluster durations, a practice that goes back to the seminal paper of EnglRuss98. Table (ref) shows descriptive statistics for two duration series, cleaned as described in Section (ref). Recorded durations are differences between two trading times, both truncated to the second. The last six columns give the percentage of durations recorded as $0s$, $1s$, $2s$, $3s$, $4s$, and $5s$. The values from $1s$ to $5s$ vary smoothly, but far more durations are recorded as $0s$ than are recorded as $1s$.
This zero inflation clearly needs to be addressed. It arises because of nearly simultaneous related trades. These often occur when a market order is matched against, and filled with, several limit orders on the opposite side of the market. It may also occur if many traders submit limit orders to be executed at a round price, as suggested by VereRodrEspa02, or if important news synchronizes a flurry of trading, or if traders use algorithmic trading strategies that can be triggered by another trade. We do not try to distinguish different causes of related trades.
In fact, the amount of zero inflation is greater than it might appear: the duration between two trades will be recorded as $0s$ if the second trade occurs during the {\em remainder} of the same clock second as the first, but as $1s$ if it occurs at {\em any} time during the next clock second. Extrapolating percentages from $1s$ to $5s$ back to $0s$ and dividing by two gives us a rough idea of the percentage of all durations that are regular and recorded as $0s$: about 5% for the RY series and 7% for the POT series. Of course, these represent a much larger percentage of regular durations, around 14% and 22%, respectively. While it is likely that a large majority of durations recorded as $0s$ are cluster durations, it would be very surprising if all were.
The same-second aggregation rule amounts to removing all durations recorded as $0s$ from the sample. The result is a truncated sample of regular durations. The fact that the truncation is at both the mode and lower bound of the distribution is particularly unfortunate. Trading intensity is understated, and the understatement varies with trading intensity: when it is high, there are more short regular durations and more spurious aggregation of unrelated trades.
Another aggregation rule, proposed by GramWell02, aggregates any sequence of transactions within the same clock second where prices are non-decreasing or non-increasing. We will call this the {\em GW} aggregation rule. Figure (ref) shows histograms of durations that are classified as regular by this rule. Not enough of the durations recorded as $0s$ are classified as regular for compatibility with a smooth density near zero. This is not surprising, as many pairs of unrelated trades will feature price changes of the same sign by coincidence.
While the spurious aggregation of unrelated trades is our main concern with these rules, we also question whether even error-free aggregation would be desirable. If a market order matches nine limit orders, five times as many traders are getting their orders filled than if it matches a single limit order. Treating these cases as equivalent may understate liquidity, as perceived by traders.
Our model is an example of what Engl02 calls a {\em multiplicative error} model, where the scale of the conditional duration distribution depends on the history of the process, as well as observables such as time of day, but the shape does not.\footnote{ For recent surveys on the analysis of high-frequency financial durations using multiplicative error models see Pacu08, Haut12 and BhogVari18. } We will call the conditional distribution of duration divided by scale the {\em normalized} conditional distribution; typically, it has unit mean. Two basic models are the data-driven Autoregressive Conditional Duration model (ACD), introduced by EnglRuss98, where the scale depends deterministically on past durations; and the parameter-driven Stochastic Conditional Duration model (SCD) introduced by BauwVere04, where it is a latent stochastic process. BauwGiot00 propose a logarithmic version of the ACD model, avoiding parameter restrictions; SCD models usually feature a similar logarithmic specification. BauwVere04 provide empirical evidence favouring the SCD model over a similar ACD model in the four data sets they analyzed.
The most commonly used conditional distributions in SCD and ACD models are the exponential, and two generalizations: the gamma and the Weibull EnglRuss98, BauwGiot00, BauwVere04, FengJianSong04, StriForbMart06, MenKolkWirj15. EnglRuss98 and others conclude that the exponential is too inflexible and favour the Weibull or gamma. However, these conclusions favouring the Weibull or gamma come largely from studies in which durations recorded as $0s$ are removed from the sample, and the density and hazard functions of these distributions feature extreme behaviour near zero. Except for the knife-edge special case where they reduce to an exponential, which has a constant hazard function, the hazard function is either zero at a duration of zero and increasing; or infinite at a duration of zero and decreasing.\footnote{ The generalized gamma and the Burr distribution generalize the gamma and Weibull distribution. They were proposed as conditional distribution for ACD models by Lund99 and GramMaue00. These distributions, unlike the gamma and the Weibull, allow for non-monotonic hazard functions, but they retain the property that their hazard function is bounded away from zero and infinity only for very special cases.}
We argue that extreme variation of the hazard function of regular durations near zero is implausible. In queueing theory, a simple model for arrival times (of, say, customers at an ATM) is the Poisson process. It is reasonable when there are a large number of potential customers, acting independently and homogenously in time, and the probability of any {\em given} customer arriving in a given time interval is much smaller than the probability of {\em some} customer arriving in that interval. In a Poisson process, durations between arrivals are exponentially distributed. The constant hazard function of the exponential makes it {\em memoryless}: the probability of an arrival in the next minute does not depend on how long you have been waiting.
Of course, trading intensity in financial markets varies over time. But after conditioning on relevant predictors and latent states measuring trading intensity, we would expect the distribution of regular durations to be not far from an exponential---its hazard rate a function of this conditioning information---due to the large number of unrelated potential traders and the fact that a regular duration is, by definition, the time interval between unrelated trades.
The hazard functions of mixtures of exponentials are bounded away from zero and infinity,\footnote{ DeLuGall04 used a mixture of two exponentials in ACD models and found that this specification provides a better fit than a Weibull distribution. DeLuGall09 again use a mixture of two exponentials but allow mixture weights to depend on observable market activity. } but they are decreasing---see BarlMarsPros63---which is a restrictive feature. Moreover, not all decreasing hazard functions can be easily captured by mixtures of exponentials. In simulations we do not report, we find that adding mixture components after the second yields little: only components with the largest and smallest hazards are important, in the sense that the posterior distributions of the weights of other components are highly concentrated near zero.\footnote{ Other mixture distributions have been proposed for the normalized conditional density. WirjKolkMen13 did a Bayesian analysis of the SCD model, with leverage, using three types of two-component mixtures: two exponentials, two Weibulls and two gammas. }
The rest of the paper is organized as follows. We describe the transaction data we analyze in Section (ref), our SCD model in Section (ref) and our methods for posterior simulation in Section (ref). In Section (ref), we conduct an artificial data experiment to test for the correctness of our posterior simulator and then illustrate our methods in an application featuring two equities traded on the Toronto Stock Exchange. We conclude in Section (ref).
Our data comes from the Tick Data historical database of securities listed on major stock exchanges. We use transaction data for two equities traded on the Toronto Stock Exchange (TSX): the Royal Bank of Canada (RY), and the Potash Corporation of Saskatchewan (POT). For each equity, we observe transactions over five consecutive trading days, from March 17 to March 21, 2014.
The trading hours of the TSX are 9:30 am to 4:00 pm, Monday to Friday. We ignore the {\em pre-open} trading session and the {\em extended} trading session after closing, and only use data from the {\em continuous} trading sessions. These data include the time of each transaction (the time stamp), the price, the volume (in number of shares), the type of trading session (pre-open, continuous, or extended), and various indicator variables specifying, among other things, whether the transaction was delayed or subsequently corrected. Trading times are truncated to the second; in other data sets, transaction times are recorded to higher levels of precision.
BrowGall06 argue that careful data cleaning is an important part of high-frequency duration analysis. We follow standard practice to remove observations that are obviously erroneous or result from atypical market conditions, such as at opening and closing. We delete entries identified by the exchange as incorrect, corrected, delayed or canceled. Then we remove transactions with aberrant prices using the method suggested by BrowGall06. Descriptive statistics of the cleaned data are reported in Table (ref). For each duration series, we report the number of observations, followed by the sample mean, standard deviation, maximum, and coefficient of variation. An exponential has a coefficient of variation of 1, so the empirical distributions are considerably overdispersed relative to the exponential. The second part of the table reports the percentage of durations recorded as $0s$, $1s$, $2s$, $3s$, $4s$, and $5s$.
The left panel of Figure (ref) shows the cumulative number of trades of POT during each of five consecutive days of trading, from March 17 to 21, 2014. The right panel zooms in on a typical half-hour period to give finer detail; it shows the cumulative number of trades on March 17, 2014 between 11:30am and noon. The figure illustrate many widely known stylized facts about durations. Trading intensity varies over time. Part of the variation is predictable given time of day: there are more trades near the beginning and end of the trading session than in the middle. The remaining, stochastic, part is highly persistent. Trades often arrive in clusters, within a very short interval and without a marked difference in trading intensity before and after the cluster.
Our model builds on the multiplicative-error stochastic volatility model of BauwVere04. It differs in three important ways. First, it has a second state variable, one which indicates which durations are cluster durations and which are regular. Since this replaces the usual practice of aggregating trades into clusters before analysis, our model is intended for unaggregated data. For the purpose of comparison, we also define a special case of our model where cluster durations have probability zero, for samples in which durations have been aggregated into clusters. We will call the special case the {\em regular-duration} model and the general model the {\em all-duration} model. Placing the classification of trades within the model makes it possible to infer that the number of regular durations that happen to be recorded as $0s$, while small, is not zero. We will see that this overcomes the problems arising from the spurious aggregation of unrelated trades. Second, we use a flexible distribution for the scale-normalized duration distribution. Our approach makes it easy to choose a prior distribution for this distribution that shrinks towards an exponential distribution, which we have argued is a sensible first-order conditional model for regular durations. Third, we jointly estimate the diurnal (time-of-day) pattern along with everything else, which gives a better accounting of the uncertainty about how trading intensity evolves. We describe a family of prior distributions to complete the model.
As with most models in the literature, time is continuous, which may not be suitable for data where trading intensity is high and transaction times are truncated to the second. In Section (ref) below, we describe a discrete-time version of the model that is more suitable for such data.
We observe transaction times over $D$ trading days in the interval $[t_\mathrm{open}, t_\mathrm{close}]$, where $t_\mathrm{open}$ and $t_\mathrm{close}$ are the opening and closing times. All times of day are measured in seconds after midnight. For each day $d=1,\dots,D$, denote the sequence of transaction times by $t_{d,0},t_{d,1},\dots,t_{d,n_d}$ and construct the durations $y_{d,i} \equiv t_{d,i} - t_{d,i-1}$, $i=1,\dots,n_d$.
The conditional distribution of each $y_{d,i}$ depends on the values $s_{d,i} \in \{0,1\}$ and $x_{d,i} \in \mathbb{R}$, at time $t_{d,i-1}$, of two latent states. The state $s_{d,i}$ is a mixture component indicator: $s_{d,i} = 0$ indicates that $y_{d,i}$ is a cluster duration; $s_{d,i} = 1$, a regular duration. The state $x_{d,i}$ gives the trading intensity at $t_{d,i-1}$, but only if the duration is regular. The conditional density of $y_{d,i}$, given $s_{d,i}$ and $x_{d,i}$, is
where $p_0$ and $p_1$ are density functions on $[0, \infty)$.
At each day $d$, the indicator process $s_{d,i}$ is first-order Markov and stationary, with
where $\xi_0 \equiv (1-\xi_{11}) / (2-\xi_{00}-\xi_{11})$ and $\xi_1 \equiv (1-\xi_{00}) / (2-\xi_{00}-\xi_{11})$ by stationarity. The trading intensity process $x_d(t)$, $t \in [t_\mathrm{open},t_\mathrm{close}]$, is the sum of a deterministic function $m(t)$ describing a diurnal pattern---equation (ref) below---and a zero-mean OU process. Sampling the $x_d(t)$ process at all transaction times gives the values $x_{d,i} \equiv x_d(t_{d,i-1})$; then the discrete time process $x_{d,i}$ is first order autoregressive, but not homogenous due to irregularly spaced trading times:
where $\phi \geq 0$ is the mean reversion parameter, and $\sigma$ is the marginal standard deviation parameter of the OU process.
To model the diurnal pattern, we follow EileMarx96, LangBrez04 and BrowVann13 and specify the deterministic process $m(t)$ as a cubic B-spline function, a piecewise polynomial defined on a set of $M$ equally spaced knots $t_\mathrm{open} = \kappa_1 < \cdots < \kappa_M = t_\mathrm{close}$. Knots $\kappa_1$ and $\kappa_M$ have multiplicity 4 and the rest have multiplicity 1. This gives an expansion that is the following linear combination of $L=M+2$ B-spline basis functions:
where $B_l(\cdot)$ denotes the $l$-th B-spline basis function, a local cubic polynomial, and $\delta_l$ its coefficient. The basis functions depend on the location and multiplicity of the knots. Setting the multiplicity of the knots $\kappa_1$ and $\kappa_M$ to 4 makes $m(t_\mathrm{open}) = \delta_1$ and $m(t_\mathrm{close}) = \delta_L$. This makes $\delta_1$ and $\delta_L$ easier to interpret. See Boor_1978 or Dierckx_1993 for more on B-splines.
We now specify the component densities $p_0(\cdot)$ (cluster) and $p_1(\cdot)$ (regular). We specify $p_0(\cdot)$ as the following mixture of two exponentials, with hazard parameters $\lambda_1$ and $\lambda_2$:
Given the nature of cluster durations, $p_0(\cdot)$ should concentrate most of its probability near zero; the prior distribution over $\lambda_1$, $\lambda_2$ and $\pi$ should be chosen so that this is true with high prior probability.
Unlike cluster durations, regular durations depend on market conditions, as captured by the intensity state $x_{d,i}$. Our specification of $p_1(\cdot)$ conforms to the multiplicative error structure proposed by BauwVere04, which gives the regular duration as $y_{d,i} = e^{x_{d,i}} \epsilon_{d,i}$, where $\epsilon_{d,i}$ is the contemporaneous value of a non-negative iid process with $E[\epsilon_{d,i}] = 1$. The $\epsilon_{d,i}$ and $x_{d,i}$ processes are independent. We denote the density of $\epsilon_{d,i}$ by $p_\epsilon(\cdot)$ and call it the normalized duration density. Thus, the conditional density of regular duration $y_{d,i}$, given $x_{d,i}$, is
Given $s_{d,i}=1$, the quantity $e^{x_{d,i}}$ is the conditional mean of $y_{d,i}$. It gives only the scale of the distribution, the shape being determined by $p_\epsilon(\cdot)$. We argued that some commonly used duration distributions are unsuitable because of their restrictive or implausible hazard functions. To complete the model specification, we will now propose a new flexible family of normalized densities.
We adopt an approach similar to one described in FerrStee06, where the cdf $P_\epsilon(\cdot)$ of a flexible density $p_\epsilon(\cdot)$ is constructed as $ P_\epsilon(\epsilon) \equiv G\big(F(\epsilon)\big), $ where $F(\cdot)$ is a parametric continuous cdf, with density $f(\cdot)$, having the same support as $P_\epsilon(\cdot)$; and $G(\cdot)$ is a flexible continuous cdf on $[0,1]$, with density $g(\cdot)$. The cdf $P_\epsilon(\cdot)$ can be viewed as a perturbed version, depending on $G(\cdot)$, of the parametric cdf $F(\cdot)$. When $G(\cdot)$ is uniform on $[0,1]$, there is no perturbation and $P_\epsilon(\cdot) = F(\cdot)$. As noted by FerrStee06, distributions defined in this way cover the entire class of continuous distributions, since any such $P_\epsilon(\cdot)$ can be constructed for a suitable choice of $G(\cdot)$. The construction implies the normalized density
We argued that the exponential distribution was a theoretically promising first order approximation of a conditional duration distribution, and so we specify $F(\cdot)$ as an exponential, in the hope of capturing realistic hazard functions using a distortion $G(\cdot)$ with a small number of terms. The hazard rate of $F(\cdot)$, which we denote $\tilde \lambda$, will be substituted out using the scale normalization condition $E[\epsilon] = 1$. At the same time, we want $P_\epsilon(\cdot)$ to be flexible, allowing for large departures from the exponential if the data warrant it. For this reason, we choose a flexible functional form for $g(\cdot)$, a $J$'th order Bernstein polynomial with positive coefficients adding to one. This is a $J$-component mixture of beta densities, each with integer-valued shape parameters adding to $J+1$; coefficients of the first and last terms determine $g(0)$ and $g(1)$, respectively. Bernstein polynomials can approximate any continuous density on $[0,1]$ arbitrarily closely (in $\sup$ norm) for $J$ sufficiently large.\footnote{ Let $F(\cdot)$ be a bounded continuous function on $[0,1]$. The Bernstein polynomial of order $n$ of the function $F(\cdot)$ is defined as \[ B_n^F(x) = \sum_{k=0}^n F\left(\frac{k}{n}\right)\binom{n}{k} x^k (1-x)^{n-k}. \] Lore53 shows that $\lim_{n\rightarrow\infty} \sup_{x \in [0,1]} |B_n^F(x) - F(x)| = 0$. } The component probabilities are the coefficients of the Bernstein basis polynomials. Specifically, $ g(z) = \sum_{j=1}^J \beta_j \operatorname{Beta}(z \,\vert\, j,J-j+1), $ where $\sum_{j=1}^J \beta_j = 1$, $\beta \equiv (\beta_1,\ldots,\beta_J) \geq 0$, and $\operatorname{Beta}(z \,\vert\, a, b)$ denotes the beta density with shape parameters $a$ and $b$, for $a,b > 0$. The $J$ Bernstein basis polynomials of order $J$ form a partition of unity, so that if $\beta = (1/J,\ldots,1/J)$, then $g$ is uniform and we get back the original exponential density, $p_e(\cdot) = f(\cdot)$. This makes it easy to choose a prior distribution for $\beta$ so as to centre the induced prior for $p_\epsilon(\cdot)$ around the exponential distribution and to control the amount of shrinkage towards it. We treat the order $J$ of the Bernstein polynomial as fixed, not as a parameter to estimate, due to the difficulty of computing marginal likelihoods for models with two different state variables. Instead of estimating $J$, we will compare results over different values of $J$. For a discussion of Bayesian nonparametric density estimation using Bernstein polynomials with unknown $J$, see Petr99a,Petr99b and PetrWass02.
Substituting the exponential distribution function $F(\epsilon) = 1-e^{-\tilde \lambda \epsilon}$ and the above expression for $g(z)$ into equation (ref) gives $p_\epsilon(\cdot)$ as a polynomial in $e^{-\tilde \lambda \epsilon}$, which we write explicitly as the following linear combination of exponential densities:
where $\alpha \equiv (\alpha_1,\dots,\alpha_J) = T\beta$ and $T$ is a $J\times J$ matrix with $T_{j,J-i+1} = \binom{J}{j}\binom{j}{i}(-1)^{j-i} \, (i/j)$ if $j\leq i$ and $0$ otherwise. Non-negativity of $\beta$ ensures the non-negativity of $p_\epsilon(\cdot)$, but as the $\alpha_j$ are not necessarily all non-negative, $p_\epsilon(\cdot)$ is not necessarily a mixture of exponentials. It is easy to check that the hazard function for $p_\epsilon(\cdot)$ is bounded away from zero and infinity, and its limiting value as $\epsilon$ goes to infinity is the hazard parameter $\tilde \lambda$ of the exponential distribution with cdf $F(\cdot)$. The scale normalization condition $\mathrm{E}[\epsilon]=1$ gives \( \tilde \lambda = \sum_{j=1}^J \frac{\alpha_j}{j}, \) which we use to substitute out $\tilde \lambda$ from the expression for $p_\epsilon(\cdot)$, freeing us from having to impose restrictions on $\alpha$ or $\beta$.
We now illustrate the flexibility we achieve with just a few terms. Figure (ref) shows examples of decreasing, increasing and non-monotonic hazard functions that can be captured with $J=3$ components. The solid lines are the density and (constant) hazard functions of an exponential with mean equal to one. The other density and corresponding hazard functions are for various values of $(\beta_1,\beta_2,\beta_3)$. Taking into account the adding-up constraint, there are two degrees of freedom, the same as a unit-mean mixture of two exponentials. The dashed lines show pairs of density and hazard functions where the hazard is monotonic; dash-dotted lines, pairs where it is not monotonic. Recall that all mixtures of exponentials have a decreasing hazard function.
To complete the model, we describe a prior distribution for the parameter vectors $(\phi,\sigma)$, $\delta$, $\beta$, $\xi$, $\pi$ and $\lambda$. These vectors are {\em a priori} independent.
We specify a multivariate log-Normal prior distribution for the vector $(\phi,\sigma)$ of OU parameters: $ \theta \equiv (\log \phi,\log \sigma) \sim \mathcal N(\bar \theta, \bar \Sigma). $ The transformation eliminates the need for parameter restrictions.
We induce a Gaussian prior for the vector $\delta$ of coefficients defining the diurnal pattern by specifying the following Gaussian prior on an invertible linear function of $\delta$:
where $v$ and $\Delta$ are the mean and first difference operators defined by \[ v \equiv \frac{1}{L}
, \quad \Delta \equiv
, \] $0_{L-1}$ is an $(L-1) \times 1$ vector of zeros, $I_{L-1}$ is the $(L-1) \times (L-1)$ identity matrix, $\bar{\delta}$ is a scalar location parameter and $\tau$ and $\bar{h}$ are scalar precision parameters. Thus, first differences $\delta_l - \delta_{l-1}$ are iid $\operatorname{\mathcal N}(0,\tau^{-1})$ and the arithmetic mean $v\delta$, is independent of the first differences, with $v\delta \sim \operatorname{\mathcal N}(\bar{\delta}, \bar{h}^{-1})$. This is the same as the random walk prior proposed by LangBrez04, except that they use $v = (1, 0, \ldots, 0)$. Both priors favour smoothness and are agnostic with respect to the signs of derivatives; since the derivatives of B-splines are linear combinations of the first differences $\delta_l - \delta_{l-1}$, we can interpret $\tau$ as a smoothing parameter for the diurnal pattern. Higher values of $\tau$ imply more shrinkage towards a flat diurnal pattern. Unlike the prior proposed by LangBrez04, where the prior variances of the $\delta_l$ increase with $l$, our prior has the symmetry property that $(\delta_1,\ldots,\delta_L)$ and $(\delta_L,\ldots,\delta_1)$ have the same distribution, and prior uncertainty about trading intensity is the same at the opening and at the closing of the market.
Equation (ref) induces the prior $\delta \sim N(\bar{\mu}, \bar{H})$, where $\bar{\mu} = \bar{\delta} L v$ and $\bar{H} = \bar{h} vv' + \tau \Delta' \Delta$. The prior covector, which we use in Section (ref) below, is $\bar{c} \equiv \bar{H}\bar{\mu} = \bar{h} v$. We will estimate $\tau$, and specify its prior as the following scaled chi-square: $\bar s \tau \sim \chi^2(\bar \nu)$.
For the beta mixture weights indexing the normalized duration density, we specify $\beta \sim \operatorname{Dirichlet}(\bar M \bar \beta)$, where $\bar \beta=(\bar \beta_1,\dots,\bar \beta_J) > 0$, $\sum \bar \beta_j=1$, and $\bar M>0$. The prior mean of $\beta$ is $\bar \beta$; $\bar M$ is a concentration parameter. When $\bar \beta = (1/J,\dots,1/J)$, the prior mean of $\beta$ corresponds to $g(\cdot)$ being uniform on $[0,1]$ and therefore $p_\epsilon(\cdot)$ being an exponential density. In this sense, $\bar \beta = (1/J,\dots,1/J)$ centres the prior distribution for $p_\epsilon(\cdot)$ around the exponential density, with $\bar M$ controlling the amount of shrinkage towards it; higher values of $\bar M$ imply more shrinkage.
The diagonal elements of the Markov transition matrix $\xi$ of the latent indicator $s_{d,i}$ are independent and beta distributed: $ \xi_{kk} \sim \operatorname{Beta}(\bar a_k,\bar b_k)$, $k=0,1$. The hazard parameters of the cluster duration density are independent gammas: $\lambda_k \sim \operatorname{Gamma}(\bar a_{\lambda_k}, \bar b_{\lambda_k})$, $k=1,2$, and the weight $\pi$ of the first component is beta, with $\pi \sim \operatorname{Beta}(\bar a_\pi, \bar b_\pi)$.
We conclude the exposition of the model by giving the joint density of all parameters, latent variables and observations, making explicit all conditional independence relationships. We refer to the model as the all-duration flexible SCD model (FSCD). Let $s$, $x$ and $y$ be the flat vectors of all indicators, states and durations, respectively. Then the joint density is
Here, the densities for the initial values $s_{d,1}$ and $x_{d,1}$ on each day $d$ are understood to be $ p(s_{d,1} \,\vert\, s_{d,0}, \xi) \equiv p(s_{d,1} \,\vert\, \xi) $ and $ p(x_{d,1} \,\vert\, x_{d,-1}, t_{d,0}, t_{d,-1}, \phi, \sigma, \delta) \equiv p(x_{d,1} \,\vert\, t_{d,0}, \sigma, \delta). $
We also define the regular-duration model (R-FSCD), the special case of our model where all durations are regular, which is suitable for use with data where trades are aggregated into clusters. The joint density of the regular-duration model is obtained by removing factors related to $\xi$, $\lambda$, $\pi$ and $s$ in the above expressions.
Here we describe posterior simulation methods for Bayesian inference in our flexible SCD model. We do this for the all-duration model; methods for the regular-duration model require only straightforward modifications. We end this section by presenting some recommended adjustments to the model and simulation methods, when transaction times are truncated to the second.
We use Markov chain Monte Carlo (MCMC) to sample the joint posterior distribution of parameters, latent indicators and state variables. For the all-duration model, there are six Gibbs blocks, updating $(x,\phi,\sigma)$, $(x, \beta)$, $(\delta,\tau)$, $s$, $\xi$ and $(\lambda,\pi)$. For the regular-duration model, we only require the first three blocks. Note that the state sequence $x$ is updated twice, once jointly with $\phi$ and $\sigma$ and once jointly with $\beta$. We now describe each of the Gibbs blocks in turn.
It is widely known that when there is strong posterior dependence between state variables (here, $x$) and the parameters of their dynamics (here, $\phi$ and $\sigma$), updating both in a single block improves numerical efficiency. However, in non-linear non-Gaussian state-space models, it is difficult to draw the state sequence as a single block, even when conditioning on the parameters of the state dynamics. Several methods have been proposed, and many of them have been applied to draw latent volatilities in stochastic volatility models KimShepChib98, ChibNardShep02, RichZhan07. Here, an additional difficulty is the non-parametric nature of the measurement distribution, which appears to rule out methods based on auxiliary mixture models. In this paper, we use the HESSIAN method of McCa12, a procedure to draw the state sequence as a single block that does not require data augmentation. McCa12 shows how to adapt it to jointly draw the state sequence and associated parameters, which we do here.
Recall that we specify a prior for $\theta \equiv (\log \phi, \log \sigma)$. Let $\theta_{-}$ be the vector of all parameters except $\theta$. Our joint proposal $(x^*, \theta^*)$ consists of a proposal $\theta^*$ drawn from a proposal density $q(\theta \,\vert\, y)$ followed by a conditional proposal $x^*\,\vert\,\theta^*$ drawn from a proposal density $q(x \,\vert\, \theta, \theta_{-}, s, y)$. We accept the pair $(x^*, \theta^*)$ with probability \[ \min \left \{ 1,\, \frac { p(y \,\vert\, \theta^*, \theta_{-} ,s, x^*) p(x^* \,\vert\, \theta^*, \theta_{-}) p(\theta^*)} { p(y \,\vert\, \theta, \theta_{-}, s, x) p(x \,\vert\, \theta, \theta_{-}) p(\theta) } \times \frac { q(\theta \,\vert\, y) q(x \,\vert\, \theta, \theta_{-},s, y) } { q(\theta^* \,\vert\, y) q(x^* \,\vert\, \theta^*, \theta_{-},s, y) } \right \}. \] The proposal density $q(\theta \,\vert\, y)$ is a multivariate Student's $t$ density $t_\nu(\hat \theta, \hat \Sigma_\theta)$ with $\nu=15$ degrees of freedom, where $\hat \theta$ and $\hat \Sigma_\theta$ are approximations of the mean and covariance of the marginal posterior distribution of $\theta$, computed during a burn-in period. During the burn-in period, we use the adaptive random walk Metropolis approach described in Viho12. At the end of the burn-in period, we compute $\hat \theta$ and $\hat \Sigma_\theta$ as the sample mean and covariance for the second half of the burn-in period.
We draw $x \,\vert\, \theta^*,\theta_{-},s, y$ using the HESSIAN method. In general terms, and suppressing notation for any parameters there may be, this gives a close approximation $q(x \,\vert\, y)$ of the conditional density $p(x \,\vert\, y)$ of the state sequence $x$ given the observed sequence $y$, for state-space models with univariate states in which $p(y \,\vert\, x) = \prod_{i=1}^n p(y_i \,\vert\, x_i)$ and $x \sim \mathcal N(\bar \Omega^{-1}\bar c, \bar \Omega^{-1})$, with $\bar \Omega$ tridiagonal. Tridiagonality of $\bar \Omega$ corresponds to $x$ being Markov but not necessarily homogenous. The method is generic, as the only model-specific code required consists of a routine to evaluate $\log p(y_i \,\vert\, x_i)$, and its first five derivatives with respect to $x_i$, at a given point. For our stochastic duration model, we compute exact values of these derivatives without deriving analytic expressions for them; instead, we exploit automatic routines to combine evaluations of derivatives of primitive functions using Fa\'a di Bruno's rule, which is much easier. Details are provided in Appendix (ref).
In computational experiments not reported here, we discovered that numerical efficiency is improved by drawing $(x, \beta)$ as a block, compared to drawing $\beta$ alone, even accounting for the additional computation required to draw $x$ again in a second block. Recall that $x$ and $\beta$ govern the scale and shape, respectively, of the conditional distribution of regular durations.
We update $(x, \beta)$ in almost the same way as we update $(x, \theta)$. For the purposes of drawing proposals, we use the logistic transformation $\vartheta(\beta) = (\log(\beta_1/\beta_J),\dots,\log(\beta_{J-1}/\beta_J))$, which maps the $J-1$ dimensional simplex to $\mathbb{R}^{J-1}$. The absence of positivity and adding-up constraints is convenient, as is the fact that the posterior distribution of $\vartheta$ is more nearly Gaussian than that of $\beta$, especially when some elements of $\beta$ have high posterior density close to zero.
Let $\beta_{-}$ be the vector of all parameters except $\beta$. The joint proposal $(x^*, \beta^*)$ consists of a proposal $\vartheta^*$ drawn from a proposal density $q_\vartheta(\vartheta \,\vert\, y)$, which is transformed back to $\beta^*$ using the inverse transformation $\beta(\vartheta) = (1+\sum_{j=1}^{J-1} \exp(\vartheta_j))^{-1} (\exp(\vartheta_1), \ldots, \exp(\vartheta_{J-1}), 1)$, followed by a conditional proposal $x^*|\beta^*$ drawn from a proposal density $q(x \,\vert\, \beta^*, \beta_{-}, s, y)$. We accept $(x^*,\beta^*)$ with probability \[ \min \left \{ 1,\, \frac { p(y \,\vert\, \beta^*, \beta_{-} ,s, x^*) p(x^* \,\vert\, \beta^*, \beta_{-}) p(\beta^*)} { p(y \,\vert\, \beta, \beta_{-}, s, x) p(x \,\vert\, \beta, \beta_{-}) p(\beta) } \times \frac { q_\beta(\beta \,\vert\, y) q(x \,\vert\, \beta, \beta_{-}, s, y) } { q_\beta(\beta^* \,\vert\, y) q(x^* \,\vert\, \beta^*, \beta_{-}, s, y) } \right \}, \] where the implied proposal density $q_\beta(\beta \,\vert\, y)$ for $\beta$ is given by $q_\beta(\beta \,\vert\, y) = q_\vartheta(\vartheta(\beta) \,\vert\, y) / \prod_{i=1}^J \beta_j$ and the proposal density $q_\vartheta(\vartheta \,\vert\, y)$ for the transformed parameter $\vartheta$ is a $(J-1)$-variate Student's $t$ with $\nu=15$ degrees of freedom, and mean and variance parameters determined during a burn-in period in exactly the same way as the proposal density $q(\theta \,\vert\, y)$, described previously.
Since the $\beta_j$ can be interpreted as mixture weights, it might seem that data augmentation might be useful here. However, $\beta$ does not only give the mixture weights defining the non-parametric distribution $G(\cdot)$. It also determines the hazard parameter of the exponential distribution $F(\cdot)$, through the unit-mean normalization of the conditional duration density for regular durations.
We update $(\delta,\tau)$ using two sub-blocks. We first draw $\tau$ from its conditional posterior distribution: $\dbar s \tau \,\vert\, \delta \sim \chi^2(\dbar \nu),$ where $\dbar s = \bar s + \delta' \Delta' \Delta \delta$ and $ \dbar \nu = \bar \nu + L - 1$. We then draw $\delta$ from its conditional posterior distribution: $\delta \sim N(\dbar \mu, \dbar H)$, where the posterior precision is $\dbar H = \bar H + W'W$, the posterior mean is $\dbar \mu = \dbar H^{-1} \dbar c$ and the posterior covector is $\dbar c = \bar c + W'\tilde{v}$. The vector $\tilde{v}$ and matrix $W$ come from writing the state equation (ref) as $\tilde{v} \sim \mathcal N(W\delta,I_N)$ where $v$ and $W$ are organized in blocks \[
, \quad
, \] with, for $d = 1,\ldots,D$ and $l = 1,\ldots,L$, \[ \tilde{v}_d =
\] and \[ W_{dl} =
. \]
Latent indicators are updated via a single-move, or one-at-a-time, sampler. As Fruh06 points out, single-move sampling is faster than drawing states as a block, as no filtering is required. Since there is little posterior autocorrelation (most probabilities are close to zero or one, regardless of past and future values) there is little loss of numerical precision. The relative conditional probabilities of drawing $s_{d,i}=0$ and $s_{d,i}=1$, given $s_{-(d,i)}$, the rest of the indicators, are given by
for $d=1,\dots,D$ and $i=1,\ldots,n_d$. Here, the first and the last conditional transition probabilities are understood to be $\xi_{j,s_{d,0}} \equiv \xi_j$ and $\xi_{j,s_{d,n_d+1}} \equiv 1$.
Given the latent indicators, we update the transition probabilities using a Metropolis-Hasting step. The prior for $(\xi_{00},\xi_{11})$ ($\xi_{00}$ and $\xi_{11}$ are independent betas) is nearly conditionally conjugate, but since the first indicator of each day comes from the marginal distribution, not exactly so. The target density can be written as
where $N_{lk} = \sum_{d,i} \mathbf 1\{s_{d,i} =l, s_{d,i+1} = k\}$ is the number of transitions from $l$ to $k$ over all $D$ days. We draw a proposal using two independent beta distributions, $\xi_{00}^* \,\vert\, s \sim \operatorname{Beta}(N_{00} + \bar a_0, N_{01} + \bar b_0)$ and $\xi_{11}^* \,\vert\, s \sim \operatorname{Beta}(N_{11} + \bar a_1, N_{10} + \bar b_1)$. This would be an exact draw from the conditional posterior if we were conditioning on the first indicator in each day. We correct for the approximation by accepting the proposal with probability \[ \min \left\{ 1, \prod_{d=1}^D \left(\frac{1-\xi^*_{11}}{1-\xi_{11}}\right)^{1-s_{d,1}} \left(\frac{1-\xi^*_{00}}{1-\xi_{00}}\right)^{s_{d,1}} \left(\frac{2-\xi_{00}-\xi_{11}}{2-\xi^*_{00}-\xi^*_{11}}\right) \right\}. \]
Recall that cluster durations have a mixture distribution governed by a vector $(\pi,1-\pi)$ of mixture weights and a vector $\lambda = (\lambda_1, \lambda_2)$ of exponential rate parameters. To update $(\lambda, \pi)$, we first draw component indicators $z_{d,i} \in \{1,2\}$ for each $(d,i)$ for which the duration $y_{d,i}$ is a cluster duration; that is, for which $s_{d,i} = 0$. The indicators $z_{d,i}$ are conditionally independent, with probability mass function given, up to a multiplicative factor, by
We then draw $(\lambda,\pi)$ from its conditional distribution given all the $z_{di}$ and $y_{di}$ where $s_{d,i}=0$. The priors for $\lambda_1$, $\lambda_2$ and $\pi$ are conditionally conjugate, and the conditional posterior distributions are $ \lambda_k \,\vert\, z, s, y \sim \operatorname{Gamma}(\bar a_{\lambda_k} + \tilde N_k, \bar b_{\lambda_k} + \tilde N_k \tilde y_k), $ $k=1,2$, and $ \pi \,\vert\, z, s, y \sim \operatorname{Beta}(\bar a_\pi + \tilde N_1, \bar b_\pi + \tilde N_2), $ where $\tilde{N_k} = \sum_{d,i} \mathbf 1\{z_{d,i} = k\}$ and $\tilde y_k = (1/\tilde{N_k})\sum_{d,i} y_{d,i} \mathbf1\{z_{d,i} = k\}$.
Discreteness of duration data arising from limited recording precision is usually ignored, and models mostly feature continuous distributions. For data recorded to millisecond precision or finer, this is innocuous, but for second precision, the truncation error can be a large fraction of many regular durations, leading to estimation bias and making continuous-time models unreliable for hypothesis testing GMMT05,ZhanLiuBai10,SchnKomlAhma10. For this reason, we recommend some minor adjustments to the model when using data truncated to the second. Discreteness is a property of the recording technology and not the underlying process, so we treat the data as a censored continuous-time process. BlasHolyToma20 recently proposed an alternative approach based on parametric discrete distributions, with and without zero inflation.
We first replace the conditional density function of regular durations in equation (ref) by a probability mass function obtained by integrating over the range of possible values,
For the probability of $y_{d,i}=0$, the integral is from zero to one.
We then replace the distribution of cluster durations by a Bernoulli distribution with parameter $\zeta$; the values of $0s$ and $1s$ have probabilities $\zeta$ and $(1-\zeta)$, respectively. We specify a conditionally conjugate beta prior: $\zeta \sim \operatorname{Beta}(\bar a_\zeta, \bar b_\zeta)$.
We require some straightforward adjustments to our posterior simulator. The block updating $(\lambda,\pi)$ is replaced by a block updating $\zeta$. The conditional posterior distribution for $\zeta$ is $\zeta \,\vert\, s, y \sim \operatorname{Beta}(\bar a_\zeta + N_{00}, N_{11} + \bar b_\zeta)$. The computation of the derivatives required by the HESSIAN method needs straightforward adjustment.
Here, we first report results from an artificial data experiment meant to test for the correctness of our posterior simulators. We then illustrate the use of the flexible SCD model with an empirical application, using transaction data for two equities traded on the Toronto Stock Exchange.
When posterior simulation methods are based on incorrect analysis or when there are coding errors in their implementation, we cannot rely on the results they produce. The tests for program correctness described here are similar to those described in Gewe04 and the title of this section comes from the title of that paper.
The idea of the exercise is to simulate a Markov chain whose stationary distribution is the {\em joint} distribution of parameters, latent state variables {\em and data}, making use of the same simulation methods that will be used later for posterior simulation, together with an additional Gibbs block to draw data from their conditional distribution given parameters and state variables. If the posterior simulation methods are correct in concept and implementation, the marginal distribution of the parameter vector, with respect to this stationary distribution, is identical to its (known) prior distribution. This is a strong condition with many easily testable implications.
We will need to simplify part of our model in order to proceed. The problem is that we record the values of the underlying continuous-time OU process $x(t)$ only at the times that trades occur. Redrawing duration data changes the trading times, which requires conditioning on the entire path of the continuous-time process, which is impractical. For the GIR simulations only, we modify the latent state process described in (ref), replacing the OU process with a homogenous autoregressive process where the transition distribution from $x_i$ to $x_{i+1}$ depends only on $x_i$ and not the duration $y_i$. We parameterize the process in a way that resembles the sampled OU process, giving
The additional Gibbs block updating the data $y$ from its conditional distribution given parameters and latent variables is described in Appendix (ref). Sampling from the posterior distribution requires some minor modifications. Since we draw a new sample of artificial data at each iteration, the adaptive schemes implemented to approximate the mean and the covariance of $p(\theta \,\vert\, y)$ and $p(\beta \,\vert\, y)$ during the burn-in period do not work well. Instead, we use independence Metropolis-Hastings updates, where the proposal distribution has the same mean and covariance as the parameter vector in question.\footnote{ We approximate the mean and covariance matrix of $p(\vartheta)$ for $\vartheta(\beta) = (\log(\beta_1/\beta_J),\dots,\log(\beta_{J-1}/\beta_J))$by
where $\psi$ is the digamma function and $\psi'$ the trigamma functions. Those relations are the solution of the approximation of a Dirichlet distribution by a logistic normal distribution minimizing their Kullback-Leibler divergence. } We set the number of components of the normalized distribution to $J=3$ and use a B-spline function defined on two knots, $t_\mathrm{open}$ and $t_\mathrm{close}$, giving a diurnal pattern that is an expansion with $L=4$ cubic polynomials. To avoid trades in simulations occurring after $t_\mathrm{close}$, where the diurnal pattern is undefined, we choose the sample size $n$ (the number of durations, not the size of the simulation sample) and prior distributions such that the probability that the last transaction of the day occurs after $t_\mathrm{close}$ is extremely small. We fix $D=1$ and $t_0=t_\mathrm{open}$, choose a sample size of $n=50$ observations, and set the length of the trading session to $600$ seconds. Values of the prior hyper-parameters are shown in Table (ref). The tighter prior distribution, and much smaller number of observations, compared with our empirical application, ensure high numerical precision with a moderate amount of computation.
We verify the correctness of both the continuous-time model and the discrete-time version that accommodates the truncation of trading times to the second. We generate a sample of size $10^6$ in each case. Table (ref) shows the results, for the continuous-time model on the left and the discrete-time model on the right. In each case, the first column gives, for selected parameters $\vartheta$, the difference between the prior mean $E[\vartheta]$ and the simulation sample mean $\bar\vartheta$; the second, the numerical standard error (i.e.\ the simulation standard deviation quantifying error in finite simulation) of the sample mean; and the third, the $t$-statistic for the test of the hypothesis that the mean of the parameter, with respect to the stationary distribution of the Markov chain, equals the (known) prior mean. Numerical standard errors are computed using the overlapping batch means method FlegJone10. Rather than reporting results for all elements of $\delta$, we report results for their mean $v'\delta$ and the smoothing parameter $\tau$. Sample means are close to the true prior means, relative to the numerical standard error. Each of the hypotheses is a necessary condition for the correctness of our simulation methods. The results fail to cast doubt on this correctness: one hypothesis (out of 22) is rejected at the 5% level and none at the 1% level.
We demonstrate our FSCD model and posterior simulation methods using the transaction data described in Section (ref). The data, from March 2014, are for two equities traded on the Toronto Stock Exchange: the Royal Bank of Canada (RY) and the Potash Corporation (POT). We analyze the full samples directly using the all-duration model. To compare our probabilistic approach to classifying durations as regular or cluster durations with deterministic aggregation rules, we also construct subsamples of the durations classified as regular by the GW aggregation rule. We will refer to these subsamples as GW-filtered subsamples, and analyze them using the regular-duration model (R-FSCD). We remind the reader that the same-second aggregation rule of common practice is cruder than the GW aggregation rule; we would expect distortions arising from the same-second aggregation rule to be even more serious than those reported here.
For both specifications (all-duration and regular-duration) we report results for four models, each with a fixed value of $J$, the number of terms in the normalized conditional density; those values are $J=2,3,4$ and $5$. Rather than estimate $J$, which is difficult in models with two different state variables, we present results for each value of $J$ and compare them. In all models, the diurnal pattern is specified as a B-spline function defined on knots set on each half-hour, giving an expansion with $L=16$ piecewise polynomials.\footnote{ Trade durations are often standardized using a cubic spline specification with knots set at each hour, with extra knots in the first and last half hour to better capture rapid changes of trading intensity at the beginning and end of the trading day. } Values of the prior hyper-parameters are shown in the third column of Table (ref). We select a fairly diffuse prior distribution for the log-transformed parameters of the latent intensity state $x_{d,i}$, the mean of the coefficients of the B-spline function $m(t)$ describing diurnal patterns and the transition probabilities of the latent indicator $s_{d,i}$. The values of $(\bar s, \bar \nu)$ give a more diffuse prior distribution for the smoothing parameter $\tau$ than the one suggested by LangBrez04. For reasons given in the introduction, we center the prior distribution for the flexible conditional duration density around the exponential distribution. We set the concentration parameter for moderate shrinkage towards the exponential distribution: most of posterior variances of the $\beta_j$ coefficients are less than a quarter of the prior variances. The prior distribution for $\zeta$ favours cluster durations of $0s$ over $1s$, without being too informative.
Table (ref) and (ref) show the results for the RY series and the POT series, respectively. For each parameter, we report the posterior mean and standard deviation, and the relative numerical efficiency (RNE) for the posterior mean. Defined in Gewe89, the relative numerical efficiency is a variance ratio that quantifies the numerical precision of the sample mean of a ergodic process, relative to that of a (hypothetical) iid sample. RNE times sample size gives the size of an iid sample with the same numerical standard error. Numerical standard errors are computed using the overlapping batch mean method FlegJone10. The posterior samples consist of 50,000 retained draws recorded after a burn-in period of 15,000 draws.
It is well known that it is difficult to sample efficiently the persistence and variance parameters of latent states in non-Gaussian state-space models. We obtain a numerical efficiency for $\phi$ and $\sigma$ that is considerably higher than that reported for the analogous parameters using the block sampling method in StriForbMart06.\footnote{ The latent intensity state in their analysis follows a Gaussian AR(1) process with fixed autocorrelation and fixed innovation variance rather than an OU process sampled at irregular intervals. } The lowest numerical efficiency we obtain for these parameters is more than ten times the highest numerical efficiency they obtained for analogous parameters, using parametric duration distributions. Posterior sample means of the mixture weights have numerical efficiency ranging from 0.03 to 0.47 for the all-duration models and from 0.31 to 0.66 for the regular-duration models. Although numerical efficiencies are lower for the all-duration models, especially for specifications with a more flexible conditional duration density, they are quite good as far as methods for non-linear non-Gaussian state-space models go.
Table (ref) reports posterior quantiles and moments of the half-life $t_{1/2}$, measured in seconds, of the OU process $x_d(t)$: the length of time it takes for $x_d(t)$ and $x_d(t+t_{1/2})$ to have a correlation of 1/2 between them. This quantity is more easily interpretable than the mean reversion parameter $\phi$. Persistence of the latent intensity process is fairly high, and more so for RY. There is a fair degree of posterior uncertainty about $t_{1/2}$, and the posterior distribution is somewhat sensitive to how regular durations are classified and to the number of terms, $J$, in the normalized duration density. The marginal standard deviation $\sigma$ is estimated more precisely (see Tables (ref) and (ref)). Its distribution is somewhat sensitive to how regular durations are classified, but very little to $J$.
Figure (ref) shows the posterior mean (dashed line) and 25 posterior draws of the diurnal pattern $m(t)$, for the full sample and the all-duration model with $J=5$ terms. The left panel is for the RY series and the right panel for POT. The diurnal pattern for RY is flatter, indicating less predictable variation in trading intensity. In both cases, we obtain the usual inverted U-shaped diurnal pattern found in most studies, with more trading intensity near the opening and closing times. The posterior variation in the diurnal patterns is fairly small compared to the variation in average intensity through the day. The posterior distribution for $\delta$ is not very sensitive to $J$ or to the choice between deterministic and probabilistic classification. For this reason, we do not show illustrations similar to Figure (ref) for other specifications.
For all-duration models, we record, at each posterior draw, whether each duration recorded as $0s$ or $1s$ is a regular or a cluster duration; classifications vary from draw to draw. We summarize this in two ways: first, we compute, for each of these durations, the probability that it is regular and illustrate the variation of this probability over durations; second, we describe the posterior distribution of the number of these durations that are regular. Figure (ref) shows histograms of classification probabilities for durations recorded as $0s$ and as $1s$ under the all-duration model with $J=5$, for both series. The horizontal axes give the posterior probability that a duration is regular, and the vertical height of the bar at a given histogram bin gives the proportion of durations recorded as $0s$ or $1s$ whose posterior probability of being regular is within that bin. Every duration recorded as $0s$ has a posterior probability less than $0.3$ of being a regular duration; for a large majority, it is less than $0.05$. In contrast, most durations recorded as $1s$ have a posterior probability of more than $0.5$ of being regular. Histograms for other values of $J$ are similar and we do not report them.
While durations recorded as $0s$ are each quite unlikely to be regular, they are very numerous, and so the probability that many of them are regular is nonetheless very high. Table (ref) shows posterior quantiles and moments of the number of durations recorded as $0s$ or $1s$ classified as regular. Results for durations recorded as $0s$ are reported in the first four rows; those for durations recorded as $1s$, in the last four rows. For comparison, the GW rule classifies 386 durations recorded as $0s$ as regular for the RY series and 187 for the POT series. The GW rule does not apply to durations recorded as $1s$ and so all 2346 RY and all 1252 POT durations recorded as $1s$ are considered regular. Our probabilistic approach classifies as regular many more durations recorded as $0s$ and slightly fewer durations recorded as $1s$. For both series, the posterior standard deviation of their number increases with the number of terms $J$ in the normalized duration density; more terms allows more flexibility in the shape of the normalized conditional density near zero.
To resume: for each duration recorded as $0s$, classifying it as a cluster duration is the right choice, under symmetric loss. Collectively, these zero/one decisions lead one to severely underestimate the number of these durations that are regular. With our probabilistic approach, we get a very good idea about how many of these durations are regular, although of course we do not magically discover which ones. This approach is exactly how we avoid artifacts arising from the spurious aggregation of unrelated, but nearly simultaneous, trades.
Figure (ref) show normalized conditional density functions at the posterior mean of the $\beta$ coefficients. Upper panels are for the all-duration models and the full sample; lower panels for the regular-duration models and the GW-filtered subsamples; RY on the left and POT on the right. We see that the conditional density at the posterior mean of $\beta$ varies much more with $J$ at shorter durations than it does at longer durations. For regular-duration models, the density at zero decreases appreciably in $J$, especially for the RY series. With extra flexibility, the density increasingly fits the spurious scarcity of durations recorded as $0s$ in the GW-filtered sample; the resulting distortions are seen well away from zero. For the all-duration models (upper panels), there is less variation of the density at zero, despite the uncertainty about the number of durations recorded as $0s$ that are regular, and little difference between the two series. For the POT series, the mean densities obtained using $3$, $4$ and $5$ terms are barely distinguishable. For the all-duration models, the density at zero increases, rather than decreases in $J$, but the variation is not very pronounced, even near zero, where uncertainty about the number of regular durations recorded as $0s$ or $1s$ comes into play. While we are not able to estimate $J$, results do not change much between $J=4$ and $J=5$, and we suggest that the $J=5$ specification is able to capture well the normalized conditional density. Figure (ref) illustrates the posterior distribution of the hazard function for the normalized conditional duration density, for the full samples and the all-duration model with $J=5$. The dashed line shows the hazard function at the posterior mean of $\beta$. Solid lines give 25 posterior draws. In both cases, the hazard functions have a smoothly decreasing shape with high posterior probability. Although posterior precision is fairly high, there is still a non-negligible posterior probability that the hazard function is non-monotonic, at least for the POT series. The smooth variation in the hazard function, together with the possibility of it being non-monotonic, would have been impossible to capture with commonly used parametric conditional duration densities.
Models in the literature are designed to capture regular durations, those between unrelated trades. They are not intended, nor well suited, to capture the observed clustering of related trades. Common practice is to aggregate seemingly related trades into clusters and model only the “regular” durations between clusters. Even if trades could be classified as related or not without error, it is not clear that this would be desirable since it involves discarding information relevant to liquidity measurement and market microstructure. Furthermore, since it is not easy to tell related trades from unrelated trades that just happen to occur within the same second, errors of classification are inevitable. The most common aggregation rule, the same-second rule, amounts to calling all durations recorded as $0s$ cluster durations, and all others, regular durations. It is clear, however, that many of the durations recorded as $0s$ must be regular by happenstance---we just don't know which ones---and they are erroneously classified as cluster durations by the same-second rule. One consequence is to understate trade intensity and liquidity, especially at times of high intensity. Another is that the abrupt change in the number of durations between $0s$ and $1s$ makes it difficult to fit a normalized conditional duration distribution. The rule suggested by GramWell02 clearly mitigates the problem, but it does not eliminate it.
The solution we proposed is to make our model a mixture model, with a binary state variable indicating which durations are cluster durations and which are regular. Identification of the two states comes from the very tight distribution of cluster durations and, more subtly, the shrinkage of the normalized conditional duration density towards an exponential density, which varies slowly near zero. This probabilistic, rather than deterministic, classification, allows us to learn that any given pair of trades recorded in the same second are very probably related, at the same time as we learn that many such pairs are not in fact pairs of related trades.
Despite not learning which durations recorded as $0s$ are regular, we are able to estimate quite precisely the normalized conditional duration density. Its hazard function does not exhibit the large changes near zero that occur using the GramWell02 rule, which we claim is an artifact of the misclassification of many unrelated but nearly simultaneous trades as being related.
We introduced a flexible distribution for regular durations. Appealing to queueing theory, we argued that a good first-order model for durations between unrelated trades is an exponential distribution. We introduce a normalized conditional distribution for regular durations that is flexible, and also expressible as a perturbation of an exponential distribution. This allows us to shrink towards the exponential distribution.
Due in part to efficient draws of the latent trade intensity state, and despite the flexible distribution, numerical efficiency of posterior simulation is considerably better than that of previous studies where duration distributions are parametric.
In the empirical application, we found that the conditional hazard function for regular durations varies much less than what is found in many studies. We attribute this to better, probabilistic, classification of trades as related or not, and using flexible duration distributions instead of parametric distributions whose hazard functions have implausible behaviour near zero.