Bayesian Vector Zero-Inflated Count Model; A Time Series Model for Covid-19 Variables

Bayesian Vector Zero-Inflated Count Model: A Time Series Model for Covid-19 Variables

Amali. E. Victoria, Nwaosu. C. Sylvester, Kuhe A. David, Udomouh . E. Francis, Tamber J. Abraham

Department Physical Sciences, Joseph Sarwuam Tarka University, Makurdi, Benue State.

DOI: https://doi.org/10.51584/IJRIAS.2025.100800130

Received: 13 August 2025; Accepted: 19 August 2025; Published: 22 September 2025

ABSTRACT

This study addresses the challenges of modeling multivariate time series data arising from the COVID-19 pandemic, characterized by overdispersion, zero inflation, and dynamic interdependencies among key indicators such as new cases, death cases, and recoveries. The aim is to develop and evaluate advanced statistical models that can accurately capture these complexities for improved surveillance and forecasting. The methodology adopts a Bayesian framework to estimate both Vector Linear and Log-Linear Zero-Inflated Negative Binomial INGARCH (ZINB-INGARCH) models. The data consist of daily reported COVID-19 cases, deaths, and recoveries from 2020 to 2023. Bayesian Markov Chain Monte Carlo (MCMC) methods were employed to estimate model parameters, using prior distributions tailored for overdispersed and zero-inflated count data. Results showed that the Log-Linear ZINB-INGARCH model performed better than the linear variant in terms of parameter significance, stationarity, and dynamic stability. Particularly, stationarity parameters in the log-linear model (e.g., d₃=−0.507, 95% HDI: [-0.990, -0.020]) demonstrated clear evidence of mean reversion, while autoregressive coefficients such as A[3,3]=0.150, (HDI: [0.074, 0.200]) confirmed strong self-excitation. The dispersion parameters ϕ₁=2.793, ϕ₂=2.903, and ϕ₃=2.116 further confirmed the presence of significant overdispersion. In contrast to the linear model’s higher zero-inflation rates ( ), the log-linear model revealed minimal zero-inflation probabilities around π_i ​≈0.019, indicating better fit to structural zeros in the data. Residual diagnostics supported these results, with Ljung-Box p-values of 1.000 indicating successful modeling of temporal dependencies. The study provides robust empirical evidence that log-linear ZINB-INGARCH models are more effective in capturing the statistical properties of COVID-19 data, offering improved tools for policymakers and epidemiologists in managing future public health crises.

Keywords; : Overdispersion, Zero inflated distribution, Count data, Bayesian inference, Zinb-Ingarch, Covid-19, Multivariate Time Series.

INTRODUCTION

The COVID-19 pandemic has highlighted the urgent need for statistical models capable of handling complex interactions among public health indicators such as new infections, recoveries, and fatalities. Existing multivariate time series models including Multivariate Count Autoregression (Fokianos, 2020), multivariate INGARCH (Heinen et al., 2007), and state-space approaches (Shapovalova et al., 2022) have been extensively applied to discreet multivariate count data. However, their effectiveness diminishes when dealing with discrete data that exhibit overdispersion and a high frequency of zeros, features commonly observed in COVID-19 time series (Zhang, 2022). Time series models like Autoregressive (AR), Moving Average (MA), and ARMA frameworks are fundamentally designed for univariate, continuous, and stationary data, thus limiting their applicability to epidemiological count data, which are inherently non-Gaussian, overdispersed, and often zero-inflated (Zhu et al., 2020).

To better handle such characteristics, time series methodologies have evolved through the integration of Generalized Linear Models (GLMs), giving rise to models such as the Generalized Autoregressive Moving Average (GARMA) model (Benjamin et al., 2003) and GLARMA (Generalized ARMA) a framework of Davis et al., 2002. The Inter-Generalized Autoregressive Conditional Hertroskedasticity (INGARCH) model was introduced by Ferland et al. (2006), with specification based on Poisson assumptions. Subsequent advancements allowed for alternative count distributions, including the negative binomial (Zhu, 2011), generalized Poisson (Zhu, 2012), and generalized compound Poisson (Gonçalves et al., 2023), enabling better accommodation of overdispersion. Despite these developments, standard Poisson and negative binomial-based INGARCH models often fail to adequately model structural zeros, a limitation addressed by Zero-Inflated INGARCH (ZI-INGARCH) models. These models enhance the INGARCH framework by incorporating mixture components to explicitly model excess zeros (Fokianos et al., 2020; Gonçalves et al., 2023). Notably, the Zero-Inflated Generalized Poisson (ZIGP)-INGARCH and Zero-Inflated Negative Binomial (ZINB)-INGARCH models extend this capacity by accommodating a wide range of dispersion scenarios and zero inflation, making them well suited for real-world applications such as COVID-19 case tracking, hospitalization monitoring, and mortality estimation (Zhu, 2012)

LITERATURE REVIEW

Recent advancements in time series modeling have increasingly focused on extending methodologies to multivariate settings to jointly capture overdispersion, zero inflation, and both serial and cross-series interactions critical features in pandemic data involving new infections, deaths, recoveries, vaccine effects and so on. Most research in this area has focused on univariate cases, while multivariate studies have largely been limited to Poisson-related models. Multivariate zero-inflated negative binomial models provide a modern and flexible alternative to Poisson-based approaches, which are constrained by the assumption of equidispersion (Gonçalves et al., 2016; Tsamtsakiri, 2023). Research has evolved along two primary paths: univariate and multivariate approaches. In the univariate case, foundational works of Li (1994) pave way for the GLARMA models proposed by Davis et al. (2001) and GARMA models developed by Benjamin et al. (2003). These were further advanced by Andrade (2017), who applied Bayesian estimation via MCMC and Box-Cox transformations. Theoretical foundations for stationarity were established by Woodard et al. (2011), justifying the use of maximum likelihood estimation in GARMA model of Benjain (et al 2003), the works of Aryuyuen et al. (2014) that propose the Zero Inflated Negative Binomial-Generalized Exponential (ZINB-GE) distribution as a new model for count data with excess zeros and overdispersion. Using both simulated datasets under varying parameter settings and two real-world applications hospital stays among older U.S. residents and household purchases of consumer goods the study evaluates the model’s performance. The ZINB-GE is formulated as a mixture of a Bernoulli process for structural zeros and a Negative Binomial-Generalized Exponential process for counts. Parameters are estimated via maximum likelihood using numerical methods in R. Simulation results show that estimation improves with larger sample sizes, though higher zero proportions increase mean squared error. Empirical findings reveal that the ZINB-GE model fits the real datasets significantly better than traditional ZIP and ZINB models, supported by lower AIC and BIC values and favorable goodness-of-fit tests. The study concludes that the ZINB-GE distribution provides a flexible and effective alternative for modeling overdispersed zero-inflated count data and so on.

In the multivariate case, complex dependence structures have been adopted to better capture inter-series relationships. Soyer et al (2021) set out to review Bayesian approaches for modeling multivariate time series of counts, motivated by the growing need to capture both serial dependence and cross-sectional dependence in applications such as disease surveillance, accidents, and consumer behaviour. They focus on Poisson, negative binomial, INAR, and state-space models, with estimation based on Bayesian inference using MCMC, particle learning, and related methods. Their findings show that Bayesian modeling provides flexible tools for handling overdispersion, zero inflation, and forecasting, though computational challenges remain in high-dimensional settings. Bárdossy et al (2017) focus on drought analysis, aiming to construct multivariate count time series models capable of describing dependence across drought-related variables and over time. Using hydrological and climatological data, they develop models based on INAR processes and multivariate stochastic structures. Bayesian estimation methods are employed, yielding results that capture the persistence and correlation of drought events more realistically than univariate methods. The study concludes that these models improve drought risk assessment and provide valuable insights for water resource management.

