Bayesian modelling of integer-valued transfer function models

Abstract

External events are commonly known as interventions that often affect times series of counts. This research introduces a class of transfer function models that include four different types of interventions on integer-valued time series: abrupt start and abrupt decay (additive outlier), abrupt start and gradual decay (transient shift), abrupt start and permanent effect (level shift) and gradual start and permanent effect. We propose integer-valued transfer function models incorporating a generalized Poisson, log-linear generalized Poisson or negative binomial to estimate and detect these four types of interventions in a time series of counts. Utilizing Bayesian methods, which are adaptive Markov chain Monte Carlo (MCMC) algorithms to obtain the estimation, we further employ deviance information criterion (DIC), posterior odd ratios and mean squared standardized residual for model comparisons. As an illustration, this study evaluates the effectiveness of our methods through a simulation study and application to crime data in Albury City, New South Wales (NSW) Australia. Simulation results show that the MCMC procedure is reasonably effective. The empirical outcome also reveals that the proposed models are able to successfully detect the locations and type of interventions.

1 Introduction

One of the major challenges in time series modelling is the presence of outliers as the series is often affected by external events commonly known as interventions. These various external events can be major corporate, political or economic policy initiatives or changes, technological advancements, work stoppages, sales promotions, advertising and so forth. These interventions influence the time series data either on a single or a few datapoints, or they impact the whole process from some specific time $T$ . (Box and Tiao (1975)) pioneer intervention detection and estimation analysis to solve the problem of Los Angeles pollution. Notable contributions and extensions concerning outliers in time series analysis include those by (Chang et al. (1988)), (Chen and Liu (1993)), (Chareka et al. (2006)) and others. Both known and unknown effects have practical aspects. When the change point is known, it is an intervention effect; otherwise, it is a transfer function for the unknown locations.

Modelling time series of counts is an important task in research today and has garnered significant attention for decades with applications in various scientific areas, such as health science, economics, finance, environment, criminology, epidemiology, etc. (McKenzie (1985)) and (Al-Osh and Alzaid (1987)) introduce the first-order non-negative integer-valued autoregressive (INAR) model that is based on a binomial thinning operator. Most count data models make use of the Poisson distribution, but it is unable to describe over-dispersion. (Ferland et al. (2006)) propose the integer-valued generalized autoregressive conditional heteroscedastic (GARCH) model to handle the over-dispersion problem, and subsequently (Fokianos and Tjøstheim (2011)) consider a log-linear Poisson model for time series of counts. Our study therefore introduces a class of transfer function models for time series of counts that covers interventions in integer-valued GARCH processes as special cases.

(Fokianos and Fried (2010)) and (Fokianos and Fried (2012)) present three types of intervention in their proposed models, which are additive outliers, transient shift and level shift for Poisson and log linear Poisson integer-valued GARCH processes. (Liboschik et al. (2014)) also look into these three types of intervention for the Poisson integer-valued GARCH process both for a known time and an unknown time of intervention. (Chen and Lee (2016)) target a class of zero-inflated generalized Poisson (GP) integer-valued GARCH models with a structural break. We introduce herein a class of transfer function models that incorporate four different types of intervention on integer-valued time series processes: abrupt start and abrupt decay (additive outlier), abrupt start and gradual decay (transient shift), abrupt start and permanent effect (level shift) and gradual start and permanent effect (see Chapter 14 in Wei, [2006], for the ARIMA process). Our study includes intervention effects from the integer-valued GARCH models of (Fokianos and Fried (2010)) and (Liboschik et al. (2014)). The novelty of this study's contribution is that we detect the four types of intervention effects by incorporating a ratio of two finite (rational) polynomials in the proposed models, which treat other models as special cases in our models.

The traditional GARCH(1,1) model does not show good performance for frequentist estimations when run over a sample size of 500 or 1,000 (see Karmakar and Roy (2021)). However, Bayesian formulation alleviates this problem as it does not require a huge sample size. A similar context appears for an integer-valued GARCH model. We offer Bayesian methods to detect and estimate those four types of intervention effects for time series of counts with GP or negative binomial (NB) distributions. GP is a mixture of Poisson distribution (Joe and Zhu (2005)) and a very versatile discrete distribution with applications in various areas of study. The advantage of the GP integer-valued GARCH model is that it accommodates the over-dispersion problem (Chen and Lee (2016)). The NB distribution frequently appears in several research studies to describe over-dispersed count data, such as (Zhu (2011)), (Chen and Lee (2017)), (Chen and Khamthong (2020)) and (Chen et al. (2021)), among others. We target to estimate and detect the locations of interventions following GP, log-linear GP and NB integer-valued transfer function models for time series count data in a Bayesian framework.

Many papers in the literature have successfully employed Bayesian techniques to make inferences or to forecast the class of an integer-valued GARCH family; see (Fried et al. (2015)), (Chen and Lee (2016)), (Chen and Lee (2017)), (Chen et al. (2019)), (Chen et al. (2021)), etc. We also employ Bayesian techniques for detection and estimation based on adaptive Markov chain Monte Carlo (MCMC) methods, which have many advantages as follows: (a) Bayesian methods do not rely on the asymptotic property that can be an obstacle when employing frequentist methods in small-sample situations; (b) these methods allow for simultaneous estimation of all unknown parameters and locations of interventions; (c) they enable efficient and flexible handling of complex models; and (d) they properly impose parameter constraints as part of the prior distribution.

We employ deviance information criterion (DIC) by (Spiegelhalter et al. (2002)) and mean squared standardized residual (MSSR) for model comparisons. In the context of Bayesian inference, one may frame hypothesis testing as a special case of model comparison. This study further presents the use of the posterior odds ratio (POR) to test one or multiple interventions at unknown locations, which shows the posterior relevance of the two hypotheses as a single number.

The rest of the article runs as follows. Section 2 introduces the GP, log-linear GP and NB transfer function models. Section 3 presents Bayesian MCMC methods to estimate the parameters and addresses model selection for the specified models. Section 4 performs a simulation study for illustration. In particular, we study misspecification of the distribution and the autoregressive order. Section 5 implements the proposed transfer function models to three categories of crime data in Albury City, New South Wales Australia (NSW). Section 6 provides concluding remarks.

2 Transfer function models

This section presents the development of the proposed transfer function models. Here, we incorporate the four types of intervention effects to GP, log-linear GP and NB integer-valued transfer function models. A random variable $X$ is said to have a GP distribution denoted by $P_x(\lambda ,\psi )$ with parameters $\lambda$ and $\psi$ if:

$$\begin{array}{ccc}\hfill P_x(\lambda ,\psi )=\left\{\begin{array}{cc}\lambda (\lambda +\psi x{)}^{x-1}{e}^{-\lambda -\psi x}/x!,\hfill & \mathit{for}x=0,1,2,\dots \hfill \\ 0,\hfill & \mathit{for}x>m\mathit{if}\psi <0,\hfill \end{array}\right.& \hfill \hfill & \hfill \end{array}$$

where $\text{max}(-1,-\lambda /m)<\psi <1,$ and $m(\ge 4)$ is the largest positive integer for which $\lambda +\psi m>0.$ Let $\left\{X_t\right\}$ be an intervention-free integer-valued GARCH(1,1) process such that:

$$X_t|{\mathit{X}}_{t-1}\sim \mathit{GP}({\lambda }_t,\psi ),:{\lambda }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\lambda }_{t-1},$$

(2.1)

where ${\mathit{X}}_{t-1}$ denotes the past information available up to time $t-1$ , and the parameters $\lambda >0$ , ${\alpha }_0>0,:{\alpha }_1,{\beta }_1\ge 0$ and ${\alpha }_1+{\beta }_1<1$ . We assume $0\le \psi <1$ to emphasize the over-dispersion situation. According to (Lee et al. (2015)), $\left\{X_t\right\}$ is strictly stationary and ergodic with all finite moments.

2.1 GP transfer function model

The first two intervention effects on conditional intensity are ‘abrupt start abrupt decay’ and ‘abrupt start gradual decay’, which are respectively additive outlier and transient shifts.

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{GP}({\kappa }_t,\psi ),\hfill \\ \hfill & \hfill \hfill & {\kappa }_t={\lambda }_t+\frac{{\omega }_s(B)}{{\delta }_q(B)}{P}_t,\hfill \\ \hfill & \hfill \hfill & {\lambda }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\lambda }_{t-1},\hfill \\ \hfill & \hfill \hfill & \hfill \end{array}$$

(2.2)

