Non-linear INAR(1) processes under an alternative geometric thinning operator

We propose a novel class of first-order integer-valued AutoRegressive (INAR(1)) models based on a new operator, the so-called geometric thinning operator, which induces a certain non-linearity to the models. We show that this non-linearity can produce better results in terms of prediction when compared to the linear case commonly considered in the literature. The new models are named non-linear INAR(1) (in short NonLINAR(1)) processes. We explore both stationary and non-stationary versions of the NonLINAR processes. Inference on the model parameters is addressed and the finite-sample behavior of the estimators investigated through Monte Carlo simulations. Two real data sets are analyzed to illustrate the stationary and non-stationary cases and the gain of the non-linearity induced for our method over the existing linear methods. A generalization of the geometric thinning operator and an associated NonLINAR process are also proposed and motivated for dealing with zero-inflated or zero-deflated count time series data.


Introduction
A first-order integer-valued autoregressive (INAR(1)) process is defined by a sequence of non-negative integer-valued random variables {X t } t∈N satisfying where "•" is the binomial thinning operator by Steutel and van Harn (1979), which is defined by α • X t−1 ≡ X t−1 k=1 ζ t,k , {ζ t,k } t,k∈N is a doubly infinite array of independent and identically distributed (iid) Bernoulli random variables with P (ζ t,k = 1) = 1 − P (ζ t,k = 0) = α ∈ (0, 1), and { t } t∈N is assumed to be a sequence of iid non-negative integer-valued random variables, with t independent of X s−1 and ζ s,k , for all k ≥ 1 and for all s ≤ t. INAR processes have been commonly used to fit count time series data. In the context of branching process, the random variable X t can be seen as the total population at time t, α • X t−1 is the number of survivors at time t − 1, while t stands for the immigration at time t. The random variable ζ t,k denotes the survival or not of the k-th individual of the population at time t − 1. Observe that, under the assumption that μ ≡ E( t ) < ∞, the conditional expectation of X t given X t−1 is linear on X t−1 , that is, E(X t |X t−1 ) = α X t−1 + μ . (2) The INAR process given in (1) has been studied by Alzaid and Al-Osh (1987), McKenzie (1988), and Dion et al. (1995). Conditional least squares estimation for this model was explored, for instance, by Wei and Winnicki (1990), Ispány et al. (2003), Freeland and McCabe (2005), and Rahimov (2008). For a comprehensive review on thinning-based INAR processes and some of their generalizations, see Scotto et al. (2015).
Alternative integer-valued processes based on non-additive innovation through maximum and minimum operations have been proposed in the literature. Littlejohn (1992) and Littlejohn (1996) considered a discrete minification processes based on a thickening operator; see also Kalamkar (1995) for an alternative class of minification processes. Scotto et al. (2016) defined a class of max-INAR processes by X t = max{α • X t−1 , t }, for t ≥ 1, with "•" being the binomial thinning operator, while Aleksić and Ristić (2021) used the minimum function and a modified negative binomial INAR-operator to define their processes as X t = min{α X t−1 , t }, where α X = 1+X i=1 G i , for a non-negative random variable X , with {G i } i∈N being a sequence of iid geometric distributed random variables with mean α > 0, where { t } t∈N is defined as before.
For the count processes {X t } t∈N considered in these works, a certain non-linearity is induced in the sense that the conditional expectation E(X t |X t−1 ) is non-linear on X t−1 (and also the conditional variance) in contrast with (2). We refer to these models as "non-linear" throughout this paper. On the other hand, the immigration interpretation in a populational context is lost due to the non-additive innovation assumption.
Our chief goal in this paper is to propose a novel class of non-linear INAR(1) processes that preserve the additive innovation, while still having a practical interpretation, in contrast with the existing non-linear INAR models, where this interpretation is lost. With that purpose we introduce a new INAR-operator, the so-called geometric thinning operator, which is of own interest. The new models are named non-linear INAR(1) (in short NonLINAR(1)) processes. We show that the proposed NonLINAR(1) processes can produce better results in terms of prediction when compared to the linear case, which is commonly considered in the literature. We now highlight other contributions of the present paper: (i) development of inferential procedures and numerical experiments, which are not well-explored for the existing non-linear models aforementioned; (ii) properties of the novel geometric thinning operator are established; (iii) a particular NonLINAR(1) process with geometric marginals is investigated in detail, including an explicit expression for the autocorrelation function; (iv) both stationary and non-stationary cases are explored, being the latter important for allowing the inclusion of covariates, a feature not considered by the aforementioned papers on non-linear INAR models; (v) a generalization of the geometric thinning operator and an associated NonLINAR process are also proposed and motivated for dealing with zero-inflated or zerodeflated count time series data.
The paper is organized as follows. In Sect. 2, we introduce the new geometric thinning operator and explore its properties. Section 3 is devoted to the development of the NonLINAR processes based on the new operator, with a focus on the case where the marginals are geometrically distributed. Two methods for estimating the model parameters are discussed in Sect. 4, including Monte Carlo simulations to evaluate the proposed estimators. In Sect. 5, we introduce a non-stationary NonLINAR process allowing for the inclusion of covariates and provide some Monte Carlo studies. Section 6 is devoted to two real data applications. Finally, in Sect. 7, we develop a generalization of the geometric thinning operator and an associated NonLINAR model.