Abatan (2013) develops Bayesian methodologies for count time series with the goal of addressing overdispersion, serial dependence, and zero inflation in real-world applications. Using both simulated data and empirical examples, the dissertation advances generalized INAR processes, regression models, and hierarchical state-space formulations. Estimation relies on MCMC algorithms with Gibbs sampling and Metropolis–Hastings, often supported by data augmentation. Results demonstrate robust parameter recovery and strong predictive performance, leading to the conclusion that Bayesian modeling is an effective and coherent framework for analyzing count time series, though improvements in computational efficiency remain an area for future work. Also worthy to mention is the Time-varying correlation models introduced by Tse et al. (2002), while the copula-based models were developed by Heinen and Rengifo (2007) and further extended by Alzaatreh et al. (2022) and Safari-Katesari et al. (2020). Tsionas et al. (2019) introduced Bayesian copula-VARs for high-dimensional settings using MCMC. Comparative evaluations of copula-based and state-space models using particle MCMC and QMLE were presented by Shapovalova et al. (2021). More recently, Zheng et al. (2022) proposed hybrid GARMA-GARCH models, and Qiu et al. (2018) introduced Bayesian Structural Time Series (BSTS) models with spike-and-slab priors. Tsamtsakiri (2023) further advanced this line of work with Bayesian multivariate INGARCH and CARR models using Sarmanov copulas, highlighting persistent computational challenges and the ongoing need for flexible zero-inflated modeling frameworks.

A significant gap emerges in multivariate Negative Binomial INGARCH models capable of handling zero inflation, particularly for applications like COVID-19 data where excess zeros, overdispersion, and complex inter-series dependencies coexist. Existing multivariate count models either neglect zero inflation (e.g., standard INGARCH) or fail to combine it with Negative Binomial distributions’ dispersion flexibility (e.g., copula approaches). The nature of COVID-19 surveillance data with its spatial-temporal patterns and abundance of zero counts (whether from surveillance gaps or true absence of cases) calls for analytical methods that can concurrently account for (1) multivariate dependence via parsimonious INGARCH structures, (2) overdispersion through Negative Binomial marginals, (3) zero inflation via a degenerate mixture component, and (4) computational feasibility for high-dimensional public health surveillance. No current framework integrates all these features, presenting a clear opportunity for research. To address this, the present study applies Bayesian Vector Linear and Log-Linear Zero-Inflated Negative Binomial INGARCH (ZINB-INGARCH) models to COVID-19 data.

MATERIALS AND METHODS

The data consists of daily COVID-19 new cases, death cases, and recoveries in Nigeria. key epidemiological indicators with complex temporal and cross-variable dynamics. The proposed

Zero- Inflated Negative Binomial INGARCH model

The Vector Zero-Inflated Negative Binomial INGARCH model was introduced by Zhu (2011), Chen et al. (2016) examined a Negative Binomial overdispersed model with a time-varying conditional autoregressive mean, treating both the shape and scale parameters as unknown, and estimated them using Bayesian MCMC methods.