$$\begin{array}{ccc}\hfill & \hfill \hfill & P_t=\left\{\begin{array}{cc}1,\hfill & \mathit{if}t=T\hfill \\ 0,\hfill & \mathit{if}t\ne T,\hfill \end{array}\right.\hfill \\ \hfill & \hfill \hfill & \hfill \end{array}$$

(2.3)

$$\begin{array}{ccc}\hfill & \hfill \hfill & {\omega }_s(B)={\omega }_0-{\omega }_{1}B-\dots -{\omega }_s{B}^s\mathit{and}{\delta }_q(B)=1-{\delta }_{1}B-\dots -{\delta }_q{B}^q.\hfill \end{array}$$

(2.4)

We call ${\omega }_s(B)/{\delta }_q(B)$ a transfer function.

The set-ups for Inv-1 and Inv-2 are the same as (Liboschik et al. (2014)) where the intervention effects are not propagated via the feedback mechanism of the conditional intensity, but only via the contaminated observations. Let $f_t(.)={\omega }_s(B)/{\delta }_q(B)$ . The impact from Inv-2, $f_t(.)$ , can be stated as follows:

$$\begin{array}{ccc}\hfill f_t(.)=\left\{\begin{array}{cc}0,\hfill & t<T\hfill \\ {\omega }_0{\delta }^{t-T},\hfill & t\ge T.\hfill \end{array}\right.& \hfill \hfill & \hfill \end{array}$$

Inv-1 corresponds to the additive outlier in conditional intensity at time $T$ . It only affects a single observation. Inv-2 corresponds to transient shifts where $\delta \in (0,1)$ . It results in effects that decay with rate $\delta$ . On the other hand, following the set-up of (Fokianos and Fried (2010)) we consider intervention effects added to the underlying conditional intensity process ${\kappa }_t$ .

Let $X_t=Y_t+X_t^{*}$ be the contaminated process, $Y_t$ is the intervention-free process with conditional intensity ${\lambda }_t$ and $X_t^{*}$ is a sequence of Poisson random variables with conditional intensity depending on the transfer function. We replace Eq. (2.2) by the following:

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{GP}({\kappa }_t,\psi ),\hfill \\ \hfill \mathit{Inv}-2\mathit{A}& \hfill \hfill & {\kappa }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\kappa }_{t-1}+{\omega }_0{P}_t,\hfill \\ \hfill & \hfill \hfill & \hfill \end{array}$$

(2.5)

$$\begin{array}{ccc}\hfill \mathit{Inv}-2\mathit{B}& \hfill \hfill & {\kappa }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\kappa }_{t-1}+\frac{{\omega }_0}{(1-\delta B)}{P}_t.\hfill \end{array}$$

(2.6)

(Inv-2A): For Inv-2A, κ _t can be expressed as follows:

