Abstract

Various models for time series of count data account for discreteness, overdispersion and serial dependence. In addition to these, accounting for covariates incorparation pattern are the complexities that arise while dealing with data which involve seasonality aspects. Specifying a model that can handle such kind of time series of count data is very important in several real-life application. However in this paper, we present a non-stationary autoregressive model where covariates information are incorporated in the Branching Negative Binomial (bNB) autoregressive model in order to assess the seasonality in the process of time series event. A simulation study is done to evaluate how well the proposed strategy performs, and inference is based on maximum likelihood estimate. The model is used to analyse a real-world dataset, which is an infectious disease with covariates information, including temperature and rainfall.

Keywords

Non-Stationary Time Series, Covariates, Branching Process, Time Series of Count Data, Overdispersion

Share and Cite:

1. Introduction

Various models for time series of count data account for discreteness, overdispersion and serial dependence. Among those applications that require these aspects, road accident counts (see e.g. [1]), sandstorm counts (see e.g. [2]), stock market application (see e.g. [3]; [4]), crime analysis (see e.g. [5]; [6]), among others. Several models have been presented by scholars to handle the analysis of such kind of data. Those models include Zero truncated Poison INAR, INARGH, NGINAR, GLM, Zero-modified geometric INAR, bNB, INAR etc., (see [7]-[10]). In addition, complexities arise while specifying a modeling strategy that is able to deal with count time series data which are not stationary. Some of the nonstationary data involve trend, seasonal, cyclical patterns. It is common to encounter time series of counts that are observed together with covariate information. These kinds of data commonly arise in fields of application such as epidemiology, accident analysis, ecology, forest fires among others. A large number of studies revealed that, these nonstationary aspects are due to some natural phenomena in time series event process such as temperature, rainfall, whether, humidity, etc. In particular, weather variables may affect accident counts (see e.g. [11]), environmental factors influence species counts recorded over time (see e.g. [12]), occurrence of wildfires over time is influenced by meteorological variables (see e.g. [13]), etc. In such case non-stationarity needs to be incorporated into the modeling framework as well as mechanisms to capture dependence and other distributional characteristics. There are various modeling strategies that have been proposed to deal with the issues arising when handling time series of counts observed with covariate information.

Studying the effect of weather conditions on daily crash counts using a discrete time-series model is proposed in this paper [14]. The model was applied to daily car crash data, metereological data and traffic exposure data from the Netherlands to analyze the risk impact of weather conditions on the observed counts. According to the authors of [15], count data time series may have non-stationarity from trends or covariates. While modelling count data time series with Markov processes based on binomial thinning, they proposed an extension of stationary time series based on binomial thinning such that the univariate marginal distributions are always in the same parametric family, such as negative binomial. In order to conduct a statistical analysis of data collected by the Massachusetts Water Resources Authority (MWRA) between 1996 and 2002 to evaluate the effects of court-mandated improvements in sewage treatment, [16] proposed a nonstationary negative binomial time series with time-dependent covariates. Utilizing financial and economic variables as exogenous covariates, [17] developed a class of Poisson autoregressive models with exogenous covariates (PARX) that can be used to model and forecast monthly corporate defaults in the US from 1982 to 2011 using these variables as time series attributes. The results showed the effectiveness of their model in capturing the time series dynamics of corporate defaults, including the well-known clustering of default counts in the data. Furthermore, they discovered that while generally speaking, current defaults do affect the likelihood that other firms will experience default in the future, in recent years, economic and financial factors at the macro level have been able to explain a significant portion of the correlation of US firm defaults over time. [18] presented parameter time series models with covariates for integer Z. The used real data applications on accidents and financial returns are given. [19] presented a Bayesian analysis of a time series of counts to assess its dependence on an explanatory variable. As application, they analyse the incidence of an infectious disease with an explanatory variable which is he number of grams of antibiotic (third generation) cephalosporin. Statistically, they realized that there is significant relationship between disease occurrence and use of the antibiotic, lagged by three months. A popular approach entails developing regression models by introducing covariates in the marginal distribution, dependence component or both. Typical examples where the parameters are models as a function of covariates include Integer Autoregressive and Generalized Linear Modeling frameworks. However, we extend the Branching Negative Binomial Autoregressive model presented by the authors in [10] to a NonStationary model, where covariate information is incorporated to the dependence parameter.