\(Y_{i,t} = \left\{ \begin{array}{r}
\pi_{i,t} + \left( 1 – \pi_{i,t} \right).NB\left( 0;r_{i},\lambda_{i,t} \right)\ \ \ \ \ \ \ \ \ \ if\ y = 0 \\
\left( 1 – \pi_{i,t} \right).NB\left( y;r_{i},\lambda_{i,t} \right)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ if\ y > 0
\end{array} \right.\ \) (1)

where in \(\mathbf{Y}_{\mathbf{i,t}}\) a time series, \(\mathbf{r}_{\mathbf{i,t}}\mathbf{=}\frac{\mathbf{r}_{\mathbf{i}}}{\mathbf{r}_{\mathbf{i}}\mathbf{+}\mathbf{\lambda}_{\mathbf{i,t}}}\) ensures\(\mathbf{\ E\lbrack}\mathbf{Y}_{\mathbf{i,t}}\mathbf{\rbrack =}\mathbf{\lambda}_{\mathbf{i,t}}\).

\[P\left( Y_{i,t} = y \right) = \binom{y + r – 1}{y}\left( \frac{r}{\lambda + r} \right)^{r}{\left( \frac{\lambda}{r + \lambda} \right)}^{y}\]

The Mean and Variance for ZINB

\[E\lbrack XY\rbrack = \left( 1 – \pi_{i,t} \right)\mathbf{\lambda}_{\mathbf{i},\mathbf{t}},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Var\lbrack Y\rbrack = \left( 1 – \pi_{i,t} \right)\mathbf{\lambda}_{\mathbf{i},\mathbf{t}}\left( \frac{1}{r} + \pi_{i,t} \right)\]

Hence the conditional density for Linear Vector-INGARCH Model structures is that of VARMA model is represented as

\[\mathbf{Y}_{it}/\mathcal{F}_{t – 1} \sim ZINB(\mathbf{\pi}_{t},\mathbf{r}_{i},\mathbf{\lambda}_{i,t})\ \ \]

\(\mathbf{\lambda}_{\mathbf{t}}\mathbf{\ = \ d\ + \ A}\mathbf{Y}_{\mathbf{t – 1}}\mathbf{\ + \ B}\mathbf{\lambda}_{\mathbf{t – 1}}\) (2)

Given three variables \({\ \ Y}_{1t},\ {\ \ Y}_{2t},\ {and\ \ \ Y}_{t2}\) as New Cases, Death Cases and Recovery respectivelyof COVID-19 variable the

linear model is given as
\[\lambda_{i,t} = d_{i} + \sum_{j = 1}^{3}{A_{i,j}Y_{j,t – i} + \sum_{i = 1}^{3}{B_{i.j}\lambda_{j,t – i}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ i,\ j\ = \ 1,2,3\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3)}\]

\[\lambda_{1,t} = d_{1} + A_{11}Y_{1,t – 1} + A_{12}Y_{2,t – 1} + A_{13}Y_{3,t – 1} + B_{11}\lambda_{1,t – 1} + B_{12}\lambda_{2,t – 1} + B_{13}\lambda_{3,t – 1}\]

\[\lambda_{2,t} = d_{2} + A_{21}Y_{1,t – 1} + A_{22}Y_{2,t – 1} + A_{23}Y_{3,t – 1} + B_{21}\lambda_{1,t – 1} + B_{22}\lambda_{2,t – 1} + B_{23}\lambda_{3,t – 1}\]

\[\lambda_{3,t} = d_{3} + A_{31}Y_{1,t – 1} + A_{32}Y_{2,t – 1} + A_{33}Y_{3,t – 1} + B_{31}\lambda_{1,t – 1} + B_{32}\lambda_{2,t – 1} + B_{33}\lambda_{3,t – 1}\]

Matrix form is given as;

\[\begin{pmatrix}
\mathbf{\lambda}_{\mathbf{1,t}} \\
\mathbf{\lambda}_{\mathbf{2,t}} \\
\mathbf{\lambda}_{\mathbf{3,t}}
\end{pmatrix}\mathbf{=}\begin{pmatrix}
\mathbf{d}_{\mathbf{1}} \\
\mathbf{d}_{\mathbf{2}} \\
\mathbf{d}_{\mathbf{3}}
\end{pmatrix}\mathbf{+}\begin{pmatrix}
\mathbf{A}_{\mathbf{11}} & \mathbf{A}_{\mathbf{12}} & \mathbf{A}_{\mathbf{13}} \\
\mathbf{A}_{\mathbf{21}} & \mathbf{A}_{\mathbf{22}} & \mathbf{A}_{\mathbf{23}} \\
\mathbf{A}_{\mathbf{31}} & \mathbf{A}_{\mathbf{32}} & \mathbf{A}_{\mathbf{33}}
\end{pmatrix}\begin{pmatrix}
\mathbf{Y}_{\mathbf{1,t – 1}} \\
\mathbf{Y}_{\mathbf{2,t – 1}} \\
\mathbf{Y}_{\mathbf{3,t – 1}}
\end{pmatrix}\mathbf{+}\begin{pmatrix}
\mathbf{B}_{\mathbf{11}} & \mathbf{B}_{\mathbf{12}} & \mathbf{B}_{\mathbf{13}} \\
\mathbf{B}_{\mathbf{21}} & \mathbf{B}_{\mathbf{22}} & \mathbf{B}_{\mathbf{23}} \\
\mathbf{B}_{\mathbf{31}} & \mathbf{B}_{\mathbf{32}} & \mathbf{B}_{\mathbf{33}}
\end{pmatrix}\begin{pmatrix}
\mathbf{\lambda}_{\mathbf{1,t – 1}} \\
\mathbf{\lambda}_{\mathbf{2,t – 1}} \\
\mathbf{\lambda}_{\mathbf{3,t – 1}}
\end{pmatrix}\]

The linear model assumes that the linear dependence of \(\mathbf{\lambda}_{\mathbf{i,t}}\) on \(\mathbf{\lambda}_{\mathbf{j,t – 1}}\)and \(\mathbf{Y}_{\mathbf{j,t – 1}}\mathbf{\ \ i,j\ = \ 1,2,3}\) where \({\mathbf{(}\mathbf{A}_{\mathbf{j}}\mathbf{)}}_{\mathbf{q}}\mathbf{,\ \ j = 1,\ \ (}{\mathbf{B}_{\mathbf{i}}\mathbf{)}}_{\mathbf{p}}\mathbf{\ \ i = 1}\) are \(\mathbf{d \times d}\) unknown matrices and all the elements of \(\mathbf{\pi}\), \({\mathbf{(}\mathbf{A}_{\mathbf{j}}\mathbf{)}}_{\mathbf{q}}\) \(\mathbf{j = 1,\ (}{\mathbf{B}_{\mathbf{i}}\mathbf{)}}_{\mathbf{p}}\mathbf{\ i = 1}\) are positive such that\(\mathbf{\ \lambda}_{\mathbf{i,t}}\mathbf{\ }\)> 0. for all i and t, Tsamtsakiri et al. (2023) for the matrix form of the equation.

The Log-linear relationship is given as:

\[\mathbf{Y}_{it}\ |\mathcal{F}_{t – 1} \sim ZINB({\pi_{t},r}_{i},\lambda_{i,t})\ \]

\(\mathbf{v}_{t}\ = \ \mathbf{d}\ + \ \mathbf{A}\mathbf{v}_{t – 1}\ + \ \mathbf{B}{log(\mathbf{Y}}_{t – 1} + \mathbf{1}_{d})\)

\[v_{i,t} = d_{i} + \sum_{j = 1}^{3}\begin{array}{r}
A_{i,j}{log(Y}_{j,t – i} + 1) + \sum_{i = 1}^{3}{B_{i.j}v_{j,t – i}}\ \ \ \ \ \ \ \ \ \ \ i,\ j\ = \ 1,2,3\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (4)\ \ \ \ \ \ \ \ \\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\end{array}\]

where\(\ \nu_{t}\ \equiv \ log\lambda_{t}\) is defined component wise (i.e.\(\mathbf{v}_{i,t} = \ log\ \mathbf{\lambda}_{i,t}\)) and \(\mathbf{1}_{p}\) denotes the p-dimensional vector which consists of ones (Tsamtsakir,.2023). \(\mathbf{d} > 0\) is a p-dimensional vector and A, B are d × d unknown matrices. The elements of d, A and B are assumed to be positive such that \(\mathbf{\lambda}_{i,t} > \ 0\)

\[log(\lambda_{1,t}) = d_{1} + A_{11}log(Y_{1,t – 1} + 1) + A_{12}log(Y_{2,t – 1} + 1) + A_{13}log(Y_{3,t – 1} + 1) + B_{11}log(\lambda_{1,t – 1}) + B_{12}log(\lambda_{2,t – 1}) + B_{13}log(\lambda_{3,t – 1})\ \]

\[log(\lambda_{2,t}) = d_{2} + A_{21}log(Y_{1,t – 1} + 1) + A_{22}log(Y_{2,t – 1} + 1) + A_{23}log(Y_{3,t – 1} + 1) + B_{21}log(\lambda_{1,t – 1}) + B_{22}log(\lambda_{2,t – 1}) + B_{23}log(\lambda_{3,t – 1})\ \]

\[log(\lambda_{3,t}) = d_{3} + A_{31}log(Y_{1,t – 1} + 1) + A_{32}log(Y_{2,t – 1} + 1) + A_{33}log(Y_{3,t – 1} + 1) + B_{31}log(\lambda_{1,t – 1}) + B_{32}log(\lambda_{2,t – 1}) + B_{33}log(\lambda_{3,t – 1})\ \]

therefore the log-linear Vector ZIGP model (Fokianos et.al., 2020) in matrix form is given as;

\(\begin{pmatrix}
\mathbf{log(}\mathbf{\lambda}_{\mathbf{1,t}}\mathbf{)} \\
\mathbf{log(}\mathbf{\lambda}_{\mathbf{2,t}}\mathbf{)} \\
\mathbf{log(}\mathbf{\lambda}_{\mathbf{3,t}}\mathbf{)}
\end{pmatrix}\mathbf{=}\begin{pmatrix}
\mathbf{d}_{\mathbf{1}} \\
\mathbf{d}_{\mathbf{2}} \\
\mathbf{d}_{\mathbf{3}}
\end{pmatrix}\mathbf{+}\begin{pmatrix}
\mathbf{A}_{\mathbf{11}} & \mathbf{A}_{\mathbf{12}} & \mathbf{A}_{\mathbf{13}} \\
\mathbf{A}_{\mathbf{21}} & \mathbf{A}_{\mathbf{22}} & \mathbf{A}_{\mathbf{23}} \\
\mathbf{A}_{\mathbf{31}} & \mathbf{A}_{\mathbf{32}} & \mathbf{A}_{\mathbf{33}}
\end{pmatrix}\begin{pmatrix}
\mathbf{log(Y}_{\mathbf{1,t – 1}}\mathbf{+ 1)} \\
\mathbf{log(}\mathbf{Y}_{\mathbf{2,t – 1}}\mathbf{+ 1)} \\
\mathbf{log(}\mathbf{Y}_{\mathbf{3,t – 1}}\mathbf{+ 1)}
\end{pmatrix}\mathbf{+}\begin{pmatrix}
\mathbf{B}_{\mathbf{11}} & \mathbf{B}_{\mathbf{12}} & \mathbf{B}_{\mathbf{13}} \\
\mathbf{B}_{\mathbf{21}} & \mathbf{B}_{\mathbf{22}} & \mathbf{B}_{\mathbf{23}} \\
\mathbf{B}_{\mathbf{31}} & \mathbf{B}_{\mathbf{32}} & \mathbf{B}_{\mathbf{33}}
\end{pmatrix}\begin{pmatrix}
\mathbf{log(\lambda}_{\mathbf{1,t – 1}}\mathbf{)} \\
\mathbf{log(\lambda}_{\mathbf{2,t – 1}}\mathbf{)} \\
\mathbf{log(}\mathbf{\lambda}_{\mathbf{3,t – 1}}\mathbf{)}
\end{pmatrix}\)

ANALYSIS

Bayesian statistics provides a framework for estimating parameters by updating prior beliefs with observed data. This is done using Bayes’ theorem, which allows us to derive the posterior density of the parameter vector θ given the data y. hence the prior is given as

Prior distribution for the parameter to be estimated \(\phi_{i,t}\sim Gama\ (2,1) = = (\frac{b^{a}}{\Gamma(a)x^{a – 1}})e^{- bx}\)

\(\pi\sim Beta\ (2,2) = \ \frac{1}{B(a,b)}x^{a – 1}{1 – x}^{b – 1}\), \(A_{i}\sim halfN(0,1) = \frac{1}{\sigma\sqrt{2\pi}}e^{(\frac{- x^{2}}{2\sigma^{2}})}\)

\(B_{i}\sim\ halfN(0,1) = \frac{1}{\sigma\sqrt{2\pi}}e^{(\frac{- x^{2}}{2\sigma^{2}})}\), \(d_{i}\sim\ halfN(0,1) = \frac{1}{\sigma\sqrt{2\pi}}e^{(\frac{- x^{2}}{2\sigma^{2}})}\)

The joint likelihood for all time points \(t = 1,\ldots,T\) and series\(i = 1,2,3\) is:

\[L{(Y}_{it}/\theta) = \prod_{t = 1}^{T}{\prod_{i = 1}^{3}{p(Y_{it}/\lambda_{it}(\theta),r_{i}{,\ \pi}_{i})}} = \ \prod_{t = 1}^{T}{\prod_{i = 1}^{3}{\lbrack\pi_{i}^{I\left( Y_{it} = 0 \right)}.\ (1 – \pi)^{I\left( Y_{it} > 0 \right)}.NB(Y_{it}/\lambda_{it},r_{i})}}\]

Posterior Distribution is given as

\[p\left( \Theta\mid Y_{it} \right)\alpha\ p(\theta) \times L{(Y}_{it}/\theta)\]

The Hierarchical structure

Observation layer \(\mathbf{Y}_{t} \sim ZINB(\ \ \lambda_{t},\phi_{k},\ \pi_{k})\), The latent process

\[\mathbf{\lambda}_{t} = exp(\mathbf{\ d\ + \ A}{\log\mathbf{(Y}}_{t – 1}\mathbf{+}\mathbf{1}_{d}\mathbf{)\ + \ B}Log\ \mathbf{\lambda}_{t – 1})\]

\[L(\theta/Y) = \prod_{t = 1}^{3}{{\pi_{it}I}_{y_{kt = 0}} + (1 – {\ \pi}_{it})NB(Y_{i,t}/{\phi_{i,t},\ \pi}_{it},\lambda_{i,t}) = \frac{\Gamma(y + \phi_{i,t})}{\Gamma(\phi_{i,t})y!}\left( \frac{(\phi_{i,t})}{\phi_{i,t} + \lambda_{i,t}} \right)^{\phi_{i,t}}\left( \frac{\lambda_{i,t}}{\phi_{i,t} + \lambda_{i,t}} \right)^{Y_{i,t}})}\]