Geometric thinning operator: definition and properties
In this section, we introduce a new thinning operator and derive its main properties. We begin by introducing some notation. For two random variables X and Y , we write min{X , Y } := X ∧ Y to denote the minimum between X and Y . The probability generating function (pgf) of a non-negative integer-valued random variable Y is denoted by for all values of s for which the right-hand side converges absolutely. The n-th derivative of Y (x) with respect to x and evaluated at x = x 0 is denoted by (n) Y (x 0 ). Let Z be a geometric random variable with parameter α > 0 and probability function assuming the form In this case, the pgf of X is and the parameter α has the interpretation α = E(Z ) > 0. The shorthand notation Z ∼ Geo(α) will be used throughout the text. We are ready to introduce the new operator and explore some of its properties.
Definition 1 (Geometric thinning operator) Let X be a non-negative integer-valued random variable, independent of Z ∼ Geo(α), with α > 0. The geometric thinning operator is defined by Remark 1 The operator defined in (4) satisfies α X ≤ X , like the classic binomial thinning operator •. Therefore, is indeed a thinning operator.
In what follows, we present some properties of the proposed geometric thinning operator. We start by obtaining its probability generating function.
Proposition 1 Let X be a non-negative integer-valued random variable with pgf X . Then, the pgf of α X is given by Proof By the independence assumption between X and Z , it holds that Hence, The second term on the last equality can be expressed as The result follows by rearranging the terms.
The next result gives us the moments of α X , which will be important to discuss prediction and forecasting in what follows.
Proposition 2 Let be the geometric thinning operator in (4). It holds that the n-th factorial moment of α X is given by Proof The result follows by using the pgf given in Proposition 1 and the generalized Leibniz rule for derivatives, namely In what follows, the notation X ⇒ Y means X weakly converges to Y .
Proposition 3 Let be the geometric thinning operator in (4). Then, Proof The proof follows immediately from Proposition 1 and the Continuity Theorem for pgf's.
We now show a property of the operator of own interest.
Proposition 4 Let Z 1 , . . . , Z n be independent geometric random variables with parameters α 1 , . . . , α n , respectively. Assume that X 1 , . . . , X n are non-negative integer-valued random variables independent of the Z 's, and let α i X i = min (X i , Z i ). Then, Proof We prove (5) by induction on n. For n = 2, it holds that the proof is complete.
We finish this section by discussing the zero-modified geometric (ZMG) distribution and some of its properties. Such a distribution will play an important role in the construction of our model in Sect. 3. We say that a random variable Y follows a ZMG distribution with parameters μ > 0 and π ∈ (−1/μ, 1) if its probability function is given by We denote Y ∼ ZMG(π, μ). The geometric distribution with mean μ is obtained as a particular case when π = 0. For π < 0 and π > 0, the ZMG distribution is zero-deflated or zero-inflated with relation to the geometric distribution, respectively.
For π ∈ (0, 1), Y satisfies the following equality in distribution: Y d = B Z, where B and Z are independent random variables with Bernoulli (with success parameter π ) and geometric (with mean μ) distributions, respectively. The associated pgf can be computed as Now, assume that X ∼ Geo(μ), with μ > 0. We have that In the next section, we introduce our class of non-linear INAR processes and provide some of their properties.

