Modeling seasonal epidemic data using integer autoregressive model based on binomial thinning

Abstract

The epidemic surveillance data are always in the form of counts observed weekly, monthly or yearly. Integer Autoregressive (INAR) models are the most suitable models for modeling such data. As most of the epidemic data has inherent seasonality in it, the INAR models need to be modified accordingly to take care of such seasonal behavior of the data. In this paper a seasonal geometric INAR(1) model based on binomial thinning is proposed with a seasonal period ‘s’ (GINAR(1) ${}_{s}$ ). The thinning models based on binomial thinning are much easier to work with, than those based on negative binomial thinning, in terms of mathematical and computational complexity. Various inferential and probabilistic properties of the model are studied. The forecasting ability of the GINAR(1) ${}_{s}$ model has been compared with that of the non seasonal counterparts. Extensive simulation study has been carried out to validate the coherent forecasting ability of the model. The model performs well for overdispersed low count time series data. The analysis of an epidemic data has been carried out to examine the performance of the proposed model.

Keywords

Binomial thinning coherent forecasting geometric distribution INAR models seasonality

1. Introduction

There are number of diseases which reoccur seasonally and affect the human life, which in turn brings burden to the economy of a country. Controlling epidemic outbreak and mass vaccination programs to counter these recurrence need lots of money and man power. Most of the epidemic surveillance data are count data and hence researchers use integer-valued autoregressive time series models for modeling such types of data. When there is an inherent seasonality in the series, the usual INtegar AutoRegressive (INAR) models will not give good forecasts. Hence, one has to take into consideration the seasonality in the data while proposing models.

The literature on seasonality of epidemics dates back to Soper (1929), who studied the periodicity in disease prevalence in the case of local epidemics like Measles. He also claimed that infectivity changes seasonally. Seasonality in gastroenteritis and for influenza can be found in the work of Lossli et al. (1943). Even though the seasonal epidemic outbreaks were there since centuries, such types of data were recorded only from later part of 19 ${}^{\rm th}$ century (see Serfling, 1963). A detailed discussion about the seasonal nature of emergence of various diseases can be found in Cook et al. (1990). The applicability of INAR models for modeling epidemic data has been discussed in Cardinal et al. (1999). Survival behavior of host, behavior of pathogen and host immune system are related to the seasonal variation (see Grassly & Fraser, 2006). Various contagious diseases exhibit seasonal patterns, a recent review on seasonality in Influenza and Respiratory Syncytial Virus (RSV) can be found in Bloom-Feshbach et al. (2013).

The smallest time period for which a phenomenon repeats itself is called a seasonal period. The spread of the disease depends on factors such as temperature, humidity, population density, communication facilities etc., in that particular zone or region. The climatic conditions may bring some seasonal pattern in the disease progression over time and hence the seasonality needs to be incorporated in the model.

Integer autoregressive models have been extensively used in modeling various count time series data that occur in health, engineering, insurance etc., see Ristić et al. (2009). Monteiro et al. (2010) have proposed a seasonal INteger AutoRegressive model of order one (INAR(1)). However, in this model, the number of parameters to be estimated is very large and hence not of much practical use. Bourguignon et al. (2016) proposed a seasonal INAR(1) model with Poisson marginal distribution. Recently, Tian et al. (2018) have proposed seasonal INAR(1) model with geometric marginal distribution and based on negative binomial thinning. Awale et al. (2017a) have carried out a comparative study on coherent forecasting ability of non seasonal INAR(1) model with geometric marginal distribution and based on both the binomial and negative binomial thinning, called as Geometric INAR(1) (GINAR(1)) and New Geometric INAR(1) (NGINAR(1)) models respectively. Models based on binomial thinning are much handy than the models based on negative binomial thinning.

In this paper we have proposed a seasonal geometric INAR(1) model with geometric marginal distribution and based on binomial thinning. From the simulation study and data analysis it can be observed that the model with binomial thinning is much easier for coherent forecasting than the model based on negative binomial thinning.

The paper is organized as follows. Section 2 deals with the introduction and basic properties of the model. Estimation of the parameters is discussed in Section 3. Extensive simulation study has been reported in Section 4. An epidemic data has been modeled with the help of the proposed seasonal GINAR(1) ${}_{s}$ model and the parameter estimation and forecasting has been carried out using the method discussed. The details of the analysis are reported in Section 5. Section 6 concludes the paper.

2. Seasonal geometric INAR(1) model based on binomial thinning

The integer-valued auto-regressive process of order one with geometric marginal distribution and seasonal period ‘ $s$ ’, (GINAR(1) ${}_{s}$ ) is defined as,

$\displaystyle X_{t}=\phi\circ X_{t-s}+Z_{t},∼{}∼{}t\geqslant s,$ (1)

where, ‘ $\circ$ ’ is a binomial thinning operator (Stuetal & Van Harn, 1979) , $\phi\circ X=\sum_{i=0}^{X}W_{i}$ , $W_{i}$ are i.i.d. as

$\displaystyle P(W_{i}=0)=1-\phi=1-P(W_{i}=1),∼{}∼{}\phi\in(0,1).$

Here, $Z_{t}=U_{t}∼{}M_{t},\forall∼{}t$ , with $U_{t}$ independent of $M_{t}$ ,

$\displaystyle P[U_{t}=0]=\phi=1-P[U_{t}=1],$

and $\{M_{t}\}$ an i.i.d. geometric sequence with

$\displaystyle P[M_{t}=j]=(1-\theta)\theta^{j},∼{}∼{}j\geqslant 0.$

Then, the marginal distribution of $\{X_{t}\}$ is given by

