Stochastic Conditional Duration Model with Intraday Seasonality and Limit Order Book Information

Abstract

It is a widely known fact that the intraday seasonality of trading intervals for financial transactions such as stocks is short at the beginning of business hours and long in the middle of the day. In this paper, we extend the stochastic conditional duration (SCD) model to capture the pattern of intraday trading intervals and propose a new Markov chain Monte Carlo method to estimate this intraday seasonality simultaneously. To efficiently generate the Monte Carlo sample, we used a hybrid of the Gibbs/Metropolis–Hastings (MH) sampling scheme and also applied generalized Gibbs sampling. In addition to capturing this intraday seasonality, this paper also considers limit order book information. Three-day tick data for three stocks obtained from Nikkei NEEDS are used for estimation, and model selection is performed on smooth parameters, Weibull distribution and Gamma distribution. The typical intraday regularity of frequent trading immediately after the start of trading is confirmed, and the spread of the limit order book information is also found to affect the trading time interval.

1. Introduction

Recently, intraday transaction data have come to attract more and more researcher attention with the increased availability of high frequency data, development in computer technology and the increasing market share of high-frequency trading (HFT). According to Hosaka (2014), in the Tokyo Stock Exchange (TSE), HFT accounted for about 55% of total orders in May 2013, compared with about 23% in September 2012. In such a situation, the importance of the analysis of intraday data has been increasing.

The first characteristics of these intraday transaction data are that they are irregularly spaced. Thus, the conventional models for regularly spaced time series, such as the autoregressive conditional heteroskedasticity (ARCH) model by Engle (1982), the generalized autoregressive conditional heteroskedasticity (GARCH) model by Bollerslev (1986) and the stochastic volatility (SV) model by Taylor (1982), cannot be applied to the transaction data. As an alternative to these models for fixed interval data, Engle and Russell (1998) propose the autoregressive conditional duration (ACD) model. The word “duration” means the interval between two transaction times. In the ACD model, this duration is treated as the product of the conditional duration, which depends on past durations and noise.

Bauwens and Veredas (2004) proposed the stochastic conditional duration (SCD) model and assumed this duration to be driven by a latent variable. The state equation of the basic SCD model is a stationary first order autoregressive process (AR(1)) on the logarithm of the latent variable, and the observation equation is a product of the latent variable and noise, which follows either a Gamma or a Weibull distribution. The similarity of the forms between the SCD model and the SV model motivates some researchers to use an MCMC method for estimation. Strickland et al. (2006) developed a Bayesian MCMC method which employs a hybrid Gibbs/Metropolis–Hastings (MH) sampling scheme for estimating the SCD model. In this method, as proposed in Shephard and Pitt (1997), the state variables are divided into several blocks and are drawn in each block. Watanabe and Omori (2004) suggested that the blocking scheme proposed by Shephard and Pitt (1997) yields an estimation bias and Omori and Watanabe (2008) proposed a new efficient method for asymmetric stochastic volatility models.

The second characteristic of intraday transaction data is that they have an intraday seasonality. Trading frequency is known to be higher just after the start and just before the end of the market (i.e., shorter trading time intervals), which is called intraday seasonality. Engle and Russell (1998) estimated the ACD model using “diurnal adjusted data” in which diurnal seasonality was removed in advance by cubic spline. The log-ACD model proposed by Bauwens and Giot (2000) also estimates in two steps using cubic spline. Veredas et al. (2002) estimated the representation of intraday seasonality with a quartic kernel simultaneously in the log-ACD model. Brownlees and Vannucci (2013) proposed the mixed ACD model that incorporates the ACD model, intraday seasonality and daily random effect. The parameters of the mixed ACD model are estimated jointly by MCMC, and intraday seasonality is represented by a cubic B-Spline. In the SCD model, Bauwens and Veredas (2004) used a quartic kernel, and Feng et al. (2004) used a piecewise cubic spline to remove intraday seasonality in advance. Men et al. (2019) proposed the threshold stochastic conditional duration (TSCD) model. The model assumes that the AR(1) process switches between the two regimes, allowing us to consider the behavior of two types of traders with different behavioral criteria. Gingras and McCausland (2020) proposed the all-duration flexible SCD (FSCD) model. This is an extension of the SCD model that attempts to model transactions occurring at the same second with a mixture distribution and uses the B-Spline to capture intraday seasonality.

The main contributions of this study are twofold. First, we extend the simultaneous intraday seasonality estimation algorithm using MCMC within the SCD model. In the previous study by Strickland et al. (2006), the intraday seasonality was simultaneously estimated within the SCD model by MCMC, however, there are no parameters to be estimated by the MCMC loop because cubic smoothing spline is used. and the roughness parameter is estimated by generalized cross-sectional verification. In this study, intraday seasonality is simultaneously estimated via parameters by including the coefficients of the basis spline (B-Spline) in the MCMC loop. B-Spline is a piecewise polynomial function driven by control points. This property allows for the expression of smooth curves without the use of higher-order functions. In addition, since the curve is expressed as a linear combination, it can be estimated in the same framework as the estimation of the effect of external information, such as the limit order information. The second contribution is the treatment of the limit order information in the SCD model. Although the limit order information plays a very important role in capturing recent trends in financial markets, there have been only a few analyses using limit order information in the SCD model, such as Sugiura et al. (2015). In this paper, we estimate the effect of the limit order book information on the trading time interval using the variables of limit order price difference (spread), order volume, and the balance between buy and sell orders.

The structure of this paper is as follows. In Section 2, we introduce the SCD model. Section 3 describes a MCMC sampling method for this model. Section 4 provides estimates for three stocks listed on the TSE. Section 5 presents the conclusions.

2. Stochastic Conditional Duration Model

2.1. A Proposed Model

Defining ${\tau }_i$ as the time of the ith transaction, the ith duration $y_i$ can be written as $y_i={\tau }_{i+1}-{\tau }_i$. The SCD model is a non-Gaussian state space model and considers this duration $y_i$ as a product of a latent valuable ${\alpha }_i$ and a positive random variable ${\epsilon }_i$. The state equation is a stationary AR(1) process on the logarithm of the latent variable with noise ${\eta }_i$. Both ${\epsilon }_i$ and ${\eta }_i$ are mutually and serially independent in this model. In this paper, we consider the SCD model as

$$\begin{array}{cc}\hfill y_i& =exp\left(x_i^{\prime }\beta +{\alpha }_i\right){\epsilon }_i,{\epsilon }_i>0,i\in \{1,\dots ,N\},\hfill \end{array}$$

(1)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\alpha }_i& =\varphi {\alpha }_{i-1}+{\eta }_i,{\eta }_i\sim \mathcal{N}(0,{\sigma }^2),i\in \{2,\dots ,N\},\hfill \\ \hfill {\alpha }_1& \sim \mathcal{N}(0,{\sigma }^2/(1-{\varphi }^2)),\left|\varphi \right|<1,\hfill \end{array}\end{array}$$

(2)

where a $\left(1\times k\right)$ vector $x_i^{\prime }$ is the ith row of the basis spline (B-Spline) function $\left(N\times k\right)$ matrix X and a $\left(k\times 1\right)$ vector $\beta$ is the control points of B-Spline. The random variable ${\epsilon }_i$ is supposed to follow a positive-valued continuous distribution. In the literature,

$$\begin{array}{cc}\hfill \mathrm{Weibull}(\gamma ,1):& p\left({\epsilon }_i\right|\gamma )=\gamma {\epsilon }_i^{\gamma -1}exp(-{\epsilon }_i^{\gamma }),\hfill \\ \hfill \mathrm{Gamma}(\gamma ,\gamma ):& p\left({\epsilon }_i\right|\gamma )=\frac{{\gamma }^{\gamma }}{\mathsf{\Gamma }\left(\gamma \right)}{\epsilon }_i^{\gamma -1}exp(-\gamma {\epsilon }_i),\hfill \end{array}$$

are often used1 as the distribution of ${\epsilon }_i$. In our study, we will follow the precedent. Thus, the probability density function (pdf) of $y_i$ is