$${\kappa }_t=\{\begin{array}{c}{\lambda }_t,t<T\\ {\lambda }_t+{\displaystyle {\sum }_{i=0}^{t-1}{\beta }_1^i{X}_{t-(i+1)}^{\ast }+{\omega }_0{\beta }_1^{t-T},}t\ge T\end{array}$$

(Inv-2B): For Inv-2B, κ _t can be derived as follows:

$${\kappa }_t=\{\begin{array}{c}{\lambda }_t,t<T\\ {\lambda }_t+{\displaystyle {\sum }_{i=0}^{t-1}{\beta }_1^i{X}_{t-(i+1)}^{\ast }+{\omega }_0{\displaystyle {\sum }_{i=0}^t{\beta }_1^i{\delta }^{t-i-T},}}t\ge T.\end{array}$$

Note that $X_t^{*}=0$ when $t<T$ . The impact from Inv-2A depends on the value of ${\beta }_1$ , in which it decreases rapidly for a smaller value of ${\beta }_1$ . Moreover, the impact from Inv-2B becomes gradually smaller as time grows, and the size of impact decreases according to $\delta$ and ${\beta }_1$ . Each subplot in Figure 1 displays the time plot of the intervention-free integer-valued GARCH(1,1) process and the contaminated process after $T=50$ under the GP integer-valued transfer function model. The impact from Inv-2A affects 7 time points after the time of intervention at $T=50$ , while Inv-2B affects 12 time points after the intervention in which its effect is an abrupt increase and subsequently a gradual decrease of the time series starting at $T=50$ .

We can add another type of intervention at a conditional intensity, known as level shift and gradual start with permanent effect based on a step function, $S_t$ .

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{GP}({\kappa }_t,\psi ),\hfill \\ \hfill & \hfill \hfill & {\kappa }_t={\lambda }_t+\frac{{\omega }_s(B)}{{\delta }_q(B)}{S}_t,\hfill \end{array}$$

(2.7)

$$\begin{array}{ccc}\hfill & \hfill \hfill & {\lambda }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\lambda }_{t-1},\hfill \\ \hfill & \hfill \hfill & \hfill \end{array}$$

(2.8)

$$\begin{array}{ccc}\hfill & \hfill \hfill & S_t=\left\{\begin{array}{cc}0,\hfill & \mathit{if}t<T\hfill \\ 1,\hfill & \mathit{if}t\ge T,\hfill \end{array}\right.\hfill \end{array}$$

(2.9)

where the definitions of ${\omega }_s(B)$ and ${\delta }_q(B)$ are the same as in Eq. (2.4).

(Inv-3): Abrupt start permanent effect: (s, q) = (0, 0).

(Inv-4): Gradual start permanent effect: (s, q) = (0, 1).

Let $g_t(.)={\omega }_s(B)/{\delta }_q(B)$ for Inv-4. The impact from Inv-4, $g_t(.)$ , is as follows:

$$\begin{array}{ccc}\hfill g_t(.)=\left\{\begin{array}{cc}0,\hfill & t<T\hfill \\ {\omega }_0\left(\frac{1-{\delta }^{t+1-T}}{1-\delta }\right),\hfill & t\ge T,\hfill \end{array}\right.& \hfill \hfill & \hfill \end{array}$$

where the impact displays a gradual start and a permanent effect. When this intervention occurs in the conditional intensity, we have:

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{GP}({\kappa }_t,\psi ),\hfill \\ \hfill \mathit{Inv}-4\mathit{A}& \hfill \hfill & {\kappa }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\kappa }_{t-1}+{\omega }_0{S}_t,\hfill \\ \hfill & \hfill \hfill & \hfill \end{array}$$

(2.10)

$$\begin{array}{ccc}\hfill \mathit{Inv}-4\mathit{B}& \hfill \hfill & {\kappa }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\kappa }_{t-1}+\frac{{\omega }_0}{(1-\delta B)}{S}_t.\hfill \end{array}$$

(2.11)

(Inv-4A): For Inv-4A, κ _t can be obtained as follows:

$${\kappa }_t=\{\begin{array}{c}{\lambda }_t,t<T\\ {\lambda }_t+{\displaystyle {\sum }_{i=0}^{t-1}{\beta }_1^i{X}_{t-(i+1)}^{\ast }+{\omega }_0\left(\frac{1-{\beta }_1^{t+1}}{1-{\beta }_1}\right),}t\ge T.\end{array}$$

(Inv-4B): For Inv-4B, κ _t can be derived as follows:

$${\kappa }_t=\{\begin{array}{c}{\lambda }_t,t<T\\ {\lambda }_t+{\displaystyle {\sum }_{i=0}^{t-1}{\beta }_1^i{X}_{t-(i+1)}^{\ast }+{\omega }_0{\displaystyle {\sum }_{i=0}^t\left(\frac{{\beta }_1^{r+1-(i+T)}-{\delta }^{t+1-(i+T)}}{{\beta }_1-\delta }\right)},}t\ge T.\end{array}$$

The impact from both Inv-4A and Inv-4B displays a gradual start and a permanent effect, but at different speeds. The rate of level shift from Inv-4A is faster than that from Inv-4B. Figure 1 illustrates both intervention effects. One can clearly see the impact from Inv-4A in Figure 1 with an abrupt increase and permanent effect starting at $T=50$ , while there is a gradual increase and permanent effect after $T=50$ in Inv-4B.

Figure 1

Simulated data under the GP integer-valued transfer function with $({\alpha }_0,{\alpha }_1,{\beta }_1,\psi ,T)=$ $(1.0,0.3,0.3,0.18,50)$ and ${\omega }_0=10$ for Inv-2A; $({\omega }_0,\delta )=(10,0.4)$ for Inv-2B; ${\omega }_0=2$ for Inv-4A; and $({\omega }_0,\delta )=(2,0.4)$ for Inv-4B

2.2 Log-linear GP transfer function model

We now turn to the log-linear GP integer-valued transfer function models with a general form of intervention.

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{GP}({\kappa }_t,\psi ),:{\nu }_t^{*}\equiv \text{log}({\kappa }_t)\hfill \\ \hfill & \hfill \hfill & {\nu }_t^{*}={\nu }_t+\frac{{\omega }_s(B)}{{\delta }_q(B)}{I}_t,\hfill \end{array}$$

(2.12)

$$\begin{array}{ccc}\hfill & \hfill \hfill & {\nu }_t={\alpha }_0+{\alpha }_1\text{log}(X_{t-1}+1)+{\beta }_1{\nu }_{t-1}.\hfill \end{array}$$

(2.13)

We alternatively have the following:

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{GP}({\kappa }_t,\psi ),:{\nu }_t\equiv \text{log}({\kappa }_t)\hfill \\ \hfill & \hfill \hfill & {\nu }_t={\alpha }_0+{\alpha }_1\text{log}(X_{t-1}+1)+{\beta }_1{\nu }_{t-1}+\frac{{\omega }_s(B)}{{\delta }_q(B)}{I}_t,\hfill \end{array}$$

(2.14)

where $I_t$ is an indicator function, and the definitions of ${\omega }_s(B)$ and ${\delta }_q(B)$ are the same as in Eq. (2.4). The model in Eq. (2.12) includes the following interventions.

(Inv-2A):Abrupt start abrupt decay (additive outlier) at ν _T : (s, q) = (0, 0) and I _t = P _t .

(Inv-2B): Abrupt start gradual decay (transient shift): (s, q) = (0, 1) and I _t = P _t .

(Inv-4A): Gradual start permanent effect: (s, q) = (0, 0) and I _t = S _t .

(Inv-4B): Gradual start permanent effect: (s, q) = (0, 1) and I _t = S _t .

To ensure stability conditions in the original process, we restrict the following conditions from (Fokianos and Tjøstheim (2011)).

$$\begin{array}{ccc}\hfill & \hfill \hfill & |{\beta }_1|<1, {\alpha }_1>0, |{\alpha }_1+{\beta }_1|<1,\hfill \\ \hfill & \hfill \hfill & \hfill \end{array}$$

(2.15)

$$\begin{array}{ccc}\hfill & \hfill \mathit{or}\hfill & |{\beta }_1|<1, {\alpha }_1<0, |{\beta }_1||{\alpha }_1+{\beta }_1|<1.\hfill \end{array}$$

(2.16)

(Chen et al. (2019)) employ the condition in (2.13), because there is stronger dependence under it as discussed by (Fokianos and Tjøstheim (2011)).

2.3 NB transfer function model

The NB distribution has become popular in research studies due to its flexibility. We also adopt a NB integer-valued transfer function model herein. We assume that $X_t|{\mathit{X}}_{t-1}$ follows a NB integer-valued transfer function model (see Zhu (2011)).

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{NB}(r,p_t),:\frac{1-p_t}{p_t}:={\kappa }_t\hfill \\ \hfill & \hfill \hfill & {\kappa }_t={\lambda }_t+\frac{{\omega }_s(B)}{{\delta }_q(B)}{I}_t,\hfill \end{array}$$

(2.17)

$$\begin{array}{ccc}\hfill & \hfill \hfill & {\lambda }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\lambda }_{t-1}.\hfill \end{array}$$

(2.18)

We alternatively have the following:

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{NB}\left(r,\frac{1}{1+{\kappa }_t}\right)\hfill \\ \hfill & \hfill \hfill & {\kappa }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\kappa }_{t-1}+\frac{{\omega }_s(B)}{{\delta }_q(B)}{I}_t,\hfill \end{array}$$

(2.19)

where $I_t$ is an indicator function, and the definitions of ${\omega }_s(B)$ and ${\delta }_q(B)$ are the same as in Eq. (2.4). Again, the proposed model in Eq. (2.16) can also cover the following interventions.

(Inv-2A):Abrupt start abrupt decay (additive outlier) at κ _T : (s, q) = (0, 0) and I _t = P _t .

(Inv-2B): Abrupt start gradual decay (transient shift): (s, q) = (0, 1) and I _t = P _t .

(Inv-4A): Gradual start permanent effect: (s, q) = (0, 0) and I _t = S _t .

(Inv-4B): Gradual start permanent effect: (s, q) = (0, 1) and I _t = S _t .

(Zhu (2011)) states the stationary condition for a NB integer-valued GARCH model as:

$$r({\alpha }_1{)}^2+(r{\alpha }_1+{\beta }_1{)}^2<1.$$

(2.20)

We incorporate this condition into our prior.

3 Bayesian inference

Bayesian methods generally require the specification of a prior distribution $p(\theta )$ . After observing data $\text{X}=(X_1\dots ,X_n{)}^{\prime }$ , we update our beliefs and calculate the posterior distribution $p(\theta |\text{X})$ . There are two scenarios in the study: when $T$ is known as an intervention effect, and when $T$ is unknown for a transfer function. We focus on a much complicated situation: $T$ and $({\omega }_0,\delta )$ are unknown parameters in the transfer function model.

Let ${\mathbf{\theta }}_1$ denote all unknown parameters in the integer-valued transfer function model with GP distribution, where ${\theta }_1=({\beta }^{\prime },{\omega }_0,\delta ,\psi ,T{)}^{\prime }$ with $\beta =({\alpha }_0,{\alpha }_1,{\beta }_1{)}^{\prime }$ . The conditional log-likelihood function for the proposed model is as follows:

$$\text{ln}\mathcal{L}(\mathit{X}|{\mathbf{\theta }}_1)=\sum _{t=2}^n(\text{ln}{\kappa }_t+(X_t-1)\text{ln}({\kappa }_t+\psi X_t)-({\kappa }_t+\psi X_t)-\text{ln}{X}_t!),$$

(3.1)

where ${\kappa }_t$ includes a specified intervention effect. For a NB distribution, let ${\mathbf{\theta }}_2=({\beta }^{\prime },{\omega }_0,\delta ,r,T{)}^{\prime }$ be all parameters in (2.16), in which the conditional log-likelihood function is as follows:

$$\begin{array}{ccc}\hfill \text{ln}\mathcal{L}(\mathit{X}|{\mathbf{\theta }}_2)=\sum _{t=2}^n\text{ln}\left\{\frac{\Gamma (X_t+r)}{\Gamma (X_t+1)\Gamma (r)}\right\}+r\sum _{t=2}^n\text{ln}\left(\frac{1}{1+{\kappa }_t}\right)+\sum _{t=2}^n{X}_t\text{ln}\left(\frac{{\kappa }_t}{1+{\kappa }_t}\right).& \hfill \hfill & \hfill \\ \hfill & \hfill \hfill & \hfill \end{array}$$

(3.2)

We make use of the following parameter groups for both GP and log-linear GP integer-valued transfer function models: (a) $\mathbf{\beta }$ , (b) $({\omega }_0,\delta )$ , (c) $\psi$ and (d) $T$ . We replace (c) by $r$ when we consider a NB integer-valued transfer function model. We choose the prior to be uninformative over the possible region for $\mathbf{\beta }$ , the estimation is dominated by the likelihood, and so the results are robust to the choice of prior. A constrained uniform prior on the parameter $\mathbf{\beta }$ for the possible region goes as follows.

1. GP integer-valued transfer function: p $(\mathbf{\beta })I(A_1)$ , where $A_1=\{{\alpha }_0>0,:{\alpha }_1,{\beta }_1\ge 0,:{\alpha }_1+{\beta }_1<1\}$ .

2. Log-linear GP integer-valued transfer function: p $(\mathbf{\beta })I(A_2)$ , where $A_2$ is in Eq. (2.13).

3. NB integer-valued transfer function: p $(\mathbf{\beta })I(A_3)$ , where $A_3$ is in Eq. (2.17).

We now discuss the prior set-up for the second group, $({\omega }_0,\delta )$ . For ${\omega }_0$ , we employ a gamma prior when we have a prior belief ${\omega }_0$ is positive. In other words, ${\omega }_0\sim G(b_1,b_2)$ , for GP and NB integer-valued transfer function models, where $b_1$ and $b_2$ are the shape and rate parameters, respectively. We employ a flat and non-informative prior for ${\omega }_0$ in the log-linear GP integer-valued transfer function model since there is no positive restriction in this case. For the parameter $\delta$ , we assume a beta prior for $\delta$ , $\mathit{Beta}(a_1,a_2)$ . Regarding the parameter $r$ in the NB distribution, we choose a truncated gamma distribution, $r\sim G(k_1,k_2)I(r>1)$ , to be greater than 1 on $r$ . For the parameter $\psi$ in the GP distribution, we assume $\psi \sim \mathit{Unif}(0,1)$ .

When the location of $T$ is unknown, there are several choices for the prior of $T$ in the literature. (Carlin et al. (1992)) manage each point as a potential change point and choose a uniform prior from 1 to n when they deal with the British coal mining disasters dataset. (Chib (1998)) treats each point as a possible outlier and chooses the point with the highest posterior probability as the change point of the series data, which (Gamerman and Lopes (2006)) later adopt. Alternatively, (Chen et al. (2011)) employ a continuous but constrained uniform prior for the change point and subsequently discretize the estimates so that they become an actual time index. (Chen and Lee (2021)) utilize a similar prior setting, take data-dependent choices of the hyperparameters on the multiple breakpoint parameters and afterwards discretize the estimates to be the actual time indices. (Raftery (1996)), (Richardson and Green (1997)) and (Wasserman (2000)) propose data-dependent priors to achieve weakly informative priors. As people consider increasingly complicated problems, the demand for data-dependent priors will become more common; see (Wasserman (2000)). We follow the method of (Chen et al. (2011)) to employ a continuous and data-dependent prior on time $T$ ; that is, $T\sim \mathit{Unif}({\mathcal{P}}_a,{\mathcal{P}}_{1-a})$ , where ${\mathcal{P}}_a$ is the $a$ th percentile of the time index.

For a multiple transfer function model, without loss of generality, we consider two interventions in the integer-valued transfer function model. We assume the prior of ( $T_1$ , $T_2$ ) as follows.

$$T_1\sim \mathit{Unif}({\mathcal{P}}_a,{\mathcal{P}}_b),T_2|T_1\sim \mathit{Unif}({\mathcal{P}}_c,{\mathcal{P}}_d),$$

(3.3)

where ${\mathcal{P}}_a$ , ${\mathcal{P}}_b$ and ${\mathcal{P}}_d$ are respectively the $h$ th, $(1-2h)$ th and $(1-h)$ th percentiles of the set of integers $\{1,2,\dots ,n\}$ . Moreover, ${\mathcal{P}}_c=T_1+\mathit{hn}$ , and this set-up assures at least $h\%$ of observations are in-between $T_1$ and $T_2$ . Accordingly, our priors for $T_1$ and $T_2$ are flat over the region and guarantee that $T_1<T_2$ and at least $h\%$ of observations are in-between $T_1$ and $T_2$ .

3.1 Posterior distributions

By the Bayes theorem, the conditional posterior distribution is proportional to the likelihood function multiplied by the prior density for each group. The conditional posterior distribution for parameter group ${\theta }_{\mathit{j}}$ is as follows:

$$\begin{array}{ccc}\hfill p({\mathbf{\theta }}_j|\mathit{X},{\mathbf{\theta }}_{\setminus j})& \hfill \hfill & \mathcal{L}(\mathit{X}|\mathbf{\theta })p({\theta }_{\mathit{j}}|{\theta }_{\setminus \mathit{j}}),\hfill \end{array}$$

(3.4)

where $p({\theta }_{\mathit{j}})$ is its prior density, and ${\theta }_{\setminus \mathit{j}}$ denotes the vector of $\mathbf{\theta }$ without the component of ${\theta }_{\mathit{j}}$ . Since the conditional posterior distributions are non-standard forms, we employ the random-walk Metropolis algorithm for parameters $\psi$ (or $r$ ) and $T$ . As we assign a continuous but constrained uniform prior on time $T$ , afterwards we discretize the estimates so that they become an actual time index. To speed up the convergence, we follow the adaptive MCMC algorithm of (Chen and So (2006)) to sample MCMC iterates for the two parameter groups, $\mathbf{\beta }$ and $({\omega }_0,\delta )$ , combining a random-walk Metropolis algorithm and an independent kernel Metropolis-Hastings (MH) algorithm to accelerate convergence and to allow for optimal mixing. We carry out a total of $N$ MCMC iterations, discard the first $M$ as a burn-in period and use the remaining $(N-M)$ iterations for Bayesian inference. We provide the MCMC procedure for the integer-valued transfer function with GP distribution as follows.

• (1)

  Set up initial values of ${\mathbf{\theta }}^{[0]}=({\mathbf{\beta }}^{[0]},{\omega }_0^{[0]},{\delta }^{[0]},{\psi }^{[0]},T^{[0]})$ .

• (2a)

  When $i\le M$ , we apply the random walk Metropolis algorithm for sampling ${\mathbf{\beta }}^{[i]}$ :

• (2b)

  When $i>M$ , we apply the independent kernel MH algorithm for sampling ${\mathbf{\beta }}^{[i]}$ :

• (3a)

  When $i\le M$ , we apply the random walk Metropolis algorithm for sampling $\left({\omega }_0^{[i]},{\delta }^{[i]}\right)$ .

• (3b)

  When $i>M$ , we apply the independent kernel MH algorithm for sampling $\left({\omega }_0^{[i]},{\delta }^{[i]}\right)$ .

• (4)

  Sample ${\psi }^{[i]}$ from $P(\psi |\mathit{X},{\mathbf{\beta }}^{[i]},{\omega }_0^{[i]},{\delta }^{[i]},T^{[i-1]})$ .

• (5)

  Sample $T^{[i]}$ from $P(T|\mathit{X},{\mathbf{\beta }}^{[i]},{\omega }_0^{[i]},{\delta }^{[i]},{\psi }^{[i]})$ , and we discretize $T^{[i]}$ as an integer value.

• (6)

  Go back to Step 2.

Note that we choose $c$ in Step (2a)(i) in order to control the updating (acceptance) rate in-between 25% to 50%.

Model selection and diagnostic checking both play crucial roles in empirical applications. To select the best-fitted model among all competing models, we employ DIC, POR and MSSR. First, we express DIC as follows.

$$\begin{array}{ccc}\hfill \mathit{DIC}& \hfill =\hfill & \overline{D(\mathbf{\theta })}+P_D,p_D=\overline{D(\mathbf{\theta })}-D(\bar{\mathbf{\theta }}),\hfill \\ \hfill D(\mathbf{\theta })& \hfill =\hfill & -2\text{ln}\mathcal{L}(\mathit{X}|\mathbf{\theta })+2\text{ln}p(\mathit{X}),\hfill \end{array}$$

where $\overline{D(\mathbf{\theta })}=E[D(\mathbf{\theta })]$ , and $\bar{\mathbf{\theta }}$ is the posterior mean (or posterior mode for the locations of interventions) of $\mathbf{\theta }$ ; $\text{ln}\mathcal{L}(\mathit{X}|\mathbf{\theta })$ is the log-likelihood function in (3.1) or (3.2). We obtain DIC as a by-product of the MCMC procedure.

Bayesian hypothesis testing can represent a special case of model comparison where a model refers to a likelihood function and a prior distribution. We offer a comparison of ${\mathcal{H}}_0:M_1$ versus ${\mathcal{H}}_1:M_2$ by POR:

$$\mathit{POR}=\frac{p(M_1|\mathit{X})}{p(M_2|\mathit{X})}=\frac{p(\mathit{X}|M_1)}{p(\mathit{X}|M_2)}\times \frac{\mathit{Pr}(M_1)}{\mathit{Pr}(M_2)},$$

(3.5)

where $p(\mathit{X}|M_i)$ denotes the marginal likelihood under model $M_i$ , and $\mathit{Pr}(M_i)$ is the prior probability of $M_i$ , for which we can choose an equal weight to present no preference for each model. We express the posterior probability of ${\mathcal{H}}_i$ based on the idea of Berger and Delampady (1987) as follows:

$$\mathit{Pr}({\mathcal{H}}_i|\mathit{X})={\left[1+\frac{p(\mathit{X}|M_j)}{p(\mathit{X}|M_i)}\times \frac{\mathit{Pr}(M_j)}{\mathit{Pr}(M_i)}\right]}^{-1}={\left[1+\frac{1}{\mathit{POR}}\right]}^{-1},$$

(3.6)

where $p(\mathit{X}|M_i)$ denotes the marginal likelihood under the null hypothesis ${\mathcal{H}}_i$ .

We favour the model with the smallest DIC value. To check the adequacy of the model, we calculate the standardized Pearson residuals proposed by (Jung et al. (2006)):

$$z_t=\frac{X_t-E(X_t|{\mathit{X}}_{t-1})}{\sqrt{\mathit{Var}(X_t|{\mathit{X}}_{t-1})}}.$$

(3.7)

If the model is correctly specified, then these residuals should have mean zero and variance one with no significant serial correlation in $z_t$ and $z_t^2$ . We also employ MSSR for model comparison, defined below as follows:

$$\begin{array}{ccc}\hfill \mathit{MSSR}& \hfill =\hfill & \frac{1}{n-1}\sum _{t=2}^n{\tilde{z}}_t^2\hfill \\ \hfill {\tilde{z}}_t^2& \hfill =\hfill & \frac{1}{N-M}\sum _{i=M+1}^N{z}_t^2\approx \frac{1}{N-M}\sum _{i=M+1}^N{\left(\frac{X_t-E(X_t|{\mathit{X}}_{t-1},{\mathbf{\theta }}^{[i]})}{\sqrt{\mathit{Var}(X_t|{\mathit{X}}_{t-1},{\mathbf{\theta }}^{[i]}})}\right)}^2,\hfill \end{array}$$

(3.8)

where the last equation is a simulation-based estimate, and ${\mathbf{\theta }}^{[i]}$ is the $i$ th MCMC iterate.

4 Simulation study

We conduct a simulation study to examine the effectiveness of MCMC procedures for the following scenarios: (a) integer-valued transfer functions with GP or NB distributions; (b) one or multiple interventions with known or unknown locations; (c) the four types of interventions based on Inv-2A, Inv-2B, Inv-4A and Inv-4B; and (d) model misspecification.

Table 1

Simulation results for integer-valued transfer function models with known location obtained from 200 replications

Log-Linear GP Integer-Valued Transfer Function
Intervention	Par.	True	Mean	Median	Std.	P 2.5	P 97.5
Inv-2A	α 0	1.0	1.1352	1.1327	0.1453	0.8578	1.4273
	α 1	0.3	0.3079	0.3078	0.0301	0.2494	0.3669
	β 1	0.3	0.2406	0.2411	0.0637	0.1141	0.3639
	ψ	0.18	0.1802	0.1803	0.0188	0.1431	0.2169
	ω 0	1.0	0.9745	0.9814	0.1779	0.6052	1.3029
Inv-2B	α 0	1.0	1.0896	1.0879	0.1186	0.8620	1.3273
	α 1	0.3	0.3243	0.3243	0.0302	0.2655	0.3837
	β 1	0.2	0.1294	0.1300	0.0507	0.0280	0.2270
	ψ	0.3	0.2987	0.2987	0.0178	0.2635	0.3335
	ω 0	1.0	0.9949	0.9992	0.2014	0.5879	1.3782
	δ	0.6	0.6020	0.6112	0.1116	0.3594	0.7936
Inv-4A	α 0	0.5	0.4830	0.4822	0.0469	0.3931	0.5763
	α 1	0.3	0.2990	0.2989	0.0314	0.2380	0.3610
	β 1	0.3	0.3104	0.3115	0.0426	0.2247	0.3910
	ψ	0.3	0.2987	0.2988	0.0202	0.2590	0.3383
	ω 0	0.3	0.2948	0.2942	0.0036	0.2909	0.3007
Inv-4B	α 0	0.5	0.5536	0.5523	0.0805	0.3997	0.7154
	α 1	0.3	0.3003	0.3002	0.0342	0.2337	0.3674
	β 1	0.2	0.1539	0.1543	0.0702	0.0148	0.2905
	ψ	0.18	0.1786	0.1786	0.0243	0.1311	0.2265
	ω 0	0.3	0.3445	0.3440	0.0587	0.2315	0.4604
	δ	0.3	0.2535	0.2494	0.1010	0.0711	0.4612

4.1 Intervention effect (known location)

We organize simulation studies by using a sample of size $n=500$ with 200 replications. For each simulated dataset, we use 20, 000 MCMC iterations and discard 8, 000 burn-in iterations. We set the initial values for unknown parameters: $({\alpha }_0,{\alpha }_1,{\beta }_1)=(0.2,0.1,0.1)$ , $({\omega }_0,\delta ,\psi ,r)=(0.1,0.1,0.1,1.1)$ ; and choose the hyperparameters to be $G(3,1)$ for ${\omega }_0$ , Beta(2,2) for $\delta$ and $G(7,1)I(r>1)$ for $r$ in the NB distribution. We further consider a multiple-intervention model with sample size 700, which is as follows:

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathcal{D}(.)\hfill \\ \hfill & \hfill \hfill & {\kappa }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\kappa }_{t-1}+{\omega }_{01}{P}_{1,t}+\frac{{\omega }_{02}}{(1-\delta B)}{P}_{2,t},\hfill \end{array}$$

