# Integer-valued autoregressive models based on quasi Pólya thinning operator

Jean Peyhardi

► To cite this version:

Jean Peyhardi. Integer-valued autoregressive models based on quasi Pólya thinning operator. Statistical Inference for Stochastic Processes, In press, 27 (3), pp.813-838. 10.1007/s11203-024-09316-3. hal-04691912

HAL Id: hal-04691912

https://hal.science/hal-04691912v1

Submitted on 11 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Integer-valued autoregressive models based on quasi Pólya thinning operator

Jean Peyhardi¹

¹IMAG, Université de Montpellier, CNRS, Montpellier, 34090, France.

Contributing authors: jean.peyhardi@umontpellier.fr;

## Abstract

Autoregressive models adapted to count time series have received less attention than their classical counterparts for continuous time series. The main approach is based on thinning stochastic operation that preserves the discrete nature of the variable between successive times. The binomial thinning operator is the most popular and the Poisson distribution emerges as the natural choice for the error distribution of such an autoregressive counting process. The present paper introduces the quasi Pólya thinning operator, that includes the binomial thinning operator as a special case. The family of additive modified power series distribution is defined and is shown to be the natural choice for the error distribution of such a counting process. We obtain the most general class of integer-valued autoregressive models of order 1 (INAR(1) models) with margins having analytic form and the property of closure under convolution introduced by Joe (1996). It includes the usual cases of Poisson and generalized Poisson margins, but also the less usual cases of binomial and negative binomial margins and the new case of generalized negative binomial margin. These models cover a high range of dispersion that are strictly ordered from the binomial case to the generalized negative binomial case. Asymptotic normality of the maximum likelihood estimator (MLE) for such INAR(1) models is obtained. Finally, the proposed INAR(1) models are applied on simulated and real datasets.

**Keywords:** Binomial thinning operator, Quasi Pólya distribution, INAR model, modified power series distribution
# 1 Introduction

A popular class of models for count time series is the class of INAR(1) models based on thinning operator. A thinning operator is a stochastic operator that shrinks a random count variable into a smaller one. This kind of operator is thus useful to define $X_t$ from $X_{t-1}$, ensuring that $X_t$ is still an integer. The most popular thinning operator is the binomial thinning operator, usually defined as follows. Let $X$ denote a random count variable and $\rho \in (0, 1)$ a given parameter, then $X' \stackrel{d}{=} \rho \circ X$ denotes the thinned version of $X$. It is defined as the random sum $X' = \sum_{i=1}^{X} Z_i$ where $Z_i \sim \mathcal{B}(\rho)$ are independent and identically distributed random Bernoulli variables with parameter $\rho$. An INAR(1) model based on binomial thinning operator is defined by the equation

$$X_t = \rho \circ X_{t-1} + \varepsilon_t, \quad (1)$$

where the residual $\varepsilon_t$ is assumed to be independent of $(X_s)_{s \le t-1}$ and all the residuals $(\varepsilon_t)_{t \in \mathbb{N}}$ are independent and identically distributed according to the distribution $\mathcal{L}_\varepsilon$. The most natural choice of error distribution $\mathcal{L}_\varepsilon$ associated with the binomial thinning operator is the Poisson distribution, since it is closed under this thinning operation and under convolution. It implies that equation (1) defines a stationary process with Poisson marginal distribution, as shown by McKenzie (1985). But count time series often present over dispersion and the Poisson distribution becomes inappropriate. In the case of negative binomial margin, McKenzie (1987) explicitly derived an expression for the error distribution, but for other margins it becomes more difficult.

It is possible to use other thinning operator to obtain the desired marginal distribution. Several extensions of the binomial thinning operator have been proposed in the literature. Latour (1998) introduced a generalized thinning operator by relaxing the assumption of Bernoulli distribution for the latent variables $Z_i$. For instance, the specific case of geometric distribution $Z_i \sim \mathcal{G}(\rho)$, has been studied by Ristić et al (2009) leading to the negative binomial operator denoted by $\rho *$. But assuming a full range for $Z_i$ in $\mathbb{N}$ instead of $\{0, 1\}$ implies that the constraint $\rho * X \le X$ does not hold anymore, and the operation $\rho \circ$ could not been viewed as a thinning operation. Another extension of the binomial thinning operator is the random coefficient thinning operator studied by Zheng et al (2007), where the variables $Z_i$ given $\rho$ are assumed to follow independent Bernoulli $\mathcal{B}(\rho)$ and $\rho$ is assumed to be a random variable. Assume for instance a beta distribution for $\rho$ leads to the beta binomial thinning operator. It could be shown that the natural choice for the error distribution in this case is the negative binomial distribution (Joe, 1996).

The present paper proposes to view the binomial thinning operation as a compound distribution operation, i.e., $\rho \circ X$ given $X = n$ follows a binomial distribution, denoted by $\mathcal{B}_n(\rho)$. In other words, if $\mathcal{L}$ denotes the distribution of $X$ then $\rho \circ X$ follows the compound distribution $\mathcal{B}_n(\rho) \wedge \mathcal{L}$, where the notation $\wedge_n \mathcal{L}$ means that the parameter $n$ is considered as a random variable following the distribution $\mathcal{L}$. It is then possible to replace the binomial distribution by another distribution supported on $\{0, \dots, n\}$, such as the hypergeometric distribution for instance or the quasi binomial distribution. These two examples have been studied by Al-Osh and Alzaid (1991) and Alzaid and
Al-Osh (1993) respectively. They have shown that if a binomial (resp. generalized Poisson) distribution is assumed for the residuals, then the process is stationary with the same distribution. These models, such as the binomial thinning with Poisson margin or the beta binomial thinning with negative binomial margin, fall into the class of INAR(1) models with margins in the convolution-closed infinitely divisible class described by Joe (1996). Another INAR(1) model with negative binomial margin and negative binomial error distribution has been proposed by Al-Osh and Aly (1992), using an operator based on a sum of geometrically distributed variables with a number of terms given by a binomial thinning operation.

In fact, the hypergeometric, binomial and beta binomial thinning operators have previously been studied in the framework of damage models; see (Patil and Ratnaparkhi, 1975) for a review. In this framework, a random variable *N* is split into two random variables $Y_1$ and $Y_2$, i.e., such that $N = Y_1+Y_2$. The conditional distribution of $Y_1$ given $N = n$ is called damage distribution, e.g., the binomial distribution. The variable $Y_1$ can be viewed as a damaged version of the original variable *N*. Two properties of such a model have been intensively studied: the closure under damage operation and the characterization of the original distribution of *N* through the condition of independence between $Y_1$ and $Y_2$. For instance, Rao (1965) showed that binomial, Poisson and negative binomial distribution are closed under the binomial thinning operation. It was generalized by Janardan and Rao (1986) which demonstrated the closure of generalized Pólya Eggenberger distributions under quasi Pólya thinning operation.

Based on this closure property, we propose a general class of INAR(1) models using the quasi Pólya distribution as thinning operator; see (Peyhardi, 2023) for a study of this stochastic operator. The present paper relates the closure properties under thinning operator to the marginal distribution of the stationary INAR(1) process. The proposed class of INAR(1) models includes the cases of binomial and quasi binomial thinning operators, already studied, the less usual cases of hypergeometric and beta binomial thinning operators and the new case of quasi beta binomial thinning operator. All properties of these models will be studied in a general way based on the class of quasi Pólya distributions. For instance, it will be shown that the natural choice for the error distribution is the class of additive modified power series distribution that lead to a stationary process with the same marginal distribution. Recently, another general approach was proposed by Sellers et al (2020) based on the Conway-Maxwell-Binomial (CMB) thinning operator with Conway-Maxwell-Poisson (CMP) as marginal distribution. It includes the three cases of hypergeometric, binomial and beta binomial thinning operators. But the first case is obtained as a limit, not really as a subcase. Moreover, the CMB and CMP distributions do not have closed form and thus need the computation of a normalizing constant. The present paper introduces the most general class of INAR(1) models with margins having analytic form and the property of closure under convolution introduced by Joe (1996).

The paper is organized as follows. Section 2 introduces some preliminary notations and definitions about the class of quasi Pólya distributions and its counterpart, the class of additive modified power series distributions. Section 3 presents the quasi Pólya thinning operator and its properties. Section 4 defines INAR(1) models based on quasi Pólya thinning operator. It focuses on the natural choice for the error distribution that
is preserved at stationary state for the marginal distribution. The five specific INAR(1) models based on quasi Pólya thinning operator are described and the standardized Pearson residuals are given in order to evaluate the model adequacy. Section 5 gives regularity conditions under which the asymptotic normality of the MLE is obtained. Finally, Section 6 presents some applied results of INAR(1) models on simulated and real datasets.

## 2 Preliminary notations and definitions

Let us introduce the parametric sequence

$$a_{\theta}^{[c,d]}(n) = \frac{\theta(\theta + dn)^{(n;c)}}{(\theta + dn)n!}, \qquad (2)$$

where $n \in \mathbb{N}$, $\theta \in \Theta$ (with $\Theta = \mathbb{N}^*$ when $c = -1$, $d = 0$ and $\Theta = \mathbb{R}_+^*$ otherwise), $c \in \{-1, 0, 1\}$, $d \in \mathbb{R}_+$ and $m^{(n;c)} := \prod_{k=0}^{n-1}(m+ck)$, where $\mathbb{N}^*$ (resp. $\mathbb{R}_+^*$) denotes the set of positive integers (resp. positive reals). Let $\binom{\theta}{n}$ denote the binomial coefficient using the gamma function, i.e., $\binom{n+\theta-1}{n} = \frac{\Gamma(n+\theta)}{\Gamma(n+1)\Gamma(\theta)}$ that is well-defined for continuous values of $\theta \in \mathbb{R}_+^*$. One can distinguish five equivalence classes for the parametric sequence $a_{\theta}^{[c,d]}(n)$ according to the $(c,d)$ values, with the following representative elements:

1. $a_{\theta}^{[-1,0]}(n) = \binom{\theta}{n}$

2. $a_{\theta}^{[0,0]}(n) = \theta^n/n!$

3. $a_{\theta}^{[1,0]}(n) = \binom{n+\theta-1}{n}$

4. $a_{\theta}^{[0,d]}(n) = \frac{\theta(\theta+dn)^{n-1}}{n!}$ with $d > 0$

5. $a_{\theta}^{[1,d]}(n) = \frac{\theta}{\theta+dn} \binom{\theta+dn+n-1}{n}$ with $d > 0$

Remark that the second case is included in the fourth case if $d = 0$ and the first and third cases are included in the fifth case if $d = -1$ and $d = 0$ respectively. The case $c = -1$ and $d \neq 0$ is omitted since $a_{\theta}^{[-1,d+1]}(n) = a_{\theta}^{[1,d]}(n)$ for all $d \ge -1$. In each case, it can be shown that the parametric sequence $a_{\theta}^{[c,d]}(n)$ respects the two following properties:

(P1) Additive closure under convolution:

$$a_{\theta}^{[c,d]} * a_{\gamma}^{[c,d]} = a_{\theta+\gamma}^{[c,d]}$$

(P2) Modified power series convergence: there exists a function $g$ such that

$$\sum_{k=0}^{\infty} a_{\theta}^{[c,d]}(k)\{g(\alpha)\}^k = h_{\theta}(\alpha)$$

for all $\theta \in \Theta$ and $\alpha \in (0, R)$. Note that $R$ is the radius of convergence of the power series when $g$ is the identity function.
The property (P1) is exactly the convolution identity of Abel when $c = 0$ and those of Hagen-Rothe when $c \neq 0$; see (Chu, 2010) for demonstration. The symbol * here denotes the convolution of sequences, i.e., $(a_\theta * a_\gamma)(n) = \sum_{k=0}^n a_\theta(k)a_\gamma(n-k)$. The property (P2) corresponds to usual power series when $c = 0$ and their Lagrangian expansion when $c \neq 0$: see (Consul and Famoye, 2006) for the demonstration of the case $c \neq 0$ and see Table A2 in Appendix A for details about the function $g$ and the normalizing constant $h$ in the five cases. These properties allow defining the following two classes of distributions.