Posterior distribution

\[P(\theta/y)\ \alpha\ L(\theta/y_{kt}) \times P(A_{k}) \times P(B_{k}) \times P(\phi_{k}) \times P(\pi_{k}) \times P(d_{k})\]

Where \(\theta = (A,B,d,\ \phi,\pi)\)

The posterior Estimate for log linear vector ZINB-INGARCH model is described as fellows

Log-likelihood for MCMC

\[Log\ L = \sum_{t,i}^{}\left\lbrack \mathbb{I}\left( Y_{it} = 0 \right)\log\left( {\ \pi}_{i} + \left( 1 – {\ \pi}_{i} \right)NB(0) + {\mathbb{I(}Y}_{it} > 0 \right)\ Log\left( 1 – {\ \pi}_{i} \right) + \log{NB(Y_{it})}) \right\rbrack\]

For \(Y_{it} = 0\), \(\ \ \ \ \ \ \ \ \ Log(\pi + (1 – \pi) \cdot \exp(\phi\log\phi – \phi\log(\lambda_{t} + \phi))\)

For \(Y_{it} > 0,\ \ \ \ \ \log(1 – \pi) + \log\Gamma(y + \phi) – \log\Gamma(\phi) – \log\Gamma(y + 1) + \phi\log\phi – \phi\log(\lambda_{t} + \phi) + y\log\lambda_{t} – y\log(\lambda_{t} + \phi)\)

Log priors

Log prior for \((A,B,d)\sim N(0,0.5)\)

\(\log P(A,B,d) = – \frac{1}{2}Log(2\pi\ .0.25) – \frac{{(d – 0)}^{2}}{2\ .0.25}\). \(A,B \geq 0\ and\ \ d\mathbb{\ \in R\ }hence\ d\ retain\ normal\ prior\)

For \(\phi\sim\ Gamma(1,1)\), Log P(\(\phi) = – \phi + 0.log\phi\)

For \(\pi\) Beta(1,1), Log P(\(\pi) = 0\)

Log Posterior Distribution

Log\(P(\theta/y)\ = \log L(\theta/y_{kt}) + \log P(A_{k}) + \log P(B_{k}) + \log P(\phi_{k}) + \log P(\pi_{k}) + \log P(d_{k}) + constant\)

\[LogP(\theta/y) = \sum_{t = 1}^{T}{\sum_{k = 1}^{3}{LogP\left( Y_{kt} = y_{kt}/\lambda_{kt},\phi_{k},\ \pi_{k} \right) + \ \sum_{k,i,j = 1}^{3}\left\lbrack Logp\left( A_{ki} \right) + \ Log\left( B_{ki} \right) + logP(d_{k}) \right\rbrack +}}\sum_{k = 1}^{3}{\left\lbrack Log\ p\left( \phi_{k} \right) + logP(\pi_{k}) \right\rbrack + C}\ \]

\[LogP(\theta/y) = \sum_{t = 1}^{T}{\sum_{k = 1}^{3}{\left\{ \begin{array}{r}
Log(\pi + (1 – \pi) \cdot \exp(\phi\log\phi – \phi\log(\lambda_{t} + \phi))\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ if\ y_{kt} = 0 \\
\log(1 – \pi) + \log\Gamma(y + \phi) – \log\Gamma(\phi) – \log\Gamma(y + 1) + \phi\log\phi – \phi\log\left( \lambda_{t} + \phi \right) + y\log\lambda_{t} – y\log\left( \lambda_{t} + \phi \right),ify_{kt} > 0
\end{array} \right.\ \ +}}\sum_{k,i,j = 1}^{3}\left\lbrack \frac{A_{ki}^{2}}{2.0.25} – \log{(0.5\sqrt{2\pi})} \right\rbrack\mathbb{I}_{A_{ki} > 0} + \sum_{k,i,j = 1}^{3}\left\lbrack \frac{B_{ki}^{2}}{2.0.25} – \log{(0.5\sqrt{2\pi})} \right\rbrack\mathbb{I}_{B_{ki} > 0} + \sum_{k,i,j = 1}^{3}\left\lbrack \frac{d_{ki}^{2}}{2.0.25} – \log{(0.5\sqrt{2\pi})} \right\rbrack + \sum_{k = 1}^{3}{\left\lbrack – \phi_{k} + 0.Log\phi_{k} – Log\Gamma(1) \right\rbrack +}\sum_{k = 1}^{3}\lbrack 0\rbrack + C\]

Posterior Sampling via MCMC

We use the No-U-Turn Sampler (NUTS) to draw samples from the posterior:

Step 1: Initialize parameters \(\theta^{(0)} = (\alpha^{(0)},\beta^{(0)},d^{(0)},\phi^{(0)},\pi^{(0)})\)

Step 2: For each MCMC iteration \(\mathbf{s =}1,\ldots,S\)

Step 3: Propose new parameters \(\theta^{*}\) using Hamiltonian dynamics:

Step 4: Simulate trajectories in parameter space, guided by gradients of \(logP(\theta \mid y).\)

Step 5: Accept/Reject \(\theta^{*}\) based on the Metropolis-Hastings ratio:

\[r = \frac{p(\theta^{*} \mid y)}{p(\theta^{(s – 1)} \mid y)}\]

Step 6: Store \(\theta^{(s)} = \theta^{*}\ \) if accepted, else \(\theta^{(s)} = \theta^{(s – 1)}\)

Output: Posterior samples \(\{\theta^{(1)},\ldots,\theta^{(s)}\}.\)

Gradient Calculation:

Compute \(\nabla\theta logp(\theta \mid y)\) for NUTS using automatic differentiation (JAX in Python).

1. Example gradient for α:

\[\frac{\partial logp(\theta \mid y)}{\partial y_{kt}} = \sum_{t = 0}^{T}\frac{\partial logP(Y_{t} = y_{t}/\theta)}{\partial\lambda_{t}}.\ \frac{\partial\lambda_{t}}{\partial y_{kt}}\]

2. Conditional Mean (\(\mathbf{\lambda}_{t}\) ​) Gradient:

\[\frac{\partial\mathbf{\lambda}_{t}}{\partial A} = \frac{\mathbf{\lambda}_{t} \cdot log(1 + y_{t – 1})}{1 + Alog(1 + y_{t – 1}) + Blog(\lambda_{t – 1})}.\]

MCMC (e.g., Random Walk Metropolis) suffers in high dimensions but HMC uses gradient information to propose more efficient moves

HMC Algorithm (Hoffman & Gelman, 2014):

Augment with momentum\(\mathbf{r} \sim N(0,\mathbf{M})\)

Hamiltonian dynamics is given as

\[H(\theta,r) = – log(\theta/\mathbf{Y}) + \frac{1}{2}r^{T}M^{- 1}r\]

Leapfrog integrator (discretized dynamics)

\[\mathbf{r}^{(t + \epsilon/2)} = \mathbf{r}^{(t)} + \frac{\epsilon}{2}\nabla_{\theta}\ log(\theta^{(t)}/\mathbf{Y})\]

\[\theta^{(t + \epsilon)} = \theta^{(t)} + \epsilon\mathbf{M}^{- 1}r^{(t – \epsilon/2)}\]

\[\mathbf{r}^{(t + \epsilon)} = \mathbf{r}^{(t + \epsilon/2)} + \frac{\epsilon}{2}\nabla_{\theta}\ log(\theta^{t + \epsilon)}/\mathbf{Y})\]

Metropolis acceptance:

\[\alpha = min(1,\ exp(\theta^{(t)},\mathbf{r}^{(t)} – H(\theta^{*},\mathbf{r}^{*})))\]

RESULTS

Table 1: Descriptive Statistics of COVID-19 Variables

The descriptive statistics in Table 1 demonstrate that new cases exhibit a mean of (229) substantially higher than the median (79), indicating strong positive skewness. This is further supported by the skewness coefficient (3.1) and high kurtosis (15.2), suggesting a heavily right-skewed distribution with extreme outliers (maximum value of 4006 cases). These patterns reflect periodic spikes in COVID-19 cases, likely corresponding to outbreak waves.

Figure 1: Plot of New Cases of COVID-19 2020 – 2023 from 2020 – 2023

Figure 2: Time Plot of Death Cases of COVID-19 From 2020 – 2023 from 2020 – 2023

Figure 3: Time Plot for Recovery of COVID-19 from 2020 – 2023

Table 3: Posterior Estimate for Vector Linear ZINB INGARCH

HDI (Highest Density Interval), R-hat (̂R) (Gelman-Rubin statistic)

1. Estimated Parameters

\(\phi = \lbrack 0.89,\ \ 2.46,\ \ 0.71\rbrack\) \(\pi = \lbrack 0.48\ ,\ \ \ \ 0.45,\ \ \ 0.46\)]

Matrix form is given as;

\[\begin{pmatrix}
\lambda_{1,t} \\
\lambda_{2,t} \\
\lambda_{3,t}
\end{pmatrix} = \begin{pmatrix}
0.67 \\
0.25 \\
0.12
\end{pmatrix} + \begin{pmatrix}
0.06 & 0.10 & 0.06 \\
0.01 & 0.06 & 0.01 \\
0.09 & 0.11 & 0.08
\end{pmatrix}\begin{pmatrix}
Y_{1,t – 1} \\
Y_{2,t – 1} \\
Y_{3,t – 1}
\end{pmatrix} + \begin{pmatrix}
0.67 & 0.12 & 0.44 \\
0.02 & 0.02 & 0.07 \\
0.45 & 0.10 & 0.21
\end{pmatrix}\begin{pmatrix}
\lambda_{1,t – 1} \\
\lambda_{2,t – 1} \\
\lambda_{3,t – 1}
\end{pmatrix}\]

\[\lambda_{1,t} = \ 0.67\ + \ 0.06 Y_{1,t – 1} + 0.10Y_{2,t – 1}\ + 0.06Y_{3,t – 1}\ + \ 0.67\lambda_{1,t – 1}\ + 0.12\lambda_{2,t – 1}\ + 0.44\lambda_{3,t – 1}\]

\[\lambda_{2,t} = \ 0.25\ + \ 0.01Y_{1,t – 1}\ + 0.06Y_{2,t – 1}\ + 0.01Y_{3,t – 1}\ + \ 0.02\lambda_{1,t – 1}\ + 0.07\lambda_{2,t – 1}\ + 0.02\lambda_{3,t – 1}\]

\[\lambda_{3,t} = \ 0.12\ + \ 0.09Y_{1,t – 1}\ + 0.11Y_{2,t – 1}\ + 0.08Y_{3,t – 1}\ + \ 0.43\lambda_{1,t – 1}\ + 0.10\lambda_{2,t – 1}\ + 0.21\lambda_{3,t – 1}\]

In Table 2, the Vector Linear ZINB-INGARCH model reveals strong short-term dependence for New Cases, especially B[1,1]=0.67 and B[1,3]= 0.44, indicating high influence from its own past and from Recovery. Recovery also shows notable dependence: B[3,1]=0.43 and B[3,3]=0.21. Vector Autoregressive effects are relatively weak across all series. Overdispersion is most in Death Rate (ϕ₂=2.46), suggesting high variability. Zero-inflation is consistently high across series (π₁=0.48, π₂=0.45, π₃=0.46), confirming frequent structural zeros. All R-hat values equal 1.00, reflecting excellent MCMC convergence. Overall, the model captures strong MA dynamics and overdispersion, particularly for Death Rate

Figure 4: Residual Analysis for Linear Vector ZINB INGARCH Model

The Q-Q plots for New Cases, Death Rate, and Recovery shows deviations from normality, particularly in the tails, suggesting non-Gaussian residuals plot. The Residuals vs. Fitted plots show potential heteroscedasticity, with variance increasing at higher fitted values. The ACF plots reveal lingering autocorrelation in New Cases and Death Rate, signaling some unaccounted temporal dependencies. For all variables, ACF decays, but Recovery decays faster, suggesting a better model fit for that series. These findings are supported by Albarracín et al. 2019. Mathew et al. (2022) and However, the results contrast with Andreassen (2013) – who assumed stronger autoregressive influence in multivariate negative binomial frameworks without allowing for ZINB structure or heavy MA dependence.

Table 3: Posterior Estimate for Log-Linear Vector ZINB INGARCH Model

HDI (Highest Density Interval), R-hat (̂R) (Gelman-Rubin statistic)

2. Estimated Parameters

\(\phi = \lbrack 2.793,\ 3.993,\ 2.116\rbrack\), \(\pi = \lbrack 0.019,\ 0.016,\ \ 0.019\)]