The rest of the paper is organized as follows; Section 2 outlines the approach adopted in developing the non-stationary branching negative binomial (bNB) autoregressive model as well the procedure for inference, Section 3 provides the results for the simulation of stationary bNB autoregressive model, Section 4 presents the results from the application of the model to real-life data, Section 5 performs the goodness-of-fit tests of models. Finally, Section 6 concludes the work.

2. Methods and Materials

The existing stationary model presented in [10] is defined by

$Y_t=N_{t-1}+E_t$ (1)

where

$N_{t-1}~Bin\left(Y_{t-1},\frac{\rho \pi }{1-\rho +\rho \pi }\right)\text{and}{E}_t~NB\left(m+N_{t-1},\frac{\pi }{1-\rho +\rho \pi }\right)$

They denote the model as $bNB\left(m,\pi ,\rho \right)$ .

The branching concept is that, since $N_{t-1}$ is generated from $Y_{t-1}$ , which is considered as previous observation. The generating process follows binomial law with probability $\frac{\rho \pi }{1-\rho +\rho \rho }$ . In terms of dependence parameter $\rho$ , it’s the rate to determine how many offsprings are generated from the previous $Y_{t-1}$ which also specifies the dependence in the autoregressive process from $Y_t$ .

The marginal distribution of the arising process is Negative Binomial because the error component $E_t$ has a Negative Binomial distribution.

In this framework, we extend the model (Equation 1) proposed in [10] by incorporating covariates information in the dependence parameter to obtain a non-stationary model capable to handle seasonal pattern in times series of count data.

2.1. Model Specification

Incorporating covariates into the dependence parameter, we considered the logit function in order to preserve the nature of the parameter ($\rho \in \left[0,1\right]$ ) from the stationary model. Hence, it will be a time varying function, denoted by ${\rho }_t$ ; and it follows

${\rho }_t=\frac{1}{1+\text{exp}\left(-X^{\text{T}}\beta \right)}$ (2)

where $\beta =\left(\begin{array}{l}{\beta }_0\hfill \\ {\beta }_1\hfill \\ {\beta }_2\hfill \\ ⋮\hfill \\ {\beta }_p\hfill \end{array}\right)$ , $X=\left(\begin{array}{l}1\hfill \\ X_1\hfill \\ X_2\hfill \\ ⋮\hfill \\ X_p\hfill \end{array}\right)$

where ${\beta }_1,{\beta }_2,\cdots ,{\beta }_p$ are regression coefficients and $x_1,x_2,\cdots ,x_p$ the covariates.

The derived non-stationary $bNB$ model is given in (Equation 3) with parameters $\theta =\left(m,\pi ,\beta \right)$ .

$Y_t=N_{t-1}+E_t$ (3)

where

$N_{t-1}~Bin\left(Y_{t-1},\frac{{\rho }_t\pi }{1-{\rho }_t+{\rho }_t\pi }\right)\text{and}{E}_t~NB\left(m+N_{t-1},\frac{\pi }{1-{\rho }_t+{\rho }_t\pi }\right)$

We denote the obtained model as $bNB\left(m,\pi ,{\rho }_t\right)$ for ${\rho }_t$ defines in (Equation 2).

In this process (Equation 1), distribution of the two components are given bellow;

$P\left(N_{t-1}=k\right)=\left(\begin{array}{c}{Y}_{t-1}\\ k\end{array}\right){\left(\frac{{\rho }_t\pi }{1-{\rho }_t+{\rho }_t\pi }\right)}^k{\left(1-\frac{{\rho }_t\pi }{1-{\rho }_t+{\rho }_t\pi }\right)}^{Y_{t-1}-k}$ (4)

$P\left(E_t=r\right)=\left(\begin{array}{c}\left(m+N_{t-1}\right)+r-1\\ r\end{array}\right){\left(\frac{\pi }{1-{\rho }_t+{\rho }_t\pi }\right)}^{\left(m+N_{t-1}\right)}\times {\left(1-\frac{\pi }{1-{\rho }_t+{\rho }_t\pi }\right)}^r$ (5)

2.2. Model Estimation

2.2.1. Likelihood Function

