Bayesian modelling of nonlinear negative binomial integer-valued GARCHX models

Abstract

This study focuses on modelling dengue cases in northeastern Thailand through two meteorological covariates: cumulative rainfall and average maximum temperature. We propose two nonlinear integer-valued GARCHX models (Markov switching and threshold specification) with a negative binomial distribution, as they take into account the stylized features of weekly dengue haemorrhagic fever cases, which contain nonlinear dynamics, lagged dependence, overdispersion, consecutive zeros and asymmetric effects of meteorological covariates. We conduct parameter estimation and one-step-ahead forecasting for two proposed models based on Bayesian Markov chain Monte Carlo (MCMC) methods. A simulation study illustrates that the adaptive MCMC sampling scheme performs well. The empirical results offer strong support for the Markov switching integer-valued GARCHX model over its competitors via Bayes factor and deviance information criterion. We also provide one-step-ahead forecasting based on the prediction interval that offers a useful early warning signal of outbreak detection.

Keywords

Dengue fever Integer-valued time series Markov chain Monte Carlo method Markov switching consecutive zeros meteorological covariates

1 Introduction

Dengue is a mosquito-borne viral disease that is widespread throughout the tropics. Female mosquitoes of the species Aedes aegypti are the main vector of dengue, and Aedes albopictus is the secondary vector of dengue virus in Asia (De Melo et al., 2012). Dengue fever causes flu-like symptoms lasting from two days to one week and usually occurs after an incubation period of 4–10 days after infected by mosquito bites. When severe dengue is developing (also known as dengue haemorrhagic fever [DHF]), the critical phase occurs 3–7 days after the first onset of the disease (A detailed clinical guidance can be found on the website of Centers for Disease Control and Prevention [CDC] USA.).

The literature in recent decades has extensively studied meteorological variables such as temperature and rainfall for their potential as early warning tools against dengue, an infectious disease (e.g., Hii et al., 2012). Numerous studies have reported the incidences of dengue in northeastern Thailand, with the tendency for high dengue transmission to rise during the rainy season (Chareonsook et al., 1999; Limkittikul et al., 2014) since the amount of rainfall may offer breeding sites for the mosquitos. This study models weekly DHF cases in Thailand associated with environmental conditions. We suspect that a nonlinear relationship exists between dengue cases and cumulative rainfall/average maximum temperature, which would allow us to look at dengue cases and two meteorological variables rather than just a simple linear relationship.

Many studies focus on analysing the dengue cases associated with meteorological variables via diverse statistical models, for example, regression analysis, autoregressive integrated moving average models generalized linear models and Poisson regression model (Halide and Ridd, 2008; Silawan et al., 2008; Wongkoon et al., 2012; Lowe et al., 2011, 2013).

In practice, such data are usually non-normally distributed and overdispersed (the mean is less than the variance). There are also a large proportion of zeros and consecutive zeros. Several articles model time series of counts using a Poisson distribution. For example, Ferland et al. (2006) propose an integer-valued GARCH (INGARCH) model to deal with overdispersion. Xu et al. (2012) apply a dispersed INARCH model to analyse weekly dengue cases in Singapore. Wang et al. (2014) suggest a self-excited threshold Poisson autoregression. Chen and Lee (2016) introduce the generalized threshold Poisson INGARCH model to cope with overdispersion, dynamic nonlinearity and structural change. Troung et al. (2017) and Chen et al. (2019) further propose some nonlinear Poisson INGARCH models on various time series of counts.

The Poisson model provides the main instrument for describing time series of counts, but other distributional assumptions can be employed instead—the most natural among other candidates being the negative binomial distribution (Fokianos, 2012). Some studies recommend using a negative binomial distribution to cope with overdispersion such as Zhu (2011), Chen et al. (2016) and Chen and Lee (2017). Since we frequently observe dynamic nonlinearity of infectious diseases, we introduce two nonlinear negative binomial INGARCH models (Markov switching and threshold specifications) along with the exogenous covariates in the conditional mean equation (named ‘INGARCHX’ models hereafter) to handle asymmetric effects, overdispersion, consecutive zeros and lagged dependence for the weekly time series of DHF cases. Chen et al. (2019) propose two types of Markov switching (MS) INGARCHX models with the Poisson distribution to analyse weekly DHF cases. The novelty of this study's contribution with respect to that of Chen et al. (2019) is that we incorporate a negative binomial distribution with a greater variance than its Poisson mean to describe time series of counts. The negative binomial with two parameters (shape and scale) makes our proposed model very flexible.

The Markov switching model of Hamilton (1989) and the threshold model of Tong (1978, 1983), known as the regime switching models, are two of the most popular nonlinear time-series models in the literature. The Markov switching models are considered to be dynamic nonlinear with a distinct state of disease incidences, while a threshold model captures the dynamic behaviour of time series by switching the regimes. The general idea for both models is that a process behaves differently when applying a different model. These two classes of models share many common features, such as nonlinear effects. Chen et al. (2011) and Chan et al. (2017) emphasize the similarities and differences for threshold and Markov switching models in the economics and finance domain.

This study conducts parameter estimation and one-step-ahead forecasting for two proposed models based on Bayesian methods. Several authors successfully employ Bayesian approaches to estimate parameters for time series of counts (see Malyshikina et al., 2009; Manitz and Höhle, 2013; Chen and Lee, 2016; Martínez-Bello et al., 2017). Furthermore, we consider a numerical approximation to Bayesian model selection and forecasting, in which we heavily utilize Markov chain Monte Carlo (MCMC) techniques to overcome the difficulty arising in the frequentist method. We provide one-step-ahead forecasting via the predictive distribution, which is a by-product of the sampling algorithm. Our study develops a dengue forecasting model that provides early warning of a dengue outbreak one week in advance via the Bayesian method.

The rest of this article runs as follows. Section 2 introduces negative binomial two-state Markov switching and threshold INGARCHX models, which take into account the Bayesian inference and MCMC methods to evaluate model parameters. We consider model comparison, prediction and diagnostic checking for the proposed models. Section 3 presents simulation studies to evaluate and compare with the two models. Section 4 applies the analytical results and discussion based on the proposed models to the four provinces’ weekly DHF cases. We then make comparisons between competing models. Section 5 gives concluding remarks.

2 Nonlinear negative binomial INGARCH models

We consider two nonlinear specifications: Markov switching and threshold mechanism. The major difference between Markov switching models and threshold models is that the former assume that the underlying state process resulting in the nonlinear dynamics is latent, whereas the latter usually allow the nonlinear dynamics to be driven by observable variables. Both nonlinear specifications are widely applied across diverse fields of endeavour through time-series analysis applicable for capturing regime variation (Hamilton, 1989; Chen et al., 2009).

More precisely, we introduce two nonlinear INGARCHX models with negative binomial variables. One is the MS INGARCHX model, and the other is the threshold INGARCHX model where ‘X’ refers to two meteorological covariates: weekly cumulative rainfall and weekly mean maximum temperature. Suppose that $s_{t}$ follows a first-order Markov chain with the following transition matrix:

℘ = [\begin{matrix} \Pr (s_{t} = 1 ∣ s_{t - 1} = 1) & \Pr (s_{t} = 2 ∣ s_{t - 1} = 1) \\ \Pr (s_{t} = 1 ∣ s_{t - 1} = 2) & \Pr (s_{t} = 2 ∣ s_{t - 1} = 2) \end{matrix}] = [\begin{matrix} p_{11} & p_{12} \\ p_{21} & p_{22} \end{matrix}],

(2.1)

where $p_{ij}$ $(i, j = 1, 2)$ denotes the transition probabilities of $s_{t} = j$ given that $s_{t - 1} = i$ . Clearly, the transition probabilities satisfy $p_{i 1} + p_{i 2} = 1$ . A challenging task with the MS INGARCH model is that it may not easily obtain estimations, because the state variables are unobservable. We overcome this difficulty using the Bayesian approach.

2.1 Markov switching INGARCHX model

Let ${Y_{t}}$ represent a univariate time series of DHF cases at time $t$ . Denote two observed meteorological covariates $X_{t} = {X_{1 t}, X_{2 t}}$ , where $X_{1 t}$ and $X_{2 t}$ represent the weekly cumulative rainfall and the weekly average maximum temperate at time $t$ , respectively. We here consider modelling the time-series weekly DHF cases using lag-one meteorological covariates (see Martínez-Bello et al., 2017). Assume that the conditional distribution follows a negative binomial distributed $NB (r, p_{t})$ with mean $r λ_{t}$ (for a real positive $r$ ). In practice, we present $λ_{t} = (1 - p_{t}) / p_{t}$ . We denote the MS INGARCHX model with state variables in (2.1) by

p (Y_{t} | ℱ_{t - 1}, X_{t - 1}) = (\frac{Y_{t} + r - 1}{r - 1}) {(\frac{1}{1 + λ_{t}})}^{r} {(\frac{λ_{t}}{1 + λ_{t}})}^{Y_{t}}, Y_{t} \geq 0

(2.2)