(4.1)

where D denotes a distribution, and $P_{i,t}$ is defined in Eq. (2.3). Due to limited space, we only list partial results in tables of this article. Concerning the MCMC diagnostic, we observe that all ACF plots decay very quickly for all estimates, indicating a fair mixing of MCMC iterates for each parameter. Furthermore, the traceplots exhibit, on average, that the MCMC iterates after the burn-in period are close to the true parameter values. The true parameter values for the entertained models with interventions are described in Tables 1–3, which display the summary statistics from 200 replications for the average of posterior means, medians, standard deviations and 95% credible intervals. We summarize the findings as follows.

Table 2

Simulation results for multiple interventions obtained from 200 replications

D	Par.	True	Mean	Median	Std.	P 2.5	P 97.5
1) GP Inv-2A 2) GP Inv-2B	α 0	1	1.0483	1.0460	0.1183	0.8231	1.2865
	α 1	0.3	0.2942	0.2941	0.0225	0.2504	0.3387
	β 1	0.3	0.2944	0.2948	0.0491	0.1970	0.3894
	ψ	0.18	0.1759	0.1757	0.0203	0.1364	0.2158
	ω 01	10	9.9944	9.9542	1.7584	6.6709	13.5534
	ω 02	8	7.7511	7.7173	1.4199	5.0780	10.6262
	δ 1	0.5	0.4852	0.4860	0.0731	0.3405	0.6258