Let $T=\left\{t_0<t_1<\cdots <t_n\right\}$ be an increasing set of times and $y=\left\{y_0,y_1,\cdots ,y_n\right\}\subset \mathbb{N}$ an arbitrary set of non-negative integers with covariates $X=\left\{X_1,X_2,\cdots ,X_n\right\}\in ℝ$ . Then the joint pmf for $Y~bNB\left(m,\pi ,{\rho }_t\right)$ is the product of the marginal and the conditionals.

$P\left[Y=y|X,m,\pi ,\beta \right]=p\left(y_o\right){\displaystyle \prod }_{t=1}^{n}P\left(y_t|x_t,y_{t-1}\right)$ (6)

where $m>0$ , $\pi \in \left[0,1\right]$ , and $\beta \in {ℝ}^{N+1}$ .

Conditional or transition probability mass function (pmf)

$\begin{array}{c}p\left(y_t|y_{t-1}\right)=P\left[Y_t=y_t|Y_{t-1}=y_{t-1}\right]\\ ={\displaystyle \sum _{N_{t-1}=0}^{\infty }}\left(\begin{array}{c}{y}_{t-1}\\ N_{t-1}\end{array}\right){\left(\frac{{\rho }_t\pi }{1-{\rho }_t+{\rho }_t\pi }\right)}^{N_{t-1}}\times {\left(\frac{1-{\rho }_t}{1-{\rho }_t+{\rho }_t\pi }\right)}^{y_{t-1}-N_{t-1}}\\ \times \left(\begin{array}{c}m+N_{t-1}+E_t-1\\ E_t\end{array}\right){\left(\frac{\pi }{1-{\rho }_t+{\rho }_t\pi }\right)}^{m+N_{t-1}}\times {\left(\frac{\left(1-{\rho }_t\right)\left(1-\pi \right)}{1-{\rho }_t+\pi }\right)}^{E_t}\end{array}$

$\begin{array}{c}=\frac{y_{t-1}!\Gamma \left(m+y_{t-1}\right){\pi }^m{\left(1-\pi \right)}^{y_{t-1}}{\left(1-{\rho }_t\right)}^{y_t+y_{t-1}}}{{\left(1-{\rho }_t+{\rho }_t\pi \right)}^{m+y_{t-1}+y_t}}\\ \times {\displaystyle \sum _{N_{t-1}=0}^{\infty }}\{\frac{1}{N_{t-1}!\left(y_t-N_{t-1}\right)!\left(y_{t-1}-N_{t-1}\right)!\Gamma \left(m+N_{t-1}\right)}\\ \times {\left[\frac{{\rho }_t{\pi }^2}{{\left(1-{\rho }_t\right)}^2\left(1-\pi \right)}\right]}^{N_{t-1}}\}\\ =\frac{\Gamma \left(m+y_t\right)}{\Gamma \left(m\right)y_t!}\frac{{\pi }^m{\left(1-{\rho }_t\right)}^{y_t+y_{t-1}}{\left(1-\pi \right)}^{y_t}}{{\left(1-{\rho }_t+{\rho }_t\pi \right)}^{m+y_t+y_{t-1}}}\\ \times {}_2{F}_1\left(-y_{t-1},-y_t;m;z\right)\end{array}$ (7)

where $z=\frac{{\rho }_t{\pi }^2}{{\left(1-{\rho }_t\right)}^2\left(1-\pi \right)}$ and ${}_2{F}_1\left(a,b;c;k\right)$ is Gauss’ hypergeometric function with parameters $a,b,c$ and $k$ .

Marginal distribution

The marginal distribution for $y_o$ is:

$p\left(y_o\right)=\frac{\text{Γ}\left(m+y_o\right)}{\text{Γ}\left(m\right)y_o!}{\pi }^m{\left(1-\pi \right)}^{y_o}$ (8)

Substituting (Equation 7) and (Equation 8) in (Equation 6), the likelihood is defined as follow