λ_{t} = \{\begin{matrix} α_{0}^{(1)} + α_{1}^{(1)} Y_{t - 1} + β_{1}^{(1)} λ_{t - 1} + \sum_{i = 1}^{2} γ_{i}^{(1)} X_{it - 1}, & if : s_{t} = 1 \\ α_{0}^{(2)} + α_{1}^{(2)} Y_{t - 1} + β_{1}^{(2)} λ_{t - 1} + \sum_{i = 1}^{2} γ_{i}^{(2)} X_{it - 1}, & if : s_{t} = 2, \end{matrix}

(2.3)

where $F_{t} = σ (Y_{t - 1}, Y_{t - 2}, \dots)$ denotes all information available up to time $t - 1$ . The Markovian state variable in (2.1), $s_{t}$ , switches the mechanism of the model structure at time $t$ , and its transition probabilities indicate the persistence of each state, because the state variables are unobservable. Here, $α_{i} = (α_{0}^{(i)}, α_{1}^{(i)}, β_{1}^{(i)}, γ_{1}^{(i)}, γ_{2}^{(i)})^{'}$ is a non-negative parameter vector in state $i$ .

2.2 Threshold INGARCHX model

We next conduct the threshold INGARCH model by incorporating covariates. Assume that $Y_{t}$ and $X_{t}$ are observed to be the same as in the MS INGARCHX model aforementioned. Therefore, we define the threshold INGARCHX as follows:

p (Y_{t} | F_{t - 1}, X_{t - 1}) = (\begin{matrix} Y_{t} + r - 1 \\ r - 1 \end{matrix}) {(\frac{1}{1 + λ_{t}})}^{r} {(\frac{λ_{t}}{1 + λ_{t}})}^{Y_{t}}, Y_{t} \geq 0

(2.4)

λ_{t} = \{\begin{matrix} α_{0}^{(1)} + α_{1}^{(1)} Y_{t - 1} + β_{1}^{(1)} λ_{t - 1} + \sum_{i = 1}^{2} γ_{i}^{(1)} X_{it - 1}, & if : Y_{t - d} \leq c \\ α_{0}^{(2)} + α_{1}^{(2)} Y_{t - 1} + β_{1}^{(2)} λ_{t - 1} + \sum_{i = 1}^{2} γ_{i}^{(2)} X_{it - 1}, & if : Y_{t - d} > c, \end{matrix}

where $Y_{t - d}$ is the threshold variable determining the dynamic switching mechanism of the model, $d$ is a delay lag and $c$ is the threshold value. Based on Zhu (2011), the conditional mean and variance of $Y_{t}$ in (2.3) and (2.4) are

\begin{matrix} E (Y_{t} | F_{t - 1}, X_{t - 1}) = r λ_{t}, : var (Y_{t} | F_{t - 1}, X_{t - 1}) = r λ_{t} (1 + λ_{t}) . \end{matrix}

This indicates that $var (Y_{t} | F_{t - 1}, X_{t - 1}) > E (Y_{t} | F_{t - 1}, X_{t - 1})$ . When $γ_{1}^{(i)} = γ_{2}^{(i)} = 0$ , the second-order stationary condition for a linear negative binomial INGARCH model is

r {(α_{1}^{(i)})}^{2} + {(r α_{1}^{(i)} + β_{1}^{(i)})}^{2} < 1, i = 1 .

(See Remark 3 of Zhu, 2011; Chen and Lee, 2017). Both models in Equations (2.3) and (2.4) suppose that the initial observation $Y_{1}$ is given. We consider the initial value $λ_{1}$ of the intensity process and $r$ (the number of successful trials) to be unknown parameters. We conduct Bayesian estimation, model comparison and forecasting for the proposed models.

2.3 Bayesian inference

We next discuss how to make inferences of the model parameters and join together unobserved state components for the Markov switching model with a Bayesian framework.

Assume that $Y = (Y_{1}, \dots, Y_{n})^{'}$ and $X = (X_{1}, X_{2})$ , and we define the conditional likelihood function of the MS INGARCHX model as

p (Y | X, s, θ) = \prod_{t = 2}^{n} \frac{Γ (Y_{t} + r)}{Γ (Y_{t} + 1) Γ (r)} {(\frac{1}{1 + λ_{t}})}^{r} {(\frac{λ_{t}}{1 + λ_{t}})}^{Y_{t}},

(2.5)

where $θ = (α_{1}^{'}, α_{2}^{'}, p_{11}, p_{22}, r, λ_{1})^{'}$ represents the vector of all unknown parameters, $s = (s_{2}, \dots, s_{n})^{'}$ , and $λ_{1}$ denotes the initial intensity. In addition to the parameter vector $θ$ , we also treat unobserved state variables ${s_{t}, t = 2, \dots, n}$ as parameters.

In practice, we believe that the coefficients of two climatic factors should have a positive relationship with weekly DHF cases, that is, $γ_{i} \geq 0, : i = 1, 2$ . To assign the required constraints on $α_{i}$ for each regime $i$ , we use a constrained uniform prior on the parameters in $α_{i}, : i = 1, 2$ imposed by an indicator $I (A)$ , that is, $p (α_{i}) \propto I (A)$ , where A is the set that satisfies

α_{0}^{(i)} > 0, α_{1}^{(i)}, β_{1}^{(i)}, γ_{j}^{(i)} \geq 0, r (α_{1}^{(i)})^{2} + (r α_{1}^{(i)} + β_{1}^{(i)})^{2} < 1, i, j = 1, 2 .

(2.6)

We consider $λ_{1}$ and $r$ to follow the gamma distribution as priors. For the intensity $λ_{1}$ , the prior is $Gamma (a_{1}, a_{2}),$ where $a_{1}$ is the shape parameter and $a_{2}$ is a rate parameter. For the number of successful trials $r$ , the prior follows a truncated gamma distribution to be greater than 1 with hyper-parameters $(k_{1}, k_{2})$ , where $k_{1}$ and $k_{2}$ are the shape and rate parameters, respectively. For the Markov switching model, we further utilize a conjugate prior distribution, $Beta (u_{i 1}, u_{i 2})$ , for $p_{ii}$ by

p (p_{ii}) \propto p_{ii}^{u_{i 1} - 1} (1 - p_{ii})^{u_{i 2} - 1}, i = 1, 2 .

(2.7)

To solve a label switching problem for Markov switching models, we use the condition of the first-order stationary for an INGARCH(1,1) as a restriction to relabel the MCMC outputs to have a higher intensity in state 2 as

\frac{r α_{0}^{(1)}}{1 - r α_{1}^{(1)} - β_{1}^{(1)}} \leq \frac{r α_{0}^{(2)}}{1 - r α_{1}^{(2)} - β_{1}^{(2)}},

(2.8)

where $r α_{0} / (1 - r α_{1} - β_{1})$ is the first-order stationary for the INGARCH(1,1) model (see Remark 2 of Zhu, 2011).

We subsequently use parameter groups (i) $α_{i}, : i = 1, 2$ ; (ii) ${p_{11}, p_{22}}$ ; (iii) $r$ ; (iv) $λ_{1}$ ; and (v) $s_{t}$ for inference. Suppose that the parameter groups are a priori independent. The posterior conditional distribution for each group is proportional to the likelihood function in (2.6) times the prior density for a stated group as

p (θ_{ℓ} | Y, X, s, θ_{\neq ℓ}) \propto p (Y | X, s, θ) p (θ_{ℓ} | θ_{\neq ℓ}),

(2.9)

where $θ_{ℓ}$ denotes the parameter(s) of each group, $p (θ_{ℓ})$ is its prior distribution and $θ_{\neq ℓ}$ is the vector of all model parameters, excluding $θ_{ℓ}$ . We describe the conditional posterior distributions as follows.

The posterior conditional distribution of $α_{i}, : i = 1, : 2$ , is

p (α_{i} | Y, X, s, α_{j \neq i}, p_{11}, p_{22}, r, λ_{1}) \propto p (Y | X, s, θ) p (α_{i}) .

(2.10)

The posterior conditional distribution of the transition probability $p_{ii}$ follows a Beta distribution, which is

p_{ii} | Y, X, s, θ_{\neq p_{ii}} \sim Beta (u_{i 1} + n_{ii}, u_{i 2} + n_{i} - n_{ii}), : i = 1, 2,

where $n_{ii}$ denotes the number of observations in state $i$ at time $t - 1$ that remain at state $i$ at time $t$ , and $n_{i}$ denotes the number of observations in state $s = i$ .

The posterior conditional distribution for $r$ is

p (r | Y, X, s, α_{1}, α_{2}, p_{11}, p_{22}, λ_{1}) \propto p (Y | X, s, θ) p (r) .

(2.11)

The posterior conditional distribution of $λ_{1}$ is

p (λ_{1} | Y, X, s, α_{1}, α_{2}, p_{11}, p_{22}, r) \propto p (Y | X, s, θ) p (λ_{1}) .

(2.12)

The posterior conditional distribution of state $s_{t}, : t = 2, \dots, n - 1$ , is

\begin{matrix} p (s_{t} | Y, X, s_{- t}, θ) & \propto & \prod_{t = 2}^{n} p (Y_{t} | Y, X, s, θ) \\ \times : p (s_{t} | s_{t - 1}, p_{11}, p_{22}) p (s_{t + 1} | s_{t}, p_{11}, p_{22}), \end{matrix}

(2.13)

where $s_{- t}$ denotes the state vector $s$ excluding the $t$ th element. Furthermore, the posterior conditional distribution of state $s_{n}$ , when $t = n$ , is the same as (2.14), but the last term $p (s_{t + 1} | s_{t}, p_{11}, p_{22})$ disappears. Consequently, we indicate the posterior conditional probability of $s_{t} = i$ by

\Pr (s_{t} = i | Y, X, s_{- t}, θ) = \frac{p (Y | X, s_{t} = i, s_{- t}, θ) \Pr (s_{t} = i | s_{- t}, θ)}{\sum_{j = 1}^{2} p (Y | X, s_{t} = j, s_{- t}, θ) \Pr (s_{t} = j | s_{- t}, θ)} .

(2.14)

Herein, if a random number generated is less than $\Pr (s_{t} = 1 | Y, X, s_{- t}, θ)$ , then we set $s_{t} = 1$ and 2 otherwise.

All conditional posterior distributions are non-standard forms except for ${p_{11}, p_{22}}$ . We employ the random walk Metropolis method for the parameters $r$ and $λ_{1}$ . To accelerate convergence, we use an adaptive Metropolis—Hastings (MH) algorithm for $α_{i}$ that incorporates a random walk Metropolis algorithm and an independent-kernel MH algorithm (e.g., Chen and So, 2006). We further monitor the MCMC convergence of each parameter using ACF and trace plots. An adaptive MCMC method can help to guarantee adequate coverage of an appropriate acceptance rate of the posterior distribution to be between 25% and 50% (see Gelman et al., 1996).

To account for the threshold INGARCHX model, we denote the vector of all unknown parameters in Equation (2.4) by $ϑ = (α_{1}', α_{2}', c, d, r, λ_{1})'$ , where a parameter vector is in each regime $j$ , $α_{j} = (α_{0}^{(j)}, α_{1}^{(j)}, β_{1}^{(j)}, γ_{1}^{(j)}, γ_{2}^{(j)})',$ $j = 1, 2$ . Assume that the initial observations $(Y_{1}, \dots, Y_{d_{0}})$ are given, where $d_{0}$ is a maximum delay. Given that $Y = (Y_{d_{0 + 1}}, \dots, Y_{n})$ , the conditional likelihood function of the threshold INGARCHX model is

p (Y | X, θ) = \prod_{t = d_{0} + 1}^{n} \{\sum_{j = 1}^{2} [\frac{Γ (Y_{t} + r)}{Γ (Y_{t} + 1) Γ (r)} {(\frac{1}{1 + λ_{t}})}^{r} {(\frac{λ_{t}}{1 + λ_{t}})}^{Y_{t}} I_{jt}]\},

(2.15)

where $I_{1 t}$ represents an indicator variable $I (Y_{t - d} < c)$ , and $I_{2 t} = 1 - I_{1 t}$ . We use a continuous yet constrained uniform prior on $c$ , that is, $c \sim Unif (l, u)$ and afterward discretized the estimation as to be an integer. The prior for the delay lag is a discrete uniform with maximum delay $d_{0}$ , that is, $p (d) = 1 / d_{0}$ for $d = 1, \dots, d_{0}$ . We assume that $α_{j}, : j = 1, 2$ , follows the uniform prior $p (α_{j}) \propto I (A)$ that satisfies Equation (2.6). Next, we adopt the same priors for $λ_{1}$ and $r$ as in Markov switching models. We divide the parameters of interest into the following parameter groups: (i) $α_{j}, : j = 1, 2$ , (ii) $c$ , (iii) $r$ , (iv) $λ_{1}$ and (v) $d$ for the threshold INGARCHX model.

2.4 Forecasting, model selection and checking model adequacy

It is very important to evaluate forecasting performance based on the fitted model. We can forecast ${Y_{t + i}, i \geq 1}$ by adding one extra step to each iteration of the MCMC sampler described in the previous sub-session. The proposed procedure is very similar to the idea of Brockwell (2007), as a by-product of the sampling algorithm, whereby the predictive distribution for $Y_{t + 1}$ takes into account model uncertainty and can easily be obtained. Herein, we use current observation $Y_{t}$ to predict one week ahead, and an alarm is triggered by an indicated triangle whenever the observation $Y_{t + 1}$ exceeds the threshold upper boundary of the 95% credible intervals based on a posterior predictive distribution. The predictive likelihood $p (Y_{t + 1} | F_{t}, X_{t})$ is the one-step-ahead predictive density for model $M_{1}$ evaluated at count $Y_{t + 1}$ :

\begin{matrix} p (Y_{t + 1} | F_{t}, X_{t}) & = & \int p (Y_{t + 1} | F_{t}, X_{t}, θ) p (θ | F_{t}, X_{t}) d θ \end{matrix}

(2.16)

\approx \frac{1}{N - M} \sum_{i = M + 1}^{N} p (Y_{t + 1} | F_{t}, X_{t}, θ^{[i]}),

(2.17)

where ${θ^{[M + 1 : N]}}$ are the MCMC draws from the posterior distributions for the model given the information ${F_{t}, X_{t}}$ .

For comparison of different models, we use a formal Bayesian approach. Let there be two models $M_{1}$ and $M_{2}$ with the parameter vectors $θ_{1}$ and $θ_{2}$ , respectively. Denote the marginal likelihood under model $M_{j}$ by $p ({Y | X, M}_{j})$ . The ratio of posterior probabilities of the models, $p (M_{1} | Y, X)$ and $p (M_{2} | Y, X)$ , is the posterior odds ratio (POR) given by

\frac{p (M_{1} | Y, X)}{p (M_{2} | Y, X)} = \frac{p (Y | X, M_{1})}{p (Y | X, M_{2})} \times \frac{\Pr (M_{1})}{\Pr (M_{2})} .

(2.18)

The RHS of Equation (2.18) expresses the Bayes factor (BF) and the prior odds ratio, respectively. In this case, assuming that we have equal preferences of these models, their prior probabilities are $\Pr (M_{1}) = \Pr (M_{2}) = 0.5$ . If the ratio in (2.18) is larger than 1, then model $M_{1}$ is favoured; if the ratio is less than 1, then model $M_{2}$ is favoured. In a Bayesian framework, a wildly used statistics for model selections is a deviance information criterion (DIC), as recommended by Spiegelhalter et al. (2002):

DIC = E [- 2 log p (Y | X, θ)] + p_{D},

where $p_{D} = E [- 2 log p (Y | X, θ)] + 2 log p (Y | X, \tilde{θ})$ , and $\tilde{θ}$ denotes the posterior mean of $θ$ . DIC directly uses the output from the MCMC samples. The model with the smallest DIC is favoured among all competing models.

A key part for the analysis of the time series is to check the adequacy of the specified model. Jung et al. (2006) offer the standardized Pearson residuals in order to check a dynamic structure of the mean and variance by

a_{t} = \frac{Y_{t} - E (Y_{t} | F_{t - 1}, X_{t - 1})}{\sqrt{var (Y_{t} | F_{t - 1}, X_{t - 1})}}, t = 1, \dots, n .

Clearly, the mean and variance of these residual should be close to zero and unit, respectively, without any significant serial correlation in $a_{t}$ and $a_{t}^{2}$ if a specified model is correct.

3 A simulation study

We now conduct a simulation study to illustrate the performance of the MCMC methods, in terms of parameter estimation and model selection. We simulate samples of size $n = 700$ from two specific models, specified as

\begin{array}{l} M o d e l 1 : Y_{t} | F_{t - 1}, X_{t - 1} ~ N B (5, \frac{1}{1 + λ_{t}}), \\ λ_{t} = {\begin{array}{l} 0.2 + 0.02 Y_{t - 1} + 0.20 λ_{t - 1} + 0.02 X_{1 t - 1}, & i f s_{t} = 1 \\ 0.5 + 0.02 Y_{t - 1} + 0.25 λ_{t - 1} + 0.01 X_{1 t - 1}, & i f s_{t} = 2, \end{array} \\ \Pr (s_{t} = 1 | s_{t - 1} = 1) = 0.92, \Pr (s_{t} = 2 | s_{t - 1} = 2) = 0.89, \end{array}

(3.1)

\begin{array}{l} Model 2 : Y_{t} | F_{t - 1}, X_{t - 1} ~ NB (5, \frac{1}{1 + λ_{t}}), \\ λ_{t} = {\begin{array}{l} 0.3 + 0.04 Y_{t - 1} + 0.20 λ_{t - 1} + 0.02 X_{1 t - 1}, & Y_{t - d} \leq 5 \\ 0.6 + 0.02 Y_{t - 1} + 0.15 λ_{t - 1} + 0.01 X_{1 t - 1}, & Y_{t - d} > 5, \end{array} \end{array}

(3.2)

where $X_{1 t}$ generates from a normal distribution, $N (20, 1)$ , which is then discretized as an integer. For Model 1, we randomly select the state process, $s$ , according to a first-order Markov process in (2.1) to generate values ${1, 2}$ as an initial guess based on $p_{ii} = 0.5, i = 1, 2$ . We set up the initial values of intensity parameters $α_{i} = (0.1, 0.1, 0.1, 0.1)^{'}$ . The priors of $p_{ii}, i = 1, 2$ follow a $Beta (90, 10)$ distribution for transition probabilities. Next, we select the hyper-parameters for the prior of $λ_{1}$ and $r$ to be $Gamma (5, 1)$ and $Gamma (25, 5)$ , respectively.

Model 2 is the same as the previous set-up except for the true value of $d = 1$ and $d_{0} = 3$ is the maximum delay. We employ maximum a posteriori probability (MAP) estimation for $(c, d)$ since both parameters are integers. Note that we choose other different hyper parameters, and as a result, the estimates are very similar. Therefore, we only present this current set-up. Finally, we take a total sample of iterations $N = 30 000$ using a burn-in sample $M = 8 00$ and store every fifth iteration within simulating 200 replications for each model.

Table 1 reports simulation results containing the averages of posterior means, medians and 95% credible intervals over the 200 replications for Models 1 and 2. The coverage probabilities indicate the percentages of the 95% posterior intervals that cover the true values generated by MCMC iterations. Overall, the posterior means of Models 1 and 2 are near the true values with coverage rates of at least 83%.

In order to make a fair comparison between two models, the conditional likelihood functions start with 4. For model selection using BFs to choose the best model, we provide simulation results in Table 2 from 200 replications from Models 1 and 2 in Table 1 where 700 datasets are tested. The performance of model selection is sound under the correct model selection rates being 88% and 86% for the two scenarios, respectively.

Table 1:

Simulation results for Models 1 and 2 obtained from 200 replications

Par.	Model 1					Model 2
Par.	True	Mean	2.5%	97.5%	Coverage	True	Mean	2.5%	97.5%	Coverage
$α_{0}^{(1)}$	0.20	0.217	0.045	0.418	0.900	0.30	0.306	0.034	0.593	1.000
$α_{1}^{(1)}$	0.02	0.014	0.001	0.039	0.980	0.04	0.036	0.002	0.100	1.000
$β_{1}^{(1)}$	0.20	0.222	0.032	0.465	0.950	0.20	0.224	0.025	0.481	0.980
$γ_{1}^{(1)}$	0.02	0.025	0.007	0.041	0.995	0.02	0.021	0.003	0.043	1.000
$α_{0}^{(2)}$	0.50	0.444	0.213	0.699	0.830	0.60	0.564	0.285	0.865	1.000
$α_{1}^{(2)}$	0.02	0.015	0.001	0.041	0.980	0.02	0.017	0.001	0.057	1.000
$β_{1}^{(2)}$	0.25	0.240	0.036	0.498	0.950	0.15	0.141	0.009	0.330	0.960
$γ_{1}^{(2)}$	0.01	0.013	0.001	0.029	0.985	0.01	0.008	0.001	0.024	0.900
$p_{11}$	0.92	0.905	0.835	0.957	1.000
$p_{22}$	0.89	0.885	0.812	0.942	1.000
$r$	5	4.943	3.250	6.882	0.995	5	4.919	3.169	7.005	1.000
$c^{★}$						5	5.045	4.695	5.890	1.000
$d^{★}$						1	1	1	1	0.995

Note: signifies MAP estimate for each iterate.

Table 2:

Comparisons of Models 1 and 2 based on 200 replications

Criteria		Bayes factors
	Fitted model
True model		Model 1	Model 2	No. of simulations
Model 1		176	24	200
Model 2		28	172	200

Note: The number in the table indicates ln(BF) > 0.

4 Analytical results

This study designs retrospective data analysis of time series for the weekly DHF cases in the four provinces of Maha Sarakham, Kalasin, Sakon Nakhon and Mukdahan in northeastern Thailand, covering approximately 12 years from January 2006 to October 2017 (612 weeks) as collected by the Bureau of Epidemiology, Department of Disease Control (DDC), Ministry of Public Health (MOPH), Thailand. For meteorological covariates, we take weekly cumulative rainfall (mm) and weekly average maximum temperature (C) from Thai Meteorological Department (TMD), Thailand. Figure 1 specifies the geographic locations in order of the four provinces.

We consider the time series of the weekly DHF cases (total sample size 612) for the four provinces, denoted by DS1, DS2, DS3 and DS4, in northeastern Thailand as shown in Table 3. All datasets of weekly DHF cases appear to be overdispersed. Table 3 reports the summary statistics for all DS1–DS4, including the percentage of zeros and consecutive zeros (at least two consecutive zeros), and also presents weekly cumulative rainfall and weekly mean maximum temperature. We observe that DS1–DS4 have at least 7% and 5% zeros and consecutive zeros, respectively. Figure 2 shows time-series plots for DS1–DS4. Obviously, these datasets include consecutive zeros, especially Sakon Nakhon (DS3) and Mukdahan (DS4).

Figure 1:

Geographic map of four provinces in northeastern Thailand

Table 3:

Summary statistics for the weekly DHF cases, cumulative rainfall (mm) and mean maximum temperature (C)

	Province	Mean	Variance	Min	Q1	Median	Q3	Max	Percentage of 0
DS1	Maha Sarakham	15.77	486.66	0	3	8	21	157	8.50	(5.39)
DS2	Kalasin	14.08	300.69	0	3	9	17	149	7.35	(5.23)
DS3	Sakon Nakhon	9.25	232.02	0	1	3	11	97	23.69	(18.79)
DS4	Mukdahan	7.41	130.10	0	1	4	9	92	17.81	(12.91)

	Zone No.	Rainfall					Temperature
		Mean	Median	Min	Max		Mean	Median	Min	Max
DS1		27.63	6.40	0	317.7		33.47	33.47	22.17	41.33
DS2	27.13	7.20	0	303.0		32.28	32.19	22.73	39.96
DS3	8	34.15	12.05	0	456.1		32.00	32.07	22.76	40.61
DS4	10	28.14	6.95	0	293.9		32.76	32.94	23.68	40.39

Note: The number in brackets of the last column means at least 2 consecutive zeros (%).

Figure 2:

Weekly DHF cases in northeastern Thailand, DS1: Maha Sarakham; DS2: Kalasin; DS3: Sakon Nakhon; DS4: Mukdahan

Figure 3:

ACF plots of weekly DHF cases, DS1–DS4

Figure 3 depicts the ACF plots for the DHF cases in the four provinces. It shows that there exists series dependence in each DHF cases. Here ${X_{1 t}, X_{2 t}}$ are the exogenous variables, which are weekly cumulative rainfall and weekly mean maximum temperature of the target province, respectively. In this scheme, the information about ${X_{1 t}, X_{2 t}}$ helps predict $Y_{t + 1}$ , as do the previous values of $Y_{t}$ themselves.

Figure 4 plots the meteorological covariates’ time series of cumulative rainfall and temperature for each dataset, respectively, which exhibit a repeated pattern. For computational convenience, we transform $X_{it}$ to be $x_{it} = [X_{it} - min (X_{i})] / S_{x}, : i = 1, 2$ , where $S_{x}$ stands for the standard deviation of $X_{i}$ , since magnitudes of some $X_{it}$ are truly large. We fit linear, Markov switching and threshold INGARCHX models, respectively denoted by Models 1, 2 and 3, for each province and provide Bayesian model comparisons. Table 4 presents the natural logarithm of the BFs and DIC values of these three candidate models. The natural logarithm values of the BFs are all positive in Model 2 versus Model 3 for the four provinces, while they are all negative for Model 1 versus Model 2 (and Model 3). This indicates that the BFs INGARCHX model is superior among three models for describing the DHF cases. The smallest value of DIC indicates the best-fitting model among the three models for each province. All datasets are in favour of the MS INGARCHX model as well as are based on DIC criterion. In summary, the Markov switching negative binomial INGARCHX model is a favourable model for all datasets.

Figure 4:

Time-series plots of meteorological covariates. Weekly cumulative rainfall (upper panel) and weekly mean maximum temperature (lower panel) for DS1–DS4, respectively

Table 4:

Model comparison between the three models for DS1-DS4

Criterion	Model	DS1	DS2	DS3	DS4
	$M_{1}$ vs. $M_{2}$	$-$ 37.27	$-$ 78.02	$-$ 77.84	$-$ 68.87
$log (BF)$	$M_{1}$ vs. $M_{3}$	$-$ 5.48	$-$ 5.94	$-$ 16.56	$-$ 12.14
	$M_{2}$ vs. $M_{3}$	31.79	72.08	61.28	56.71
	$M_{1}$	3614.86	3640.69	2942.26	2900.69
DIC	$M_{2}$	3601.81	3552.53	2856.82	2787.75
	$M_{3}$	3608.79	3630.50	2911.18	2881.70

Note: $M_{1}$ is negative binomial INGARCHX model

$M_{2}$ is Markov switching negative binomial INGARCHX model; and

$M_{3}$ is Threshold negative binomial INGARCHX model.

Table 5:

Bayesian estimates of the two INGARCHX models for DS1

Par.	DS1: MS INGARCHX				DS1: Threshold INGARCHX
Par.	Mean	Median	2.5%	97.5%	Mean	Median	2.5%	97.5%
$α_{0}^{(1)}$	0.018	0.015	0.002	0.065	0.040	0.040	0.017	0.065
$α_{1}^{(1)}$	0.078	0.078	0.062	0.090	0.074	0.074	0.062	0.086
$β_{1}^{(1)}$	0.289	0.286	0.207	0.384	0.341	0.341	0.294	0.386
$γ_{1}^{(1)}$	0.028	0.024	0.002	0.080	0.071	0.071	0.029	0.116
$γ_{2}^{(1)}$	0.011	0.011	0.002	0.024	0.009	0.008	0.001	0.017
$α_{0}^{(2)}$	0.101	0.010	0.016	0.199	0.095	0.095	0.051	0.138
$α_{1}^{(2)}$	0.076	0.076	0.063	0.089	0.097	0.097	0.088	0.105
$β_{1}^{(2)}$	0.344	0.346	0.237	0.424	0.292	0.292	0.248	0.334
$γ_{1}^{(2)}$	0.057	0.057	0.006	0.117	0.040	0.040	0.009	0.074
$γ_{2}^{(2)}$	0.018	0.017	0.002	0.042	0.056	0.057	0.015	0.097
$p_{11}$	0.816	0.819	0.737	0.881
$p_{22}$	0.804	0.807	0.723	0.874
$r$	8.033	7.978	7.225	9.116	6.667	6.655	6.149	7.251
$c^{★}$					5
$d^{★}$					1

Note: signifies MAP estimate.

Table 6:

Bayesian estimates of the two INGARCHX models for DS2

Par.	DS2: MS INGARCHX				DS2: Threshold INGARCHX
Par.	Mean	Median	2.5%	97.5%	Mean	Median	2.5%	97.5%
$α_{0}^{(1)}$	0.019	0.019	0.010	0.029	0.115	0.114	0.088	0.141
$α_{1}^{(1)}$	0.059	0.058	0.049	0.069	0.071	0.071	0.053	0.090
$β_{1}^{(1)}$	0.415	0.416	0.361	0.465	0.393	0.393	0.348	0.438
$γ_{1}^{(1)}$	0.020	0.020	0.010	0.030	0.020	0.019	0.002	0.039
$γ_{2}^{(1)}$	0.005	0.005	0.001	0.008	0.007	0.007	0.001	0.015
$α_{0}^{(2)}$	0.334	0.334	0.263	0.407	0.071	0.071	0.017	0.130
$α_{1}^{(2)}$	0.095	0.095	0.085	0.106	0.118	0.118	0.106	0.130
$β_{1}^{(2)}$	0.118	0.118	0.045	0.193	0.215	0.215	0.164	0.267
$γ_{1}^{(2)}$	0.100	0.099	0.050	0.149	0.117	0.116	0.049	0.186
$γ_{2}^{(2)}$	0.053	0.052	0.012	0.095	0.055	0.054	0.012	0.104
$p_{11}$	0.850	0.852	0.793	0.899
$p_{22}$	0.781	0.783	0.700	0.854
$r$	7.948	7.932	7.263	8.774	5.928	5.918	5.471	6.447
$c^{★}$					4
$d^{★}$					1

Note: signifies MAP estimate.

Table 7:

Bayesian estimates of the two INGARCHX models for DS3

Par.	DS3: MS INGARCHX				DS3: Threshold INGARCHX
Par.	Mean	Median	2.5%	97.5%	Mean	Median	2.5%	97.5%
$α_{0}^{(1)}$	0.020	0.020	0.013	0.027	0.030	0.030	0.015	0.045
$α_{1}^{(1)}$	0.074	0.074	0.068	0.081	0.109	0.109	0.087	0.128
$β_{1}^{(1)}$	0.385	0.385	0.363	0.408	0.454	0.452	0.395	0.513
$γ_{1}^{(1)}$	0.008	0.008	0.003	0.013	0.014	0.014	0.002	0.028
$γ_{2}^{(1)}$	0.004	0.004	0.001	0.008	0.007	0.007	0.001	0.014
$α_{0}^{(2)}$	0.363	0.362	0.301	0.424	0.080	0.080	0.019	0.145
$α_{1}^{(2)}$	0.127	0.127	0.118	0.135	0.172	0.172	0.153	0.193
$β_{1}^{(2)}$	0.115	0.115	0.053	0.175	0.154	0.154	0.100	0.210
$γ_{1}^{(2)}$	0.129	0.129	0.089	0.170	0.053	0.052	0.011	0.095
$γ_{2}^{(2)}$	0.245	0.245	0.187	0.305	0.177	0.177	0.087	0.265
$p_{11}$	0.906	0.907	0.871	0.935
$p_{22}$	0.803	0.805	0.723	0.867
$r$	6.018	6.013	5.474	6.612	4.371	4.366	3.950	4.805
$c^{★}$					4
$d^{★}$					1

Note: signifies MAP estimate.

Table 8:

Bayesian estimates of the two INGARCHX models for DS4

Par.	DS4: MS INGARCHX				DS4: Threshold INGARCHX
Par.	Mean	Median	2.5%	97.5%	Mean	Median	2.5%	97.5%
$α_{0}^{(1)}$	0.020	0.020	0.011	0.029	0.026	0.026	0.014	0.039
$α_{1}^{(1)}$	0.080	0.080	0.071	0.089	0.083	0.083	0.074	0.093
$β_{1}^{(1)}$	0.404	0.404	0.371	0.437	0.475	0.475	0.440	0.508
$γ_{1}^{(1)}$	0.008	0.008	0.003	0.013	0.008	0.008	0.001	0.017
$γ_{2}^{(1)}$	0.007	0.007	0.003	0.012	0.011	0.011	0.003	0.019
$α_{0}^{(2)}$	0.450	0.450	0.380	0.518	0.120	0.120	0.042	0.199
$α_{1}^{(2)}$	0.114	0.114	0.100	0.128	0.129	0.129	0.115	0.144
$β_{1}^{(2)}$	0.119	0.119	0.060	0.177	0.110	0.110	0.030	0.193
$γ_{1}^{(2)}$	0.137	0.137	0.087	0.187	0.045	0.043	0.006	0.094
$γ_{2}^{(2)}$	0.170	0.171	0.112	0.225	0.243	0.241	0.108	0.381
$p_{11}$	0.910	0.911	0.872	0.941
$p_{22}$	0.761	0.762	0.676	0.837
$r$	6.082	6.072	5.485	6.737	5.703	5.694	5.202	6.253
$c^{★}$					8
$d^{★}$					1

Note: signifies MAP estimate.

Figure 5:

Estimated $\Pr (s_{t} = 2 | Y, X, s_{- t}, θ)$ is based on the Markov switching INGARCHX model for DS1, DS2 (upper panel), DS3 and DS4 (lower panel), respectively

Tables 5 –8 report Bayesian estimates of the Markov switching and threshold INGARCHX models for DS1–DS4, including posterior means, medians and 95% credible intervals (25 and 97.5 percentiles), respectively. We provide some evidence of convergence of the Markov chains in the Appendix, which includes some trace plots and ACF plots of parameter iterates by fitting the Markov switching negative binomial INGARCHX model. The results are convincing becuase convergence is clear in each dataset. We illustrate via ACF plots that the adaptive method for $α_{i}$ indeed leads to reduced autocorrelations for MCMC iterates.

Moreover, ${γ_{1}^{(i)}, γ_{2}^{(i)}, i = 1, 2}$ appears in different regimes, which shows $γ_{1}^{(i)} < γ_{2}^{(i)}$ for both proposed models. The estimated transition probabilities $p_{22}$ and $p_{11}$ are at least 0.76, with $p_{22}$ less than $p_{11}$ , indicating that the transition probabilities reveal the persistence of state 1.

To assess model adequacy, we employ standardized residuals through the posterior estimates of the target model. By using the Ljung—Box $Q$ statistics based on the standardized residuals from lags 1 to 12, the result confirms no significant autocorrelation patterns and henceforth the adequacy of the fitted models.

Figure 5 provides the average probabilities of state 2, $\Pr (s_{t} = 2 | Y, X, s_{- t}, θ)$ , for DS1–DS4 based on the MS INGARCHX model. These figures deliver insight as to the state change captured in the modelled dataset. When the estimates of posterior probabilities are greater than 0.5, they can capture the high peaks in state 2. We observe that the estimated $\Pr (s_{t} = 2 | Y, X, s_{- t}, θ)$ are analogous to the pattern of weekly DHF time-series plots in Figure 2. For outbreak detection of DS1–DS4, we utilize the MS INGARCHX model to predict one-step-ahead using current information. Figures 6 and 7 display the observed and predicted DHF cases, including the blue line that indicates the upper bound of the 95% Bayesian credible intervals. The vertical line denotes $Y_{t + 1}$ and triangles (red colour) indicate alarms when the observed cases are higher than the corresponding upper bound of the 95% credible interval at time $t$ . We suggest that practitioners can increase the percentage of credible intervals if they want to reduce false alarms.

Figure 6:

DS1–DS2: One–step ahead forecast from the Markov switching INGARCHX models

Figure 7:

DS3–DS4: One–step ahead forecast from the Markov switching INGARCHX models

5 Conclusion

This study provides two nonlinear integer-valued GARCHX (Markov switching and threshold) models with negative binomial distribution for modelling weekly DHF cases that consist of meteorological covariates. Applying Bayesian inferences of these models demonstrates some fascinating statistical features, consisting of overdispersion, lagged dependence, consecutive zeros, nonlinear dynamics and switching coefficients of meteorological covariates.

This article demonstrates the benefits of forecasting using the posterior predictive distribution via the joint modelling of unknown parameters and state variables. We employ a Bayesian framework through the MCMC method for making inferences. Parameter estimation methods of simulation studies illustrate that they are appropriately estimated with sufficient coverage properties. Comparisons between the three competing models based on the BF and DIC confirm that the MS INGARCHX model outperforms the threshold INGARCHX model. The evidence herein offers strong support for the proposed MS INGARCHX model over its competitors. Nevertheless, some fundamental properties exists for the threshold INGARCHX model, which might properly describe other datasets. We use predictive credible intervals that take into account ‘upper threshold’ for monitoring and providing early warning signals of outbreaks. The posterior probabilities convey clearly the captured state changing in the modelled datasets.

Finally, the proposed models allow for a warning several weeks in advance of dengue epidemics, which involve forecasting weekly cumulative rainfall and average maximum temperature as well.

Appendix

We provide some evidence of convergence of the Markov chains in Section 4. Some trace plots and ACF plots of parameter iterates by fitting the Markov switching negative binomial INGARCHX model.

Figure A1:

DS2: Trace and ACF plots of MCMC estimates

Figure A2:

DS3: Trace and ACF plots of MCMC estimates

Figure A3:

DS4: Trace and ACF plots of MCMC estimates

Footnotes

Acknowledgments

We thank the editor, the associate editor and the anonymous referees for their valuable time and careful comments on our article, which have led to an improved version of it.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

Cathy WS Chen’s research is funded by the Ministry of Science and Technology, Taiwan (MOST107-2118-M-035-005-MY2).

References

Brockwell

(2007) Likelihood-based analysis of a class of generalized long-memory time series models. Journal of Time Series Analysis 28 386–407.

Chan

Hansen

Timmermann

(2017) Guest editors’ introduction: Regime switching and threshold models. Journal of Business & Economic Statistics 35 159–61.

Chareonsook

Foy

Teeraratkul

Silarug

(1999) Changing epidemiology of dengue hemorrhagic fever in Thailand. Epidemiology and Infection 122 161–66.

Chen

CWS

Khamthong

Lee

(2019) Markov switching integer-valued GARCH models for dengue counts. Journal of the Royal of Statistical Society Series C: Applied Statistics . doi: 10.1111/rssc.12344.

Chen

CWS

Lee

(2016) Generalized Poisson autoregressive models for time series of counts. Computational Statistics and Data Analysis 99 51–67.

Chen

CWS

Lee

(2017) Bayesian causality test for integer-valued time series models with applications to climate and crime data. Journal of the Royal Statistical Society Series C , 66 797–814.

Chen

CWS

MKP

(2006) On a threshold heteroscedastic model. International Journal of Forecasting 22 73–89.

Chen

CWS

MKP

Sriboonchitta

(2016) Autoregressive conditional negative binomial model applied to over-dispersed time series of counts. Statistical Methodology 31 73–90.

Chen

CWS

MKP

Lin

EMH

(2009) Volatility forecasting with double Markov switching GARCH models. Journal of Forecasting 28 681–97.

10.

Chen

CWS

MKP

Liu

(2011) A review of threshold time series models in finance. Statistics and Its Interface 4 167–82.

11.

De Melo

DPO

Scherrer

Eiras

ÁE

(2012) Dengue fever occurrence and vector detection by larval survey ovitrap and MosquiTRAP: A space-time clusters analysis. Public Library of Science One 7 e42125.

12.

Ferland

Latour

Oraichi

(2006) Integer-valued GARCH processes. Journal of Time Series Analysis 27 923–42.

13.

Fokianos

(2012) Count time series models. In Handbook of Statistics 30—Time Series Analysis: Methods and Applications edited by Rao

T Subba

Rao

S Subba

Rao

. Pages 315–347. Amsterdam: Elsevier B.V.

14.

Gelman

Roberts

Gilks

(1996) Efficient metropolis jumping rules. In Bayesian Statistics 5 edited by Bernardo

Berger

Dawid

Smith

AFM

. Pages 599–607. Oxford: Oxford University Press.

15.

Halide

Ridd

(2008) A predictive model for dengue hemorrhagic fever epidemics. International Journal of Environmental Health Research 18 253–65.

16.

Hamilton

(1989) A new approach to the economic analysis of nonstationary time series and the business cycle. Econometrica 57 357–84.

17.

Hii

Zhu

Rocklöv

(2012) Forecast of dengue incidence using temperature and rainfall. PLOS Neglected Tropical Diseases 6 e1908.

18.

Jung

Kukuk

Liesenfeld

(2006) Time series of count data: Modelling and estimation and diagnostics. Computational Statistics and Data Analysis 51 2350–64.

19.

Limkittikul

Brett

L’Azou

(2014) Epidemiological trends of dengue disease in Thailand (2000–2011): A systematic literature review. Public Library of Science Neglected Tropical Diseases 8 e3241.

20.

Lowe

Bailey

Stephenson

Graham

Coelho

CAS

Carvalho

Barcellos

(2011) Spatio-temporal modelling of climate-sensitive disease risk: Towards an early warning system for dengue in Brazil. Computers and Geosciences 37 371–81.

21.

Lowe

Bailey

Stephenson

Jupp

Graham

Barcellos

Carvalho

(2013) The development of an early warning system for climate-sensitive disease risk with a focus on dengue epidemics in Southeast Brazil. Statistic in Medicine 32 864–83.

22.

Malyshikina

Mannering

Tarko

(2009) Markov switching negative binomial models: An application to vehicle accident frequencies. Accident Analysis and Prevention 41 217–26.

23.

Manitz

Höhle

(2013) Bayesian outbreak detection algorithm for monitoring reported cases of campylobacteriosis in Germany. Biometrical Journal 55 509–26.

24.

Martínez- Bello

López-Quílez

Torres-Prieto

(2017) Bayesian dynamic modeling of time series of dengue disease case counts. PLOS Neglected Tropical Diseases 11 e0005696.

25.

Silawan

Singhasivanon

Kaewkungwal

Nimmanitya

Suwonkerd

(2008) Temporal patterns and forecast of dengue infection in northeastern Thailand. Southeast Asian Journal of Tropical Medicine and Public Health 39 90–8.

26.

Spiegelhalter

Best

Carlin

Van der Linde

(2002) Bayesian measures of model complexity and fit (with discussion). Journal of the Royal Statistical Society: Series B , 64 583–616.

27.

Tong

(1978) On a threshold model. In Pattern Recognition and Signal Processing, edited by Chen

. Amsterdam: Sijthoff & Noordhoff pp. 575–86.

28.

Tong

(1983). Threshold Models in Nonlinear Time Series Analysis. Lecture Notes in Statistics . New York, NY: Springer-Verlag.

29.

Truong

Chen

CWS

Sriboonchitta

(2017) Hysteretic Poisson INGARCH model for integer-valued time series Statistical Modelling 17 401–22.

30.

Wang

Liu

Yao

Davis

(2014) Self-excited threshold Poisson autoregression. Journal of the American Statistical Association 109 777–87.

31.

Wongkoon

Jaroensutasinee

(2012) Assessing the temporal modeling for prediction of dengue infection in northern and northeastern Thailand. Tropical Biomedicine 29 339–48.

32.

H-Y

Xie

Goh

(2012) A model for integer-valued time series with conditional overdispersion. Computational Statistics and Data Analysis 56 4229–42.

33.

Zhu

(2011) A negative binomial integer-valued GARCH model. Journal of Time Series Analysis 32 54–67.

Bayesian modelling of nonlinear negative binomial integer-valued GARCHX models

Abstract

Keywords

1 Introduction

2 Nonlinear negative binomial INGARCH models

Figure 1:

Geographic map of four provinces in northeastern Thailand

Weekly DHF cases in northeastern Thailand, DS1: Maha Sarakham; DS2: Kalasin; DS3: Sakon Nakhon; DS4: Mukdahan

ACF plots of weekly DHF cases, DS1–DS4

Time-series plots of meteorological covariates. Weekly cumulative rainfall (upper panel) and weekly mean maximum temperature (lower panel) for DS1–DS4, respectively

Estimated Pr ( s t = 2 | Y , X , s − t , θ ) is based on the Markov switching INGARCHX model for DS1, DS2 (upper panel), DS3 and DS4 (lower panel), respectively

DS1–DS2: One–step ahead forecast from the Markov switching INGARCHX models

DS3–DS4: One–step ahead forecast from the Markov switching INGARCHX models

Appendix

DS2: Trace and ACF plots of MCMC estimates

DS3: Trace and ACF plots of MCMC estimates

DS4: Trace and ACF plots of MCMC estimates

Footnotes

Acknowledgments

Declaration of conflicting interests

Funding

References

Estimated $\Pr (s_{t} = 2 | Y, X, s_{- t}, θ)$ is based on the Markov switching INGARCHX model for DS1, DS2 (upper panel), DS3 and DS4 (lower panel), respectively