Non-linear INAR(1) processes
In this section, we introduce a novel class of non-linear INAR(1) processes based on the new geometric thinning operator defined in Sect. 2 and explore a special case when the marginals are geometrically distributed.
Definition 2 A sequence of random variables {X t } t∈N is said to be a non-linear INAR(1) process (in short NonLINAR(1)) if it satisfies the stochastic equation with α X t−1 = min (X t−1 , Z t ), {Z t } t∈N being a sequence of iid random variables with Z 1 ∼Geo(α), { t } t∈N being an iid non-negative integer-valued random variables called innovations, where t is independent of X t−l and Z t−l+1 , for all l ≥ 1, and X 0 is some starting non-negative value/random variable.

Remark 2
The random variable Z t in Definition 2 determines the number of survivals at time t. Either the previous population is reduced to Z t individuals if Z t < X t−1 or everybody survives if Z t ≥ X t−1 . As argued in Remark 1, is a thinning operator. We remark that the minimum operation to construct count processes is also known in the literature as minimization; for instance, see Littlejohn (1992). Non-linearity for count time series models is achieved by using minimum or maximum operations, but differently from the existing literature, our proposed methodology induces nonlinearity to the model and keeps the additive innovation assumption at the same time. Therefore, the populational interpretation (with survivals and immigration processes) is still valid under our framework.
We now obtain some properties of the NonLINAR(1) processes.
Theorem 6 Let {X t } t∈N be a NonLINAR(1) process. Then, {X t } t∈N is stationary and ergodic.
Proof Note that {X t } t∈N is a homogeneous discrete time Markov chain (see Proposition (5)). Therefore, using the Markov property, it is not hard to see that (X 1 , . . . , X m ) and (X k , . . . , X k+m ) have the same distribution, for all m ∈ N and for all k ∈ N. This gives the stationarity of the process. To prove ergodicity, let σ (X t , X t+1 , . . . ) be the sigma-field generated by the random variables X t , X t+1 . . . . Definition (7) leads to Since T is a tail sigma-field, it follows by Kolmogorov's 0-1 Law that every event A ∈ T has probability 0 or 1. It follows that the process is ergodic (see Shiryaev (2019), Definition 2, pg. 43).

Proposition 7
The joint pgf of the discrete random vector (X t , X t−1 ) is given by with X (·) being the pgf of X and (·) as in (11), where s 1 and s 2 belong to some intervals containing the value 1.
Proof We have that Therefore, Proposition 8 The 1-step ahead conditional mean and conditional variance are given by

respectively.
Proof From the definition of the NonLINAR(1) processes, we obtain that for all x = 0, 1, . . . . The conditional expectation above can be obtained from Proposition 2 with X being a degenerate random variable at x (i.e. P(X = x) = 1). Then, it follows that The conditional variance can be derived analogously, so details are omitted.

Remark 3
Note that the conditional expectation and variance given in Proposition 8 are non-linear on X t−1 in contrast with the classic INAR processes where they are linear.
From now on, we focus our attention on a special case from our class of non-linear INAR processes when the marginals are geometrically distributed. From (7), it follows that a NonLINAR(1) process with geometric marginals is well-defined if the function is a proper pgf, with s belonging to some interval containing the value 1, where X (s) and α X (s) are the pgf's of geometric distributions with means μ and αμ 1 + α + μ , respectively. More specifically, we have which corresponds to the pgf of a zero-modified geometric distribution with parameters μ and α/(1 + μ + α); see (6). This shows that a NonLINAR(1) process with geometric marginals is well-defined.

Definition 3 The stationary geometric non-linear INAR (Geo-NonLINAR) process
{X t } t∈N is defined by assuming that (7) Remark 4 Note that imposing a geometric distribution for the marginals of the NonLI-NAR process implies that the innovations are ZMG distributed. Conversely, assuming a ZMG distribution as above for the innovations implies that the marginals are geometrically distributed. Therefore, Definition 3 ensures that the process has geometric marginals.
From (11), we have that the mean and variance of the innovations { t } t≥1 are given by respectively. Additionally, the third and forth moments of the innovations are