$\begin{array}{c}P\left[Y=y,X|m,\pi ,\beta \right]\\ =\frac{\Gamma \left(m+y_o\right)}{\Gamma \left(m\right)y_o!}{\pi }^m{\left(1-\pi \right)}^{y_o}\times {\displaystyle \prod _{t=1}^n}\left\{\frac{\Gamma \left(m+y_t\right)}{\Gamma \left(m\right)y_t!}\frac{{\pi }^m{\left(1-{\rho }_t\right)}^{y_{t-1}+y_t}{\left(1-\pi \right)}^{y_t}}{{\left(1-{\rho }_t+{\rho }_t\pi \right)}^{m+y_{t-1}+y_t}}\right\}\\ \underset{}{\overset{}{}}{\times }_2{F}_1\left(-y_{t-1},-y_t;m;z\right)\}\\ ={\displaystyle \prod _{t=0}^n}\left\{\frac{\Gamma \left(m+y_t\right){\pi }^m}{\Gamma \left(m\right)y_t!}\right\}\times \frac{{\left(1-\pi \right)}^{{\displaystyle {\sum }_{t=0}^n{y}_t}}{\left(1-{\rho }_t\right)}^{{\displaystyle {\sum }_{t=1}^n\left(y_{t-1}+y_t\right)}}}{{\left(1-{\rho }_t+{\rho }_t\pi \right)}^{nm+{\displaystyle {\sum }_{t=1}^n\left(y_{t-1}+y_t\right)}}}\\ \times {\displaystyle \prod _{t=1}^n}{}_2{F}_1\left(-y_{t-1},-y_t;m;z\right)\\ ={\displaystyle \prod _{t=0}^n}\left\{\frac{\Gamma \left(m+y_t\right){\pi }^m}{\Gamma \left(m\right)y_t!}\right\}{\displaystyle \prod _{t=1}^n}{}_2{F}_1\left(-y_{t-1},-y_t;m;z\right)\\ \times \frac{{\left(1-\pi \right)}^{y_{+}}{\left(1-{\rho }_t\right)}^{2y_{+}-y_o-y_n}}{{\left(1-{\rho }_t+{\rho }_t\pi \right)}^{nm+2y_{+}-y_o-y_n}}\end{array}$ (9)

where $y_{+}={\displaystyle \sum }_{t=0}^n{y}_t$ , ${\rho }_t=\frac{1}{1+\text{exp}\left(-X^{\text{T}}\beta \right)}$ , $z=\frac{{\rho }_t{\pi }^2}{{\left(1-{\rho }_t\right)}^2\left(1-\pi \right)}$

2.2.2. Log-Likelihood Function

Taking log (natural logarithm) of the previous equation (Equation 9), we obtain the log-likelihood (Equation 10)

$\begin{array}{l}\log L\left(m,\pi ,\beta |Y=y\right)\\ ={\displaystyle \sum }_{t=0}^n\log \Gamma \left(m+y_t\right)+\left(n+1\right)m\log \pi +{\displaystyle \sum }_{t=1}^n{\log }_2{F}_1\left(-y_{t-1},-y_t;m;z\right)\\ -\left(n+1\right)\log \Gamma \left(m\right)+y_{+}\log \left(1-\pi \right)+\left(2y_{+}-y_o-y_n\right)\log \left(1-{\rho }_t\right)\\ -\left(nm+2y_{+}-y_o-y_n\right)\log \left(1-{\rho }_t+{\rho }_t\pi \right)\end{array}$ (10)

The parameters to be estimated are the marginal, the dependence and regression parameters. The parameter vector is $\theta =\left(m,\pi ,{\beta }_0,\cdots ,{\beta }_p\right)$ , which are obtained by optimizing the log-likelihood function using the optim() routine in R software.

In statistical modeling, the maximum likelihood estimation (MLE) method is commonly used to estimate the parameters of a probability distribution based on a given set of data. The MLE method involves finding the values of the parameters that maximize the log-likelihood function, which is the logarithm of the probability of observing the data given the parameter values.

However, in this case it is not easy to differentiate the log-likelihood function with respect to every parameter because of the presence of special functions such as the Gauss hypergeometric function. The derivative of the log-likelihood function with respect to a parameter may not have a closed-form solution, making it challenging to find the parameter value that maximizes the log-likelihood function.

Furthermore, the presence of the Gamma function in the expression of the derivative of the log-likelihood function can make it difficult to obtain a closed-form solution for the derivative of the log-likelihood function with respect to the parameter. The derivative of the Gamma function involves the digamma function, which does not have a closed-form solution in general.