\[\begin{pmatrix}
log(\lambda_{1,t}) \\
log(\lambda_{2,t}) \\
log(\lambda_{3,t})
\end{pmatrix} = \begin{pmatrix}
– 0.416 \\
– 0.399 \\
– 0.507
\end{pmatrix} + \begin{pmatrix}
0.051 & 0.051\ & 0.054 \\
0.052 & 0.050 & 0.053 \\
0.055 & 0.055 & 0.056
\end{pmatrix}\begin{pmatrix}
{log(Y}_{1,t – 1} + 1) \\
log(Y_{2,t – 1} + 1) \\
log(Y_{3,t – 1} + 1)
\end{pmatrix} + \begin{pmatrix}
0.211 & 0.161 & 0.219 \\
0.142 & 0.123 & 0.137 \\
0.221 & 0.166 & 0.247
\end{pmatrix}\begin{pmatrix}
{log(\lambda}_{1,t – 1}) \\
{log(\lambda}_{2,t – 1}) \\
log(\lambda_{3,t – 1})
\end{pmatrix}\]

Log-Linear ZINB INGARCH Model Equations

Conditional Mean Equations

\[log\left( \lambda_{1,t} \right) = \ – 0.416\ + \ 0.112log\left( Y_{1,t – 1} + 1 \right) + 0.116log\left( Y_{2,t – 1} + 1 \right) + 0.141log\left( Y_{3,t – 1} + 1 \right) + \ 0.051\ log(\lambda_{1,t – 1})\underbrace{} + 0.051log(\lambda_{2,t – 1})\ + 0.054log(\lambda_{3,t – 1})\]