Proposition 9
The autocovariance and autocorrelation functions at lag 1 of the Geo-NonLINAR process are respectively given by After some algebra, the result follows.
In the following proposition, we obtain an expression for the conditional expectation E(X t |X t−k = ). This function will be important to find the autocovariance function at lag k ∈ N and to perform prediction and/or forecasting.
for all integer k ≥ 2.
Proof Let F t = σ (X 1 , . . . , X t ) denote the sigma-field generated by the random variables X 1 , . . . , X t . By the Markov property it is clear that for all k ≥ 1. The proof proceeds by induction on k. Equation (14) and Proposition 8 give us that with α * = α 1+α . Using (10), we obtain that Assume that (13) is true for k = n − 1. Using (14), we have From the definition of H n , we obtain Note that Proposition 11 Let h j , g j be as in Proposition 10 and writeh j = μh j , for j ∈ N. It holds that Proof Note that where the third equality follows by (13). A thorough inspection of the definition of H k gives Note that the argument of the function in (16) is just a constant times the derivative of X 1 (s) with respect to s and evaluated at α * . More specifically, . (17) The second equality follows from (3). Plugging (17) in (16), we obtain The result follows by plugging (18) in (15).

Parameter estimation
In this section, we discuss estimation procedures for the geometric NonLINAR process through conditional least squares (CLS) and maximum likelihood methods. We assume that X 1 , . . . , X n is a trajectory from the Geo-NonLINAR model with observed values x 1 , . . . , x n , where n stands for the sample size. We denote the parameter vector by θ ≡ (μ, α) . For the CLS method, we define the function Q n (θ) as The CLS estimators are obtained as the argument that minimizes Q n (θ), i.e.
Since we do not have an explicit expression forθ cls , numerical optimization methods are required to solve (20). This can be done through optimizer packages implemented in softwares such as R (R Core Team 2021) and MATLAB. The gradient function associated with Q n (·) can be provided for these numerical optimizers and is given by A strategy to get the standard errors of the CLS estimates based on bootstrap is proposed and illustrated in our empirical illustrations; please see Sect. 6.
We now discuss the maximum likelihood estimation (MLE) method. Note that our proposed Geo-NonLINAR process is a Markov chain (by definition) and therefore the likelihood function can be expressed in terms of the 1-step transition probabilities derived in Proposition 5. The MLE estimators are obtained as the argument that maximizes the log-likelihood function, that is,θ mle = arg max θ n (θ ), with where the conditional probabilities in (21) are given by (8) and P (X 1 = x 1 ) is the probability function of a geometric distribution with mean μ. There is no closed-form expression available forθ mle . The maximization of (21) can be accomplished through numerical methods such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm implemented in the R package optim. The standard errors of the maximum likelihood estimates can be obtained by using the Hessian matrix associated with (21), which can be evaluated numerically.
In the remaining of this section, we examine and compare the finite-sample behavior of the CLS and MLE methods via Monte Carlo (MC) simulation with 1000 replications per set of parameter configurations, with the parameter estimates computed under both approaches. All the numerical experiments presented in this paper were carried out using the R programming language.
We consider four simulation scenarios with different values for θ = (μ, α) , namely: (I) θ = (2.0, 1.0) , (II) θ = (1.2, 0.5) , (III) θ = (0.5, 1.5) , and (IV) θ = (0.3, 0.5) . To illustrate these configurations, we display in Fig. 1 simulated trajectories from the Geo-NonLINAR process and their associated autocorrelation function (ACF) and partial autocorrelation function (PACF) under Scenarios I and IV. In Table 1, we report the empirical mean and root mean squared error (RMSE) of the parameter estimates obtained from the MC simulation based on the MLE and CLS methods. We can observe that both approaches produce satisfactory results and also a slight advantage of the MLE estimators over the CLS for estimating α, mainly in terms of RMSE, which is already expected. This advantage can also be seen from Fig. 2, which presents boxplots of the parameter estimates for μ and α under the Scenarios I and IV with sample sizes n = 100 and n = 500. In general, the estimation procedures considered here produce estimates with bias and RMSE decreasing towards zero as the sample size increases, therefore giving evidence of consistency. We also present in Fig. 3 the histograms of the standardized Monte Carlo estimates of μ and α under MLE and CLS approaches along with standard normal density curves. From this figure, we can observe a good normal approximation of the proposed estimators.