Therefore, when there is a difficulty in deriving the log-likelihood function with respect to every parameter because of the Gauss hypergeometric function in the expression and its derivative has difficulty exhibiting a parameter because of the Gamma function, we can use numerical optimization techniques such as the optim() function in R software to estimate the parameters that maximize the log-likelihood function. The optim() function uses numerical methods to find the minimum or maximum of a given function, which in this case is the negative log-likelihood function.

By using the optim() function in R, we can estimate the parameter values that maximize the log-likelihood function, even when closed-form solutions for the derivative of the log-likelihood function with respect to the parameter are not available. This allows us to perform statistical inference and make predictions based on the estimated parameter values.

In the next section, the R software package is used in the simulation process to estimate parameters of the non-stationary $bNB$ models.

3. Simulation Study

In this simulation study, our target is bear on one covariate variable where $X=\left(1,X_1\right)$ and $\beta =\left({\beta }_0,{\beta }_1\right)$ . Base on this restriction, the Monte Carlo experiments have been conducted several time to evaluate the performance of the Branching Negative Binomial estimators for the parameters. The experiment entails generating $n$ time series of size $N$ from the models and then estimating the parameter vector $\theta =\left(m,\pi ,{\beta }_0,{\beta }_1\right)$ . For each combination ($\theta$ , $N$ ), we compute the mean, bias, and the mean-squared error (MSE) given by,

