INARMA Modeling of Count Time Series

Abstract

While most of the literature about INARMA models (integer-valued autoregressive moving-average) concentrates on the purely autoregressive INAR models, we consider INARMA models that also include a moving-average part. We study moment properties and show how to efficiently implement maximum likelihood estimation. We analyze the estimation performance and consider the topic of model selection. We also analyze the consequences of choosing an inadequate model for the given count process. Two real-data examples are presented for illustration.

1. Introduction

When dealing with stationary real-valued time series, the autoregressive moving-average (ARMA) models constitute a popular baseline model [1]. Besides the special case of purely autoregressive models, also full ARMA models are commonly used in practice, i.e., ARMA models that include a moving-average part (of typically low order). The integer-valued counterpart to the ARMA model, the INARMA model for ARMA-like count time series, dates back to McKenzie [2], Al-Osh and Alzaid [3], see Chapters 2 and 3 in Weiß [4] for a recent survey. In contrast to the ordinary ARMA models, research and applications nearly exclusively concentrate on the purely autoregressive INAR models (a few works also consider pure moving-average-type models, i.e., INMA models, like Brännäs and Hall [5], Brännäs and Quoreshi [6], Aleksandrov and Weiß [7]), but autoregressive INARMA models with an additional MA-component are rarely used in practice. This might be caused by the fact that there is only little work on stochastic properties of full INARMA models which, in turn, complicates the application of such models in practice. However, it seems to be mainly due to practical issues, e.g., it is not clear how to efficiently estimate the model parameters or to select the model order. While maximum likelihood (ML) estimation is easily done for the Markovian INAR$\left(p\right)$ models, it is assessed until now that “the inclusion of a moving average component renders maximum likelihood estimation infeasible” [8] (p. 1). The only likelihood-related approach discussed until now in the literature seems to be the Markov chain Monte Carlo (MCMC) method for INARMA processes developed by Neal and Subba Rao [9], Enciso-Mora et al. [10], which has been used for conducting inference on model parameters and forecast distributions as well as for model selection, also see Alzahrani et al. [11]. Furthermore, Enciso-Mora et al. [10] propose to use the Expectation Maximization (EM) algorithm for ML estimation, but only for the purely autoregressive INAR models. Wheatley et al. [12] conjecture the asymptotic equivalence between an INARMA process and an ARMA point process, which, in turn, would allow to adapt their EM scheme for ARMA point processes to the INARMA case. In the present article, however, we show that a direct numerical maximization of the log-likelihood function is tractable by considering the proposed efficient implementation of likelihood computation.

In Section 2, we give the definition of INARMA models and a more detailed discussion of those models, which are relevant for our study, namely INAR$\left(p\right)$ and INARMA(1,1) models. Then we turn to the question of model fitting. In Section 3, we show that it is possible to efficiently implement ML estimation also for the INARMA(1,1) model. Furthermore, we compare the performance of the ML approach for INAR$\left(p\right)$ and INARMA(1,1) models with simulations. Section 4 discusses the task of model selection across different types of INARMA models, and Section 5 analyzes the consequences of fitting the wrong type of INARMA model to the actual data-generating process (DGP). Our findings are also illustrated with two real-data examples in Section 6. Finally, we conclude in Section 7.

2. INARMA Models for Count Time Series

The basic idea behind INARMA models is to modify the ordinary ARMA recursion by replacing multiplications by so-called ‘binomial thinning’ operations $\{\rho \circ Z\}$ with $\rho \in [0,1]$, see [13], where

$$\rho \circ Z|Z\sim \mathrm{Bin}(Z,\rho ).$$

(1)

The conditional binomial distribution immediately implies that

$$E(\rho \circ Z)=\rho E\left(Z\right),\mathrm{Var}(\rho \circ Z)=\rho (1-\rho )E\left(Z\right)+{\rho }^2\mathrm{Var}\left(Z\right).$$

(2)

Using this integer-valued operation as a substitute for the ordinary multiplication, the formal definition of INARMA$(p,q)$ models is given by the recursion [14]

$$Y_t={\alpha }_1\circ Y_{t-1}+\dots +{\alpha }_p\circ Y_{t-p}+R_t+{\beta }_1\circ R_{t-1}+\dots +{\beta }_q\circ R_{t-q},$$

(3)

where $\left(R_t\right)$ constitutes a sequence of independent and identically distributed (i.i.d.) random variables, commonly referred to as the ‘innovations’. We denote the innovations’ mean by $E\left(R_t\right)=\lambda$ and the variance by Var($R_t)=\nu \lambda$. So $\nu =\mathrm{Var}\left(R_t\right)/E\left(R_t\right)$ expresses the dispersion ratio, which equals 1 in the case of Poisson innovations (equidispersion).