\[\log\left( \lambda_{2,t} \right) = \ – 0.399\ + \ 0.109\log\left( Y_{1,t – 1} + 1 \right) + 0.115\log\left( Y_{2,t – 1} + 1 \right)\ + 0.136\log\left( Y_{3,t – 1} + 1 \right) + \ 0.052\log\left( \lambda_{1,t – 1} \right) + 0.050\log\left( \lambda_{2,t – 1} \right)\ + 0.053\log\left( \lambda_{3,t – 1} \right)\]

\[log(\lambda_{3,t}) = \ – 0.507\ + \ 0.119log(Y_{1,t – 1} + 1)\ + 0.125log(Y_{2,t – 1} + 1) + 0.150log(Y_{3,t – 1} + 1) + \ 0.055log(\lambda_{1,t – 1}) + 0.055log(\lambda_{3,t – 1}) + 0.056log(\lambda_{3,t – 1})\]

Table 3 shows that parameter estimates reveal distinct dynamic patterns across the three epidemiological series. The diagonal dominance in autoregressive parameters (A[0,0]= 0.112, A[1,1]= 0.115, A[2,2}= 0.150) indicates strong self-excitation effects, consistent with the transmission dynamics described in Held et al. (2005). Cross-component parameters show meaningful, though weaker, interactions—particularly between New Cases and Recovery (A[0,2]= 0.141, B[0,2] = 0.054) mirroring the lagged New cases – Recovery relationships observed by Chen et al. (2021) in their U.S. COVID-19 time-series analysis.