$$p\left(y_i\right|{\alpha }_i,\beta ,\gamma )=\left\{\begin{array}{c}{\displaystyle \gamma y_i^{\gamma -1}exp\left(-\gamma (x_i^{\prime }\beta +{\alpha }_i)-y_i^{\gamma }exp(-\gamma (x_i^{\prime }\beta +{\alpha }_i))\right),}\hfill \\ (\mathrm{Weibull}\mathrm{distribution}),\hfill \\ {\displaystyle \frac{{\gamma }^{\gamma }}{\mathsf{\Gamma }\left(\gamma \right)}{y}_i^{\gamma -1}exp\left(-\gamma \left(x_i^{\prime }\beta +{\alpha }_i+y_{i}exp(-x_i^{\prime }\beta -{\alpha }_i)\right)\right),}\hfill \\ (\mathrm{gamma}\mathrm{distribution}).\hfill \end{array}\right.$$

(3)

2.2. Basis Spline (B-Spline)

To imply the intraday seasonality by the model, we follow Eilers and Marx (1996), Lang and Brezger (2004), Brownlees and Vannucci (2013) and Gingras and McCausland (2020) and use the B-Spline in the observation equation. B-Spline is a piecewise polynomial function driven by control points. Owing to the properties, we can express a smooth curve without using the higher order function. In our model, the B-Spline is implemented as

$$x_i^{B^{\prime }}=(B_{1,n}\left({\tau }_i\right),\dots ,B_{k,n}\left({\tau }_i\right)),k=m-n-1$$

(4)

where m is the number of the knots $t_1\le t_2\le \cdots \le t_m$, and n is the order of B-Spline. $B_{j,n}\left({\tau }_i\right)$ is the basis function defined by means of the Cox–de Boor recursion formula by De Boor and De Boor (1978) as

$$\begin{array}{cc}\hfill B_{j,0}\left({\tau }_i\right)& =\left\{\begin{array}{c}1\mathrm{if}{t}_j\le {\tau }_i<t_{j+1}\hfill \\ 0\mathrm{otherwise},\hfill \end{array}\right.\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill B_{j,n}\left({\tau }_i\right)& =\frac{{\tau }_i-t_j}{t_{j+n}-t_j}{B}_{j,n-1}\left({\tau }_i\right)+\frac{t_{j+n+1}-{\tau }_i}{t_{j+n+1}-t_{j+1}}{B}_{j+1,n-1}\left({\tau }_i\right)\hfill \end{array}$$

(6)

2.3. Limit Order Book Information

By the same method, we added the order book information as

$$x_i^{\prime }=(x_i^{B^{\prime }},d_{i1},d_{i2},d_{i3})$$

(7)

where $d_1$ is the price spread, $d_2$ is the sum of the volume at the best bid and ask price, and $d_3$ is the bid-ask volume ratio, which means the value of the best bid volume divided by the best ask volume minus one for centering. This design allows us to simultaneously estimate the regression coefficients $\beta$ on the diurnal seasonality and the effect of external information in the SCD model.

3. Estimation Method

3.1. Joint Posterior Distribution

The SCD model (1)–(2) is a non-linear non-Gaussian state–space model where (1) is the measurement equation, (2) is the transition equation, and ${\alpha }_i$, $i\in \{1,\dots ,n\}$ are the state variables. Since (2) is a stationary AR(1) model, the joint probability distribution of $\alpha =[{\alpha }_1;\cdots ;{\alpha }_n]$ is $\mathcal{N}(0,{\sigma }^2{V}^{-1})$ where

$$V=\left[\begin{array}{ccccccc}1& -\varphi & 0& \dots & \dots & \dots & 0\\ -\varphi & 1+{\varphi }^2& -\varphi & \ddots & & & ⋮\\ 0& -\varphi & 1+{\varphi }^2& -\varphi & \ddots & & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & \ddots & -\varphi & 1+{\varphi }^2& -\varphi & 0\\ ⋮& & & \ddots & -\varphi & 1+{\varphi }^2& -\varphi \\ 0& \dots & \dots & \dots & 0& -\varphi & 1\end{array}\right],$$

(8)

is a tridiagonal matrix, and it is positive definite as long as $\left|\varphi \right|<1$. The prior of $\beta$, $\gamma$, $\varphi$ and ${\sigma }^2$ in our study are

$$\begin{array}{cc}\hfill & \beta \sim \mathcal{N}({\overline{\mu }}_{\beta },{\overline{A}}_{\beta }^{-1}),\gamma \sim \mathrm{Gamma}(a_{\gamma },b_{\gamma }),\hfill \\ \hfill & \frac{\varphi +1}{2}\sim \mathrm{Beta}(a_{\varphi },b_{\varphi }),{\sigma }^2\sim \mathrm{Inv}.\mathrm{Gamma}(a_{\sigma },b_{\sigma }),\hfill \end{array}$$

(9)

where

$$\begin{array}{cc}\hfill \mathcal{N}(\mu ,\sigma ):& p\left(x\right|\mu ,\sigma )=\frac{1}{\sqrt{2\pi {\sigma }^2}}exp\left\{-\frac{{(x-\mu )}^2}{2{\sigma }^2}\right\},\hfill \\ \hfill \mathrm{Gamma}(\alpha ,\beta ):& p\left(x\right|\alpha ,\beta )=\frac{{\beta }^{\alpha }}{\mathsf{\Gamma }\left(\alpha \right)}{x}_i^{\alpha -1}exp(-\beta x),x>0,\alpha ,\beta >0,\hfill \\ \hfill \mathrm{Beta}(\alpha ,\beta ):& p\left(x\right|\alpha ,\beta )=\frac{\mathsf{\Gamma }(\alpha +\beta )}{\mathsf{\Gamma }\left(\alpha \right)\mathsf{\Gamma }\left(\beta \right)}{x}^{\alpha -1}{(1-x)}^{\beta -1},0<x<1,\alpha ,\beta >0\hfill \\ \hfill \mathrm{Inv}.\mathrm{Gamma}(\alpha ,\beta ):& p\left(x\right|\alpha ,\beta )=\frac{{\beta }^{\alpha }}{\mathsf{\Gamma }\left(\alpha \right)}{x}_i^{-\alpha -1}exp(-\frac{\beta }{x}),x>0,\alpha ,\beta >0.\hfill \end{array}$$

Therefore, the joint posterior density of $(\alpha ,\beta ,\gamma ,{\sigma }^2)$ for the SCD model (1)–(2) is

$$p(\alpha ,\beta ,\gamma ,{\sigma }^2|y)\propto \prod _{i=1}^{n}p\left(y_i\right|{\alpha }_i,\beta ,\gamma )·p\left(\alpha \right|\varphi ,{\sigma }^2)·p\left(\beta \right)p\left(\gamma \right)p\left(\varphi \right)p\left({\sigma }^2\right),$$

(10)

where $p\left(y_i\right|{\alpha }_i,\beta ,\gamma )$ is the pdf of the Weibull or Gamma distribution in (3), $p\left(\alpha \right|\varphi ,{\sigma }^2)$ is the pdf of $\mathcal{N}(0,{\sigma }^2{V}^{-1})$, and $p\left(\beta \right),\dots ,p\left({\sigma }^2\right)$ are the pdf’s of the priors in (9). Since analytical evaluation of the joint posterior distribution (10) is impractical, we apply an MCMC method.

In the MCMC method, we generate a random sample of $(\alpha ,\beta ,\gamma ,\varphi ,{\sigma }^2)$, say $\left\{\right({\alpha }^{\left(r\right)},{\beta }^{\left(r\right)},$ ${\gamma }^{\left(r\right)},{\varphi }^{\left(r\right)},{\sigma }^{2\left(r\right)}{\left)\right\}}_{r=1}^m$, from the joint posterior distribution (10), and numerically evaluate the posterior statistics necessary for Bayesian inference with Monte Carlo integration. The outline of the sampling scheme is given as follows:

Step 1: 

Initialize $({\alpha }^{\left(0\right)},{\beta }^{\left(0\right)},{\gamma }^{\left(0\right)},{\varphi }^{\left(0\right)},{\sigma }^{2\left(0\right)})$ and set the counter $r=1$.

Step 2: 

Generate ${\alpha }^{\left(r\right)}$ from $p\left(\alpha \right|{\beta }^{(r-1)},{\gamma }^{(r-1)},{\varphi }^{(r-1)},{\sigma }^{2(r-1)},y)$.

Step 3: 

Generate ${\beta }^{\left(r\right)}$ from $p\left(\beta \right|{\alpha }^{\left(r\right)},{\gamma }^{(r-1)},{\varphi }^{(r-1)},{\sigma }^{2(r-1)},y)$.

Step 4: 

Generate ${\gamma }^{\left(r\right)}$ from $p\left(\gamma \right|{\alpha }^{\left(r\right)},{\beta }^{\left(r\right)},{\varphi }^{(r-1)},{\sigma }^{2(r-1)},y)$.

Step 5: 

Generate ${\varphi }^{\left(r\right)}$ from $p\left(\varphi \right|{\alpha }^{\left(r\right)},{\beta }^{\left(r\right)},{\gamma }^{\left(r\right)},{\sigma }^{2(r-1)},y)$.

Step 6: 

Generate ${\sigma }^{2\left(r\right)}$ from $p\left({\sigma }^2\right|{\alpha }^{\left(r\right)},{\beta }^{\left(r\right)},{\varphi }^{\left(r\right)},y)$.

Step 7: 

Let $r=r+1$, and go to Step 2 until the burn-in iterations are completed.

Step 8: 

Reset the counter $r=1$, and repeat Step 2–6 m times to obtain $\left\{\right({\alpha }^{\left(r\right)},{\beta }^{\left(r\right)},{\gamma }^{\left(r\right)},{\varphi }^{\left(r\right)},$ ${\sigma }^{2\left(r\right)}{\left)\right\}}_{r=1}^m$.

3.2. State Variables

The conditional posterior density of the state variables $\alpha$ is

$$p\left(\alpha \right|\beta ,\gamma ,\varphi ,{\sigma }^2,y)\propto \prod _{i=1}^{n}p\left(y_i\right|{\alpha }_i,\beta ,\gamma )·p\left(\alpha \right|\varphi ,{\sigma }^2).$$

(11)

First we consider the second-order Taylor expansion of the log likelihood

$$\begin{array}{cc}\hfill & \ell \left(\alpha \right)\equiv \sum _{i=1}^{n}logp\left(y_i\right|{\alpha }_i,\beta ,\gamma ),\hfill \\ \hfill & logp\left(y_i\right|{\alpha }_i,\beta ,\gamma )=\hfill \\ \hfill & \left\{\begin{array}{c}log\gamma -log{y}_i+\gamma (log{y}_i-x_i^{\prime }\beta -{\alpha }_i)-exp\left(\gamma (log{y}_i-x_i^{\prime }\beta -{\alpha }_i)\right),\hfill \\ (\mathrm{Weibull}\mathrm{distribution}),\hfill \\ \gamma log\gamma -log\mathsf{\Gamma }\left(\gamma \right)-log{y}_i+\gamma \left(log{y}_i-x_i^{\prime }\beta -{\alpha }_i-exp(log{y}_i-x_i^{\prime }\beta -{\alpha }_i)\right),\hfill \\ (\mathrm{Gamma}\mathrm{distribution}),\hfill \end{array}\right.\hfill \end{array}$$

(12)

in the neighborhood of ${\alpha }^{*}=[{\alpha }_1^{*};\cdots ;{\alpha }_n^{*}]\in {\mathbb{R}}^n$:

$$\ell \left(\alpha \right)\approx \ell \left({\alpha }^{*}\right)+g{\left({\alpha }^{*}\right)}^{\prime }(\alpha -{\alpha }^{*})-\frac{1}{2}{(\alpha -{\alpha }^{*})}^{\prime }Q\left({\alpha }^{*}\right)(\alpha -{\alpha }^{*}),$$

(13)

where

$$\begin{array}{cc}\hfill g\left({\alpha }^{*}\right)& \equiv {\nabla }_{\alpha }\ell \left({\alpha }^{*}\right)=[g_1\left({\alpha }^{*}\right);\cdots ;g_n\left({\alpha }^{*}\right)],\hfill \\ \hfill Q\left({\alpha }^{*}\right)& \equiv -{\nabla }_{\alpha }{\nabla }_{\alpha }^{\prime }\ell \left({\alpha }^{*}\right)=\mathrm{diag}\{q_1\left({\alpha }^{*}\right),\dots ,q_n\left({\alpha }^{*}\right)\},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill g_i\left({\alpha }^{*}\right)& =\left\{\begin{array}{cc}-\gamma +\gamma exp\left(\gamma (log{y}_i-x_i^{\prime }\beta -{\alpha }_i^{*})\right),\hfill & (\mathrm{Weibull}\mathrm{distribution}),\hfill \\ -\gamma +\gamma exp(log{y}_i-x_i^{\prime }\beta -{\alpha }_i^{*}),\hfill & (\mathrm{Gamma}\mathrm{distribution}),\hfill \end{array}\right.\hfill \\ \hfill q_i\left({\alpha }^{*}\right)& =\left\{\begin{array}{cc}{\gamma }^{2}exp\left(\gamma (log{y}_i-x_i^{\prime }\beta -{\alpha }_i^{*})\right),\hfill & (\mathrm{Weibull}\mathrm{distribution}),\hfill \\ \gamma exp(log{y}_i-x_i^{\prime }\beta -{\alpha }_i^{*}),\hfill & (\mathrm{Gamma}\mathrm{distribution}),\hfill \end{array}\right.\hfill \\ \hfill & i\in \{1,\dots ,n\}.\hfill \end{array}$$

Note that $Q\left({\alpha }^{*}\right)$ is always positive definite.

Since the log prior density of $\alpha$ is

$$logp\left(\alpha \right)\equiv -\frac{n}{2}log\left(2\pi {\sigma }^2\right)+\frac{1}{2}log(1-{\varphi }^2)-\frac{1}{2{\sigma }^2}{\alpha }^{\prime }V\alpha ,$$

(14)

the conditional posterior density of $\alpha$ can be approximated by

$$\begin{array}{cc}\hfill & p\left(\alpha \right|\beta ,\gamma ,\varphi ,{\tau }^2,y)\hfill \\ \hfill & =Cexp\left[\ell \left(\alpha \right)+logp\left(\alpha \right)\right]\hfill \\ \hfill & \approx Cexp\left[\ell \left({\alpha }^{*}\right)+g{\left({\alpha }^{*}\right)}^{\prime }(\alpha -{\alpha }^{*})-\frac{1}{2}{(\alpha -{\alpha }^{*})}^{\prime }Q\left({\alpha }^{*}\right)(\alpha -{\alpha }^{*})+logp\left(\alpha \right)\right]\hfill \\ \hfill & =Cexp\left[\ell \left({\alpha }^{*}\right)-\frac{n}{2}log\left(2\pi {\tau }^2\right)+\frac{1}{2}log(1-{\varphi }^2)\right]\hfill \\ \hfill & \times exp\left[g{\left({\alpha }^{*}\right)}^{\prime }(\alpha -{\alpha }^{*})-\frac{1}{2}{(\alpha -{\alpha }^{*})}^{\prime }Q\left({\alpha }^{*}\right)(\alpha -{\alpha }^{*})-\frac{1}{2{\tau }^2}{\alpha }^{\prime }V\alpha \right],\hfill \end{array}$$

(15)

where C is the normalizing constant. By completing the square in (15), we have

$$\begin{array}{cc}\hfill & g{\left({\alpha }^{*}\right)}^{\prime }(\alpha -{\alpha }^{*})-\frac{1}{2}{(\alpha -{\alpha }^{*})}^{\prime }Q\left({\alpha }^{*}\right)(\alpha -{\alpha }^{*})-\frac{1}{2{\sigma }^2}{\alpha }^{\prime }V\alpha \hfill \\ \hfill & =-\frac{1}{2}{(\alpha -{\mu }_{\alpha }\left({\alpha }^{*}\right))}^{\prime }{\mathsf{\Sigma }}_{\alpha }{\left({\alpha }^{*}\right)}^{-1}(\alpha -{\mu }_{\alpha }\left({\alpha }^{*}\right))\hfill \\ \hfill & +\frac{1}{2}g{\left({\alpha }^{*}\right)}^{\prime }Q{\left({\alpha }^{*}\right)}^{-1}g\left({\alpha }^{*}\right)\hfill \\ \hfill & -\frac{1}{2}{(g\left({\alpha }^{*}\right)+Q\left({\alpha }^{*}\right){\alpha }^{*})}^{\prime }\left({\mathsf{\Sigma }}_{\alpha }\left({\alpha }^{*}\right)+Q{\left({\alpha }^{*}\right)}^{-1}\right)(g\left({\alpha }^{*}\right)+Q\left({\alpha }^{*}\right){\alpha }^{*}),\hfill \end{array}$$

(16)

where

$${\mathsf{\Sigma }}_{\alpha }\left({\alpha }^{*}\right)={\left(\frac{1}{{\tau }^2}V+Q\left({\alpha }^{*}\right)\right)}^{-1},{\mu }_{\alpha }\left({\alpha }^{*}\right)={\mathsf{\Sigma }}_{\alpha }\left({\alpha }^{*}\right)\left(g\left({\alpha }^{*}\right)+Q\left({\alpha }^{*}\right){\alpha }^{*}\right).$$

Therefore, the right-hand side of (15) is proportional to the pdf of the following normal distribution:

$$\alpha \sim \mathcal{N}\left({\mu }_{\alpha }\left({\alpha }^{*}\right),{\mathsf{\Sigma }}_{\alpha }\left({\alpha }^{*}\right)\right).$$

(17)

Hoping that the approximation (15) is sufficiently accurate, we apply the Metropolis–Hastings (MH) algorithm to generate $\alpha$ from the conditional posterior distribution (11) by using (17) as the proposal distribution.

In practice, however, we need to address two issues:

• The choice of ${\alpha }^{*}$ is crucial to make the approximation (15) workable.
• The acceptance rate of the MH algorithm tends to be low when $\alpha$ is a high-dimensional vector.

We address the former issue by using the mode of the conditional posterior density as ${\alpha }^{*}$. The search of the mode is performed by the following recursion:

Step 1: 

Initialize ${\alpha }^{*\left(0\right)}$, and set the counter $r=1$.

Step 2: 

Update ${\alpha }^{*\left(r\right)}$ by ${\alpha }^{*\left(r\right)}={\mu }_{\alpha }\left({\alpha }^{*(r-1)}\right)$.

Step 3: 

Let $r=r+1$ and go to Step 2 unless ${max}_{i=1,\dots ,n}|{\alpha }_i^{*\left(r\right)}-{\alpha }_i^{*(r-1)}|$ is less than the preset tolerance level.

It turns out that the above algorithm is equivalent to the Newton–Raphson method, and in our experience, it mostly attains convergence in a few iterations.

We address the latter issue by applying a so-called block sampler. In the block sampler, we randomly partition $\alpha$ into several sub-vectors (blocks), generate each block ${\alpha }_b$ from its conditional distribution given the rest of the blocks ${\alpha }_r$ and apply the MH algorithm to each generated block.

$$\begin{array}{c}\hfill \alpha =\left(\begin{array}{c}{\alpha }_b\\ {\alpha }_r\end{array}\right),{\alpha }^{*}=\left(\begin{array}{c}{\alpha }_b^{*}\\ {\alpha }_r^{*}\end{array}\right),{\mathsf{\Sigma }}_{\alpha }{\left({\alpha }^{*}\right)}^{-1}=\left(\begin{array}{cc}{\mathsf{\Omega }}_{bb}& {\mathsf{\Omega }}_{br}\\ {\mathsf{\Omega }}_{rb}& {\mathsf{\Omega }}_{rr}\end{array}\right)\end{array}$$

(18)

Using (17), the kernel of the conditional distribution $p({\alpha }_b\mid {\alpha }_r)$ is introduced by

$$\begin{array}{cc}\hfill & {(\alpha -{\alpha }^{*})}^{\prime }{\mathsf{\Sigma }}_{\alpha }{\left({\alpha }^{*}\right)}^{-1}(\alpha -{\alpha }^{*})\hfill \\ \hfill =& {\left(\begin{array}{c}{\alpha }_b-{\alpha }_b^{*}\\ {\alpha }_r-{\alpha }_r^{*}\end{array}\right)}^{\prime }\left(\begin{array}{cc}{\mathsf{\Omega }}_{bb}& {\mathsf{\Omega }}_{br}\\ {\mathsf{\Omega }}_{rb}& {\mathsf{\Omega }}_{rr}\end{array}\right)\left(\begin{array}{c}{\alpha }_b-{\alpha }_b^{*}\\ {\alpha }_r-{\alpha }_r^{*}\end{array}\right)\hfill \\ \hfill =& {({\alpha }_b-{\alpha }_b^{*})}^{\prime }{\mathsf{\Omega }}_{bb}({\alpha }_b-{\alpha }_b^{*})+2{({\alpha }_r-{\alpha }_r^{*})}^{\prime }{\mathsf{\Omega }}_{rb}({\alpha }_b-{\alpha }_b^{*})\hfill \\ \hfill & +{({\alpha }_r-{\alpha }_r^{*})}^{\prime }{\mathsf{\Omega }}_{rr}({\alpha }_r-{\alpha }_r^{*})\hfill \end{array}$$

(19)

By completing the square in (19), we have

$$\begin{array}{cc}\hfill & {({\alpha }_b-{\alpha }_b^{*})}^{\prime }{\mathsf{\Omega }}_{bb}({\alpha }_b-{\alpha }_b^{*})+2{({\alpha }_r-{\alpha }_r^{*})}^{\prime }{\mathsf{\Omega }}_{rb}({\alpha }_b-{\alpha }_b^{*})+{({\alpha }_r-{\alpha }_r^{*})}^{\prime }{\mathsf{\Omega }}_{rr}({\alpha }_r-{\alpha }_r^{*})\hfill \\ \hfill =& {({\alpha }_b-{\mu }_{{\alpha }_b}\left({\alpha }^{*}\right))}^{\prime }{\mathsf{\Sigma }}_{{\alpha }_b}{\left({\alpha }^{*}\right)}^{-1}({\alpha }_b-{\mu }_{{\alpha }_b}\left({\alpha }^{*}\right))\hfill \\ \hfill & -2{({\alpha }_r-{\alpha }_r^{*})}^{\prime }{\mathsf{\Omega }}_{rb}{\alpha }_b^{*}+{({\alpha }_r-{\alpha }_r^{*})}^{\prime }({\mathsf{\Omega }}_{rr}-{\mathsf{\Omega }}_{rb}{\mathsf{\Omega }}_{bb}^{-1}{\mathsf{\Omega }}_{br})({\alpha }_r-{\alpha }_r^{*})\hfill \end{array}$$

(20)

where

$${\mu }_{{\alpha }_b}\left({\alpha }^{*}\right)={\alpha }_b^{*}-{\mathsf{\Omega }}_{bb}^{-1}{\mathsf{\Omega }}_{br}({\alpha }_r-{\alpha }_r^{*}),{\mathsf{\Sigma }}_{{\alpha }_b}\left({\alpha }^{*}\right)={\mathsf{\Omega }}_{bb}^{-1}$$

Therefore, the right-hand side of (19) is the kernel of the pdf of the following normal distribution:

$${\alpha }_b\sim \mathcal{N}\left({\mu }_{{\alpha }_b}\left({\alpha }^{*}\right),{\mathsf{\Sigma }}_{{\alpha }_b}\left({\alpha }^{*}\right)\right).$$

(21)

3.3. Regression Coefficients

The simulation strategy for the regression coefficients $\beta$ is almost identical to the one for the state variables $\alpha$. The second-order Taylor expansion of the log likelihood (12) with respect to $\beta$ in the neighborhood of ${\beta }^{*}\in {\mathbb{R}}^m$ is

$$\ell \left(\beta \right)\approx \ell \left({\beta }^{*}\right)+g{\left({\beta }^{*}\right)}^{\prime }(\beta -{\beta }^{*})-\frac{1}{2}{(\beta -{\beta }^{*})}^{\prime }Q\left({\beta }^{*}\right)(\beta -{\beta }^{*}),$$

(22)

where

$$\begin{array}{cc}\hfill g\left({\beta }^{*}\right)& \equiv {\nabla }_{\beta }\ell \left({\beta }^{*}\right)=\sum _{i=1}^n{g}_i\left({\beta }^{*}\right)x_i,\hfill \\ \hfill Q\left({\beta }^{*}\right)& \equiv -{\nabla }_{\beta }{\nabla }_{\beta }^{\prime }\ell \left({\beta }^{*}\right)=\sum _{i=1}^n{q}_i\left({\beta }^{*}\right)x_i{x}_i^{\prime },\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill g_i\left({\beta }^{*}\right)& =\left\{\begin{array}{cc}-\gamma +\gamma exp(\gamma (log{y}_i-x_i^{\prime }{\beta }^{*}-{\alpha }_i),\hfill & (\mathrm{Weibull}\mathrm{distribution}),\hfill \\ -\gamma +\gamma exp(log{y}_i-x_i^{\prime }{\beta }^{*}-{\alpha }_i),\hfill & (\mathrm{Gamma}\mathrm{distribution}),\hfill \end{array}\right.\hfill \\ \hfill q_i\left({\beta }^{*}\right)& =\left\{\begin{array}{cc}{\gamma }^{2}exp\left(\gamma (log{y}_i-x_i^{\prime }{\beta }^{*}-{\alpha }_i)\right),\hfill & (\mathrm{Weibull}\mathrm{distribution}),\hfill \\ \gamma exp(log{y}_i-x_i^{\prime }{\beta }^{*}-{\alpha }_i),\hfill & (\mathrm{Gamma}\mathrm{distribution}),\hfill \end{array}\right.\hfill \\ \hfill & i\in \{1,\dots ,n\}.\hfill \end{array}$$

Note that $Q\left({\beta }^{*}\right)$ is always positive definite.

With the prior $\beta \sim \mathcal{N}({\overline{\mu }}_{\beta },{\overline{A}}_{\beta }^{-1})$, the conditional posterior density of $\beta$ can be approximated by

$$\begin{array}{cc}\hfill & p\left(\beta \right|\alpha ,\gamma ,\varphi ,{\tau }^2,y)\hfill \\ \hfill & =Cexp\left[\ell \left(\beta \right)+logp\left(\beta \right)\right]\hfill \\ \hfill & \approx Cexp\left[\ell \left({\beta }^{*}\right)-\frac{1}{2}log\left(2\pi \right)+\frac{1}{2}log\left|{\overline{A}}_{\beta }\right|\right]\hfill \\ \hfill & \times exp\left[g{\left({\beta }^{*}\right)}^{\prime }(\beta -{\beta }^{*})-\frac{1}{2}{(\beta -{\beta }^{*})}^{\prime }Q\left({\beta }^{*}\right)(\beta -{\beta }^{*})-\frac{1}{2}{(\beta -{\overline{\mu }}_{\beta })}^{\prime }{\overline{A}}_{\beta }(\beta -{\overline{\mu }}_{\beta })\right].\hfill \end{array}$$

(23)

By completing the square as in (16), the proposal distribution for the MH algorithm is derived as

$$\beta \sim \mathcal{N}\left({\mu }_{\beta }\left({\beta }^{*}\right),{\mathsf{\Sigma }}_{\beta }\left({\beta }^{*}\right)\right),$$

(24)

where

$${\mathsf{\Sigma }}_{\beta }\left({\beta }^{*}\right)={\left({\overline{A}}_{\beta }+Q\left({\beta }^{*}\right)\right)}^{-1},{\mu }_{\beta }\left({\beta }^{*}\right)={\mathsf{\Sigma }}_{\beta }\left({\beta }^{*}\right)\left({\overline{A}}_{\beta }{\overline{\mu }}_{\beta }+g\left({\beta }^{*}\right)+Q\left({\beta }^{*}\right){\beta }^{*}\right).$$

The search algorithm for ${\beta }^{*}$ is the same as ${\alpha }^{*}$. Since the dimension of $\beta$ is considerably smaller than $\alpha$, it is not necessary to apply the block sampler in our experience.

3.4. Shape Parameter

The sampling strategy for the shape parameter $\gamma$ is also the same as $\alpha$ and $\beta$. Since the prior of $\gamma$ is not Gaussian, instead of the log likelihood (12), we consider the second-order Taylor expansion of the log conditional posterior density of $\gamma$:

$$f\left(\gamma \right)\equiv \sum _{i=1}^{n}logp\left(y_i\right|{\alpha }_i,\beta ,\gamma )+logp\left(\gamma \right)+\mathrm{constant},$$

with respect to $\gamma$ in the neighborhood of ${\gamma }^{*}>0$, i.e.,

$$f\left(\gamma \right)\approx f\left({\gamma }^{*}\right)+g\left({\gamma }^{*}\right)(\gamma -{\gamma }^{*})-\frac{1}{2}q\left({\gamma }^{*}\right){(\gamma -{\gamma }^{*})}^2,$$

(25)

where

$$\begin{array}{cc}\hfill g\left({\gamma }^{*}\right)& \equiv {\nabla }_{\gamma }f\left({\gamma }^{*}\right)\hfill \\ \hfill & =\left\{\begin{array}{c}{\displaystyle \frac{n}{{\gamma }^{*}}+\sum _{i=1}^n\left(u_i-u_i{e}^{{\gamma }^{*}{u}_i}\right)+\frac{a_{\gamma }-1}{{\gamma }^{*}}-b_{\gamma },}\hfill \\ (\mathrm{Weibull}\mathrm{distribution}),\hfill \\ {\displaystyle n+nlog{\gamma }^{*}-n{\psi }^{\left(0\right)}\left({\gamma }^{*}\right)+\sum _{i=1}^n(u_i-e^{u_i})+\frac{a_{\gamma }-1}{{\gamma }^{*}}-b_{\gamma },}\hfill \\ (\mathrm{Gamma}\mathrm{distribution}),\hfill \end{array}\right.\hfill \\ \hfill q\left({\gamma }^{*}\right)& \equiv -{\nabla }_{\gamma }^{2}f\left({\gamma }^{*}\right)\hfill \\ \hfill & =\left\{\begin{array}{cc}{\displaystyle \frac{n}{{\gamma }^{*2}}+\sum _{i=1}^n{u}_i^2{e}^{{\gamma }^{*}{u}_i}+\frac{a_{\gamma }-1}{{\gamma }^{*2}},}\hfill & (\mathrm{Weibull}\mathrm{distribution}),\hfill \\ {\displaystyle -\frac{n}{{\gamma }^{*}}+n{\psi }^{\left(1\right)}\left({\gamma }^{*}\right)+\frac{a_{\gamma }-1}{{\gamma }^{*2}},}\hfill & (\mathrm{Gamma}\mathrm{distribution}),\hfill \end{array}\right.\hfill \\ \hfill u_i& \equiv log{y}_i-x_i^{\prime }\beta -{\alpha }_i,i\in \{1,\dots ,n\},\hfill \end{array}$$

and ${\psi }^{\left(s\right)}$ is the polygamma function of order s. $q\left({\gamma }^{*}\right)$ for the Weibull distribution is positive for any ${\gamma }^{*}>0$ if $n+a_{\gamma }>1$ which is satisfied in any applications. To prove $q\left({\gamma }^{*}\right)>0$ for the Gamma distribution, we use the following formula (Equation (6.1.38) in Abramowitz et al. (1988)):

$$log\mathsf{\Gamma }\left(\gamma \right)=\frac{1}{2}log\left(2\pi \right)+\left(\gamma -\frac{1}{2}\right)log\gamma -\gamma +\frac{\theta }{12\gamma },\gamma >0,0<\theta <1.$$

(26)

Taking the differential of both side of (26) twice, we have

$${\psi }^{\left(1\right)}\left(\gamma \right)={\nabla }_{\gamma }^{2}log\mathsf{\Gamma }\left(\gamma \right)=\frac{1}{\gamma }+\frac{1}{2{\gamma }^2}+\frac{\theta }{6{\gamma }^3}.$$

(27)

Replacing ${\psi }^{\left(1\right)}$ in $q\left({\gamma }^{*}\right)$ with (27), we have

$$q\left({\gamma }^{*}\right)=\frac{n+2a_{\gamma }-2}{2{\gamma }^{*2}}+\frac{n\theta }{6{\gamma }^{*3}},$$

(28)

which is positive if $n+2a_{\gamma }>2$. This condition is always satisfied in practice.

We directly apply the completing-the-square technique to (25) and obtain the proposal distribution for the MH algorithm which is derived as

$$\gamma \sim \mathcal{N}\left({\mu }_{\gamma }\left({\beta }^{*}\right),{\sigma }_{\gamma }^2\left({\gamma }^{*}\right)\right),$$

(29)

where

$${\sigma }_{\gamma }^2\left({\gamma }^{*}\right)=\frac{1}{q\left({\gamma }^{*}\right)},{\mu }_{\gamma }\left({\gamma }^{*}\right)={\gamma }^{*}+\frac{g\left({\gamma }^{*}\right)}{q\left({\gamma }^{*}\right)}.$$

If we use the mode of $f\left(\gamma \right)$ as ${\gamma }^{*}$, $g\left({\gamma }^{*}\right)=0$ always holds due to the global concavity of $f\left(\gamma \right)$. Thus ${\mu }_{\gamma }\left({\gamma }^{*}\right)$ is effectively identical to ${\gamma }^{*}$. The search algorithm for ${\gamma }^{*}$ is identical to those for ${\alpha }^{*}$ and ${\beta }^{*}$ except for the fact that ${\gamma }^{*}$ is a scalar. The implementation of the MH algorithm with the proposal distribution (29) is straightforward.

3.5. AR(1) Coefficient and Variance

Once the state variables $\alpha$ are generated, the conditional posterior density of $\varphi$ is given by

$$\begin{array}{cc}\hfill p\left(\varphi \right|y,\alpha ,{\sigma }^2)& \propto \sqrt{1-{\varphi }^2}exp\left[-\frac{(1-{\varphi }^2){\alpha }_1^2+{\sum }_{i=2}^n{({\alpha }_i-\varphi {\alpha }_{i-1})}^2}{2{\sigma }^2}\right]\hfill \\ \hfill & \times {(1+\varphi )}^{a_{\varphi }-1}{(1-\varphi )}^{b_{\varphi }-1}.\hfill \end{array}$$

(30)

By completing the square, we have

$$\begin{array}{cc}\hfill & (1-{\varphi }^2){\alpha }_1^2+\sum _{i=2}^n{({\alpha }_i-\varphi {\alpha }_{i-1})}^2\hfill \\ \hfill & =(1-{\varphi }^2){\alpha }_1^2+\sum _{i=2}^n{\alpha }_i^2-2\varphi \sum _{i=2}^n{\alpha }_i{\alpha }_{i-1}+{\varphi }^2\sum _{i=2}^n{\alpha }_{i-1}^2\hfill \\ \hfill & ={\alpha }_1^2+\sum _{i=2}^n{\alpha }_i^2-2\varphi \sum _{i=2}^n{\alpha }_i{\alpha }_{i-1}+{\varphi }^2\sum _{i=2}^{n-1}{\alpha }_i^2\hfill \\ \hfill & =\sum _{i=2}^{n-1}{\alpha }_i^2{\left(\varphi -\frac{{\sum }_{i=2}^n{\alpha }_i{\alpha }_{i-1}}{{\sum }_{i=2}^{n-1}{\alpha }_i^2}\right)}^2+{\alpha }_1^2-\frac{{\left({\sum }_{i=2}^n{\alpha }_i{\alpha }_{i-1}\right)}^2}{{\sum }_{i=2}^{n-1}{\alpha }_i^2}.\hfill \end{array}$$

With the above expression in mind, we use the following normal distribution:

$$\varphi \sim \mathcal{N}\left(\frac{{\sum }_{i=2}^n{\alpha }_i{\alpha }_{i-1}}{{\sum }_{i=2}^{n-1}{\alpha }_i^2},\frac{{\sigma }^2}{{\sum }_{i=2}^{n-1}{\alpha }_i^2}\right),$$

(31)

as the proposal distribution in the MH algorithm for $\varphi$.

For the variance ${\sigma }^2$, the exact expression of the conditional posterior distribution is available. With the inverse Gamma prior, it is derived as

$${\sigma }^2|y,\alpha ,\varphi \sim \mathrm{Inv}.\mathrm{Gamma}\left(a_{\sigma }+\frac{n}{2},b_{\sigma }+\frac{1}{2}{\alpha }^{\prime }V\alpha \right).$$

(32)

3.6. Generalized Gibbs Sampler

To facilitate convergence of the Markov chain sampling path, we utilize a generalized Gibbs sampler by Liu and Sabatti (2000). To do so, we first apply the following reparameterization:

$$\xi =\frac{1}{\gamma },\tau =\sqrt{{\sigma }^2}.$$

to the posterior distribution of $(\alpha ,\beta ,\gamma ,\varphi ,{\sigma }^2)$. Then, the pdf of the reparameterized posterior distribution is given by

$$f(\alpha ,\beta ,\xi ,\varphi ,\sigma |y)=p(\alpha ,\beta ,{\xi }^{-1},\varphi ,{\sigma }^2|y)\left|J\right|,$$

(33)

where $p(\alpha ,\beta ,\gamma ,\varphi ,{\sigma }^2|y)$ is the pdf of the original posterior distribution, and $\left|J\right|={\xi }^{-2}\times 2\tau$ is the Jacobian. Next we introduce the scale transformation: $(c\alpha ,c\beta ,c\xi ,c\sigma )$, $c>0$. In the generalized Gibbs sampler, we generate a random scale c from its conditional posterior distribution and apply the scale transformation to the parameters so that they would be updated more frequently and bounce more freely around in the parameter space. The conditional posterior density of c is proportional to

$${\left|c\right|}^{n+k+1}p(c\alpha ,c\beta ,c\xi ,\varphi ,c\sigma |y),$$

(34)

where $n+k+1$ is the number of transformed parameters minus 1. The outline of the generalized Gibbs sampler is summarized as follows.

Step 1: 

Obtain $({\alpha }^{\left(r\right)},{\beta }^{\left(r\right)},{\gamma }^{\left(r\right)},{\sigma }^{2\left(r\right)})$ at the r-th run of the sampling scheme.

Step 2: 

Generate c from (34) with the MH algorithm.

Step 3: 

Transform the parameters as $(c{\alpha }^{\left(r\right)},c{\beta }^{\left(r\right)},{\gamma }^{\left(r\right)}/c,c^2{\sigma }^{2\left(r\right)})$.

At the end of every cycle in the sampling scheme, we execute Step 1–3. Since the pdf of c is non-standard, we use the MH algorithm to generate c in Step 2 in the same manner as $(\alpha ,\beta ,\gamma )$, and we execute Step 3 only if the newly generated c is accepted in the MH algorithm.

4. Empirical Applications

4.1. Data description

The proposed SCD model is estimated with tick data of Toyota Motor Corporation (TMC) (Aichi, Japan), Asahi Group Holdings, Ltd. (AGH) (Tokyo, Japan) and Japan Airlines Co., Ltd. (JAL) (Tokyo, Japan) for 4 January 2016, 23 June 2016 and 27 July 2016 obtained from Nikkei NEEDS. The dates are chosen because 4 January was the first trading day of the year, 23 June was the day before the Brexit vote and 27 July was the day before the Bank of Japan’s policy meeting.

The daily trading session on the TSE, where these three companies are listed, is divided into two parts: the morning session is from 9:00 to 11:30, and the afternoon session is from 12:30 to 15:00. In this study, only the morning session is treated, and 9:00 is normalized to $t=0$ and 11:30 to $t=1$. As in the study by Bauwens and Veredas (2004), we exclude data that occurred at the same time. This is based on the assumption that the simultaneously recorded trades are only executed by small orders on the board due to one large order, i.e., they are caused by a single event. Since the TSE tick-data is recorded in microseconds, the percentage of null duration is very small, as shown in

Table 1. Basic information comparison for Toyota Motor Corporation (TMC) (Aichi, Japan), Asahi Group Holdings, Ltd. (AGH) (Tokyo, Japan) and Japan Airlines Co., Ltd. (JAL) (Tokyo, Japan).

	Date	No. Durations > 0	No. Durations = 0
TMC	4 January 2016	8086	579
	23 June 2016	8263	291
	27 July 2016	8359	457
AGH	4 January 2016	1385	50
	23 June 2016	1092	33
	27 July 2016	3581	151
JAL	4 January 2016	2755	107
	23 June 2016	1308	38
	27 July 2016	1814	21

. For example, Veredas et al. (2002) find that 26.5% of the data is null duration.

4.2. Empirical Results

The results in this section are produced by 20,000 iterations and 10,000 burn-ins. The setting of hyperparameters is below

$$(a_{\varphi },b_{\varphi })=(245,5),(a_{\tau },b_{\tau })=(5,0.04),({\gamma }_a,{\gamma }_b)=(1,1).$$

Since our proposed models include the parameters of B-Spline, the smoothing parameter $\lambda$ and the number of knots k, we need to choose a model from candidate parameters $\lambda =\{5,10,100\}$ and $k=\{150,30,15,10\}$. The number of k means that we set knots every 1, 5, 10 and 30 minutes, respectively. As shown in Table 1, the trading volume of TMC is several times larger than that of the two companies AGH and JAL; thus we set k = (150, 30, 15) for TMC and k = (30, 15, 10) for AGH and JAL as candidates. We use the Deviance information criterion (DIC) proposed by Spiegelhalter et al. (2002) as a criterion for model selection.

Table 2. DIC for Weibull distribution (4 January 2016, TMC, Aichi, Japan).

	Smoothing Parameter
Num of Knots	5	10	100
150	−146,550	−146,558	−146,578
30	−146,562	−146,571	−146,560
15	−146,562	−146,546	−146,535

The best of each distribution is underlined, and the best overall is in bold.

,

Table 3. DIC for Gamma distribution (4 January 2016, TMC, Aichi, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
150	−146,387	−146,455	−146,437
30	−146,417	−146,434	−146,423
15	−146,454	−146,422	−146,389

The best of each distribution is underlined, and the best overall is in bold.

,

Table 4. DIC for Weibull distribution (23 June 2016, TMC, Aichi, Japan).

	Smoothing Parameter
Num of Knots	5	10	100
150	−152,178	−152,175	−152,148
30	−152,114	−152,102	−152,103
15	−152,093	−152,099	−152,041

The best of each distribution is underlined, and the best overall is in bold.

,

Table 5. DIC for Gamma distribution (23 June 2016, TMC, Aichi, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
150	−151,919	−150,961	−151,948
30	−151,938	−151,950	−151,896
15	−151,953	−151,961	−151,898

The best of each distribution is underlined, and the best overall is in bold.

,

Table 6. DIC for Weibull distribution (27 July 2016, TMC, Aichi, Japan).

	Smoothing Parameter
Num of Knots	5	10	100
150	−153,018	−153,018	−153,021
30	−153,024	−153,028	−153,027
15	−153,001	−153,012	−153,001

The best of each distribution is underlined, and the best overall is in bold.

and

Table 7. DIC for Gamma distribution (27 July 2016, TMC, Aichi, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
150	−152,932	−152,927	−152,921
30	−152,924	−152,925	−152,950
15	−152,901	−152,946	−152,935

The best of each distribution is underlined, and the best overall is in bold.

show the DIC for each setting of TMC on 4 January 2016, 23 June 2016 and 27 July 2016, respectively. The best of each distribution is underlined, and the best overall is in bold. Using the variable selections in bold, the MCMC result of the limit order book information and plots of intraday seasonality are shown. The estimated parameters about the limit order book information are reported in

Table 8. Estimated parameters about limit order book (4 January 2016, TMC, Aichi, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
Spread	0.8862	0.0390	[0.8098 0.9623]
Volume	$-3.2627\times {10}^{-6}$	$2.4334\times {10}^{-6}$	[$-8.0522\times {10}^{-6}1.5150\times {10}^{-6}$]
Bid-ask ratio	−0.1706	0.1067	[−0.3799 0.0404]

,

Table 9. Estimated parameters about limit order book (23 June 2016, TMC, Aichi, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
Spread	1.4483	0.0560	[1.3389 1.5553]
Volume	$-5.8621\times {10}^{-6}$	$4.5002\times {10}^{-2}$	[$-1.4395\times {10}^{-5}2.6732\times {10}^{-6}$]
Bid-ask ratio	0.8105	0.1151	[0.5866 1.0511]

and

Table 10. Estimated parameters about limit order book (27 July 2016, TMC, Aichi, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
Spread	1.2863	0.0480	[1.1951 1.3798]
Volume	$-1.8708\times {10}^{-5}$	$3.6379\times {10}^{-6}$	[$-2.5667\times {10}^{-5}-1.1645\times {10}^{-5}$]
Bid-ask ratio	0.0093	0.1152	[−0.2199 0.2328]

. It is intuitive that the spread parameter is positive. This indicates that the larger the spread, the more open the trading interval. The effect of trading volume is negative only on 27 July. The effect of bid–ask ratio is positive only on 23 June. The 95% of CIs for the other items include 0, thus no impact can be claimed. Figure 1, Figure 2 and Figure 3 display the curve of intraday seasonality on 4 January, 23 June and 27 July 2016. Figure 1 takes the typical form of the inversed U-shape; however, Figure 3 shows that the duration becomes sparse toward noon. This is likely due to a change in trading behavior ahead of the Bank of Japan’s policy meeting.

Table 11. DIC for Weibull distribution (4 January 2016, AGH, Tokyo, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
30	−22,290	−22,200	−22,289
15	−22,112	−22,135	−22,273
10	−22,059	−22,120	−22,038

The best of each distribution is underlined, and the best overall is in bold.

,

Table 12. DIC for Gamma distribution (4 January 2016, AGH, Tokyo, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
30	−22,084	−22,146	−22,111
15	−22,129	−22,140	−22,059
10	−22,125	−22,111	−22,092

The best of each distribution is underlined, and the best overall is in bold.

,

Table 13. DIC for Weibull distribution (23 June 2016, AGH, Tokyo, Japan).

	Smoothing Parameter
Num of Knots	5	10	100
30	−16,633	−16,637	−16,632
15	−16,632	−16,636	−16,635
10	−16,638	−16,639	−16,628

The best of each distribution is underlined, and the best overall is in bold.

,

Table 14. DIC for Gamma distribution (23 June 2016, AGH, Tokyo, Japan).

	Smoothing Parameter
Num of Knots	5	10	100
30	−16,604	−16,606	−16,604
15	−16,607	−16,613	−16,608
10	−16,609	−16,606	−16,601

The best of each distribution is underlined, and the best overall is in bold.

,

Table 15. DIC for Weibull distribution (27 July 2016, AGH, Tokyo, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
30	−62,228	−62,231	−62,234
15	−62,244	−62,218	−62,235
10	−62,235	−62,238	−62,216

The best of each distribution is underlined, and the best overall is in bold.

and

Table 16. DIC for Gamma distribution (27 July 2016, AGH, Tokyo, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
30	−62,267	−62,274	−62,254
15	−62,265	−62,275	−62,260
10	−62,273	−62,252	−62,261

The best of each distribution is underlined, and the best overall is in bold.

show the DIC for each setting of AGH on 4 January 2016, 23 June 2016 and 27 July 2016, respectively. In

Table 17. Estimated parameters about limit order book (4 January 2016, AGH, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
Spread	0.9452	0.1127	[0.7121 1.1289]
Volume	$-2.6640\times {10}^{-4}$	$1.0761\times {10}^{-4}$	[$-4.6876\times {10}^{-4}-5.6517\times {10}^{-5}$]
Bid-ask ratio	−0.4060	0.4020	[−1.2842 0.3290]

,

Table 18. Estimated parameters about limit order book (23 June 2016, AGH, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
Spread	0.8802	0.1100	[0.6645 1.0960]
Volume	$-3.3787\times {10}^{-4}$	$1.4132\times {10}^{-4}$	[$-6.0885\times {10}^{-4}-6.1209\times {10}^{-5}$]
Bid-ask ratio	0.2303	0.3793	[−0.5134 0.9776]

and

Table 19. Estimated parameters about limit order book (27 July 2016, AGH, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
Spread	1.0511	0.0791	[0.9201 1.1741]
Volume	$-9.9515\times {10}^{-5}$	$3.4337\times {10}^{-5}$	[$-1.6464\times {10}^{-4}2.8444\times {10}^{-5}$]
Bid-ask ratio	−0.2928	0.1532	[−0.5582 0.0043]

, as in the case of TMC, the effect of the spread is positive on all dates. In other words, the larger the spread, the longer the interval between the next transaction. The trading volume parameters are negative on 4 January and 27 July, while the bid–ask ratio does not appear to affect the interval on all dates. Negative impact of trading volume indicates that traders may tend to falter or wait and see when large volume limit orders are placed. Figure 4, Figure 5 and Figure 6 display the curve of intraday seasonality on 4 January, 23 June and 27 July 2016. The credible interval on 4 January and 23 June is wider, but this is most likely due to the small sample size, as Table 1 shows.

Table 20. DIC for Weibull distribution (4 January 2016, JAL, Tokyo, Japan).

	Smoothing Parameter
Num of Knots	5	10	100
30	−47,826	−47,808	−47,803
15	−47,840	−47,823	−47,813
10	−47,817	−47,804	−47,765

The best of each distribution is underlined, and the best overall is in bold.

,

Table 21. DIC for Gamma distribution (4 January 2016, JAL, Tokyo, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
30	−47,845	−47,861	−47,843
15	−47,850	−47,806	−47,827
10	−47,831	−47,820	−47,791

The best of each distribution is underlined, and the best overall is in bold.

,

Table 22. DIC for Weibull distribution (23 June 2016, JAL, Tokyo, Japan).

	Smoothing Parameter
Num of Knots	5	10	100
30	−16,633	−16,637	−16,632
15	−16,632	−16,636	−16,635
10	−16,638	−16,639	−16,628

The best of each distribution is underlined, and the best overall is in bold.

,

Table 23. DIC for Gamma distribution (23 June 2016, JAL, Tokyo, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
30	−16,604	−16,606	−16,604
15	−16,607	−16,613	−16,608
10	−16,609	−16,606	−16,601

The best of each distribution is underlined, and the best overall is in bold.

,

Table 24. DIC for Weibull distribution (27 July 2016, JAL, Tokyo, Japan).

	Smoothing Parameter
Num of Knots	5	10	100
30	−29,421	−29,422	−29,422
15	−29,418	−29,420	−29,416
10	−29,419	−29,418	−29,406

The best of each distribution is underlined, and the best overall is in bold.

and

Table 25. DIC for Gamma distribution (27 July 2016, JAL, Tokyo, Japan).

	Smoothing Parameters
Num of Knots	5	10	100
30	−29,481	−29,499	−29,482
15	−29,524	−29,483	−29,461
10	−29,491	−29,477	−29,468

The best of each distribution is underlined, and the best overall is in bold.

show the DIC for each setting of JAL on 4 January 2016, 23 June 2016 and 27 July 2016, respectively. The estimation results of the limit order book information parameters in

Table 26. Estimated parameters about limit order book (4 January 2016, JAL, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
Spread	1.1827	0.0529	[1.0809 1.2932]
Volume	$-6.5667\times {10}^{-5}$	$3.4338\times {10}^{-5}$	[$-1.3631\times {10}^{-4}1.9506\times {10}^{-7}$]
Bid-ask ratio	0.0616	0.2112	[−0.3058 0.5194]

,

Table 27. Estimated parameters about limit order book (23 June 2016, JAL, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
Spread	1.4040	0.1372	[1.1384 1.6733]
Volume	$-3.0252\times {10}^{-5}$	$1.2035\times {10}^{-5}$	[$-5.3968\times {10}^{-5}-7.0058\times {10}^{-6}$]
Bid-ask ratio	0.0729	0.3496	[−0.6238 0.7474]

and

Table 28. Estimated parameters about limit order book (27 July 2016, JAL, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
Spread	1.0004	0.0469	[0.9136 1.0864]
Volume	$-5.9807\times {10}^{-5}$	$2.8536\times {10}^{-5}$	[$-1.1086\times {10}^{-4}-3.9883\times {10}^{-6}$]
Bid-ask ratio	−0.4280	0.1643	[−0.7945 −0.1419]

show almost the same results as those for TMC and AGH. Specifically, the spread parameters are positive on all dates, the volume parameters are negative on 23 June and 27 July, and the bid-ask ratio parameter is positive only on 27 July. Figure 7, Figure 8 and Figure 9 display the curve of intraday seasonality on 4 January, 23 June and 27 July 2016. The credible interval on 23 June is wider, but this is most likely due to the small sample size, as Table 1 shows.

To summarize the results of the estimation of the parameters of the limit order book information for the three stocks, the spread has a positive impact in all nine patterns of the three dates for the three stocks. We argue that the fact that trades are more likely to occur when price differentials are narrow is a near-obvious result from the perspective of understanding the order generation mechanism, but it is a useful result from the perspective of analyzing and modeling the time interval of trades based on the limit order book information. As for volume, five of the nine patterns show a negative impact, indicating that large volume at the best quote may not lead to active trading. Unlike the spread, there are some patterns in which the effect of volume is not pronounced. Therefore, identifying under what conditions the effect of volume changes is a future issue. The bid–ask ratio is positive in one of the nine patterns and negative in one pattern. The ratio of buy to sell order volume does not seem to affect the time interval.

4.3. Model Comparison

For the models selected in Section 4.2, we compare their performance with existing methods. The comparison is made with the SCD model based on Bauwens and Veredas (2004) (the original SCD). Following this paper, intraday seasonality is removed in advance using the Nadaraya–Watson estimator. The kernel is the quartic, and the bandwidth h is set to $2.78s{N}^{-1/5}$ where s is the standard deviation of the data.

The DIC and the widely applicable information criterion (WAIC) proposed by Watanabe (2013) are used as the model selection criteria.

Out of a total of eight patterns shown in

Table 29. Comparison of information criteria (TMC, Aichi, Japan).

	DIC	WAIC
proposed-SCD (4 January 2016)	−146,578	−143,175
original-SCD (4 January 2016)	−137,586	−137,275
proposed-SCD (23 June 2016)	−152,178	−156,533
original-SCD (23 June 2016)	−151,017	−152,319
proposed-SCD (27 July 2016)	−153,028	−144,322
original-SCD (27 July 2016)	−143,312	−143,084

,

Table 30. Comparison of information criteria (AGH, Tokyo, Japan).

	DIC	WAIC
proposed-SCD (4 January 2016)	−22,290	−22,132
original-SCD (4 January 2016)	−25,160	−25,239
proposed-SCD (23 June 2016)	−16,639	−16,647
original-SCD (23 June 2016)	−20,711	−20,675
proposed-SCD (27 July 2016)	−62,275	−62,025
original-SCD (27 July 2016)	−62,786	−62,015

and

Table 31. Comparison of information criteria (JAL, Tokyo, Japan).

	DIC	WAIC
proposed-SCD (4 January 2016)	−47,861	−49,786
original-SCD (4 January 2016)	−49,472	−49,284
proposed-SCD (23 June 2016)	−16,639	−21,825
original-SCD (23 June 2016)	-	-
proposed-SCD (27 July 2016)	−29,524	−32,512
original-SCD (27 July 2016)	−32613	−32,402

, the proposed model outperformed in the WAIC in six patterns and the model with a priori removal of intraday seasonality in the DIC in five patterns. This may be due to the better generalization performance of the proposed model with simultaneous estimation, since the DIC tends to select overfitted models, as indicated by Chan and Grant (2016). The original model for JAL on 23 June could not be estimated because the matrix was not positive definite.

Table 32. Comparison of estimated parameters (4 January 2016, TMC, Aichi, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
${\varphi }_{proposed}$	0.8215	0.0211	[0.7772    0.8600]
${\gamma }_{proposed}$	0.5391	0.0105	[0.5202    0.5614]
${\sigma }_{proposed}$	0.8994	0.0706	[0.7708    1.0492]
${\varphi }_{original}$	0.9995	0.0001	[0.9993    0.9998]
${\gamma }_{original}$	0.5216	0.0048	[0.5122    0.5309]
${\sigma }_{original}$	0.1883	0.0088	[0.1717    0.2055]

,

Table 33. Comparison of estimated parameters (23 June 2016, TMC, Aichi, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
${\varphi }_{proposed}$	0.5812	0.0188	[0.5475    0.6195]
${\gamma }_{proposed}$	0.7672	0.0231	[0.7238    0.8093]
${\sigma }_{proposed}$	1.8313	0.0516	[1.7316    1.9234]
${\varphi }_{original}$	0.9970	0.0007	[0.9955    0.9982]
${\gamma }_{original}$	0.5088	0.0065	[0.4963    0.5219]
${\sigma }_{original}$	0.8340	0.0318	[0.7706    0.8988]

,

Table 34. Comparison of estimated parameters (27 July 2016, TMC, Aichi, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
${\varphi }_{proposed}$	0.7213	0.0233	[0.6754    0.7620]
${\gamma }_{proposed}$	0.6168	0.0172	[0.5863    0.6506]
${\sigma }_{proposed}$	1.3649	0.0755	[1.2317    1.5081]
${\varphi }_{original}$	0.9993	0.0002	[0.9988    0.9996]
${\gamma }_{original}$	0.5083	0.0047	[0.4991    0.5176]
${\sigma }_{original}$	0.2815	0.0125	[0.2587    0.3080]

,

Table 35. Comparison of estimated parameters (4 January 2016, AGH, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
${\varphi }_{proposed}$	0.9116	0.0264	[0.8566    0.9579]
${\gamma }_{proposed}$	0.3706	0.0110	[0.3497    0.3930]
${\sigma }_{proposed}$	0.5502	0.1285	[0.3093    0.8022]
${\varphi }_{original}$	0.9875	0.0033	[0.9877    0.9934]
${\gamma }_{original}$	0.3776	0.0115	[0.3557    0.4006]
${\sigma }_{original}$	1.3037	0.1065	[1.1071    1.5176]

,

Table 36. Comparison of estimated parameters (23 June 2016, AGH, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
${\varphi }_{proposed}$	0.9650	0.0150	[0.9291    0.9874]
${\gamma }_{proposed}$	0.3513	0.0087	[0.3345    0.3685]
${\sigma }_{proposed}$	0.1179	0.0366	[0.0647    0.2062]
${\varphi }_{original}$	0.9916	0.0025	[0.9862    0.9960]
${\gamma }_{original}$	0.3362	0.0099	[0.3173    0.3559]
${\sigma }_{original}$	1.0119	0.0888	[0.8462    1.1915]

,

Table 37. Comparison of estimated parameters (27 July 2016, AGH, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
${\varphi }_{proposed}$	0.2824	0.0153	[0.2494    0.3066]
${\gamma }_{proposed}$	1.3463	0.0952	[1.1919    1.5324]
${\sigma }_{proposed}$	2.6611	0.0395	[2.5852    2.7386]
${\varphi }_{original}$	0.9958	0.0021	[0.9904    0.9983]
${\gamma }_{original}$	0.3944	0.0115	[0.3692    0.4114]
${\sigma }_{original}$	0.6704	0.1291	[0.5394    0.9674]

,

Table 38. Comparison of estimated parameters (4 January 2016, JAL, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
${\varphi }_{proposed}$	0.3760	0.0180	[0.3402    0.4128]
${\gamma }_{proposed}$	1.9805	0.1875	[1.6585    2.3511]
${\sigma }_{proposed}$	2.9829	0.0447	[2.8972    3.0719]
${\varphi }_{original}$	0.9900	0.0024	[0.9850    0.9944]
${\gamma }_{original}$	0.3727	0.0092	[0.3553    0.3909]
${\sigma }_{original}$	1.1598	0.0732	[1.0233    1.3094]

,

Table 39. Comparison of estimated parameters (23 June 2016, JAL, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
${\varphi }_{proposed}$	0.9140	0.0311	[0.8452    0.9664]
${\gamma }_{proposed}$	0.3455	0.0121	[0.3239    0.3714]
${\sigma }_{proposed}$	0.6072	0.1817	[0.2827    0.9946]
${\varphi }_{original}$	-	-	[-    -]
${\gamma }_{original}$	-	-	[-    -]
${\sigma }_{original}$	-	-	[-    -]

and

Table 40. Comparison of estimated parameters (27 July 2016, JAL, Tokyo, Japan).

Parameters	Posterior Mean	Posterior SD	95%CI
${\varphi }_{proposed}$	0.3277	0.0155	[0.2911    0.3577]
${\gamma }_{proposed}$	1.1510	0.1310	[0.8505    1.3794]
${\sigma }_{proposed}$	3.1048	0.0667	[2.9681    3.2322]
${\varphi }_{original}$	0.9977	0.0007	[0.9963    0.9990]
${\gamma }_{original}$	0.2459	0.0066	[0.2332    0.2590]
${\sigma }_{original}$	0.3712	0.0248	[0.3228    0.4210]

compare the results of the parameter estimation for each of the three companies in the proposed model and the original model. In all eight patterns, AR(1) Coefficient $\varphi$ of the original model exceeds that of the proposed model, and it can be observed that the prior intraday seasonality removal does not successfully remove the persistence of the state variable. The shape parameter $\gamma$ changes significantly in Table 37, Table 38 and Table 40 where the Gamma distribution is selected, while it changes little in the remaining five patterns where the Weibull distribution is selected. The variance $\sigma$ of the proposed model increased in six out of eight patterns compared to that of the original model. This may be due to the fact that the extreme rigidity of the state variable $\alpha$ in the original model has been eliminated because $\varphi$ no longer takes values close to one. The two patterns in Table 35 and Table 36 with decreased $\sigma$ in the proposed model are consistent with the two highest when $\varphi$ is compared among the proposed models (0.9650, 0.9116), supporting the above interpretation.

5. Conclusions

This paper proposed an extension of the SCD model for estimating the intraday seasonality and duration clustering simultaneously, while previous studies remove this intraday seasonality externally. The intraday seasonality is represented by a B-Spline, and its coefficients are simulated along with the model parameters in the MCMC loop. This has the advantage not only of allowing simultaneous estimation of the intraday seasonality, but also of facilitating the extension to the estimation of the effect of external information such as the limit order information. One limitation of the method in this paper is that the number of knots and smoothing parameters must be selected based on an information criterion, which is computationally expensive. Future research should address these issues and reduce the computational cost. We also use the block-sampling for state variable and change a setting for the parameter in the Gamma distribution for a stable sampling.

In the empirical analysis, the proposed SCD model is applied to three days of microsecond tick-data for three stocks listed on the TSE stock market. The proposed method outperforms the existing method in the WAIC sense in many of the estimation results for the three stocks with different trading frequencies, consistent with previous studies. Since very few studies have incorporated the limit order book information into the SCD models, the use of the limit order book information and discussion of the corresponding parameter estimation results are also contributions of this paper. We find that the spread has a positive impact on the trading time interval for all patterns, while order size has a negative impact for about half of the patterns. Future studies will confirm the superiority of simultaneous estimation in forecasting with a wide range of data.