Yang Lu

The Predictive Distributions of Thinning-based Count Processes

Keywords: compound autoregressive process, probabilistic forecast of counts, matrix arithmetic

This paper shows that the term structure of conditional, or predictive distributions allows for closed form expression in a large family of (possibly higher-order, or infinite order) thinning-based count processes such as INAR(p), INARCH(p), NBAR(p), and INGARCH(1,1). Such predictive distributions are currently often deemed intractable by the literature and existing approximation methods are usually time consuming and induce approximation errors. In this paper, we propose a Taylor's expansion algorithm for these predictive distributions, which is both exact and fast. Through extensive simulation exercises, we demonstrate its advantages with respect to existing methods in terms of the computational gain and/or precision.

Introduction

Nonnegative integer-valued count process models are widely used in domains such as marketing [START_REF] Böckenholt | Mixed INAR (1) Poisson Regression Models: Analyzing Heterogeneity and Serial Dependencies in Longitudinal Count Data[END_REF]], economics [START_REF] Blundell | Market Share, Market Value and Innovation in a Panel of British Manufacturing Firms[END_REF]; [START_REF] Brännäs | Generalized integer-valued autoregression[END_REF]; [START_REF] Harris | Semiparametric independence testing for time series of counts and the role of the support[END_REF]], finance [START_REF] Heinen | Multivariate Autoregressive Modeling of Time Series Count Data using Copulas[END_REF]; [START_REF] Bien | An Inflated Multivariate Integer Count Hurdle Model: An Application to Bid and Ask Quote Dynamics[END_REF]; [START_REF] Kirchner | An Estimation Procedure for the Hawkes Process[END_REF] and insurance [START_REF] Gouriéroux | Heterogeneous INAR (1) Model with Application to Car Insurance[END_REF]]. The benchmark model, introduced by [START_REF] Mckenzie | Some Simple Models for Discrete Variate Time Series[END_REF], Al-Osh and [START_REF] Al-Osh | First-order Integer-valued Autoregressive (INAR(1)) Process[END_REF] in the first-order case (called INAR(1)) postulates that:

X t = α • X t-1 + t , ∀t, (1)
where the thinning operator is defined as follows: conditionally on X t-1 , variable α•X t-1 has the distribution B(X t-1 , α), that is the binomial distribution with probability parameter α and size parameter X t-1 . Moreover, this variable is conditionally independent from the i.i.d. sequence of innovation ( t ), which follows Poisson distribution.

Since this seminal work, extensions have flourished. [START_REF] Du | The Integer-Valued Autoregressive (INAR (p)) Model[END_REF] introduce higher-order INAR(p) models satisfying:

X t = p i=1 α i • X t-i + t , ∀t, (2)
where given all the past observations X t-1 = {X t-1 , X t-2 , ...}, variables α i • X t-i are mutually independent with Binomial distribution B(X t-i , α i ), and are independent from t , which is an i.i.d. sequence, and the distribution of the latter is not necessarily Poisson. In particular, if

( t )
is Poisson (resp. negative binomial) distributed, then we say that (X t ) follows the Poisson-INAR(p) [resp. NB-INAR(p)] process. Note that [START_REF] Alzaid | An Integer-valued pth-order Autoregressive Structure (INAR (p)) Process[END_REF] has also proposed an alternative INAR(p) specification. This latter is less tractable and even the conditional expectation formula is highly complicated [see [START_REF] Jung | Coherent Forecasting in Integer Time Series Models[END_REF] for its expression in the special case where p = 2]. Thus this model is not considered in our paper. This choice is consistent with the majority of recent papers [see e.g. [START_REF] Pedeli | Likelihood Estimation for the INAR (p) Model by Saddlepoint Approximation[END_REF] and references therein].

Meanwhile, it has also been proposed to relax the binomial distributional assumption of the thinning operator [see e.g. [START_REF] Latour | Existence and stochastic structure of a non-negative integer-valued autoregressive process[END_REF]]. Roughly speaking, these models trace their origin to the applied probability literature, in which they are called Galton-Watson process with immigration, or branching process, see e.g. [START_REF] Venkataraman | A limit theorem on a subcritical galton-watson process with immigration[END_REF]. They are special cases of model:

X t = p i=1 Z i * X t-i + t , (3)
where integer p could be finite or infinite; the operator * is defined as follows. Given X t-1 , each Z i * X t-i is the sum of X t-i i.i.d. copies of count variable Z i . In other words, equation (3) does not include random coefficient count processes in which the different copies of Z i are dependent [see e.g. [START_REF] Zheng | Inference for pth-order random coefficient integer-valued autoregressive processes[END_REF]]. This distribution of Z i usually depends on some unknown parameters and is not necessarily Bernoulli, which is the case in model (2). Moreover variables Z i * X t-i are mutually independent, and are independent from the i.i.d. innovation sequence t . Prominent examples of special cases of model (3) include, among others:

• [START_REF] Ristić | A New Geometric First-order Integervalued Autoregressive (NGINAR (1)) Process[END_REF] and [START_REF] Gouriéroux | Negative Binomial Autoregressive Process with Stochastic Intensity[END_REF] consider first order models (p = 1), and assume that Z 1 is negative binomial (NB) distributed. Then they make special assumptions on the distribution of t to ensure that process (X t ) is stationary with NB stationary distribution.

• [START_REF] Rydberg | A modelling framework for the prices and times of trades made on the new york stock exchange[END_REF]; [START_REF] Weiß | Inarch (1) processes: Higher-order moments and jumps[END_REF] consider the INARCH(1) model in which the conditional distribution of X t given X t-1 is Poisson P(αX t-1 + β), with parameters α, β > 0. This corresponds to model (3) with p = 1, and both Z 1 and t are Poisson distributed with parameters α and β, respectively.

• Zhu and Joe (2010) also focuses on first-order models but assume that Z 1 and t follow other non-NB and non-Binomial distributions.

• Finally, a few models can be written in the form (3) with p = ∞. For instance, Kirchner In this paper we call models of type (3) thinning-based [see [START_REF] Weiß | Thinning operations for modeling time series of countsâa survey[END_REF]; [START_REF] Ristić | A New Geometric First-order Integervalued Autoregressive (NGINAR (1)) Process[END_REF]; [START_REF] Zhu | Negative Binomial Time Series Models Based on Expectation Thinning Operators[END_REF]; [START_REF] Scotto | Thinning-Based Models in the Analysis of Integer-Valued Time Series: A Review[END_REF]; [START_REF] Kirchner | Hawkes and INAR (∞) Processes[END_REF]]. Note that some papers in the literature also use the term "generalized thinning". We choose to avoid this latter terminology since we want to emphasize the our framework does not apply to models with random coefficient thinning. Family (3) is significantly larger than (2), but up to now the literature has only been able to handle first-order models, with quite restrictive parametric assumptions on Z i 's and t .

