A Bayesian Analysis for Negative Binomial INGARCH Models with Minimal Prior Information

A Bayesian Analysis for Negative Binomial INGARCH Models with Minimal Prior Information Xiaoyin Wang Yunwei Cui Department of Mathematics, Towson University Abstract This study presents a Bayesian framework for estimating Negative Binomial Integer-valued Generalized Autoregressive Conditional Heteroskedasticity (NB-INGARCH) models, designed to address over-dispersed count time series data. Minimally informative priors such as Dirichlet prior and hierarchical prior are employed to ensure the flexibility together with the model constraints. Numerical inference is carried out via Markov Chain Monte Carlo (MCMC) methods. Applications to empirical data demonstrate that the proposed Bayesian approach delivers more robust and interpretable estimates compared to frequentist methods. Keywords: Bayesian estimation, Markov Chain Monte Carlo, NB-INGARCH model, Dirichlet distribution, time series analysis, minimal-informative priors. 1 Introduction Count time series have broadly used in the various fields such as criminology, finance, public health, and marketing. Rececenly, it was applied to a study analyzing soft drink sales (Jia et al., 2023) and COVID-19 mortality data (Palmer et al., 2021). Early work in count time series modeling largely focused on conditional Poisson distributions (Davis et al., 2000; Ferland et al., 2006; Weiß, 2008; Fokianos et al., 2009, 2020; Kong and Lund, 2025). However, the Poisson assumption of equal mean and variance often fails in practice, particularly for over-dispersed data. Empirical studies have shown that non-Poisson alternatives, such as the Negative Binomial distribution, often capture this variability more effectively (see Davis and Wu, 2009; Zhu, 2011; Christou and Fokianos, 2014; Davis and Liu, 2016; Ahmad and Francq, 2016). Another notable shortcoming of existing research is that the estimation of the dispersion parameter has received limited attention, especially in classical settings. Zhu (2012) emphasized the need for further research and suggested moment-based approaches, while Davis and Liu (2016) introduced a two-stage profile likelihood method: first estimating the parameters 𝜶,𝜷, 𝛿 as functions of the dispersion parameter 𝜙, followed by optimizing the dispersion parameter using the profile likelihood. Nonetheless, these works did not provide a comprehensive inference for the dispersion parameter. Bayesian methods have been widely applied in continuous GARCH models (Bauwens and Lubrano, 1998; Aus´ ın and Galeano, 2007) and Poisson GARCH models (Ardia, 2008, 2009; Ardia and Hoogerheide, 2010). More recently, Chen et al. (2021) demonstrated that Bayesian analysis with naive priors shows Negative Binomial INGARCH models outperform Poisson-based alternatives. A key strength of the Bayesian framework is its ability to jointly estimate dispersion and model parameters within a unified process. The choice of prior distribution is critical in Bayes inference. The gamma prior has been commonly used for dispersion (Chen and Lee, 2017; Chen and Khamthong, 2020; Chen et al., 2021; Chu and Yu, 2023), yet the selection of the hyperparameters is not explicitly stated, and often relies on subjective judgment. 1 To address these limitations, this paper develops Bayesian estimation of the NB-INGARCH process under non-informative priors, allowing inference to rely primarily on observed data. For further reviews of count time series models, see Fokianos (2012); Davis et al. (2021); Liu et al. (2023). We begin by providing an overview of NB-INGARCH models and their dynamic specifications in Section 2. Following this, we introduce Bayesian estimation techniques, elaborating on the selection of prior distributions for model parameters and the derivation of the joint posterior distribution in the Section 3. Obtaining the closed-form expression for the posterior distribution of model parameters is challenging, and Markov Chain Monte Carlo (MCMC) methods can be applied to generate numerical results. To demonstrate the practical application of Bayesian estimation for NB-INGARCH models, Section 4 presents an empirical example, on the transaction data of Ericsson B stock. Since it is widely used in literature, this dataset serves as a standardized benchmark for comparing and evaluating our research findings. Finally, we present conclusions and remarks in Section 5. 2 The NB-INGARCH(p, q) Models This section presents the NB-INGARCH model, which provides a flexible framework for overdispersed count time series. The INGARCH process specifies counts 𝑌𝑡with a time-varying mean 𝜆𝑡evolving as 𝜆𝑡=𝛿+ 𝑝 ∑︁ 𝑖=1 𝛼𝑖𝜆𝑡−𝑖+ 𝑞 ∑︁ 𝑗=1 𝛽𝑗𝑌𝑡−𝑗,(1) where 𝛿 > 0 is the baseline level, 𝛼𝑖, 𝛽𝑗≥0 are coefficients for lagged means and counts, and 𝑝, 𝑞 denote the model orders. The equation (1) describes that the dynamics of 𝜆𝑡depends on the parameters 𝜶=𝛼1,· · · , 𝛼𝑝, 𝜷=𝛽1,· · · , 𝛽𝑞, and 𝛿, as well as the lagged value of (𝜆𝑡−1,· · · , 𝜆𝑡−𝑝)and (𝑌𝑡−1,· · · ,𝑌𝑡−𝑞). The state parameters 𝜶is often referred to as the persistence parameter, and 𝜷is the variability parameter, while 𝛿sets the long-term baseline. Conditioned on 𝜆𝑡,𝑌𝑡follows a Negative Binomial distribution with mean 𝜆and dispersion 𝜙: 𝑃(𝑋=𝑘)=𝑘+𝜙−1 𝑘 𝜙 𝜙+𝜆𝜙𝜆 𝜙+𝜆𝑘 , with variance 𝜆+𝜆2/𝜙. This mean-dispersion parameterization is especially useful for modeling the evolution of counts through the conditional mean. Davis and Liu (2016) show that in order to ensure the existence of a unique strictly stationary and ergodic solution for the time series {𝑌𝑡, 𝜆𝑡}with 𝐸(𝑌𝑡)<∞and 𝐸(𝜆𝑡)<∞, it requires Í𝑖𝛼𝑖+Í𝑗𝛽𝑗<1. The likelihood for observations 𝒀=(𝑌1, . . . , 𝑌𝑇)is 𝐿(𝜙, 𝜶,𝜷, 𝛿 |𝒀)= 𝑇 Ö 𝑡=1𝑦𝑡+𝜙−1 𝑦𝑡 𝜙 𝜙+𝜆𝑡𝜙𝜆𝑡 𝜙+𝜆𝑡𝑦𝑡 . 3 Bayesian Estimation 3.1 Prior Specification Choosing priors is a critical yet challenging task in Bayesian analysis. The selection of the prior distribution in Bayesian statistics typically depends on the study’s particular objectives and the characteristics of the available data. Priors can be informative or non-informative depending on the extent of prior knowledge. In the absence of prior information, we adopt non-informative or hierarchical priors to allow data-driven inference. While the dispersion parameter 𝜙in the Negative Binomial distribution can theoretically be treated as a continuous variable, it is often modeled using a gamma prior due to its flexibility and conjugacy properties (Chen and Lee, 2017; Chen and Khamthong, 2020; Chen et al., 2021; Chu and Yu, 2023). However, the process of selecting hyperparameters for the gamma prior lacks formalized 2 guidance and tends to rely heavily on subjective judgment. In practical applications, 𝜙is frequently interpreted as an integer-valued parameter, which aligns with the discrete nature of count data. Although the Poisson distribution is a natural choice for modeling such data, it does not offer the same adaptability as the gamma distribution in expressing a wide range of prior beliefs about 𝜙. To integrate the flexibility of the gamma distribution and the natural alignment of the Poisson with count data, we construct a hierarchical Poisson-gamma prior for 𝜙as below: 𝜙∼Poisson(𝜌), 𝜌 ∼Gamma(𝜏, 𝜂), with hyperparameters 𝜏and 𝜂chosen to yield vague distributions, allowing it to vary freely across a broad range of plausible values. Notably, as the hyperparameters 𝜏and 𝜂of the gamma prior approach zero, the prior distribution for 𝜙becomes nearly flat, indicating that all feasible values of 𝜙are assigned approximately equal probability. The resulting marginal prior is 𝜋(𝜙) ∝ Γ(𝜙+𝜏) (𝜂+1)𝜙+𝜏+1𝜙!. As 𝜏, 𝜂 →0, the marginal prior for 𝜙becomes inversely proportional to 𝜙, i.e., the prior density behaves like 1/𝜙. Interestingly, this form of inverse proportionality corresponds to the Jeffreys prior for the dispersion parameter in the Negative Binomial distribution. For the state parameters (𝜶and 𝜷), uniform distributions with stationarity constraints are used as priors as discussed in Chen and Lee (2017); Chen and Khamthong (2020); Chen et al. (2021), where the individual components 𝛼𝑖and 𝛽𝑗are treated as independent parameters. Nonetheless, the constraints imposed by the model are not embedded directly in the prior distributions but are instead enforced during the computational phase of inference. We adopt a Dirichlet distribution over the (𝑝+𝑞)-simplex to enforce the stationarity constraint Í𝑖𝛼𝑖+Í𝑗𝛽𝑗<1. This distribution is equivalent to say that (𝜶,𝜷)jointly follow a multivariate uniform distribution over the open standard (𝑝+𝑞)-simplex, and that 𝛼𝑖and 𝛽𝑗(where 𝑖=1,· · · , 𝑝 and 𝑗=1,· · · , 𝑞) each have a marginal distribution of 𝐵𝑒𝑡𝑎(1, 𝑝 +𝑞). For the positive baseline parameter 𝛿, previous studies such as Chen and Lee (2017); Chen and Khamthong (2020); Chen et al. (2021) have adopted a uniform (or flat) prior over the positive real line (0,∞). While this choice is straightforward and non-informative, it has been frequently criticized for being improper, as it does not integrate to a finite value. To address this concern, we instead adopt a hierarchical truncated normal prior, which ensures properness while retaining flexibility in modeling uncertainty about 𝛿. 𝛿∼𝑁(𝜇, 𝜎2)𝐼(0,∞), 𝜇 ∼𝑁(0,105), 𝜎2∼Inv-Gamma(0.1,0.1), with hyperparameters set to produce diffuse distributions. 3.2 Posterior Distribution Let 𝜋(𝜙),𝜋(𝛿), and 𝜋(𝜶,𝜷)denote the prior densities for 𝜙,𝛿, and (𝜶,𝜷), respectively. Assuming independence, the joint posterior is 𝑝(𝜙, 𝜶,𝜷, 𝛿 |𝒀) ∝ 𝐿(𝒀|𝛿, 𝜶,𝜷)𝜋(𝛿)𝜋(𝜶,𝜷)𝜋(𝜙).(2) In theory, a closed-form joint posterior can be obtained by normalizing equation (2). However, due to the function’s complexity, integrating over all parameters is often impractical, so MCMC methods such as Gibbs sampling, implemented via JAGS, are commonly used to remedy this type of difficulties. In Bayesian analysis, the posterior mean is a common choice for point estimation of the parameter, and is known as the Bayes estimator; whereas the highest posterior density (HPD) credible interval is often preferred for interval estimation, especially for asymmetric or multimodal posteriors. This interval is characterized as the shortest interval containing a specified proportion of the posterior 3 probability density. It represents the most credible range of values for the parameter given the observed data and the prior distribution. Stationary parameters such as the mean 𝜇and standard deviation 𝜎are essential for characterizing the long-term behavior of a time series. Leveraging the posterior distribution enables more accurate estimation of these quantities, offering deeper insights into the underlying dynamics. 4 Empirical Example We illustrate the advantages of the proposed Bayesian approach using Ericsson B stock transaction data. The series of transactions per minute on July 2, 2002 (460 observations) has been widely studied in Davis and Liu (2016); Fokianos et al. (2009); Fokianos and Neumann (2013); Doukhan and Kengne (2015); Diop and Kengne (2017) and serves as a benchmark. Figure 1 shows the time series and autocorrelation function (ACF), indicating significant temporal dependence and overdispersion (mean = 9.91, variance = 32.84). We fit a NB-INGARCH(1,1) model: 𝑌𝑡| F𝑡−1∼𝑁𝐵(𝜆𝑡, 𝜙),(3) 𝑃(𝑌𝑡=𝑦𝑡| F𝑡−1)=𝑦𝑡+𝜙−1 𝑦𝑡 𝜙 𝜙+𝜆𝑡𝜙𝜆𝑡 𝜙+𝜆𝑡𝑦𝑡 ,(4) 𝜆𝑡=𝛿+𝛼𝜆𝑡−1+𝛽𝑌𝑡−1,(5) with 𝜙 > 0, 𝛼, 𝛽 ≥0, 𝛿 > 0, and F𝑡−1=𝜎{𝑌𝑡−1,𝑌𝑡−2, . . . }. The model is stationary if 0 < 𝛼 +𝛽 < 1 (Davis and Liu, 2016). As previously discussed in Section 3, we assign a Poisson-Gamma prior for 𝜙with hyperparameters 𝜏=𝜂=0.0001, a Dirichlet(1,1,1) prior for (𝛼, 𝛽), and a hierarchical truncated normal prior for 𝛿. Posterior inference is conducted via Markov Chain Monte Carlo (MCMC) using JAGS in R. Initialization, burn-in, and sampling lengths are automatically managed by functions within the JAGS package, which iteratively extends the simulation until convergence is achieved for all monitored parameters (Gelman and Rubin, 1992). Additional iterations are performed to compensate any autocorrelation. Convergence diagnostics confirm adequate mixing of the chains (Figure 2). The posterior distributions (Figure 3) show that 𝜙peaks near 8, 𝛼exhibits left skewness, while the remaining parameters are right-skewed. Due to these asymmetries, highest posterior density (HPD) intervals are preferred for inference (Gelman et al., 2013; Davidson-Pilon, 2015). A summary of the parameter estimates is listed in Table 1. Par. Mean SD Lower95 Upper95 Median 𝜙7.89 0.97 6.00 10.00 8 𝛼0.83 0.04 0.76 0.90 1 𝛽0.14 0.03 0.09 0.19 0 𝛿0.28 0.16 0.01 0.58 0 Table 1: Ericsson Stock: Bayesian Estimates Summary Statistics One advantage of the Bayesian approach is its effectiveness in estimating the dispersion parameter 𝜙in the NB-INGARCH model. Due to the difficulty associated with the likelihood base estimation method and the seperation nature of the parameters, traditional methods often fix 𝜙and estimate only the state parameters, 𝛼,𝛽, and 𝛿. The Bayesian approach allows joint inference of all parameters in a single procedure. For the dataset considered, setting 𝜙=8 is common (e.g., Davis and Wu, 2009; Davis and Liu, 2016), and the Bayesian estimate aligns with this value. Table 2 presents MLEs for 𝛼,𝛽, and 𝛿, which are similar to the Bayesian estimates. Notably, the classical 95% confidence interval for 𝛿has a negative lower bound, violating its positivity constraint, whereas the Bayesian HPD credible interval respects parameter constraints. 4 Par. Mean SD Lower95 Upper95 𝛼0.84 0.03 0.78 0.91 𝛽0.13 0.03 0.07 0.18 𝛿0.27 0.14 −0.01 0.54 Table 2: Ericsson Stock: Summary Statistics of Classical MLE For the NB-INGARCH(1,1) model, the stationary mean 𝜇and standard deviation 𝜎are 𝜇=𝛿 1−𝛼−𝛽, 𝜎 =v u u t𝜇+𝜇2 𝜙+1+1 𝜙𝛽2𝜇(1+𝜇 𝜙) (1− (𝛼+𝛽)2−𝛽2 𝜙) . The MLEs of the stationary parameters 𝜇and 𝜎2are ˆ𝜇= ˆ 𝛿 1−ˆ𝛼−ˆ 𝛽 =9,and ˆ𝜎=v u u u tˆ𝜇+ˆ𝜇2 ˆ 𝜙+ (1+1 ˆ 𝜙) ˆ 𝛽2ˆ𝜇(1+ˆ𝜇 ˆ 𝜙) (1− ( ˆ𝛼+ˆ 𝛽)2−ˆ 𝛽2 ˆ 𝜙) =5.05. Estimating standard error of the MLE is challenging due to the complex asymptotic properties of the model. In contrast, the Bayesian approach offers a comprehensive assessment of uncertainty, providing both standard errors and credible intervals directly (Bauwens and Lubrano, 1998; Aus´ ın and Galeano, 2007; Ardia, 2008). Table 3 reports the Bayesian point and interval estimates of 𝜇and 𝜎. Par. Mean SD Lower95 Upper95 Median 𝜇10.68 1.63 7.96 14.15 10 𝜎6.64 3.21 4.38 10.73 6 Table 3: Ericsson Stock: Bayesian Estimates of Stationary Parameters 5 Conclusion This paper investigates the NB-INGARCH(𝑝, 𝑞) model, a count time series model characterized by a conditional Negative Binomial distribution, making it particularly suitable for over-dispersed count data. Classical parameter estimation poses challenges due to the complexity of the likelihood function. To address this, we advocate for Bayesian estimation using minimally informative priors, which help mitigate prior influence on the joint posterior and provide a robust and unified approach for inference (Bauwens and Lubrano, 1998; Aus´ ın and Galeano, 2007; Ardia, 2008). Unlike classical two-stage approaches, Bayesian methods allow for the joint estimation of all model parameters, including the dispersion parameter 𝜙. The use of hierarchical priors enhances flexibility, and Dirichlet priors captures the dependency among the parameters. MCMC-based numerical results show that Bayesian highest posterior density (HPD) intervals are more reliable than classical confidence intervals, as they respect model constraints. Furthermore, the Bayesian framework facilitates comprehensive estimation of stationary parameters, such as the mean and variance, by providing both point estimates and credible intervals, which offers deeper insights into the underlying structure of the time series and strengthens the interpretability of model outputs. References Ahmad, A. and Francq, C. (2016). Poisson QMLE of count time series models. Journal of Time Series Analysis, 37(3):291–314. 5 Ardia, D. (2008). Bayesian estimation of the GARCH(1, 1) model with normal innovations. Financial Risk Management with Bayesian Estimation of GARCH Models: Theory and Applications, pages 17–37. Ardia, D. (2009). Bayesian estimation of a Markov-switching threshold asymmetric GARCH model with student-t innovations. The Econometrics Journal, 12(1):105–126. Ardia, D. and Hoogerheide, L. F. (2010). Bayesian estimation of the GARCH(1, 1) model with student-t innovations. The R Journal, 2(2):41–47. Aus´ ın, M. C. and Galeano, P. (2007). Bayesian estimation of the Gaussian mixture GARCH model. Computational Statistics & Data Analysis, 51(5):2636–2652. Bauwens, L. and Lubrano, M. (1998). Bayesian inference on GARCH models using the Gibbs sampler. The Econometrics Journal, 1(1):C23–C46. Chen, C. W. and Khamthong, K. (2020). Bayesian modelling of nonlinear negative binomial integervalued GARCHX models. Statistical Modelling, 20(6):537–561. Chen, C. W. and Lee, S. (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: Applied Statistics, 66(4):797–814. Chen, C. W., Lee, S., and Khamthong, K. (2021). Bayesian inference of nonlinear hysteretic integervalued GARCH models for disease counts. Computational Statistics, 36(1):261–281. Christou, V. and Fokianos, K. (2014). Quasi-likelihood inference for negative binomial time series models. Journal of Time Series Analysis, 35(1):55–78. Chu, Y. and Yu, K. (2023). Bayesian log-linear beta-negative binomial integer-valued GARCH model. Computational Statistics, pages 1–20. Davidson-Pilon, C. (2015). Bayesian Methods for Hackers: Probabilistic Programming and Bayesian Inference. Addison-Wesley. Davis, R. A., Dunsmuir, W. T. M., and Wang, Y. (2000). On autocorrelation in a Poisson regression model. Biometrika, 87:491–505. Davis, R. A., Fokianos, K., Holan, S. H., Joe, H., Livsey, J., Lund, R., Pipiras, V., and Ravishanker, N. (2021). Count time series: A methodological review. Journal of the American Statistical Association, 116(535):1533–1547. Davis, R. A. and Liu, H. (2016). Theory and inference for a class of nonlinear models with application to time series of counts. Statistica Sinica, pages 1673–1707. Davis, R. A. and Wu, R. (2009). A negative binomial model for time series of counts. Biometrika, 96(3):735–749. Diop, M. L. and Kengne, W. (2017). Testing parameter change in general integer-valued time series. Journal of Time Series Analysis, 38(6):880–894. Doukhan, P. and Kengne, W. (2015). Inference and testing for structural change in general Poisson autoregressive models. Electronic Journal of Statistics, 9:1267–1314. Ferland, R., Latour, A., and Oraichi, D. (2006). Integer-valued GARCH process. Journal of time series analysis, 27(6):923–942. Fokianos, K. (2012). Count time series models. In Handbook of statistics, volume 30, pages 315–347. Elsevier. 6 Fokianos, K. and Neumann, M. H. (2013). A goodness-of-fit test for Poisson count processes. Electronic Journal of Statistics, 7:793–819. Fokianos, K., Rahbek, A., and Tjøstheim, D. (2009). Poisson autoregression. Journal of the American Statistical Association, 104(488):1430–1439. Fokianos, K., Støve, B., Tjøstheim, D., and Doukhan, P. (2020). Multivariate count autoregression. Bernoulli, 26(1):471–499. Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian Data Analysis. CRC Press, 3rd edition. Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4):457–472. Jia, Y., Kechagias, S., Livsey, J., Lund, R., and Pipiras, V. (2023). Latent Gaussian count time series. Journal of the American Statistical Association, 118(541):596–606. Kong, J. and Lund, R. B. (2025). Poisson count time series. Journal of Time Series Analysis, 45. Liu, M., Zhu, F., Li, J., and Sun, C. (2023). A systematic review of INGARCH models for integervalued time series. Entropy, 25(6):922. Palmer, W. R., Davis, R. A., and Zheng, T. (2021). Count-valued time series models for COVID-19 daily death dynamics. Stat, 10(1):e369. Weiß, C. H. (2008). Serial dependence and regression of Poisson INARMA models. Journal of Statistical Planning and Inference, 138(10):2975–2990. Zhu, F. (2011). A negative binomial integer-valued GARCH model. Journal of Time Series Analysis, 32:54–67. Zhu, F. (2012). Modeling overdispersed or underdispersed count data with generalized poisson integer-valued GARCH models. Journal of Mathematical Analysis and Applications, 389(1):58– 71. 7 Time Series Plot Time Number of Transactions 0 100 200 300 400 0 10 30 0 5 10 15 20 25 0.0 0.4 0.8 Lag ACF Auto−Correlation Function Plot Figure 1: Time Series and ACF Plots of Ericsson B during July 2nd, 2002 8 Figure 2: Ericsson Stock: MCMC Convergence Diagnostic Plots 9