The model’s stationarity was mathematically confirmed through eigenvalue analysis, with the spectral radius \(\rho(\alpha + \beta)\) having a posterior mean of 0.438 (95% HDI: [0.372, 0.512]), satisfying the stationarity and stability conditions for INGARCH-type models derived by Fokianos et al. (2009). The dispersion parameters ϕ exhibited series-specific patterns that align with theoretical expectations. For New Cases and Recovery (ϕ1=2.793, ϕ3=2.116. Collectively, these findings validate the model’s specification against the dispersion spectrum framework proposed by Zhu (2012), which emphasizes tailoring distributional assumptions to match observed variance characteristics in count data.

Table 4: Residual Diagnostics and Model Fit Assessment for log linear ZINB INGARCH Model

Figure 5: Residual Analysis for Log-linear Vector Zigp Ingarch

However, the strong non-normality implies that: Prediction intervals should rely on the model’s native discrete distribution rather than normal approximations. The quadratic approximation used in Laplace-based inference remains appropriate, given the large sample size (n > 1000). Together, these results validate the model’s capacity to represent both the first moment (mean) and second moment (autocovariance structure) while highlighting the necessity of using flexible count distributions for accurate uncertainty quantification and credible interval estimation

DISCUSSION

Based on the comparison of the two models, the Log-Linear Vector ZINB-INGARCH model provides a better overall fit and more reliable parameter estimates than the Vector Linear ZINB-INGARCH model. The log-linear model demonstrates stronger evidence of stationarity, with negative and statistically meaningful d parameters, while the linear model exhibits non-stationarity and wider uncertainty intervals. The autoregressive structure in the log-linear model is more robust, as indicated by its A matrix coefficients having narrow credible intervals that exclude zero, whereas the linear model shows weak or uncertain dependence across many parameters. Additionally, the log-linear model better captures overdispersion in the data through higher and more varied dispersion estimates (ϕ), while its low zero-inflation probabilities (π) suggest that it avoids unnecessary complexity when excess zeros are minimal. Both models show good convergence diagnostics, but the overall statistical and structural performance favors the log-linear formulation. Therefore, the Log-Linear Vector ZINB-INGARCH model is more appropriate and effective for estimating the data.

REFERENCES

1. Abatan, A. O. (2013). Bayesian analysis of count time series models (Doctoral dissertation, University of Sheffield).
2. Aktekin, T., Polson, N., and Soyer, R. (2018). Sequential Bayesian Analysis of Multivariate Count Data. Bayesian Analysis, B(2), 385–409. https://doi.org/10.1214/17-BA1054.
3. Albarracin, O.Y.E., Airlane, P.A. and Ho, L.L. (2019). Generalized Autoregressive and Moving Average Models: Multicollinearity, Interpretation, and a New Modified Model,’ Journal of Statistical Computation and Simulation.
4. Al-Wahsh, H. and Hussein, A. (2020). A Bivariate Autoregressive Poisson Model and its Application to Asthma-related Emergency Room Visits,’ Statistics in Medicine, 39(3) 3184–3194.
5. Alzaatreh, A., Famoye, F., and Lee, C. (2022). Multivariate Count Data Regression Models and Their Applications. In: Advances in Statistical Modeling and Inference. DOI: 10.1007/978-3-031-13971-0_11​.
6. Andrade, B.S., Andrade, M.G. and Ehlers, R.S. (2015). Bayesian GARMA models for count data, Communications in Statistics: Case Studies, Data Analysis and Applications, 1(2), pp. 192–205.
7. Andrade, B.S., (2017). GARMA models: A New Perspective using Bayesian Methods and Transformation.
8. Aryuyuen, S., Bodhisuwan, W., & Supapakorn, T. (2014). Zero inflated negative binomial-generalized exponential distribution and its applications. Songklanakarin Journal of Science and Technology, *36*(4), 483–491.
9. Bárdossy, A., and Pegram, G. G. S. (2017). Multivariate count time series models for drought analysis. Stochastic Environmental Research and Risk Assessment, 31(5), 1349–1362. https://doi.org/10.1007/s00477-016-1240-4
10. Benjamin, M.A., Rigby, R.A. and Stasinopoulos, D.M. (2003). ‘Generalized Autoregressive Moving Average Models’, Journal of the American Statistical Association, 98(461), pp. 214–223.
11. Biswas, A. and Song, P.X.K. (2009) ‘Bayesian forecasting of much count-valued time series,’ Journal of Business and Economic Statistics, 38, pp. 872–887.
12. Bosowski, N. M., (2016). Linear and log-linear models for count time series analysis (Master’s thesis, Northeastern University). https://hdl.handle.net/2047/D20234153
13. Box, G. E. P., Jenkins, G. M., Reinsel, G. C., and Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.
14. Christou, V. and Fokianos, K. (2014) ‘Quasi-likelihood inference for negative binomialtime series models’, Journal of Time Series Analysis, 35, pp. 55–78.
15. Cox, D.R. (1981) ‘Statistical analysis of time series: Some recent developments’, Scandinavian Journal of Statistics, 8, pp. 93–115.
16. Cui, Y. and Zhu, F. (2018) ‘A new bivariate integer-valued GARCH model allowing for negativecross-correlation’, TEST, 27, pp. 428–452.
17. Cui, Y., and Zheng, X. (2020). Bayesian copula selection for mixed data types in high dimensions. Journal of Computational and Graphical Statistics, 29(4), 810-822. https://doi.org/10.1080/10618600.2020.1753538
18. Davis, R. A., Dunsmuir, W. T. M., and Streett, S. B. (2001). Observation Driven Models for Poisson Counts.
19. Davis, R.A., Fokianos, K., Holan, S.H., Joe, H., Livse, J., Lund, R., Pipiras, V. and Ravishanker, N. (2021) ‘Count time series: A methodological review’, Journal of the American Statistical Association, 116, pp. 1533–1547.
20. Davis, R.A., Holan, S.H., Lund, R. and Ravishanker, N. (Eds.) (2016) Handbook of discrete-valued time series. London: Chapman and Hall/CRC.
21. Dabson A.J. and Barnet A.G (2018) An Introduction to Generalized Linear Model, London: Chapman and Hall/CRC.(4)
22. Deng, D., Sun, Y., and Tian. G.-L. (2023). Multivariate zero-inflated binomial model for the analysis of correlated proportional data. Statistical Methods in Medical Research, 32(5), 891–908. https://doi.org/10.1177/09622802231172020
23. Diop, M. L., Diop, A., and Diongue, A. K. (2018). A negative binomial mixture integer-valued GARCH model. Afrika Statistika, 13(2), 1645–1666. https://doi.org/10.16929/as/1645.126
24. Doukhan, P., Fokianos, K., and Tjøstheim, D., (2012). On weak dependence conditions for Poisson autoregressions. Statistics & Probability Letters, 82(5), 942–948.
25. Durrett, R. (2019). Probability: Theory and Examples (5th ed.). Cambridge University Press.
26. Du, J.G. and Li, Y. (1991) ‘The integer-valued autoregressive INAR(p) model’, Journal of Time Series Analysis, 12, pp. 129–142.
27. Durbin, J. and Koopman, S.J. (2000) ‘Time series analysis of non-Gaussian observations based on state-space models from both classical and Bayesian perspectives’, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 62(1), pp. 3–56.
28. Faraway,J.J.,(2016). Extending the Linear Model with R: Generalized Linear, Mixed Effects and Nonparametric Regression Models (2nd ed.). CRC Press
29. Ferland, R., A. Latour, and D. Oraichi (2006). Integer-valued GARCH process. Journal of Time Series Analysis 27(6), 923–942
30. Fokianos, K., Støve, B., Tjøstheim, D. and Doukhan, P. (2020) ‘Multivariate count autoregression’, Bernoulli, 26, pp. 471–499.
31. Fokianos, K. (2021). Autoregressive models for count time series: A survey. Journal of the American Statistical Association, 116(535), 1543–1558.  https://doi.org/10.1080/01621459.2021.1922096
32. Gonçalves, E., N. Mendes-Lopes, and F. Silva (2016). Zero-Inflated Compound PoissonDistributions in Integer-Valued GARCH Models. Statistics 50(3), 558–578.
33. Gill, R. (2022). Bayesian methods for time series analysis. Journal of Time Series Analysis, 43(5), 789-812. https://doi.org/10.1111/jtsa.12642
34. Hardin, J. W., and Hilbe, J. M. (2018). Generalized linear models and extensions (4th ed.). Stata Press.
35. Harvey, A., and Kattuman, P. (2021). Time series models based on growth curves with applications to forecasting coronavirus. Journal of the Royal Statistical Society: Series A.
36. Haugh, M. (2016). The importance of copulas in finance and insurance. Wiley StatsRef: Statistics Reference Online, 1-10. https://doi.org/10.1002/9781118445112.stat07901
37. Heyde, C.C. (1997).Quasi-likelihood and its applications: A general approach to optimal parameter estimation. New York: Springer.
38. Heinen, A. (2003). Modelling time series count data: an autoregressive conditional Poissonmodel. Available at SSRN 1117187
39. Heinen, A. and E. Rengifo (2007). Multivariate autoregressive modeling of time series count data using copulas. Journal of Empirical Finance 14, 564 – 583.
40. Hilbe, J.M. (2014) Modeling count data. Cambridge: Cambridge University Press.
41. Jarque, C.M. and Bera, A.K. (1987) ‘A test for normality of observations and regression residuals’, International Statistical Review, 55, pp. 163–172.
42. Joe, H. (1997) Multivariate models and dependence concepts. London: Chapman and Hall.Kedem, B. and Fokianos, K. (2002) Regression models for time series analysis. Hoboken, NJ: Wiley.
43. Kirchner, M. (2016) ‘Hawkes and INAR(∞) processes’, Stochastic Processes and Their Applications, 126, pp. 2494–2525.
44. Kucharski, A. J., Russell, T. W., Diamond, C., Liu, Y., Edmunds, J., Funk, S., and Eggo, R. M. (2020).Early dynamics of transmission and control of COVID-19: A mathematical modelling study. The Lancet Infectious Diseases, 20(5), 553-558. https://doi.org/10.1016/S1473-3099(20)30144-4
45. Lembert, B. (2016) Student’s guide to Bayesian statistics. SAGE Publications.
46. Lambert, D. (2018). Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics, 34(1), 1-14. (Original work published 1992). https://doi.org/10.1080/00401706.1992.10485228
47. Liu, H. (2012).Some models for time series of counts. Ph.D. thesis, Columbia University, USA.
48. Lütkepohl, H. (2004). Forecasting with VARMA Models. European University Institute, Working Paper ECO No. 2004/25.
49. Lütkepohl, H. (2005) New introduction to multiple time series analysis. Berlin: Springer-Verlag.
50. Mathew J., Bhattacharya s., Sen S., and Das I., (2022). Multivariate Inflated Negative Binimial Regressive for Correlated multivariate count Data. Dependence Modelling (10) 290-307.
51. McCullagh, P. and Nelder, J.A. (1989) Generalized linear models (2nd ed.). London: Chapman and Hall.
52. Mukhopadhyay, S., & Satish, V. (2018). Predictive Likelihood for Coherent Forecasting of Count Time Series.arXiv:1807.02850.
53. Nelder, J.A. and Wedderburn, R.W.M. (1972) ‘Generalized linear models’, Journal of the Royal Statistical Society: Series A (General), 135, pp. 370–384.
54. Qiu, J., Jammalamadaka, S. R., and Ning, N. (2018). Multivariate Bayesian Structural Time Series Model. Journal of Machine Learning Research, 19, 1-33​.
55. Safari-Katesari, H., &Samadi, S. Y. (2020). Modeling count data via copulas: Spearman’s rho and application to cervical cancer. Journal of Applied Statistics, 47(5), 802–820. [DOI]
56. Serhiyenko, V. (2015), Dynamic modeling of multivariate counts: Fitting, diagnostics, and applications. Ph.D. thesis, University of Connecticut, USA.
57. Silveira de Andrade, B. (2017). GARMA Models, a New Perspective Using Bayesian Methodsand Transformations. Doctoral Dissertation, Universidade de São Paulo (USP).
58. Shephard, N. and Pitt, M.K. (1997) ‘Likelihood analysis of non-Gaussian measurement time series’, Biometrika, 84, pp. 653–667.
59. Shapovalova, Y., Baştürk, N., and Eichler, M. (2021).Multivariate Count Data Models for Time Series Forecasting. Entropy, 23(718). https://doi.org/10.3390/e23060718&#8203}.
60. Smith, M.S. and Khaled, M.A. (2012) ‘Estimation of copula models with discrete margins via Bayesian data augmentation’, Journal of the American Statistical Association, 107, pp. 290–303.
61. Shumway, R. H., and Stoffer, D. S. (2017).Time Series Analysis and Its Applications: With R Examples (4th ed.). Springer.
62. Soyer, R., and Zhang, D. (2021). Bayesian modeling of multivariate time series of counts. Wiley Interdisciplinary Reviews: Computational Statistics, 13(4), e1559. https://doi.org/10.1002/wics.1559
63. Tsionas, M. G., and Kumbhakar, S. C. (2019). Bayesian estimation of large time-varying VARs using copulas. Journal of Econometrics, 210(2), 286–306. [DOI]
64. Tsay, R.S. (2014) Multivariate time series analysis. Hoboken, NJ: Wiley.
65. Tse, Y. K. and Tsui A. K. C. (2002). A multivariate generalized autoregressive conditional heteroscedasticity model with time-varying correlations. Journal of Business & Economic Statistics 20(3), 351–362
66. Tsamtsakiri, P. (2023). Modeling Multivariate Time Series. Doctoral Dissertation, Athens University of Economics and Business​.
67. Wang, X., D. Wang, and H. Zhang (2020). Poisson autoregressive process modeling via the penalized conditional maximum likelihood procedure. Statistical Papers 61(1), 245–260
68. Wang, C., Liu, H., Yao, J.-F., Davis, R.A. and Li, W.K. (2014) ‘Self-excited threshold Poisson autoregression’, Journal of the American Statistical Association, 109, pp. 777–787.
69. Wang, F. and H. Wang (2018). Modelling non-stationary multivariate time series of counts via common factors. Journal of the Royal Statistical Society: Series B 80, 769–791.
70. Wood, S. N. (2017).Generalized Additive Models: An Introduction with R (2nd ed.). CRC
71. Woodard, D.W., Matteson, D.S. and Henderson, S.G. (2011) ‘Stationarity of count-valued and nonlinear time series models’, Electronic Journal of Statistics, 5, pp. 800–828.
72. Wu, H. and Peirus, S. (2018) ‘An introduction to vector Gegenbauer processes with long memory’, The ISIs Journal for the Rapid Dissemination of Statistics Research, 7.
73. Wu, W.B. and Shao, X. (2004) ‘Limit theorems for iterated random functions’, Journal of Applied Probability, 41, pp. 425–436.
74. Yiu, K.T. and Albert, K.C.T. (2022) ‘A multivariate generalized autoregressive conditional heteroscedasticity model with time-varying correlation’, Journal of Business and Economic Statistics, 20(3), pp. 351–362.
75. Yang, L., Frees, E.W. and Zhang, Z. (2020) ‘Nonparametric estimation of copula regression models with discrete outcomes’, Journal of the American Statistical Association, 115, pp. 707–720.
76. Jung, R. C., &Liesenfeld, R. (2006). “Zero-Inflated Count Time Series Models with Applications to High-Frequency Financial Data.” Journal of Empirical Finance, 13(1), 23–39.
77. Zeger, S.L. (1988) ‘A regression model for time series of counts’, Biometrika, 75, pp. 621–629.
78. Zeger, S. L., Liang, K.-Y., & Albert, P. S. (1988).Models for Longitudinal Data: A Generalized Estimating Equation Approach. Biometrics, 44(4), 1049-1060.
79. Zeger, S.L. and Qaqish, B. (1988) ‘Markov regression models for time series: A quasi-likelihood approach’, Biometrics, 44, pp. 1019–1031.
80. Zheng, T., Xiao, H. and Chen, R. (2015) ‘Generalized ARMA models with martingalizeddifference errors’, Journal of Econometrics, 189, pp. 492–506.
81. Zheng, X., Chen, Y., and Li, D. (2021). GARMA-GARCH models for non-Gaussian time series with conditional heteroskedasticity. Journal of Econometrics, 224(2), 352-370. https://doi.org/10.1016/j.jeconom.2021.05.003
82. Zheng, T., Xiao, H., & Chen, R. (2022). Generalized Autoregressive Moving Average Models with GARCH Errors. Journal of Time Series Analysis, 43(125–146).
83. Zhang, Y., Zhou, H., Zhou, J. and Sun, W. (2017) ‘Regression models for multivariate count data’, Journal of Computational and Graphical Statistics, 26, pp. 1–13.
84. Zhang, M. (2022). Copula-based Analysis for Dependent Count Data. Ph.D. Dissertation, The George Washington University​
85. Zhu, F., and Wang, G. (2020). “Negative Binomial INGARCH Models for Overdispersed Count Time Series.” Journal of Time Series Analysis, 41(3), 343–361.
86. Zhu, X., Pan, R., Li, G., Liu, Y. and Wang, H. (2017) ‘Network vector autoregression’, The Annals of Statistics, 45, pp. 1096–1123.