Dealing with non-stationarity
In many practical situations, stationarity can be a non-realistic assumption; for instance, see Brännäs (1995), Enciso-Mora et al. (2009), andWang (2020) for works that investigate non-stationary Poisson INAR process. Motivated by that, in this section, we propose a non-stationary version of the Geo-NonLINAR process allowing for time-varying parameters. Consider where w t and v t are p × 1 and q × 1 covariate vectors for t ≥ 1, and β and γ are p × 1 and q × 1 vectors of associated regression coefficients. We define a time-varying or non-stationary Geo-NonLINAR process by and It is also assumed that t is independent of X t−l and Z t−l+1 , for all l ≥ 1. Under these assumptions, the marginals of the process (22) are Geo(μ t ) distributed, for t ∈ N. This claim can be proved by following the same steps as for the stationary case. Remark 5 The transition probabilities and conditional mean and variance for the nonstationary Geo-NonLINAR(1) process are given by expressions in Propositions 5 and 8, respectively, just by replacing μ and α by μ t and α t .
We consider two estimation methods for the parameter vector θ = (β, γ ) . The first one is based on the conditional least squares. The CLS estimator of θ is obtained by minimizing (19) with μ t and α t instead of μ and α, respectively. According to Wang (2020), this procedure might not be accurate in the sense that non-significant covariates can be included in the model. In that paper, a penalized CLS (PCLS) method is considered. Hence, a more accurate estimator is obtained by minimizing Q n (θ) = Q n (θ )+n p+q j=1 P δ (|θ i |), where P δ (·) is a penalty function and δ is a tuning parameter. See Wang (2020) for possible choices of penalty function. This can be used as a selection criterion and we hope to explore it in a future paper. A second method for estimating the parameters is the maximum likelihood method. The log-likelihood function assumes the form (21) with μ and α replaced by μ t and α t , respectively.
For the non-stationary case, we carry out a second set of Monte Carlo simulations by considering trend and seasonal covariates in the model as follows: μ t = exp(β 0 + β 1 t/n + β 2 cos(2π t/12)) and α t = exp(γ 0 + γ 1 t/n), for t = 1, . . . , n. The above structure aims to mimic realistic situations when dealing with epidemic diseases. We here set the following scenarios: (V) (β 0 , β 1 , β 2 , γ 0 , γ 1 ) = (2.0, 1.0, 0.7, 2.0, 1.0) and (VI) (β 0 , β 1 , β 2 , γ 0 , γ 1 ) = (3.0, 1.0, 0.5, 3.0, 2.0). We  Table 2 reports the empirical mean and the RMSE (within parentheses) of the parameter estimates based on the MLE and CLS methods. We can observed that the MLE method outperforms the CLS method for all configurations considered, as expected since we are generating time series data from the "true" model. This can be also seen from Fig. 4, which presents the boxplots of MLE and CLS estimates under the Scenarios V with sample sizes n = 200, 500. Regardless, note that the bias and RMSE of the CLS estimates decrease as the sample size increases. Figure 5 displays the standardized Monte Carlo estimates under the MLE and CLS methods along with the standard normal density curve. As in the stationary case, we can observe a good normal approximation, which is more satisfactory under the MLE method since it uses the full distributional assumption and we are generating data from the correct model.

Real data applications
In this section, we discuss the usefulness of our methodology under stationary and non-stationary conditions. In the first empirical example, we consider the monthly number of polio cases reported to the U.S. Centers for Disease Control and Prevention from January 1970 to December 1983, with 168 observations. The data were obtained through the gamlss package in R. Polio (or poliomyelitis) is a disease caused by poliovirus. Symptoms associated with polio can vary from mild flu-like symptoms   disease that is caused by M. leprae. It mainly affects the skin, the peripheral nerves mucosa of the upper respiratory tract, and the eyes. According to the World Health Organization, about 208,000 people worldwide are infected with Hansen's disease. The data are displayed in Table 3.

