Aldo M Garay
email: agaray@de.ufpe.br

Francyelle L Medina

Isaac Jales
email: isaacjales@uern.br

Patrice Bertail
email: patrice.bertail@gmail.com

First-order integer valued AR processes with zero-inated innovations

Keywords: Count time series, ZINAR(1) process, Bootstrap, EM algorithm

To deal with time series process with excess of zeros, we extend the INAR(1) process by considering that the innovations follow dierent zeroinated models, called the ZI-INAR(1) model. We present some of its theoretical properties, develop an ecient EM algorithm for parameter estimation and propose several bootstrap techniques to construct condence intervals for the parameters. Finally, we present the relevance and applicability the of proposed ZI-INAR(1) model through simulation studies and an application to a real dataset.

To analyze discrete time series with excess of zeros, in this manuscript we extend the ZINAR(1) model, considering that the innovations follow a class of zero-inated (ZI) models, called the ZI-INAR(1) model. Some properties, such as mean, variance and joint distribution are developed. We propose an EM algorithm to estimate the unknown parameters by maximizing the conditional likelihood function. For this purpose, we start by introducing some denitions, notations and properties about the rst-order integer valued autoregressive processes INAR [START_REF] Abramowitz | Handbook of mathematical functions: With formulas, graphs, and mathematical tables applied mathematics series[END_REF] and the zero-inated (ZI) models.