**Definition 1.** A discrete non-negative variable $X$ is said to follow a quasi Pólya distribution if its probability mass function (pmf) is given by

$$P(X = k) = \frac{a_{\theta}^{[c,d]}(k)a_{\gamma}^{[c,d]}(n-k)}{a_{\theta+\gamma}^{[c,d]}(n)} \mathbb{1}_{k \le n},$$

for some $(\theta, \gamma) \in \Theta^2$ and $n \in \mathbb{N}$. This distribution will be denoted by $qP_n^{[c,d]}(\theta, \gamma)$.

Property (P1) ensures that we obtain a proper distribution. This class of distributions includes for instance the hypergeometric distribution when $c = -1$ and $d = 0$, or the binomial distribution when $c = 0$ and $d = 0$. More generally, this is the family of quasi Pólya distributions, introduced by Janardan (1973) as urn model with predetermined strategies. The different examples of quasi Pólya distributions are summarized in Table 1 and more details are given in Table A1 of Appendix A. The moments of quasi Pólya distributions are recalled in the following Lemma.

**Lemma 1.** Let $Y$ follow a quasi Pólya distribution $qP_n^{[c,d]}(\theta, \gamma)$ and $Z = n - Y$, then

$$\bullet E[Y] = \frac{\theta}{\theta + \gamma} n$$

$$\bullet V(Y) = \frac{\theta\gamma}{\theta + \gamma} \left[ \frac{n^2}{\theta + \gamma} - \sum_{j=1}^{n-1} \frac{n^{(j+1; -1)}(d+c)^{j-1}}{\{\theta + \gamma + dx_{t-1} + (n-j+2)c\}^{(j;c)}} \right]$$

$$\bullet \text{ If } d = 0, Y \text{ follows a simple Pólya distribution } P_n^{[c]}(\theta, \gamma), \text{ then}$$

$$V(Y) = \frac{\theta\gamma}{(\theta + \gamma)(\theta + \gamma + c)} n + c \frac{\theta\gamma}{(\theta + \gamma)^2(\theta + \gamma + c)} n^2$$
$$Cov(Y, Z) = - \frac{\theta\gamma}{(\theta + \gamma)(\theta + \gamma + c)} n - c \frac{\theta\gamma}{(\theta + \gamma)^2(\theta + \gamma + c)} n^2$$

The two first points of this Lemma have been shown by Janardan and Schaeffer (1977) and the third by Johnson et al (1997).

**Definition 2.** A discrete non-negative variable $X$ is said to follow an additive modified power series distribution if its pmf is given by

$$P(X = k) = \frac{a_{\theta}^{[c,d]}(k)\{g(\alpha)\}^k}{h_{\theta}(\alpha)}$$

for some $\theta \in \Theta$ and $\alpha \in (0, R)$. This distribution will be denoted by $APS^{[c,d]}(\theta, \alpha)$.
<table><caption>**Table 1** Equivalence classes of quasi Pólya distributions $qP_n^{[c,d]}(\theta, \gamma)$ according to $(c, d)$ values. The cell with quasi hypergeometric distribution is greyed out since it is equivalent to the quasi beta binomial distribution, more precisely we have $qH_n^{[d]}(\theta, \gamma) = q\beta B_n^{[d-1]}(\theta, \gamma)$. Note that $\rho = \frac{\theta}{\theta+\gamma}$ and $d' = \frac{d}{\theta+\gamma}$.</caption><thead><tr><th></th><th>d = 0</th><th>d > 0</th></tr></thead><tbody><tr><td>c &lt; 0</td><td>hypergeometric<br/>H<sub>n</sub>(θ, γ)</td><td>quasi hypergeometric<br/>qH<sub>n</sub><sup>[d]</sup>(θ, γ)</td></tr><tr><td>c = 0</td><td>binomial<br/>B<sub>n</sub>(ρ)</td><td>quasi binomial<br/>qB<sub>n</sub><sup>[d]</sup>(ρ)</td></tr><tr><td>c &gt; 0</td><td>beta binomial<br/>βB<sub>n</sub>(θ, γ)</td><td>quasi beta binomial<br/>qβB<sub>n</sub><sup>[d'</sup>(θ, γ)</td></tr></tbody></table>

Property (P2) ensures that we obtain a proper distribution. This is a subclass of the modified power series distributions introduced by Gupta (1974). It includes for instance the binomial distribution when $c = -1$ and $d = 0$, or the Poisson distribution when $c = 0$ and $d = 0$. The different examples of additive modified power series distributions are summarized in Table 2 and more details are given in Table A2 of Appendix A. The mean and variance of such distributions are expressed in the following Lemma.

**Lemma 2.** Let $X \sim APS^{[c,d]}(\theta, \alpha)$ then there exists two functions $f_1$ and $f_2$ such that the expectation and variance of $X$ can be decomposed as follows:

$$E[X] = \theta f_1(\alpha),$$

$$V(X) = \theta f_2(\alpha).$$

*Proof.* Let $X \sim APS^{[c,d]}(\theta, \alpha)$. According to Gupta (1974) we have

$$\mu = E[X] = \frac{g(\alpha) h'_{\theta}(\alpha)}{g'(\alpha) h_{\theta}(\alpha)},$$