Polio data analysis
We begin the analysis of the polio data by providing plots of the observed time series and the corresponding sample ACF and PACF plots in Fig. 6. These plots give us evidence that the count time series is stationary. Table 4 provides a summary of the polio data with descriptive statistics, including mean, median, variance, skewness, and kurtosis. From the results in Table 4, we can observe that counts vary between 0 and 14, with the sample mean and variance equal to 1.333 and 3.505, respectively, which suggests overdispersion of the data.
For comparison purposes, we consider the classic first-order INAR process with E(X t |X t−1 ) = κ X t−1 + μ(1 − κ), where μ = E(X t ) and κ = corr(X t , X t−1 ) ∈ (0, 1). This linear conditional expectation on X t−1 holds for the classic stationary INAR processes such as the binomial thinning-based ones, in particular, the Poisson  INAR(1) model by Alzaid and Al-Osh (1987). The aim is to evaluate the effect of the nonlinearity of our proposed models on the prediction in comparison to the classic INAR(1) processes. We consider the CLS estimation procedure, where just the conditional expectation is considered. This allows for a more flexible approach since no further assumptions are required. To obtain the standard errors of the CLS estimates, we consider a parametric bootstrap where some model satisfying the specific form for the conditional expectation holds. In this first application, for our NonLINAR process, we consider the geometric model derived in Sect. 3. For the classic INAR, the Poisson model by Alzaid and Al-Osh (1987) is considered in the bootstrap approach. This strategy to get standard errors has been considered, for example, by Maia et al. (2021) for a class of semiparametric time series models driven by a latent factor. In order to compare the predictive performance of the competing models, we compute the sum of squared prediction errors (SSPE) defined by SSPE = n t=2 (x t −μ t ) 2 , whereμ t = E(X t |X t−1 ) is the predicted mean at time t (see Proposition 8), for t = 2, . . . , n. Table 5 summarizes the fitted models by providing CLS estimates and their respective standard errors, and the SSPE values. The SSPE results in Table 5 show the superior performance of the NonLINAR process over the classic INAR process in terms of prediction. This can also be observed from Fig. 7, where the NonLINAR process shows a better agreement between the observed and predicted values.  1970 1972 1974 1976 1978 1980 1982 1984 1970 1972 1974 1976 1978 1980 1982 1984  To evaluate the adequacy of our proposed NonLINAR process, we consider the Pearson residuals defined by R t ≡ (X t −μ t )/σ t , withσ t = Var(X t |X t−1 ), for t = 2, . . . , n, where we assume that the conditional variance takes the form given in Proposition 8. Figure 8 presents the Pearson residuals against the time, its ACF, and the qq-plot against the normal quantiles. These plots show that the data correlation was well-captured. On the other hand, the qq-plot suggests that the Pearson residuals are not normally distributed. Actually, this discrepancy is not unusual especially when dealing with low counts; for instance, see Zhu (2011) and Silva and Barreto-Souza (2019). As an alternative way to check the adequacy, we use the normal pseudoresiduals introduced by Dunn and Smyth (1996), which is defined by R * t = −1 (U t ), where (·) is the standard normal distribution function and U t is uniformly distributed on the interval (Fθ (x t − 1), Fθ (x t )), where Fθ (·) is the fitted predictive cumulative distribution function of the NonLINAR process. Figure 9 shows the pseudo residuals against the time, its ACF, and qq-plot. We can observe that the pseudo-residuals are not correlated and are approximately normally distributed. Therefore, we conclude that the NonLINAR process provides an adequate fit to the polio count time series data.
We now analyze the predictive performance of the proposed model by conducting an out-of-sample forecasting exercise through a rolling estimation window  approach. More specifically, we split the data X 1 , . . . , X n into the first n 0 observations X 1 , . . . , X n 0 and the remaining time series X n 0 +1 , . . . , X n , where n 0 < n. Hence, we estimate the model parameters using the trajectory X 1 , . . . , X n 0 and forecast X n 0 +1 by using the conditional expectation given in Proposition 8. Thereafter, we update the training dataset including the observation X n 0 +1 and reestimate the model parameters using X 1 , . . . , X n 0 , X n 0 +1 . Based on this fitted model, we forecast X n 0 +2 using the conditional 1-step ahead expectation as before. This procedure is repeated until we reach the last observation.
For the polio data, we consider n 0 = 84, which corresponds to December 1976. Figure 10 displays the polio time series and the 1-step ahead predicted values. This figure reveals a satisfactory out-of-sample forecasting performance of the NonLINAR(1) process since there is a good agreement between the observed and predicted values.