1) Log GP Inv-2A 2) Log GP Inv-2B	α 0	1.0	1.1809	1.1797	0.1061	0.9772	1.3923
	α 1	0.3	0.3119	0.3118	0.0249	0.2633	0.3609
	β 1	0.2	0.1053	0.1058	0.0484	0.0091	0.1988
	ψ	0.3	0.3005	0.3005	0.0150	0.2708	0.3299
	ω 01	1.5	1.4597	1.4698	0.1947	1.0461	1.8117
	ω 02	1.0	0.9975	1.0020	0.1991	0.5952	1.3749
	δ 1	0.6	0.6084	0.6176	0.1106	0.3679	0.7978
1) NB Inv-2A 2) NB Inv-2B	α 0	0.50	0.5107	0.5094	0.0557	0.4050	0.6233
	α 1	0.80	0.0764	0.0763	0.0059	0.0651	0.0882
	β 1	0.30	0.2854	0.2856	0.0437	0.1991	0.3704
	r	5.00	5.2168	5.2135	0.1707	4.8879	5.5535
	ω 01	6.00	6.7328	6.6666	1.4407	4.1056	9.7258
	ω 02	5.00	4.9608	4.9081	1.1260	2.9096	7.3024
	δ 1	0.50	0.4824	0.4828	0.0699	0.3446	0.6179

The two interventions, P_1t and P_2t in Eq. (4.1), are known locations and are given as follows:

$P_{1t}=\{\begin{array}{c}1, if t=230\\ 0, if t\ne 230\end{array},$

$P_{1t}=\{\begin{array}{c}2, if t=450\\ 0, if t\ne 450\end{array},$

1. For different values of ${\omega }_0$ and $\psi$ for Inv-2A and Inv-4A, the MCMC estimation results are sound, in which the posterior mean and median estimates are close to the true values. A similar situation applies to ${\omega }_0$ , $\psi$ and $\delta$ for Inv-2B and Inv-4B models.

2. The results show a successful estimation for single and multiple intervention models (Tables 1 and 2). The average posterior means (medians) are overall satisfactorily close to the true values of the parameters, and all true values are within the 95% credible intervals, indicating that the MCMC method is reasonably effective.

3. The estimate of $\psi$ in a GP distribution is close to 0, when we misspecify the true Poisson distribution to be a GP distribution. The estimate of $r$ is around 5.8 in Table 3 when we misspecify the true distribution $\mathit{GP}({\lambda }_t,0.18)$ to be a NB distribution. The estimate of $\psi \cong 0.28$ in the GP distribution when we misspecify the true NB $(7,p_t)$ distribution to be a GP distribution.

Table 3

Simulation results of the mis-specified models obtained from 200 replications

		True Model
		Inv-2B & GP					Inv-2B & NB
Fitted Model	Par.	True	Mean	P 2.5	P 97.5	Fitted Model	True	Std.	P 2.5	P 97.5
GP Inv-2B	α 0	1.0	1.0937	0.8087	1.3948	GP Inv-2B		1.6818	1.2896	2.0945
	α 1	0.3	0.3005	0.2483	0.3534			0.2411	0.1961	0.2867
	β 1	0.3	0.2704	0.1498	0.3889			0.2737	0.1581	0.3871
	ψ	0.18	0.1763	0.1297	0.2237			0.2840	0.2429	0.3255
	ω 0	10.0	9.6974	5.8622	14.3607			16.0218	10.2482	22.6612
	δ	0.6	0.5951	0.3310	0.7851			0.5969	0.3769	0.7614
Log GP Inv-2B	α 0		0.2394	0.1135	0.3692	Log GP Inv-2B		0.5025	0.3561	0.6538
	α 1		0.4618	0.3944	0.5299			0.3832	0.3237	0.4431
	β 1		0.1639	0.0603	0.2634			0.1685	0.0685	0.2649
	ψ		0.1821	0.1356	0.2294			0.2797	0.2388	0.3207
	ω 0		1.3957	0.8956	1.8864			1.7874	1.3843	2.1610
	δ		0.6549	0.4347	0.8246			0.6380	0.4768	0.7754
NB Inv-2B	α 0		0.2518	0.1842	0.3259	NB Inv-2B	0.3	0.3405	0.2541	0.4340
	α 1		0.0660	0.0518	0.0814		0.05	0.0521	0.0414	0.0636
	β 1		0.2292	0.1132	0.3448		0.3	0.2616	0.1457	0.3766
	r		5.8170	5.0304	6.7068		7.0	6.8805	6.1306	7.6961
	ω 0		5.8686	3.4383	8.7209		5.0	4.7184	2.8908	6.8293
	δ		0.4012	0.1760	0.6141		0.5	0.5220	0.3122	0.7019

4.2 Transfer function (unknown location)

We treat the location of $T$ as an unknown parameter and choose an initial value $T=101$ , which is far away from the true value $T=290$ . Since $T$ is an integer value, we report Bayesian maximum a posteriori (MAP) estimate for the unknown location of $T$ . The MAP estimate corresponds to the most probable value of the posterior distribution. As in Tables 4–7, the average posterior means (medians) are fairly close to the true values when the underlined models are correct. Regarding the multiple effects, we pick up the initial values for $(T_1^{[0]},T_2^{[0]})=(140,280)$ , which are some distance away from the true values, (230, 450). The MCMC method is reasonably effective even when we consider the complicated multiple cases in Table 6. The MCMC method with continuous data-dependent priors of $T_1$ and $T_2$ provides a soundly effective estimate for unknown locations. The simulation results show that the average of MAP estimates of $T$ is close to the true value for the proposed transfer function models, even though we misspecify the distribution in Table 7.

Table 4

Simulation results for integer-valued transfer function models with unknown location of T obtained from 200 replications Intervention GP Integer-Valued Transfer Function

Intervention	GP Integer-Valued Transfer Function
	Par.	True	Mean	Median	Std.	P 2.5	P 97.5
Inv-2A	α 0	1.0	1.0899	1.0872	0.1341	0.8347	1.3593
	α 1	0.3	0.3022	0.3021	0.0231	0.2573	0.3478
	β 1	0.3	0.2744	0.2748	0.0501	0.1753	0.3715
	ψ	0.3	0.2978	0.2977	0.0213	0.2564	0.3397
	ω 0	20	20.0462	20.2344	4.5726	10.7151	28.1121
	T ∗	290	289.630
Inv-2B	α 0	2.0	1.9510	1.9469	0.2417	1.4893	2.4367
	α 1	0.3	0.3023	0.3023	0.0227	0.2580	0.3471
	β 1	0.3	0.3070	0.3073	0.0482	0.2115	0.4005
	ψ	0.3	0.3012	0.3012	0.0187	0.2646	0.3379
	ω 0	15	15.4053	15.2379	3.4663	9.1493	22.6685
	δ	0.6	0.5954	0.6085	0.1973	0.1979	0.9282
	T ∗	290	289.415
Inv-4A	α 0	2.0	2.2087	2.2023	0.3028	1.6331	2.8189
	α 1	0.3	0.3012	0.3012	0.0226	0.2573	0.3456
	β 1	0.3	0.2734	0.2738	0.0525	0.1692	0.3751
	ψ	0.3	0.2969	0.2968	0.0197	0.2582	0.3354
	ω 0	2.0	2.2046	2.1926	0.4136	1.4293	3.0482
	T ∗	290	288.5
Inv-4B	α 0	0.5	0.5049	0.5038	0.0564	0.3978	0.6185
	α 1	0.3	0.2991	0.2991	0.0225	0.2553	0.3433
	β 1	0.3	0.2958	0.2962	0.0436	0.2086	0.3801
	ψ	0.3	0.3008	0.3007	0.0228	0.2562	0.3454
	ω 0	2.0	2.0015	1.9753	0.5399	1.0299	3.1242
	δ	0.6	0.5980	0.6033	0.1098	0.3699	0.7952
	T ∗	290	289.765

Note: ∗The average of the MAP estimates of T.

Table 5

Simulation results for integer-valued transfer function models with unknown location of T obtained from 200 replications

Log GP Integer-Valued Transfer Function
Intervention	Par.	True	Mean	Median	Std.	P 2.5	P 97.5
Inv-2A	α 0	1.0	1.1105	1.1077	0.1449	0.8350	1.4026
	α 1	0.3	0.3064	0.3062	0.0299	0.2482	0.3654
	β 1	0.3	0.2513	0.2519	0.0624	0.1273	0.3715
	ψ	0.3	0.3036	0.3037	0.0166	0.2709	0.3360
	ω 0	1.0	1.1848	1.1979	0.1812	0.7935	1.5049
	T ∗	290	290
Inv-2B	α 0	1.0	1.0732	1.0717	0.1064	0.8692	1.2858
	α 1	0.3	0.3174	0.3173	0.0280	0.2628	0.3724
	β 1	0.2	0.1454	0.1461	0.0552	0.0352	0.2513
	ψ	0.3	0.3028	0.3029	0.0181	0.2672	0.3381
	ω 0	1.0	1.1697	1.1784	0.1822	0.7892	1.5033
	δ	0.6	0.6082	0.6118	0.1479	0.3108	0.8828
	T ∗	290	289.30
Inv-4A	α 0	0.5	0.4539	0.4541	0.0338	0.3875	0.5202
	α 1	0.3	0.3132	0.3128	0.0331	0.2491	0.3789
	β 1	0.3	0.3050	0.3054	0.0365	0.2324	0.3752
	ψ	0.18	0.1843	0.1844	0.0234	0.1386	0.2301
	ω 0	0.5	0.4899	0.4908	0.0062	0.4798	0.4975
	T ∗	290	289.75
Inv-4B	α 0	0.5	0.4734	0.4731	0.0546	0.3673	0.5813
	α 1	0.3	0.3006	0.3005	0.0325	0.2376	0.3648
	β 1	0.25	0.2648	0.2651	0.0456	0.1744	0.3539
	ψ	0.18	0.1820	0.1820	0.0239	0.1353	0.2290
	ω 0	0.3	0.3143	0.3148	0.0497	0.2164	0.4097
	δ	0.3	0.2482	0.2434	0.1021	0.0662	0.4598
	T ∗	290	290.29

Note: ∗The average of the MAP estimates of T .

Table 6

Simulation study of multiple interventions for unknown locations of T₁ and T₂

Interventions	Par.	True	Mean	Median	Std.	P 2.5	P 97.5
1) GP Inv-2A 2) GP Inv-2B	α 0	1	1.0856	1.0834	0.1221	0.8529	1.3308
	α 1	0.3	0.2994	0.2977	0.0228	0.2533	0.3426
	β 1	0.3	0.2809	0.2814	0.0506	0.1806	0.3784
	ψ	0.18	0.1786	0.1784	0.0208	0.1383	0.2201
	ω01	10	10.4126	10.4002	2.0313	6.4619	14.4448
	ω02	8	7.8227	7.8125	1.6791	4.5762	11.1461
	δ 1	0.5	0.5050	0.5068	0.0821	0.3403	0.6600
	T 1	230	230.26	230
	T 2	450	449.615	450
1) Log GP Inv-2A 2) Log GP Inv-2B	α 0	1.0	1.1534	1.1522	0.1063	0.9484	1.3651
	α 1	0.4	0.4106	0.4105	0.0267	0.3585	0.4633
	β 1	0.2	0.1321	0.1324	0.0457	0.0417	0.2211
	ψ	0.3	0.3045	0.3045	0.0142	0.2765	0.3324
	ω 0	1.0	0.7431	0.7475	0.1728	0.4021	1.0662
	ω 1	1.0	1.0049	1.0054	0.1048	0.7980	1.2085
	δ 1	0.6	0.6149	0.6169	0.0630	0.4861	0.7322
	T*1	230	227.61	230.0
	T*2	450	449.66	450.0
1) NB Inv-2A 2) NB Inv-2B	α 0	0.5	0.5237	0.5226	0.0588	0.4123	0.6421
	α 1	0.08	0.0782	0.0781	0.0063	0.0662	0.0908
	β 1	0.3	0.2949	0.2951	0.0451	0.2062	0.3824
	r	5	5.0282	5.0260	0.1682	4.7098	5.3622
	ω 1	6	6.7416	6.6803	1.4940	3.9965	9.8334
	ω2	5	4.5337	4.4865	1.1472	2.4229	6.9073
	δ	0.5	0.4970	0.4979	0.0807	0.3375	0.6519
	T*1	230	229.70	230
	T ∗2	450	450.03	450

Note: ∗The third and fourth cells provide the average and median of the MAP estimates of T _i , respectively.

Table 7

Simulation results of the misspecified distribution with unknown location of T obtained from 200

		True Model
		Inv-2B & GP					Inv-2B & NB
Fitted Model	Par.	True	Mean	P 2.5	P 97.5	Fitted Model	True	Std.	P 2.5	P 97.5
GP Inv-2B	α 0	2.0	2.2065	1.6894	2.7499	GP Inv-2B		1.9877	1.5285	2.408
	α 1	0.3	0.3022	0.2584	0.3465			0.0895	0.0438	0.1366
	β 1	0.3	0.2706	0.1715	0.3680			0.1681	0.0262	0.3361
	ψ	0.3	0.2991	0.2621	0.3361			0.2187	0.1727	0.2656
	ω 0	10	9.2788	4.3860	15.1423			14.2996	8.7787	20.8479
	δ	0.6	0.5365	0.1940	0.8134			0.4464	0.2048	0.6563
	T ∗	290	288.55					289.16
Log GP Inv-2B	α 0		0.4435	0.2949	0.5972	Log GP Inv-2B		0.6146	0.4553	0.7867
	α 1		0.4726	0.4123	0.5342			0.1779	0.1133	0.2426
	β 1		0.2121	0.1289	0.2928			0.1554	0.0085	0.2894
	ψ		0.3037	0.2666	0.3404			0.2192	0.1728	0.2661
	ω 0		1.1236	0.6612	1.5203			2.3072	1.9664	2.6238
	δ		0.6062	0.3051	0.8844			0.6075	0.3070	0.8846
	T ∗		289.31					289.195
NB Inv-2B	α 0		0.3928	0.2856	0.5105	NB Inv-2B	0.4	0.4471	0.3393	0.5590
	α 1		0.0520	0.0425	0.0623		0.02	0.0206	0.0097	0.0323
	β 1		0.2461	0.1376	0.3531		0.2	0.1454	0.0201	0.2990
	r		8.5456	7.5609	9.5953		6.0	5.9365	5.0070	6.9937
	ω 0		5.6787	3.3835	8.3362		5.0	5.3402	3.2454	7.7821
	δ		0.3040	0.0908	0.5353		0.3	0.3561	0.1679	0.5376
	T ∗		288.265				290	287.375

Note: ∗The average of the MAP estimates of T.

When we fit an integer-valued transfer function model to a set of data, the true order of the underlying process is often unknown; hence, it is likely to be misspecified. Table 8 now scrutinizes the effects on estimation in misspecified autoregressive order for the following model.

Table 8

Simulation results of the misspecified autoregressive order with unknown location of T obtained from 200 replications