One explanation is that, up to now extensions (2) and ( 3) have the serious downside that the term structure of predictive probability mass functions (p.m.f.'s) is believed to be intractable, rendering estimation and/or forecasting cumbersome. This difficulty is especially acute when the order p, and/or the horizon of forecasting h increases. To our knowledge, the term structure of p.m.f.'s alows for closed form expression only in two special first-order models. This paper solves the above open difficulty for model (3). We show that the term structure of conditional p.m.f.'s in model (3) has closed form expression at any horizon, and hence significantly improve the existing estimation and forecasting procedures. The method is general in the sense that it does not require the distributions of Z i 's and t to belong to the classical families such as Bernoulli, Poisson or negative binomial. Rather, it is applicable so long as some mild assumption (see Assumption 1) on the log p.g.f. of these variables is satisfied. Our method explores the Compound Autoregressive (CaR) property of process (3), and the natural link between the probability generating function (p.g.f.) and the corresponding p.m.f. of a count distribution.

This approach has the distinctive advantage that the algorithm only involves low-dimensional matrix operations, whose computational cost is insensitive to p and remains very low even for large p.

The paper is organized as follows. Section 2 computes the Taylor's expansion of the onestep-ahead conditional p.g.f. and deduces the corresponding p.m.f., which facilitates, among others, the likelihood-based estimation. Section 3 derives the multiple-step-ahead predictive p.m.f.'s for finite-order models, before extending the result to some infinite order models such as

Poisson-INGARCH(1,1). Section 4 compares our approach and existing methods, first in terms of maximum likelihood estimation, then in terms of higher-horizon forecasting. In particular, for the estmation we show that depending on model parameters, our model is generically either comparable, or faster than the other, exact method of Bu et al. (2008). For the forecasting, we show that our method is more precise and much less time consuming than competing methods proposed by [START_REF] Mccabe | Efficient Probabilistic Forecasts for Counts[END_REF] and [START_REF] Jung | Coherent Forecasting in Integer Time Series Models[END_REF]. Section 5 concludes. Technical details and proofs are gathered in Appendices.

Likelihood-based estimation

Conducting maximum likelihood (ML) estimation requires the evaluation of the conditional p.m.f. p(X t |X t-1 ). In model (3), by standard convolution we get [see e.g. [START_REF] Drost | Efficient Estimation of Auto-Regression Parameters and Innovation Distributions for Semiparametric Integer-Valued AR (p) models[END_REF]]:

p(X t |X t-1 ) = n1+n2+•••+np+np+1=Xt P[ t = n p+1 ] p i=1 P[Z i * X t-i = n i |X t-i ]. (4)
For instance, in INAR(p) models, Z i 's are Bernoulli distributed and we have P

[Z i * X t-i = n i |X t-i ] = Xt-i ni α ni i (1 -α i ) Xt-i-ni .
Except in a few special cases such as INARCH(p) where each Z i and t follows Poisson distribution, the above p-dimensional summation cannot be further simplified, and becomes increasingly cumbersome for large p. Bu et al. (2008) point out that, this computation can be partially simplified by using the fact that the convolution of p + 1 variables is the convolution between the first variable and that of the p remaining ones, which can be further decomposed recursively. However for large p, the recursive method remains computationally heavy [see also [START_REF] Pedeli | Likelihood Estimation for the INAR (p) Model by Saddlepoint Approximation[END_REF] for discussions]. While there exist convenient and consistent moment based estimators [see e.g. Al-Osh and [START_REF] Al-Osh | First-order Integer-valued Autoregressive (INAR(1)) Process[END_REF]], they typically suffer from significant efficiency loss [see e.g. Bu et al. (2008)]. Recently, [START_REF] Pedeli | Likelihood Estimation for the INAR (p) Model by Saddlepoint Approximation[END_REF] propose a saddle-point approximation of p(X t |X t-1 ) for ML estimation. Although this method is faster than the brute force convolution, it still has some drawbacks. First, the approximated p.m.f. p(•|X t-1 ) does not sum up to one across N = {0, 1, 2, ...}. Second, the method only applies to INAR(p) models. Third, this approximation itself lacks closed form expression and its computation requires numerical inversion of some nonlinear functions. This leads to extra computational complexity, as well as approximation error. For instance, [START_REF] Pedeli | Likelihood Estimation for the INAR (p) Model by Saddlepoint Approximation[END_REF] report that for certain parameter values, the relative error of the likelihood function can be as large as 2 to 5 percent. As a consequence, this approximate ML estimator is inconsistent and simulation results show that in finite sample, the bias can be significantly larger than that of the exact ML estimator (see their Table 1).

The Taylor's expansion based algorithm

Our solution is based on the property that the one-step-ahead conditional p.g.f. of model (3)

satisfies:

E[u Xt |X t-1 ] = exp p i=1 X t-i log(E[u Zi ]) + log(E[u t ]) , ( 5
)
which is exponential affine in X t-1 . Such processes (X t ) are called (p-th order) Compound Autoregressive (CaR(p)) processes [see [START_REF] Darolles | Structural Laplace Transform and Compound Autoregressive Models[END_REF]] and by Proposition 7 of [START_REF] Darolles | Structural Laplace Transform and Compound Autoregressive Models[END_REF], we can first deduce that such a process is strictly stationary if and only if:

p i=1 E[Z i ] < 1.
Then we remark that this conditional p.g.f. is linked to the conditional p.m.f. through:

E[u Xt |X t-1 ] = ∞ k=0 P[X t = k|X t-1 ]u k . ( 6
)
Thus to get P[X t = k|X t-1 ] for all k = 0, 1, ..., n, where n is a fixed integer, it suffices to compute the n-th order Taylor's expansion at u = 0 for the corresponding conditional p.g.f. To this end we make the following assumption:

Assumption 1. In model (3), the log p.g.f. of variables Z i , i = 1, ..., p, and t all have tractable

Taylor's expansions:

g i (u) := log E[u Zi ] = ∞ k=0 c i,k u k , ∀i = 1, ..., p, ∀u, (7)
g 0 (u) := log E[u t ] = ∞ k=0 c 0,k u k , ∀u. (8)
By "tractable", we require the coefficients c i,k and c 0,k , where k = 0, ..., and i = 1, ..., p, to be analytic functions of the parameters characterizing the distributions of Z i 's and of ( t ). This assumption is quite mild and is satisfied by many standard count distributions. For instance, distributions in Examples 1-4 below have tractable p.g.f. and are infinitely differentiable at u = 0, such that Assumption 1 automatically holds.

Example 1 (Bernoulli and Poisson distributions). If X follows Bernoulli distribution B(1, α), then:

log E[u X ] = log(αu + 1 -α) = log(1 -α) + log(1 + α 1 -α u) = log(1 -α) + ∞ k=1 (-1) k-1 k α 1 -α k u k .
On the other hand, if X follows Poisson distribution P(λ), then:

log E[u X ] = -λ + λu.
Thus, in particular, Assumption 1 is satisfied by the Poisson-INAR(p) model, that is equation ( 2) with Poisson distributed innovation ( t ). In this case, the coefficients in the Taylor's expansions ( 7) and ( 8) are given by: c 0,0 = -λ, c 0,1 = λ, and c 0,j = 0, ∀j ≥ 2.

c i,0 = log(1 -α i ), and
c i,k = (-1) k α i 1 -α i k , ∀k ≥ 1, i = 1, ..., p.
Example 2 (NB distribution (including geometric distribution), see [START_REF] Ristić | A New Geometric First-order Integervalued Autoregressive (NGINAR (1)) Process[END_REF]; [START_REF] Gouriéroux | Negative Binomial Autoregressive Process with Stochastic Intensity[END_REF]). If X follows NB distribution with probability parameter q ∈ (0, 1) and size parameter r [in the following we will use the notation X ∼ N B(r, q)], then:

E[u X ] = (1 -q) r (1 -qu) r , ∀u ∈ [0, 1 q ), (9)
hence log E[u X ] = r log(1 -q) -r log(1 -qu) = r log(1 -q) + r ∞ j=1 q j j u j .
Example 3 (Compound Poisson distributions, see [START_REF] Schweer | Compound Poisson INAR (1) Processes: Stochastic Properties and Testing for Overdispersion[END_REF]). We say that variable X follows compound Poisson distribution if it has the representation:

X = (d) N j=1 Z j ,
where symbol = (d) denotes equality in distribution; variable N is Poisson P(λ) distributed, and is independent from the sequence (Z j ) j , where (Z j ) j is an i.i.d. sequence of count variables. It is easily checked that the log p.g.f. of X is:

log E[u X ] = λ E[u Z1 ] -1 = λ ∞ k=0 P[Z 1 = k]u k -1 .
Thus log E[u X ] has tractable Taylor's expansion, so long as the p.m.f. of Z 1 is tractable.

Example 4 (Distributions defined through the p.g.f.). [START_REF] Zhu | Negative Binomial Time Series Models Based on Expectation Thinning Operators[END_REF] propose several firstorder, thinning-based count processes in which Z 1 has simple p.g.f. but complicated p.m.f. For instance, the non central NB distribution has the p.g.f.:

E[u X ] = (1 -α) + (α -γ)u 1 -αγ -(1 -α)γu , ∀u ∈ [0, 1 -αγ (1 -α)γ ),
where parameters α, γ lie between 0 and 1. For this distribution, the log p.g.f. is easily expanded:

log E[u X ] = log 1 -α 1 -αγ + log(1 + α -γ 1 -α u) -log[1 - (1 -α)γ 1 -αγ u] = log 1 -α 1 -αγ + ∞ k=1 (-1) k-1 k α -γ 1 -α k u k + ∞ j=1 (1 -α)γ 1 -αγ k u k .
Another example is the Hermite distribution, which is used by [START_REF] Fernández-Fontelo | Integer-valued AR processes with Hermite innovations and time-varying parameters: An application to bovine fallen stock surveillance at a local scale[END_REF] in INAR models as the distribution of the innovation t . Such a distribution has a log. p.g.f.:

log E[u t ] = a 1 (u -1) + a 2 (u 2 -1),
where parameters a 1 , a 2 ≥ 0.

Finally, the majority of other distributions considered in the count process literature who do not have simple p.g.f. (or log p.g.f.) usually have simple p.m.f. For instance, [START_REF] Heinen | Multivariate Autoregressive Modeling of Time Series Count Data using Copulas[END_REF] introduce the double Poisson distribution, Jazi et al. ( 2012) consider zero-inflated mixture distributions. In Appendix 3 we show how the log p.g.f. can be easily expanded at any given order, so long as the p.m.f. is tractable.

Under Assumption 1, we can rewrite equation ( 5) into:

E[u Xt |X t-1 ] (10) = exp c 0,0 + p i=1 X t-i c i,0 exp n j=1 (c 0,j + p i=1 X t-i c i,j )u j + O(u n+1 ) (11) = exp(A 0 ) exp n j=1 A j u j + O(u n+1 ) (12) = exp(A 0 ) n k=0 1 k! n j=1 A j u j k + O(u n+1 ), (13)
where coefficients A j are given by:

A j = c 0,j + p i=1 X t-i c i,j , ∀j = 0, ..., n. ( 14
)
From equation ( 10) to ( 11), we have expanded the log p.g.f. of Z i 's and t ; from equation ( 12)

to ( 13) we have expanded the exponential function.

Thus to deduce the coefficient in front of the term u n in expansion ( 13), it suffices to compute recursively the n + 1 first terms of the polynomial

n j=1 A j u i k
for each k. While this can be achieved using a symbolic calculation package such as Mathematica (see [START_REF] Gordy | Saddlepoint Approximation of CreditRisk +[END_REF] for a similar application to credit risk), we suggest below a matrix-based algorithm that is simpler and more suitable for statistical packages.

Proposition 1. For any coefficients γ i , i = 0, ..., n, the (n + 1) first coefficients of polynomial

n j=0 γ j u j k
, where k = 0, ...n, are given by the column vector:

            γ 0 0 0 • • • 0 γ 1 γ 0 0 • • • 0 γ 2 γ 1 γ 0 • • • 0 • • • • • • • • • • • • • • • γ n γ n-1 • • • γ 1 γ 0             k             1 0 0 • • • 0             , ∀k = 0, ..., n. ( 15
)
The proof is obvious and omitted. In equation ( 13) we have γ 0 = 0, thus the square matrix above, called the Toeplitz matrix, becomes triangular inferior. Thus we get:

            P[X t = 0|X t-1 ] P[X t = 1|X t-1 ] • • • P[X t = n -1|X t-1 ] P[X t = n|X t-1 ]             = exp(A 0 ) n k=0 1 k!             0 0 0 • • • 0 A 1 0 0 • • • 0 A 2 A 1 0 • • • 0 • • • • • • • • • • • • • • • A n A n-1 • • • A 1 0             k             1 0 0 • • • 0             . ( 16
)
Note that although the above algorithm involves a Taylor's expansion (of the conditional p.g.f.), it provides the exact value of the conditional p.m.f, because of the identity (6) (see also

Example 5 for further discussions). In other words, if we were to compute, say,