$\widehat{\theta }=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\widehat{\theta }}_i,Bias\left(\widehat{\theta }\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left({\widehat{\theta }}_i-\theta \right)\text{and}MSE\left(\widehat{\theta }\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left({\widehat{\theta }}_i-\theta \right)}^2$

where ${\widehat{\theta }}_i$ is the estimated parameter vector values from the $i^{th}$ simulated series.

Using the software package R, we generate by Monte Carlo replicates n = 1000 random samples of length N = 250, N = 500 and N = 1000 from the non-stationary bNB, then calculate the mean, bias and mean-squared error of the estimator.

Simulation Results

To assess the performance of the proposed approach of incorporating covariates in terms of estimation, we focused on the case in which the covariates are introduced in the dependence parameter as it is more promising in terms of accommodating non-stationary patterns in both the marginal distribution $m$ and the dependence parameter, $\rho$ . We take motivation from epidemiology and biostatistics where harmonic functions are commonly employed within regression models as standard tools for capturing seasonal patterns (Equation 11) where we consider sinusoidal term $X=\text{sin}\left(2\pi t/12\right)$ as covariates information in this simulation. Then the corresponding time dependence correlation parameter is given by

${\rho }_t=\frac{1}{1+\exp \left(-{\beta }_0-{\beta }_1\sin \left(2\pi t/12\right)\right)}$ (11)

Figure 1. Time series plot and autocorrelation functions from a sample size of 200 simulated data with parameter $m=15$ , $\pi =0.6$ , ${\beta }_0=0.7$ and ${\beta }_1=0.8$ .

Figure 1 shows plots of time series data generated from the nonstationary bNB model and the corresponding autocorrelation functions. Particularly, we display seasonal effects by considering parameters $m=15$ , $\pi =0.6$ ,${\beta }_0=0.7$ and ${\beta }_1=0.8$ .

Table 1. Summary statistics for the estimator for different parameter values $\theta =\left(m,\pi ,{\beta }_0,{\beta }_1\right)$ and different sample sizes $N$ . These statistics are obtained from 1000 Monte Carlo replications of the developed model.

Table 2. Summary statistics for the estimator for different parameter values $\theta =\left(m,\pi ,{\beta }_0,{\beta }_1\right)$ and different sample length of the series $N$ . These statistics are obtained from 1000 Monte Carlo replications of the non-stationary developed model by varying regression parameters.

Figure 2. Boxplots with results from the simulation experiment for $N=250$ .

Figure 3. Boxplots with results from the simulation experiment for $N=1000$ .

Figure 4. Boxplots with results from the simulation experiment for $N=250$ .

In Table 1 and Table 2, we report the mean, bias and mean-squared error with parameters $\theta =\left(m,\pi ,{\beta }_0,{\beta }_1\right)$ obtained using 1000 Monte Carlo replications of time series of length 250, 500 and 1000. As expected, the results indicate that the biases decrease as the sample size increases. We realized that in this scenario, when the parameter $m$ varies by increasing, its mean square error increases while for rest of parameters decrease; and its bias increase while for the rest of parameters decrease. Also, when parameter $\pi$ varies by increasing, its bias decreases with mean square errors of ${\beta }_0$ and ${\beta }_1$ . These results are also illustrated in the boxplots via Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7. All these show that there is an improvement in the results as the sample size increases.

Figure 5. Boxplots with results from the simulation experiment for $N=1000$ .

Figure 6. Boxplots with results from the simulation experiment for $N=250$ .

Figure 7. Boxplots with results from the simulation experiment for $N=1000$ .

The following figures are the boxplots of the simulation experiment by using parameters from Table 1, where the black line is the median and the red is true parameter.

4. Application to Real-Life Data

The purpose of this application is to assess : (1) the effect of temperature, and (2) the effect of rainfall on the Dengue fever count series. Therefore we applied the nonstationary Branching Negative Binomial model to the weekly data from Thailand over the period January 2006 to February 2013, which is 373 observations in total.

Table 3 is the summary statistic of the temperature data and Figure 8 is its time series plot. The temperature is higher in the middle of the year while lower between the end and the beginning of two consecutive years, respectively. In particular, higher from April to September, and lower from October to March.

Table 3. Summary statistics for the weekly mean temperature of Thailand from January 2006 to February 2013.

Figure 8. Weekly plot of Thailand temperature over the period of January 2006 to October 2013.

Here, Table 4 displays the summary statistic of the rainfall data and Figure 9 is its time series plot. The rainfall quantity seems to increase in the middle of the year while decreases between the end and the beginning of two consecutive years, respectively. In particular, higher from April to September, and lower from October to March.

Table 4. Summary statistics for the weekly rainfall mean of Thailand from January 2006 to February 2013.

Table 5 is the summarize statistic of the Dengue Fever Counts data and Figure 10 is its time series plot. The variance and the mean are respectively 17.68 and 3.20, which implies the existence of overdispersion in the data since the variance is greater than the mean. This is one of the requirements of our developed models. So far, in addition to the Seasonality pattern in the data (Figure 10), the partial autocorrelation functions (Figure 11) displayed that data is autoregressive of order one. Therefore our model is suitable to apply to this particular data.

Figure 9. Weekly plot of Thailand rainfall volume over the period of January 2006 to February 2013.

Table 5. Summary statistics for the weekly mean temperature of Thailand from January 2006 to February 2013.

Figure 10. Time series plot of weekly Dengue Fever counts data from Thailand, over the period of January 2006 to February 2013.

Figure 11. Plot of autocorrelation functions.

This application is to fit our developed model and the seasonality in Dengue Fever counts data, where temperature is an effect factor (Figure 12). By fitting the Data to the nonstationary bNB model, we obtained the parameters in Table 6.

Table 6. Estimate parameters obtained by fitting Dengue Fever and temperature data in the nonstationary bNB model.

Figure 12. Time series plots of Dengue Fever observations in blue and its estimates values from the model in red when covariate is temperature.

This application is to check how best our developed model can capture the seasonality in Dengue Fever counts data, where Rainfall is also an effect factor (Figure 13). By fitting the Data to the nonstationary bNB model, we obtained the parameters and their corresponding Standard error in Table 7. The fitting plot of Dengue Fever process and its estimates.

Table 7. Estimate parameters obtained by fitting Dengue Fever and rainfall data in the nonstationary bNB model.

Figure 13. Time series plots of Dengue Fever observations in blue and its Estimates values from the model in red when covariate is rainfall.

5. Goodness-of-Fit Test

The chi-square test of goodness of fit is a statistical test used to determine if there is a significant difference between the observed data and the expected data. In the case of time series of count data model, the chi-square test can be used to assess whether the observe counts follow a specified distribution or if there are any significant deviations from the expected counts. To perform a chi-square test for goodness of fit on this study, we first need to define the expected distribution. This is done by fitting bNB models to the observed data in order to determine the parameters of the models, which are displayed in Table 6 and Table 7. Having the expected counts for each time point in the time series via estimate parameters, the observed counts and expected counts are compared using the chi-square test statistic as follows

${\chi }^2={\displaystyle \sum }_{t=0}^n\frac{{\left(y_t-{\widehat{y}}_t\right)}^2}{{\widehat{y}}_t}$ (12)

where the $\sum$ is taken over all time points in the series. The test statistic follows a chi-square distribution with $\left(n-1\right)$ degrees of freedom, where $n$ is the number of time points in the time series. To determine whether the observed counts significantly differ from the expected counts, the test statistic is compared to the critical value of the chi-square distribution at a specified significance level, $\alpha =0.05$ . In this case, the null hypothesis (the observed counts follow the nonstationary bNB model) is rejected if the test statistic is greater than the critical value and it can be conclude that there are significant deviations from the expected counts.

Non-Stationary Branching Negative Binomial Autoregressive Model

Formulation of the null and alternative hypotheses

• Null hypothesis $H_0$ : Dengue Fever data follow the models

• Alternative hypothesis $H_a$ : Dengue Fever data does not follow it

Test Statistic ${\chi }^2$

Figure 14. Chi-square goodness of fit test outcome of the non-stationary bNB model to Dengue Fever count data with covariable, temperature.

Figure 15. Chi-square goodness of fit test outcome of the non-stationary bNB model to Dengue Fever count data with covariable, rainfall.

In this case, two chi-square statistics are determined, including the one where temperature is used as covariate and the one where rainfall is used as covariates, since they generate different expected data for Dengue Fever counts. Then Figure 14 and Figure 15 show the test outcomes from those two cases.

Critical Value

Knowing level of significance $\alpha =5\text{\%}$ and the degree of freedom $df=n-1=372$ , qchisq $\left(1-\alpha ,df\right)$ allowed us to get the critical value, which is ${\chi }^2=417.9735$ , by use of R software.

As result, it’s seeing in Figure 14 that the ${\chi }^2$ value ($\chi =393.83$ ) is less than the critical value (${\chi }_0=417.9735$ ), which implies that the non-stationary bNB model fits well the Dengue Fever with covariate temperature.

Bearing Figure 15, it’s seeing that the ${\chi }^2$ value ($\chi =343.83$ ) is less than the critical value. This also proves that the non-stationary bNB model fits well the Dengue Fever count data when rainfall is the covariate.

6. Conclusions and Suggestions

The findings from both the simulation study and the analysis of the real-world dataset have significant implications for the proposed nonstationary autoregressive model for time series count data with covariates. Starting with the simulation study, it provides a controlled environment to evaluate the performance of the proposed modeling strategy. By simulating count time series data with known characteristics, researchers can assess how well the model captures the underlying seasonality and incorporates the covariate information. The simulation study allows for comparisons between the true values and the estimates obtained through maximum likelihood estimation. The results of the simulation study can provide insights into the accuracy and robustness of the proposed model in capturing the complex dynamics of count data with seasonality and covariate effects. Moving on to the analysis of the real-world dataset, which focuses on an infectious disease with covariates such as temperature and rainfall, the application of the proposed model provides insights into the relationship between the disease occurrence and the covariate variables. The findings can shed light on how weather factors impact the incidence of the infectious disease. By considering nonstationarity and incorporating the covariate information into the model, the analysis can provide a more comprehensive understanding of the factors influencing the disease occurrence and improve the accuracy of predictions or assessments. However, it is important to acknowledge the limitations of the study. One limitation may be the assumptions made in the proposed model. While the nonstationary autoregressive model with covariate incorporation offers an innovative approach, the assumptions underlying the model may not fully capture the complexities of the realworld data. Additionally, the generalizability of the findings might be limited to the specific context of the dataset used in the analysis. Different datasets with varying characteristics and applications may yield different results. Furthermore, the proposed model’s performance depends on the quality and availability of the covariate information. In some cases, the covariates may not be accurately measured or may contain missing data, which can introduce uncertainties and affect the model’s performance. It is crucial to address such limitations and explore techniques to handle missing or imperfect covariate data. Overall, the findings from both the simulation study and the analysis of the realworld dataset provide valuable insights into the proposed nonstationary autoregressive model’s effectiveness in capturing seasonality and incorporating covariates in count time series data. These findings pave the way for further research and refinement of the modeling approach, considering potential limitations and exploring alternative methodologies to enhance the analysis of count data with covariate information and seasonality aspects.

Data Availability Statement

Data available on request from the authors

Conflicts of Interest

The authors declare no potential conflict of interests.

References