Regarding Equation (3), some caution is necessary. Since the binomial thinning operations are random operations, one has to carefully think about the joint distribution among the thinnings, and the joint distribution of the thinnings to the other random variables in Equation (3). To obtain feasible stochastic properties, one commonly formulates several independence assumptions. This is exemplified in the sequel by discussing the purely autoregressive INAR$\left(p\right)$ models as well as the INARMA(1,1) model.

2.1. A Brief Survey on INAR Models

After its introduction by McKenzie [2], Al-Osh and Alzaid [3], the INAR(1) model was extended to a pth-order INAR$\left(p\right)$ model by Alzaid and Al-Osh [15], Du and Li [16], but in two different ways. Both extensions follow the recursion $Y_t={\alpha }_1\circ Y_{t-1}+\dots +{\alpha }_p\circ Y_{t-p}+R_t$, but while the INAR$\left(p\right)$ model by Du and Li [16] assumes independence among all thinnings, independence to the innovations, and independence to ${\left(Y_s\right)}_{s<t}$ for the thinnings at time t, the model by Alzaid and Al-Osh [15] assumes a conditional multinomial distribution for the thinnings experienced by $Y_t$, in the sense that

$$({\alpha }_1\circ Y_t,\dots ,{\alpha }_p\circ Y_t)\sim \mathrm{Mult}(Y_t;{\alpha }_1,\dots ,{\alpha }_p).$$

Since only the INAR$\left(p\right)$ model by Du and Li [16] leads to the well-known Yule-Walker equations for the autocorrelation function (ACF) of an ordinary AR$\left(p\right)$ model, these are typically preferred in practice and also considered in the remaining article.

For the existence of the INAR$\left(p\right)$ process by Du and Li [16], ${\alpha }_{•}:={\sum }_{j=1}^p{\alpha }_j<1$ has to be assumed. Then one obtains a pth-order Markov process, where the transition probabilities are computed as a convolution between p binomial distributions and the innovations’ distribution. In particular, one obtains

$$P(Y_t=k|Y_{t-1}=l)=\sum _{j=0}^{min\{k,l\}}\left(\genfrac{}{}{0pt}{}{l}{j}\right){\alpha }^j{(1-\alpha )}^{l-j}·P(R_t=k-j)$$

(4)

for the INAR(1) model, and

$$\begin{array}{cc}\hfill P(Y_t=k|Y_{t-1}=l_1,Y_{t-2}=l_2)=& {\displaystyle \sum _{j_1=0}^{min\{k,l_1\}}\sum _{j_2=0}^{min\{k-l_1,l_2\}}\left(\genfrac{}{}{0pt}{}{l_1}{j_1}\right){\alpha }_1^{j_1}{(1-{\alpha }_1)}^{l_1-j_1}}\hfill \\ & {\displaystyle ·\left(\genfrac{}{}{0pt}{}{l_2}{j_2}\right){\alpha }_2^{j_2}{(1-{\alpha }_2)}^{l_2-j_2}·P(R_t=k-j_1-j_2)}\hfill \end{array}$$

(5)

for the INAR(2) model, see, e.g., Weiß [4]. Furthermore, the stationary marginal mean is given by $E\left(Y_t\right)=\lambda /(1-{\alpha }_{•})$, and the variance satisfies

$$\mathrm{Var}\left(Y_t\right)·\left(1-{\sum }_{i=1}^p{\alpha }_i\rho \left(i\right)\right)=E\left(Y_t\right){\sum }_{j=1}^p{\alpha }_j(1-{\alpha }_j)+\nu \lambda ,$$

(6)

also see Equation (2). The ACF $\rho \left(k\right)=\mathrm{Corr}(Y_t,Y_{t-k})$ is obtained from the ordinary AR$\left(p\right)$’s Yule-Walker equations

$$\rho \left(k\right)={\sum }_{i=1}^p{\alpha }_i\rho \left(\right|k-i\left|\right)\mathrm{for}k\ge 1.$$

(7)

In particular, the INAR(1) model satisfies $\rho \left(k\right)={\alpha }_1^k$.

2.2. Moment Properties of INARMA(1,1) Model

While the stochastic properties for INAR models as surveyed in Section 2.1 are well known among practitioners, such properties are not readily available for the INARMA(1,1) model. The INARMA(1,1) model is defined as

$$Y_t={\alpha }_1\circ Y_{t-1}+{\beta }_1\circ R_{t-1}+R_t,$$

(8)

where the ARMA parameters ${\alpha }_1,{\beta }_1\in [0,1]$ are used within the binomial thinning operations, and these are assumed to be executed independently of each other. Hence,