$$V(X) = \frac{g(\alpha)}{g'(\alpha)} \mu'(\alpha),$$

where $g$ and $h$ are given in Table A2 of Appendix A. It is easily seen that

$$\frac{h'_{\theta}(\alpha)}{h_{\theta}(\alpha)} = \begin{cases} \frac{\theta}{1+\alpha} & \text{if } c = -1 \\ \theta & \text{if } c = 0 \\ \frac{\theta}{1-\alpha} & \text{if } c = 1 \end{cases}$$

and thus with

$$f_1(\alpha) := \begin{cases} \frac{g(\alpha)}{g'(\alpha)(1+\alpha)} & \text{if } c = -1 \\ \frac{g(\alpha)}{g'(\alpha)} & \text{if } c = 0 \\ \frac{g(\alpha)}{g'(\alpha)(1-\alpha)} & \text{if } c = 1 \end{cases}$$

and $f_2(\alpha) := \frac{g(\alpha)}{g'(\alpha)} f_1'(\alpha)$ we obtain $E[X] = \theta f_1(\alpha)$ and $V(X) = \theta f_2(\alpha)$. $\square$
<table><caption>**Table 2** Equivalence classes of additive modified power series distributions APS<sup>[c,d]</sup>(θ, α) according to (c, d) values. The cell generalized binomial distribution GB<sup>[d]</sup>(θ, p) is greyed out, since it is equivalent to the generalized negative binomial distribution GNB<sup>[d-1]</sup>(θ, p). Note that $p = \frac{\alpha}{1+\alpha}$ and $\lambda = \theta\alpha$.</caption><thead><tr><th></th><th>d = 0</th><th>d > 0</th></tr></thead><tbody><tr><th>c &lt; 0</th><td>binomial<br>B<sub>&theta;</sub>(p)</td><td>generalized binomial<br>GB<sup>[d]</sup>(&theta;, p)</td></tr><tr><th>c = 0</th><td>Poisson<br>P(&lambda;)</td><td>generalized Poisson<br>GP<sup>[d]</sup>(&lambda;)</td></tr><tr><th>c &gt; 0</th><td>negative binomial<br>NB(&theta;, &alpha;)</td><td>generalized negative binomial<br>GNB<sup>[d]</sup>(&theta;, &alpha;)</td></tr></tbody></table>

It should be remarked that both distributions (definitions 1 and 2) are related by the same parametric sequence $a_{\theta}^{[c,d]}(n)$ that respects properties (P1) and (P2). This is the key point in the demonstration of the following Lemma.

**Lemma 3.** An additive modified power series distribution is additively closed under convolution, more precisely we have

$$APS^{[c,d]}(\theta, \alpha) * APS^{[c,d]}(\gamma, \alpha) = APS^{[c,d]}(\theta + \gamma, \alpha),$$

for all $(\theta, \gamma) \in \Theta^2$, $\alpha \in (0, R)$, $c \in \{-1, 0, 1\}$ and $d \in \mathbb{R}_+$ (the symbol $*$ here denotes the convolution of distributions).

The proof of this Lemma is given by Peyhardi (2023) in Proposition 1.

## 3 Quasi Pólya thinning operator

Let us first recall the definition of the classical binomial thinning operator. Let $X$ denote a count variable and $\rho \in (0, 1)$ a parameter, then $X' \stackrel{d}{=} \rho \circ X$ denote the thinned variable. It is usually defined as a random sum $X' = \sum_{i=1}^{X} Z_i$ where $Z_i \sim B(\rho)$ are independent random Bernoulli variables with parameter $\rho$. Two kinds of extension have been proposed: assume that $\rho$ is a random variable, or assume that $Z_i$ follow another count distribution than the Bernoulli.

To propose a new extension, we need to see the binomial thinning as a binomial mixed distribution. Remark that $X' \stackrel{d}{=} \rho \circ X$ is equivalent to $X' \sim B_n(\rho) \wedge \mathcal{L}$ where $\mathcal{L}$ denotes the distribution of $X$. It is thus easy to extend the binomial thinning to the quasi Pólya thinning as follows

$$X' \sim qP_n^{[c,d]}(\theta, \gamma) \wedge \mathcal{L}.$$

Note that with this parametrization, $\rho$ corresponds to the proportion $\theta/(\theta + \gamma)$. In the following, let $\rho \bullet$ denote the quasi Pólya thinning operator, i.e., $\rho \bullet X = X'$ with $X'|X = n \sim qP_n^{[c,d]}(\theta, \gamma)$ for all $n \in \mathbb{N}$. This new class of thinning operators includes five distinct examples; see Table 1 for notations.
Let us briefly recall the interpretation of (quasi) Pólya distribution in terms of urn models. The Pólya urn model is generally presented in terms of $n$ random drawings of balls from an urn, that initially contains $\theta$ white balls and $\gamma$ black balls. One ball is drawn at random and then replaced with $c$ additional balls of the same color. This procedure is repeated $n$ times and focus is made on the count $X'$ of drawn white balls. Knowing the number $n$ of draws, the count distribution for $X'$ is known as the Pólya distribution, i.e., we have $X'|X = n \sim \mathcal{P}_n^{[c]}(\theta, \gamma)$. In this framework, the parameter $\rho = \theta/(\theta + \gamma)$ corresponds to the proportion of white balls initially contained in the urn. The parameter $n$ being related to the support of the distribution, is intentionally noted as an index of the distribution, distinguishing it from other parameters $\theta$ and $\gamma$ used to define the pmf. The Pólya distribution turns out to be the hypergeometric distribution when $c = -1$ (without replacement), the binomial distribution when $c = 0$ (with replacement meaning independent draws) and the negative hypergeometric distribution when $c = 1$ (with reinforcement). Let us note that the last two cases are defined for continuous values of $\theta$ and $\gamma$. The negative hypergeometric distribution is thereby extended to the beta binomial distribution. Finally, the parameter $d$ is added to take into account a predetermined strategy, leading to the quasi Pólya distribution $q\mathcal{P}_n^{[c,d]}(\theta, \gamma)$; see Consul and Mittal (1975) and Janardan (1973) for more details. The quasi Pólya distribution turns out to the Pólya distribution when $d = 0$.

**Proposition 1.** Let $X$ denote an arbitrary random count variable, then

$$ \bullet E[\rho \cdot X] = \rho E[X] $$

• If $d = 0$, $\rho \cdot$ denotes the simple Pólya thinning operator $\mathcal{P}_n^{[c]}(\theta, \gamma) \wedge \mathcal{L}$, then

$$ V(\rho \cdot X) = \rho^2 V(X) + \rho(1-\rho) \frac{\theta + \gamma}{\theta + \gamma + c} E[X] + c\rho(1-\rho) \frac{1}{\theta + \gamma + c} E[X^2] $$

$$ \mathrm{Cov}(\rho \cdot X, X) = \rho V(X) $$

where $\rho = \theta/(\theta + \gamma)$.

*Proof.* Let $Y = \rho \cdot X$ and $Z = X - Y$. Next equalities are due to the total law of mean (variance and covariance) and moments of (quasi) Pólya given in Lemma 1.

$$ E[\rho \cdot X] = E[E[Y|X]] = E[\rho X] = \rho E[X] $$

Assuming $d = 0$ we have

$$ \begin{aligned} V(\rho \cdot X) &= V(E[Y|X]) + E[V(Y|X)] \\ &= V(\rho X) + E \left[ \frac{\theta \gamma}{(\theta + \gamma)(\theta + \gamma + c)} X + c \frac{\theta \gamma}{(\theta + \gamma)^2 (\theta + \gamma + c)} X^2 \right] \end{aligned} $$

$$ V(\rho \cdot X) = \rho^2 V(X) + \rho(1-\rho) \frac{\theta + \gamma}{\theta + \gamma + c} E[X] + c\rho(1-\rho) \frac{1}{\theta + \gamma + c} E[X^2] $$
and

$$
\begin{align*}
\mathrm{Cov}(\rho \bullet X, X) &= \mathrm{Cov}(\mathrm{E}[Y|X], \mathrm{E}[X|X]) + \mathrm{E}[\mathrm{Cov}(Y, X|X)] \\
&= \mathrm{Cov}(\rho X, X) + \mathrm{E}[\mathrm{Cov}(Y, Y + Z|X)] \\
&= \rho V(X) + \mathrm{E}[V(Y|X) + \mathrm{Cov}(Y, Z|X)] \\
\mathrm{Cov}(\rho \bullet X, X) &= \rho V(X)
\end{align*}
$$

remarking that $V(Y|X) = -\mathrm{Cov}(Y, Z|X)$ according to Lemma 1. $\square$

These results coincide with those of Weiß (2008b) about the hypergeometric thinning operator ($c = -1$ and $d = 0$), the binomial thinning operator ($c = 0$ and $d = 0$) and the beta binomial thinning operator ($c = 1$ and $d = 0$).

Now, let us describe different properties of quasi Pólya thinning operator and additive modified power series distributions, that will be used in next sections to demonstrate the stationarity of INAR(1) model based on quasi Pólya thinning operator.

**Lemma 4.** An additive modified power series distribution is closed under the corresponding quasi Pólya thinning operator if the parameters respect the sum constraint, more precisely we have

$$q\mathcal{P}_n^{[c,d]}(\theta, \gamma) \wedge_n \mathcal{APS}^{[c,d]}(\theta + \gamma, \alpha) = \mathcal{APS}^{[c,d]}(\theta, \alpha),$$

for all $(\theta, \gamma) \in \Theta^2$, $\alpha \in (0, R)$, $c \in \{-1, 0, 1\}$ and $d \in \mathbb{R}_+$.

The proof of this Lemma is given by Peyhardi (2023) in Theorem 1 (item 2).

**Lemma 5.** Let X, Y and Z denote three random count variables, then the following assertions are equivalent

$$\left\{ \begin{array}{l} X \sim \mathcal{APS}^{[c,d]}(\theta + \gamma, \alpha) \\ Y = \rho \bullet X \\ Z = X - Y \end{array} \right. \iff \left\{ \begin{array}{l} Y \sim \mathcal{APS}^{[c,d]}(\theta, \alpha) \\ Z \sim \mathcal{APS}^{[c,d]}(\gamma, \alpha) \\ Y \perp Z \\ X = Y + Z \end{array} \right.$$

*Proof.* • Assume that the left part holds. According to Lemma 4 we have $Y \sim \mathcal{APS}^{[c,d]}(\theta, \alpha)$. Since $Y|X = n \sim q\mathcal{P}_n^{[c,d]}(\theta, \gamma)$ we have $n - Y|X = n \sim q\mathcal{P}_n^{[c,d]}(\gamma, \theta)$ and thus $Z = X - Y \sim q\mathcal{P}_n^{[c,d]}(\gamma, \theta) \wedge_n \mathcal{APS}^{[c,d]}(\theta + \gamma, \alpha)$. Using Lemma 4 we obtain $Z \sim \mathcal{APS}^{[c,d]}(\gamma, \alpha)$. Finally, the joint distribution of Y and Z is given by

$$P(Y = y, Z = z) = P(Y = y, Z = z | X = y + z) P(X = y + z) \\ = \frac{a_{\theta}^{[c,d]}(y) a_{\gamma}(z) a_{\theta+\gamma}^{[c,d]}(y+z) \{g(\alpha)\}^{y+z}}{a_{\theta+\gamma}^{[c,d]}(y+z) h_{\theta+\gamma}(\alpha)}$$

$$P(Y = y, Z = z) = \frac{a_{\theta}^{[c,d]}(y) \{g(\alpha)\}^{y} a_{\gamma}^{[c,d]}(z) \{g(\alpha)\}^{z}}{h_{\theta+\gamma}(\alpha)},$$
We know that $h_{\theta+\gamma} = h_{\theta}h_{\gamma}$ and thus $Y$ and $Z$ are independent.

* Assume that the right part holds. According to Lemma 3 we have $X \sim \mathcal{APS}^{[c,d]}(\theta+\gamma, \alpha)$. Moreover, we have

$$ 
\begin{aligned}
P(Y=y|X=n) &= \frac{P(Y=y, X=n)}{P(X=n)} \\
&= \frac{P(Y=y)P(Z=n-y)}{P(X=n)} \\
&= \frac{a_{\theta}^{[c,d]}(y)\{g(\alpha)\}^y a_{\gamma}^{[c,d]}(n-y)\{g(\alpha)\}^{n-y}}{h_{\theta}(\alpha) h_{\gamma}(\alpha) a_{\theta+\gamma}^{[c,d]}(n)\{g(\alpha)\}^n} h_{\theta+\gamma}(\alpha) \\
P(Y=y|X=n) &= \frac{a_{\theta}^{[c,d]}(y)a_{\gamma}^{[c,d]}(n-y)}{a_{\theta+\gamma}^{[c,d]}(n)},
\end{aligned}
$$

i.e., $Y|X=n \sim q\mathcal{P}_n^{[c,d]}(\theta, \gamma)$ or equivalently $Y = \rho \bullet X$.

□

**Lemma 6.** If $X \sim \mathcal{APS}^{[c,d]}(\theta+\gamma, \alpha)$ then

$$ 
\begin{aligned}
E[\rho \bullet X] &= \rho E[X], \\
V(\rho \bullet X) &= \rho V(X), \\
Cov(\rho \bullet X, X) &= \rho V(X),
\end{aligned}
$$

where $\rho = \theta/(\theta + \gamma)$.

*Proof.* According to Lemma 4 we have $\rho \bullet X \sim \mathcal{APS}^{[c,d]}(\theta, \alpha)$. Then according to Lemma 2 we have

$$ 
\begin{aligned}
E[\rho \bullet X] &= \theta f_1(\alpha) \\
&= \frac{\theta}{\theta + \gamma} (\theta + \gamma) f_1(\alpha) \\
E[\rho \bullet X] &= \rho E(X),
\end{aligned}
$$

where $\rho = \frac{\theta}{\theta+\gamma}$. In the same way, we obtain $V(\rho \bullet X) = \rho V(X)$. Using Lemma 5 we have

$$ 
\begin{aligned}
\mathrm{Cov}(\rho \bullet X, X) &= \mathrm{Cov}(Y, Y + Z) \\
&= V(Y) + \underbrace{\mathrm{Cov}(Y, Z)}_{=0} \\
\mathrm{Cov}(\rho \bullet X, X) &= \rho V(X).
\end{aligned}
$$

□
4 Integer valued autoregressive models of order 1

Let $(X_t)_{t \in \mathbb{N}}$ denote a discrete time process with discrete state space. An INAR(1) model based on quasi Pólya thinning operator is defined as follows:

$$X_t = \rho \cdot X_{t-1} + \varepsilon_t$$

where $\varepsilon_t \sim \mathcal{L}_\varepsilon$ are independent and identically distributed random variables and $\rho \cdot$ denotes a quasi Pólya thinning operator (with parameters $\psi = (c, d, \theta, \gamma)$) independent of $(\varepsilon_t)_{t \in \mathbb{N}}$. Denoting by $\mathcal{L}_t$ the distribution of $X_t$ we have

$$\mathcal{L}_t = \{ q\mathcal{P}_n^{[c,d]}(\theta, \gamma) \wedge \mathcal{L}_{t-1} \} * \mathcal{L}_\varepsilon,$$

where * denotes the convolution of distributions.

**Theorem 2.** Let $(X_t)_{t \in \mathbb{N}}$ follow an INAR(1) model based on quasi Pólya thinning operator. Assume that the error distribution is the specific additive modified power series distribution $\mathcal{L}_\varepsilon = \mathcal{APS}^{[c,d]}(\gamma, \alpha)$ for some $\alpha \in (0, R)$. Then $(X_t)_{t \in \mathbb{N}}$ is an ergodic Markov chain with stationary distribution $\mathcal{APS}^{[c,d]}(\theta + \gamma, \alpha)$. Moreover, at stationary state we have

$$E[X_t] = \frac{\mu_\varepsilon}{1 - \rho},$$

$$V(X_t) = \frac{\sigma_\varepsilon^2}{1 - \rho},$$

$$\mathrm{Corr}(X_{t+h}, X_t) = \rho^h,$$

where $\mu_\varepsilon = E(\varepsilon_t)$, $\sigma_\varepsilon^2 = V(\varepsilon_t)$ and $\rho = \theta/(\theta + \gamma)$. Finally, the process $(X_t)_{t \in \mathbb{N}}$ is time reversible.

*Proof.* The conditional distribution of $X_{t+1}$ given $X_t = x_t$ is

$$p(x_{t+1}|x_t) = \sum_{k=0}^{\min(x_{t+1},x_t)} \frac{a_{\theta}^{[c,d]}(k) a_{\gamma}^{[c,d]}(x_t - k) a_{\gamma}^{[c,d]}(x_{t+1} - k) \{g(\alpha)\}^{x_{t+1}-k}}{a_{\theta+\gamma}^{[c,d]}(x_t) h_{\gamma}(\alpha)} \quad (3)$$

The support of $p(\cdot|x)$ is related to the support of $a_{\theta}^{[c,d]}(\cdot)$. A necessary condition for the ergodicity is to have $p(x+1|x) > 0$ and $p(x|x+1) > 0$ for all $x$ such that $x$ and $x+1$ belong to the support of the distribution $\mathcal{APS}^{[c,d]}(\theta + \gamma, \alpha)$. For the four cases such that $(c, d) \neq (-1, 0)$, the demonstration is straightforward since $\mathrm{Supp} \ a_{\theta}^{[c,d]} = \mathbb{N}$. For the special case $(c, d) = (-1, 0)$, we have $\mathrm{Supp} \ a_{\theta}^{[-1,0]} = \{0, \dots, \theta\}$. In this case, let show that $p(x+1|x) > 0$ for all $x \in \{0, \dots, \theta+\gamma-1\}$. It is sufficient to show that there exists $k \le x$ such that $k \le \theta$ and $x+1-k \le \gamma$ since at least one term in the sum of equation (3) would be positive. Remark that $k = \min(x, \theta)$ is a good candidate. A similar argument leads to $p(x|x+1) > 0$ for all $x \in \{0, \dots, \theta+\gamma-1\}$. The Markov chain $(X_t)_{t \in \mathbb{N}}$ is thus ergodic.
Now let us show that the stationary distribution is $\mathcal{AP}S^{[c,d]}(\theta+\gamma, \alpha)$. Assume that $X_{t-1} \sim \mathcal{AP}S^{[c,d]}(\theta+\gamma, \alpha)$ and remark that the distribution of $X_t$ is

$$
\begin{align*}
\mathcal{L}_t &= \{q\mathcal{P}_n^{[c,d]}(\theta, \gamma) \wedge \mathcal{L}_{t-1}\} * \mathcal{L}_\epsilon, \\
&= \{q\mathcal{P}_n^{[c,d]}(\theta, \gamma) \wedge \mathcal{AP}S^{[c,d]}(\theta + \gamma, \alpha)\} * \mathcal{AP}S^{[c,d]}(\gamma, \alpha) \\
&= \mathcal{AP}S^{[c,d]}(\theta, \alpha) * \mathcal{AP}S^{[c,d]}(\gamma, \alpha) \\
\mathcal{L}_t &= \mathcal{AP}S^{[c,d]}(\theta + \gamma, \alpha)
\end{align*}
$$

where the third and last equalities are obtained according to lemma 4 and Lemma 3 respectively.

Using Lemma 2 it is easily seen that $E[X_t] = \frac{\mu_\epsilon}{1-\rho}$ and $V(X_t) = \frac{\sigma_\epsilon^2}{1-\rho}$. Otherwise, the covariance is given by

$$
\begin{align*}
\mathrm{Cov}(X_{t+1}, X_t) &= \mathrm{Cov}(\rho \bullet X_t + \varepsilon_{t+1}, X_t) \\
&= \mathrm{Cov}(\rho \bullet X_t, X_t) + \underbrace{\mathrm{Cov}(\varepsilon_{t+1}, X_t)}_{=0} \\
\mathrm{Cov}(X_{t+1}, X_t) &= \rho V(X_t),
\end{align*}
$$

according to Lemma 6 and thus Corr($X_{t+1}, X_t$) = $\rho$. By recursion, we obtain Corr($X_{t+h}, X_t$) = $\rho^h$.

Finally, to demonstrate the time reversibility of the process $(X_t)_{t \in \mathbb{N}}$, remark that the joint distribution is symmetric

$$
\begin{align*}
p(x_{t+1}, x_t) &= p(x_t)p(x_{t+1}|x_t) \\
&= \frac{a_{\theta+\gamma}^{[c,d]}(x_t)\{g(\alpha)\}^{x_t}}{h_{\theta+\gamma}(\alpha)} \sum_{k=0}^{\min(x_{t+1}, x_t)} \frac{a_{\theta}^{[c,d]}(k)a_{\gamma}^{[c,d]}(x_t-k)a_{\gamma}^{[c,d]}(x_{t+1}-k)\{g(\alpha)\}^{x_{t+1}-k}}{h_{\gamma}(\alpha)} \\
p(x_{t+1}, x_t) &= \frac{\{g(\alpha)\}^{x_{t+1}+x_t}}{h_{\theta+2\gamma}(\alpha)} \sum_{k=0}^{\min(x_{t+1}, x_t)} a_{\theta}^{[c,d]}(k)a_{\gamma}^{[c,d]}(x_t-k)a_{\gamma}^{[c,d]}(x_{t+1}-k)\{g(\alpha)\}^{-k}.
\end{align*}
$$

□

## Models comparison

We have seen, in the previous section, that there are five distinct quasi-Pólya thinning operators. It implies that there are five distinct classes of INAR(1) models. They can be distinguished by their marginal distribution. Note that the number of parameters is varying between two and four according to these five cases:

1. Binomial marginal distribution: this is the case $c = -1$ and $d = 0$ that corresponds to the binomial INAR(1) model introduced by Al-Osh and Alzaid (1991). The vector of parameters is $\psi = (\theta, \gamma, \alpha)$ and the space of parameters is $\Psi = \mathbb{N}^{*2} \times \mathbb{R}_+^*$. This model will be denoted by B-INAR(1).
<table><caption>**Table 3** Unity variance function of the five additive modified power series</caption><thead><tr><th>Distribution</th><th>V(v)</th></tr></thead><tbody><tr><td>Binomial</td><td>v(1 - v)</td></tr><tr><td>Poisson</td><td>v</td></tr><tr><td>Negative binomial</td><td>v(1 + v)</td></tr><tr><td>Generalized Poisson</td><td>v(1 + v)<sup>2</sup></td></tr><tr><td>Generalized negative binomial</td><td>v(1 + v)(1 + v<sup>d+1</sup>/d)</td></tr></tbody></table>

2. **Poisson marginal distribution:** this is the case $c = 0$ and $d = 0$ that corresponds to the usual Poisson INAR(1) model introduced by McKenzie (1985); Al-Osh and Alzaid (1987). The vector of parameter is $\psi = (\rho, \alpha)$, where $\rho = \theta/(\theta + \gamma)$, and the space of parameters is $\Psi = (0, 1) \times \mathbb{R}_+^*$. This model will be denoted by P-INAR(1).

3. **Negative binomial marginal distribution:** this is the case $c = 1$ and $d = 0$ that corresponds to the negative binomial INAR(1) model discussed by Joe (1996) as a special case of convolution-closed infinitely divisible class. The vector of parameter is $\psi = (\theta, \gamma, \alpha)$ and the space of parameters is $\Psi = \mathbb{R}_+^{*2} \times (0, 1)$. This model will be denoted by NB-INAR(1).

4. **Generalized Poisson marginal distribution:** this is the case $c = 0$ and $d > 0$ that corresponds to the generalized Poisson INAR(1) model introduced by Alzaid and Al-Osh (1993). The vector of parameter is $\psi = (\rho, \alpha, d)$, where $\rho = \theta/(\theta + \gamma)$, and the space of parameters is $\Psi = (0, 1) \times \mathbb{R}_+^* \times [0, 1]$. Note that for the generalized Poisson distribution $GP^{[d]}(\lambda)$, the parameter $d$ must be included in interval $[0, 1]$ to obtain a proper distribution (Consul and Famoye, 2006). This INAR(1) model will be denoted by GP-INAR(1). It includes the P-INAR(1) model as a special case when $d = 0$.

5. **Generalized negative binomial marginal distribution:** this is the case $c = 1$ and $d > -1$ that corresponds to the generalized negative binomial INAR(1) model. This is a new INAR(1) model. The vector of parameter is $\psi = (\theta, \gamma, \alpha, d)$ and the space of parameters is $\Psi = \mathbb{R}_+^{*2} \times (0, 1) \times [0, (1 - \alpha)/\alpha]$. Note that for the generalized negative binomial distribution $GNB^{[d]}(\theta, \alpha)$, the parameter $d$ must be included in interval $[0, (1 - \alpha)/\alpha]$ to obtain a proper distribution (Consul and Famoye, 2006). This INAR(1) model will be denoted by GNB-INAR(1). It includes the B-INAR(1) model as a special case when $d = -1$ and the NB-INAR(1) model when $d = 0$.

Since additive modified power series distributions belong to the family of dispersion models (Jorgensen, 1997), their unity variance functions have closed form; see Table 3. As remarked by Jourdan and Kokonendji (2002) these unity variance functions are strictly ordered as follows,

$$V_B < V_P < V_{NB} < V_{GP} < V_{GNB}.$$

Therefore, according to the level of dispersion observed in the dataset, the autoregressive model with the appropriate marginal distribution has to be used. The different models could be compared using penalized criteria, like Akaike information criterion (AIC) as proposed by Weiß (2008b).
# 5 Model inference based on likelihood

## Maximum likelihood estimator

The log-likelihood of a count time series $x = (x_0, x_1, \dots, x_T)$ for an INAR(1) model based on quasi Pólya thinning operator with parameter $\psi = (\theta, \gamma, \alpha, d)$ is

$$L_T(\psi) = l(x_0; \psi) + \sum_{t=1}^{T} l(x_t|x_{t-1}; \psi), \quad (4)$$

where $l(x_0; \psi) = \ln p(x_0; \psi)$, $l(x_t|x_{t-1}; \psi) = \ln p(x_t|x_{t-1}; \psi)$. Since $X_0$ is assumed to follow the modified power series distribution $\mathcal{APS}^{[c,d]}(\theta + \gamma, \alpha)$, then its pmf is

$$p(x_0; \psi) = \frac{a_{\theta+\gamma}^{[c,d]}(x_0)\{g(\alpha)\}^{x_0}}{h_{\theta+\gamma}(\alpha)}.$$

Note that the log-likelihood of the first observation could be omitted, since it is asymptotically irrelevant. The right part of the log-likelihood is obtained from (3).

Let $r$ denote the dimension of the vector of parameters $\psi$ (note that $r$ is varying between 2 and 4 according to the five distinct classes of INAR(1) model). Let $l_u = \partial l/\partial \psi_u$ denote the partial derivative with respect to the $u^{th}$ component of $\psi$ for $u = 1, \dots, r$. In the same way, let $l_{uv}$ and $l_{uvw}$ denote the partial derivatives of second and third order (same notations hold for $p$). Finally, $E_\psi$ means expectation, assuming that the true parameter value is $\psi$ and $X_0$ starts with a stationary distribution. Following Joe (1997) (page 317) we introduce the following regularity conditions:

(C1) For $\psi \in \Psi$ there exists a neighborhood $N_\psi$ of $\psi$ such that for all $u, v, w$ and $x$

$$\begin{align*}
\sum_{y \in \mathbb{N}} \sup_{\psi' \in N_{\psi}} p_u(y|x; \psi) &< \infty \\
\sum_{y \in \mathbb{N}} \sup_{\psi' \in N_{\psi}} p_{uv}(y|x; \psi) &< \infty \\
E_{\psi} \left[ \sup_{\psi' \in N_{\psi}} |l_{uvw}(X_1|X_0; \psi)| \right] &< \infty
\end{align*}$$

(C2) For $u = 1, \dots, r$ $E_\psi[|l_u(X_1|X_0; \psi)|^2] < \infty$ and $\Sigma(\psi) = (\sigma_{uv}(\psi))$ is a non-singular $r \times r$ matrix where $\sigma_{uv}(\psi) = E_\psi[l_u(X_1|X_0; \psi)l_v(X_1|X_0; \psi)]$.

**Theorem 3.** Let $(X_t)_{t \in \mathbb{N}}$ follow an INAR(1) model based on quasi Pólya thinning operator. Assume that $d \neq 0$ or $c \ge 0$, avoiding the case of INAR(1) model based on hypergeometric thinning operator. Assume that the starting distribution is $\mathcal{L}_0 = \mathcal{APS}^{[c,d]}(\theta + \gamma, \alpha)$ for some $\alpha \in (0, R)$ and the error distribution is $\mathcal{L}_\epsilon = \mathcal{APS}^{[c,d]}(\gamma, \alpha)$. Under the regularity conditions (C1) and (C2), the MLE $\hat{\psi}_T$ converges
in probability to the true parameter $\psi$ with asymptotic normal distribution

$$ \frac{\hat{\psi}_T - \psi}{\sqrt{T}} \underset{T \to \infty}{\longrightarrow} \mathcal{N}(0, \Sigma^{-1}(\psi)). $$

*Proof.* The proof follows the approach of Joe (1997), that is based on the result of Billingsley (1961) about Markov chains of first order. Let us verify the following regularity conditions:

* All states of the Markov chain communicate with each other. In the proof of previous theorem, it has already been shown that there are no transient state since $p(x+1|x) > 0$ and $p(x|x+1) > 0$ for all $x \in \mathbb{N}$.

* The support of $p(\cdot|x; \psi)$ does not depend on $\psi$. According to (3) it is easily seen that $\text{Supp } a_\theta = \mathbb{N} \Rightarrow \text{Supp } p(\cdot|x; \psi) = \mathbb{N}$. Note that this condition is not verified for the B-INAR(1) model since $\text{Supp } a_\theta = \{0, 1, \dots, \theta\}$ in this case.

* $p_u, p_{uv}, p_{uvw}$ ($u, v, w \in \{1, \dots, r\}$) exist and are continuous in $\psi$. According to (3), $p(\cdot|x; \psi)$ is a sum of ratio of functions $a_\theta, a_\gamma, g(\alpha)$ and $h_\gamma(\alpha)$. Since the gamma, power and exponential functions are smooth, this condition is obtained.

* $p(\cdot|x; \psi)$ is absolutely continuous with respect to the stationary measure $\mu(\cdot; \psi)$. In our case, the stationary measure of the INAR(1) model is the modified power series distribution $\mathcal{APS}^{[c,d]}(\theta + \gamma, \alpha)$ for some $\alpha \in (0, R)$. Except for the binomial case ($c < 0$ and $d = 0$), the support of this distribution is the entire set of non-negative integers $\mathbb{N}$.

□

As remarked by Joe (1997) this asymptotic result implies, as in the i.i.d. case, that the negative inverse Hessian of $L_T(\psi)$ evaluated at the MLE $\hat{\psi}$ can be used as an estimated covariance of $\hat{\psi}$. In Section Application, the MLE is performed via the optim function contained in the stats R package. The optim argument, hessian=TRUE, allows determining the standard errors of MLE (via the Hessian pseudo-inverse when the Hessian is non-invertible).

## Likelihood ratio test

A likelihood ratio test can also be employed, since the stationary distributions of the five different INAR(1) models are nested; see Table 2 (they belong to the family of additive modified power series distributions). The case $c < 0$ is not totally separated from the case $c > 0$ since both distribution $\mathcal{APS}^{[-1,d]}$ and $\mathcal{APS}^{[1,d-1]}$ are equivalent. Therefore, the case $c < 0$ is omitted and focus is made on the four distributions: Poisson, negative binomial, generalized Poisson, generalized negative binomial.

* The NB-INAR(1) model can be compared with the P-INAR(1) model through the null hypothesis $H_0: c = 0$ against $H_1: c > 0$ (assuming $d = 0$). The log-likelihood test statistic $2(\mathcal{L}_{NB} - \mathcal{L}_P)$ follows a chi-square distribution $\chi^2(1)$ with one degree of freedom, under the null hypothesis.
* The GP-INAR(1) model can be compared with the P-INAR(1) through the null hypothesis $H'_0: d = 0$ against $H'_1: d > 0$ (assuming $c = 0$). The log-likelihood test statistic $2(\mathcal{L}_{GP} - \mathcal{L}_P)$ follows a chi-square distribution $\chi^2(1)$ with one degree of freedom, under the null hypothesis.

* The GNB-INAR(1) model can be compared with the P-INAR(1) through the null hypothesis $H''_0: (c, d) = (0, 0)$ against $H''_1: (c, d) \neq (0, 0)$. The log-likelihood test statistic $2(\mathcal{L}_{GNB} - \mathcal{L}_P)$ follows a chi-square distribution $\chi^2(2)$ with two degree of freedom, under the null hypothesis.

## Model validation

The standardized Pearson residuals are given by

$$e_t := \frac{E[X_t|x_{t-1}] - x_t}{\sqrt{V(X_t|x_{t-1})}}.$$

For an adequate model, these residuals should have a null mean and a unit variance, according to Harvey and Fernandes (1989). In our context of INAR(1) model based on quasi Pólya thinning operator, if we denote $Y_t := \rho \cdot X_{t-1}$, we have $E[X_t|x_{t-1}] = E[Y_t|x_{t-1}] + E[\varepsilon_t]$ and $V(X_t|x_{t-1}) = V(Y_t|x_{t-1}) + V(\varepsilon_t)$ where $Y_t|x_{t-1} \sim q\mathcal{P}_{x_{t-1}}^{[c,d]}(\theta, \gamma)$ and $\varepsilon_t \sim \mathcal{APS}^{[c,d]}(\gamma, \alpha)$ are independent variables. According to the result of Janardan and Schaeffer (1977) (given in Lemma 1) we have

$$E[Y_t|x_{t-1}] = \frac{\theta}{\theta + \gamma} x_{t-1}, \\ V(Y_t|x_{t-1}) = \frac{\theta\gamma}{\theta + \gamma} \left[ \frac{x_{t-1}^2}{\theta + \gamma} - \sum_{j=1}^{x_{t-1}-1} \frac{x_{t-1}^{(j+1;-1)}(d+c)^{j-1}}{\{\theta + \gamma + dx_{t-1} + (x_{t-1}-j+2)c\}^{(j;c)}} \right].$$

Moreover, according to Gupta (1974) we have

$$\mu = E[\varepsilon_t] = \frac{g(\alpha) h'_\gamma(\alpha)}{g'(\alpha) h_\gamma(\alpha)} \quad \text{and} \quad V(\varepsilon_t) = \frac{g(\alpha)}{g'(\alpha)} \mu'(\alpha).$$

Appendix B details the computation for the five different INAR(1) models.

# 6 Application

In this section we will compare the five INAR(1) models based on quasi Pólya thinning operator, described in Section 4: the B-INAR(1), P-INAR(1), NB-INAR(1), GP-INAR(1) and the GNB-INAR(1) models. The MLEs have been computed using the optim_sa function, available in the R-package optimization.

## 6.1 Simulated data

The five kinds of INAR(1) models have been simulated with a length of time series $T = 50$ and auto correlation $\rho = 0.3, \rho = 0.6$ and $\rho = 0.9$. For each of the five models,
three values of parameters $\psi$ have been simulated, $n = 1000$ times, in order to obtain different mean values: $\mu = 5$, $\mu = 20$ and $\mu = 50$. The variance has also been calibrated from under dispersion $\sigma^2 = \mu/2$ with the B-INAR(1) model to huge over dispersion $\sigma^2 = 8\mu$ with the GNB-INAR(1) model. Finally, $5 \times 3 \times 3 = 45$ scenarios have been simulated $n = 1000$ times.

For instance, focusing on the B-INAR(1) model, the parameters are given by $\theta = 2\mu\rho$, $\gamma = 2\mu(1 - \rho)$ and $\alpha = 1$. Following Theorem 2 the first value is simulated according to the binomial distribution $X_0 \sim B_{\theta+\gamma}(p)$ where $p = \alpha/(1 + \alpha)$ and thus the process is stationary with $X_t \sim B_{\theta+\gamma}(p)$ for all $t = 1, \dots, T$. Therefore, we have $E[X_t] = (\theta + \gamma)p = \mu$ and $V(X_t) = (\theta + \gamma)p(1 - p) = \mu/2$. For instance, with $\mu = 5$ we have $\psi = (\theta, \gamma, \alpha) = (6, 4, 1)$. The parameters $\psi$ of other INAR(1) models have been determined following the same way.

Then for each of the 45 scenarios, the five INAR(1) models have been inferred using the MLE on the $n = 1000$ simulated time series. The estimated marginal distributions have been compared with the true marginal distribution using the Kullback-Leibler (KL) divergence. It was computed as $\sum_i p(i) \ln \{p(i)/\hat{p}(i)\}$ for $i$ between minimum and maximum observed values, where $p$ (resp. $\hat{p}$) denotes the true (resp. estimated) pmf. The results have been summarized in Table 4 (resp. Tables C3 and C4 in Appendix C) for time series of length $T = 50$ with autocorrelation $\rho = 0.6$ (resp. $\rho = 0.3$ and $\rho = 0.9$). For each scenario, the mean (and standard deviations) of the $n = 1000$ KL divergences are given, and the smaller (among the five models) appears in bold.

The results can be interpreted according to the three kinds of dispersion: under dispersion (binomial), equi dispersion (Poisson) and over dispersion (negative binomial, generalized Poisson and generalized negative binomial). When the time series are simulated according to a B-INAR(1) model with small or middle auto correlation ($\rho = 0.3$ or $0.6$), the smallest KL divergence is obtained with the B-INAR(1) model. The other models are more divergent except the GNB-INAR(1) model when the auto correlation is small ($\rho = 0.3$). When the auto correlation is high ($\rho = 0.9$), the smallest KL divergence is obtained with the P-INAR(1) model but the results are very close between the five models. When the time series are simulated according to a P-INAR(1) model, the KL divergence is smaller with the P-INAR(1) model but is approximately the same with the GP-INAR(1) models. When the time series are simulated according to an over dispersed INAR(1) model (with NB, GP or GNB margins), the KL divergence is huge for the B-INAR(1) model or P-INAR(1) model. On the contrary, the KL divergence is small for the three other INAR(1) models, particularly with the GP and GNB-INAR(1) models. More precisely, when the time series are simulated according to a NB-INAR(1) model, the KL divergence is smaller with the GP-INAR(1) model or GNB-INAR(1) model but never with the NB-INAR(1) model. Performances of GP-INAR(1) model and GNB-INAR(1) model are very similar. Nevertheless, it should be kept in mind that the GNB-INAR(1) model has four parameters, whereas its competitor only has three. Penalized criteria could thus be used when comparing the different INAR(1) models on real data. The results are globally similar for the different mean values ($\mu = 5$, $\mu = 20$ or $\mu = 50$) except when the mean and the variance are high together ($\mu = 50$ and $\sigma^2 = 8\mu$). In this case, the GP-INAR(1) model gives very poor results and is outperformed by both NB and GNB-INAR(1) models that are stable.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Summary of soap dataset about sales counts. Plot of the count time series in top row and plots of sample probability mass function and (partial) autocorrelation function in bottom.</figcaption></figure>

Finally, the results are impact by the auto correlation values. The KL divergences are globally smaller for small auto correlation ($\rho$ = 0.3) and higher for high auto correlation ($\rho$ = 0.9).

## 6.2 Real data

As an example, let us use the dataset of soap product sales presented by MacDonald and Zucchini (2016). This is a count time series of length $T = 242$ about the weekly number of sales of a soap product in a supermarket. The counts exhibit an important over dispersion, which is also remarkable from the plot of the sample pmf in Figure 1. The partial ACF (PACF) has a significant value only at lag 1 ($\hat{\rho}(1) \approx 0.39$). It is thus natural to use INAR(1) models on this dataset. The sample mean is equal to 5.44 and the variance is equal to 15.40.

The five INAR(1) models based on quasi Pólya thinning operator have been inferred with the MLE. Summarized results are given in Table 5, allowing the comparison of these models. Using the AIC for model selection, the GP-INAR(1) model is preferred; see Figure 2 for the adjusted pmf and ACFs using this model. Note that
<table><caption>**Table 4** Means (and standard deviations) of KL divergences for different INAR(1) models based on quasi Pólya thinning operator for integer valued time series of length *T* = 50 and autocorrelation *ρ* = 0.6. The rows indicate the marginal distribution of the simulated INAR(1) model whereas the columns indicate the marginal distribution of the estimated INAR(1) model. Means of KL divergences (and standard deviations) have been computed for the 5 different estimated models, on *n* = 1000 simulations for each of the 15 scenarios. The minimal value of each scenario is written in bold.</caption><thead><tr><th>Marginal</th><th>Mean</th><th>Variance</th><th>B</th><th>P</th><th>NB</th><th>GP</th><th>GNB</th></tr></thead><tbody><tr><td rowspan="3">B</td><td>μ = 5</td><td></td><td><b>7.73e-02</b><br>(9.28e-02)</td><td>1.15e-01<br>(3.10e-02)</td><td>1.22e-01<br>(3.66e-02)</td><td>1.27e-01<br>(4.09e-02)</td><td>1.48e-01<br>(6.16e-02)</td></tr><tr><td>μ = 20</td><td>σ² = μ/2</td><td><b>7.14e-02</b><br>(7.26e-02)</td><td>1.11e-01<br>(3.20e-02)</td><td>1.23e-01<br>(4.22e-02)</td><td>1.28e-01<br>(4.48e-02)</td><td>1.63e-01<br>(6.49e-02)</td></tr><tr><td>μ = 50</td><td></td><td><b>7.37e-02</b><br>(1.13e-01)</td><td>1.12e-01<br>(3.41e-02)</td><td>1.23e-01<br>(3.91e-02)</td><td>1.36e-01<br>(6.07e-02)</td><td>1.93e-01<br>(9.72e-02)</td></tr><tr><td rowspan="3">P</td><td>μ = 5</td><td></td><td>6.11e-02<br>(8.06e-02)</td><td><b>3.87e-02</b><br>(6.04e-02)</td><td>4.20e-02<br>(6.15e-02)</td><td>4.14e-02<br>(5.84e-02)</td><td>5.80e-02<br>(8.04e-02)</td></tr><tr><td>μ = 20</td><td>σ² = μ</td><td>6.89e-02<br>(1.71e-01)</td><td><b>3.56e-02</b><br>(5.01e-02)</td><td>3.93e-02<br>(5.43e-02)</td><td>3.92e-02<br>(5.07e-02)</td><td>5.58e-02<br>(7.04e-02)</td></tr><tr><td>μ = 50</td><td></td><td>7.06e-02<br>(2.73e-01)</td><td><b>3.64e-02</b><br>(4.86e-02)</td><td>3.99e-02<br>(5.19e-02)</td><td>4.08e-02<br>(5.18e-02)</td><td>5.10e-02<br>(6.02e-02)</td></tr><tr><td rowspan="3">NB</td><td>μ = 5</td><td></td><td>2.35e-01<br>(1.27e-01)</td><td>2.19e-01<br>(1.18e-01)</td><td>7.06e-02<br>(8.09e-02)</td><td><b>6.36e-02</b><br>(7.09e-02)</td><td>6.56e-02<br>(7.37e-02)</td></tr><tr><td>μ = 20</td><td>σ² = 2μ</td><td>2.52e-01<br>(1.35e-01)</td><td>2.29e-01<br>(1.27e-01)</td><td>7.05e-02<br>(7.80e-02)</td><td><b>6.52e-02</b><br>(7.26e-02)</td><td>6.55e-02<br>(7.31e-02)</td></tr><tr><td>μ = 50</td><td></td><td>2.66e-01<br>(1.48e-01)</td><td>2.34e-01<br>(1.26e-01)</td><td>7.01e-02<br>(7.97e-02)</td><td><b>6.83e-02</b><br>(7.41e-02)</td><td>6.87e-02<br>(7.49e-02)</td></tr><tr><td rowspan="3">GP</td><td>μ = 5</td><td></td><td>7.28e-01<br>(2.22e-01)</td><td>6.99e-01<br>(2.18e-01)</td><td>7.85e-02<br>(8.88e-02)</td><td><b>6.51e-02</b><br>(7.34e-02)</td><td>6.73e-02<br>(7.84e-02)</td></tr><tr><td>μ = 20</td><td>σ² = 4μ</td><td>8.29e-01<br>(2.13e-01)</td><td>7.91e-01<br>(2.07e-01)</td><td>7.12e-02<br>(8.95e-02)</td><td>6.45e-02<br>(7.95e-02)</td><td><b>6.37e-02</b><br>(7.82e-02)</td></tr><tr><td>μ = 50</td><td></td><td>8.82e-01<br>(2.06e-01)</td><td>8.32e-01<br>(2.02e-01)</td><td>7.75e-02<br>(8.88e-02)</td><td>6.93e-02<br>(7.69e-02)</td><td><b>6.91e-02</b><br>(7.62e-02)</td></tr><tr><td rowspan="3">GNB</td><td>μ = 5</td><td></td><td>1.67e+00<br>(5.76e-01)</td><td>1.62e+00<br>(5.18e-01)</td><td>8.87e-02<br>(1.09e-01)</td><td><b>6.56e-02</b><br>(8.64e-02)</td><td>6.99e-02<br>(9.92e-02)</td></tr><tr><td>μ = 20</td><td>σ² = 8μ</td><td>2.14e+00<br>(4.52e-01)</td><td>2.07e+00<br>(4.36e-01)</td><td>7.76e-02<br>(8.44e-02)</td><td><b>6.28e-02</b><br>(7.67e-02)</td><td>6.45e-02<br>(7.48e-02)</td></tr><tr><td>μ = 50</td><td></td><td>2.39e+00<br>(5.24e-01)</td><td>9.97e+00<br>(2.94e+01)</td><td>7.89e-02<br>(8.88e-02)</td><td>7.99e+00<br>(2.99e+01)</td><td><b>6.98e-02</b><br>(7.92e-02)</td></tr></tbody></table>
<table><caption>**Table 5** Summarized results for the five INAR(1) models based on quasi Pólya thinning operator, inferred on the soap products dataset. For each model, the log-likelihood, the AIC are given. Then the estimations of the mean, variance and autocorrelation (lag from 1 to 4) of the stationary distribution are given. Finally the mean and variance of the estimated standardized Pearson residuals are given. For each column, the best value is written in bold.</caption><thead><tr><th rowspan="2">Sample</th><th rowspan="2">L</th><th rowspan="2">AIC</th><th rowspan="2">Mean</th><th rowspan="2">Var</th><th colspan="4">ACF at lag</th><th colspan="2">Pearson residuals</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>Mean</th><th>Var</th></tr></thead><tbody><tr><td>B-INAR(1)</td><td>-682.73</td><td>1371.46</td><td>5.693</td><td>5.636</td><td>0.2604</td><td>0.0678</td><td>0.0177</td><td>0.0046</td><td>0.08653</td><td>2.448</td></tr><tr><td>P-INAR(1)</td><td>-680.12</td><td>1364.25</td><td>5.340</td><td>5.340</td><td>0.2300</td><td>0.0529</td><td>0.0122</td><td>0.0028</td><td>-0.02728</td><td>2.528</td></tr><tr><td>NB-INAR(1)</td><td>-614.69</td><td>1235.39</td><td><strong>5.408</strong></td><td>12.578</td><td>0.3235</td><td>0.1047</td><td>0.0339</td><td>0.0110</td><td><strong>-0.00257</strong></td><td>1.107</td></tr><tr><td>GP-INAR(1)</td><td>-614.55</td><td><strong>1235.10</strong></td><td>5.328</td><td>13.008</td><td><strong>0.3700</strong></td><td><strong>0.1369</strong></td><td><strong>0.0507</strong></td><td><strong>0.0187</strong></td><td>-0.02426</td><td>1.120</td></tr><tr><td>GNB-INAR(1)</td><td><strong>-614.26</strong></td><td>1236.51</td><td>5.542</td><td><strong>14.856</strong></td><td>0.3050</td><td>0.093</td><td>0.0284</td><td>0.0087</td><td>0.02626</td><td><strong>0.943</strong></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Plot of the pmf and partial autocorrelation function (for Pearson residuals) for the GP-INAR(1) model.</figcaption></figure>

NB, GP and GNB-INAR(1) models have very closed values of log-likelihood, contrarily to the B and P-INAR(1) models that obtained a poor fit. That could be explained by the high over dispersion observed in the dataset. The variance is clearly underestimated with these two models. On the contrary the variance is well estimated for the three other models, with ordered values 12.578 < 13.008 < 14.856 as expected, closed to the observed variance of 15.401. The autocorrelation of first order is well captured by the GP-INAR(1) model with $\hat{\rho} = 0.37$ compared with the observed value of 0.39. For the model adequacy, we can check the mean and variance of the standardized Pearson residuals. The NB-INAR(1) model obtained the mean closest to 0 and the GNB-INAR(1) model obtained the variance closest to 1. Finally, regarding the different criteria of model adequacy, the NB, GP and GNB-INAR(1) models cannot really be distinguished.

# 7 Conclusion

The unifying class of INAR(1) models based on quasi Pólya thinning operator has been introduced. It includes two well-known INAR(1) models, based on binomial and
quasi binomial thinning operators. It also includes two less usual INAR(1) models, based on hypergeometric and beta binomial thinning operators. Finally, it includes a new INAR(1) model based on quasi beta binomial thinning operator. This unified framework allowed us to derive the different properties of these models in one way. For instance, the class of additive modified power series distributions emerges as the natural choice for the error distribution. The asymptotic normality of the MLE for such a class of models has been shown under usual regularity conditions. Finally, the proposed class of INAR(1) models have been applied on simulated and real datasets.

The huge range of variability of this class allows taking into account different situations: under dispersion with B-INAR(1) model, equi dispersion with P-INAR(1) model and over dispersion with NB, GP and GNB-INAR(1) models. In real situations, the data are often over dispersed and the GP-INAR(1) model and GNB-INAR(1) model well perform, regarding the simulations results. In the example of real dataset, the NB, GP and GNB-INAR(1) models give very similar results regarding the different criteria. The GNB-INAR(1) model is more flexible since it includes the NB-INAR(1) model as a special case and has four parameters. But the additional parameter *d* seems not sufficient to clearly distinguish it from NB or GP-INAR(1) models. It may be due to the small order of the Markov chain. Therefore, it could be interesting to compare the performances of different integer-valued autoregressive models (INAR models), with a higher order. A natural perspective is thus to define INAR(*p*) models based on quasi Pólya thinning operator and finally INARMA(*p*, *q*) models. Until now, INAR(*p*) models with closure property (i.e., same distribution for margin and residual), has been defined only for the case of Poisson margin (Alzaid and Al-Osh, 1990; Weiß, 2008a) that do not take into account under and over dispersed data. The extension to CINAR(*p*) models based on quasi Pólya thinning operator will use the properties of such an operator (Lemmas 4 and 3). The extension of the INAR(*p*) model of Alzaid and Al-Osh (1990) to non-Poisson margins will need the multivariate version of the quasi Pólya thinning operator introduced by Peyhardi (2023). The stationarity of such extensions will use the same closure properties, under thinning operator and under convolution.

## Declarations

For the purpose of Open Access, a CC-BY 4.0 public copyright licence (https://creativecommons.org/licenses/by/4.0/) has been applied by the author to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

<figure><img src="image_5.png" /></figure>

## Conflict of interest

The author declare that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
<table><caption>**Table A2** Notations and pmf of additive modified power series distributions. The cell generalized binomial distribution GB<sup>[d]</sup>(θ, p) is greyed out since it is equivalent to the generalized negative binomial distribution GNB<sup>[d-1]</sup>(θ, p). After identifiability our generalized Poisson distribution GP<sup>[d]</sup>(αθ) becomes those introduced by Consul and Jain (1973) denoted by GP<sup>[d*]</sup>(λ) with λ = θα and d* = dα.</caption><thead><tr><th>Name</th><th>Notation</th><th>Parameters space</th><th>Supp</th><th>g(α)</th><th>h<sub>θ</sub>(α)</th><th>pmf</th></tr></thead><tbody><tr><td>additive modified power series<br>(c ∈ Z, d ∈ ℝ)</td><td>APS<sup>[c,d]</sup>(θ, α)</td><td>θ ∈ Θ<sub>c</sub>, α ∈ (0, R<sub>c</sub>)</td><td></td><td></td><td></td><td><sup>a<sub>θ</sub><sup>[c,d]</sup>(k)</sup>/<sub>h<sub>θ</sub>(α)</sub><sup>k</sup></td></tr><tr><td>1. binomial<br>(c = -1, d = 0)</td><td>B<sub>θ</sub>(p)</td><td>θ ∈ ℝ<sup>*</sup>, α ∈ (0, ∞)<br>p = α/(1+α) ∈ (0, 1)</td><td>{0, ..., θ}</td><td>α</td><td>(1 + α)<sup>θ</sup></td><td>(<sup>θ</sup>/<sub>n</sub>)(<sup>α</sup>/<sub>1+α</sub>)<sup>n</sup>(<sup>1</sup>/<sub>1+α</sub>)<sup>θ-n</sup></td></tr><tr><td>2. Poisson<br>(c = 0, d = 0)</td><td>P(αθ)</td><td>θ ∈ ℝ<sub>+</sub><br>α ∈ (0, ∞)</td><td>ℝ</td><td>α</td><td>e<sup>αθ</sup></td><td><sup>θ<sup>n</sup></sup>/<sub>n!</sub>α<sup>n</sup>e<sup>-αθ</sup></td></tr><tr><td>3. Negative binomial<br>(c = 1, d = 0)</td><td>NB(θ, p)</td><td>θ ∈ ℝ<sub>+</sub><sup>*</sup><br>p = α ∈ (0, 1)</td><td>ℝ</td><td>α</td><td>(1 - α)<sup>-θ</sup></td><td>(<sup>n+θ-1</sup>/<sub>n</sub>)α<sup>n</sup>(1 - α)<sup>θ</sup></td></tr><tr><td>generalized binomial<br>(c = -1, 1 &lt; d ≤ 1+α/α)</td><td>GB<sup>[d]</sup>(θ, p)</td><td>θ ∈ ℝ<sub>+</sub><sup>*</sup>, α ∈ (0, ∞)<br>p = α/(1+α) ∈ (0, 1)</td><td>ℝ</td><td>α(1 + α)<sup>-d</sup></td><td>(1 + α)<sup>θ</sup></td><td><sup>θ</sup>/<sub>θ+dn</sub>(<sup>θ+dn</sup>/<sub>n</sub>)(<sup>α</sup>/<sub>1+α</sub>)<sup>n</sup>(<sup>1</sup>/<sub>1+α</sub>)<sup>θ+dn-n</sup></td></tr><tr><td>4. generalized Poisson<br>(c = 0, 0 &lt; d ≤ 1)</td><td>GP<sup>[d]</sup>(αθ)</td><td>θ ∈ ℝ<sub>+</sub><sup>*</sup>, α ∈ (0, ∞)</td><td>ℝ</td><td>αe<sup>-dα</sup></td><td>e<sup>αθ</sup></td><td><sup>θ(θ+dn)<sup>n-1</sup></sup>/<sub>n!</sub>(<sup>αe<sup>-dα</sup></sup>/<sub>e<sup>αθ</sup></sub>)<sup>n</sup></td></tr><tr><td>5. generalized negative binomial<br>(c = 1, 0 &lt; d ≤ 1-α/α)</td><td>GNB<sup>[d]</sup>(θ, p)</td><td>θ ∈ ℝ<sub>+</sub><sup>*</sup>, p = α ∈ (0, 1)</td><td>ℝ</td><td>α(1 - α)<sup>d</sup></td><td>(1 - α)<sup>-θ</sup></td><td><sup>θ</sup>/<sub>θ+dn</sub>(<sup>θ+dn+n-1</sup>/<sub>n</sub>)α<sup>n</sup>(1 - α)<sup>θ+dn</sup></td></tr></tbody></table>

## Appendix B Standardized Pearson residuals

The standardized Pearson residuals are given by

$$e_t := \frac{E[X_t|x_{t-1}] - x_t}{\sqrt{V(X_t|x_{t-1})}}.$$

Let us describe the specific value of E[X_t|x_{t-1}] and V(X_t|x_{t-1}) in the five cases of INAR(1) models based on quasi Pólya thinning operator.

1. For the B-INAR(1) model with parameters ψ = (θ, γ, α) we have

$$E[X_t|x_{t-1}] = \frac{\theta}{\theta + \gamma}x_{t-1} + \gamma\frac{\alpha}{1+\alpha},$$

$$V(X_t|x_{t-1}) = \frac{\theta\gamma(\theta + \gamma - x_{t-1})x_{t-1}}{(\theta + \gamma)^2(\theta + \gamma - 1)} + \gamma\frac{\alpha}{(1 + \alpha)^2}.$$

2. For the P-INAR(1) model with parameters ψ = (ρ, α) we have

$$E[X_t|x_{t-1}] = \rho x_{t-1} + (1 - \rho)\lambda,$$
$$V(X_t|x_{t-1}) = \rho(1 - \rho)x_{t-1} + (1 - \rho)\lambda.$$

3. For the NB-INAR(1) model with parameters $\psi = (\theta, \gamma, \alpha)$ we have

$$ 
\begin{aligned}
E[X_t|x_{t-1}] &= \frac{\theta}{\theta + \gamma}x_{t-1} + \gamma\frac{\alpha}{1-\alpha}, \\
V(X_t|x_{t-1}) &= \frac{\theta\gamma(\theta + \gamma + x_{t-1})x_{t-1}}{(\theta + \gamma)^2(\theta + \gamma + 1)} + \gamma\frac{\alpha}{(1-\alpha)^2}.
\end{aligned}
$$

4. For the GP-INAR(1) model with parameters $\psi = (\rho, \lambda, d)$ we have

$$ 
\begin{aligned}
E[X_t|x_{t-1}] &= \rho x_{t-1} + (1 - \rho)\frac{\lambda}{1 - d}, \\
V(X_t|x_{t-1}) &= \rho(1 - \rho)\left[x_{t-1}^2 - \sum_{j=1}^{x_{t-1}-1} \frac{x_{t-1}^{(j+1;-1)} d'^{j-1}}{(1 + d'x_{t-1})^j}\right] + (1 - \rho)\frac{\lambda}{(1 - d)^3}.
\end{aligned}
$$

5. For the GNB-INAR(1) model with parameters $\psi = (\theta, \gamma, \alpha, d)$ we have

$$ 
\begin{aligned}
E[X_t|x_{t-1}] &= \frac{\theta}{\theta + \gamma}x_{t-1} + \gamma\frac{\alpha}{1 - \alpha(1 + d)}, \\
V(X_t|x_{t-1}) &= \frac{\theta\gamma}{\theta + \gamma}\left[\frac{x_{t-1}^2}{\theta + \gamma} - \sum_{j=1}^{x_{t-1}-1} \frac{x_{t-1}^{(j+1;-1)}(d+1)^{j-1}}{\{\theta + \gamma + (d+1)x_{t-1} - j + 2\}^{(j;1)}}\right] + \gamma\frac{\alpha(1 - \alpha)}{\{1 - \alpha(1 + d)\}^3}.
\end{aligned}
$$
Appendix C Summarized results on simulated datasets

<table><caption>**Table C3** Means (and standard deviations) of KL divergences for different INAR(1) models based on quasi Pólya thinning operator for time series of length *T* = 50 and autocorrelation *ρ* = 0.3. The rows indicate the marginal distribution of the simulated INAR(1) model whereas the columns indicate the marginal distribution of the estimated INAR(1) model. Means of KL divergences (and standard deviations) have been computed for the 5 different estimated models, on *n* = 1000 simulations for each of the 15 scenarios. The minimal value of each scenario is written in bold.</caption><thead><tr><th>Marginal</th><th>Mean</th><th>Variance</th><th>B</th><th>P</th><th>NB</th><th>GP</th><th>GNB</th></tr></thead><tbody><tr><td rowspan="3">B</td><td rowspan="3">μ = 20</td><td rowspan="3">σ² = μ/2</td><td><strong>4.17e-02</strong><br>(4.45e-02)</td><td>1.05e-01<br>(1.46e-02)</td><td>1.10e-01<br>(1.71e-02)</td><td>1.05e-01<br>(1.59e-02)</td><td>6.11e-02<br>(6.14e-02)</td></tr><tr><td><strong>3.58e-02</strong><br>(3.59e-02)</td><td>1.01e-01<br>(1.44e-02)</td><td>1.08e-01<br>(1.95e-02)</td><td>1.01e-01<br>(1.44e-02)</td><td>6.49e-02<br>(4.66e-02)</td></tr><tr><td><strong>3.47e-02</strong><br>(3.77e-02)</td><td>1.00e-01<br>(1.47e-02)</td><td>1.10e-01<br>(1.90e-02)</td><td>1.00e-01<br>(1.48e-02)</td><td>5.86e-02<br>(4.32e-02)</td></tr><tr><td rowspan="3">P</td><td rowspan="3">μ = 20</td><td rowspan="3">σ² = μ</td><td>2.49e-02<br>(3.14e-02)</td><td><strong>1.59e-02</strong><br>(2.19e-02)</td><td>1.85e-02<br>(2.34e-02)</td><td>1.85e-02<br>(2.25e-02)</td><td>2.64e-02<br>(3.36e-02)</td></tr><tr><td>2.73e-02<br>(3.62e-02)</td><td><strong>1.64e-02</strong><br>(2.46e-02)</td><td>1.98e-02<br>(2.62e-02)</td><td>1.93e-02<br>(2.49e-02)</td><td>2.62e-02<br>(3.30e-02)</td></tr><tr><td>2.75e-02<br>(3.36e-02)</td><td><strong>1.64e-02</strong><br>(2.23e-02)</td><td>1.89e-02<br>(2.34e-02)</td><td>1.89e-02<br>(2.28e-02)</td><td>2.62e-02<br>(3.48e-02)</td></tr><tr><td rowspan="3">NB</td><td rowspan="3">μ = 20</td><td rowspan="3">σ² = 2μ</td><td>1.99e-01<br>(7.52e-02)</td><td>1.86e-01<br>(6.53e-02)</td><td>4.08e-02<br>(4.32e-02)</td><td>3.82e-02<br>(3.85e-02)</td><td><strong>3.81e-02</strong><br>(3.84e-02)</td></tr><tr><td>2.08e-01<br>(7.05e-02)</td><td>1.91e-01<br>(6.49e-02)</td><td>4.35e-02<br>(3.98e-02)</td><td>4.42e-02<br>(3.94e-02)</td><td><strong>4.28e-02</strong><br>(3.78e-02)</td></tr><tr><td>2.13e-01<br>(6.66e-02)</td><td>1.92e-01<br>(6.11e-02)</td><td>4.35e-02<br>(3.84e-02)</td><td>4.48e-02<br>(3.79e-02)</td><td><strong>4.27e-02</strong><br>(3.60e-02)</td></tr><tr><td rowspan="3">GP</td><td rowspan="3">μ = 20</td><td rowspan="3">σ² = 4μ</td><td>6.91e-01<br>(1.19e-01)</td><td>6.65e-01<br>(1.09e-01)</td><td>4.27e-02<br>(4.95e-02)</td><td><strong>3.58e-02</strong><br>(4.09e-02)</td><td>3.83e-02<br>(4.16e-02)</td></tr><tr><td>7.85e-01<br>(1.23e-01)</td><td>7.52e-01<br>(1.08e-01)</td><td>3.61e-02<br>(3.79e-02)</td><td>3.33e-02<br>(3.53e-02)</td><td><strong>3.33e-02</strong><br>(3.56e-02)</td></tr><tr><td>8.33e-01<br>(1.18e-01)</td><td>9.15e-01<br>(4.00e+00)</td><td>3.96e-02<br>(4.61e-02)</td><td>1.63e-01<br>(4.02e+00)</td><td><strong>3.56e-02</strong><br>(3.92e-02)</td></tr><tr><td rowspan="3">GNB</td><td rowspan="3">μ = 20</td><td rowspan="3">σ² = 8μ</td><td>1.60e+00<br>(2.93e-01)</td><td>1.56e+00<br>(2.65e-01)</td><td>4.97e-02<br>(5.16e-02)</td><td><strong>3.61e-02</strong><br>(4.28e-02)</td><td>4.00e-02<br>(4.46e-02)</td></tr><tr><td>2.07e+00<br>(2.68e-01)</td><td>2.01e+00<br>(2.53e-01)</td><td>3.73e-02<br>(3.80e-02)</td><td><strong>3.16e-02</strong><br>(3.31e-02)</td><td>3.31e-02<br>(3.41e-02)</td></tr><tr><td>2.29e+00<br>(2.97e-01)</td><td>1.07e+01<br>(2.98e+01)</td><td>3.75e-02<br>(3.86e-02)</td><td>8.73e+00<br>(3.04e+01)</td><td><strong>3.53e-02</strong><br>(3.63e-02)</td></tr></tbody></table>
<table><caption>**Table C4** Means (and standard deviations) of KL divergences for different INAR(1) models based on quasi Pólya thinning operator for time series of length *T* = 50 and autocorrelation *ρ* = 0.9. The rows indicate the marginal distribution of the simulated INAR(1) model whereas the columns indicate the marginal distribution of the estimated INAR(1) model. Means of KL divergences (and standard deviations) have been computed for the 5 different estimated models, on *n* = 1000 simulations for each of the 15 scenarios. The minimal value of each scenario is written in bold.</caption><thead><tr><th>Marginal</th><th>Mean</th><th>Variance</th><th>B</th><th>P</th><th>NB</th><th>GP</th><th>GNB</th></tr></thead><tbody><tr><td></td><td>μ = 5</td><td></td><td>1.30e+03<br>(1.31e+04)</td><td>1.70e-01<br>(1.38e-01)</td><td>1.84e-01<br>(1.50e-01)</td><td><b>1.70e-01</b><br>(1.38e-01)</td><td>3.30e-01<br>(5.07e-01)</td></tr><tr><td rowspan="2">B</td><td>μ = 20</td><td>σ² = μ/2</td><td>2.73e-01<br>(3.58e-01)</td><td><b>1.75e-01</b><br>(1.66e-01)</td><td>1.87e-01<br>(1.74e-01)</td><td>1.77e-01<br>(1.64e-01)</td><td>2.10e-01<br>(2.68e-01)</td></tr><tr><td>μ = 50</td><td></td><td>3.37e-01<br>(5.53e-01)</td><td><b>1.83e-01</b><br>(1.70e-01)</td><td>2.05e-01<br>(2.04e-01)</td><td>1.95e-01<br>(1.87e-01)</td><td>2.36e-01<br>(2.80e-01)</td></tr><tr><td rowspan="2">P</td><td>μ = 20</td><td>σ² = μ</td><td>2.19e-01<br>(3.12e-01)</td><td>1.40e-01<br>(2.03e-01)</td><td>1.47e-01<br>(2.12e-01)</td><td><b>1.38e-01</b><br>(1.97e-01)</td><td>1.89e-01<br>(5.23e-01)</td></tr><tr><td>μ = 50</td><td></td><td>2.77e-01<br>(4.42e-01)</td><td><b>1.26e-01</b><br>(1.89e-01)</td><td>1.33e-01<br>(1.96e-01)</td><td>1.27e-01<br>(1.82e-01)</td><td>1.89e-01<br>(2.99e-01)</td></tr><tr><td rowspan="2">NB</td><td>μ = 20</td><td>σ² = 2μ</td><td>2.92e-01<br>(4.50e-01)</td><td><b>1.26e-01</b><br>(1.83e-01)</td><td>1.35e-01<br>(1.91e-01)</td><td>1.31e-01<br>(1.71e-01)</td><td>2.05e-01<br>(3.58e-01)</td></tr><tr><td>μ = 50</td><td></td><td>5.25e-01<br>(4.61e+00)</td><td>3.22e-01<br>(3.16e-01)</td><td>1.98e-01<br>(2.88e-01)</td><td><b>1.62e-01</b><br>(2.31e-01)</td><td>1.71e-01<br>(2.44e-01)</td></tr><tr><td rowspan="2">GP</td><td>μ = 20</td><td>σ² = 4μ</td><td>4.71e-01<br>(6.70e-01)</td><td>3.70e-01<br>(3.70e-01)</td><td>2.49e-01<br>(2.83e-01)</td><td><b>2.03e-01</b><br>(2.44e-01)</td><td>2.25e-01<br>(2.74e-01)</td></tr><tr><td>μ = 50</td><td></td><td>5.00e-01<br>(5.84e-01)</td><td>3.80e-01<br>(3.51e-01)</td><td>2.75e-01<br>(2.88e-01)</td><td><b>2.35e-01</b><br>(2.59e-01)</td><td>2.58e-01<br>(2.80e-01)</td></tr><tr><td rowspan="2">GNB</td><td>μ = 20</td><td>σ² = 4μ</td><td>5.74e+00<br>(1.44e+02)</td><td>8.87e-01<br>(6.93e-01)</td><td>2.50e-01<br>(3.49e-01)</td><td><b>1.95e-01</b><br>(2.77e-01)</td><td>2.00e-01<br>(2.88e-01)</td></tr><tr><td>μ = 50</td><td></td><td>1.02e+00<br>(7.72e-01)</td><td>9.34e-01<br>(6.28e-01)</td><td>2.50e-01<br>(3.08e-01)</td><td><b>1.86e-01</b><br>(2.39e-01)</td><td>1.98e-01<br>(3.30e-01)</td></tr><tr><td rowspan="2">GNB</td><td>μ = 20</td><td>σ² = 8μ</td><td>1.10e+00<br>(8.67e-01)</td><td>1.13e+00<br>(4.93e+00)</td><td>3.44e-01<br>(3.32e-01)</td><td>3.95e-01<br>(4.93e+00)</td><td><b>2.44e-01</b><br>(2.96e-01)</td></tr><tr><td>μ = 50</td><td></td><td>5.47e+01<br>(7.57e+02)</td><td>1.83e+00<br>(1.96e+00)</td><td>6.52e+00<br>(1.83e+02)</td><td><b>2.39e-01</b><br>(4.48e-01)</td><td>1.23e+00<br>(2.09e+01)</td></tr><tr><td rowspan="2">GNB</td><td>μ = 20</td><td>σ² = 8μ</td><td>2.26e+00<br>(1.70e+00)</td><td>2.13e+00<br>(1.15e+00)</td><td>2.26e-01<br>(3.14e-01)</td><td><b>1.54e-01</b><br>(2.13e-01)</td><td>1.65e-01<br>(2.56e-01)</td></tr><tr><td>μ = 50</td><td></td><td>2.40e+00<br>(1.11e+00)</td><td>7.10e+00<br>(2.96e+01)</td><td>2.72e-01<br>(3.16e-01)</td><td>5.18e+00<br>(2.99e+01)</td><td><b>2.25e-01</b><br>(1.19e+00)</td></tr></tbody></table>
References

Al-Osh M, Alzaid A (1987) First-order integer-valued autoregressive (inar (1)) process. Journal of Time Series Analysis 8(3):261–275

Al-Osh M, Alzaid A (1991) Binomial autoregressive moving average models. Stochastic Models 7(2):261–282

Al-Osh MA, Aly EEA (1992) First order autoregressive time series with negative binomial and geometric marginals. Communications in Statistics-Theory and Methods 21(9):2483–2492

Alzaid A, Al-Osh M (1990) An integer-valued pth-order autoregressive structure (inar (p)) process. Journal of Applied Probability 27(2):314–324

Alzaid A, Al-Osh M (1993) Some autoregressive moving average processes with generalized poisson marginal distributions. Annals of the Institute of Statistical Mathematics 45(2):223–232

Billingsley P (1961) Statistical methods in markov chains. The annals of Mathematical Statistics pp 12–40

Chu W (2010) Elementary proofs for convolution identities of Abel and Hagen–Rothe. The Electronic Journal of Combinatorics 24(1)

Consul P, Mittal S (1975) A new urn model with predetermined strategy. Biometrische Zeitschrift 17(2):67–75

Consul PC, Famoye F (2006) Lagrangian Probability Distributions. Springer

Consul PC, Jain GC (1973) A generalization of the poisson distribution. Technometrics 15(4):791–799

Gupta RC (1974) Modified power series distribution and some of its applications. Sankhyā: The Indian Journal of Statistics, Series B pp 288–298

Harvey AC, Fernandes C (1989) Time series models for count or qualitative observations. Journal of Business & Economic Statistics 7(4):407–417

Janardan K (1973) A new four urn model with predetermined strategy. In: Technical Report 37-I. Sangamon State University Springfield IL

Janardan K, Rao RB (1986) Identifiability of a generalized markov-polya damage model. Statistics: A Journal of Theoretical and Applied Statistics 17(2):303–310

Janardan KG, Schaeffer DJ (1977) A generalization of markov-polya distribution its extensions and applications. Biometrical Journal 19(2):87–106
Joe H (1996) Time series models with univariate margins in the convolution-closed infinitely divisible class. Journal of Applied Probability 33(3):664–677

Joe H (1997) Multivariate Models and Multivariate Dependence Concepts. CRC press

Johnson N, Kotz S, Balakrishnan N (1997) Discrete Multivariate Distributions. Wiley Series in Probability and Statistics, Wiley, URL https://books.google.fr/books?id=B4F9QgAACAAJ

Jorgensen B (1997) The Theory of Dispersion Models. CRC Press

Jourdan A, Kokonendji CC (2002) Surdispersion et modèle binomial négatif généralisé. Revue de Statistique Appliquée 50(3):73–86

Latour A (1998) Existence and stochastic structure of a non-negative integer-valued autoregressive process. Journal of Time Series Analysis 19(4):439–455

MacDonald IL, Zucchini W (2016) Hidden markov models for discrete-valued time series. Handbook of Discrete-Valued Time Series pp 267–286

McKenzie E (1985) Some simple models for discrete variate time series 1. JAWRA Journal of the American Water Resources Association 21(4):645–650

McKenzie E (1987) Innovation distributions for gamma and negative binomial autoregressions. Scandinavian Journal of Statistics pp 79–85

Patil G, Ratnaparkhi M (1975) Problems of damaged random variables and related characterizations. In: A Modern Course on Statistical Distributions in Scientific Work. Springer, p 255–270

Peyhardi J (2023) On quasi pólya thinning operator. Brazilian Journal of Probability and Statistics

Rao CR (1965) On discrete distributions arising out of methods of ascertainment. Sankhyā: The Indian Journal of Statistics, Series A pp 311–324

Ristić MM, Bakouch HS, Nastić AS (2009) A new geometric first-order integer-valued autoregressive (nginar (1)) process. Journal of Statistical Planning and Inference 139(7):2218–2226

Sellers KF, Peng SJ, Arab A (2020) A flexible univariate autoregressive time-series model for dispersed count data. Journal of Time Series Analysis 41(3):436–453

Weiß CH (2008a) The combined inar (p) models for time series of counts. Statistics & Probability Letters 78(13):1817–1822

Weiß CH (2008b) Thinning operations for modeling time series of counts—a survey. AStA Advances in Statistical Analysis 92(3):319–341
Zheng H, Basawa IV, Datta S (2007) First-order random coefficient integer-valued autoregressive processes. Journal of Statistical Planning and Inference 137(1):212–229