Hansen's disease data analysis
We now analyze Hansen's disease data. A descriptive data analysis is provided in Table 6. Figure 11 presents the Hansen's count data and its corresponding sample ACF and PACF plots. This figure provides evidence that the count time series is non-stationary. In particular, we can observe a negative trend. This motivates us to use non-stationarity approaches to handle this data. We consider our non-stationary NonLINAR process with conditional mean where the following regression structure is assumed: with the term t/252 being a linear trend. For comparison purposes, we also consider the Poisson INAR(1) process allowing for covariates (Brännäs 1995 where μ t = exp (β 0 + β 1 t/252) and κ t = exp (ξ 0 + ξ 1 t/252) 1 + exp (ξ 0 + ξ 1 t/252) , for t = 1, . . . , 252.  We consider the CLS estimation procedure for both approaches considered here. Table 7 gives us the parameter estimates under the NonLINAR and PINAR(1) processes, standard errors obtained via bootstrap, and the SSPE values (we use Eq. (23) and E(X t |X t−1 ) = κ t X t−1 + μ t (1 − κ t ) to obtain the predicted values according the non-stationary Geo-NonLINAR and Poisson INAR models, respectively). To get the standard errors for the parameter estimates, we proceed similarly as done in the first application with a slight difference. Since here the counts are high, the geometric assumption cannot be valid. Therefore, we consider a non-stationary NonLINAR process with innovations following a Poisson distribution with mean μ t (1 + μ t ) 1 + μ t + α t in our bootstrap scheme. This ensures that the conditional mean is the same as in (23). From Table 7, we have that the trend is significant (using, for example, a significance level at 5%) to explain the marginal mean μ t , but not for the parameter α t , under the NonLINAR model. Furthermore, we note that the sign of the estimate of β is negative, which is in agreement with the observed negative trend. We highlight that the parameter μ t also appears in the autocorrelation structure under our approach, therefore the trend is also significant to explain the autocorrelation of the NonLINAR process. By looking at the results from the PINAR fitting, we see that the trend is significant to explain α t (parameter related to the autocorrelation) but not the marginal mean μ t . Once again, we have that the model producing the smallest SSPE is the NonLINAR process. So, our proposed methodology is performing better than the classic PINAR  Plots of Hansen's disease data (solid line) and fitted conditional means (dots) based on the nonstationary NonLINAR process (to the left) and PINAR process (to the right) model in terms of prediction. The predictive values according to both models along with the observed counts are exhibited in Fig. 12. We now check if the non-stationary NonLINAR process fits well the data. Figure 13 provides the Pearson residuals against time, its ACF plot, and the qq-plot of the residuals. By looking at this figure, we have evidence of the adequacy of the NonLINAR process to fit Hansen's disease data.
We conclude this data analysis exploring the predictive power of the NonLINAR(1) process by performing the out-of-sample forecasting exercise through a rolling estimation window as described at the end of Sect. 6.1. Here, we consider n 0 = 168 (December 2014). Figure 14 provides the plot of Hansen's data and the 1-step ahead predictions. Once again, we notice a good agreement between the observed time series and the predictions.

Generalization
In this section, we provide an extension of the geometric thinning operator and propose a non-linear INAR process based on such generalization. As we will see, alternative distributions rather than geometric for the operation in (4) can provide flexible approaches for dealing with different features on count time series. We also discuss how to handle zero-inflation or zero-deflation with respect to the geometric model.
Remark 6 Note that the ZMG operator given in (24) has the geometric thinning operator as a special case when η = 1 since Z (1,α) ∼ Geo(α). Further, we stress that the parameterization of the ZMG distribution in terms of 1 − η instead of η will be convenient in what follows. Also, we will omit the dependence of Z on (η, α) to simplify the notation.
Based on the ZMG operator, we can define a non-linear INAR process {X t } t∈N (similarly as done in Sect. 3) by where (η, α) X t−1 = min (X t−1 , Z t ), with {Z t } t∈N iid ∼ ZMG(1 − η, α), { t } t≥1 is a sequence of iid non-negative integer-valued random variables, called innovations, with t independent of X t−l and Z t−l+1 , for all l ≥ 1, with X 0 being some starting value/random variable. This is basically the same idea as before; we are just replacing the geometric assumption by the zero-modified geometric law in the thinning operation.
We now show that it is possible to construct a stationary Markov chain satisfying (25) and having marginals ZMG-distributed; this could be seen as an alternative model to the zero-modified geometric INAR(1) process proposed by Barreto-Souza (2015). Furthermore, we argue that such construction is not possible under the geometric thinning operator defined in Sect. 2 (see Remark 7 below), which motivates the ZMG thinning introduced here.

Remark 7
Note that we are excluding the case η = 1 (which corresponds to the geometric thinning operator) since the required inequality 1−π 1−πη 1 + 1+μ α < 1 does not hold in this case (1 + 1+μ α > 1). This shows that a NonLINAR process with ZMG marginals cannot be constructed based on the geometric thinning operator defined previously and therefore motivates the ZMG operator. We would like to highlight the importance of the ZMG thinning operator (which is an extension of the geometric thinning operator) since it permits us to construct a Non-LINAR(1) process with ZMG marginals. As a consequence, this model can handle inflation or deflation of zeros that cannot be accounted by the geometric Non-LINAR(1) model introduced in Sect. 3.