Denition 1 Let X be a non-negative integer-valued random variable (r.v.) and α ∈ [ 0, 1]. Then, the thinning operator `•' [START_REF] Steutel | Discrete analogues of self-decomposability and stability[END_REF] is dened as:

α • X = X i=1 Z i , (1)
where {Z i } i≥1 is a sequence of independent identically distributed (iid) Bernoulli random variables independent of X, with P (Z i = 1) = α.

Thus, considering Denition 1, the INAR(1) processes have the following stochastic structure:

Y t = α • Y t-1 + V t , t ∈ Z, (2)
where {V t } t∈Z is a sequence of non-negative integer-valued random iid variables, called innovations, with E(V t ) = µ < ∞, V ar(V t ) = σ 2 < ∞ and independent of Y t-1 , for all t.

It is important to note that given Y t-1 = y t-1 > 0, the r.v. α • Y t-1 follows a binomial distribution, with parameters y t-1 and α, and given Y t-1 = 0, we have that α • Y t-1 is a r.v. degenerate at zero. As discussed by [START_REF] Al-Osh | First-order integer-valued autoregressive (INAR (1)) process[END_REF] and [START_REF] Du | The Integer Valued Autoregressive (INAR(p)) Model[END_REF], if α ∈ (0, 1), then {Y t } t∈Z is a stationary process, whereas α = 0 and α = 1 imply, respectively, independence and non-stationarity. Also, if 0 ≤ α < 1, the INAR(1) process is second-order stationary.

Several zero-inated models have been proposed in the literature to accommodate simultaneously both overdispersion and excess zeros: (i) the zeroinated Poisson (ZIP) model, [START_REF] Lambert | Zero-inated Poisson regression, with an application to defects in manufacturing[END_REF], in which the zero counts can come from two sources: from the Poisson distribution (sampling zeros) or from Bernoulli distribution (structural zeros); (ii) the zero-inated negative binomial distribution (ZINB), [START_REF] Hall | Zero-inated Poisson and Binomial regression with random eects: a case study[END_REF], [START_REF] Ridout | A score test for testing a zero-inated Poisson regression model against zero-inated negative binomial alternatives[END_REF], [START_REF] Garay | On estimation and inuence diagnostics for zero-inated negative binomial regression models[END_REF], which is obtained by mixing a Bernoulli distribution with a baseline negative binomial (NB) distribution over the ZIP model, among others.

In the following, we dene the zero-inated distributions, through their hierarchical formulation and introduce some further properties and three particular cases. Denition 2 The discrete r.v. V follows a zero-inated (ZI) distribution, with parameters ρ and λ, if it has the following stochastic representation:

V = BU, B⊥U,
where B is a Bernoulli r.v., with P (B = 1) = 1 -ρ and 0 ≤ ρ < 1. U is a nonnegative discrete r.v., with probability mass function (pmf ) h U (u|λ). λ can be a scalar or vector parameter, indexing the distribution of U . B⊥U indicates that the r.v. B and U are independent.

As consequence of Denition 2, we can obtain the pmf of V , given by:

P (V = v) = ρ + (1 -ρ) h U (0|λ) v = 0 (1 -ρ) h U (v|λ) v ≥ 1, (3)
where h U (v|λ) = P (U = v). We denote V ∼ ZI(ρ, λ; h U (•)).

E [V ] = (1 -ρ) E [U ] and V ar [V ] = (1 -ρ) V ar [U ] + ρE 2 [U ] . (4)
The distribution of the r.v. U determines the form of the ZI distribution. Thus, we describe briey three particular cases of the exible ZI models:

The zero-inated Poisson (ZIP) model:

In this case, we assume that the r.v. U follows a Poisson distribution, with mean λ. Thus, a pmf of the r.v. V , dened in Equation ( 3), takes the form:

P (V = v) = ρ + (1 -ρ) e -λ , v = 0 (1 -ρ) e -λ λ v v! v ≥ 1. (5)
We use the notation V ∼ ZIP (ρ, λ) so that:

E [V ] = (1 -ρ) λ and V ar [V ] = (1 -ρ) λ (1 + ρλ) .
The zero-inated negative binomial (ZINB) model:

The ZINB distribution is the result of considering that the r.v. U follows a negative binomial distribution, as in Equation ( 3). We denote by V ∼ ZINB (ρ, µ, φ) and its pmf is given by:

P (V = v) =          ρ + (1 -ρ) φ µ+φ φ , v = 0 (1 -ρ) Γ (φ + v) Γ (v + 1)Γ (φ) µ µ + φ v φ µ + φ φ , v ≥ 1
where µ ≥ 0 and φ > 0 is the dispersion parameter. Γ (•) represents the gamma function. We have that:

E [V ] = (1 -ρ)µ and V ar [V ] = (1 -ρ) µ 1 + µ φ + ρµ .
When ρ = 0, the r.v. V follows a negative binomial distribution, with mean µ and dispersion parameter φ, denoted by V ∼ NB(µ, φ). For more details, see [START_REF] Hall | Zero-inated Poisson and Binomial regression with random eects: a case study[END_REF], [START_REF] Ridout | A score test for testing a zero-inated Poisson regression model against zero-inated negative binomial alternatives[END_REF], [START_REF] Garay | On estimation and inuence diagnostics for zero-inated negative binomial regression models[END_REF], [START_REF] Martinez-Flores | Asymmetric regression models with limited responses with an application to antibody response to vaccine[END_REF], [START_REF] Garay | Bayesian estimation and case inuence diagnostics for the zero-inated Negative Binomial regression model[END_REF].

The zero-inated Poisson inverse Gaussian (ZIPIG) model:

This case arises when we consider that U , given in Equation ( 3), follows a Poisson inverse Gaussian (PIG) distribution. The pmf of the r.v. V , denoted by V ∼ ZIPIG(ρ, λ, φ), is:

P (V = v) =          ρ + (1 -ρ)e φ- √ φ(φ+2λ) , v = 0; (1 -ρ) 2 ρ [φ(φ + 2λ)] -(v-1/2) 2 e φ (λφ) v v! K v-1/2 φ(φ + 2λ) , v ≥ 1 . . .
where, µ > 0 is the mean, φ is a dispersion parameter and

K λ (t) = 1 2 ∞ 0 u λ-1 e -t 2 (u+ 1 u
) du is the modied Bessel function of the third kind, [START_REF] Abramowitz | Handbook of mathematical functions: With formulas, graphs, and mathematical tables applied mathematics series[END_REF]. When ρ = 0, the r.v. V follows a PIG distribution, with mean λ and dispersion parameter φ, denoted by V ∼ PIG(µ, φ), and we have that:

E [V ] = (1 -ρ)µ and V ar [V ] = (1 -ρ) µ 1 + µ φ + ρµ .
The ZIPIG distribution is a particular case of the mixed Poisson (MP) distribution, with hierarchical representation:

V |Z = z ∼ Poisson(µz),
where Z follows an Inverse Gaussian (IG) distribution, with mean 1 and dispersion parameter φ, denoted by Z ∼ IG(1, φ). More details are presented and discussed by [START_REF] Karlis | A general em approach for maximum likelihood estimation in mixed poisson regression models[END_REF] and [START_REF] Barreto-Souza | General mixed poisson regression models with varying dispersion[END_REF].

This manuscript is organized as follows. Section 2 outlines the proposed ZI-INAR(1) model and discusses some mathematical properties, including the likelihood function. Section 3 presents the implementation of an EM-type algorithm to estimate the parameters of the model and proposes several bootstrap methods to construct condence intervals for the parameters. The suitability and applicability of the process are illustrated in Sections 4 and 5, through simulation studies and analysis of real dataset, respectively. Finally, Section 6 concludes with short remarks and some possible avenues for future research.

The ZI-INAR(1) process

The ZI-INAR(1) process is an integer-valued rst-order autoregressive process, with ZI innovations, as presented in Equation (2), given by:

Y t = α • Y t-1 + V t , t ∈ Z, (6)
where V t ∼ ZI(ρ, λ; h U (•)).

Mathematical properties

In this subsection, we present some mathematical and structural properties of the corresponding marginal distributions of the process:

Proposition 1 Let {Y t } t∈Z be a stationary ZI-INAR(1) process; then:

E[Y t ] = (1 -ρ)E[U t ] 1 -α and V ar[Y t ] = (1 -ρ) αE[U t ] + ρE[U t ] 2 + V ar[U t ] 1 -α 2 ,
where U t denotes a r.v. with density function (or pmf) h U (u|λ), for all t ∈ Z.

Proof By using the expectation obtained in Eq. ( 4), we have that:

E[Y t |Y t-1 ] = E[α • Y t-1 + V t |Y t-1 ] = E[α • Y t-1 |Y t-1 ] + E[V t |Y t-1 ] = αY t-1 + (1 -ρ) E [U t ] . (7)
Thus,

E[Y t ] = E [E (Y t |Y t-1 )] = E [αY t-1 + (1 -ρ) E (U t )] = αE[Y t-1 ] + (1 -ρ) E [U t ] . (8)
Eq. ( 8) is a consequence of the fact that the process is stationary, that is,

E[Y t ] = µ Y , for all t ∈ Z.
To obtain the variance of the process, we use:

V ar[Y t |Y t-1 = y t-1 ] = V ar[α • Y t-1 + V t |Y t-1 = y t-1 ] = V ar[α • Y t-1 |Y t-1 = y t-1 ] + V ar[V t |Y t-1 = y t-1 ] = α(1 -α)y t-1 + (1 -ρ) V ar [U t ] + ρE 2 [U t ] . (9)
Thus, by [START_REF] Politis | The stationary bootstrap[END_REF] and ( 9) we obtain:

V ar[Y t ] = V ar[E(Y t |Y t-1 )] + E[V ar(Y t |Y t-1 )] = α 2 V ar[Y t-1 ] + α(1 -α)E[Y t-1 ] + (1 -ρ)(V ar[U t ] + ρE 2 [U t ]) = α 2 V ar[Y t-1 ] + α(1 -ρ)E[U t ] + (1 -ρ)(V ar[U t ] + ρE 2 [U t ]), (10)
where [START_REF] Garay | On estimation and inuence diagnostics for zero-inated negative binomial regression models[END_REF] is obtained by the expectation given in Eq. ( 8).

The likelihood function

Let y = (y 1 , y 2 , . . . , y n ) be realizations of the ZI-INAR(1) process. Then the likelihood function of the unknown parameters θ = (α, ρ, λ) , given y, can be written as:

L (θ|y) = P (Y 1 = y 1 , Y 2 = y 2 , . . . , Y n = y n ) = P (Y 1 = y 1 ) P (Y 2 = y 2 |Y 1 = y 1 ) . . . P (Y n = y n |Y 1 = y 1 , . . . , Y n-1 = y n-1 ) = P (Y 1 = y 1 ) n t=2 P (Y t = y t |Y t-1 = y t-1 ) , (11)
where

P (Y t = y t |Y t-1 = y t-1 ) = min{yt-1,yt} k=0 y t-1 k α k (1 -α) yt-1-k × ρI {0} (y t -k) + (1 -ρ)h (y t -k|λ) . ( 12
)
I A (•) denotes the indicator function, i.e., I A (y) = 1, if y ∈ A and I A (y) = 0 otherwise.

Note that Equation [START_REF] Hall | Zero-inated Poisson and Binomial regression with random eects: a case study[END_REF] represents the probability transition of a stationary Markov chain, of state y t-1 to y t . Thus, the marginal probability function can be dened by:

P (Y t = y t ) = ∞ yt-1=0 P (Y t-1 = y t-1 ) P (Y t = y t |Y t-1 = y t-1 ) . ( 13
)
From Equations (11)(13), we have that:

L (θ|y) = P (Y 1 = y 1 ) n t=2 min{yt-1,yt} k=0 y t-1 k α k (1 -α) yt-1-k × ρI {0} (y t -k) + (1 -ρ)h (y t -k|λ) . (14)
The marginal distribution is intractable, so a simple approach to deal with this is to condition it on the observed Y 1 , and estimate the parameters by conditional maximum likelihood (CML); see [START_REF] Al-Osh | First-order integer-valued autoregressive (INAR (1)) process[END_REF] and [START_REF] Jazi | First-order integer valued AR processes with zero inated Poisson innovations[END_REF]. To estimate the parameters of the ZI-INAR (1) process by maximizing this conditional likelihood function directly prevents the possibility of analytical solutions. One alternative is to maximize the complete conditional likelihood using the expectationmaximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF], which is stable and straightforward to implement, since the iterations converge monotonically and no second derivatives are required. In the next section, we discuss a technique to nd the ML estimates of the parameters vector θ, based on the EM algorithm.

3 Maximum likelihood estimation and bootstrap resampling methods

Parameter estimation via the EM algorithm

In this section we develop an EM-type algorithm for maximum likelihood estimation of the parameters of the ZI-INAR(1) process. The key to the development of our EM-type algorithm is to consider the presence of latent variables and treat the problem as if these variables were in fact observed.

As suggested by [START_REF] Hall | Zero-inated Poisson and Binomial regression with random eects: a case study[END_REF] and [START_REF] Neal | MCMC for integer-valued ARMA processes[END_REF], we dene the latent variables W = (W 1 , . . . , W n ) and S = (S 1 , . . . , S n ) for all t ≥ 1, where:

S t is dened by S t = α • Y t-1 . Thus S t |Y t-1 = y t-1 , α ∼ Bin (y t-1 ; α i ) , if y t-1 >
0 and it is a degenerate distribution at zero if y t-1 = 0, where Bin (y; α) represents the binomial distribution with parameters y and α ∈

[0, 1].
W t is a latent dichotomous variable, that is:

W t = 1 if V t is from the zero state 0 if V t ∼ h (•|λ), with P (W t = 1) = ρ and h (•|λ) as presented in Denition 2.
Let Y c = (Y, W, S) be the complete data vector, where Y and {W, S} represent the observed data and the missing data, respectively, with Y = (Y 1 , . . . , Y n ) . Then the joint probability function for Y ct = (Y t , S t , W t ) is given by:

P (Y ct = y ct ) = P (Y t = y t , S t = s t , W t = w t |Y t-1 = y t-1 ) = P (Y t = y t , W t = w t |S t = s t , Y t-1 = y t-1 ) P (S t = s t |Y t-1 = y t-1 ) = ρ wt ((1 -ρ) h(y t -s t |λ)) 1-wt × y t-1 s t α st (1 -α) yt-1-st , for all t ≥ 1.
Thus, the complete likelihood function is dened by:

L c (θ|y c ) = P (Y 1 = y 1 ) n t=2 P (Y t = Y t , S t = s t , W t = w t |Y t-1 = y t-1 ) ∝ n t=2 ρ wt ((1 -ρ) h (y t -s t |λ)) 1-wt y t-1 s t α st (1 -α) yt-1-st .
Hence, the complete log-likelihood function is given by:

c (θ|y c ) = log (L c (θ|y c )) ∝ n t=2 w t log (ρ) + n t=2 (1 -w t ) log (1 -ρ) + n t=2 (1 -w t ) log (h (y t -s t |λ)) + n t=2 s t log (α) + n t=2 (y t-1 -s t ) log (1 -α) .
The EM algorithm has several appealing properties relative to other iterative algorithms such as Newton-Raphson and Fisher's scoring method for nding MLEs, [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]. One of them is that the sequence of estimates from the EM algorithm increases the likelihood function (θ|y) at each iteration, and under standard regularity conditions the sequence converges to a stationary point of the likelihood.

Thus, the EM-algorithm proceeds in two steps:

E-step: Let θ (k)
be the current k-th step estimate of θ. By using the property of conditional expectation, we compute the Q(θ| θ

(k)
) function

given by:

Q(θ| θ (k) ) = E c (θ|y c )|y, θ (k)
.

(

) M-step: Maximize Q(θ| θ (k) 15
) with respect to θ, obtaining θ

.

Observe that the expression of the Q-function, in Eq. ( 15), is determined by the knowledge of the following expectations:

s (k) t = E S t |y, θ (k)
w

(k) t = E W t |y, θ (k)
and

Q * t (λ|θ (k) ) = E (1 -W t ) log h(y t -S t |λ)|y, θ (k)
.

Thus, the Q-function can be written in more simple form as:

Q(θ| θ (k) ) ∝ n t=2 w (k) t log (ρ) + n t=2 1 -w (k) t log (1 -ρ) + n t=2 Q * t (λ|θ (k) ) + n t=2 s (k) t log (α) + n t=2 y t-1 -s (k) t log (1 -α) .
At each step, s

(k) t can be obtained by using the following results:

P S t = s t |Y t-1 = y t-1 , θ (k) = P S t = s t , Y t = y t |Y t-1 = y t-1 , θ (k)
P Y t = y t |Y t-1 = y t-1 , θ (k) = P Y t = y t |S t = s t , Y t-1 = y t-1 , θ (k)
× P S t = s t |Y t-1 = y t-1 , θ (k) min{yt-1,yt} st=0 P Y t = y t , S t = s t |Y t-1 = y t-1 , θ (k) = y t-1 s t (α (k) ) st (1 -α (k) ) yt-1-st ρ (k) I {0} (y t -s t ) + (1 -ρ (k) )h(y t -s t |λ (k) ) min{yt-1,yt} st=0 P Y t = y t , S t = s t |Y t-1 = y t-1 , θ (k) (16)
and

P W t = w t |Y t-1 = y t-1 , θ (k) = P W t = w t , Y t = y t |Y t-1 = y t-1 , θ (k) P Y t = t t |Y t-1 = y t-1 , θ (k) = P W t = w t , Y t = y t |Y t-1 = y t-1 , θ (k) min{yt-1,yt} st=0 P Y t = y t , S t = s t |Y t-1 = y t-1 , θ (k)
.

Then, for all t ≥ 1, we have that:

s (k) t = min{yt-1,yt} st=0
s t P S t = s t |Y t-1 = y t-1 , θ (k) (17)
w

(k) t = 1 × P W t = 1|y, θ (k)
=

P W t = 1, Y t = y t |Y t-1 = y t-1 , θ (k) min{yt-1,yt} st=0 P Y t = y t , S t = s t |Y t-1 = y t-1 , θ (k)
.

Therefore, our EM algorithm for the ZI-INAR(1) process can be summarized in the following way:

E-step: Given θ (k) , for t ≥ 1 we compute s (k) t , w (k)
t , as given in Eq. ( 17)- [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF], respectively, and

Q * t (λ|θ (k) ). M-step: Update θ (k) by maximizing Q(θ| θ (k)
) over θ, which leads to the following expressions:

α (k+1) = n t=2 s t n t=2 y t-1 , ρ (k+1) = n t=2 w t n -1 , λ (k+1) = arg max λ n t=2 Q * t (λ|θ (k) ) . (19)
In the following, we develop the procedure to obtain the expressions

Q * t (λ|θ (k) ) and λ (k+1)
, considering the three particular cases of the ZI models seen before, i.e., when V ∼ ZI(ρ, λ). For this, we dene the expectation b t s t with B t = 1 -W t , given by:

b t s t = E B t S t |y, θ (k) = 1 -ρ (k) min{yt-1,yt} st=0 s t y t-1 s t α (k) st 1 -α (k) yt-1-st h y t -s t | λ (k) P Y t = y t |Y t-1 = y t-1 , θ (k)
.

Thus,

-If V t ∼ ZIP(ρ,λ), then from Eq. ( 5) h(•|λ) represents a pmf of the Poisson distribution, with parameter λ. Consequently,

Q * t (λ|θ (k) ) ∝ -λ(1 -w (k) t ) + log(λ)(1 -w (k) t )y t + log(λ) b t s t (k) (20)
and from Eq. ( 19) and ( 20):

λ (k+1) = n t=2 (1 -w (k) t )y t - n t=2 b t s t (k) n t=2 (1 -w (k) t )
.

-If V t ∼ ZINB(ρ,µ,φ), then h(•|λ) represents the pmf of the negative binomial distribution with parameters λ = (µ, φ). Thus,

Q * t (λ|θ (k) ) ∝ g t (φ) (k) + log(µ) -log(µ + φ) (1 -w (k) t )y t -b t s t (k) + -log (Γ (φ)) + φ (log(φ) -log(µ + φ)) 1 -w (k) t (21)
with

g t (φ) (k) = E B t log Γ (y t -S t + φ)|y, θ (k) = 1 -ρ (k) min{yt-1,yt} st=0 log Γ (y t -s t + φ)P(S t = s t |Y t-1 = y t-1 , θ (k)
)h(y t -s t | θ

)

P Y t = y t |Y t-1 = y t-1 , θ (k)
.

(

) 22
When the M-step turns out to be analytically intractable, it can be replaced by a sequence of conditional maximization (CM) steps. The procedure is known as the ECM algorithm [START_REF] Meng | Maximum likelihood estimation via the ECM algorithm: A general framework[END_REF]. Thus, from Eq. ( 19), ( 21) and ( 22), we have that for the ZINB-INAR(1) process, λ

= ( µ (k+1) , φ (k+1) ) are given by:

µ (k+1) = n t=2 (1 -w (k) t )y t - n t=2 b t s t (k) n t=2 (1 -w (k) t )
.

φ (k+1) = arg max φ n t=2 Q * t (µ (k+1) , φ| θ (k) )
φ (k+1) is obtained using the optim routine in the R software ( [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]).

-If V t ∼ ZIPIG(ρ,µ,φ), then λ = (µ, φ) and Its hierarchical representation is:

V t |Z t = z t , W t = 0 ∼ Poisson(µz t ) V t |Z t = z t , W t = 1 follows a degenerate at zero; (23)
Z t |W t = 0 ∼ IG(1, φ), (24)
W t ∼ Bin(1, ρ), t = 2, 3, . . . , (25)
the density function of Z t |W t = w t is:

f (z t |w t ) = φ 1/2 z -3/2 t exp - φ 2 (z t -1) 2 z t 1-wt . (26)
In order to develop the EM algorithm for the ZIPIG-INAR model, we add the latent variable Z t . Thus, the complete-data is dened by Y c = (Y, W, S, Z). Thus,

Q * t (µ, φ|θ (k) ) ∝ -µ b t z t (k) -log(µ) b t s t (k) + log(µ)(1 -w (k) t )y t + log(φ) 2 (1 -w (k) t ) + φ(1 -w (k) t ) - φ 2 b t z t (k) + b t /z t (k) , (27)
where the expectations b t z t (k)

and b t /z t (k)
are dened by:

b t z t (k) = E E B t Z t |S t , y, θ (k) = min{yt-1,yt} st=0 1 -ρ (k) (y t -s t + 1)h y t -s t + 1| λ (k) × P S t = s t |Y t = y t , Y t-1 = y t-1 , θ (k)
× 1

µ (k) ρ (k) I {yt=st} + 1 -ρ (k) h y t -s t | λ (k) (28)
and

b t /z t (k) = E E B t Z -1 t |S t , , y, θ (k) = min{yt-1,yt} st=0 µ (k) h(y t -s t -1| λ (k) )P(S t = s t |Y t = y t , Y t-1 = y t-1 , θ (k)
)

(y t -s t )h(y t -s t | λ (k)
)

I {yt>st} + ρ (k) ( ρ (k) + 2 µ (k) ) + 1 1 -ρ (k) h(0| λ (k)
)

× P S = y t |Y t = y t , Y t-1 = y t-1 , θ (k)
× 1

ρ (k) ρ (k) + 1 -ρ (k) h 0| λ (k) I {yt=st} , (29)
where

P S t = s t |Y t = y t , Y t-1 = y t-1 , θ (k)
is given in Eq. ( 16). Thus, from Eq. ( 19) and ( 27), we have that:

µ (k+1) = n t=2 (1 -w (k) t )y t - n t=2 b t s t (k) n t=2 b t z t . φ (k+1) = n t=2 (1 -w (k) t ) n t=2 b t z t (k) + n t=2 b t /z t (k) -2 n t=2 (1 -w (k) t )
.

This process is iterated until some convergence rule is satised. Here, we use the Aitken acceleration-based stopping criterion [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF] as a convergence rule.

This criterion is based on the fact that the limit of the sequence (k+1) = ( θ (k+1) |y), denoted by ∞ , can be approximated by

(k+1) ∞ = (k) + ( (k+1) - (k) )/(1 -c (k) ), where c (k) = ( (k+1) -(k) )/( (k) -(k-1)
). As suggested by [START_REF] Zeller | Finite mixture of regression models for censored data based on scale mixtures of normal distributions[END_REF], we decided to stop the algorithm when (k+1) ∞ -(k+1) < ε = 10 -5 .

Bootstrap resampling methods

In this section, we propose several bootstrap alternatives to construct condence intervals for the parameter. Actually, we can indierently use the posterior mean of the Bayesian approach or the standard maximum likelihood estimators (mle), or any other ecient estimators. In any case, the standard error can be dicult to compute, so several bootstrap approaches can be used to solve this issue. Here we propose two types of approaches, based respectively on some strong parametric assumptions, on the mixing properties of the INAR(p) process and on the Markov property of this process.

(a) Parametric bootstrap approach:

Let θ = (α, ρ, λ) be the parameters of the ZI-INAR(1) process, where we denote θ n = ( α n , ρ n , λ n ) as the mle estimator obtained using the EM algorithm developed in Section 3.1. The properties of mle for general INAR(p) models are studied, for instance, in [START_REF] Du | The Integer Valued Autoregressive (INAR(p)) Model[END_REF], who showed strong consistency of mle, and [START_REF] Latour | Existence and Stochastic Structure of a Non Negative Integer-Valued Autoregressive Process[END_REF] who studied the eciency and asymptotic normality in Proposition 6.1. Other estimators, based on saddle point approximation of the likelihood, which may be easier to implement, were proposed in [START_REF] Pedeli | Likelihood Estimation for the INAR(p) Model by Saddlepoint Approximation[END_REF]. In the following, we denote by θ n any kind of asymptotically ecient estimator of the parameter θ.

The parametric bootstrap method simply consists of generating new data in the model with estimated parameters. That is:

Step 1: Generate the random variable V * t from the ZI distribution, with pa- rameters θ n .

Step 2: Generate a ZI-INAR(1) model recursively as follows:

Y * 1 = y 1 Y * t = α n • Y * t-1 + V * t , t = 2, ..., n, meaning that if Y * t-1 > 0, then α n • Y * t-1 |Y * t-1 d = Y * t-1 i=1 Z * i with Z * i ∼ Bin(1, α n ) or α n • Y * t-1 |Y * t-1 is degenerate at zero if Y * t-1 = 0.
Notice that the rst observation of the bootstrap ZI-INAR(1) process is set to the rst observation of the observed data: asymptotically, even if the generated process is not stationary, this will not perturb the asymptotic properties of the bootstrap process. Now the bootstrap counter- The centered and normalized bootstrap distribution is then given by, for u = (u 1 , . . . , u 4 ) ∈ R 4 :

parts of θ n = ( α n , ρ n , λ n ), say θ * n = ( α * n , ρ * n , λ * n ),
K * (B) n (u) = 1 B B b=1 I { √ n( θ * (b) n -θn)≤u} .
If the bootstrap distribution is asymptotically valid, that is, it is a convergent distribution (at least in probability) of P(

√ n( θ n -θ) ≤ u) then
it is easy to use the quantile of the Monte Carlo bootstrap distribution to construct (simultaneous) condence intervals for the parameter θ or any regular (at least dierentiable) functional of the parameter.

Concerning the choice of B, it is known from the work by [START_REF] Hall | On the Number of Bootstrap Simulations Required to Construct a Condence Interval[END_REF] that one needs to choose B so that it at least has the same order as n and such that (B + 1)(1 -γ)/2 is an integer (when constructing asymptotic two-sided condence intervals of level γ. So typically for γ = 95% (for size of n smaller than 1000) we choose B = 999.

Proposition 2 Assume that the parameter space of θ = (α, ρ, with λ = (µ, φ) is given by

Θ =]0, 1[×]0, 1[×]0, ∞[×]0, ∞[
and consider θ n to be the mle of the ZI-INAR(1) process. Then the parametric bootstrap in asymptotically correct, meaning that almost surely along the sample, notice that any ecient estimator has the LAE property, so that one may use indierently the mle in the construction or the estimators proposed in [START_REF] Pedeli | Likelihood Estimation for the INAR(p) Model by Saddlepoint Approximation[END_REF].

sup u∈R 4 P( √ n( θ * n -θ n ) ≤ u) -P( √ n( θ n -θ) ≤ u) → 0 when n → ∞ Proof: In [5], theorem

(b) Moving block bootstrap approach

When no specic assumption is made about the distribution of the residuals, it is still possible to implement semiparametric estimators of the ZI-INAR process, as described in [START_REF] Pedeli | Likelihood Estimation for the INAR(p) Model by Saddlepoint Approximation[END_REF]. In that case, a parametric bootstrap process cannot be used. Moreover, to assess the robustness of the parametric assumptions, it may also be interesting to implement a more robust version of the bootstrap process. The more general method is based on splitting the original time series into overlapping blocks which are then resampled to reconstruct the original time series. The procedure in our case is as follows:

Step 1: Choose a length b (which will typically be of size b = o ( √ n)) . Dene the overlapping blocks

B 1 = (Y 1 , ..., Y b ), B 2 = (Y 2 , ..., Y b+1 ), ...., B n-b+1 = (Y n-b+1 , ..., Y n )
Use circular block bootstrap, or even better, apply stationary block bootstrap (see [START_REF] Politis | The stationary bootstrap[END_REF]), which allows simulating a stationary version of the moving block bootstrap, in place of the moving block bootstrap method.

Step 2: Draw without replacement [ n b ] + 1 blocks, which are bound together (and possible truncated at the end) to form a new time series of size n.

After this, compute the statistics of interest, either the mle or estimators based on approximations (see [START_REF] Pedeli | Likelihood Estimation for the INAR(p) Model by Saddlepoint Approximation[END_REF]) of the corresponding time series.

Step 3: Just like in the Step 3 of the parametric bootstrap, use Monte Carlo method to obtain an approximation of the bootstrap distribution.

It has been shown in [START_REF] Kunsch | The jackknife and the bootstrap for general stationary observations[END_REF] that under some strong mixing conditions and provided that b n → ∞ and bn √ n → 0, the moving block bootstrap is asymptotically valid. A lot of variation and modications have been proposed in the literature to obtain valid second-order approximations (see for instance [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF] for a complete overview and references). However, it should be noticed that according to the distribution of V t , the INAR(p) process may not strongly mix. The process can be shown to strongly mix when the residuals V t have Poisson distribution. However, it can be shown that the process is only weakly mixing, or ψ -weak mixing for more general distributions; see [START_REF] Doukhan | On weak dependence conditions: The case of discrete valued processes[END_REF]. This includes the ZI cases studied here, when 0 < α < 1, since the moments of all residuals are nite. The validity of the stationary moving block bootstrap method (as well as variations including the circular block bootstrap) in the weak mixing case has been studied in [START_REF] Hwang | Strong consistency of the stationary bootstrap under ψ-weak dependence[END_REF], see theorem 3.2. As noticed by the authors, their theorem can be easily extended to functionals of means, Frechet dierentiable functionals (which can be simply linearized) or even Hadamard dierentiable functionals (including M-estimators). As a consequence, this method will be asymptotically valid for the mle, which is LAE (it is less obvious for the estimators proposed in [START_REF] Pedeli | Likelihood Estimation for the INAR(p) Model by Saddlepoint Approximation[END_REF]).

Condence intervals

In the second simulation study, we use and compare two dierent types of asymptotically valid condence intervals, mainly bootstrap variance based asymptotic intervals and the standard percentile method. We briey describe these methods in our context. Consider θ k the components of the mle of θ = (θ 1 , . . . , θ 4 ) = (α, ρ, µ, φ) and denote by θ * (b) k,l , b = 1, ..., B, the corresponding Bootstrap replications obtained at the Monte-Carlo step, either by the parametric (l = 1) or moving block (l = 2) techniques. In the simulation we will choose B such that (B + 1)(1 -γ)/2 is an integer, as suggested by [START_REF] Hall | On the Number of Bootstrap Simulations Required to Construct a Condence Interval[END_REF].

Variance based asymptotic interval

The bootstrap variance of θ k is simply given by:

V * k,l = 1 B B b=1 θ * (b) k,l - 1 B B i=1 θ * (i) k,l 2
. Tables 4 and5 show the mean value of the bootstrap estimator and the proportion of times that the condence intervals contain the true value of the

Conclusion

This article presents a study of the ZI-INAR(1) processes, which are an extension of the ZINAR(1) process, introduced by [START_REF] Jazi | First-order integer valued AR processes with zero inated Poisson innovations[END_REF]. The aim is to propose models with more exibility and that can deal with excess zeros. To t the models, we develop an ecient EM algorithm to estimate the parameters and apply several alternative bootstrap techniques to construct condence intervals for the parameters. The performance of our proposals was evaluated through simulation studies and application to a real dataset of monthly number of drug oenses recorded. In this illustrative example, there is strong evidence of a rst-order dependence structure and zero-inated counts. We conjecture that our methods can be extended by adding a moving average structure to the process. An in-depth investigation of such extension is beyond the scope of the present manuscript, but is an interesting topic for further research.

  are the mle obtained using our proposed EM algorithm, considering the bootstrap process Y * = (Y * 1 , . . . , Y * n ). Step 3: Since the exact bootstrap distribution of θ * n may be dicult (and time consuming) to compute, it is replaced by a Monte Carlo approximation that is repeats the procedure in Step 2, B times, by generating times series y * (b) , b=1,. . . , B , for realization of the bootstrap process Y * . Then we can compute the corresponding θ * (b) n , b = 1, . . . , B estimators.

4. 2 3 .

 23 Simulation 2: Bootstrap condence intervals In this simulation study we analyze the eciency of the bootstrap methods presented in Section 3. We simulate N = 100 samples of size n = 196 from the following processes with parameters α = 0.25 and ρ = 0.4: (i) ZIP-INAR(1), with λ = 4; and (ii) ZINB-INAR(1) and ZIPIG-INAR(1), with µ = 4 and φ = For each sample, we generate B = 999 bootstrap replicates and the length of the circular block bootstrap adopted is b = 14.

Figure 1

 1 Figure 1 contains the panels (a) and (b), which represent the time evolution and the bar plot of the monthly count of drug oenses recorded, respectively. The two plots indicate a signicant presence of zero value observations (43%). In panel (c), the partial autocorrelation function (PACF) suggests a rst-order dependence structure. For comparative purposes, we consider the INAR(1) process with Poisson (Po), NB and PIG innovations and ZI-INAR(1) processes.Table6presents the parameter estimates and the Akaike information criterion (AIC) values,[START_REF] Akaike | A new look at the statistical model identication[END_REF], of all the considered processes. This criterion indicates that the ZIPIG-INAR(1) process is more appropriate for this dataset. Table 7 shows the parameter estimates of the ZIPIG-INAR(1) process and the condence interval, using parametric bootstrap (B-PAR) and circular block bootstrap (CBB).

Fig. 1 :

 1 Fig. 1: Drugs oence dataset. (a) Time series plot. (b) Bar plot. (c) PACF.

Table 1 :

 1 RB and RRSE of the parameter estimates of ZIP-INAR(1) processes, with α = 0.3 and λ = 2.

	ρ	n	RB	α	RRSE	RB	ρ	RRSE	RB	λ	RRSE
		100	-0.0668 0.3141 -0.0764 0.3436 -0.0035 0.1163
	0.3	300 500	-0.0166 0.1635 -0.0090 0.1943 -0.0229 0.1188 -0.0128 0.1372	0.0058 0.0033		0.0777 0.0579
		1000 -0.0020 0.0870 -0.0150 0.0995 -0.0024 0.0395
		100	-0.0433 0.2362 -0.0212 0.1337	0.0007		0.1582
	0.6	300 500	-0.0118 0.1289 -0.0068 0.0707 -0.0058 0.0895 -0.0117 0.0985 -0.0050 0.0488 0.0021 0.0650
		1000 -0.0095 0.0750 -0.0031 0.0388 -0.0040 0.0472

Table 2 :

 2 RB and RRSE of the parameter estimates of ZINB-INAR(1) processes, with α = 0.3 and µ = 2.

				α		ρ		µ		φ
	φ	ρ	n	RB	RRSE	RB	RRSE	RB	RRSE	RB	RRSE
			100	0.0051	0.2080	-0.0001	0.7010	0.0971	0.3614	0.3417	1.8875
		0.3	300 500	0.0042 0.0022	0.1117 0.0842	-0.0365 -0.0238	0.5423 0.4571	0.0344 0.0143	0.2355 0.1873	0.1165 0.0696	0.6861 0.4845
	0.75		1000 100	0.0019 0.0052	0.0595 0.1893	-0.0424 -0.1443	0.3500 0.4576	-0.0048 0.0306	0.1376 0.4736	0.0111 0.2758	0.2943 3.2307
		0.6	300 500	-0.0148 0.0083	0.1091 0.0830	-0.0688 -0.0684	0.2908 0.2552	0.0106 -0.0158	0.2980 0.2584	0.0990 0.0125	0.9289 0.6428
			1000	-0.0021	0.0612	-0.0242	0.1472	-0.0076	0.1692	0.0602	0.3940
			100	-0.0157	0.2241	-0.0666	0.6120	0.0183	0.2399	0.2144	1.4996
		0.3	300 500	-0.0015 -0.0019	0.1219 0.0899	-0.0834 -0.0185	0.4179 0.3245	-0.0147 0.0084	0.1604 0.1294	0.0298 0.0585	0.6910 0.4348
	1.5		1000 100	0.0020 -0.0200	0.0687 0.2105	-0.0333 -0.1236	0.2436 0.3800	-0.0043 -0.0118	0.0950 0.3557	-0.0039 0.0959	0.2898 2.1109
		0.6	300 500	0.0005 -0.0109	0.1204 0.0949	-0.0407 -0.0256	0.1930 0.1351	-0.0130 -0.0167	0.2056 0.1607	0.0401 0.0213	1.1545 0.6549
			1000	-0.0004	0.0603	-0.0218	0.0883	-0.0184	0.1161	-0.0314	0.3832
			100	-0.0388	0.2513	-0.1162	0.5348	0.0030	0.2104	0.0319	0.9771
		0.3	300 500	-0.0050 -0.0212	0.1370 0.1089	-0.0642 -0.0502	0.3803 0.2699	-0.0077 -0.0061	0.1390 0.1037	0.1027 0.0178	0.9715 0.5092
	2.5		1000 100	-0.0066 -0.0169	0.0764 0.2278	-0.0281 -0.1465	0.1829 0.3334	-0.0082 -0.0959	0.0755 0.2975	-0.0169 -0.1273	0.3569 1.3087
		0.6	300 500	-0.0125 -0.0030	0.1194 0.0889	-0.0165 -0.0088	0.1262 0.0954	-0.0176 -0.0022	0.1730 0.1286	0.0525 0.0730	1.0700 0.8590
			1000	-0.0006	0.0697	-0.0065	0.0690	-0.0025	0.0977	< 0.001	0.4792

Table 3 :

 3 RB and RRSE of the parameter estimates of ZIPIG-INAR(1) processes, with α = 0.3 and µ = 2.

				α		ρ		µ		φ
	φ	ρ	n	RB	RRSE	RB	RRSE	RB	RRSE	RB	RRSE
			100	0.0010	0.2044	-0.0046	0.5836	0.0391	0.2632	0.1688	1.8383
		0.3	300 500	-0.0056 -0.0002	0.1175 0.0887	-0.0665 -0.0243	0.3940 0.2975	-0.0078 0.0004	0.1739 0.1293	0.0431 0.0237	0.5878 0.3576
	0.75		1000 100	-0.0016 -0.0205	0.0614 0.2034	-0.0117 -0.0825	0.2135 0.3159	-0.0004 0.0369	0.0918 0.3589	0.0301 0.2246	0.2578 2.6273
		0.6	300 500	-0.0007 -0.0073	0.1125 0.0871	-0.0345 -0.0121	0.1567 0.1040	-0.0168 -0.0096	0.2096 0.1539	-0.0025 0.0717	0.7321 0.5155
			1000	-0.0027	0.0677	-0.0095	0.0790	-0.0110	0.1204	0.0187	0.3109
			100	-0.0340	0.2193	-0.0740	0.5251	0.0183	0.2156	0.1622	1.3835
		0.3	300 500	-0.0068 -0.0050	0.1367 0.1032	-0.0243 -0.0162	0.3257 0.2326	0.0016 0.0007	0.1346 0.0959	0.0398 0.0268	0.7542 0.4768
	1.5		1000 100	-0.0065 -0.0320	0.0702 0.2046	-0.0273 -0.0452	0.1705 0.2319	-0.0001 0.0160	0.0762 0.2772	-0.0163 0.1455	0.2293 1.6193
		0.6	300 500	-0.0052 -0.0051	0.1066 0.0913	-0.0212 -0.0096	0.1103 0.0885	-0.0034 -0.0015	0.1659 0.1253	0.0623 0.0543	0.9166 0.9582
			1000	0.0022	0.0640	0.0025	0.0573	0.0051	0.0942	0.0481	0.3459
			100	-0.0437	0.2472	-0.0750	0.4556	0.0187	0.1978	0.1961	1.3687
		0.3	300 500	-0.0095 -0.0040	0.1382 0.1188	-0.0194 0.0002	0.2708 0.2095	-0.0006 0.0061	0.1141 0.0938	0.0972 0.0628	0.9488 0.4610
	2.5		1000 100	-0.0042 -0.0486	0.0720 0.2160	-0.0017 -0.0850	0.1384 0.2530	0.0070 -0.0336	0.0621 0.2658	0.0338 -0.0715	0.3133 1.1974
		0.6	300 500	-0.0107 -0.0088	0.1299 0.0918	-0.0161 -0.0143	0.0976 0.0777	-0.0046 -0.0097	0.1420 0.1051	0.0965 0.0423	0.9962 0.7257
			1000	0.0025	0.0653	-0.0031	0.0503	-0.0027	0.0776	0.0077	0.3785

parameter (Cov.), as well as the percentage of times that the parameter value is less than the lower bound (θ k < LL) and greater than the upper bound (U L < θ k ). In general, the proportion of the intervals containing the true value of the parameters is close to the nominal level (95%). However, these proportions are smaller than 90% for the intervals based on percentiles of the ZINB-INAR(1) process.

Table 4 :

 4 Parametric bootstrap approach results using condence intervals.

	Model ZIP-INAR(1) ZINB-INAR(1) ZIPIG-INAR(1)	α ρ µ α ρ µ φ α ρ µ φ	boots. value mean θ k < LL U L < θ k Cov. θ k < LL U L < θ k Cov. Asymptotic interval Standard percentile method 0.25 0.2504 3% 4% 93% 3% 6% 91% 0.40 0.4038 3% 1% 96% 3% 1% 96% 4.00 3.9864 1% 4% 95% 1% 4% 95% 0.25 0.2430 1% 2% 97% 0% 3% 97% 0.40 0.3772 3% 1% 96% 3% 3% 94% 4.00 3.9360 2% 6% 92% 2% 6% 92% 3.00 4.4043 0% 4% 96% 5% 2% 93% 0.25 0.2441 4% 3% 93% 4% 3% 93% 0.40 0.3863 0% 1% 99% 0% 2% 98% 4.00 4.0022 7% 1% 92% 7% 1% 92% 3.00 5.0738 0% 2% 98% 4% 0% 96%

Table 5 :

 5 Circular block bootstrap results using condence intervals.To further demostrate the usefulness of our approach, we analyze a real dataset of the monthly number of drug oenses recorded from January 1990 to December 2001, with 144 observations, in Pittsburgh's census tract 2206. The data were obtained from Forecasting Principles site: http://www.forecastingprinciples.com.

	Model ZIP-INAR(1) ZINB-INAR(1) ZIPIG-INAR(1)	α ρ µ α ρ µ φ α ρ µ φ	boots. value mean θ k < LL U L < θ k Cov. θ k < LL U L < θ k Cov. Asymptotic interval Standard percentile method 0.25 0.2375 2% 5% 93% 2% 5% 93% 0.40 0.3945 3% 3% 94% 3% 3% 94% 4.00 3.9912 2% 4% 94% 2% 5% 93% 0.25 0.2249 0% 9% 91% 0% 11% 89% 0.40 0.3575 2% 4% 94% 1% 9% 90% 4.00 3.9057 3% 9% 88% 3% 9% 88% 3.00 4.1189 0% 7% 93% 4% 7% 89% 0.25 0.2258 0% 6% 94% 0% 6% 94% 0.40 0.3693 0% 8% 92% 0% 9% 91% 4.00 3.9847 3% 1% 96% 5% 1% 94% 3.00 4.8409 0% 5% 95% 4% 3% 93%
	5 Real dataset: Drug oenses

Table 6 :

 6 Parameter estimates

	Processo	α	ρ	µ	φ	AIC
	Po-INAR(1)	0.212	-	1.679	-	590.74
	ZIP-INAR(1)	0.181 0.512 3.577	-	565.77
	NB-INAR(1)	0.071	-	1.977 0.471	550.43
	ZINB-INAR(1)	0.070 0.138 2.296 0.630	552.20
	PIG-INAR(1)	0.072	-	1.973 0.336	554.53
	ZIPIG-INAR(1) 0.065 0.325 2.946 0.903 549.41

Acknowledgements

Aldo M. Garay would like to acknowledge the support of the National Council for Scientic and Technological Development, CNPq-Brazil, Grant 420082/2016-6. This research has also been conducted as part of the project Labex MME-DII (ANR11-LBX-0023-01).

Thus, the asymptotic normal approximation leads to an asymptotic γ condence interval of the type

, where u 1+γ 2 is a quantile of order 1+γ 2 of a standard normal distribution.

Standard percentile method

Dene respectively θ * k,l(

2 ) .

Simulation Study

In this section we show the performance of the estimation procedure for the ZI-INAR(1) process, using the EM algorithm proposed and bootstrap resampling methods described before. Thus, we present examples considering articial datasets, in dierent scenarios. All the computational codes were implemented using the R software ( [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]) and the program codes are available from us on request.

We present two simulation studies. The rst one investigate the asymptotic properties of the maximum likelihood estimation. In the second one, we use and compare dierent type of asymptotically valid condence intervals.

Simulation 1: Asymptotic properties

We consider the set of sample sizes n ∈ {100, 300, 500, 1000} to analyze the consistency of the mle of parameters. We x α = 0.

where θ ij is the mle of parameter θ i , computed in the j-th sample.

Tables 13 show that both RB and RRSE decrease when the sample size increases, for all the parameter estimates.


Patrice Bertail

Prof Aldo M Garay
email: agaray@de.ufpe.br

Francyelle L Medina

Isaac Jales

Moraes Rego

A maximum likelihood and regenerative bootstrap approach for estimation and forecasting of INAR(p) processes with zero-inated innovations

Keywords: Count time series, ZI-INAR(p) process, Bootstrap, EM algorithm, Zero-inated distribution

In this work, we study a class of p-order integer-valued non-negative autoregressive (INAR(p)) processes, with innovations following zero-inated (ZI) distributions called ZI-INAR(p) processes. Based on the EM algorithm, we present an estimation procedure of the parameters of the model. We also develop a regenerative bootstrap method to construct condence intervals for the parameters, and more generally Fréchet dierentiable functionals, as well as forecasting distributions for future values. We study the performance of the proposed method and demonstrate the utility of our approaches through a simulation study and an application to a real dataset.

Introduction

To deal with discrete time series models, a common approach is to consider the nonnegative integer-valued autoregressive (INAR) process proposed by [START_REF] Mckenzie | Some simple models for discrete variate time series[END_REF] and [START_REF] Al-Osh | First-order integer-valued autoregressive (INAR(1)) process[END_REF]. This process is based on the binomial thinning operator and a Poisson distribution for the innovations. But in some cases, in which the dataset presents excess of zeros, the Poisson distribution may be inappropriate. In this context, [START_REF] Jazi | First-order integer valued AR processes with zero inated Poisson innovations[END_REF] introduced the ZINAR(1) process, assuming that the innovations follow a zero-inated Poisson distribution. [START_REF] Barreto-Souza | Zero-Modied Geometric INAR(1) process for modelling count time series with deation or ination of zeros[END_REF] and [START_REF] Risti¢ | Zero-inated NGINAR(1) process[END_REF] developed the INAR(1) process with zero-inated geometric (ZMG) marginal distribution. In order to analyze a p-order counting process with excess of zeros, [START_REF] Maiti | Modelling and coherent forecasting of zeroinated count time series[END_REF] proposed the ZIPPAR(p) process, which assumes a zero-inated Poisson marginal distribution and uses a Pegram operator; [START_REF] Garay | Bayesian analysis of the p-order integer-valued ar process with zero-inated poisson innovations[END_REF] extended the ZINAR(1) process and proposed the ZINAR(p) process based on the binomial thinning operator, with zero-inated Poisson innovations, and developed a Bayesian approach for parameter estimation and prediction. These models are based on dierent operators and assume the zero-inated Poisson distribution. In this paper, to model discrete time series with excess of zeros and order dependence larger than one, we propose the INAR(p) process based on the binomial thinning operator, considering a class of zero-inated (ZI) innovations. We develop an ecient EM algorithm for estimating the parameters by maximum likelihood. However, it is known that the estimators may be asymptotically Gaussian with a complicated variance. Thus, in order to construct condence intervals for the parameters and forecasting distribution for future values, we propose specic bootstrap resampling techniques.

It is important to recall that in the context of the INAR processes, dierent bootstrap methods have been developed, among which we can mention [START_REF] Cardinal | On the application of integer-valued time series models for the analysis of disease incidence[END_REF], [START_REF] Jung | Coherent forecasting in integer time series models[END_REF], [START_REF] Kim | A non-stationary integer-valued autoregressive model[END_REF], [START_REF] Kim | Coherent forecasting in binomial AR(p) model[END_REF], [START_REF] Jentsch | Bootstrapping INAR models[END_REF] and [START_REF] Bisaglia | Model-based INAR bootstrap for forecasting INAR(p) models[END_REF]. In this manuscript we use the regenerative bootstrap, studied by [START_REF] Bertail | Regenerative block bootstrap for markov chains[END_REF], extending the approaches of [START_REF] Athreya | Bootstrapping Markov chains: countable case[END_REF] and [START_REF] Datta | On the First-Order Edgeworth Expansion for a Markov Chain[END_REF], which exploit the regeneration properties of Markov chains when a (recurrent) state is innitely often visited by the chain. The main idea consists of resampling a deterministic number of data blocks corresponding to regeneration cycles, involving a successive visit to a well-chosen observation. Because of some inadequate standardization, the regeneration-based bootstrap method proposed in [START_REF] Datta | On the First-Order Edgeworth Expansion for a Markov Chain[END_REF] is not always second-order correct (its rate is O P (n -1/2 ) only). However, [START_REF] Bertail | Regenerative block bootstrap for markov chains[END_REF] proposed a sequential modication which is second order correct, even if the chain is not stationary. They extended this approach to general Harris Markov chains (using the Nummelin splitting approach), but their simulation studies showed that the method works the best when the chain is atomic or takes a nite number of values. Since the INAR process, and more generally integer-valued processes, can be approximated by Markov chains with atoms (each visited point or sequence of visited points can be seen as an atom), the regenerative approach appears natural in this context. We develop the asymptotic validity of the method in the context of p-order Markov chains for general dierentiable functionals. This allows us to get both the asymptotic normality of the maximum likelihood estimation (m.l.e.) and the asymptotic validity of the regenerative bootstrap of the m.l.e. This bootstrap method is also robust to the model specication (since we only use the Markov assumption and not the specic form of the model).

This manuscript is organized as follows. Section 2 introduces some denitions, notations and characteristic features of the p-order non-negative integer-valued autoregressive process (INAR(p)) and the zero-inated (ZI) models. Section 3 outlines the proposed ZI-INAR(p) model and discusses some statistical properties, including the description of the likelihood function. Section 4 presents the implementation of an EM-type algorithm to estimate the model's parameters. We then propose a regenerative bootstrap method to construct condence intervals for the parameters and for estimating the forecasting distributions for future values. In order to illustrate the usefulness of the proposed method, several articial and real datasets are analyzed in Sections 5 and 6, respectively. Finally, Section 7 concludes with short remarks and some possible avenues for future research.

Preliminaries

The INAR(p) process, presented by [START_REF] Jin-Guan | The integer-valued autoregressive (INAR(p)) model[END_REF], is a model to deal with count time series. Suppose that Y is a non-negative integer-valued random variable (r.v.) and α ∈ [0, 1], then the thinning operator • is dened as:

α • Y = Y i=1 Z i ,
where {Z i } i∈Z is a sequence of independent and identically distributed (iid) random variables, independent of Y , which follows a Bernoulli distribution with parameter α.

Denition 1. The INAR(p) process is dened as follows:

Y t = p i=1 α i • Y t-i + V t , t ∈ Z, (1)
where {V t } t∈Z is an iid sequence of non-negative integer r.v.'s, dened as innovations, with mean µ V and variance σ 2 V < ∞, and α i ∈ [0, 1] for all i = 1, . . . , p. As commented by [START_REF] Jin-Guan | The integer-valued autoregressive (INAR(p)) model[END_REF], given {Y t-1 = y t-1 , . . . , Y t-p = y t-p }, the r.v.'s α i •Y t-i , are conditionally independent. On the other hand, the sequence of innovations

{V t } t∈Z is independent of Y t-1 , Y t-2 , . . .. If y t-i > 0, then α i • Y t-i |Y t-i = y t-i ∼ Bin(y t-i , α)
, where Bin(a,b) denotes the binomial distribution, with parameters a and b.

If y t-i = 0, then α i • Y t-i |Y t-i = y t-i
follows a degenerate distribution at zero.

In the following, we present briey two important results about the stationarity condition of the INAR(p) process, which will be useful for our model.

-Jin-Guan and Yuan (1991) showed that the stationarity conditions of INAR(p) and AR(p) are the same, as stated in the following theorem: Theorem 1. Let {V t } t∈Z be i.i.d. non-negative integer-valued random variables with

E [V t ] = µ V , V ar [V t ] = σ 2
V < ∞ and suppose that α i ∈ [0, 1] , for i = 1, . . . , p. If the roots of the polynomial

x p -α 1 x p-1 -• • • -α p-1 x -α p = 0
are inside the unit circle, then there exists a unique stationary non-negative integervalued random series {Y t } satisfying

Y t = α 1 • Y t-1 + • • • + α p • Y t-p + V t ; Cov (Y s , V t ) = 0 for s < t.
Proof. See Theorem 2.1 in Jin-Guan and Yuan (1991).

- [START_REF] Doukhan | On weak dependence conditions: The case of discrete valued processes[END_REF] proved that for discrete valued processes, the weak dependence mixing conditions are satised under natural assumptions (whereas most of the times strong mixing assumptions are not satised). In particular, the authors further showed that p i=1 α i < 1 is a sucient condition for stationarity and ergodicity of the INAR(p) process.

Several zero inated models have been proposed to deal with datasets with excess of zeros. In the following, we present the zero-inated distributions and their hierarchical formulation, as well as some particular cases and properties. Denition 2. A non-negative discrete r.v. V , with parameters π and λ, follows a zeroinated distribution, denoted by V ∼ ZI(π, λ), if it can be represented by the stochastic form:

V = BU, with B⊥U, (2)
where the r.v. B follows a Bernoulli distribution with mean 1 -π, for π ∈ [0, 1). U is a non-negative discrete r.v. with nite variance and parameter vector λ. B⊥U indicates that the r.v.'s B and U are independent.

From Denition 2, we have that the probability function (pf ) of V is given by:

P (V = v) = π + (1 -π) h U (0|λ) , if v = 0; (1 -π) h U (v|λ) , if v ≥ 1,
where h U (•|λ) denotes the pf of U . We have that:

E (V ) = (1 -π) E (U ) and V ar (V ) = (1 -π) V ar (U ) + πE 2 (U ) .
Hereafter, we describe three particular cases of the ZI models.

The zero-inated Poisson (ZIP) model:

When U follows a Poisson distribution with mean λ, we say that V follows a ZIP distribution, denoted by V ∼ ZIP(π, λ). In this case, the pf of the r.v. V is given by:

P (V = v) = π + (1 -π) e -λ , if v = 0; (1 -π) e -λ λ v v! , if v ≥ 1,
and

E (V ) = (1 -π) λ and V ar (V ) = (1 -π) λ + πλ 2 .
The zero-inated negative binomial (ZINB) model:

The ZINB model is obtained when U follows a negative binomial distribution, with mean µ ≥ 0 and dispersion parameter φ > 0. In this case, λ = (µ, φ) and we have that:

P (V = v) =      π + (1 -π) φ µ+φ φ , if v = 0; (1 -π) Γ(φ+v) Γ(v+1)Γ(φ) µ µ+φ v φ µ+φ φ , if v ≥ 1,
where Γ(•) denotes the gamma function. We use the notation V ∼ ZINB(π, µ, φ). Note that if π = 0, then V follows a negative binomial distribution, denoted by V ∼ NB(µ, φ)

E (V ) = (1 -π)µ and V ar (V ) = (1 -π) µ 1 + µ φ + πµ . (3)
For more details, see [START_REF] Hall | Zero-inated Poisson and Binomial regression with random eects: a case study[END_REF], [START_REF] Ridout | A score test for testing a zero-inated Poisson regression model against zero-inated negative binomial alternatives[END_REF] and [START_REF] Garay | On estimation and inuence diagnostics for zero-inated negative binomial regression models[END_REF], among others.

The zero-inated Poisson inverse Gaussian (ZIPIG) model:

In this case, U follows a Poisson inverse Gaussian (PIG) distribution with mean µ ≥ 0 and dispersion parameter φ > 0. Thus, λ = (µ, φ) and the pf of the r.v. V is given by:

P (V = v) =        π + (1 -π)e φ- √ φ(φ+2λ) , if v = 0; (1 -π) 2 π [φ(φ + 2λ)] -(v-1/2) 2 e φ (λφ) v v! K v-1/2 φ(φ + 2λ) , if v ≥ 1, where K λ (t) = 1 2 ∞ 0 u λ-1 e -t 2 (u+ 1 u )
du is the modied Bessel function of the third kind [START_REF] Abramowitz | Handbook of mathematical functions: With formulas, graphs, and mathematical tables applied mathematics series[END_REF]. We denote V ∼ ZIPIG(π, µ, φ). The expectation and variance of the ZIPIG model coincide with the expectation and the variance of the ZINB model, given in Eq. (3). When π = 0, we have that r.v. V ∼ PIG(µ, φ). We recall that the PIG distribution belongs to the class of the mixed Poisson distributions, that is, if X ∼ PIG(µ, φ), then X has the following stochastic structure: X|Z = z ∼ Poisson(µz), where Z ∼ IG(1, φ). The notation IG(1, φ) represents the inverse Gaussian distribution with mean 1 and dispersion parameter φ > 0. More details, about the mixed Poisson distributions are presented and discussed by [START_REF] Karlis | A general em approach for maximum likelihood estimation in mixed poisson regression models[END_REF] and [START_REF] Barreto-Souza | General mixed poisson regression models with varying dispersion[END_REF].

In the analysis of discrete time-series data with excess of zeros, the INAR(1) process with marginals or innovations following the zero-inated Poisson distribution, as proposed by [START_REF] Maiti | Modelling and coherent forecasting of zeroinated count time series[END_REF] and [START_REF] Garay | Bayesian analysis of the p-order integer-valued ar process with zero-inated poisson innovations[END_REF], may be inadequate. Thus, to provide a more exible process, with dierent zero-inated models, we propose and study in this manuscript the INAR(p) process with zero-inated innovations, namely ZI-INAR(p) processes.

The ZI-INAR(p) processes

The ZI-INAR(p) processes are integer-valued p-order autoregressive time series with ZI innovations, as dened by Eq. ( 1), where {V t } t∈Z follows a zero-inated distribution. Since the ZI-INAR(p) processes have the same stochastic structure as the INAR(p) process and considering the results obtained by Jin-Guan and Yuan (1991) and [START_REF] Doukhan | On weak dependence conditions: The case of discrete valued processes[END_REF] presented above, we assume p i=1 α i < 1 to ensure the stationarity of the ZI-INAR(p) process. We now present some mathematical properties of these processes.

Mathematical properties

Let {Y t } t∈Z be a stationary ZI-INAR(p) process; then the mean and the covariance function γ(•) of Y t are given by:

E (Y t ) = (1 -π)E(U t ) 1 - p i=1 α i and γ (h) = p i=1 α i γ (h -i) + (1 -π) V ar(U t ) + πE 2 (U t ) I {h} (0),
where I A (•) denotes the indicator function, i.e., I A (y) = 1, if y ∈ A and I A (y) = 0 otherwise.

These expressions can be obtained by extending the results presented by Jin-Guan and Yuan (1991), who showed that the INAR(p) process has the same correlation structure as the standard AR(p) process.

The Likelihood function

Let y = (y 1 , . . . , y n ) be realizations of the ZI-INAR(p) process. The likelihood function of the vector of parameters θ = (α, π, λ) , with α = (α 1 , . . . , α p ) , given y, is dened by:

L (θ|Y) = P (Y 1 = y 1 , . . . , Y p = y p ) n t=p+1 P (Y t = y t |Y t-1 = y t-1 , . . . , Y t-p = y t-p ) , (4)
where

P (Y t = y t |Y t-1 = y t-1 , . . . , Y t-p = y t-p ) = min{y t-1 ,yt} k 1 =0 y t-1 k 1 α k 1 1 (1 -α 1 ) y t-1 -k 1 × min{y t-2 ,yt-k 1 } k 2 =0 y t-2 k 2 α k 2 2 (1 -α 2 ) y t-2 -k 2 × . . . × min{y t-p ,yt-p-1 i=1 k i} kp=0 y t-p k p α kp p (1 -α p ) y t-p -kp × πI {yt-p i=1 k i} (0) + (1 -π)h Ut y t - p i=1 k i |λ .
Eq. ( 4) is expressed in terms of the joint probability P (Y 1 = y 1 , . . . , Y p = y p ), which is not available. Thus, the exact likelihood function is intractable. An approach to estimate the parameters is to use the conditional log-likelihood function, given by:

(θ|y) ∝ n t=p+1 log [P (Y t = y t |Y t-1 = y t-1 , . . . , Y t-p = y t-p )] , (5)
which is obtained by ignoring the initial values P (Y 1 = y 1 , . . . , Y p = y p ). See [START_REF] Hamilton | Time Series Analysis[END_REF], [START_REF] Bu | Maximum likelihood estimation of higher-order integervalued autoregressive processes[END_REF], Prado and West (2010, Sec 2.3.5), [START_REF] Jazi | First-order integer valued AR processes with zero inated Poisson innovations[END_REF], [START_REF] Weiss | An Introduction to Discrete-Valued Time Series[END_REF] and [START_REF] Garay | Bayesian analysis of the p-order integer-valued ar process with zero-inated poisson innovations[END_REF] for more details and references. For ZI-INAR(p) processes, the direct maximization of the expression ( 5) is not easy due to its form. An alternative is to consider numerical optimization using the EM algorithm. Its properties ensure the monotone convergence to a stationary point of the log-likelihood function, in contrast to direct maximization. Recall that one may converge to a local solution, according to the choice of the initial values.

Maximum likelihood estimation and bootstrap resampling methods

The EM algorithm

In order to develop the EM algorithm for maximum likelihood estimation (m.l.e.) of the parameters of the ZI-INAR(p) process, we consider the augmented data representation of the process, by combining the observed data and some latent variables described below.

Thus, as suggested by [START_REF] Hall | Zero-inated Poisson and Binomial regression with random eects: a case study[END_REF], [START_REF] Neal | MCMC for integer-valued ARMA processes[END_REF] and [START_REF] Garay | Bayesian analysis of the p-order integer-valued ar process with zero-inated poisson innovations[END_REF], considering t = p + 1, . . . , n and i = 1, . . . , p, we dene the following latent variables:

-S t,i = α i • Y t-i ; where S t,i |Y t-i = y t-i follows a binomial distribution, with parameters y t-i and α i , when y t-i > 0 and if y t-i = 0, S t,i follows a degenerate distribution at zero.

-From Denition 2, the distribution of V t is a mixture with components h U (•|λ) and the degenerate distribution at zero, with weights 1 -π and π, respectively. Thus, a latent binary random variable W t exists such that:

W t ∼ Bern(π),
where Bern(π) denotes the Bernoulli distribution with parameter π, and

V t |W t = 0 ∼ h U (•|λ) , V t |W t = 1 follows a degenerate distribution at zero. (6)
-From Eq. ( 1), we have that V t = Y t -p i=1 S t,i and as commented in [START_REF] Garay | Bayesian analysis of the p-order integer-valued ar process with zero-inated poisson innovations[END_REF], considering y t(-p) = {y t-1 , y t-2 , . . . , y t-p } for t = p + 1, . . . , n, we have that S t,1 , S t,2 , . . . , S t,p are independent when y (t-p) is given; the same is true for (Z t , W t ) and (S t,1 , S t,2 , . . . , S t,p ).

Let Y = (Y 1 , . . . , Y n ) , W = (W p+1 , . . . , W n ) , S = (S p+1 , . . . , S n ) , with S t = (S t,1 , . . . , S t,p ) and Y c = Y c p+1 , . . . , Y cn . Then, the joint probability function of Y ct = (Y t , W t , S t ) is given by: P Y ct = y ct |Y t(-p) = y t(-p) = P Y t = y t , W t = w t , S t = s t |Y t(-p) = y t(-p) = P S t = s t |Y t(-p) = y t(-p) × P Y t = y t , W t = w t |S t = s t , Y t(-p) = y t(-p) = p i=1 y t-i s t,i α s t,i i (1 -α i ) y t-i -s t,i π wt (1 -π)h Ut y t - p i=1 s t |λ 1-wt
.

Consequently, the complete log-likelihood function can be expressed as:

c (θ|y c ) = log   P (Y 1 = y 1 , . . . , Y p = y p ) n t=p+1 P Y ct = y ct |Y t(-p) = y t(-p)   ∝ p i=1      n t=p+1 s t,i   log (α i )    + p i=1      n t=p+1 (y t-i -s t,i )   log (1 -α i )    + n t=p+1 w t log (π) + n t=p+1 (1 -w t ) log (1 -π) + n t=p+1 (1 -w t ) log h Ut (y t - p i=1 s t,i |λ) .
The EM-type algorithm for parameter estimation contains two steps. For the E-step, let

θ (k)
be the k-th step estimate of θ. The Q-function is calculated as:

Q θ| θ (k) = E c (θ|y c )|y, θ (k)
.

In the M-step, the Q-function Q θ| θ

(k) is maximized with respect to θ = (α, π, λ) , obtaining θ (k+1)
.

For the parameter estimation of the ZI-INAR(p) processes, the expression of the Qfunction is completely determined by the knowledge of the following expectations:

s (k) t,i = E S t,i |y, θ (k)
, w

(k) t = E W t |y, θ (k)
, as well as

Q * t λ| θ (k) = E (1 -W t ) log h Ut y t - p i=1 S t,i |y, θ (k)
.

Finally, the Q-function can be written as:

Q θ| θ (k) ∝ p i=1      n t=p+1 s (k) t,i   log (α i )    + p i=1      n t=p+1 y t-i -s (k) t,i   log (1 -α i )    + n t=p+1 w (k) t log (π) + n t=p+1 1 -w (k) t log (1 -π) + n t=p+1 Q * t λ| θ (k)
.

The maximization of Q θ| θ

(k)
over θ must be obtained under restrictions on the parameter α in order to ensure the standard condition for stationarity and ergodicity of the ZI-INAR(p) process ( see [START_REF] Doukhan | On weak dependence conditions: The case of discrete valued processes[END_REF]). Recall that this process will be stationary if and only if α i ∈ [0, 1), for i = 1, . . . , p, and p i=1 α i < 1. In order to satisfy these restrictions, we establish the following proposition, which proposes the transformation in ∆ : (α 1 , . . . , α p ) → (β 1 , . . . , β p ) :

Proposition 1. Let ∆ : (α 1 , . . . , α p ) → (β 1 , . . . , β p ) be a transformation in α = (α 1 , . . . , α p ) , such that

β i = ∆ i (α) = α i /(1 -p j =i α j ), i = 1, . . . , p. If α i ∈ [0, 1) and p i=1 α i < 1.
We then have that the ∆ transformation admits an inverse given by

α i = ∆ -1 i (β) = (1 - p i=1 β i 1-β i 1 + p i=1 β i 1-β i ) β i 1 -β i , i = 1, ..., p,
where β = (β 1 , . . . , β p ) .

Proof . Considering p i=1 α i = c ∈ [0, 1) we have that β i = α i 1-c+α i and it follows by direct inversion that

α i = (1 -c) β i 1 -β i .
Now by summing this expression of the α i values, we get

c = (1 -c) p i=1 β i 1 -β i yielding to the expression of c in terms of β i given by c = p i=1 β i 1-β i 1 + p i=1 β i 1-β i leading to the result.
β (k+1) = arg max β∈[0,1) p Q ∆ -1 (β), π, λ| θ (k) α (k+1) i = ∆ -1 i β (k+1) π (k+1) = n t=p+1 w t n -p λ (k+1) = arg max λ    n t=p+1 Q * t λ| θ (k)    , (7)
where

∆ -1 (β) = (∆ -1 1 (β), . . . , ∆ -1 p (β)) .
It is important to note that to obtain the expressions

Q * t λ| θ (k) and λ (k+1) , it is necessary to compute the expectation b t s t,i = E B t S t,i |y, θ (k)
, where B t = 1 -W t .

Thus, we present the expressions Q * t λ| θ

(k)
and λ

(k+1)

considering the particular case of ZI models:

-If V t ∼ ZIP(π,λ), then:

Q * t λ| θ (k) ∝ -λ(1 -w (k) t ) + log(λ)(1 -w (k) t )y t -log(λ) p i=1 b t s t,i (k)
and by maximizing the function

Q * t λ| θ (k)
, over λ, we have that:

λ (k+1) = n t=p+1 (1 -w (k) t )y t - n t=p+1 p i=1 b t s t,i (k) n t=p+1 (1 -w (k) t )
.

-If V t ∼ ZINB(π,λ), then λ = (µ, φ) and

Q * t µ, φ| θ (k) ∝ g t (φ) (k) + [log(µ) -log(µ + φ)] (1 -w (k) t )y t - p i=1 b t s t,i (k) + -log (Γ(φ)) + φ [log(φ) -log(µ + φ)] 1 -w (k) t , (8)
where g t (φ)

(k) = E B t log Γ (y t -p i=1 S t,i + φ) |y, θ (k)
.

In order to circumvent the computation of g t (φ)

(k)

, we adopt the ECM algorithm, proposed by [START_REF] Liu | The ECME algorithm: a simple extension of EM and ECM with faster monotone convergence[END_REF], which is an extension of the EM algorithm that is obtained by replacing the M-step by a sequence of conditional maximization steps. Thus, we have that:

µ (k+1) = n t=p+1 (1 -w (k) t )yt - n t=p+1 p i=1 bts t,i (k) n t=p+1 (1 -w (k) t )
and k+1) , φ|y ,

φ (k+1) = arg max φ α (k) , π (k) , µ (
where (•|y) is given by Equation ( 5).

-If V t ∼ ZIPIG(π,µ,φ), then U t ∼ PIG(µ,φ) and λ = (µ, φ). Consider now the stochastic representation of the mixed Poisson model for U t , that is:

U t |Z t = z t ∼ Poisson(µz t ) (9)
and Z t ∼ IG(1, φ), for t = p + 1, . . . , n.

Thus, the complete-data for the ZIPIG cases is now dened by Y c = (Y, W, S, Z) and the joint probability function of Y ct = (Y t , W t , S t , Z t ) is given by:

f (y ct |y t(-p) ) = P S t = s t |Y t(-p) = y t(-p) × f (y t , w t , z t |s t , y t(-p) ) = p i=1 y t-i s t,i α s t,i i (1 -α i ) y t-i -s t,i π wt (1 -π)h Ut,Zt y t - p i=1 s t , z t |λ 1-wt
, where h Ut,Zt (•|λ) is the joint density function of (U t , Z t ), that is:

h Ut,Zt (u t , z t |λ) = P(U t = u t |Z t = z t ) × f (z t ). Consequently, the Q-function Q * t λ| θ (k)
is given by:

Q * t λ| θ (k) = E (1 -W t ) log P U t = y t - p i=1 S t,i |Z t = z t × f (z t ) |y, θ (k) ∝ -µ b t z t (k) -log(µ) n i=p b t s t,i (k) + log(µ)(1 -w (k) t )y t + log(φ) 2 (1 -w (k) t ) + φ(1 -w (k) t ) - φ 2 b t z t (k) - φ 2 b t /z t (k) , ( 10
)
where b t z t (k) = E B t Z t |y, λ (k) and b t /z t (k) = E B t Z -1 t |y, λ (k)
.

Thus, from Equations ( 7) and (10), we have that:

µ (k+1) = n t=p+1 (1 -w (k) t )y t - n t=p+1 p i=1 b t s t,i (k) n t=p+1
(1 -w

(k) t )
and

φ (k+1) = n t=p+1 (1 -w (k) t ) n t=p+1 b t z t (k) + n t=p+1 b t /z t (k) -2 n t=p+1 (1 -w (k) t )
.

The procedure is iterated until some convergence criterion is satised. In this manuscript, |y , where (θ|y) is given in

(5), the Aitken acceleration-based stopping criterion uses the fact that the limit of the sequence (1) , (2) , . . ., denoted by ∞ , can be aproximated by

(k+1) ∞ = (k) + ( (k+1) - (k) )/(1 -c (k) ), where c (k) = ( (k+1) -(k) )/( (k) -(k-1)
). As suggested by Zeller et al.

(2019), we stop the algorithm when

(k+1) ∞ -(k+1) < ε = 10 -5 .
All the expressions and details of our EM algorithm are presented in Appendix A of the supplementary material.

Regenerative bootstrap method

In what follows, we present and discuss some denitions and asymptotic properties of the regenerative bootstrap method that justify the use of this technique for the ZI-INAR(p) processes.

The regenerative bootstrap

We briey recall the method developed by [START_REF] Bertail | Regenerative block bootstrap for markov chains[END_REF], assuming that the process is a Markov chain of order 1. It is assumed that the process is an aperiodic irreducible recurrent Markov chain with an atom A. For instance, for the ZI-INAR(1) process, the atom may be chosen to be A = {0}. However, if the chain is of length p = 1 (for instance for ZI-INAR(p) processes), one can vectorize the chain and choose A to be a sequence of p integers appearing a great number of times in the chain. In fact, the method is asymptotically independent of the choice of the atom, but in practice it is important to choose an atom which is visited a great number of times. For that, in an observed times series, one may count how many times each state is visited, then how many times a sequence of states, for instance {(0, 1)}, {(0, 2)}, ...{(0, 1, 2)}, {(0, 1, 3)}, is visited to determine empirically the order of the chain and the adequate atom. To simplify the exposition, we assume here that A = {a} is reduced to a point. Denote by τ 0 = τ 0 (1) = inf {n ≥ 1, Y n = a} the hitting time on a. Then the hitting time times are dened sequentially by τ 0 (j) = inf {n > τ 0 (j -1), Y = a} for j ≥ 2. The sequence {τ 0 (j)} j 1 denes successive times called regeneration times, at which the chain forgets its past. It follows from the strong Markov property that, for any initial law ν of Y 0 , the sample paths can be divided into i.i.d. blocks of random length corresponding to consecutive visits to A, also called cycles:

B 1 = (Y τ 0 (1)+1 , . . . , Y τ 0 (2) ), . . . , B j = (Y τ 0 (j)+1 , . . . , Y τ 0 (j+1) ), . . . , . . . , B ln-1 = (Y τ 0 (ln-1)+1 , . . . , Y τ 0 (ln) ), . . . If the chain does not start at Y 0 = a, then the rst block B 0 = (Y 1 , ..., Y τ 0 (1) ) is independent of the next ones but has a dierent distribution, depending on the law of the initial distribution. Similarly, the last block B ln = Y τ 0 (ln)+1 , ..., Y n , in a stretch of length n, is independent of the proceeding ones but may be incomplete, in that it does not nish at the atom {a}. This explains why these blocks should be removed in the sequel.

Consider now a general statistic T n of interest. For sake of generality, we also assume there eventually exists an adequate standardization S n . (typically an estimator of the standard deviation of the statistics of interest). For the reason invoked before regarding B 0 and B ln , we assume that the statistics T n and the estimate of the variance are constructed using the data blocks B 1 , . . . , B ln-1 only. We are thus interested in estimating accurately its sampling distribution under P ν (.), say H (n)

Pν (y) = P ν (S -1 n (T n -θ) y).
The regenerative procedure is performed in ve steps as follows.

Step 1. Count the number of visits l n to the atom A up to time n. Divide the observed sample path

Y (n) = (Y 1 , . . . , Y n ) into l n + 1 blocks, B 0 , B 1 , . . . , B ln-1 , B (n)
ln . Drop the rst and last (non-regenerative) blocks.

Step 2. Draw sequentially bootstrap data blocks B * 1,n , . . . , B * k,n independently from the em-

pirical distribution F n = (l n -1) -1 ln-1 j=1 δ B j of the blocks {B j } 1 j ln-1 conditional on Y (n) , until the length l * (k) = k j=1 l(B * j,n ) of the bootstrap data series is larger than n. Let l * n = inf{k 1, l * (k) > n}.
The fact that the number of blocks in the reconstructed bootstrap path is random is crucial for the good properties of the method.

Step n -1,n ). The RBB distribution is then given by:

H * n (x) = P * (S * -1 n (T * n -T n ) y | Y (n) )
where P * (. | Y (n) ) denotes the conditional probability given Y (n) .

Step 5. Repeat independently the procedure above in Steps 2 to Step 4 B times to compute successive values S * -1 b,n (T * b,n -T n ) , b = 1, ..., B and a Monte-Carlo approximation of H * n (y) , say:

H * B,n (y) = 1 B B b=1 I {S * -1 b,n (T * b,n -Tn) y}.
B can be chosen so that it does not aect the good properties of the bootstrap distribution: typically of order B = O( √ n), see [START_REF] Bertail | Regenerative block bootstrap for markov chains[END_REF]. When one is interested in computing two-sided condence intervals of coverage γ, it has been recognized since the work of [START_REF] Hall | On the Number of Bootstrap Simulations Required to Construct a Condence Interval[END_REF] that choosing B such that (1-γ)/2(B +1) is an integer will be more ecient. This is due to the fact that the quantile of order (1-γ)/2 of the bootstrap distribution is exactly the largest value of (B + 1)(1 -γ)/2 among the order values of the quantities S * -1 b,n (T * b,n -T n ). Most of the time B = 999 (for time series smaller than 1000) is a good choice. We choose B = 999 in our applications. This procedure is asymptotically correct and can be applied, for instance, to maximum likelihood estimation of the parameters of interest for the ZI-INAR(p) processes.

Asymptotic properties

In what follows, P ν (respectively P a for a a given atom) will denote the probability measure of the underlying probability space such that Y 0 ∼ ν (resp. Y 0 = a), E µ (.) is the P ν -expectation (resp. E a (.) being P a -expectation). Recall that if the chain is Markov of order p, we can always vectorize it, so as to get a chain of order 1. In that case the atom is a point a ∈ R p . Assume that the chain is irreducible, aperiodic, positive recurrent, such that the return time to the atom starting from a has nite expectation

E a τ 0 < ∞.
Then it is known that the chain admits a stationary measure say µ. Let f : E → be a µ-integrable function and consider the empirical mean

µ n (f ) = n -1 n i=1 f (Y i ) of the unknown mean µ(f ) = E µ (f (Y 1
)) constructed on a data segment of size n, say X (n) . [START_REF] Bertail | Edgeworth expansions of suitably normalized sample mean statistics for atomic Markov chains[END_REF] (see Proposition 3.1) showed that if ν = µ, the chain is not stationary and the rst data block B 0 induces a signicant bias depending on the the unknown measure ν. This is well understood in the Bayesian literature; whereby when starting a Markov chain from an initial distribution ν = µ, one has to wait for a long time before reaching stationarity. This is known as the burn-in period. The last block also induces some bias. This suggests dropping the rst and last blocks and instead using:

µ n (f ) = (τ A (l n ) -τ A ) -1 ln-1 i=1 S j (f )
(with the convention that empty sums are equal to 0), where we put

f (B j ) = τ A (j+1) i=1+τ A (j) {f (X i ) -µ(f )}
for the sum over a block (of random length l(B j ) = τ 0 (j + 1) -τ 0 (j)).

We will only focus on asymptotic properties and not second-order theory. They do not hold (directly) since the processes under consideration here are not continuous and essentially take discrete values (so that the usual Cramer condition for Edgeworth expansions may fail).

In this framework, we may use indierently µ n (f ) or µ n (f ). Their asymptotic behavior and the asymptotic validity of the bootstrap hold for both estimators. It is then easy to prove the following result (see Meyn and Tweedie, 1996) for the mean.

Theorem 2. Central limit theorem Assume that E ν τ 0 < ∞ and E a τ 2 0 < ∞, and further assume that we have the block moment conditions

E a   τ 0 i=1 |f |(Y i ) 2   < ∞ and dene σ 2 (f ) = E a S 1 (f ) 2 E a τ 0
which is is assumed to be non degenerate and nite

0 < σ 2 (f ) < ∞
Then for any starting distribution ν, we have

√ n (µ n (f ) -µ (f )) d - → N 0, σ 2 (f ) as n → ∞.
In the case of INAR(p), the block moment conditions may be replaced by drift conditions which are satised under more easily veriable conditions involving the parameters of the time series: this is done in [START_REF] Drost | Note on integer-valued bilinear time series model[END_REF], Theorem 2.2. These conditions trivially reduce to p i=1 α i < 1 and existence of the moment of order 2 of the noise.

In light of the expression of the asymptotic variance and recalling that n ln-1 → E a τ 0 a.s., it is easy to see that the following plug-in estimator of the asymptotic variance,

σ 2 n (f ) = n -1 ln-1 j=1 (f (B j ) -µ n (f )l(B j )) 2 , ( 11
)
is a convergent estimator of the asymptotic variance. The bootstrap versions of µ n (f ) and σ 2 n (f ) are dened by:

µ * n (f ) = n * -1 a l * n -1 j=1 f (B * j ) and σ * 2 n (f ) = n * -1 a l * n -1 j=1 f (B * j ) -µ * n (f )l(B * j ) 2 , ( 12
)
with n * a = l * n -1 j=1 l(B * j ) being the total size of the full blocks in the bootstrap reconstruction. Then, the results of [START_REF] Bertail | Regenerative block bootstrap for markov chains[END_REF] apply. In the following, d -→ * means convergence in distribution conditionally to the data, and it holds either in probability along the data or a.s. along the data, meaning that the cdf of the bootstrap distribution converges either in probability or almost surely to the true distribution.

Theorem 3. Assume that the conditions of Theorem 2 are satised and that in addition,

E ν τ 2 0 < ∞, E a τ 6 0 < ∞ and E a   τ 0 i=1 |f |(X i ) 6   < ∞.
Then, for any starting distribution ν, we have:

√ n (µ * n (f ) -µ n (f )) d - → * N 0, σ 2 (f ) almost surely as n → ∞.
This result ensures the validity of the bootstrap method for the mean, for the ZIN AR(p) processes. If one is interested in µ n additional moment conditions are needed for the moment of order 2 under ν of the rst block to ensure that it does not overly perturb the statistics of interest. For instance, we may have problems if the initial distribution has a Pareto type tail with no moments. Then intuitively in that case, the return time to the atom may be so long that we cannot observe any complete blocks. Now, the parameter θ in the m.l.e procedure can be seen as a regular M-estimator in terms of the distribution of a block (see [START_REF] Bertail | Bootstrapping robust statistics for markovian data applications to regenerative r-statistics and l-statistics[END_REF]). Each block can be seen as a random variable in the torus T = ∪ ∞ k=1 R k . Indeed, the M-estimator dening the m.l.e., θ n , can be seen as the solution of the block-estimating equation

ln-1 j=1 ∂ ∂θ l(B j , θ) = 0 where l(B i , θ) = log P a (Y τ 0 (i)+1 = y τ 0 (i)+1 , ...., Y τ 0 (i+1) = a )
is the log likelikood corresponding to the i-th block.

We recall a few denitions in this framework. Let P T denote the set of all probability measures of the torus T. A Markov functional is a function T : P T -→ (V, ||.||) with its value in some separable Banach space V endowed with a norm ||.||. If, for all L in P T ,

T (1) (b, L) = lim t→0 T ((1 -t)L + tδ b ) -T (L) t
exists for any b ∈ T, it is called the inuence function of T based on blocks.

We now recall the notion of Frechet dierentiability of functionals based on blocks.

Denition 3. A functional T : P T → ( V,||.||) is said to be Fréchet dierentiable at L 0 ∈ P T , for a norm ||.||, if there exists a continuous linear operator DT L 0 (from the set of signed measures of the form L -L 0 ) and a function (1) (•, L 0 ) : R → (, • ), which is continuous at 0 with (1) (0, L 0 ) = 0, such that

∀ L ∈ P, T (L) -T (L 0 ) = DT L 0 (L -L 0 ) + R (1) (L, L 0 ), where R (1) (L, L 0 ) = ||L -L 0 || (1) (||L -L 0 ||, L 0 ).
In this framework, its is also assumed that T has an inuence function T (1) (•, L 0 ), and that the following integral representation holds for DT L 0 :

∀ L 0 ∈ P, DT L 0 (L -L 0 ) = E T (1) (x, L 0 )L(dx).
Now, by the Kac representation (see [START_REF] Bertail | Edgeworth expansions of suitably normalized sample mean statistics for atomic Markov chains[END_REF]), θ is the value T (L) depending on the distribution L of one block such that

E a ∂ ∂θ l(B 1 , θ) = 0.
Just like in the i.i.d. case, a simple estimator of the distribution of the blocks is given by the empirical distribution of the blocks constructed for the data, that is:

L n = (l n -1) -1 ln-1 j=1 δ B j .
It follows that the plug-in-estimator θ(L n ) is simply the value such that

ln-1 j=1 ∂ ∂θ l(B j , θ) = 0,
which is the m.l.e.

Assuming l(B i , θ) is twice dierentiable in θ and that E a [ ∂ 2 l(B i ,θ) ∂θ 2 ] is dened and positive, then under some regularity assumptions and with an adequate choice of ||.|| (we refer to [START_REF] Bertail | Bootstrapping robust statistics for markovian data applications to regenerative r-statistics and l-statistics[END_REF] for discussion of the choice of the metric), θ may be seen as a Fréchet dierentiable functional with inuence function (in terms of the i.i.d. blocks) given by:

T (1) g (b, θ) = -E a ∂ 2 l(B i , θ) ∂θ 2 -1 ∂ ∂θ l(b, θ).
Then, under any starting distribution ν, as n → ∞, we have n 1/2 ( θ n -θ) ⇒ N (0, σ 2 (θ)),

where:

σ 2 (θ) = -E a ∂ 2 l(B i , θ) ∂θ 2 -1 V ar a ( ∂ ∂θ l(B i , θ)) E a [τ 0 ] -E a ∂ 2 l(B i , θ) ∂θ 2 -1
.

Unfortunately this variance is very complicated to estimate because of the dependence on θ within each block and because the log likelihood involving blocks of random length can not be explicitly computed. Nevertheless, the following results obtained in [START_REF] Bertail | Bootstrapping robust statistics for markovian data applications to regenerative r-statistics and l-statistics[END_REF] ensure the validity of the bootstrap method for a large class of dierentiable functionals, including M and L estimators based on blocks.

Theorem 4. If T :

P T → R is Fréchet dierentiable at L and ||L n -L|| = O Pν (n -1/2 ) (or R (1) (L n , L) = o Pν (n -1/2 )) as n → ∞, and if E a [τ 0 ] < ∞ and 0 < V ar a (T (1) (B 1 , L)) < ∞,
then under P ν , we have the convergence in distribution

n 1/2 (T (L n ) -T (L)) d → N (0, E a [τ 0 ]V ar a (T (1) (B 1 , L)), as n → ∞.

If in addition ||L

* n -L n || = O P (n -1/2 ) or the remainder (or R (1) (L * n , L n ) = o Pν (n -1/2 ))as
n → +∞, then we also have:

n 1/2 (T (L * n ) -T (L n )) d * → N (0, E a [τ 0 ]V ar a (T (1) (B 1 , L)), as n → ∞.
Moreover, the bootstrap variance converges at least in P-probability to the correct one.

V ar * (n 1/2 (T (L * n ) -T (L n ))) Pr → E a [τ 0 ]V ar a (T (1) (B 1 , L)), as n → ∞.
This results justify the use of the regenerative bootstrap to obtain convergent estimators of the variance as well as condence intervals for the m.l.e estimators in the ZINAR(p) case.

Prediction in a regenerative framework.

In this section, we present and discuss some strategies to predict and evaluate the predictive capability of our model. Having observed {y 1 , . . . , y n }, we want to predict Y n+h , for some h ≥ 1.

As mentioned by [START_REF] Weiss | An Introduction to Discrete-Valued Time Series[END_REF], for real-valued processes, the most common type of point forecasting is the conditional mean, as this is known to be optimal in the sense of the mean squared error. Thus, the h-step-ahead conditional mean is given by the recurrence equations:

y n+h = E(Y n+h |Y n+h-1 = y n+h-1 , . . . , Y n+h-p+1 = y n+h-p+1 ) = p i=1 α i • y n+h-i + (1 -π)E (U h+1 ) , h ≥ 1
with the convention that y n+h-i = y n+h-i when i ≥ h.

However, the main disadvantage of the mean forecast is that it will usually lead to non-integer value predictions, while Y n+h will certainly take an integer value. A possible solution is of course to round the predictions to the closest integer value, but this will create an uncontrolled bias.

In this manuscript, we are more interested in the full distribution of the predicted value, that is:

p n+h (j) = P (Y n+h = j|Y n+h-1 = y n+h-1 , . . . , Y n+h-p+1 = y n+h-p+1 )
for j = 0, . . . , ∞ (or some large H) and the successive h ≥ 1. Thus, in order to obtain the successive h-step-ahead predictive distribution of Y n+h , h ≥ 1, we instead propose to use a mixed procedure with both a model based prediction procedure and the regenerative bootstrap ideas, as described before, to obtain condence intervals for the estimated distribution of the future observations.

Let θ = (α, π, λ) , with α = (α 1 , . . . , α p ) , the parameter vector of the ZI-INAR(p) process. We update the values of {Y n+1 , Y n+2 , . . .} one component at a time, using a procedure similar to that described by [START_REF] Neal | MCMC for integer-valued ARMA processes[END_REF] and [START_REF] Garay | Bayesian analysis of the p-order integer-valued ar process with zero-inated poisson innovations[END_REF] for predictive inference, but including an additional level of the regenerative procedure to get the distribution of the predictor:

Step 1 Repeat the following procedure B times:

(a) Generate a bootstrap sample using the regenerative bootstrap method (by resampling complete regenerative blocks).

n+h,i + v (m) n+h , (13)
for h ≥ 1 and m = 1, . . . , M .

Step 3 For b = 1, . . . , B, compute:

p b n+h (j) = # m, y (m) n+h = j M , h ≥ 1.
One can also compute the mean predictor:

y b n+h = 1 M M m=1 y (m)
n+h , h ≥ 1 but this will not necessarily be an integer.

The empirical distribution of p b n+h (j) , b = 1, . . . , B clearly yields an approximation of the distribution for the estimator of the predictive distribution of y n+h , h ≥ 1. Thus, for a xed j, we can construct a condence interval of the true probability p n+h (j).

We might construct the bootstrap distribution of the mean predicted value and obtain condence intervals for it, but it is less informative.

We choose B = 999 and M = 199, which will allow determining uniquely the (1 -γ)/2 and (1+γ)/2 quantiles of the bootstrap distributions for the usual risks γ = 90%, 95%, 99%.

The regenerative phase only relies on the hypothesis that the times series under study has a Markov chain structure and thus does not depend on our chosen model. This has the advantage of obtaining a valid approximation of the distribution of the m.l.e. even if the original model is mispecied.

Criteria for comparing models

The idea behind information criteria is to balance the goodness of t against the model size by adding an appropriate penalty term to (θ|y). Thus, in order to compare the dierent models considered in this manuscript, we use the Akaike information criterion (AIC) see [START_REF] Akaike | A new look at the statistical model identication[END_REF] and the Bayesian information criterion (BIC) see [START_REF] Schwarz | A new look at the statistical model identication[END_REF] computed based on the CML approach. However, as mentioned by [START_REF] Weiss | An Introduction to Discrete-Valued Time Series[END_REF], since the number of terms in (θ|y), dened in Eq. ( 5), varies with changes of p, one may insert the factor n/(n -p) before (θ|y) to account for this distortion. Thus, the AIC M od and BIC M od are given by:

AIC M od = -2 n (n -p) (θ|y) + 2κ and BIC M od = -2 n (n -p) (θ|y) + p × log(n),
where κ represents the number of parameters of the model.

Simulation Study

The goal of this section is to evaluate the nite-sample performance of the parameter estimates for the ZI-INAR(p) model. We generated articial samples from this model, considering p ∈ {1, 2, 3} , π ∈ {0.3, 0.6} and α = 0.3, α = (α 1 , α 2 ) = (0.3, 0.2) and α = (α 1 , α 2 , α 3 ) = (0.25, 0.2, 0.15) for ZI-INAR(1) to ZI-INAR(3) processes, respectively. The sample sizes were xed to n ∈ {100, 300, 500, 1000} and the ZI models considered for the innovation were respectively: (i) ZIP, with λ = 2; and (ii) ZINB and ZIPIG with µ = 2 and φ = 1.5. We denote the processes INAR(p) with ZIP, ZINB and ZIPIG innovations by ZIP-INAR(p), ZINB-INAR(p) and ZIPIG-INAR(p), respectively. For each of these combinations of model, p order and π, we generated N = 300 replicates of size n and analyzed the relative bias (RB) and root relative mean square error (RRMSE), given by:

RB( θ i ) = 1 N N j=1 θ ij -θ i θ i and RRMSE( θ i ) = 1 N N j=1 θ ij -θ i θ i 2 ,
where θ ij is the estimate of parameter θ i , in the j-th replicate. 4, namely: mean, standard deviation (SD) and 2.5% and 97.5% quantiles. These condence intervals aim at stressing that the probability to have one or two assaults may be very high. Indeed, for two assaults, this probability can be bigger than 5%, which may be unacceptable, and indeed the true realizations are equal to 2 in both cases (for observations #421 and #422). Even if the chance of having zero assaults is high, notice that the probability of having one assault is between 0.16 and 0.30, which is very high and should denitely be taken into account. Figure 3-4 presents the histograms of the relative frequency of the predicted value for the data points #421 and #422. These histograms show that for 0 to 3 assaults, the distribution is close to a normal one but becomes very erratic for k=4. Indeed, in that case there is a very small probability of having more than 4 assaults, and the asymptotic distribution is expected to look more like a Poisson distribution (which is the case here). innovative EM-type algorithm to obtain ML parameter estimates computationally. We also present a regenerative bootstrap method to construct condence intervals for the parameters and construct the forecasting distribution for future values. We implemented all the programs in the R software (R Development Core Team, 2019) and the program codes are available from us on request. A Monte-Carlo simulation study was conducted to investigate the performance of the model in the presence of dierent levels of zero proportions. Finally, a real data analysis demonstrated the applicability and benet of the proposed approach for practical cases, where we showed strong evidence of high-order dependence and inated zero counts. The utility of our approach can be further extended, for example, by adding a moving average structure. Another possible extension of the work is to consider a full Bayesian approach as a basis for inference and prediction, as in [START_REF] Garay | Bayesian analysis of the p-order integer-valued ar process with zero-inated poisson innovations[END_REF].

As shown by

  we use the Aitken acceleration-based stopping criterion (McLachlan and Krishnan, 2008) as a convergence rule. Dening (k+1) = θ (k+1)

  3. From the resampled data blocks, construct a pseudo-trajectory by binding the blocks together Y * (n) = (B * 1,n , . . . , B * l * n -1,n ) and truncating the joint blocks to get a time series of length n. Then recompute the value of the statistics of interest T * n = T n (Y * (n) ).

Step 4 .

 4 If S n = S(B 1 , ..., B ln-1 ) is an estimator of the variance (otherwise set it to 1) of the original statistic T n . Similarly compute S * n,bn = S(B * 1,n , . . . , B * l *

  These hypotheses are clearly satised for stationary INAR(p) processes which hit 0 (see proposition 2.1 in Drost et al. (2010)).

  (b) Estimate the parameters θ, of the model, using the EM procedure described before for this boostrap sample. These parameters are denoted by θb = α b , π b , λ b , with b = 1, . . . , B.Step 2 For b = 1, . . . , B, repeat the following steps M times: (a) Draw the value of w (m) n+h from the distribution Bern π b . (b) Draw the value of v (m) n+h from the distribution of V n+h |w

  (c) For 1 ≤ i ≤ p, obtain: s (m) n+h,i = 0, if y n+h-i = 0 Value drawn from Bin y n+h-i ; α b i , otherwise.

Figure 2 :

 2 Figure 2: Assaults dataset. Bar plot of the predictive distributions of data points #421 and #422 in ZIP-INAR(2) model.

Figure 3 :Figure 4 :

 34 Figure 3: Assaults dataset. Histograms of the relative frequency for the predicted values of data point #421 in the ZIP-INAR(2) model.

Table 1 ,

 1 for all the parameters of the ZIP-INAR(p) model, the RB and RRMSE decrease when the sample size increases indicating that the estimates based on the proposed EM-type algorithm do satisfy good asymptotic properties. The results for the cases ZINB-INAR(p) and ZIPIG-INAR(p) are similar, and are presented in Appendix B of the supplementary material.

Table 1 :

 1 RB and RRMSE (in parentheses) of the parameter estimates of the ZIP-INAR(p) process. π b and λ b are the m.l.e of the parameters, respectively, obtained using the proposed EM algorithm considering the b-th regenerative bootstrap sample.

	p π n	α 1	α 2	α 3	π	λ
		100 -0.067 (0.314)	-	-	-0.076 (0.344) -0.004 (0.116)
	0.3	300 -0.017 (0.164) 500 -0.023 (0.119)	--	--	-0.009 (0.194) -0.013 (0.137)	0.006 (0.078) 0.003 (0.058)
	1	1000 -0.002 (0.087) 100 -0.043 (0.236)	--	--	-0.015 (0.010) -0.002 (0.040) -0.021 (0.134) 0.001 (0.158)
	0.6	300 -0.012 (0.129) 500 -0.012 (0.099)	--	--	-0.007 (0.071) -0.006 (0.090) -0.005 (0.049) 0.002 (0.065)
		1000 -0.010 (0.075)	-	-	-0.003 (0.039) -0.004 (0.047)
		100 -0.060 (0.428) -0.106 (0.399)	-	-0.142 (0.476)	0.037 (0.176)
	0.3	300 -0.003 (0.240) -0.058 (0.252) 500 -0.004 (0.173) -0.018 (0.183)	--	-0.054 (0.279) -0.026 (0.210) -0.001 (0.070) 0.013 (0.096)
	2	1000 -0.010 (0.120) -0.005 (0.138) 100 -0.020 (0.332) -0.103 (0.398)	--	-0.017 (0.150) -0.056 (0.178) -0.014 (0.158) 0.006 (0.053)
	0.6	300 -0.004 (0.195) -0.032 (0.203) 500 0.004 (0.146) -0.029 (0.157)	--	-0.018 (0.097) -0.010 (0.101) -0.009 (0.076) -0.003 (0.084)
		1000 -0.009 (0.105) -0.017 (0.115)	-	-0.009 (0.048) -0.006 (0.052)
		100 -0.098 (0.580) -0.128 (0.506) -0.001 (0.529) -0.159 (0.612)	0.057 (0.218)
	0.3	300 -0.020 (0.291) -0.010 (0.305) -0.041 (0.289) -0.009 (0.349) 500 -0.012 (0.242) -0.014 (0.215) -0.044 (0.237) -0.046 (0.271)	0.032 (0.123) 0.010 (0.092)
	3	1000 -0.010 (0.168) -0.002 (0.166) -0.015 (0.151) -0.019 (0.192) 100 -0.019 (0.506) -0.099 (0.522) -0.055 (0.469) -0.094 (0.281) -0.016 (0.215) 0.007 (0.059)
	0.6	300 -0.002 (0.263) -0.043 (0.292) -0.037 (0.285) -0.020 (0.121) 500 0.003 (0.198) -0.024 (0.216) -0.035 (0.216) -0.018 (0.095) -0.001 (0.095) 0.008 (0.120)
		1000 0.008 (0.150)	-0.029 (0.144) -0.014 (0.151) -0.010 (0.067) -0.004 (0.063)
	6. Real Dataset: Assaults
	This section considers the Assaults dataset, previously analyzed by Garay et al. (2020)

from a Bayesian approach. This dataset concerns the weekly number of assaults, recorded from January 2008 to December 2015 at Federal University of Pernambuco (UFPE) -Brazil.

Originally, It has 422 observations, but here we only use the rst 420 observations to estimate the parameters and we do not consider the last two for prediction comparison purposes.

where

Table 3

 3 In order to obtain the forecasting distribution for the last two values of the dataset (observations #421 and #422), we use the strategy presenting in Section 4.3, in the ZIP-INAR(2) model. The bar plots of the predictive distributions of the data points #421 and #422 are presenting in Figure2, showing close agreement with the bar plot of the data given in Figure1.

	: m.l.e., SE-Boots and condence intervals (CI Boots) for the ZIP-INAR(2) and ZIP-INAR(3) pro-
	cesses.
	Param			ZIP-INAR(2)				ZIP-INAR(3)
		EM	Mean Boots SE Boots	CI Boots	EM	Mean Boots SE Boots 95% CI Boots
	α 1	0.181	0.180	0.035	(0.112, 0.248) 0.173	0.172	0.036	(0.102, 0.242)
	α 2	0.101	0.105	0.042	(0.023, 0.187) 0.084	0.089	0.043	(0.004, 0.173)
	α 3					0.082	0.072	0.039	[0.000, 0.148)
	π	0.517	0.516	0.057	(0.404, 0.628) 0.549	0.542	0.060	(0.424, 0.660)
	λ	1.167	1.173	0.119	(0.939, 1.407) 1.151	1.168	0.122	(0.928, 1.407)

Table 4 :

 4 Assaults dataset. Summaries of the predictive distributions of Y421 and Y422 in the ZIP-INAR(2) model.

	Predict Values		# 421			# 422
	j	Mean	SD	Q 0.025	Q 0.975	Mean	SD	Q 0.025	Q 0.975
	0	0.595 0.048	0.497	0.688	0.596 0.048	0.492	0.688
	1	0.229 0.036	0.161	0.302	0.220 0.034	0.161	0.291
	2	0.110 0.025	0.065	0.161	0.116 0.026	0.070	0.171
	3	0.047 0.017	0.020	0.080	0.048 0.017	0.020	0.085
	4	0.015 0.009	0.000	0.035	0.016 0.009	0.000	0.035
		0.004 0.005	0.000	0.015	0.004 0.005	0.000	0.015

Acknowledgements

Aldo M. Garay would like to acknowledge the support of the National Council for Scientic and Technological Development (CNPq-Brazil), Grant 420082/2016-6. This study was also conducted as part of the project Labex MME-DII (ANR11-LBX-0023-01).

As can be seen in the time series plot and bar plot, presented in panels (a) and (b) of Fig. 1 respectively, there is a signicant presence of zeros (52.4%). The partial autocorrelation function (PACF), of panel (c), indicates that a model with an order dependence structure greater than 1 should be suitable. The value of AIC M od and BIC M od for these 18 models are presented in Table 2, where we observe that the ZIP-INAR(2) and ZIP-INAR(3) processes seem to produce more accurate estimates.