$$\mathrm{Cov}({\alpha }_1\circ W,{\beta }_1\circ Z)={\alpha }_1{\beta }_1\mathrm{Cov}(W,Z).$$

(9)

This model is a special case of the GINARMA$(p,q)$ model by Dion et al. [14]. For ${\beta }_1=0$, model (8) reduces to the INAR(1) model by McKenzie [2], Al-Osh and Alzaid [3], and for ${\alpha }_1=0$, it reduces to the INMA(1) model recently surveyed by Aleksandrov and Weiß [7]. But model (8) differs from the INARMA models considered by Dungey et al. [8], McKenzie [17], Bracher [18]. The INARMA(1,1)-like model by McKenzie [17] includes the AR-component only with lag 2 (i.e., with $Y_{t-2}$ instead of $Y_{t-1}$), and the one by Dungey et al. [8] applies the ${\beta }_1$-thinning to $R_t$ instead of $R_{t-1}$. In Bracher [18], the MA-component is constructed in a different way.

The relation between the counting series $Y_t$ and the innovation terms $R_t$ is summarized through the following lemma.

Lemma 1.

For the INARMA(1,1) model (8),

(a) 

$$Cov(Y_t,R_{t+h})=\left\{\begin{array}{cc}Var\left(R_t\right)\hfill & h=0,\hfill \\ ({\alpha }_1+{\beta }_1)Var\left(R_t\right)\hfill & h=-1,\hfill \\ 0\hfill & otherwise.\hfill \end{array}\right.$$

(b) 

$Cov({\alpha }_1\circ Y_{t-1},{\beta }_1\circ R_{t-1})={\alpha }_1{\beta }_{1}Cov(Y_{t-1},R_{t-1})={\alpha }_1{\beta }_{1}Var\left(R_t\right)$.

Using Lemma 1, it is proven in Appendix A, under stationary moments condition, that

$$E\left(Y_t\right)=\frac{(1+{\beta }_1)\lambda }{1-{\alpha }_1},$$

(10)

also see Corollary 2 in Dion et al. [14]. Analogously, we obtain

$$\mathrm{Var}\left(Y_t\right)=\frac{\left[{\alpha }_1(1+{\beta }_1)+{\beta }_1(1-{\beta }_1)\right]\lambda +\left[{\beta }_1({\beta }_1+2{\alpha }_1)+1\right]\nu \lambda }{1-{\alpha }_1^2}.$$

(11)

The lag-h autocovariance with $h\ge 1$ is given by

$$\mathrm{Cov}(Y_t,Y_{t+h})={\alpha }_1^h\mathrm{Var}\left(Y_t\right)+{\alpha }_1^{h-1}{\beta }_1\nu \lambda ,$$

(12)

see Appendix A for the derivation.

3. ML Estimation for INARMA(1,1) Models

ML estimation is easily implemented for the INAR models surveyed in Section 2.1, because these are Markov models such that the (conditional) likelihood function

$$L\left(\mathit{\theta }\right)=P(Y_T=y_T,\dots ,Y_2=y_2|Y_1=y_1)$$

factorizes as

$$L\left(\mathit{\theta }\right)={\prod }_{t=p+1}^{T}P(Y_t=y_t|Y_{t-1}=y_{t-1},\dots ,Y_{t-p}=y_{t-p}),$$

see Appendix B.2.1 in Weiß [4]. The required transition probabilities are computed like in Equations (4) and (5). An estimate for the parameter vector $\mathit{\theta }$ (which contains the thinning parameters ${\alpha }_1,\dots ,{\alpha }_p$ plus all parameters related to the innovations) is computed by numerically maximizing (the logarithm of) $L\left(\mathit{\theta }\right)$. For the non-Markovian INARMA(1,1) model, in contrast, an efficient likelihood computation is much more demanding.

3.1. Efficient Implementation of ML Estimation

According to (conditional) ML estimation, the estimate of the parameter vector $\mathit{\theta }$ (which contains the two thinning parameters ${\alpha }_1,{\beta }_1$ plus all parameters related to the innovations) is obtained as a maximizer of

$$L\left(\mathit{\theta }\right)=P(Y_T=y_T,\dots ,Y_2=y_2|Y_1=y_1).$$

For an efficient recursive computation of $L\left(\mathit{\theta }\right)$ (or the logarithm thereof, i.e., of $\ell \left(\mathit{\theta }\right)$), we adapt an approach known from ML estimation for Hidden-Markov models, see Zucchini et al. [19], Weiß [4]. Let us consider the probabilities

$$b_{kl}\left(t\right)=P(R_t=k,R_{t-1}=l;Y_t=y_t,\dots ,Y_2=y_2|Y_1=y_1).$$

Then $L\left(\mathit{\theta }\right)={\sum }_{k,l}{b}_{kl}\left(T\right)$ holds. At a first glance, it seems that this sum has to be taken for $k,l=0,\dots ,\infty$. However, because of the model recursion $Y_t={\alpha }_1\circ Y_{t-1}+{\beta }_1\circ R_{t-1}+R_t$, it is clear that $R_t\le Y_t$ for all t. So if $M=max\{y_T,\dots ,y_1\}$, then the $b_{kl}\left(t\right)$ only have to be computed for $k,l=0,\dots ,M$. Furthermore, $L\left(\mathit{\theta }\right)={\sum }_{k=0}^{y_T}{\sum }_{l=0}^{y_{T-1}}{b}_{kl}\left(T\right)$ holds.

The likelihood function is now computed recursively. We have

$$\begin{array}{cc}\hfill b_{kl}(t+1)=& \sum _{i=0}^{y_{t-1}}P(R_{t+1}=k,R_t=l,R_{t-1}=i;Y_{t+1}=y_{t+1},\dots |Y_1=y_1)\hfill \\ \hfill =& P(Y_{t+1}=y_{t+1}|R_{t+1}=k,R_t=l,Y_t=y_t)P(R_{t+1}=k)\hfill \\ & ·\sum _{i=0}^{y_{t-1}}P(R_t=l,R_{t-1}=i;Y_t=y_t,\dots |Y_1=y_1)\hfill \\ \hfill =& P({\alpha }_1\circ y_t+{\beta }_1\circ l+k=y_{t+1})P(R_{t+1}=k)\sum _{i=0}^{y_{t-1}}{b}_{li}\left(t\right)\hfill \\ \hfill =& P(R_{t+1}=k)\left(\sum _{i=0}^{y_{t-1}}{b}_{li}\left(t\right)\right)\hfill \\ & ·\sum _{y=0}^{y_{t+1}-k}\left(\genfrac{}{}{0pt}{}{y_t}{y}\right){\alpha }_1^y{(1-{\alpha }_1)}^{y_t-y}\left(\genfrac{}{}{0pt}{}{l}{y_{t+1}-k-y}\right){\beta }_1^{y_{t+1}-k-y}{(1-{\beta }_1)}^{l-y_{t+1}+k+y}.\hfill \end{array}$$

For initialization, we approximate

$$\begin{array}{cc}\hfill b_{kl}\left(2\right)=& P(R_2=k,R_1=l;Y_2=y_2|Y_1=y_1)\hfill \\ \hfill =& P(Y_2=y_2|R_2=k,R_1=l,Y_1=y_1)P(R_2=k)P(R_1=l|Y_1=y_1)\hfill \\ \hfill \approx & P(R_2=k)P(R_1=l|R_1\le y_1)\hfill \\ & ·\sum _{y=0}^{y_2-k}\left(\genfrac{}{}{0pt}{}{y_1}{y}\right){\alpha }_1^y{(1-{\alpha }_1)}^{y_1-y}\left(\genfrac{}{}{0pt}{}{l}{y_2-k-y}\right){\beta }_1^{y_2-k-y}{(1-{\beta }_1)}^{l-y_2+k+y}.\hfill \end{array}$$

From these derivations, an algorithmic scheme for log-likelihood computation can be derived. First note that the update step $t\mapsto t+1$ only requires the sums ${\sum }_{i=0}^{y_{t-1}}{b}_{li}\left(t\right)$ but not the individual $b_{li}\left(t\right)$. Hence, it suffices to define the $(M+1)$-dimensional vectors ${\mathit{a}}_t$ with entries $a_{t,l}={\sum }_{i=0}^{y_{t-1}}{b}_{li}\left(t\right)$, which may have non-zero entries only for $l=0,\dots ,y_t$. Next, the initialization step for $t=2$ can be included in the update step $t\mapsto t+1$ by defining $a_{1,l}=P(R_1=l|R_1\le y_1)$. Finally, computations can be simplified by defining the following matrices:

$$\mathbf{D}=\mathrm{diag}\left(P(R=0),\dots ,P(R=M)\right)$$

and ${\mathbf{Q}}_t={\left(q_{t,kl}\right)}_{k,l=0,\dots ,M}$ for $t=2,\dots ,T$ with

$$q_{t,kl}=\sum _{y=0}^{y_t-k}\left(\genfrac{}{}{0pt}{}{y_{t-1}}{y}\right){\alpha }_1^y{(1-{\alpha }_1)}^{y_{t-1}-y}\left(\genfrac{}{}{0pt}{}{l}{y_t-k-y}\right){\beta }_1^{y_t-k-y}{(1-{\beta }_1)}^{l-y_t+k+y}.$$

Then our algorithm becomes

$${\mathit{a}}_1,{\mathit{a}}_t=\mathbf{D}{\mathbf{Q}}_t{\mathit{a}}_{t-1}\mathrm{for}t=2,\dots ,T,L\left(\mathit{\theta }\right)={\sum }_{k=0}^{y_T}{a}_{T,k}={\mathbf{1}}^{\top }{\mathit{a}}_T.$$

Finally, to make computations numerically more stable, also see Weiß [4], Zucchini et al. [19], we define

$$\begin{array}{c}{w}_1={\mathbf{1}}^{\top }{\mathit{a}}_1,{\mathit{\varphi }}_1={\mathit{a}}_1/w_1;\hfill \\ {\mathit{u}}_t:=\mathbf{D}{\mathbf{Q}}_t{\mathit{\varphi }}_{t-1},\frac{w_t}{w_{t-1}}={\mathbf{1}}^{\top }{\mathit{u}}_t,{\mathit{\varphi }}_t={\mathit{u}}_t/\frac{w_t}{w_{t-1}}\mathrm{for}t=2,\dots ,T.\hfill \end{array}$$

The log-likelihood function is obtained as

$$\ell \left(\mathit{\theta }\right)=ln{w}_T=ln{w}_1+{\sum }_{t=2}^{T}ln\frac{w_t}{w_{t-1}}.$$

The log-likelihood function $\ell \left(\mathit{\theta }\right)$ is numerically maximized by using a standard optimization routine.

Remark 1.

It should be noted that the same recursive scheme could also be used for implementing ML estimation for the INMA$\left(1\right)$ model (i.e., where the AR-part is missing). Up to now, moment estimation or conditional least squares estimation are used for this model, see [5,7], and also the MCMC framework might be applied for this purpose, see [9,10,12]. If one wants to do ML estimation instead, the above scheme can be used together with one modification: the computation of $b_{kl}\left(t\right)$ and thus ${\mathbf{Q}}_t$ has to be simplified by setting ${\alpha }_1=0$. This implies that

$$\begin{array}{cc}\hfill q_{t,kl}=& P(Y_t=y_t|R_t=k,R_{t-1}=l,Y_{t-1}=y_{t-1})\hfill \\ \hfill =& P({\beta }_1\circ l+k=y_t)=\left(\genfrac{}{}{0pt}{}{l}{y_t-k}\right){\beta }_1^{y_t-k}{(1-{\beta }_1)}^{l-y_t+k}.\hfill \end{array}$$

The presented scheme is also easily modified to capture seasonality or trend. If these are incorporated by time-dependent thinning parameters or innovations parameters, respectively, as suggested by Freeland and McCabe [20] (p. 704), then one just has to modify ${\mathbf{Q}}_t$ or ${\mathbf{D}}_t=\mathrm{diag}\left(P(R_t=0),\dots ,P(R_t=M)\right)$ (this diagonal matrix would then also depend on time t) accordingly.

Furthermore, it can also be adapted to fit higher-order INARMA models. If the AR-order p is $>1$, we have to compute

$$\begin{array}{cc}\hfill q_{t,kl}=& P(Y_t=y_t|R_t=k,R_{t-1}=l,Y_{t-1}=y_{t-1},\dots ,Y_{t-p}=y_{t-p})\hfill \\ \hfill =& P({\alpha }_1\circ y_{t-1}+\dots +{\alpha }_p\circ y_{t-p}+{\beta }_1\circ l=y_t-k)\hfill \end{array}$$

as a convolution of $p+1$ binomial distributions. An MA-order $q>1$ is more cumbersome as it goes along with an increase of dimensionality. We then have to define

$$b_{k_0{k}_1\dots k_q}\left(t\right)=P(R_t=k_0,R_{t-1}=k_1,\dots ,R_{t-q}=k_q;Y_t=y_t,\dots ,Y_2=y_2|Y_1=y_1)$$

with $q+1$ subscripts $k_0,k_1,\dots ,k_q=0,\dots ,M$.

3.2. Performance of ML Estimation

We simulated INARMA processes of orders (1,0), (2,0) and (1,1), and with the innovations being either Poisson (Poi) or negative binomially (NB) distributed. In the latter case, the dispersion factor $\nu$ was chosen equal to 1.5. For the marginal means, we considered $\mu \in \{3,6\}$, and the dependence parameters were always chosen such that $\rho \left(1\right)\in \{0.35,0.70\}$. We set ${\alpha }_2=0.25$ for the INAR(2) models, and ${\beta }_1=0.25$ for the INARMA(1,1) models. The considered sample sizes are $T\in \{100,250,500,1000\}$. For each scenario, we simulated 1000 replications, and the model parameters were estimated with the ML approach described before (always choosing the appropriate model type).

Boxplots of the simulation results (and also tables with the means of the obtained estimates) are summarized in Appendix B.1. As can be seen from the boxplots, the respective ML estimates are certainly less biased with increasing T and also show decreasing dispersion for any of the considered models. However, the final sample properties differ a lot between different types of models and different types of parameters. One general observation is that the additional dispersion introduced by NB innovations deteriorates the estimation performance. We do not only observe increased bias for the innovations’ mean $\lambda$ and the dependence parameters ${\alpha }_1,{\alpha }_2,{\beta }_1$, also the estimation of $\nu$ itself suffers from rather large dispersion. Furthermore, the distribution of $\widehat{\nu }$ is positively skewed and sometimes causes a strong overestimation of the true dispersion level $\nu$. In addition, for the innovations’ mean $\lambda$, we observe positive bias and skewness for small sample sizes, whereas the AR(1) parameter ${\alpha }_1$ is negatively biased. The same holds for the AR(2) parameter ${\alpha }_2$ (but with stronger bias than for ${\alpha }_1$), whereas ${\beta }_1$ is positively biased for small T. Compared across models, we observe much more dispersion for the estimates of the dependence parameters in the INARMA(1,1) case than in the purely autoregressive cases. It is also worth noting that the dispersion of ${\widehat{\alpha }}_1$ for INARMA(1,1) processes decreases with increasing mean $\mu$.

4. Model Selection for INARMA Processes

In the simulations of Section 3.2, we always fitted the true model type to the given count time series. In practice, however, the true model behind the DGP is not known and therefore has to be identified based on the available data. A widely used approach is to apply information criteria (IC) for this purpose, and especially Akaike’s and the Bayesian IC (AIC and BIC, respectively) are routinely used for this purpose see [21,22]. These criteria are computed together with the ML estimation of each candidate model, and that model is selected as the final one which minimizes the value of AIC or BIC, respectively. More details on these and further ICs can be found in Neath and Cavanaugh [21], Cavanaugh and Neath [22], Weiß and Feld [23]. In Weiß and Feld [23], the performance of these criteria was analyzed for count time series mainly generated by regression-type DGPs. They confirmed the consistency of the BIC in their study, but the actual rates of correct model identifications for smaller sample sizes T were often best for the AIC. Another related study is that of Zhu et al. [24]. They used AIC and BIC for selecting the order of the components of their mixture autoregressive Poisson regression model, but found out that these ICs “do not give a very satisfactory result” in this context.

Since model selection across INARMA models was not considered in Weiß and Feld [23], we also analyzed the AIC’s and BIC’s performance in our simulation study. For this purpose, we used all of the six considered models as possible candidate models for any of the simulated count time series. The obtained numbers of model selections (out of 1000 replications) are tabulated in Appendix B.2. Some of the conclusions found by Weiß and Feld [23] are confirmed also here. The BIC’s ability for identifying the correct model always improves with increasing T, whereas the AIC for the smallest model, i.e., for the INAR(1) model in our comparison, stabilizes at a rate below 80%. On the other hand, it does most often better than the BIC for smaller sample sizes such as $T\le 250$ (with exceptions only for INAR(1) DGPs). Besides looking at the correct identifications, it is also interesting to study the possible mis-identifications. It becomes clear that the purely autoregressive INAR models are mainly confused among themselves (by either choosing the wrong model order or the wrong distribution family), but we rarely observe an erroneous identification as INARMA(1,1). For an INARMA(1,1) DGP, in contrast, there is a very large risk of being mis-identified as an INAR(1) process. For $T=100$ and also for larger T if $\rho \left(1\right)=0.35$, such a mis-identification is even the most frequent decision, especially if using the BIC. This clearly differs from the INAR case, where the model order of the INAR(2) model is correctly identified in the majority of cases for any scenario. So not only the parameter estimation of an adequately chosen INARMA(1,1) model requires rather large sample sizes, these are also required for being able to correctly identify an INARMA(1,1) model at all. This shows that model selection should not solely rely on an information criterion, but further diagnostics should be done (e.g., a comparison of properties of the fitted models with the corresponding sample properties, see Section 6 below). At this point, it deems appropriate to recall the “two-units rule” for interpreting AIC and BIC, see Tables 1 in [21,22]. It says that if some candidate model’s AIC (BIC) differs from the smallest AIC (BIC) by a difference $\le 2$, it should also be considered as a “viable candidate” [22] (p. 6).

5. On Properties of (Mis-)Fitted Models

The previous Section 4 showed that there is a considerable risk of choosing the wrong model for the actual count DGP. Such a (possibly falsely) chosen model is then used for interpreting the data, for forecasting future values, or for setting up a control chart for progressive process monitoring, see Weiß [4]. Hence, it is important to ask for the consequences of such a possible misfit of the DGP’s true model. To answer this question (to some extent), we studied important stochastic properties of the fitted models and compared them to the respective properties of the true DGP. More precisely, we focused on the marginal mean $\mu$ and dispersion ratio ${\sigma }^2/\mu$ as well as on the ACF values $\rho \left(1\right),\rho \left(2\right),\rho \left(3\right)$, see Section 2 for the required formulae. Appendix B.3 provides tables of the means of these properties.

The consequences of overfitting a Poi-INAR(1) DGP are quite moderate, a slight effect on $\rho \left(2\right)$ has to be noted. Things change if the INAR(1) process has NB innovations. If then falsely fitting a Poi model, we are not only unable to capture the apparent overdispersion ${\sigma }^2/\mu >1$, also the ACF values are clearly undervalued. If erroneously fitting an INAR(1) or INARMA(1,1) model to a Poi-INAR(2) DGP, we severely underrate the actual ACF values, and this becomes even worse in the presence of overdispersion (i.e., NB-INAR(2) DGP). While analogous deviations are also observed for INARMA(1,1) DGPs, they appear to be less pronounced than in the INAR(2) case. Here, we usually underrate $\rho \left(1\right)$ but overrate $\rho \left(2\right),\rho \left(3\right)$. For the NB-INARMA(1,1) DGP, it is even worse to falsely choose a Poi-INARMA(1,1) model than an NB-INAR(1) or NB-INAR(2) model. So altogether, it appears that an inappropriate choice of the dispersion structure has the worst effect on the quality of the fitted model.

6. Real-Data Examples

Let us discuss two real-data examples, which both are time series $y_1,\dots ,y_T$ of monthly counts. The first time series consists of monthly counts of crime offense reports regarding burglaries from the 43th police car beat in Pittsburgh (1990–2001, so $T=144$). It is available at the Forecasting Principles website, http://www.forecastingprinciples.com/index.php/crimedata. The second time series consists of monthly counts of claims caused by burn related injuries in the heavy manufacturing industry (1987–1994, so $T=96$), see Example 2.5.1 in Freeland [25] for further details. Plots of both time series are provided by Figure 1.

Both time series exhibit a quickly decaying ACF and a moderate extent of empirical overdispersion (about 30%), see the first row in

Table 1. Upper part: Sample properties (italic font) and properties of maximum likelihood (ML)-fitted models for monthly counts of burglaries (left) and claims (right). Lower part: Results from parametric bootstrap experiment with 10,000 replications.

	Burglaries Counts:						Claims Counts:
	AIC	$\mathit{\mu }$	${\mathit{\sigma }}^{\mathbf{2}}/\mathit{\mu }$	$\mathit{\rho }\left(\mathbf{1}\right)$	$\mathit{\rho }\left(\mathbf{2}\right)$	$\mathit{\rho }\left(\mathbf{3}\right)$	AIC	$\mathit{\mu }$	${\mathit{\sigma }}^{\mathbf{2}}/\mathit{\mu }$	$\mathit{\rho }\left(\mathbf{1}\right)$	$\mathit{\rho }\left(\mathbf{2}\right)$	$\mathit{\rho }\left(\mathbf{3}\right)$
Sample	–	4.319	1.271	0.255	0.014	0.040	–	8.604	1.320	0.452	0.362	0.213
Poi-INAR(1)	643.7	4.311	1.000	0.210	0.044	0.009	490.5	8.667	1.000	0.396	0.157	0.062
NB-INAR(1)	641.9	4.312	1.264	0.238	0.057	0.013	490.4	8.670	1.236	0.426	0.181	0.077
Poi-INAR(2)	646.7	4.309	1.000	0.208	0.043	0.009	487.2	8.610	1.063	0.419	0.318	0.183
NB-INAR(2)	644.6	4.309	1.273	0.236	0.056	0.013	488.1	8.610	1.250	0.444	0.339	0.202
Poi-INARMA(1,1)	643.9	4.316	1.030	0.248	0.024	0.002	492.5	8.666	1.007	0.398	0.155	0.061
NB-INARMA(1,1)	642.7	4.319	1.273	0.266	0.019	0.001	492.4	8.682	1.264	0.418	0.175	0.073
Bootstrap for	Fitted NB-INARMA(1,1):						Fitted NB-INAR(2):
lower quartile		4.153	1.149	0.200	−0.052	−0.070		8.188	1.048	0.331	0.218	0.069
median		4.313	1.256	0.254	0.006	−0.009		8.594	1.190	0.407	0.295	0.151
upper quartile		4.479	1.371	0.306	0.065	0.049		9.031	1.357	0.480	0.368	0.233
standard error		0.247	0.167	0.077	0.086	0.086		0.634	0.235	0.111	0.112	0.121

, such that INAR(1), INAR(2) and INARMA(1,1) models with either Poisson or NB innovations are reasonable candidate models. These were fitted to both time series using ML estimation, as outlined in Section 3. Since both time series are rather short, we use the AIC for (initial) model selection; our analyses in Section 4 showed that the BIC most often has very low rates of correct model identification in such a case. The obtained AIC values as well as relevant stochastic properties of the fitted models (in analogy to Section 5) are summarized in the upper part of Table 1.

For the burglaries counts, the AIC selects the NB-INAR(1) model, but the AIC of the NB-INARMA(1,1) model is also quite low. Actually, also the corresponding Poi-counterparts have AICs that deviate from the lowest AIC by not more than two units, only the INAR(2) models violate this rule. If we now compare the tabulated properties of the fitted models with their corresponding sample counterparts (the latter are printed in italic font), the NB-INARMA(1,1) model shows the least deviations. The NB-INAR(1) model, in contrast, does not only show a smaller mean $\mu$ and dispersion ratio ${\sigma }^2/\mu$ than observed in the sample, also the ACF values show more deviations (and the Poi-models do even worse). So actually, the NB-INARMA(1,1) model provides the best fit to the burglaries counts data. As a further check for model adequacy, we computed an acceptance envelope in the spirit of Tsay [26] for the properties considered in Table 1. For this purpose, a parametric bootstrap with 10,000 replications for the fitted NB-INARMA(1,1) model was done. From the sample properties obtained for each simulation run, quartiles and standard errors for the sample statistics have been computed. It can be seen that the values for $\widehat{\mu }$, ${\widehat{\sigma }}^2/\widehat{\mu }$, $\widehat{\rho }\left(1\right)$, $\widehat{\rho }\left(2\right)$, $\widehat{\rho }\left(3\right)$ as computed from the time series (first line of Table 1) are always within the respective lower and upper quartile in the lower part of Table 1, so they do not contradict the fitted NB-INARMA(1,1) model.

But why did the AIC select the NB-INAR(1) model? Here, it is important to recall the discussion in Section 4. If the true DGP would be NB-INARMA(1,1), then the AIC has a very low rate of correct model identification for small T. If we consider the scenario $\mu =3$, $\rho \left(1\right)=0.35$ and $T=100$ in Appendix B.2, which is reasonably close to our data example, then the AIC falsely selects an NB-INAR(1) model in the majority of cases. So the AIC (and even more the BIC) is not a reliable tool for model selection for such small T. It gives rough indication in the sense that Poisson models do not work well, or that an INAR(2) model does not improve over an INAR(1) model. However, the decision for the final model should also consider further aspects like the stochastic properties of the fitted models.

An analogous conclusion applies to the claims counts data. This time, the AIC prefers the Poi-INAR(2) model against the remaining candidate models, but the NB-INAR(2) model also has an AIC satisfying the two-units rule. However, both the dispersion ratio and the ACF values are much smaller for the fitted Poi-INAR(2) model than for the sample itself, whereas the fitted NB-INAR(2) model gives a clearly better agreement (Table 1 also provides some results from a parametric bootstrap for the fitted NB-INAR(2) model with 10,000 replications, which do not contradict the adequacy of this model). Hence, it appears that the AIC was misleading also for this data example. Looking into Appendix B.2, DGP NB-INAR(2) with $\mu =6$, $\rho \left(1\right)=0.35$ and $T=100$, we see a strong tendency of the AIC for falsely selecting a Poi-INAR(2) model (and this becomes even worse with increasing $\rho \left(1\right)$). So we conclude again that for short count time series, the AIC (and especially the BIC) should be used with caution. It gives a rough orientation regarding the correct model type (the choice of an INAR(2) model appears to be justified), but the final model selection should also use further criteria like a comparison of model properties.

7. Conclusions

We derived an efficient scheme for ML estimation of INARMA(1,1) models, and we also discussed possible extensions to higher-order models, or to INARMA models with seasonality or trend. Then we compared the INARMA(1,1) model to the INAR(1) and INAR(2) model (these three models constitute a reasonable set of candidate models for many applications) in several respects. The performance of ML estimation is generally rather good, but a small-sample bias might be observed especially in the presence of overdispersion. The BIC is consistent for model selection across these INARMA models, but it shows a worse performance than the AIC for small sample sizes (such as $T\le 250$). In addition, the AIC has to be treated with caution for very short time series (such as $T=100$): It gives a rough orientation but should always be complemented by further selection criteria. The consequences of fitting the wrong model to the actual DGP might be particularly severe in the presence of overdispersion, where a misfit of the dispersion behavior also causes a misfit of the ACF values.