True Model in Eq. (4.2)							True Model in Eq. (4.2)
Fitted Model	Par.	True	Mean	P 2.5	P 97.5	Fitted Model	True	Std.	P 2.5	P 97.5
GP Inv-2A	α 0	1.0	1.2230	0.8377	1.6373	GP Inv-2B	1.0	1.2746	0.8796	1.7029
	α 1	0.3	0.2820	0.2331	0.3322		0.3	0.2878	0.2392	0.3374
	β 1	0.3	0.4915	0.3977	0.5813		0.3	0.4798	0.3860	0.5689
	β2	0.2					0.2
	ψ	0.18	0.1822	0.1411	0.2237		0.18	0.1809	0.1398	0.2227
	ω 0	10.0	9.2569	4.5804	14.8425		10	9.1740	4.8718	14.4429
	δ						0.6	0.5197	0.1922	0.7976
	T∗	290	281.95				290	289.405
Log GP Inv-4A	α 0	0.3	0.2671	0.1991	0.3342	Log GP Inv-4B	0.3	0.3060	0.2250	0.3870
	α 1	0.1	0.0988	0.0346	0.1646		0.1	0.0962	0.0291	0.1641
	β 1	0.2	0.3634	0.2928	0.4302		0.2	0.2965	0.2043	0.3877
	β 2	0.1					0.1
	ψ	0.18	0.1780	0.1292	0.2275		0.18	0.1801	0.1301	0.2305
	ω 0	0.9	0.8946	0.8811	0.9129		0.7	0.7483	0.5421	0.9454
	δ						0.3	0.2552	0.0815	0.4448
	T ∗	290	289.885				290	290.045
NB Inv-4A	α 0	0.3	0.3532	0.2602	0.4501	NB Inv-4B	0.3	0.3256	0.2575	0.3984
	α 1	0.05	0.0514	0.0387	0.0649		0.05	0.0502	0.0405	0.0605
	β 1	0.2	0.2752	0.1539	0.4014		0.2	0.2709	0.1835	0.3592
	β2	0.1					0.1
	r	6.0	5.5996	5.1009	6.1138		7.0	7.0429	6.4562	7.6571
	ω 0	1.0	1.2088	0.8629	1.5907		1.0	1.0272	0.4878	1.6531
	δ						0.5	0.5202	0.2371	0.7700
	T ∗	290	290.145				290	289.995

Note: ∗Average of the MAP estimates of T.

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{GP}({\kappa }_t,\psi ),\hfill \\ \hfill & \hfill \hfill & {\kappa }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\kappa }_{t-1}+{\beta }_2{\kappa }_{t-2}+\frac{{\omega }_s(B)}{{\delta }_q(B)}{I}_t,\hfill \end{array}$$

(4.2)

where $I_t$ is either $P_t$ or $S_t$ . The estimate of ${\beta }_1$ tends to have a positive bias when we ignore the high-order term ${\beta }_2{\kappa }_{t-2}$ in (4.2). In other words, underfitting the autoregressive order always comes with some cost, in which the estimate ${\beta }_1$ is larger than the true value.

Figure 2

Time plots of Data 1, Data 2 and Data 3 from January 1995 to March 2020

5 Real examples

To demonstrate the proposed methodology, we use monthly crime counts of Albury City in NSW Australia from the NSW Bureau of Crime Statistics and Research (https://www.bocsar.nsw.gov.au/) covering from January 1995 to March 2020. Specifically, we use the number of against justice procedures (breach apprehended violence order), sexual offenses (indecent assault, act of indecency and other sexual offenses), and disorderly conduct (trespassing) and name them Data 1, 2 and 3, respectively. Figure 2 exhibits the time plot of Data 1–3. Table 9 shows the descriptive statistics of Data 1–3. We observe that the variance is greater than the mean in these datasets, thus presenting overdispersion. From a visual inspection of Data 1, we notice an intervention at $T=77$ (May 2001).

Table 9

Descriptive statistics for data examples

Data	Crime Category	Mean	Variance	Min	Q1	Median	Q3	Max
Data 1	Against justice procedure	11.25	40.8789	0	7	10	15	46
Data 2	Sexual offense	4.908	8.806	0	3	5	7	24
Data 3	Disorderly conduct	6.835	13.6813	0	4	7	9	20

Data 2 lists monthly counts of sexual offenses, and we specify that the intervention starts at $T=265$ , which is January 2017. One possible reason for this abrupt increase of sexual assaults is the relaxation of the Liquor Amendment Act 2014 announced by the Premier of NSW in December 2016. According to the 2016 Australian Bureau of Statistics (ABS) personal safety survey, $50\%$ (or 321, 000) of women believe that alcohol or another substance contributed to their most recent incident of sexual assault. For Data 3, we estimate the intervention locations later. We check the convergence diagnostics of the MCMC iterates by providing the trace and ACF plots according to the fitted model for each dataset. Due to space limitation, we do not provide any diagnostic plots here. We conclude that all ACF plots die down quickly, indicating no autocorrelation, and that all MCMC samples seem to mix well.

Table 10 presents results of model selection for Data 1 and 2 when the locations of interventions are known. Results show that the NB integer-valued transfer function model (Inv-1) and NB integer-valued transfer function model (Inv-2A) are the favoured models for Data 1 and 2, respectively, in terms of the lowest DIC. Bayesian estimation for Data 1 and 2 when time $T$ is unknown appears in Table 11. We use the MAP estimate for the unknown location of $T$ . For Data 1, all models indicate the intervention happens at $T=77$ (May 2001). For Data 2, almost all models are in agreement that intervention occurs at $T=265$ , except for the GP cases, which detect an intervention at $T=270$ . Table 12 summarizes the DIC and MSSR results for all candidate models with an unknown location of intervention. We prefer the lowest DIC and MSSR. The favoured models are in agreement with the favoured models with the known location of $T$ .

Table 10

MCMC results for Data 1–2 for the known location of intervention

Fitted Model		α 0	α 1	β 1	ω 0	δ	ψ	r	DIC
Data 1	T = 77
GP	Inv-1	2.3756	0.3851	0.1003	11.3553		0.3369		1789.946
GP	Inv-2	2.5328	0.3676	0.1023	9.3087	0.2924	0.3398		1793.695
GP	Inv-2A	2.4991	0.3746	0.0997	10.6478		0.3363		1790.941
GP	Inv-2B	2.5324	0.3693	0.1000	8.8567	0.2878	0.3389		1795.052
Log GP	Inv-1	0.4400	0.5659	0.0989	1.4473		0.3228		1783.632
Log GP	Inv-2	0.4646	0.5527	0.0996	1.3043	0.2037	0.3276		1794.579
Log GP	Inv-2A	0.4566	0.5573	0.1001	1.4010		0.3246		1784.274
Log GP	Inv-2B	0.4370	0.5600	0.1000	1.1732	0.1986	0.3335		1792.862
NB	Inv-1	0.0941	0.0416	0.6040	5.6916			8.0163	1769.389
NB	Inv-2	0.1885	0.0528	0.4763	5.0281	0.2639		7.5696	1775.168
NB	Inv-2A	0.1794	0.0562	0.4413	4.2078			7.7961	1779.023
NB	Inv-2B	0.2445	0.0660	0.3974	3.9888	0.2586		6.9208	1785.130
Data 2	T = 265
GP	Inv-1	0.0978	0.5807	0.1033	8.4511		0.3226		1560.893
GP	Inv-2	0.1048	0.5766	0.1087	7.3267	0.3083	0.3212		1560.495
GP	Inv-2A	0.1040	0.5667	0.0987	8.0794		0.3327		1562.543
GP	Inv-2B	0.2096	0.5579	0.1012	7.0760	0.3138	0.3192		1547.861
Log GP	Inv-1	0.0992	0.5982	0.0852	1.4944		0.2647		1459.210
Log GP	Inv-2	0.1001	0.5881	0.1008	1.2895	0.2858	0.2645		1460.649
Log GP	Inv-2A	0.0993	0.5924	0.0995	1.3898		0.2631		1453.189
Log GP	Inv-2B	0.1000	0.5852	0.1022	1.2086	0.2815	0.2662		1456.199
NB	Inv-1	0.3696	0.0996	0.1042	5.7815			5.3208	1454.827
NB	Inv-2	0.3835	0.1001	0.1040	4.9690	0.3244		5.2605	1455.702
NB	Inv-2A	0.5304	0.1004	0.0767	5.9518			4.5897	1450.531
NB	Inv-2B	0.3995	0.0993	0.1377	4.7914	0.3233		5.0041	1454.665

Note: Boxes indicate the favoured models.

Table 11

MCMC estimation for the candidate models with unknown location of intervention

Data 1	T = 77
GP	Inv-1	2.4297	0.3779	0.0978	11.0933		0.3404		77
GP	Inv-2	2.4819	0.3693	0.1024	8.8387	0.2924	0.3418		77
GP	Inv-2A	2.5343	0.3728	0.0997	10.7178		0.3357		77
GP	Inv-2B	2.5089	0.3667	0.1002	7.6416	0.3195	0.3450		77
Log GP	Inv-1	0.4380	0.5648	0.1017	1.4438		0.3220		77
Log GP	Inv-2	0.4458	0.5562	0.1009	1.2932	0.2076	0.3323		77
Log GP	Inv-2A	0.4350	0.5638	0.1005	1.4126		0.3271		77
Log GP	Inv-2B	0.3923	0.5740	0.0997	1.0494	0.2264	0.3405		77
NB	Inv-1	0.1154	0.0440	0.5673	6.0444			8.0476	77
NB	Inv-2	0.2238	0.0551	0.4420	5.1032	0.2672		7.4790	77
NB	Inv-2A	0.1686	0.0561	0.4497	4.1971			7.7945	77
NB	Inv-2B	0.2190	0.0614	0.3946	3.8008	0.2605		7.5282	77
Data 2
GP	Inv-1	0.0971	0.5809	0.1030	6.5194		0.3254		270
GP	Inv-2	0.1021	0.5769	0.1031	5.6062	0.3426	0.3272		270
GP	Inv-2A	0.0958	0.5733	0.1065	5.9332		0.3303		270
GP	Inv-2B	0.2290	0.5517	0.0988	5.3719	0.3450	0.3244		270
Log GP	Inv-1	0.1013	0.5850	0.0967	1.1417		0.2762		265
Log GP	Inv-2	0.1000	0.5774	0.1209	0.9090	0.3429	0.2690		265
Log GP	Inv-2A	0.1001	0.5903	0.1001	0.9159		0.2703		265
Log GP	Inv-2B	0.1003	0.5861	0.1011	0.7095	0.3517	0.2758		265
NB	Inv-1	0.3691	0.0989	0.1514	5.5208			5.0920	265
NB	Inv-2	0.4234	0.1006	0.1196	3.9150	0.3088		4.9435	265
NB	Inv-2A	0.6823	0.1007	0.0601	6.1924			4.0450	265
NB	Inv-2B	0.3720	0.0997	0.1584	4.3086	0.3136		4.9794	265

Note: ∗denotes MAP estimate.

Table 12

Model comparison for the unknown location of intervention

Fitted Model		DICData 1	MSSR	DICData 2	MSSR
GP	Inv-1	1790.665	0.9648	1563.149	1.9416
GP	Inv-2	1794.885	0.9746	1567.980	1.9462
GP	Inv-2A	1790.971	0.9752	1568.164	1.9337
GP	Inv-2B	1802.899	1.0243	1551.923	1.5720
Log GP	Inv-1	1783.303	0.9830	1463.485	1.0323
Log GP	Inv-2	1795.07	0.9708	1467.808	1.0543
Log GP	Inv-2A	1784.501	0.9748	1462.615	1.0489
Log GP	Inv-2B	1795.21	0.9761	1467.470	1.0596
NB	Inv-1	1769.785	0.8530	1457.186	0.8788
NB	Inv-2	1777.176	0.8698	1469.238	0.8822
NB	Inv-2A	1779.238	0.8656	1451.772	0.7089
NB	Inv-2B	1784.348	0.8582	1458.875	0.8709

Note: Boxes indicate the favoured models.

Concerning Data 3, we suspect there are multiple abrupt start and abrupt decay outliers that are sparsely spread throughout the data. We further consider multiple interventions, $M_2$ , as a gradual start and a permanent effect (Inv-4B) at $T_1$ and abrupt decay outlier in ${\kappa }_t$ (Inv-2A) at $T_2$ . When the underlying process is GP, we can express $M_2$ as follows.

$$\begin{array}{ccc}\hfill & \hfill \hfill & X_t|{\mathit{X}}_{t-1}\sim \mathit{GP}\left({\kappa }_t,\psi \right),\hfill \\ \hfill & \hfill \hfill & {\kappa }_t={\alpha }_0+{\alpha }_1{X}_{t-1}+{\beta }_1{\kappa }_{t-1}+\frac{{\omega }_{01}}{(1-\delta B)}{S}_t+{\omega }_{02}{P}_t,\hfill \\ \hfill \mathit{where}& \hfill \hfill & S_t=\left\{\begin{array}{cc}0,\hfill & \mathit{if}t<T_1\hfill \\ 1,\hfill & \mathit{if}t\ge T_1,\hfill \end{array}\right.P_t=\left\{\begin{array}{cc}1,\hfill & \mathit{if}t=T_2\hfill \\ 0,\hfill & \mathit{if}t\ne T_2.\hfill \end{array}\right.\hfill \end{array}$$

(5.1)

We provide parameter estimation and model comparisons for the candidate models in Tables 13 and 14 for Data 3. We detect that the locations of interventions occur at $T_1=121$ (January 2005) and $T_2=214$ (October 2012) built on both NB and GP integer-valued transfer function models. We favour a single intervention ( $M_1$ ) if POR is greater than 1 as $M_1$ is placed in the numerator. All PORs in Table 14 are less than 1, for which the multiple interventions are superior to the single intervention, no matter which distribution. The posterior probability of $H_0:M_1$ is less than 0.5, which also confirms that multiple interventions are preferable. In summary, the GP integer-valued transfer function model in Eq. (5.1) outperforms any candidate models.

Table 13

Bayesian estimates of the multiple intervention model for Data 3

Fitted Model	Par.	Mean	Median	Std.	P 2.5	P 97.5
1) GP Inv-4B2) GP Inv-2A	α 0	0.7898	0.7872	0.1110	0.5752	1.0182
	α 1	0.1683	0.1680	0.0295	0.1109	0.2286
	β 1	0.5713	0.5737	0.0508	0.4671	0.6648
	ψ	0.0989	0.0992	0.0354	0.0307	0.1679
	ω 1	0.6664	0.6626	0.1299	0.4316	0.9332
	δ	0.3481	0.3474	0.0766	0.2013	0.5019
	ω2	5.9879	5.9843	1.0236	4.0710	7.9699
	T*1	121
	T*2	214
1) Log GP Inv-4B2) Log GP Inv-2A	α 0	0.7018	0.7002	0.1082	0.4693	0.9112
	α 1	0.2209	0.2208	0.0360	0.1492	0.2925
	β 1	0.1821	0.1812	0.0701	0.0483	0.3218
	ψ	0.1046	0.1057	0.0341	0.0377	0.1720
	ω 1	0.2976	0.2961	0.0508	0.2085	0.4037
	δ	0.3613	0.3600	0.0806	0.2080	0.5223
	ω2	0.7305	0.7325	0.1893	0.3365	1.0957
	T*1	121
	T*2	214