P[X t = 5|X t-1 ],
then an expansion of order 5 or 20 would lead to the same result.

On the contrary to formula (4), which becomes increasingly cumbersome when p increases, the computational cost of ( 15) is independent of p. Indeed, it only involves square matrices and vectors of dimension (n + 1), where n does not depend on p. The order p only impacts the computation of the scalars A 0 , ..., A n , whose cost is negligible. Thus this method is extremely convenient, especially for higher-order models.

One referee rightly points out that the computational cost of our algorithm is still (inevitably) increasing in n. Thus, if the values of the counts we deal with are very large, then this method can also become cumbersome. This downside is nevertheless partially mitigated by the fact that the successive powers of the Toeplitz matrix are still Toeplitz and thus can be computed very efficiently. More precisely, we can check that:

            0 0 0 • • • 0 γ 1 0 0 • • • 0 γ 2 γ 1 0 • • • 0 • • • • • • • • • • • • • • • γ n γ n-1 • • • γ 1 0             k =             0 0 0 • • • 0 γ 1,k 0 0 • • • 0 γ 2,k γ 1,k 0 • • • 0 • • • • • • • • • • • • • • • γ n,k γ n-1,k • • • γ 1 0             (17)
for any integer k, where the sequence (γ i,k ) satisfies the recursion:

γ 1,k+1 = 0 γ 2,k+1 = γ 1,k γ 1 • • • = • • • γ n,k+1 = γ n-1,k γ 1 + γ n-2,k γ 2 + • • • + γ 1,k γ n-1 .
In other words, the multiplication between Toeplitz matrices involves at most the computation of n coefficients for a total of O( 1 2 n 2 ) operations, where the symbol O(•) should be understood as being equivalent to 1 2 n 2 when n goes to infinity. . As a comparison, ordinary matrix multiplications usually requires a cost of O(n 3 ) operations. As a consequence, in total, computing the conditional p.m.f. of X t given X t-1 involves roughly O( 1 2 X 3 t ) operations since we have to raise the Toeplitz matrix up to power X t .

To summarize, the Taylor's expansion approach makes the ML estimation feasible for higherorder models without inducing bias. Alternatively, it also facilitates Bayesian inference [see [START_REF] Mccabe | Bayesian Predictions of Low Count Time Series[END_REF]; Bisaglia and Canale ( 2016)] of the same models, the ease of which usually depending (among others) on the tractability of the likelihood function.

Examples

As an illustration, let us derive the conditional p.m.f. for some popular specifications of the distributions of Z i and t .

Example 5 (INARCH(p) model). In this example we assume that, in equation (3), variables Z 1 , ..., Z p and t follow Poisson distributions with parameters λ 1 , ..., λ p and λ 0 , respectively.

This can be viewed as the higher-order extension of the (Poisson-)INARCH(1) model studied by [START_REF] Rydberg | A modelling framework for the prices and times of trades made on the new york stock exchange[END_REF]; [START_REF] Weiß | Inarch (1) processes: Higher-order moments and jumps[END_REF]. For this model, the one-step-ahead conditional distribution is Poisson with parameter:

λ t := λ 0 + p i=1 X t-i λ i .
While it is possible to directly deduce the conditional p.m.f.:

P[X t = k|X t-1 ] = e -λt λ k t k! , ( 18
)
to support our claim that equation ( 16) provides an exact expression [see the comment below ( 16)], let us now check that we can recover equation ( 18) from ( 16). We first remark that for this model, A 0 , A 1 , ..., A n defined in ( 14) are given by:

A 0 = -λ 0 - p i=1 X t-i λ i = -λ t , A 1 = λ t , and A 2 = A 3 = • • • = 0.
Then by some simple algebra, we can rewrite the RHS of equation ( 16) into:

exp(-λ t ) n k=0 λ k t k M             0 0 0 • • • 0 1 0 0 • • • 0 0 1 0 • • • 0 • • • • • • • • • • • • • • • 0 0 • • • 1 0             k             1 0 0 • • • 0             = exp(-λ t )             1 λ t λ 2 t 2! • • • λ n t n!             .
Thus equation ( 16) provides indeed the exact expression of the conditional p.m.f. 2), under the assumption that t is Poisson P(λ) distributed. The conditional p.g.f. is:

E[u Xt |X t-1 ] = exp p i=1 X t-i log(α i u + 1 -α i ) + λ(u -1) = exp(A 0 ) n k=0 1 k! n j=1 A j u j k + O(u n+1 ), with A 0 = -λ + p i=1 X t-i log(1 -α i ), A 1 = λ + p i=1 X t-i α i 1 -α i ,
and

A j = (-1) j-1 j p i=1 X t-i ( α i 1 -α i ) j , ∀j = 2, ..., n.
Example 7 (NB-INAR(p) model). Let us now give the corresponding formulas, when the innovation term in the INAR(p) model ( 2) follows the N B(r, q) distribution instead of the Poisson distribution. By (5), the conditional p.g.f. is:

E[u Xt |X t-1 ] = exp p i=1 X t-i log(α i u + 1 -α i ) + r log(1 -q) -r log(1 -qu) = exp r log(1 -q) + p i=1 X t-i log(1 -α i ) n k=0 1 k! n j=1 A j u j k + O(u n+1 ),
where

A 0 = r log(1 -q) + p i=1 X t-i log(1 -α i ), A j = rq j j + (-1) j-1 j p i=1 X t-i ( α i 1 -α i ) j , ∀j = 1, ..., n.
Example 8 (NBAR(p) model). Let us finally assume that in equation ( 3), variables t and Z i , i = 1, ..., p follow N B(r 0 , q 0 ) and N B(r i , q i ) distributions, respectively. This model can be viewed as a higher-order extension of the Negative Binomial Autoregressive [NBAR(1)] process considered in [START_REF] Ristić | A New Geometric First-order Integervalued Autoregressive (NGINAR (1)) Process[END_REF]; [START_REF] Gouriéroux | Negative Binomial Autoregressive Process with Stochastic Intensity[END_REF] and thus will be called NBAR(p).

Its conditional p.g.f. is:

E[u Xt |X t-1 ] = exp p i=1 X t-i r i [log(1 -q i ) -log(1 -q i u)] + r 0 [log(1 -q 0 ) -log(1 -q 0 u)] = exp(A 0 ) n k=0 1 k! n j=1 A j u j k + O(u n+1 ), with A 0 = r 0 log(1 -q 0 ) + p i=1 X t-i r i log(1 -q i ), A j = r 0 q j 0 j + p i=1 X t-i r i q j i j , ∀j = 1, ..., n.
Finally, in all the above four examples, once the coefficients A j are computed, it suffices to insert them into (16) to get the vector of conditional p.m.f.

Multi-step forecasting

One of the major challenges of count process models is that even though the predictive means