$\displaystyle P[X_{t}=x]=(1-\theta)\theta^{x},∼{}x\geqslant 0,∼{}0<\theta<1.$

The probability distribution of $Z_{t}$ is

$\displaystyle P(Z_{t}=z)=\left\{\begin{array}[]{ll}(1-\theta+\phi\theta),&% \text{if}\ z=0,\\ (1-\theta)(1-\phi)\theta^{z},&\text{if}\ z\geqslant 1.\end{array}\right.$ (2)

.

Following Bourguignon et al. (2016) and Tian et al. (2018) it is easy to show that the process $\{X_{t}\}_{t\geqslant 1}$ is a stationary process.

The sample path of the process along with autocorrelation function (ACF) and partial autocorrelation function (PACF) are shown in Fig. 1. The process in Eq. (1) has two components, viz., autoregressive and innovation. The number of survivors from the elements at time $t-s$ are given by $\phi\circ X_{t-s}$ and the new arrivals during the time interval $(t-s,t]$ are denoted by the innovation term $Z_{t}$ .

Figure 1.

Time series, sample ACF and PACF plots for various combination of parameters for GINAR(1) ${}_{s}$ processes, first row: $\phi=$ 0.3, $\theta=$ 0.5; second row: $\phi=$ 0.5, $\theta=$ 0.8; third row: $\phi=$ 0.6, $\theta=$ 0.4; fourth row: $\phi=$ 0.8, $\theta=$ 0.3, $s=$ 52.

The marginal mean, variance and pgf of $\{X_{t}\}$ are $E(X_{t})=\theta/(1-\theta)$ , $V(X_{t})=\theta/(1-\theta)^{2}$ and $\Phi_{X_{t}}(\upsilon)=(1-\theta)/(1-\upsilon\theta)$ , $|\upsilon|<\frac{1}{\theta}$ , respectively. The conditional mean and the variance are

$\displaystyle E(X_{t}|X_{t-s})=\phi X_{t-s}+(1-\phi)\frac{\theta}{1-\theta}$

and

$\displaystyle V(X_{t}|X_{t-s})=\phi(1-\phi)X_{t-s}+(1-\phi)\theta(1+\phi\theta% )\frac{1}{(1-\theta)^{2}}.$

The conditional pmf is derived in the following theorem. The proof of which is deferred to Appendix.

.

Let $\{X_{t}\}$ be a GINAR(1) ${}_{s}$ process defined as in Eq. (1). Then,

$\displaystyle P(X_{t}=y|X_{t-s}=x)=\left\{\begin{array}[]{l}(1-\theta)\theta^{% y-x}\sum\limits_{r=0}^{y}{{x}\choose{r}}\theta^{x-r}\phi^{r}(1-\phi)^{x-r+1}\\ +{{x}\choose{y}}\phi^{y+1}(1-\phi)^{x-y},∼{}∼{}∼{}∼{}∼{}∼{}∼{}∼{}∼{}∼{}∼{}∼{}∼% {}∼{}∼{}\mbox{if}∼{}∼{}y\leqslant x,\\ \\ (1-\theta)\theta^{y-x}(1-\phi)\{\phi+(1-\phi)\theta\}^{x},∼{}∼{}∼{}∼{}\mbox{if% }∼{}∼{}y>x.\end{array}\right.$

From Eq. (1) and using the recursion, we can write,

$\displaystyle X_{t+k}=\phi^{q}\circ X_{t-r}+\sum_{j=0}^{q-1}\phi^{j}\circ Z_{t% +k-js},$ (3)

where, $q=[\frac{k}{s}]$ , $r=qs-k$ and $[\cdot]$ denotes the upper integer part of the ‘ $\cdot$ ’.

The $k$ -step ahead conditional pgf of the process can be derived as,

$\displaystyle G_{X_{t+k}|X_{t-r}}(\upsilon)=E({\upsilon}^{X_{t+k}}|X_{t-r})=E(% {\upsilon}^{\phi^{q}\circ X_{t-r}}|X_{t-r})\prod_{j=0}^{q-1}E({\upsilon}^{\phi% ^{j}\circ Z_{t+k-js}}).$

Using the property of binomial thinning operator (Al-osh & Alzaid, 1987), one can write

$\displaystyle\phi^{q}\circ X_{t-r}|X_{t-r}\sim Bin(X_{t-r},\phi^{q})$ (4)

and therefore,

$\displaystyle E({\upsilon}^{\phi^{q}\circ X_{t-r}}|X_{t-r})={(1-\phi^{q}+% \upsilon\phi^{q})}^{X_{t-r}}.$ (5)

As $Z_{t}$ ’s are i.i.d, we can write,

$\displaystyle\prod_{j=0}^{q-1}E({\upsilon}^{\phi^{j}\circ Z_{t}})=\prod_{j=0}^% {q-1}E(E({\upsilon}^{\phi^{j}\circ Z_{t}}|Z_{t}))=\prod_{j=0}^{q-1}E(1-\phi^{j% }+\upsilon\phi^{j})^{Z_{t}}=\prod_{j=0}^{q-1}G_{Z_{t}}(1-\phi^{j}+\upsilon\phi% ^{j}).$

But,

$\displaystyle G_{Z_{t}}(\upsilon)=\frac{1-\theta+\phi\theta(1-\upsilon)}{1-% \theta\upsilon},∼{}∼{}\theta\upsilon\neq 1.\ \text{Therefore,}$

$\displaystyle\prod_{j=0}^{q-1}E({\upsilon}^{\phi^{j}\circ Z_{t}})=\prod_{j=0}^% {q-1}\bigg{(}\frac{1-\theta+\phi\theta(1-(1-\phi^{j}+\upsilon\phi^{j}))}{1-% \theta(1-\phi^{j}+\upsilon\phi^{j})}\bigg{)}=\prod_{j=0}^{q-1}\bigg{(}\frac{1-% \theta+\theta\phi^{j+1}(1-\upsilon)}{1-\theta+\theta\phi^{j}(1-\upsilon)}\bigg% {)}=\frac{1-\theta+\theta\phi^{q}(1-\upsilon)}{1-\theta\upsilon},$

which implies that

$\displaystyle G_{\sum\limits_{j=0}^{q-1}\phi^{j}\circ Z_{t+k-js}}(\upsilon)=% \frac{1-\theta(1-\phi^{q}+\phi^{q}\upsilon)}{1-\theta\upsilon}.$ (6)

Hence, the conditional pgf is given by,

$\displaystyle G_{X_{t+k}|X_{t-r}}(\upsilon)={\big{(}1-\phi^{q}+\upsilon\phi^{q% }\big{)}}^{X_{t-r}}\bigg{(}\frac{1-\theta+\theta\phi^{q}(1-\upsilon)}{1-\theta% \upsilon}\bigg{)},q=\bigg{[}\frac{k}{s}\bigg{]}.$ (7)

From Eq. (7) it can be shown that,

$\displaystyle G_{X_{t+k}|X_{t-r}}(\upsilon)\rightarrow G_{X_{t}}(\upsilon)=% \frac{1-\theta}{1-\theta\upsilon}\ \text{as}\ k\rightarrow\infty.$

Using Eq. (6 ), the pmf of the second term in Eq. (3) is given as

$\displaystyle P\bigg{(}\sum\limits_{j=0}^{q-1}\phi^{j}\circ Z_{t+k-js}=y-r% \bigg{)}=\left\{\begin{array}[]{ll}(1-\theta+\theta\phi^{q}),&\text{if}\ r=y,% \\ (1-\theta)(1-\phi^{q})\theta^{y-r},&\text{if}\ r=0,1,\ldots,y-1.\end{array}\right.$ (8)

From Eqs (3), (4) and (8), the $k$ -step ahead conditional pmf of the process can be derived as

$\displaystyle P(X_{t+k}=y|X_{t-r}=x)=\left\{\begin{array}[]{ll}(1-\theta-{\phi% }^{q}+\theta{\phi}^{q})\theta^{y-x}\sum\limits_{r=0}^{y}{{x}\choose{r}}\theta^% {x-r}\phi^{qr}(1-\phi^{q})^{x-r}\\ \quad+{{x}\choose{y}}\phi^{q(y+1)}(1-\phi^{q})^{x-y},&\text{if}∼{}∼{}y% \leqslant x,\\ \\ (1-\theta-{\phi}^{q}+\theta{\phi}^{q})\theta^{y-x}({\phi}^{q}+(1-\phi^{q})% \theta)^{x},&\text{if}∼{}∼{}y>x,\end{array}\right.$ (9)

where, $q=\big{[}\frac{k}{s}\big{]}$ . The $k$ -step ahead conditional mean and variance are respectively,

$\displaystyle E(X_{t+k}|X_{t-r})=\phi^{q}X_{t-r}+(1-\phi^{q})\frac{\theta}{1-\theta}$ (10)

and

$\displaystyle V(X_{t+k}|X_{t-r})=\phi^{q}(1-\phi^{q})X_{t-r}+\Big{(}\frac{\phi% }{1-\phi^{2}}(1-\phi^{q-1}-\phi^{q}+\phi^{2q-1})\Big{)}\mu_{z}+\frac{1-\phi^{2% q}}{1-\phi^{2}}\sigma_{z}^{2},$ (11)

where,

$\displaystyle\mu_{z}=\frac{\theta(1-\phi)}{1-\theta}\ \text{and}\ \sigma_{z}^{% 2}=\frac{(1-\phi)\theta(1+\phi\theta)}{(1-\theta)^{2}}$

are the mean and variance of $Z_{t}$ . From the Eqs (10) and (11), we observe that, as $k\rightarrow\infty$ the conditional mean and variance converge to the marginal mean and variance respectively. In the following lemma, we obtain the autocorrelation function of the process $\{X_{t}\}$ defined in Eq. (1). The proof of the lemma is deferred to the Appendix.

.

If $\{X_{t}\}$ is a process defined as in Eq. (1), then the autocorrelation function of the process $\{X_{t}\}$ is given by

$\displaystyle\rho(h)=\left\{\begin{array}[]{ll}\rho^{h/s},&\text{if}\ h∼{}% \text{is a multiple of}∼{}s,\\ 0,&\text{otherwise}.\end{array}\right.$

3. Estimation of the parameters of GINAR(1)

{}_{s}

model

In this section we consider the conditional maximum likelihood and conditional least squares estimation of the model parameters. Conditional maximum likelihood estimators (CML) can be obtained by maximizing the conditional log-likelihood function

$\displaystyle\log L(x_{1},\ldots,x_{n};\phi,\theta)=\sum_{t=s+1}^{n}\log P(X_{% t}|X_{t-s}),$

where, $P(X_{t}|X_{t-s})$ is given in Theorem 1. The conditional least squares (CLS) estimators of the parameters can be obtained by minimizing the function

$\displaystyle S_{n}(\phi,\theta)=\sum_{t=s+1}^{n}(X_{t}-E(X_{t}|X_{t-s}))^{2},$

with respect to $\phi$ and $\theta$ . The differentiation results in the following estimating equations,

$\displaystyle\sum\limits_{t=s+1}^{n}\bigg{(}X_{t}-\phi X_{t-s}-(1-\phi)\frac{% \theta}{1-\theta}\bigg{)}X_{t-s}=0$ (12)

and

$\displaystyle\sum\limits_{t=s+1}^{n}\bigg{(}X_{t}-\phi X_{t-s}-(1-\phi)\frac{% \theta}{1-\theta}\bigg{)}=0.$ (13)

Solving Eqs (12) and (13) for $\phi$ and $\theta$ , we get the following estimators,

$\displaystyle\hat{\phi}_{\textit{CLS}}=\frac{(n-s)\sum_{t=s+1}^{n}X_{t}X_{t-s}% -\big{(}\sum_{t=s+1}^{n}X_{t}\big{)}\big{(}\sum_{t=s+1}^{n}X_{t-s}\big{)}}{(n-% s)\sum_{t=s+1}^{n}{X}_{t-s}^{2}-\big{(}\sum_{t=s+1}^{n}X_{t-s}\big{)}^{2}}$

and

$\displaystyle\hat{\theta}_{\textit{CLS}}=\frac{\sum_{t=s+1}^{n}X_{t}-\hat{\phi% }_{\textit{CLS}}\sum_{t=s+1}^{n}X_{t-s}}{(n-s)(1-\hat{\phi}_{\textit{CLS}})+(% \sum_{t=s+1}^{n}X_{t}-\hat{\phi}_{\textit{CLS}}\sum_{t=s+1}^{n}X_{t-s})}.$

The asymptotic normality of these estimators is proved in the following theorem.

.

Let $\hat{\Theta}_{\textit{CLS}}=(\hat{\phi}_{\textit{CLS}},\hat{\theta}_{\textit{% CLS}})^{\prime}$ be the conditional least squares (CLS) estimator of the parameter $\Theta=(\phi,\theta)^{\prime}$ of a seasonal INAR(1) model with geometric marginal distribution as defined in Eq. (1). Then,

$\displaystyle\sqrt{n}(\hat{\Theta}_{\textit{CLS}}-\Theta)\lx@stackrel{{% \scriptstyle d}}{{\rightarrow}}N_{2}(\underline{0},\Sigma_{\phi,\theta}),$

where,

$\displaystyle\Sigma_{\phi,\theta}=S^{-1}\textit{VS}^{-1^{\prime}},$

$\displaystyle S=\left(\begin{array}[]{cc}\displaystyle\frac{-\theta}{(1-\theta% )^{2}}&\frac{-\theta(1-\phi)}{(1-\theta)^{3}}\\ 0&-\frac{1-\phi}{(1-\theta)^{2}}\end{array}\right)$

and

$\displaystyle V=\left(\begin{array}[]{cc}V_{11}&V_{12}\\ V_{21}&V_{22}\end{array}\right).$

The elements of the matrix $V$ are

$\displaystyle V_{11}=\frac{(1-\phi)\theta\left(1+\theta-\phi\left(-1-4\theta+2% \theta^{2}+\theta^{3}\right)\right)}{(1-\theta)^{4}},$ $\displaystyle V_{12}=V_{21}=\frac{(1-\phi)\theta\left(\theta+\phi\left(1-% \theta+\theta^{3}\right)\right)}{(1-\theta)^{4}},$

and

$\displaystyle V_{22}=\frac{\left(1-\phi^{2}\right)\theta}{(1-\theta)^{2}}.$

3.1 Coherent forecasting and forecast accuracy measures

Even though various extensions of the INAR models have been developed, very little effort has been made to address the forecasting issues. The traditional minimum mean square error (MMSE) forecast, the conditional mean, $\hat{X}_{t+k}=E(X_{t+k}|\mathcal{F}_{t})$ , where, $\mathcal{F}_{t}$ is the available information up to time $t$ , need not result as a count in the case of INAR models. Therefore, to address this problem, coherent forecasts are introduced. Coherent forecasts are the forecasts which take only integer values, i.e., forecasts that are in the form of the underlying data. The coherent forecasts are usually obtained from the $k$ -step ahead conditional distribution. The idea of coherent forecasting for INAR models was first proposed by McCabe (2004). These authors have considered the coherent forecasting issues related to Poisson integer autoregressive model of order one. Jung and Tremayne (2006) extended the idea of coherent forecasting to INAR(2) Poisson models. Maiti et al. (2015) considered the coherent forecasting in count time series, using Box and Jenkins’s AR( $p$ ) model. Maiti and Biswas (2015) studied coherent forecasting in geometric INAR(1) model based on binomial thinning.

The forecast accuracy is measured using measures such as Predictive Root Mean Square Error (PRMSE), Prediction Mean Absolute Error (PMAE) and Percentage of True Prediction (PTP). Let the observed data set $\{X_{1},\ldots,X_{n},X_{n+1},\ldots,X_{n+m}\}$ be partitioned into two sets. The first $n$ observations are used for the estimation of the parameters of the model and remaining $m$ observations called the test set is used for the computation. The measures are given below,

$\displaystyle\textit{PRMSE}(k)=\sqrt{\frac{1}{m-k+1}\sum\limits_{i=1}^{n+m-k}% \Big{(}X_{i+k}-\hat{X}_{i+k}^{\textit{Mean}}\Big{)}^{2}},$ $\displaystyle\textit{PMAE(k)}=\frac{1}{m-k+1}\sum\limits_{i=1}^{n+m-k}|X_{i+k}% -\hat{X}_{i+k}^{\textit{Med}}|$

and

$\displaystyle\textit{PTP(k)}=\frac{1}{m-k+1}\sum\limits_{i=1}^{n+m-k}I\Big{(}X% _{i+k}=\hat{X}_{i+k}\Big{)}\times 100\%,$

where, $I(\cdot)$ is the indicator function. For computation of PTP, we have used $\hat{X}_{n+k}=\hat{X}_{n+k}^{\textit{Med}},\hat{X}_{n+k}^{\textit{Mode}}$ and the $k$ -step ahead conditional mean, after rounding it to the nearest integer. The quantities $\hat{X}_{n+k}^{\textit{Mode}}$ , $\hat{X}_{n+k}^{\textit{Median}}$ and $\hat{X}_{n+k}^{\textit{Mean}}$ are obtained from the $k$ -step ahead conditional pmf.

The standard prediction intervals assume the predictive probability distribution to be symmetric. However, for these models, it is uni-modal and positively skewed. Hence, we find the $100(1-\gamma)\%$ highest predictive probability (HPP) interval for $X_{t+k}$ as, $C_{k}=(X_{L},X_{U}),$ with $C_{k}=\{y:p_{k}(y|x)\geqslant K_{\gamma}\}$ and $K_{\gamma}$ is the largest number such that,

$\displaystyle P(X_{L}\leqslant X_{t+k}\leqslant X_{U}|X_{t}=x)=\sum\limits_{y=% X_{L}}^{X_{U}}p_{k}(y|x)\geqslant(1-\gamma),$

where,

$\displaystyle p_{k}(y|x)=P(X_{t+k}=y|X_{t}=x).$

4. Simulation study

The simulation results in Table 1 present the parameter estimates, their standard errors and coverage probability associated with seasonal GINAR(1) (GINAR(1) ${}_{s}$ ) model. For this study, we have simulated 1000 series each of size 300, 500 and 1000 with parameter values as given in Table 1, with seasonal period $s=$ 52. The average values of the estimates over 1000 simulations, estimate of standard errors and the coverage probability for 95% confidence interval are given in row one, row two and row three respectively for each sample size and for all parameter combinations. From this table, it can be seen that both CML and CLS estimators perform well and get stabilized for larger sample sizes. CML estimates have smaller standard error (SE) than the CLS estimators. Therefore, we have used CML estimates for the data analysis.

Table 1
Parameter estimates, their standard errors and coverage probability associated with GINAR(1) ${}_{s}$ model

$n$	$\hat{\phi}_{cml}$	$\hat{\theta}_{cml}$	$\hat{\phi}_{cls}$	$\hat{\theta}_{cls}$
	$\phi=$ 0.2, $\theta=$ 0.3
300	0.1964	0.2997	0.1904	0.2985
	0.0032	0.0016	0.0040	0.0015
	0.949	0.946	0.945	0.943
500	0.1999	0.3000	0.1962	0.2992
	0.0019	0.0009	0.0023	0.0009
	0.947	0.953	0.955	0.945
1000	0.1988	0.2999	0.1978	0.2991
	0.0009	0.0004	0.0023	0.0009
	0.945	0.947	0.943	0.947
	$\phi=$ 0.5, $\theta=$ 0.6
300	0.4971	0.5979	0.4846	0.5983
	0.0018	0.0017	0.0032	0.0018
	0.945	0.953	0.945	0.954
500	0.4977	0.5982	0.4940	0.5989
	0.0011	0.0010	0.0020	0.0010
	0.945	0.953	0.953	0.950
1000	0.5001	0.5990	0.4983	0.5994
	0.0005	0.0005	0.0010	0.0005
	0.953	0.950	0.954	0.943
	$\phi=$ 0.8, $\theta=$ 0.9
300	0.7996	0.8971	0.7898	0.8978
	0.0005	0.0008	0.0021	0.0010
	0.944	0.946	0.0946	0.948
500	0.7796	0.8983	0.7973	0.4940
	0.0003	0.0004	0.0012	0.0030
	0.947	0.948	0.943	0.949
1000	0.7997	0.8991	0.7958	0.8989
	0.0001	0.0002	0.0005	0.0002
	0.947	0.945	0.944	0.954
	$\phi=$ 0.9, $\theta=$ 0.5
300	0.8983	0.4868	0.8878	0.4896
	0.0009	0.0048	0.0018	0.0073
	0.942	0.944	0.925	0.984
500	0.8987	0.4919	0.8922	0.4940
	0.0005	0.0027	0.0011	0.0030
	0.942	0.955	0.933	0.949
1000	0.8985	0.4937	0.8970	0.4982
	0.0002	0.0014	0.0005	0.0015
	0.947	0.943	0.942	0.958

First row – estimate, second row – standard error, third row – coverage probability for respective sample sizes.

A simulation study has been carried out to examine the performance of the seasonal as well as the non seasonal models for seasonal data. We have carried out 1000 simulations of the series with size 600 from GINAR(1) ${}_{s}$ model for various parameter combinations mentioned in the Table 2. Out of these 600 observations first 400 were used for the parameter estimation and remaining for computation of forecast accuracy measures.

Table 2

PRMSE and PMAE comparison of GINAR(1) and GINAR(1) ${}_{s}$ models with HPP coverage

	GINAR(1)		GINAR(1) ${}_{s}$			GINAR(1)		GINAR(1) ${}_{s}$
	PRMSE	PMAE	PRMSE	PMAE	HPP Cov.	PRMSE	PMAE	PRMSE	PMAE	HPP Cov.
$k$	$\phi=$ 0.2, $\theta=$ 0.3					$\phi=$ 0.4, $\theta=$ 0.5
1	0.7776	0.4282	0.7623	0.4199	0.973	1.4047	1.0001	1.2882	0.7667	0.949
2	0.7776	0.4281	0.7626	0.4200	0.973	1.4036	0.9989	1.2876	0.7666	0.949
3	0.7774	0.4281	0.7624	0.4199	0.973	1.4033	0.9985	1.2873	0.7663	0.949
4	0.7776	0.4282	0.7624	0.4201	0.973	1.4033	0.9985	1.2873	0.7662	0.949
5	0.7774	0.4280	0.7624	0.4200	0.973	1.4026	0.9982	1.2865	0.7657	0.949
	$\phi=$ 0.6, $\theta=$ 0.7					$\phi=$ 0.8, $\theta=$ 0.9
1	2.7882	1.9505	2.2244	1.2330	0.855	9.3811	6.6127	5.6197	2.4265	0.607
2	2.7871	1.9488	2.2244	1.2328	0.855	9.3782	6.6074	5.6204	2.4273	0.607
3	2.7867	1.9489	2.2235	1.2323	0.855	9.3785	6.6081	5.6213	2.4282	0.607
4	2.7864	1.9491	2.2235	1.2321	0.855	9.3769	6.6066	5.6184	2.4269	0.607
5	2.7861	1.9492	2.2225	1.2321	0.855	9.3789	6.6082	5.6213	2.4277	0.607
	$\phi=$ 0.8, $\theta=$ 0.3					$\phi=$ 0.9, $\theta=$ 0.2
1	0.7758	0.4330	0.4650	0.1618	0.870	0.5387	0.2527	0.2347	0.0493	0.811
2	0.7750	0.4315	0.4649	0.1618	0.870	0.5380	0.2511	0.2348	0.0494	0.811
3	0.7750	0.4315	0.4649	0.1618	0.870	0.5381	0.2512	0.2348	0.0494	0.811
4	0.7751	0.4315	0.4648	0.1618	0.870	0.5378	0.2510	0.2347	0.0494	0.811
5	0.7751	0.4317	0.4649	0.1619	0.870	0.5376	0.2510	0.2346	0.0494	0.811

Table 2 presents the comparison of seasonal GINAR(1) model and non seasonal GINAR(1) model. The estimate of coverage probability for 95% HPP interval is also given in the same table for GINAR(1) ${}_{s}$ model. Here it can be observed that the range of $\phi$ and $\theta$ is $(0,1)$ . When $\phi$ takes value zero, the process $\{X_{t}\}$ becomes independently and identically distributed (i.i.d.) and when it takes value one, the model becomes non stationary. Most of the INAR(1) models do not perform well when the parameter values are near to the boundary of the non-stationary region of parameter space. Similar phenomenon is observed in the case of the test for randomness in $\phi$ (see Awale et al., 2019) and the Fisher’s dispersion test (see Lee et al., 2017), i.e., the tests do not maintain the level and the power, when the thinning parameter values are near to one. In the current model, seasonal component is added to the basic GINAR(1) model and hence when the parameter values near one, the HPP intervals based on seasonal forecast distribution do not perform well.

It can be observed that the PRMSE and PMAE for the seasonal model are smaller than that of the non-seasonal model and this is observed for all the parameter combinations. Here, we conclude that the GINAR(1) ${}_{s}$ model performs better than the GINAR(1) model in terms of PRMSE and PMAE. We have obtained the coherent forecasts such as median and mode, using the $k$ -step ahead conditional distribution given in Eq. (9). The traditional mean forecast is obtained using the Eq. (10).

Table 3

PTP comparison of GINAR and GINAR(1) ${}_{s}$ models

	GINAR(1)			GINAR(1) ${}_{s}$			GINAR(1)			GINAR(1) ${}_{s}$
	Med	Mode	Mean	Med	Mode	Mean	Med	Mode	Mean	Med	Mode	Mean
$k$	$\phi=$ 0.2, $\theta=$ 0.3						$\phi=$ 0.4, $\theta=$ 0.5
1	70.02	70.04	65.87	69.52	70.08	63.65	38.37	49.77	25.12	52.56	53.93	26.66
2	70.04	70.04	66.31	69.51	70.07	63.66	38.34	49.79	25.13	52.56	53.92	26.67
3	70.03	70.03	66.29	69.51	70.07	63.66	38.39	49.80	25.14	52.57	53.93	26.67
4	70.03	70.03	66.29	69.50	70.06	63.65	38.37	49.79	25.15	52.57	53.93	26.68
5	70.03	70.03	66.30	69.51	70.07	63.66	38.38	49.80	25.15	52.59	53.94	26.68
	$\phi=$ 0.6, $\theta=$ 0.7						$\phi=$ 0.8, $\theta=$ 0.9
1	18.52	29.74	13.47	45.33	46.06	19.45	5.42	9.99	3.98	38.36	38.36	7.75
2	18.51	29.89	13.48	45.34	46.07	19.46	5.42	10.02	3.98	38.36	38.35	7.74
3	18.52	29.89	13.48	45.35	46.07	19.46	5.42	10.04	3.98	38.36	38.35	7.75
4	18.51	29.89	13.47	45.35	46.07	19.46	5.42	10.04	3.99	38.36	38.36	7.75
5	18.50	29.89	13.47	45.34	46.06	19.46	5.41	10.04	3.98	38.37	38.36	7.76
	$\phi=$ 0.8, $\theta=$ 0.3						$\phi=$ 0.9, $\theta=$ 0.2
1	69.74	69.90	58.47	86.78	86.83	86.66	79.78	79.90	78.98	95.61	95.62	95.59
2	69.94	69.94	58.25	86.78	86.83	86.66	79.95	79.95	79.42	95.61	95.61	95.59
3	69.94	69.94	58.39	86.78	86.83	86.66	79.94	79.94	79.33	95.61	95.61	95.59
4	69.95	69.95	58.39	86.78	86.83	86.66	79.95	79.95	79.34	95.61	95.61	95.58
5	69.94	69.94	58.38	86.78	86.82	86.65	79.95	79.95	79.34	95.61	95.61	95.58

The PTP values have been computed for the six parameter combinations covering the range of $\phi$ and $\theta$ for both the models (see Table 3). The PTP for mean is computed using the rounded part of the $k$ -step ahead conditional mean. Here also, it can be observed that the PTP for mean, median and mode are higher in the seasonal model than those in the non seasonal model. This implies that the proportion of forecasting exact values is higher in the seasonal model than in the non-seasonal model. It can also be observed that the seasonal model also performs much better in terms of forecast accuracy measures (see Tables 2 and 3).

Table 4

Model selection using AIC and BIC

Model	Parameter estimates	AIC	BIC
GINAR(1)	$\hat{\phi}=$ 0.1197 (0.0014), $\hat{\theta}=$ 0.5508 (0.0011)	969.18	970.17
GINAR(1) ${}_{S}$	$\hat{\phi}=$ 0.0050 (0.0015), $\hat{\theta}=$ 0.5595(0.0013)	814.74	815.73
NGINAR(1)	$\hat{\alpha}=$ 0.2432(0.0030), $\hat{\mu}=$ 1.2211 (0.0060)	964.52	965.51
NGINAR(1) ${}_{S}$	$\hat{\alpha}=$ 0.0058 (0.0017), $\hat{\mu}=$ 1.2703 (0.0063)	814.74	815.73
POINAR(1)	$\hat{\phi}=$ 0.0535 (0.0028), $\hat{\lambda}=$ 1.2037 (0.0048)	981.42	988.90
POINAR(1) ${}_{S}$	$\hat{\phi}=$ 0.1679 (0.0035), $\hat{\lambda}=$ 1.0616 (0.0054)	829.19	836.67

Figure 2.

Time series, ACF and PACF plots for Campylobacteriosis data.

5. Application to Campylobacter infection data

In this section, we consider the problem of modeling a seasonal epidemic data using the seasonal geometric INAR(1) model. The data are from Rottweil County of Baden-Württemberg state in Germany, and consist of 312 weekly observations pertaining to the year 2001 to 2006. The data are available at the website of the Robert Koch Institute, Germany (https://survstat.rki.de). From the the ACF plot in Fig. 2, it can be seen that the series has seasonality and from PACF plot the AR(1) structure is evident. The seasonal behaviour of this data was also discussed in Louis et al. (2005) and Meldrun et al. (2005). The mean and variance are 1.2756 and 2.0781 respectively, implying overdispersion. From mean, variance and Fig. 3 it can be confirmed that the seasonal INAR(1) with geometric marginal distribution would be a suitable model for the data.

Table 5
Point forecast and HPP intervals for Compylobacteriosis data $\gamma=$ 0.95

$k$	$X_{t+k}$	Mean	Mode	Median	HPP	Mean	Mode	Median	HPP
		NGINAR(1) ${}_{S}$				GINAR(1) ${}_{S}$
1	0	1.26	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$
2	0	1.47	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$
3	2	1.20	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$
4	2	1.20	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$
5	1	1.15	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$
6	1	1.20	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$
7	3	1.26	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$
8	7	1.20	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$
9	4	1.15	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$
10	0	1.26	0	1	$\{0,..,5\}$	1.25	0	1	$\{0,..,5\}$

Figure 3.

Frequency plot for Campylobacteriosis data.

The suitability of geometric model can be seen from the bar-plot and AIC, BIC computed for the similar models in Table 4. Parameter estimates and their standard errors are given in the same table. It appears that seasonal NGINAR(1) ${}_{s}$ and GINAR(1) ${}_{s}$ provide a better fit for the data as compared to the other models.

From Table 5, it can be observed that the median and mode forecasts are similar for both the seasonal as well as non seasonal model. The mean forecast varies with lag in the case of seasonal model, but no variation found in the case of non seasonal model. Similarly, it can be seen that the HPP intervals remain the same for all the steps in non-seasonal model, but are varying in the case of seasonal model. It may be noted that as the periodicity in this model is $s=$ 52, we may not see significant changes in the forecasts up-to 10 steps. The changes in the forecasts obtained using seasonal model may be more prominent for the series with smaller periods.

From the Table 5, it can be observed that the GINAR(1) ${}_{s}$ model forecasts are as good as that of NGINAR(1) ${}_{s}$ model.

6. Conclusions

In this paper, a seasonal GINAR(1) (GINAR(1) ${}_{s}$ ) model for overdispersed data based on binomial thinning is proposed. CLS and CML estimates of the parameters are studied, CML estimators were found to be better than the CLS. The coherent forecast methodology has been used for forecasting. The forecasts are validated by extensive simulation study. The coherent forecasts obtained using GINAR(1) ${}_{s}$ model for overdispersed data performed equally well with that of the model based on negative binomial thinning and with low cost of computational time and mathematical complexity.

Footnotes

Acknowledgments

Authors are thankful to anonymous referee and editor for their comments and suggestions. Manik Awale would like to acknowledge the ASPIRE research grant of Savitribai Phule Pune University. The research of Akanksha Kashikar was supported in part by a grant from the Board of College and University Development, Savitribai Phule Pune University. T. V. Ramanathan’s research was partially supported by a grant from the Department of Science and Technology (DST), Government of India, SR/S4/MS-866/13.

Appendix

References

Al-osh

M. A.

, & Alzaid

A. A.

(1987). First order integer-valued autoregressive (INAR(1)) processes. Journal of Time Series Analysis, 8, 261-275.

Awale

Ramanathan

T. V.

, & Kale

(2017a). Coherent forecasting in integer-valued AR(1) models with geometric marginals. Journal of Data Science, 15(1), 95-114.

Awale

Balakrisha

, & Ramanathan

T. V.

(2019). Testing the constancy of the thinning parameter in a random coefficient integer autoregressive model. Statistical Papers, 60(5), 1515-1539.

Bloom-Feshbach

Alonso

W. J.

Charu

Tamerius

Simonsen

Miller

M. A.

, & Viboud

(2013). Latitudinal variations in seasonal activity of influenza and respiratory syncytial virus (RSV): A global comparative review. PLOS ONE, 8(2), e54445.

Bourguignon

Vasconcellos

K. L. P.

Reisen

V. A.

, & Ispány

(2016). A Poisson INAR(1) process with a seasonal structure. Journal of Statistical Computation and Simulation, 86(2), 373-387.

Cardinal

Roy

, & Lambert

(1999). On the application of integer-valued time series models for the analysis of diseases incidence. Statistics in Medicine, 18, 2025-2039.

Cook

S. M.

Glass

R. I.

LeBaron

C. W.

, & Ho

M.-S.

(1990). Global seasonality of rotavirus infections. Bulletin of the World Health Organization, 68(2), 171-177.

Freeland

R. K.

, & McCabe

B. P. M.

(2004). Forecasting discrete low count time series. International Journal of Forecasting, 20, 427-434.

Freeland

R. K.

, & McCabe

(2005). Asymptotic properties of CLS estimators in the Poisson AR(1) model. Statistics and Probability Letters, 73, 147-153.

10.

Godambe

V. P.

, & Heyde

C. C.

(1987). Quasi-likelihood and optimal estimation. International Statistical Review, 55, 231-244.

11.

Grassly

N. C.

, & Fraser

(2006). Seasonal infectious disease epidemiology. Proceedings of The Royal Society-B, 273, 2541-2550.

12.

Jung

R. C.

, & Tremayne

A. R.

(2006). Coherent forecasting in integer time series models. International Journal of Forecasting, 22(2), 223-238.

13.

Klimko

L. A.

, & Nelson

P. I.

(1978). On conditional least squares estimation for stochastic processes. The Annals of Statistics, 6(3), 629-642.

14.

Lee

Park

, & Cathy

W. S. Chen

(2017). On Fisher’s dispersion test for integer-valued autoregressive Poisson models with applications, Communications in Statistics – Theory and Methods, 46(20), 9985-9994.

15.

Loosli

C. G.

Lemon

H. M.

Robertson

O. H.

, & Appel

(1943). Experimental airborne influenza infection. I. influence of humidity on survival of virus in air. Proceedings of the Society for Experimental Biology and Medicine, 53(2), 205-206.

16.

Louis

V. R.

Gillespie

I. A.

O’Brien

S. J.

Russek-Cohen

Pearson

A. D.

, & Colwell

R. R.

(2005). Temperature-driven Campylobacter seasonality in England and Wales. Applied and Environmental Microbiology, 71(1), 85-92.

17.

Maiti

, & Biswas

(2015). Coherent forecasting for over-dispersed time series of counts data. Brazilian Journal of Probability and Statistics, 29(4), 747-766.

18.

Maiti

Biswas

, & Das

(2015). Coherent forecasting for count time series using Box-Jenkins’s A⁢R⁢(p) model. Statistica Neerlandica, 70(2), 123-145.

19.

McKenzie

(1986). Autoregressive moving-average processes with negative-binomial and geometric marginal distributions. Advances in Applied Probability, 18, 679-705.

20.

Meldrum

R. J.

Griffiths

J. K.

Smith

R. M. M.

, & Evans

M. R.

(2005). The seasonality of human campylobacter infection and campylobacter isolates from fresh, retail chicken in Wales. Epidemiology and Infection, 133, 49-52.

21.

Monteiro

Scotto

M. G.

, & Pereira

(2010). Integer-valued autoregressive processes with periodic structure. Journal of Statistical Planning and Infer ence, 140(6), 1529-1541.

22.

Ristić

M. M.

Bakouch

H. S.

, & Nastić

A. S.

(2009). A new geometric first or der integer-valued autoregressive (NGINAR(1)) process. Journal of Statistical Planning and Inference, 139, 2218-2226.

23.

Serfling

R. E.

(1963). Methods for current statistical analysis of excess pneumonia-influenza deaths. Public Health Reports, 78(6), 494-506.

24.

Soper

H. E.

(1929). The interpretation of periodicity in disease prevalence. Journal of the Royal Statistical Society, 92(1), 34-73.

25.

Stuetal

F. W.

, & Van Harn

(1979). Discrete analogues of self-decomposability and stability. Annals of Probability, 7, 893-899.

26.

Tian

Wang

, & Cui

(2018). A seasonal geometric INAR process based on negative binomial thinning operator. Statistical Papers, doi: 10.1007/s00362-018-1060-7.

Modeling seasonal epidemic data using integer autoregressive model based on binomial thinning

Abstract

Keywords

1. Introduction

2. Seasonal geometric INAR(1) model based on binomial thinning

.

.

.

.

4. Simulation study

Table 1 Parameter estimates, their standard errors and coverage probability associated with GINAR(1) s model

Table 5 Point forecast and HPP intervals for Compylobacteriosis data γ = 0.95

Footnotes

Acknowledgments

Appendix

References

Table 1
Parameter estimates, their standard errors and coverage probability associated with GINAR(1) ${}_{s}$ model

Table 5
Point forecast and HPP intervals for Compylobacteriosis data $\gamma=$ 0.95