1) NB Inv-4B2) NB Inv-2A	α 0	0.3098	0.3101	0.0420	0.2288	0.3907
	α 1	0.0217	0.0218	0.0053	0.0114	0.0322
	β 1	0.0502	0.0459	0.0319	0.0034	0.1225
	r	9.3827	9.3895	0.5522	8.2942	10.4811
	ω 1	0.2665	0.2651	0.0361	0.1956	0.3383
	δ	0.2741	0.2732	0.0582	0.1635	0.3933
	ω2	4.3796	4.3572	0.8867	2.8127	6.1676
	T*1	121
	T*2	214

Note: ∗denotes the MAP estimate.

Table 14

Model comparisons for Data 3: The candidate models with unknown locations of interventions

POR
Model		DIC	MSSR	M1 vs M2	Prob(M1|X)
M 1	GP Inv-4B	1470.262	1.0299	0.2811	0.2194
M 2	GP Inv-4B & Inv-2A	1464.734	1.0233
M 1	Log GP Inv-4B	1473.537	1.0265	0.0281	0.0273
M 2	Log GP Inv-4B & Inv-2A	1467.054	1.0219
M 1	NB Inv-4B	1498.633	0.7061	0.0206	0.0202
M 2	NB Inv-4B & Inv-2A	1490.483	0.7671

Figure 3

Diagnostic checking for standardized residuals. Upper panel: Data 1 is based on the NB integer-valued transfer function (Inv-1); lower panel: Data 2 is based on the NB integer-valued transfer function (Inv-2A)

Figure 4

Data 3: The estimated locations of two interventions and diagnostic checking for standardized residuals based on the GP integer-valued transfer function (Inv-4B) and (Inv-2A)

We further examine the residuals based on the lowest DIC value as the best fitted model for each dataset. Figures 3 and 4 present the time plots of the standardized residuals as well as the sample ACFs of the residuals based on the best fitted models. All of the standardized residuals are uncorrelated. On the basis of the diagnostic checking plots, we conclude that the proposed models are adequate. Figure 4 also demonstrates the estimated locations of two interventions for Data 3. It is interesting to note that the average of counts after December 2004 is double than before. Our analysis indicates a gradual start permanent effect beginning at $T_1=121$ (January 2005).

6 Conclusion

This research sets up a transfer function model for time series of counts, specifically focusing on the four types of intervention effects. We herein propose a Bayesian MCMC method based on GP, log-linear GP and NB integer-valued transfer function models to estimate and detect these intervention effects. The adaptive MCMC algorithms give reliable and accurate estimates for all unknown parameters both for known and unknown locations of interventions. The simulation results also indicate that the estimate of location is not sensitive to distribution misspecification. We apply the proposed method to crime datasets and select the favoured models by using three model selection criteria (DIC, MSSR and POR). The empirical outcome also reveals that the proposed models are able to successfully detect the locations and type of interventions. Regarding the GP or NB integer-valued transfer function model, each has its own merits. As a final remark, one can obtain a one-step-ahead prediction, $X_{n+1}$ , based on the proposed MCMC sampling scheme via ${\kappa }_{n+1}$ .