E[X t+h |X t-
p h (X t+h |X t-1 ) = ∞ X t+h-1 =0 ∞ X t+h-2 =0 • • • ∞ Xt=0 p(X t+h |X t+h-1 )p(X t+h-1 |X t+h-2 ) • • • p(X t |X t-1 ). (19)
This is a h-dimensional, infinite summation, and thus is extremely cumbersome for any h ≥ 1, even if the previous section has provided a simpler method of computing the individual terms of the form p(X t |X t-1 ).

So far there are two major solutions in the literature. [START_REF] Jung | Coherent Forecasting in Integer Time Series Models[END_REF] the probability of the process taking values larger than some integer threshold τ by 0. Then the

p-dimensional vector Y t = (X t , X t-1 , ..., X t-p+1
) is approximated by a first-order, (τ + 1) p states Markov chain, and the approximate p h (•|X t-1 ) satisfies a recursive formula (in h) involving the (τ + 1) p × (τ + 1) p transition matrix Π of the chain (Y t ). Its drawbacks is that, when τ and/or p are not small enough, the dimension of Π is too high for the computation to be feasible. This necessity to consider, for some applications, a larger threshold has been highlighted by [START_REF] Zhu | Modelling Heavy-Tailed Count Data Using a Generalised Poissoninverse Gaussian Family[END_REF] and references therein.

The algorithm

Let us now propose a simpler algorithm using the same, Taylor's expansion approach. As in Section 2, we first derive the conditional p.g.f. E[u X t+h |X t-1 ]. One nice property of the CaR count processes is that such higher-order p.g.f. is still exponential affine. In other words, the following proposition is a higher-order extension of equation ( 4).

Proposition 2. For any integer h ≥ 0, we have:

E[u X t+h-1 |X t-1 ] = exp B h,0 (u) + p i=1 B h,i (u)X t-i , ( 20
)
where functional (in u) coefficients B h,i (u), i = 0, ..., p, satisfy the recursive formula (in h):

B 1,0 (u) = log E[u t ], (21)
B 1,i (u) = log E[u Zi ], ∀i = 1, ..., p, ( 22
)
B h+1,0 (u) = B(h, 0)(u) + g 0 (e B h,1 (u) ), (23)
B h+1,i (u) = 1 i<p B h,i+1 (u) + g i (e B h,1 (u) ), ∀i = 1, ..., p. ( 24
)
Proof. See Appendix 1.

Although B h,i (u) is typically intractable for large h, their Taylor's expansions up to a finite order (at u = 0), and hence that of E[u X t+h |X t-1 ], can be easily obtained by recursion, and we will see later on that this is sufficient for probabilistic forecasting purpose. Indeed, given a large integer threshold τ and suppose that the τ -th order Taylor's expansions of B h,i (u) have already been obtained for each i = 0, ..., p, then those of B h+1,i (u) are available through:

Proposition 3. We have:

B h+1,0 (u) = B h,0 (u) + c 0,0 + τ k=0 B k h,1 (u) k! τ j=1 c 0,j j k + O(u τ +1 ), B h+1,i (u) = 1 i<p B h,i+1 (u) + c i,0 + τ k=0 B k h,1 (u) k! τ j=1 c i,j j k + O(u τ +1 ), i = 1, ..., p.
Proof. See Appendix 2.

Thus it suffices to apply Proposition 1 onceto get the τ -th order Taylor's expansions of the successive powers of B h,1 (u), to derive that of B h+1,i (u), i = 0, ..., p. More precisely, if we write the Taylor's expansion of B(h, i)(u) as:

B h,i (u) = τ k=0 b h,i,k u k + O(u τ +1 ),
then for each i = 0, ..., p, the vector of coefficients

(b h,i,0 , b h,i,1 , • • • , b h,i,τ ) satisfies the recursion:          b h+1,i,0 b h+1,i,1 • • • b h+1,i,τ          = 1 i=0          b h,0,0 b h,0,1 • • • b h,0,τ          + 1 0<i<p          b h,i+1,0 b h,i+1,1 • • • b h,i+1,τ          +          c i,0 0 • • • 0          + τ k=0 1 k! τ j=1 c i,j j k          b h,i,0 0 • • • 0 b h,i,1 b h,i,0 • • • 0 • • • • • • • • • • • • b h,i,τ b h,i,τ -1 • • • b h,i,0          k          1 0 • • • 0         
.

Finally, we can re-arrange equation ( 20) into:

E[u X t+h+1 |X t ] = exp A h+1,0 + τ j=1 A h+1,j u j + O(u τ +1 ) = exp(A h+1,0 ) τ k=0 1 k! τ j=1 A h+1,j u j k + O(u τ +1 ), (25)
where each A h+1,j , j = 0...τ is linear in X t-1 :

A h+1,j = b h+1,0,j + p i=1 b h,i,j X t-i . ( 26
)
Thus the Taylor's expansion of the RHS of equation ( 25) can be obtained by applying Proposition 1 again and the values of P[X t+h+1 = j|X t-1 ], j = 0, ...τ, are obtained by coefficient matching:

            p h+1 (0|X t-1 ) p h+1 (1|X t-1 ) • • • p h+1 (τ -1|X t-1 ) p h+1 (τ |X t-1 )             = τ k=0 exp(A h+1,0 ) k!             0 0 0 • • • 0 A h+1,1 0 0 • • • 0 A h+1,2 A h+1,1 0 • • • 0 • • • • • • • • • • • • • • • A h+1,τ A h+1,τ -1 • • • A h+1,1 0             k             1 0 0 • • • 0             . ( 27
)
In terms of the computational effort, the Taylor's expansions of B h,i (u) do not depend on the values of X t and thus only need to be computed once in an online forecasting framework. They necessitate only the multiplication between matrices and vectors of dimension (τ + 1), instead of (τ + 1) p in the Markov chain approach. Thus the method is expected to be much faster than the Markov chain approach, especially for large p. Note also that although the computation of the Taylor's expansion of B h,i is recursive, that of the h-step-ahead conditional p.m.f. does not necessitate a recursion. That is, if one is interested in the forecasting performance at a specific horizon H [see e.g. [START_REF] Chevillon | Direct multi-step estimation and forecasting[END_REF]], then she can use the above formula directly once B H,i , i = 0, ..., p are obtained, without computing the conditional p.m.f. at shorter horizons.

Examples

Let us now give the explicit forecasting formula in the four examples considered in Section 2.2.

Example 9 (INARCH(p) continued.). We first consider the INARCH(p) model introduced in

Example 5. We have:

B 1,0 (u) = λ 0 (u -1), B 1,i (u) = λ i (u -1), ∀i = 1, ..., p, B h+1,0 (u) = B(h, 0)(u) + λ 0 [e B h,1 (u) -1], (28)
B h+1,i (u) = 1 i<p B h,i+1 (u) + λ i (e B h,1 (u) -1), ∀i = 1, ..., p. (29)
By expanding ( 28) and ( 29) at u = 0, we get:

         b h+1,0,0 b h+1,0,1 • • • b h+1,0,τ          = 1 i=0          b h,0,0 b h,0,1 • • • b h,0,τ          +          -λ 0 0 • • • 0          + λ 0 τ k=0 1 k!          b h,i,0 0 • • • 0 b h,i,1 b h,i,0 • • • 0 • • • • • • • • • • • • b h,i,τ b h,i,τ -1 • • • b h,i,0          k          1 0 • • • 0          ,          b h+1,i,0 b h+1,i,1 • • • b h+1,i,τ          = 1 0<i<p          b h,i+1,0 b h,i+1,1 • • • b h,i+1,τ          +          -λ i 0 • • • 0          + τ k=0 λ i k!          b h,i,0 0 • • • 0 b h,i,1 b h,i,0 • • • 0 • • • • • • • • • • • • b h,i,τ b h,i,τ -1 • • • b h,i,0          k          1 0 • • • 0         
.

Example 10 (Poisson-INAR(p) continued). Let us now consider the Poisson-INAR(p) model

introduced in Example 6. First, the recursion for functions B h,i is:

B 1,0 (u) = λ(u -1), B 1,i (u) = log(α i u + 1 -α i ), ∀i = 1, ..., p, B h+1,0 (u) = B h,0 (u) + λ[e B h,1 (u) -1], B h+1,i (u) = 1 i<p B h,i+1 (u) + log(α i e B h,1 (u) + 1 -α i ), ∀i = 1, ..., p.
their Taylor's expansions at zero are:

         b h+1,0,0 b h+1,0,1 • • • b h+1,0,τ          = 1 i=0          b h,0,0 b h,0,1 • • • b h,0,τ          +          -λ 0 • • • 0          + λ τ k=0 1 k!          b h,i,0 0 • • • 0 b h,i,1 b h,i,0 • • • 0 • • • • • • • • • • • • b h,i,τ b h,i,τ -1 • • • b h,i,0          k          1 0 • • • 0          ,          b h+1,i,0 b h+1,i,1 • • • b h+1,i,τ          = 1 0<i<p          b h,i+1,0 b h,i+1,1 • • • b h,i+1,τ          +          log(1 -α i ) 0 • • • 0          + τ k=0 1 k! τ j=1 (-1) j-1 α j i (1 -α i ) j j k-1          b h,i,0 0 • • • 0 b h,i,1 b h,i,0 • • • 0 • • • • • • • • • • • • b h,i,τ b h,i,τ -1 • • • b h,i,0          k          1 0 • • • 0         
Example 11 (NB-INAR(p) continued). Let us now give the formula in the NB-INAR(p) model considered in Example 7. First, the recursion for functions B h,i is:

B 1,0 (u) = r 0 log(1 -q 0 ) -r 0 log(1 -q 0 u), B 1,i (u) = log(α i u + 1 -α i ), ∀i = 1, ..., p, B h+1,0 (u) = B h,0 (u) + r 0 log(1 -q 0 ) -r 0 log(1 -q 0 e B h,1 (u) ), B h+1,i (u) = 1 i<p B h,i+1 (u) + log(α i e B h,1 (u) + 1 -α i ), ∀i = 1, ..., p.
Their Taylor's expansions at zero are:

         b h+1,0,0 b h+1,0,1 • • • b h+1,0,τ          = 1 i=0          b h,0,0 b h,0,1 • • • b h,0,τ          +          r 0 log(1 -q 0 ) 0 • • • 0          + r 0 τ k=0 1 k! τ j=1 q j 0 j k-1          b h,i,0 0 • • • 0 b h,i,1 b h,i,0 • • • 0 • • • • • • • • • • • • b h,i,τ b h,i,τ -1 • • • b h,i,0          k          1 0 • • • 0                   b h+1,i,0 b h+1,i,1 • • • b h+1,i,τ          = 1 0<i<p          b h,i+1,0 b h,i+1,1 • • • b h,i+1,τ          +          r i log(1 -q i ) 0 • • • 0          + τ k=0 1 k! τ j=1 (-1) j-1 α j i (1 -α i ) j j k-1          b h,i,0 0 • • • 0 b h,i,1 b h,i,0 • • • 0 • • • • • • • • • • • • b h,i,τ b h,i,τ -1 • • • b h,i,0          k          1 0 • • • 0          .
Example 12 (NBAR(p) continued). Finally, in the Poisson-INAR(p) model, the recursion for functions B h,i is:

B 1,0 (u) = r 0 log(1 -q 0 ) -r 0 log(1 -q 0 u), B 1,i (u) = r i log(1 -q i ) -r i log(1 -q i u), ∀i = 1, ..., p, B h+1,0 (u) = B h,0 (u) + r 0 log(1 -q 0 ) -r 0 log(1 -q 0 e B h,1 (u) ), B h+1,i (u) = 1 i<p B h,i+1 (u) + log(α i e B h,1 (u) + 1 -α i ), ∀i = 1, ..., p.,
whereas their Taylor's expansion at zero are:

         b h+1,0,0 b h+1,0,1 • • • b h+1,0,τ          = 1 i=0          b h,0,0 b h,0,1 • • • b h,0,τ          +          r 0 log(1 -q 0 ) 0 • • • 0          + r 0 τ k=0 1 k! τ j=1 q j 0 j k-1          b h,i,0 0 • • • 0 b h,i,1 b h,i,0 • • • 0 • • • • • • • • • • • • b h,i,τ b h,i,τ -1 • • • b h,i,0          k          1 0 • • • 0                   b h+1,i,0 b h+1,i,1 • • • b h+1,i,τ          = 1 0<i<p          b h,i+1,0 b h,i+1,1 • • • b h,i+1,τ          +          r i log(1 -q i ) 0 • • • 0          + τ k=0 1 k! τ j=1 r i q j i j k-1          b h,i,0 0 • • • 0 b h,i,1 b h,i,0 • • • 0 • • • • • • • • • • • • b h,i,τ b h,i,τ -1 • • • b h,i,0          k          1 0 • • • 0          .
Finally, the conditional p.m.f. at horizon h is given by equation ( 27) with the same definition of coefficients A h+1,0 as in ( 26).

The case of some infinite order models

Although our paper is mainly interested in models with finite order p, it can be straightforwardly extended to some infinite order, GARCH-type models. For instance, the Poisson-INGARCH(1,1) model of [START_REF] Ferland | Integer-Valued GARCH Process[END_REF] postulates that:

X t |X t-1 ∼ P(λ t ), where λ t = ρλ t-1 + λ 0 + aX t-1 , ( 30
)
with a > 0, λ 0 > 0 and ρ ∈ [0, 1[. This model can be viewed as an infinite extension of model

(3), since by the infinite divisibility of the Poisson distribution, we have the representation:

X t = ∞ i=1 Z i * X t-i + t ( 31
)
where Z i , i = 1, ... and t follow Poisson distributions with parameters aρ i-1 and λ0 1-ρ , respectively.

The above example can be extended by replacing the Poisson distribution by some other infinitely divisible distributions [see discussions in [START_REF] Gonçalves | Infinitely divisible distributions in integervalued garch models[END_REF]]. For instance, Xu et al. ( 2012) define a NB-INGARCH(1,1) process by:

X t |X t-1 ∼ N B(r t , q),
where

r t = ρr t-1 + r 0 + aX t-1 , ( 32
)
and this process still has the infinite thinning representation (31) with Z i and t following N B(aρ i-1 , q) and N B( r0 1-ρ , q) distributions, respectively. Note that this definition is different from that given by [START_REF] Zhu | A negative binomial integer-valued garch model[END_REF], which assumes instead that the conditional negative binomial distribution has a fixed degree of freedom parameter r but a time-varying probability parameter q t . In particular, Zhu's model does not admit the thinning representation and its multi-step ahead conditional distribution is unknown.

The advantage of models of type ( 30) and ( 32) is that the one-step-ahead conditional distribution is tractable. Nevertheless, for horizon h ≥ 2, the predictive distribution is (up to now) unknown. Let us now explain how the approach developed in our paper can be adapted to these models.

A first, universal solution to tackle infinite order models is to approximate an infinite order model by a finite order one, by truncating the summation in (31) at a large, deterministic order p [see [START_REF] Kirchner | Hawkes and INAR (∞) Processes[END_REF] for a formal discussion of the Poisson-INGARCH(1,1) case]. Since the computational cost of our approach is rather insensitive to the lag p, this method would provide an excellent approximation, especially when the number of observations is relatively small and the infinite order model should anyway be truncated for estimation. The advantage of this method is that it also applies to any other infinite order thinning-based models with, for instance, Bernoulli distributed variables Z i , for which the one-step-ahead conditional p.m.f. need not necessarily be simple [this has motivated the assumption of Poisson distribution on Z i 's in [START_REF] Kirchner | Hawkes and INAR (∞) Processes[END_REF]].

However, given the popularity of INGARCH type models in the applied literature, it seems beneficial if an exact forecasting method could be derived. This is indeed the case and is due to the property that in both models (30) and ( 32), the augmented process (X t , λ t ) is first-order Markov. As a consequence, for the Poisson-INGARCH(1,1) model, say, we have the following analogy of Proposition 2: Proposition 4. In the Poisson-INGARCH(1,1) model defined by (30), the multi-step-ahead

In other words, in either of the two cases where p = 2 and p = 5, we simulate three models who differ only in terms of the value of λ. This will allow us to investigate how our algorithm performs when the average value of the observations becomes larger. For each of the simulated processes, the sample size is set to T = 500. For illustration purpose, Figure 1 reports the empirical marginal p.m.f.'s of the two INAR(2) processes. As we can see, when λ = 1, the whole process takes small values (at most 6), whereas when λ = 10, we can observe a large proportion of "large" counts. Then we compute the likelihood function using our method to compare its speed with respect to the brute force convolution. For the latter approach we have adopted the recursive algorithm suggested by Bu et al. (2008), Thm 1. Both algorithms are implemented in R using the same laptop (intel i5, 3.0 GHz, 8GB RAM) and the program is available from the author upon request. When order p is small (i.e. p = 2), both methods are roughly equivalent and sufficiently fast.

When p increases to 5, however, our method becomes more competitive. Indeed, its CPU time remains nearly unchanged, whereas the CPU of the direct method is increasing in n. As for the impact of λ (and hence the impact of the expectation of the process) on the CPU time, we can see that for both methods, when λ increases, the computation time increases. Nevertheless the speed of increase is (slightly) faster for our method than for the convolution method. An intuitive explanation of these comparison results is as follows. In Section 2, we have shown that computing the conditional p.m.f. (X t |X t-1 ) involves a complexity of O( 1 2 X 3 t ) using our method, or, roughly speaking, O( 1 2 λ 3 ). On the other hand, Bu et al. (2008)'s method involves the convolution of α 1 • X t-1 , ..., α p • X t-p , and t , by a p-step recursion. Each step involves the convolution of two integer-valued variables for O(X 2 t ) operations. Hence Bu et al. (2008)'s method has a complexity of O(pX 2 t ), or roughly O(pλ 2 ). Thus when λ is fixed and p goes to infinity, our method becomes dominant; whereas when p is fixed but λ goes to infinity, the method of Bu et al. (2008) is more suitable. Note, however, that since count time series models are specifically designed for low counts, the case where λ is very large is somehow scarce and in that case, we might just as well approximate the count time series by a continuously valued time series and use standard continuous time series models.

Multi-step forecasting

We now illustrate how our method fares against the Markov chain approach of McCabe et al.

(2011) and the Monte-Carlo approach [START_REF] Jung | Coherent Forecasting in Integer Time Series Models[END_REF]; [START_REF] Bisaglia | Bayesian Nonparametric Forecasting for INAR Models[END_REF] in terms of both the computational time and precision, when it comes to higher-horizon forecasting.

To this end we use the simulated processes (37). Given observations up to time T , we will compute the conditional distribution of X T +1 , ..., and X T +h given X T .

Clearly, the quality of approximation of the Markov chain approach depends on the choice of the threshold τ . In the following we consider three such values: τ ∈ {6, 10, 20}. As a comparison, by looking at the empirical p.m.f.'s of Figure 1, we can see that the probability that either of the two processes takes values larger than 6 is quite low.

Similarly, by the law of large numbers, the larger the number of replications N in the Monte-Carlo approach, the smaller the difference between the exact predictive distributions and their Monte-Carlo approximations. Thus we will consider three values of N , that are 5000, which is the number recommended by [START_REF] Jung | Coherent Forecasting in Integer Time Series Models[END_REF], as well as N = 50000 and N = 5000000. 2) model. For expository purpose we only report their values at points 0, 1, 2, and 6.

The CPU times correspond to the calculation of p h (k|X t ) for all horizons h = 1, 2, ..., 20 and all arguments k = 0, ..., τ . The Monte-Carlo approach is computationally the most intensive,with a CPU time that is roughly proportional to the number of draws. Our approach, on the other hand, is the quickest, with a computational time that is slowly increasing in τ . As for the Markov chain approach, its CPU time is small when τ = 6, but increases quite significantly in τ . Indeed, the Markov chain approach requires the computation of a matrix of dimension (τ + 1) p × (τ + 1) p , whose computational cost is quadratic in τ for p = 2. Note that in this comparison, we have fixed the values of the parameters. In practice, it might also be preferable to take into account the uncertainty of the parameter forecast. In other words, if the simulation approach is used, one should first draw a large number of possible values of the parameter; then for each parameter value, a large number of draws of future path of the process are drawn. As a consequence, the computational cost would be even much larger than the 6 seconds reported in the table.

Let us now compare the p.m.f.'s obtained from the three approaches. As expected, for different values of τ , our method produces the same conditional p.m.f. at various horizons. For the Monte-Carlo approach, as N increases, the simulated conditional p.m.f.'s become closer and closer to the results obtained by our method, and = 5000000, both approaches provide roughly the same results. Finally, the p.m.f.'s obtained by the Markov chain approach are nearly (resp. exactly, up to a precision of 3 decimal places) the same as the Taylor's expansion. This is expected, since as τ increases, the quality of approximation in the Markov chain approach becomes better and better. Overall, in terms of precision, the Taylor's expansion and the Markov chain approaches strongly dominate the Monte-Carlo approach.

Next, let us illustrate that the advantage of our approach still exists for other thinning-type models satisfying (3). To this end we report a similar experiment for a NBAR(2) process. For illustration we assume that Z 1 , Z 2 and t in equation (3) all follow a special NB distribution: the geometric distribution, which has a degree of freedom parameter r = 1. Their probability parameters q 1 , q 2 are set such that their respective expectations match the corresponding Bernoulli and Poisson parameters in equation ( 2), that is:

q 1 1 -q 1 = q 2 1 -q 2 = 0.2.
This makes the marginal mean of this NBAR(2) process to be same as the above Poisson-INAR(2) process (and equal to

E[ t] 1-E[Z1]-E[Z2] under the stationarity condition E[Z 1 ] + E[Z 2 ] < 1),
hence facilitating the comparison between the two tables. The p.m.f.'s are reported in Table 3 and the conclusions we can draw are similar as for For expository purpose we only report their values at points 0, 1, 2, and 6.

h h = 1 h = 5 h = 10 h = 1 h = 5 h = 10 h = 1 h = 5 h = 10 p h (0|X
In the next exercise, we compare the computational cost of forecasting up to horizon 20 We can see that for the Markov chain approach, the CPU time is often missing when τ and p are large, indicating a failure of the computer program in these cases. This is mainly due to the necessity to compute a matrix of dimension (τ + 1) p × (τ + 1) p . For instance, if p = 5, τ = 10, the algorithm requires a flash memory of 25GB, which exceeds the capacity of a standard PC.

As a consequence, while for small values of p, the Markov chain approach is an acceptable (but more costly and less precise) alternative to our approach, it becomes infeasible for higher-order models.

To conclude this section, Table 5 summarizes the above comparison in three major criteria. We can see that our method combines the strength of its two competitors, while at the same time mitigating their downsides. Finally, note that while for expository purpose we have limited our comparison to models with order p non larger than 5, there is little doubt that for larger p, the advantage of our method becomes more and more

Method

Conclusion

We have solved the open computational difficulty concerning the predictive distributions in thinning-based, higher-order count processes models. This family is quite large and includes, • being specific to the one-step-ahead predictive distribution [see Bu et al. (2008)]

• being specific to INAR(p) models [see [START_REF] Pedeli | Likelihood Estimation for the INAR (p) Model by Saddlepoint Approximation[END_REF]]

Our approach is exact, generically faster than competing methods as demonstrated by simulations, and is unified in the sense that it can be used for forecasting at any horizons, for a large family of models. Our method is based on i) the simple relationship between the p.g.f. and p.m.f. for a count distribution; ii) the CaR property of these processes. Because of its reduced computational cost and wide applicability, our method significantly facilitates the practical implementation of many of these higher-order models.

Finally, is it possible to extend the result of this paper to bivariate count processes? Recently [START_REF] Darolles | Bivariate integer-autoregressive process with an application to mutual fund flows[END_REF] show that for a specific bivariate extension of the INAR(p) models, this is indeed the case. However a systematic treatment of this question requires substantially more investigation and is left for future research.

and the proof of the second equality for B h+1,i (u) is similar and omitted.

Appendix 3: Expanding the log p.g.f. using the p.m.f.

In this section we discuss how the Taylor's expansion of the log p.g.f. can be obtained from the p.m.f. We have:

log E[u X ] = log p 0 + p 1 u + p 2 u 2 + • • • + p τ u τ + O(u τ +1 ) = log p 0 + log 1 + p 1 p 0 u + p 2 p 0 u 2 + • • • + p τ p 0 u τ + O(u τ +1 ) = log p 0 + n k=1 (-1) k k p 1 p 0 u + p 2 p 0 u 2 + • • • + p τ p 0 u τ k + O(u τ +1 ), (40)
for a fixed integer n. Thus by applying Proposition 1, the vector of the first n + 1 coefficients of the Taylor's expansion of log E[u X ] are equal to:

            log p 0 0 • • • 0 0             + τ k=1 (-1) k k             0 0 0 • • • 0 p1 p0 0 0 • • • 0 p2 p0 p1 p0 0 • • • 0 • • • • • • • • • • • • • • • pτ p0 pτ-1 p0 • • • p1 p0 0             k             1 0 0 • • • 0            
.

Remark 1. Clearly, the above series expansion requires p 0 = P[X = 0] to be non zero. Let us discuss the case where this condition is not satisfied by either Z i or t in model (3).

• If either of the Z i 's is such that Z i ≥ 1 almost surely, then we have: X t+1 ≥ X t almost surely for each t, with a nonzero probability that the inequality holds strictly. Thus the process is explosive, nonstationary. This is out of the scope of the present paper.

• If the distribution of t is such that P[ t = 0] = • • • = P[ t = K -1] = 0 and P[ t = K] > 0, then any random variable following this distribution admits the stochastic representation:

X = K + Z,
where Z is again a nonnegative count variable. Thus in this case it suffices to consider the process (Z t ) = (X t ) -K instead of the initial process (X t ). Then the above

(

  2016) introduces a Poisson-INARCH(∞) model (although the author terms this model INAR(∞), in his framework all the variables Z i , as well as t , follow Poisson distributions. Thus his model can also be called Poisson-INARCH(∞)), whereas the Poisson-INGARCH(1,1) model of Ferland et al. (2006) can also be viewed as a special case of Poisson-INARCH(∞) model (see section 3.3 for a more detailed analysis).

Example 6 (

 6 Poisson-INAR(p) model). Let us now consider the benchmark Poisson-INAR(p) model (

Figure 1 :

 1 Figure 1: Empirical marginal p.m.f. of the simulated INAR(2) process with λ = 1 (left panel) and λ = 10 process (right panel).

  using different approaches, when we increase progressively the lag p of the Poisson-INAR(p) model, from p = 2 to 5. The individual models are obtained recursively and backwardly from the INAR(5), by dropping each time the last term in equation (3), while keeping the other parameters unchanged. We have not reported the comparison of the precision of the p.m.f.'s, as these results are rather similar across different values of p. Method Taylor's expansion Markov chain Monte-Carlo, N = 50000 τ = 6 τ = 10 τ = 20 τ = 6 τ = 10 τ = 20 N = 5000 N = 50000 N = 500000 p = 2 0.

  among others, Poisson-INARCH(p), Poisson-INAR(p), NB-INAR(p), NBAR(p), as well as some infinite order models such as Poisson-INGARCH(1,1). Existing methods usually suffer from one or several of the following downsides, such as: • introducing significant approximation errors [see e.g. Pedeli et al. (2015)] • being time consuming [see e.g. McCabe et al. (2011); Jung and Tremayne (2006)]

  1 ] are often straightforward to compute, they are not informative enough to characterize the predictive count distributions. Moreover,[START_REF] Mccabe | Efficient Probabilistic Forecasts for Counts[END_REF] argue that since the predictive means are non-integer, they are incompatible with the discrete sample space. A proper, probabilistic forecast thus necessitates the multiple-step-ahead conditional p.m.f. of X t+h given

	X t-1 :

Table 1

 1 compares their respective CPU time. Note that we do not compare the likelihood values, since both methods provide the same, exact likelihood.

	Value of p		p = 2			p = 5
	Value of λ	λ = 1	λ = 5	λ = 10	λ = 1	λ = 5 λ = 10
	Simulated empirical mean	1.59	7.9	16.5	2.28	12.3	19.6
	Simulated empirical variance	1.72	7.8	17.9	2.44	14.0	25.6
	CPU of our method	0.006 s 0.015 s 0.033 s 0.007 s 0.02 s 0.05 s
	CPU of convolution method 0.008 s 0.016 s 0.031 s 0.02 s	0.06s	0.11 s

Table 1 :

 1 CPU time of the two methods for different values of p and λ. The time unit is second.

Table 2 :

 2 Table 2 reports the CPU time and the conditional p.m.f.'s at horizons h = 1, 5 and 10 obtained from the three methods, given observations up to time T of the INAR(2) process, with λ = 1. Conditional p.m.f. at different horizons obtained by different methods for the Poisson-INAR(

	Method	Taylor's expansion, τ = 6	Markov chain, τ = 6	Monte-Carlo, N = 5000
	CPU time	0.003 s			0.008 s		0.6 s
	horizon h	h = 1 h = 5	h = 10	h = 1 h = 5 h = 10 h = 1 h = 5	h = 10
	p h (0|X T )	0.235 0.198	0.193	0.236 0.198	0.193	0.231 0.210	0.200
	p h (1|X T )	0.353 0.317	0.315	0.353 0.317	0.314	0.352 0.312	0.314
	p h (2|X T )	0.250 0.258	0.259	0.25	0.258	0.259	0.258 0.255	0.248
	p h (6|X T )	0.002 0.006	0.006	0.002 0.006	0.006	0.002 0.005	0.007
	Method	Taylor's expansion, τ = 10	Markov chain, τ = 10	Monte-Carlo, N = 50000
	CPU time	0.005 s			0.06 s		6 s
	horizon h	h = 1 h = 5	h = 10	h = 1 h = 5 h = 10 h = 1 h = 5	h = 10
	p h (0|X T )	0.235 0.198	0.193	0.235 0.198	0.193	0.232 0.199	0.191
	p h (1|X T )	0.353 0.317	0.314	0.353 0.317	0.314	0.356 0.321	0.314
	p h (2|X T )	0.250 0.258	0.259	0.250 0.261	0.271	0.249 0.259	0.264
	p h (6|X T )	0.002 0.006	0.006	0.002 0.006	0.006	0.002 0.006	0.006
	Method	Taylor's expansion, τ = 20	Markov chain, τ = 20	Monte-Carlo, N = 500000
	CPU time	0.005 s			1.2 s		60 s
	horizon h	h = 1 h = 5	h = 10	h = 1 h = 5 h = 10 h = 1 h = 5	h = 10
	p h (0|X T )	0.235 0.198	0.193	0.235 0.198	0.193	0.232 0.199	0.191
	p h (1|X T )	0.353 0.317	0.314	0.353 0.317	0.314	0.356 0.321	0.314
	p h (2|X T )	0.250 0.258	0.259	0.250 0.258	0.259	0.253 0.257	0.257
	p h (6|X T )	0.002 0.006	0.006	0.002 0.006	0.006	0.002 0.006	0.006

  Table 2 and are thus omitted.

	Method	Taylor's expansion, τ = 6	Markov chain, τ = 6	Monte-Carlo, N = 5000
	CPU time	0.003 s		0.008 s		0.6 s
	horizon h	h = = 5	h = 10	h = 1 h = 5 h = 10 h = 1 h = 5	h = 10
	p h (0|X T )	0.379 0.344	0.339	0.379 0.344	0.339	0.375 0.357	0.332
	p h (1|X T )	0.265 0.250	0.250	0.265 0.250	0.250	0.270 0.247	0.260
	p h (2|X T )	0.159 0.162	0.163	0.159 0.162	0.163	0.163 0.162	0.158
	p h (6|X T )	0.014 0.020	0.021	0.014 0.020	0.021	0.015 0.020	0.019
	Method	Taylor's expansion, τ = 10	Markov chain, τ = 10	Monte-Carlo, N = 50000
	CPU time	0.005 s		0.06 s		0.6 s
	horizon h	h = 1 h = 5	h = 10	h = 1 h = 5 h = 10 h = 1 h = 5	h = 10
	p h (0|X T )	0.379 0.344	0.339	0.379 0.344	0.339	0.378 0.343	0.340
	p h (1|X T )	0.265 0.250	0.250	0.265 0.250	0.250	0.263 0.252	0.246
	p h (2|X T )	0.159 0.162	0.163	0.159 0.162	0.163	0.160 0.161	0.167
	p h (6|X T )	0.014 0.020	0.021	0.014 0.020	0.021	0.014 0.021	0.020
	Method	Taylor's expansion, τ = 20	Markov chain, τ = 20	Monte-Carlo, N = 500000
	CPU time	0.005 s		1.2 s		60 s
	horizon

Table 3 :

 3 Conditional p.m.f. at different horizons obtained by different methods for the NBAR(2) model.

	T )	0.379 0.344	0.339	0.379 0.344	0.339	0.38	0.344	0.340
	p h (1|X T )	0.265 0.250	0.250	0.265 0.250	0.250	0.265 0.250	0.250
	p h (2|X T )	0.159 0.162	0.163	0.159 0.162	0.163	0.159 0.163	0.163
	p h (6|X T )	0.014 0.020	0.021	0.014 0.020	0.021	0.015 0.020	0.021

Table 4 :

 4 CPU time for forecasting in a Poisson-INAR(5) model, using different approaches.

		003	0.005	0.006	0.008	0.06	1.2	0.6	6	60
	p = 3	0.004	0.006	0.006	0.20	0.4	5.2	0.8	9	58
	p = 4	0.004	0.006	0.006	6.0	50	NA	0.11	10	90
	p = 5	0.004	0.006	0.007	180	NA	NA	0.13	12	116

Table 5 :

 5 Comparison of the three approaches in terms of three major criteria. The sensitivity concerns the computational cost for forecasting.

		Taylor's expansion Markov chain Monte-Carlo
	Computation accuracy	Very High	High	Low
	Sensitivity of CPU time w.r.t. threshold τ	Low	Very High	NA
	Sensitivity of CPU time w.r.t. lag p	Low	Very High	High

I thank the two referees and the Associate Editor for their insightful comments and suggestions. I thank the ACPR chair in Regulation and Systemic Risk, the CNRS (délégation grant), the Society of Actuaries (2018 individual grant) as well as the ANR (Labex MME-DII grant) for financial support.

conditional p.g.f. is:

with B 1,1 (u) = u -1, B 1,0 (u) = 0, (34) B h+1,1 (u) = e aB h,1 (u) -1 + ρB h,1 (u), (35)

Proof. See Appendix 4.

Using this lemma, we can mimic the procedure described in Section 3.1 to obtain recursively the Taylor's expansion of B h,1 (u) and B h,0 (u) at zero, and then expand the exponential function to recover the exact value of the h-step-ahead conditional p.m.f. Since the conditional p.g.f.

given X t-1 is completely characterized by λ t , the resulting computational cost is comparable to model (3) with p = 1.

Finally, Proposition 4 can also be adapted to the NB-INGARCH(1,1) model ( 32). Due to space constraint this analysis is relegated to Appendix 5.

4 Comparison with existing methods using simulated data

Likelihood function

To illustrate the computational gain, we first implement the conditional likelihood function, that is the likelihood function given the first p initial values of the process, using both our method and the direct application of (4). We consider a Poisson-INAR(2) model, whose parameters and initial values are fixed as:

as well as a Poisson-INAR(5) model with parameters and initial values set as:

Appendix 1: Proof of Proposition 2

Let us proceed by induction. case where h = 1 is a direct consequence of Equation ( 5).

Assume that (20) holds for a certain horizon h, then:

Thus (20) holds for any h.

Appendix 2: Proof of Proposition 3

We have:

algorithm continues to be applicable.

Appendix 4: Proof of Proposition 4

We proceed by induction. The initial condition h = 1 is obvious. Let us assume that equations ( 35) and ( 36) are satisfied for a certain h ≥ 1. Then we have:

= e ρB h,1 (u)λt+B h,1 (u)λ0+B h,0 (u)+λt(e aB h,1 (u) -1) .

Hence ( 35) and ( 36) are also satisfied for h + 1.

Appendix 5: The predictive distribution in NB-INGARCH(1,1) model

In this section we establish the analogue of Proposition 4 for NB-INGARCH(1,1) model defined by (32). We have:

Proposition In the NB-INGARCH(1,1) model defined by (32), the multi-step-ahead conditional p.g.f. is:

with B 1,1 (u) = log( 1 -q 1 -qu ), B 1,0 (u) = 0, B h+1,1 (u) = log( 1 -q 1 -qe aB h,1 (u) ) + ρB h,1 (u), B h+1,0 (u) = B h,1 (u)λ 0 + B h,0 (u), ∀h ≥ 0.

The proof (by induction) is similar to the proof of Proposition 4 and is thus omitted.