Serge Darolles

Gaëlle Le Fol

Yang Lu

Ran Sun

Bivariate INAR Processes with an Application to Mutual Fund Flows

Keywords: Multivariate low event count process, non-linear forecasting, memory persistence, liquidity risk, mutual funds. MSC code: 62-15, JEL code: C32, C53

We propose a new family of bivariate nonnegative integer-autoregressive (BI-NAR) models for count process data. We first generalize the existing BINAR(1) model by allowing for dependent thinning operators and arbitrary innovation distribution. The extended family allows for intuitive interpretation, as well as tractable aggregation and stationarity properties. We then introduce higher order BINAR(p) and BINAR(∞) dynamics to accommodate more flexible serial dependence patterns. So far, the literature has regarded such models to be computationally intractable. Our contribution is to show that the extended BINAR family allows for closed-form predictive distributions at any horizons and for any values of p, which significantly facilitates non-linear forecasting and likelihood based estimation. Finally, a BINAR(∞) model with memory persistence is applied to open-ended mutual fund purchase and redemption order counts.

Introduction

Nonnegative low count processes has been 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]], finance [START_REF] Heinen | Multivariate Autoregressive Modeling of Time Series Count Data using Copulas[END_REF]], insurance [START_REF] Gouriéroux | Heterogeneous INAR (1) Model with Application to Car Insurance[END_REF]] and beyond, ever since the seminal work of [START_REF] Mckenzie | Some Simple Models for Discrete Variate Time Series[END_REF]. Our interest in this paper lies in the monitoring of the liquidity risk of an open-ended mutual fund (MF). A MF channels investors' cash investment into less liquid assets, and is thus structurally vulnerable to liquidity risk. This risk has recently received much attention of the regulators [see, e.g., Securities and Exchange Commission (2015); [START_REF] Darolles | Liquidity Risk and Investor Behavior: Issues, Data and Models[END_REF]], but its quantification and management remains difficult. Indeed, from the modelling point of view, the liquidity risk is quite different from traditional market risks in that they involve the daily counts of redemption and purchase orders, which are i) most of the time low integers or zero, but also have a non-null probability of taking mildly large values; ii) both cross-sectionally and serially dependent, with significant heteroscedasticity. Recently, the MF industry has started to record purchase and redemption order count data separately. This allows to distinguish auto-correlation effects and cross-effects between the two count processes, which have different economic interpretations. For instance, the clustering of the redemption counts corresponds to fund run, whereas a fund manager usually reacts to past redemptions by seeking new investors in order to stabilize the fund size, leading to a positive feedback effect between past redemption and current purchase counts.

Therefore, a bivariate count analysis can be of great interest to understand clients' behaviour and the manager's reaction to exogenous liquidity shock.

Yet the literature on bivariate count processes is still in its infancy. The benchmark approach is the Bivariate INteger-valued AutoRegressive (BINAR) model [START_REF] Latour | The Multivariate GINAR (p) Process[END_REF]; [START_REF] Pedeli | Some Properties of Multivariate INAR (1) Processes[END_REF]], which assumes that for each t:

X 1,t = α 11 • X 1,t-1 + α 12 • X 2,t-1 + 1,t , X 2,t = α 21 • X 1,t-1 + α 22 • X 2,t-1 + 2,t , (1.1)
where given X t-1 = (X 1,t-1 , X 2,t-1 ) , the binomial thinning operators are defined as follows:

for each i, j = 1, 2, variable α i,j • X j,t-1 follows the binomial distribution with size X j,t-1 and success probability α i,j ∈ [0, 1]. Moreover these variables are conditionally independent, and are unbounded count data. In this paper we show that one of these extensions has indeed tractable properties, even after extension to higher-order cases.

More precisely, we contribute to the BINAR literature in two ways. First, we extend model (1.1) (called independent BINAR(1) henceforth) by introducing (positively or negatively) dependent thinning operators and arbitrary innovation distribution. We also show that the process belongs to the compound autoregressive (CaR) family, and possesses intuitive aggregation and stationarity properties. More importantly, we clarify that in this large family of models, the predictive distributions at various horizons are easily computable via a matrix-based algorithm.

This largely facilitates likelihood-based inference and non-linear forecasting, especially when it comes to the prediction of extreme events. Second, we extend our model to higher-order dependent BINAR(p) and BINAR(∞) processes, which can better capture slowly decaying serial correlation patterns.

The paper is organized as follows. The dependent BINAR(1) model is introduced in Section 2 and extended in Section 3 to higher-order BINAR(p) and BINAR(∞) models. Section 4 computes the predictive distributions. Section 5 applies the model to forecast the counts of share purchase and redemption of a MF. Section 6 concludes. Proofs and technical details are gathered in Appendices.

2 Dependent BINAR(1) process 2.1 The dynamic specification Definition 1. We say that the bivariate count process (X t ) = (X 1t , X 2t ) with domain N 2 ≥0 is dependent BINAR(1) if it has the stochastic representation:

   X 1,t X 2,t    = X1,t-1 i=1    Z 1,i,t Z 2,i,t    + X2,t-1 j=1    Z 3,j,t Z 4,j,t    +    1,t 2,t    , ∀t, (2.1)
where given X t-1 :

• for i, j, t varying, random vectors (Z 1,i,t , Z 2,i,t ) , (Z 3,j,t , Z 4,j,t ) , ( 1,t , 2,t ) are mutually independent copies of (Z 1 , Z 2 ) , (Z 3 , Z 4 ) , and ( 1 , 2 ) , respectively.

• Couples (Z 1 , Z 2 ) and (Z 3 , Z 4 ) are independent and bivariate Bernoulli distributed. That is, they have Bernoulli marginal distribution with success probability parameters given by:

   P[Z 1 = 1], P[Z 3 = 1] P[Z 2 = 1], P[Z 4 = 1]    =    α 11 , α 12 α 21 , α 22    := A,
where all the entries of matrix A are nonnegative, whereas the two joint probabilities are respectively:

P[Z 1 = Z 2 = 1] = q 1 , P[Z 3 = Z 4 = 1] = q 2 .
• Innovation t is independent of X t-1 , X t-2 , ..., and is i.i.d. across t. Moreover, it is nonnegative and has finite variance.

The bivariate Bernoulli distribution is first introduced by [START_REF] Wicksell | Some theorems in the theory of probability, with special reference to their importance in the theory of homograde correlations[END_REF] [see also [START_REF] Marshall | A Family of Bivariate Distributions Generated by the Bivariate Bernoulli Distribution[END_REF]], and is recently used by [START_REF] Scotto | Bivariate Binomial Autoregressive Models[END_REF] to model bounded count processes.In order for it to be well defined, its parameters have to satisfy the following constraint:

Lemma 1 [START_REF] Joe | Multivariate Models and Multivariate Dependence Concepts[END_REF], page 210). Parameters (α 11 , α 12 , α 21 , α 22 , q 1 , q 2 ) of the bivariate Bernoulli distribution can take any values satisfying the following two inequalities: max(α 11 + α 21 -1, 0) ≤ q 1 ≤ min(α 11 , α 21 ), (2.2)

max(α 12 + α 22 -1, 0) ≤ q 2 ≤ min(α 12 , α 22 ). (2.
3)

The covariance between Z 1 and Z 2 (resp. Z 3 and Z 4 ) is q 1 -α 11 α 21 . In particular, if q 1 = max(α 11 + α 21 -1, 0), then the covariance is nonpositive; if q 1 = min(α 11 , α 21 ), the covariance is nonnegative; if q 1 = α 11 α 21 and q 2 = α 12 α 22 , we get the independent BINAR(1) model (1.1).

Birth-death-immigration interpretation. Let us temporarily assume that X 1,t and X 2,t count individuals of type 1 and 2 at time t, respectively. A type j 1 (j 1 = 1, 2) individual produces a type j 2 (j 2 = 1, 2) off-spring with marginal probability α j1,j2 , and joint probability q j1 . Then the population of type j (j = 1, 2) at the next period t + 1 is composed of such off-springs, plus j,t+1 immigrants of type j. In particular, if q 1 = α 11 α 21 and q 2 = α 12 α 22 , the productions of the two types of off-springs are independent events. If instead q 1 = q 2 = 0, the productions are mutually exclusive, i.e., each individual can only produce up to one off-spring.

Conditional distribution

First two conditional moments

Since we have:

E[X t |X t-1 ] = AX t-1 + E[ t ],
(2.4) we abbreviate model (2.1) into:

X t = A(q 1 , q 2 ) • X t-1 + t , (2.5)
where the operator A(q 1 , q 2 )• is linear and will be called the dependent matrix thinning operator.

The conditional covariance matrix is:

   V[X 1,t |X t-1 ], * Cov[X 1,t , X 2,t |X t-1 ], V[X 2,t |X t-1 ]    = γ + X 1,t-1 γ 12 + X 2,t-1 γ 34 , (2.6)
where γ , γ 12 , γ 34 are the covariance matrices of ( 1t , 2t ), (Z 1,t , Z 2,t ) and (Z 3,t , Z 4,t ), respectively:

γ 12 =    α 11 (1 -α 11 ), * q 1 -α 11 α 21 , α 21 (1 -α 21 )    , γ 34 =    α 12 (1 -α 12 ), * q 2 -α 12 α 22 , α 22 (1 -α 22 )    , (2.7)
and the conditional correlation between X 1,t and X 2,t is:

corr[X 1,t , X 2,t |X t-1 ] = Cov[ 1,t , 2,t ] + X 1,t-1 Cov[Z 1 , Z 2 ] + X 2,t-1 Cov[Z 3 , Z 4 ] V[ 1t ] + X 1,t-1 V[Z 1 ] + X 2,t-1 V[Z 3 ] V[ 2t ] + X 1,t-1 V[Z 2 ] + X 2,t-1 V[Z 4 ] .
In the above equation, the denominator of the right hand side does not depend on q 1 or q 2 , whereas the numerator is increasing in q 1 and q 2 . In the independent BINAR(1) model with q 1 = α 11 α 21 , q 2 = α 22 α 12 , the above correlation coefficient is (in absolute value) non larger than corr( 1t , 2t ), and becomes small whenever X 1,t-1 and/or X 2,t-1 are large, thus the conditional heteroscedasticity cannot be well captured. This downside exists also in several other competing bivariate count process models. For instance, in the bivariate Poisson autoregression of [START_REF] Liu | Some Models for Time Series of Counts[END_REF], the conditional correlation coefficient goes to zero when X 1t and X 2t are large. In copula-based bivariate count processes [START_REF] Heinen | Multivariate Autoregressive Modeling of Time Series Count Data using Copulas[END_REF]], this conditional correlation coefficient can only be computed numerically and it is not clear how it behaves when components of X t-1 are large; in our model, when both X 1,t-1 and X 2,t-1 are large, the conditional correlation is approximately:

corr[X 1,t , X 2,t |X t-1 ] ≈ X 1,t-1 Cov[Z 1 , Z 2 ] + X 2,t-1 Cov[Z 3 , Z 4 ] X 1,t-1 V[Z 1 ] + X 2,t-1 V[Z 3 ] X 1,t-1 V[Z 2 ] + X 2,t-1 V[Z 4 ] ,
which can be closer to 1 (resp. -1) if q 1 and q 2 are close to their upper (resp. lower) bounds.

Finally, if 1t , 2t have Poisson distributions, then both components are conditionally underdispersed:

V[X jt |X t-1 ] ≤ E[X jt |X t-1 ], j = 1, 2.
This differs from the models of [START_REF] Cui | A new bivariate integer-valued garch model allowing for negative cross-correlation[END_REF]; [START_REF] Gouriéroux | Negative Binomial Autoregressive Process with Stochastic Intensity[END_REF], which are conditionally equi-dispersed and over-dispersed, respectively. In general, by leaving the distribution of t unconstrained, the BINAR(1) family allows for flexible conditional dispersion.

Conditional p.g.f.

The dynamics of process (X t ) is characterized by its conditional probability generating function (p.g.f.), which is equal to:

E[u X1,t v X2,t |X t-1 ] = E[u Z1 v Z2 ] X1,t-1 E[u Z3 v Z4 ] X2,t-1 E[u 1,t v 2,t ] = a 1 (u, v) X1,t-1 a 2 (u, v) X2,t-1 b(u, v), ∀u, v ≥ 0, (2.8)
where b(u, v) is the p.g.f. of ( 1,t , 2,t ) and

a 1 (u, v) = q 1 uv + (α 11 -q 1 )u + (α 21 -q 1 )v + (1 + q 1 -α 11 -α 21 ), ∀u, v ≥ 0, (2.9) a 2 (u, v) = q 2 uv + (α 12 -q 2 )u + (α 22 -q 2 )v + (1 + q 2 -α 12 -α 22 ), ∀u, v ≥ 0 (2.10)
are the p.g.f. of (Z 1 , Z 2 ) [resp. (Z 3 , Z 4 )] respectively. This conditional p.g.f. is exponential affine in X t-1 . Such a process is called compound autoregressive (CaR) [see [START_REF] Darolles | Structural Laplace Transform and Compound Autoregressive Models[END_REF]]. A remarkable property of such processes is that the multiple-step-ahead conditional p.g.f. remains exponential affine in X t-1 :

Proposition 1. In model (2.1), we have, for any horizon h ≥ 0:

E[u X 1,t+h-1 v X 2,t+h-1 |X t-1 ] = a (h) 1 (u, v) X1,t-1 a (h) 2 (u, v) X2,t-1 b (h) (u, v), ∀u, v ≥ 0, (2.11)
where functions a v) are defined by the recursion:

(h) 1 (u, v), a (h) 2 (u, v) and b (h) (u,
a (h+1) 1 (u, v) = 1 + α 11 [a (h) 1 (u, v) -1] + α 21 [a (h) 2 (u, v) -1] + q 1 [a (h) 1 (u, v) -1][a (h) 2 (u, v) -1], a (h+1) 2 (u, v) = 1 + α 12 [a (h) 1 (u, v) -1] + α 22 [a (h) 2 (u, v) -1] + q 2 [a (h) 1 (u, v) -1][a (h) 2 (u, v) -1], b (h+1) (u, v) = b (h) (u, v)b(a (h) 1 (u, v), a (h) 2 (u, v)), ∀u, v ≥ 0, h ∈ N, with initial conditions a (0) 1 (u, v) = u, a (0) 2 (u, v) = v, and b (0) (u, v) = 1.
The proof (by induction) is straightforward and thus omitted. This proposition implies that

for each h ≥ 1, a (h) 1 (u, v) and a (h)
2 (u, v) are polynomials of degree 2 h-1 in both arguments, except when q 1 = q 2 = 0. More precisely:

Corollary 1. If q 1 = q 2 = 0, then:    a (h) 1 (u, v) -1 a (h) 2 (u, v) -1    = A    a (h-1) 1 (u, v) -1 a (h-1) 2 (u, v) -1    = (A ) h    u -1 v -1    .
(2.12)

Thus in this case a

(h) 1 (u, v) and a (h)
2 (u, v) are affine instead of polynomial. In other words, the conditional p.g.f. of [A(0, 0)•] (h) X t-1 given X t-1 at any horizon h has the same functional form as that of A h (0, 0) • X t-1 , or equivalently, the operator [A(0, 0)•] (h) = A h (0, 0)• is still a dependent matrix thinning operator. As a consequence, in terms of temporal aggregation, when observed at a lower frequency of h, the process (X th ) t is still BINAR(1).

Remark 1. Corollary 1 is easily explained by the birth-immigration interpretation: if q 1 = q 2 = 0, each individual produces at most one off-spring in the next period, and thus at most one offspring at horizon h ≥ 2. Hence the identity [A(0, 0)•] (h) = A h (0, 0) Remark 2. Note that the literature [see, e.g., [START_REF] Boudreault | Multivariate Integer-Valued Autoregressive Models Applied to Earthquake Counts[END_REF], Thm. 2.8, as well as Pedeli and Karlis (2013b), eq. 15] usually claims that in the independent BINAR model where q 1 = α 11 α 21 , q 2 = α 21 α 22 , the composite operator [A(q 1 , q 2 )•] (h) is equal to the bivariate thinning operator A h (q 1,h , q 2,h )•, where matrix A h = A h , and q 1,h (resp. q 2,h ) is the products of the two entries of the first (resp. second) column of A h . Then these authors deduce that the h-step-ahead conditional p.g.f. still has the CaR form of equation ( 1), but with functions a (h) 1

and a (h) 2

being affine rather than higher-order polynomial. From the above analysis we can see that this assertion is incorrect.

Notice also that equation (2.11) corresponds to the decomposition:

X t+h-1 = [A(q 1 , q 2 )•] (h) X t-1 + [A(q 1 , q 2 )•] (h-1) t + • • • + A(q 1 , q 2 ) • t+h-2 + t+h-1 ,
where the additive terms on the right hand side are conditionally independent given X t-1 , with conditional p.g.f.'s a ..., and b(u, v) for t+h-1 , respectively.

(h) 1 (u, v) X1,t-1 a (h) 2 (u, v) X2,t-1 for [A(q 1 , q 2 )•] (h) X t-1 , b(a (h-1) 1 (u, v), a (h-1) 2 (u, v)) for [A(q 1 , q 2 )•] (h-1) t ,
Example 1. Let us consider the case where 1t , 2t are mutually independent, Poisson P(λ 1 ), P(λ 2 ) distributed, respectively. We have b(u, v) = exp λ 1 (u -1) + λ 2 (v -1) , and the p.g.f. of

A(q 1 , q 2 ) • t+h-2 is: b(a 1 (u, v), a 2 (u, v)) = exp λ 1 (a 1 (u, v) -1) + λ 2 (a 2 (u, v) -1)
= exp λ 1 (α 11 -q 1 )(u -1) + (α 21 -q 1 )(v -1) + q 1 (uv -1) (2.13)

+ λ 2 (α 12 -q 2 )(u -1) + (α 22 -q 2 )(v -1) + q 2 (uv -1) = exp m 1 (u -1) + m 2 (v -1) + m 3 (uv -1) , (2.14)
where

m 1 = λ 1 (α 11 -q 1 )+λ 2 (α 12 -q 2 ), m 2 = λ 2 (α 21 -q 1 )+λ 2 (α 22 -q 2 ), m 3 = (λ 1 q 1 +λ 2 q 2 ).
This is the p.g.f. of a bivariate (dependent) Poisson distribution BP (m 1 , m 2 , m 3 ) with trivariate reduction [see, e.g., [START_REF] Marshall | A Family of Bivariate Distributions Generated by the Bivariate Bernoulli Distribution[END_REF]].Its correlation coefficient is:

ρ 1 = m 3 (m 1 + m 3 )(m 2 + m 3 ) = λ 1 q 1 + λ 2 q 2 (λ 1 α 11 + λ 2 α 21 )(λ 1 α 12 + λ 2 α 22 ) ,
which is nonnegative, and increasing in q 1 and q 2 . As a consequence, the conditional distribution

X t+1 |X t-1 is BP (m 1 + λ 1 , m 2 + λ 2 , m 3 ).
Under the same assumption, the p.g.f. of

[A(q 1 , q 2 )•] (2) t+h-3 is: b(a (2) 1 (u, v), a (2) 2 (u, v)) = exp λ 1 [a (2) 1 (u, v) -1] + λ 2 [a (2) 2 (u, v) -1] , (2.15)
which is exponential quadratic in u and v. The family of distributions with exponential quadratic p.g.f. is called bivariate Hermite (BH). It nests the bivariate Poisson as special case and is closed under convolution [see [START_REF] Kemp | Bivariate Hermite Distributions[END_REF]]. In particular, the conditional distribution X t+1 |X t-1 still belongs to the BH family.

Cross-sectional aggregation

In the previous subsection we have analyzed the frequency aggregation property of the BINAR(1) process. [START_REF] Brännäs | A New Approach to Modelling and Forecasting Monthly Guest Nights in Hotels[END_REF] and α 12 + α 22 , respectively, if and only if they only take values 0 and 1, that is when q 1 = q 2 = 0.

Nevertheless, given X 1,t-1 + X 2,t-1 , the sum

X1,t-1 j=1 (Z 1,j,t + Z 2,j,t ) + X2,t-1 j=1 (Z 3,j,t + Z 4,j,t ) is
generically not Binomial, except when α 11 + α 21 = α 12 + α 22 , in which case (2.16) becomes: (2.17) where Z 1,j,t + Z 2,j,t , j varying are independent of X t-1 , X t-2 , ..., and the innovation 1,t + 2,t is independent of the variables Z 1,j,t + Z 2,j,t . To summarize, we have the following property:

X 1,t + X 2,t = X1,t-1+X2,t-1 j=1 (Z 1,j,t + Z 2,j,t ) + 1,t + 2,t ,
Proposition 2. In the BINAR(1) model with q 1 = q 2 = 0 and

α 11 + α 21 = α 12 + α 22 , the sum process (X 1,t + X 2,t ) is also univariate INAR(1) with autocorrelation coefficient α 11 + α 21 .
Thus the sum process is Markov with respect to its own history. Moreover, eq. ( 2.17) says that:

(X 1,t + X 2,t |X 1,t-1 , X 2,t-1 ) depends on X 1,t-1 , X 2,t-1 only via the sum X 1,t-1 + X 2,t-1 .
In other words, the sum process can be viewed as an exogenous, common Markov factor.

It is also possible to interpret condition α 11 + α 21 = α 12 + α 22 in terms of the second order dynamics of the process. To this end we exclude two degenerate cases where matrix A is diagonal or anti-diagonal. Then condition α 11 + α 21 = α 12 + α 22 implies that

(1, 1) A = (α 11 + α 21 )(1, 1) ,
that is, α 11 + α 21 is an eigenvalue of matrix A associated with eigenvector (1, 1) . Since matrix A and vector (1, 1) are positive, by the Perron-Frobenius theorem, α 11 + α 21 is, in modulus, the simple largest eigenvalue of A. Thus, among all the linear combinations of components of X t , process (X 1,t + X 2,t ) has the largest autocorrelation coefficient α 11 + α 21 . This justifies its interpretation as a common factor.

Stationarity

The stationarity condition

The strict stationarity condition of the BINAR(1) process is given in the next proposition:

Proposition 3. Process (X t ) is both strictly and mean-variance stationary if and only if:

(1 -α 11 )(1 -α 22 ) > α 12 α 21 , (2.18)
or equivalently, if and only if the eigenvalues of A are smaller than 1 in modulus.

Note also that under condition (2.18), inequalities α 21 + α 11 > 1 and α 22 + α 12 > 1 cannot hold simultaneously. Thus in inequalities (2.2) and (2.3), at least one of the lower bounds are effectively zero.

The stationary distribution

Let us denote by b ∞ the p.g.f. of the stationary distribution. By taking expectation in (2.8), we get:

b ∞ (u, v) = b ∞ (a 1 (u, v), a 2 (u, v))b(u, v) ⇐⇒ b ∞ (u, v) = ∞ i=1 b(a (h) 1 (u, v), a (h) 2 (u, v)).
This latter expression is generically complicated, but can be simplified in the special case considered in Corollary 1:

Proposition 4. If:

• Processes ( 1,t ) and ( 2,t ) are mutually independent and Poisson P(λ 1 ), P(λ 2 ) distributed;

• and q 1 = q 2 = 0 (which is possible only if α 11 + α 21 < 1 and α 12 + α 22 < 1), then:

• the stationary distribution of process (X t ) is such that X 1t and X 2t are independent, Poisson P(λ 1,∞ ), P(λ 2,∞ ) distributed, with parameters given by:

(λ 1,∞ , λ 2,∞ ) = (Id -A) -1 (λ 1 , λ 2 ) ,
• the sum process

(Y t ) = (X 1,t + X 2,t ) is Poisson INAR(1)
, with representation:

Y t = Yt-1 j=1 B jt + η t , (2.19)
where given Y t-1 , variables (B jt ) j are i.i.d. Bernoulli with probability α = 1 -λ1+λ2 λ1,∞+λ2,∞ , and andY t-1 , andPoisson P(λ 1 + λ 2

(η t ) = ( 1,t + 2,t ) is i.i.d., independent of B jt
) distributed. Proof. See Appendix A.1.
Proposition 3 is the bivariate analogue of the well known result that the stationary distribution of a univariate Poisson-INAR(1) process is Poisson [see [START_REF] Mckenzie | Some Simple Models for Discrete Variate Time Series[END_REF]]. It is interesting to compare it with Proposition 2, since in both cases we have assumed q 1 = q 2 = 0 and the sum process is INAR(1). On the one hand, Proposition 2 requires condition α 11 + α 21 = α 12 + α 22 , but leaves the distribution of the innovation ( t ) unconstrained. On the other hand, Proposition 4 does not restrict matrix A, but is based on the independent Poisson assumption of ( t ).

The marginal moments

The simple conditional expectation allows us to derive the first two marginal moments of the process. We have:

Proposition 5. The marginal expectation of the process is given by:

E[X t ] = (Id -A) -1 E[ t ],
(2.20)

the covariance matrix is:

Γ(0) :=    V[X 1,t ], * Cov[X 1,t , X 2,t ], V[X 2,t ]    = ∞ h=0 A h γ + E[X 1t ]γ 12 + E[X 2t ]γ 34 (A ) h , (2.21)
whereas the autocovariance matrices are:

Γ(h) :=    Cov[X 1,t , X 1,t-h ], Cov[X 1,t , X 2,t-h ] Cov[X 2,t , X 1,t-h ], Cov[X 2,t , X 2,t-h ]    = A h Γ(0), ∀h ∈ N. (2.22) Proof. See Appendix A.2. Thus, since matrix A is nonnegative, the variances V[X 1t ], V[X 2t ] are increasing in q 1 , q 2 ,
ceteris paribus, whereas the marginal expectations do not depend on these two probabilities. Thus the marginal over-dispersion coefficients

V[X1,t] E[X1,t] , V[X2,t] E[X2,t]
has a wider range than the independent BINAR model.

To illustrate this larger flexibility, we compare three models which differ only by q 1 , q 2 . We

set A =    0.5, 0.3 0.4, 0.5  
 and assume 1,t , 2t to be independent, P(1) distributed. Finally, q 1 , q 2 are specified as follows:

• In Model 1, q 1 = min(α 1,1 , α 1,2 ) = 0.4, q 2 = min(α 2,1 , α 2,2 ) = 0.3, i.e. both bivariate Bernoulli variables (Z 1 , Z 2 ) and (Z 3 , Z 4 ) have maximal correlation. • In the [independent BINAR(1)] Model 2, q 1 = α 1,1 α 1,2 = 0.2, q 2 = α 2,1 α 2,2 = 0.15, i.e.,
both bivariate Bernoulli variables have zero correlation.

• In Model 3, q 1 = q 2 = 0, i.e. both bivariate Bernoulli variables have minimal correlation.

Table 1 reports the over-dispersion coefficients and correlation coefficients computed using Proposition 4, under the three above models. Note that in Model 3, The two over-dispersion coefficients as well as the correlation coefficient are equal to 1, 1, 0, respectively. Indeed in this case the marginal stationary distribution is bivariate independent Poisson. As expected, these coefficients are largest (resp. smallest) in Model 1 (resp. Model 2). In other words, by letting q 1 , q 2 vary, the dependent BINAR model allows for a much larger range of over-dispersion and coerrelation coefficients than the benchmark Model 2.

V[X 1,t ]/E[X 1,t ] V[X 1,t ]/E[X 1,t ] corr[X 1,t , X 2,t ]
3 Higher-order BINAR processes

BINAR(p) process

Similar as the INAR(p) process introduced by [START_REF] Du | The Integer-Valued Autoregressive (INAR (p)) Model[END_REF], we define the dependent BINAR (p) process as follows:

Definition 2. We say that process (X t ) is dependent BINAR(p) if it has the representation:

X t = p i=1 A i (q 1,i , q 2,i ) • X t-i + t , (3.1)
where given

X t-1 = {X t-1 , X t-2 , ...}, bivariate count variables A i (q 1,i , q 2,i ) • X t-i , i = 1, ..., p
are mutually independent, and are independent of t . Moreover, A i (q 1,i , q 2,i )•X t-i is the sum of X 1,t-i independent copies of bivariate Bernoulli variable with marginal (resp. joint) probabilities α 11,i , α 21,i (resp. q 1,i ), as well as X 2,t-i independent copies of bivariate Bernoulli variable with marginal (resp. joint) probabilities α 12,i , α 22,i (resp. q 2,i ).

Thus compared with the BINAR(1) process, the extended model (3.1) has a slightly different interpretation since each individual can produce off-springs of both types at the next p periods, and these production outcomes are independent across these periods.

The stationarity condition of the BINAR(p) process is given below.

Proposition 6. Process (3.1) is both strictly and mean-variance stationary if and only if:

p i=1 α 11,i < 1, p i=1 α 22,i < 1, (3.2) and 1 - p i=1 α 11,i 1 - p i=1 α 22,i > p i=1 α 21,i p i=1 α 12,i . (3.3)
or equivalently, if and only if the eigenvalues of p i=1 A i are smaller than 1 in modulus.

The BINAR(p) process has a weak VAR(p) representation since:

E[X t+1 |X t ] = p i=1 A i X t+1-p + E[ t+1 ], (3.4)
and its conditional p.g.f. is:

E[u X1,t v X2,t |X t-1 ] = exp log b(u, v) + p i=1 X 1,t-i log a i,1 (u, v) + p i=1 X 2,t-i log a i,2 (u, v) , (3.5)
where

a i,1 (u, v) = q 1,i uv + (α 11,i -q 1,i )u + (α 21,i -q 1,i )v + (1 + q 1,i -α 11,i -α 21,i ), a i,2 (u, v) = q 2,i uv + (α 12,i -q 2,i )u + (α 22,i -q 2,i )v + (1 + q 2,i -α 12,i -α 22,i ).
Thus the conditional p.g.f. is exponential affine in X t-1 , ..., X t-p , i.e., process (X t ) is CaR of order p [CaR(p)]. Similar as the BINAR(1) process, its h-step-ahead conditional p.g.f. is still exponential affine.

Corollary 2. We have:

E[u X 1,t+h v X 2,t+h |X t ] = exp B h,0 (u, v) + p i=1 B h,i (u, v)X t+1-i , ∀h ≥ 1, (3.6)
where B h,0 and B h,i are univariate and bivariate functions, respectively. For h = 1, their values are given by equation (3.5):

B 1,0 (u, v) = log b(u, v) , B 1,i (u, v) =    log q 1,i uv + (α 11,i -q 1,i )u + (α 21,i -q 1,i )v + (1 + q 1,i -α 11,i -α 21,i ) log q 2,i uv + (α 12,i -q 2,i )u + (α 22,i -q 2,i )v + (1 + q 2,i -α 12,i -α 22,i )    ,
whereas for h > 1, we have the following recursions:

B h+1,0 (u, v) = B h,0 (u, v) + log b B h,1 (u, v) , B h+1,i (u, v) = 1 i =p B h,i+1 (u, v) + B 1,i (B h,i (u, v)), ∀i = 1, ..., p.
The proof is a direct consequence of equation (3.5) and is omitted.

Finally, the marginal mean and variance-covariance matrix of the BINAR(p) also have closed form. Their formulas are derived in Appendix A.3.

BINAR(∞) process

3.2.1 Definition, stationarity, and memory persistence A natural extension of the BINAR(p) model is to let the order p go to infinity. More precisely:

Definition 3. We say that process (X t ) is dependent BINAR(∞) if it has the representation:

X t = ∞ i=1 A i (q 1,i , q 2,i ) • X t-i + t , ∀t, (3.7)
where variables A i (q 1,i , q 2,i ) • X t-i are defined in the same way as in Definition 2.

In the above definition, the partial sum

p i=1 A i (q 1,i , q 2,i )•X t-i given X t-1
converges almost surely to a non degenerate limit when p goes to infinity. Moreover, since all the terms on the right hand side are integer-valued, variables A i (q 1,i , q 2,i )•X t-i are null for i larger than a certain stochastic threshold τ t .

To our knowledge, this is the first infinite order INAR-type model in the literature [see [START_REF] Kirchner | Hawkes and INAR (∞) Processes[END_REF] for infinite order, univariate INARCH process]. Let us first provide its stationarity condition.

Proposition 7. Process (3.1) is both strictly and mean-variance stationary if and only if:

∞ i=1 α 11,i < 1, ∞ i=1 α 22,i < 1, (3.8) and ∞ i=1 α 12,i < ∞, ∞ i=1 α 21,i < ∞, (3.9) and 1 - ∞ i=1 α 11,i 1 - ∞ i=1 α 22,i > ∞ i=1 α 21,i ∞ i=1 α 12,i . (3.10)
or equivalently, if and only if matrix ∞ i=1 A i is finite and its eigenvalues are all smaller than 1 in modulus.

This proposition nests Proposition 3 and 6. Its proof is provided in Appendix A.4.

This stationarity condition has also implications for the ranges of q 1,i , q 2,i , which have to satisfy the analogue of inequalities (2.2) and (2.3), for each i. Since sequences (α 11,i ) i and (α 21,i ) i converge to zero, for large i we have: α 11,i + α 21,i < 1. Then the ranges become:

0 ≤ q 1,i ≤ min(α 11,i , α 21,i ), 0 ≤ q 2,i ≤ min(α 12,i , α 22,i ), for large i. (3.11)
Finally, the conditional expectation and conditional p.g.f. of a BINAR(∞) process can be deduced from (3.4) and (3.6) by replacing p by ∞. Under the stationary condition, these infinite summation are all finite.

Exact simulation

While simulating a BINAR(p) process is straightforward for finite p, this is no longer the case if p = ∞. Let us now derive an exact simulation method for this latter. The basic idea is that the infinite summation in (3.7) is almost surely finite, and thus it suffices to sample the stochastic threshold τ t ∈ N such that the infinite summation actually stops at order τ t , i.e., A τt • X t-τt is non null but A i • X t-i = 0 is zero for any i ≥ τ t + 1. First we define:

δ i (X t-1 ) := (1 + q 1,i -α 11,i -α 21,i ) X1,t-i (1 + q 2,i -α 12,i -α 22,i ) X2,t-i , ∀i ≥ 1, which is equal to P[A i • X t-i = 0|X t-1 ].
We also assume without loss of generality that:

Assumption 1.

α 11,i + α 21,i < 1, α 12,i + α 22,i < 1, ∀i ≥ 1.
Indeed if the first few terms q i,1 , q i,2 do not satisfy this condition, we can leave the corresponding variables A i • X t-i out of the infinite sum and simulate them separately. This assumption implies that δ i (X t-1 ) > 0, for any i. Then the conditional CDF of the count variable τ is:

P[τ t ≤ i|X t-1 ] = ∞ j=i+1 δ j (X t-1 ) := F (i + 1|X t-1 ), ∀i ∈ N.
This CDF has the following property:

Lemma 2. Function i → F (i|X t-1 ) is nondecreasing. Its upper limit is lim i→∞ F (i|X t-1 ) = 1, whereas its lower limit F (0|X t-1 ) is strictly positive.

Proof. See Appendix A.5.

Since F (i + 1|X t-1 ) can be easily computed, we can simulate X t given X t-1 as follows:

• Draw U from the uniform distribution on [0, 1].

• Find the unique integer

τ t ≥ 0 such that F (τ t -1|X t-1 ) < U ≤ F (τ t |X t-1 )
, where by convention we set F (-1|X t-1 ) = 0. In particular, by definition X t-τt cannot be zero since

δ τt (X t-1 ) = F (τt|Xt-1) F (τt-1|Xt-1) > 1.
• Sample a certain number of independent copies of A τt • X t-τt until we obtain the first non null sample. This is possible due to Assumption 1.

• Sample random vectors A i • X t-i for i = 1, ..., τ t -1, as well as t . Then a sample of X t is given by τt i=1 A i • X t-i + t .

A constrained specification with persistent memory

To avoid the curse of dimensionality, in the application of this paper, we will focus on the following constrained BINAR(∞) specification:

α 11,i = α 11 i d , α 21,i = α 21 i d , α 12,i = α 12 i d , α 22,i = α 22 i d , q 1,i = q 1 i d , q 2,i = q 2 i d , ∀i = 1, ... (3.12)
where the power index d > 1 to ensure that i A i is finite, and for each i, probabilities q 1,i and q 2,i satisfy the constraint:

max(α 11,i + α 21,i -1, 0) ≤ q 1,i ≤ min(α 11,i , α 21,i ), max(α 12,i + α 22,i -1, 0) ≤ q 2,i ≤ min(α 12,i , α 22,i ).
It is easily checked that this is true if and only if these two inequalities hold for i = 1.

While in a BINAR(p) model, the autocovariance function decays geometrically in the lag h, a distinct feature of the BINAR(∞) model is that it allows the autocovariance to have a hyperbolic decay rate. More precisely we have the following result:

Lemma 3. In model (3.12), the autocovariance matrix

Γ(h) = E[(X t -E[X t ])(X t-h -E[X t ]) ]
of process (X t ) decays also at the hyperbolic rate d.

Proof. See Appendix A.6.

This model has a similar spirit as the (univariate) ARCH(∞) model for asset returns [see [START_REF] Giraitis | Stationary ARCH models: Dependence Structure and Central Limit Theorem[END_REF]; [START_REF] Zaffaroni | Stationarity and Memory of ARCH (∞) Models[END_REF]]. Since d > 1, the autocovariance matrix Γ(h) is summable, ruling out the possibility that [START_REF] Giraitis | Stationary ARCH models: Dependence Structure and Central Limit Theorem[END_REF] for a similar difficulty in the ARCH(∞) literature].

∞ h=0 Γ(h) = ∞ [see
We plot in the following figure the simulated path of a BINAR(∞) process with:

A =    0.12, 0.03 0.06, 0.15    , d = 1.3, q 1 = 0.015, q 2 = 0.03, ( 1t , 2t ) ∼ BP (1, 1, 0.5).
In this model, the eigenvalues of ( ∞ i=1 1/i d )A 1 are 0.7 and 0.35, respectively, thus the persistence of the process is quite strong. The simulated means and variances of the two processes are:

Ê[X 1t ] = 3.4, Ê[X 2t ] = 4.6, V[X 1t ] = 3.7, V[X 2t ] =

Predictive distributions

McCabe and Martin (2005) argue that one of the essential properties of a count process model is the tractability of the predictive p.m.f. of X t+h |X t , for both estimation and forecasting purposes.

Indeed, in terms of estimation, there are two natural approaches for the BINAR(p) model.

The first one is the Generalized Method of Moments (GMM), based on moment restrictions derived from the conditional p.g.f. [see e.g. [START_REF] Gouriéroux | Negative Binomial Autoregressive Process with Stochastic Intensity[END_REF]]. While this approach is computationally simple, it usually induces an efficiency loss. The second approach is the maximum likelihood estimation. Although more efficient, the difficulty of the latter lies in the computation of the conditional p.m.f. (X t |X t-1 ). Indeed, since this distribution is the convolution of (p+1) bivariate discrete distributions, its expression is highly cumbersome if we apply brute-force convolution [see e.g. [START_REF] Marshall | A Family of Bivariate Distributions Generated by the Bivariate Bernoulli Distribution[END_REF] for the expression of n i=1 (Z 1,j , Z 2,j ) ]. This has been identified by [START_REF] Pedeli | Some Properties of Multivariate INAR (1) Processes[END_REF] as the major downside of higher-order (B)INAR models.

As for forecasting, the conditional expectations of count processes are non integer-valued and thus incompatible with count data. While some approximation methods [see, e.g., [START_REF] Jung | Coherent Forecasting in Integer Time Series Models[END_REF]; [START_REF] Mccabe | Efficient Probabilistic Forecasts for Counts[END_REF]] have been proposed in the univariate INAR framework, they are generically time consuming and induce significant errors [see [START_REF] Lu | The Term Structure of Predictive Distributions is Solvable for Thinning-based Count Processes[END_REF] for such an analysis].

In this section we clarify that in the dependent (and, a fortiori, independent) BINAR(p) model, the conditional p.m.f. (X t+h-1 |X t-1 ) at any horizon h = 1, ... can be computed exactly.

This makes the likelihood-based estimation feasible, and eliminates the approximation error induced for the forecasting. For expository purpose we focus on BINAR(p) models with a finite, potentially large p. We first consider the horizon h = 1, before explaining how to adapt the algorithm to higher horizons.

One-step-ahead predictive distribution

Our aim in this subsection is to compute the probabilities , n|], where the bounds m, n are chosen such that the residual

P[X t = (x 1 , x 2 ) |X t-1 ] simultaneously for all couples (x 1 , x 2 ) ∈ [|0, m|]×[|0
probability P[X 1t > m or X 2t > n|X t-1 ] is negligible.
Let us first remark that for any count process, the conditional p.g.f. and p.m.f. are linked via:

E[u X1,t v X2,t |X t-1 ] = ∞ i=0 ∞ j=0 P[X t = (i, j) |X t-1 ]u i v j , ∀u, v ≥ 0. (4.1) Thus P[X t = (x 1 , x 2 ) |X t-1
] is equal to the (x 1 , x 2 )-th order coefficient in the above Taylor expansion with respect to (u, v), at (0, 0). Let us now make use of the simple conditional p.g.f.

to compute the Taylor expansion up to order (m, n). First, we rewrite (3.5) into:

E[u X1,t v X2,t |X t-1 ] = exp log b(0, 0) + p i=1 X 1,t-i log a 1 (0, 0) + p i=1 X 2,t-i log a 2 (0, 0) × exp log b(u, v) b(0, 0) + p i=1 X 1,t-i log a i,1 (u, v) a i,1 (0, 0) + p i=1 X 2,t-i log a i,2 (u, v) a i,2 (0, 0) . (4.2)
Then we perform the (m, n)-th order Taylor expansion of log[b(u, v)/b(0, 0)] and log[a i,j (u, v)/a i,j (0, 0)] at (0, 0), where i = 1, ..., p, and j = 1, 2, respectively. For most standard bivariate count distri- with p.g.f.:

b(u, v) = (1 -b 1 -b 2 ) θ /(1 -b 1 u -b 2 v) θ , ∀u, v ≥ 0 such that b 1 u + b 2 v < 1, where b 1 , b 2 , θ > 0, 1 -b 1 -b 2 > 0.
As for log ai,j (u,v) ai,j (0,0) , we have:

log a i,j (u, v) a i,j (0, 0) = ∞ k=1 (-1) k-1 k P k i,j (u, v) = m+n k=1 (-1) k-1 k P k i,j (u, v) + o m,n (u, v), ( 4.3)
where

P i,1 (u, v) = a i,1 (u, v) a i,1 (0, 0) -1 = q 1,i uv + (α 11,i -q 1,i )u + (α 21,i -q 1,i )v 1 + q 1,i -α 11,i -α 21,i , i = 1, ..., p, P i,2 (u, v) = a i,2 (u, v) a i,2 (0, 0) -1 = q 2,i uv + (α 12,i -q 2,i )u + (α 22,i -q 2,i )v 1 + q 2,i -α 12,i -α 22,i , i = 1, ..., p,
are polynomials in u and v without constant term, whereas o m,n (u, v) represents the omitted higher-order terms in the expansion. Let us explain why the truncation in (4.3) stops at order m+n. The polynomial

P k i,j (u, v) is a linear combination of terms u k1 v k2 , where k 1 +k 2 ≥ k. Thus if k > m + n, then either k 1 > m or k 2 > n,
and u k1 v k2 is omitted in the (m, n)-th order Taylor expansion. Therefore, we only need to expand recursively each P k i,j (u, v) for k = 1, ..., m + n, i = 1, ..., p and j = 1, 2 and truncate these polynomials at order (m, n). This latter can be achieved by the following algorithm:

Proposition 8. If we represent the (m, n)-th order truncation of a polynomial P :

P (u, v) = m k1=0 n k2=0 c k1,k2 u k1 v k2 + o m,n (u, v),
by the column vector:

(c 0,0 , c 0,1 , ..., c 0,n n+1 terms , c 1,0 , c 1,1 , ..., c 1,n n+1 terms , • • • , c m,0 , c m,1 , ..., c m,n n+1 terms ) ∈ R (m+1)(n+1) ,
then the (m, n)-th order truncation of polynomial P k (u, v) is represented by the column vector:

            M 0 0 0 0 • • • M 1 M 0 0 • • • • • • M 2 M 1 M 0 • • • • • • • • • • • • • • • • • • • • • M m M m-1 • • • • • • M 0             k :=M             e 1 0 • • • • • • 0            
, where column vector e 1 = (1, 0, 0,

• • • , 0) ∈ R n+1 ; square matrix M ∈ M (m+1)(n+1) (R) is (m + 1) × (m + 1) block lower triangular Toeplitz, i.e., M = 1(i ≥ j)M i-j 0≤i,j≤m
, and the block matrices M i ∈ M n+1 (R) are themselves lower triangular Toeplitz:

M i = 1(k 1 ≥ k 2 )c i,k1-k2 0≤k1,k2≤n
, ∀i = 0, ..., m.

As an illustration, if m = n = 2 we have: After Taylor-expanding log b(u, v) and log[1 + P i,j (u, v)], i = 1, ..., p, j = 1, 2, we sum up these expansions and get:

M =                          
E[u X1,t v X2,t |X t-1 ] = c exp m k1=0 n k2=0 f k1,k2 (X t-1 )u k1 v k2 + R m (u, v) = c m+n k=0 1 k! m k1=0 n k2=0 f k1,k2 (X t-1 )u k1 v k2 k + o m,n (u, v), (4.4)
where coefficients f k1,k2 (X t-1 ) are affine in X t-1 ; whereas constant c is equal to:

c = exp p i=1 X 1,t-i log(1+q 1,i -α 11,i -α 21,i )+ p i=1 X 2,t-i log(1+q 2,i -α 12,i -α 22,i )+log b(0, 0) .
In equation (4.4), the expansion stops at order m + n for the same reason as in (4.3). Then we apply Proposition 8 and obtain the (m, n)-th order truncation of polynomials

m k1=0 n k2=0 f k1,k2 (X t-1 )u k1 v k2 k , k = 1, ..., m + n.
Finally by coefficient matching we get the p.m.f. of X t |X t-1 .

In terms of the computational cost, the Taylor expansions of log b(u, v) and log[1 + P i,j (u, v)], i = 1, ..., p, j = 1, 2 are conducted only once when t varies. Thus for each t, the computation of (X t |X t-1 ) involves essentially the computation of the right hand side of equation (4.4), whose cost is independent of p. Thus this method is applicable even for large p.

The tractability of conditional distribution has several important implications. Firstly, it allows for efficient maximum likelihood estimation. In Appendix 7, we propose a small comparison between the MLE and a GMM estimator via Monte-Carlo experiments. Secondly, the likelihood function can also be used for model selection via the information criteria [see e.g. [START_REF] Weiss | An Introduction to Discrete-valued Time Series[END_REF]],

as an alternative to Box-Jenkins type approach [see e.g. [START_REF] Bu | Model Selection, Estimation and Forecasting in INAR (p) Models: A Likelihood-based Markov Chain Approach[END_REF] for an application on univariate INAR(p) models]. Thirdly, it allows to conduct likelihood ratio type tests for statistical significance of the parameters.

Comparison with the simulation-based forecasting approach

Let us now illustrate how the exact forecasting approach fares against the state-of-the-art simulationbased method [START_REF] Mccabe | Bayesian Predictions of Low Count Time Series[END_REF]; [START_REF] Jung | Coherent Forecasting in Integer Time Series Models[END_REF]] when it comes to the computation of the one-step-ahead predictive distribution. While this latter method is general and applies to other non-BINAR count process model, it will be shown that one of the advantages of BINAR models is that the exact approach outperforms significantly the simulation approach, both in terms of computational time and forecasting accuracy (note that [START_REF] Mccabe | Efficient Probabilistic Forecasts for Counts[END_REF] propose another method to approximate the predictive p.m.f. They approximate the a univariate INAR(p) process by a Markov chain with S states, where S is a large integer. This method involves the computation of matrices of dimension S 2p × S 2p , which is extremely cumbersome when p ≥ 2).

To this end we consider the BINAR(p) model with:

A i = 1 i d    0.12, 0.06 0.03, 0.15    , i = 1, ..., p, q 1 = 0.015, q 2 = 0.03,
whereas the innovation ( t ) follows the bivariate Poisson distribution BP (2, 2, 2). Given the past observation X T , the simulation-based method consists in drawing a large number of possible future values X

(n)

T +1 , n = 1, ..., N . Then the conditional p.m.f. is approximated by:

P[X T +1 = (i, j) |X T ] ≈ 1 N N n=1 1 X (n) T +1 = (i, j) , i, j.
We first report below the run time of the two methods (for the simulation-based method we consider two values for the number of draws). Both methods are implemented in R using the same laptop (intel i5, 3.0 GHz, 8GB RAM) and the program is available from the authors upon request. Table 2: CPU time of the two methods applied to BINAR(1) and BINAR(4) models. For the exact approach, we compute the value of the conditional p.m.f. ((i, j) |X T ) for any i, j ranging from 0 to 15, although in the table we only report their values for i, j non larger than 8.

We see that the exact method is 10 to 100 times faster than the simulation method. Moreover this ratio becomes larger when the order p increases. This is expected since the run time of the simulation-based method is roughly proportional to p, whereas we have argued in the previous section that it remains roughly the same in our approach.

Let us now illustrate the accuracy of the simulation-based approach. We focus on the above BINAR(1) model, and express the approximated conditional p.m.f.'s as a percentage of the corresponding exact values.

p.m.f. computed using the exact method

X 1,T +1 = 0 X 1,T +1 = 2 X 1,T +1 = 4 X 1,T +1 = 6 X 1,T +1 = 8 X 2,T +1 = 0 0.
, N = 10000 X 1,T +1 = 0 X 1,T +1 = 2 X 1,T +1 = 4 X 1,T +1 = 6 X 1,T +1 = 8 X 2,T +1 = 0 72 % 132 % 105 % 115 % 0 % X 2,T +1 = 2 111 % 101 % 104 % 104 % 93 % X 2,T +1 = 4 87 % 96 % 98 % 99 % 109 % X 2,T +1 = 6 85 % 103 % 106 % 108 % 91 % X 2,T +1 = 8 0 % 89 % 91 % 108 % 101 % Relative accuracy of the simulation-based method, N = 100000 X 1,T +1 = 0 X 1,T +1 = 2 X 1,T +1 = 4 X 1,T +1 = 6 X 1,T +1 = 8 X 2,T +1 = 0 118 % 99 % 105 % 92 % 65 % X 2,T +1 = 2 101 % 97 % 101 % 96 % 115 % X 2,T +1 = 4 96 % 98 % 100 % 100 % 100 % X 2,T +1 = 6 99 % 96 % 101 % 97 % 98 % X 2,T +1 = 8 80 % 95 % 97 % 105 % 100 %
Table 3: Relative accuracy of the simulation approach compared to the exact approach. Due to the space constraint we have only reported the results for odd values of X 1,T +1 and X 2,T +1 .

We can see that the approximation error of the simulation approach is substantial, even with a huge number of draws (N = 100000). This is particularly the case for the probabilities of "extreme events", that is when either X 1,t+1 or/and X 2,t+1 is large. This is a serious downside of the standard simulation approach, since in finance, predicting extreme events is key to the risk management.

Multiple-step-ahead predictive distributions

Let us now adapt the above algorithm for the computation of (X t+h |X t ). For expository purpose, we focus on the BINAR(1) process. To Taylor-expand the conditional p.g.f. given in Proposition 1, we use the following procedure:

• First, we use Proposition 8 to compute the (m, m)-th order Taylor expansion of log a

(h) 1 (u, v), log a (h)
by Proposition 1, a (h) 1 and a (h) 2 are 2 h -th order polynomial in (u, v), which can be large for large h, only the terms of degree lower than min(m, 2 h ) contribute to the (m, m)-th Taylor expansion. Thus this step involves roughly the same computational effort for different values of h.

• Then we compute the (m, m)-th order Taylor expansion of

exp X 1,t log a (h) 1 (u, v) + X 2,t log a (h) 2 (u, v) + log b (h) (u, v) = c exp X 1,t log a (h) 1 (u, v) a (h) 1 (0, 0) + X 2,t log a (h) 2 (u, v) a (h) 2 (0, 0) + log b (h) (u, v) = c 2m k=0 1 k! X 1,t log a (h) 1 (u, v) a (h) 1 (0, 0) + X 2,t log a (h) 2 (u, v) a (h) 2 (0, 0) + log b (h) (u, v) b(0, 0) k + o m,m (u, v), where c = exp X 1,t log a (h) 1 (0, 0) + X 2,t log a (h)
2 (0, 0) + log b(0, 0) . This is conducted using the same method as for h = 1. Finally we deduce P[X t+h = (i, j) |X t ] for any (i, j) ∈ [|0, m|] 2 by coefficient matching.

Application

The mutual fund industry

Mutual funds (MF) are investment vehicles who invest in a wide range of assets ranging from liquid ones such as stocks, bonds, to highly illiquid ones such as hedge funds. Their clients include, for instance, insurance companies, private banks, large corporations as well as retail investors (In some funds, including the one we study, retail investors' orders are first centralized by a broker before being transferred to the fund. Thus from the fund's point of view, its client is the broker). They are traditionally much less regulated than commercial banks, but this potential loophole has recently received much attention of the regulator.

Most MF are open-ended, i.e., they allow investors to purchase new shares, or redeem their shares on a daily basis. Thus the size of the MF can feature important short-term fluctuations, making them vulnerable to liquidity risk. In particular, during a market turmoil, investors' redemption decisions tend to cluster. If the fund manager's cash holding is insufficient to meet the redemption requests, he might be foced to sell its illiquid assets, whose market liquidity would have also plunged due to the crisis. Such fire selling usually leads to significant investment loss, which in turn creates panics and triggers further redemptions. This phenomenon is called fund run. On the other hand, while inflow, that is the purchase orders of the MF, or cash holding can offset the outflow due to redemption, a sudden large inflow also dilutes the investment performance of the fund due to the lack of immediate investment opportunities. Thus they can also trigger subsequent (large) outflows. As a consequence, it is essential for the fund manager to monitor parallely the outflow and the inflow pattern of its clients, on a daily basis. The current MF literature usually focuses on the outflow only, or the net outflow, that is the difference between the outflow volume and inflow volume [see, e.g., [START_REF] Schmidt | Runs on Money Market Mutual Funds[END_REF]; [START_REF] Desmettre | Generalized Pareto Processes and Liquidity[END_REF]]. Moreover, many of these studies are based on weekly data only.

Finally, while prior studies focus on the volume of the outflow/inflow, our attention is on the number of purchase X 1,t and redemption orders X 2,t . These variables are closely related since on each trading day, the volume Y 1,t and Y 2,t of the outflow/inflow have the compound representation:

Y 1,t = X2,t j=1 S j,t , Y 2,t = X2,t j=1 B j,t , (5.1)
where S j,t (resp. B j,t ) denotes the volume of the j-th redemption (resp. purchase). Moreover, when j and t vary, (S j,t ) and (B j,t ) can be reasonably assumed to i.i.d.

In this framework, our preference for studying the counts rather than the volume is motivated by the following reasons. First, very often a large outflow volume is due to the redemption order of a large investor. These "VIP" clients usually have privileged relationship with the fund manager and in the case of a large redemption, they also tend to (privately) inform the fund manager sufficiently in advance so that the latter can avoid a massive fire selling. On the other hand, a large X 2,t spells a collective withdrawal, that is the fund run, which is the most dangerous scenario. More importantly, the non-linear forecasting of the net outflow, that is Y 2,t -Y 1,t , can be easily deduced from representation (5.1). Indeed, the conditional Laplace transform of

Y 2,t -Y 1,t given X t-1 is: E[e -u(Y2,t-Y1,t) |X t-1 ] = E (E[e uB1,t ]) X2,t (E[e -uS1,t ]) X1,t |X t-1 , ∀u ∈ R. (5.2)
The right hand side is the conditional p.g.f. of the count process evaluated at (E[e uS1,t ], E[e -uB1,t ]), which has closed form under suitable distributional assumptions of S 1,t and B 1,t (such as gamma).

Thus the left hand side of (5.2) is readily available. Then the conditional Value-at-Risk of the net outflow, say, can be accurately approximated without simulation [see [START_REF] Gordy | Saddlepoint Approximation of CreditRisk +[END_REF]]. As a consequence, in the paper we will only focus on counts.

Data description

Our dataset comes from a French equity-focused MF. We observe the daily number of purchase orders X 1t and redemption orders X 2t , during about 1000 trading days. The following three figures plot the trajectory of the two count processes during a sub-sample of 100 trading days, as well as the histograms of X 1t and X 2t . We can see that the probability of X 1,t and X 2,t taking large values is extremely small. We report below the empirical marginal moments of the two processes. Both processes feature mild unconditional over-dispersion, a typical feature of BINAR processes. Figure 5 plots their autocorrelation patterns. These ACF's and CCF decay rather slowly, resembling the simulated patterns of section 3.2.3 for a model with hyperbolic decaying coefficients A i , q 1,i , q 2,i . Thus in this section, we will focus on the estimation of this latter model. Note that we can interpret the thinning part of model (2.1) as trades by clients who also made purchase/sell orders in the previous periods. Then the binomial distributional assumption can be explained by the fact that usually clients only trade on a daily basis (even if they place several orders within a day, usually the fund only executes the orders once each day at the market closure, and thus only one aggregate order is counted for each client). Thus this interpretation suggests model (2.1) with the further constraint that q 1 = q 2 = 0, that is, clients do not make simultaneous buy and sell orders as these two will partially cancel out. In practice, in the application, although this constrained version of model (2.1) is easier to interpret, we have considered a BINAR(∞) as this latter reflects better the persistence of the ACF function.

Ê[X 1t ] Ê[X 2t ] V[X 1t ] V[X 2t ] ĉorr[X 1t , X 2t ] 2.
Since the CCF indicates a positive (marginal) correlation between the two processes, it is reasonable to assume that the innovation process ( t ) features also positive correlation between its two components. Moreover, the degree of unconditional over-dispersion of the two count variables are rather weak (see Table 4), suggesting the Poisson distribution for the marginal distribution of t . Thus we assume the distribution of the latter to be bivariate Poisson BP(λ 1 , λ 2 , λ 3 ). Then the set of parameters is: θ = (α 11 , α 22 , α 12 , α 21 , q 1 , q 2 , λ 1 , λ 2 , λ 3 , d).

Let us now interpret the regression parameters α i,j . Parameter α 11 measures how purchase decisions of investors are (positively) correlated, due to the so-called reputation effect. Moreover, since here the investors' benchmark is past yearly performance and our observations concern daily inflow/outflow movements, it is not surprising that this reputation effect decays rather slowly when daily data is used, as is shown by the ACF of X 1t . Similarly, parameter α 22 measures the panic effect among investors. Parameter α 12 captures the propensity of redemption following large recent inflows, as such inflows might dilute the fund's performance due to lack of sufficient investment opportunities. This is consistent with the CCF given in Figure 5, which indicates a positive cross correlation between X 2t and the lagged values of X 1,t . Finally, parameter α 21 captures the propensity of purchase following large outflows. This can be interpreted as the fund manager's capability of attracting new investment in order to stabilize the fund size.

Model estimation

We estimate the parameters by maximum likelihood. A practical difficulty of applying an infinite order model is that our number of observations (which is around 1000 trading days) is finite.

Thus the "∞" in ∞ i=1 A i • X t-i should be replaced by a finite, but large p. In the estimation we take p = 300, and regard the first p values of X t as initial values rather than observations. We have checked that for reasonable values of parameters, further increasing p only marginally impacts the value of the likelihood function. Table 5: Parameter estimate along with the standard errors of the estimators.

The parameter estimate α22 (= 0.285) is larger than α11 (= 0.174). In other words, the panic effect is more important than the reputation effect. Second, α12 (= 0.12) is much larger than α21 (= 0.0355), which means that existing investors are quite sensitive to large inflows and tend to redeem their shares for fear of performance drop after a large inflow. This highlights the importance of monitoring the purchase counts separately from the redemption counts. On the other hand, it seems to be difficult for the fund to attract new investors after large outflows.

Finally, the two eigenvalues of the matrix ( ∞ i=1 1/i d )A 1 are 0.58 and 0.26, respectively, which are both smaller than 1. Thus the joint process seems to be stationary.

Let us now compute the conditional p.m.f. We focus on horizon 1, since the fund manager usually adjusts the positions on a daily basis. For each trading day t in our observation period, we compute P[X 1,t+1 = m, X 2,t+1 = n|X t ] for all m, n ≤ M , where we set M to be the largest past observation M = max t,j X j,t = 17.Then we follow [START_REF] Mccabe | Bayesian Predictions of Low Count Time Series[END_REF] and compute the conditional mode ( X1,t+1 , X2,t+1 ) defined by:

( X1,t+1 , X2,t+1 ) = arg max 0≤m,n≤M P[X 1,t+1 = m, X 2,t+1 = n|X t ].
Figure 6 displays the evolution of the conditional mode X2,t against the corresponding realized value X 2,t for all the past dates. Globally, we can see that the mode forecast can quite well capture the local tendency of the count process, although the realized paths tends to be more erratic.

Diagnostic check Let us finally conduct some adequacy checks of the estimated model. We first compute the ACF/CCF of the estimated model, as well as some summary statistics of the empirical Pearson residuals [see e.g. [START_REF] Weiss | An Introduction to Discrete-valued Time Series[END_REF] for their definitions]. For the first one, since in the paper we have only derived the ACF/CCF for finite p, we resort to Monte-Carlo simulation to obtain approximations of these functions for this BINAR(∞) model. Due to space constraint this figure is provided in Appendix 8. We can remark that globally, the ACF/CCF of the estimated process are quite similar to their empirical counterparts reported in Figure 5. Moreover, the Pearson residuals seem to be well uncorrelated across different lags.

Conclusion

We have extended the BINAR(1) model to allow for dependent thinning, arbitrary errors, and higher-order dynamics. This family has intuitive interpretations, tractable stationarity properties, and are rather flexible compared to existing models. More importantly, we have derived tractable expressions for the predictive distributions, allowing for likelihood based estimation and non-linear forecasting. The model has been applied to a new application area, i.e., fund liquidity risk.

In the paper we have followed the literature by focusing on bivariate models. Is it possible

   V[X 1,t ], * Cov[X 1,t , X 2,t ], V[X 2,t ]    =    V[α 11 X 1,t-1 + α 12 X 2,t-1 ], * Cov[α 11 X 1,t-1 + α 12 X 2,t-1 , α 21 X 1,t-1 + α 22 X 2,t-1 ], V[α 21 X 1,t-1 + α 22 X 2,t-1 ]    + E    V[X 1,t |X t-1 ], * Cov[X 1,t , X 2,t |X t-1 ], V[X 2,t |X t-1 ]    = A    V[X 1,t ], * Cov[X 1,t , X 2,t ], V[X 2,t ]    A + γ + E[X 1t ]γ 12 + E[X 2t ]γ 34 .
Solving this linear matrix equation yields solution (2.21). Finally, equation (2.22) is a direct consequence of equation (2.4).

A.3. The first two marginal moments of BINAR(p) model

Under the stationarity condition, the marginal expectation of a BINAR(p) model satisfies:

E[X t ] = p i=1 A i E[X t ] + E[ t ] ⇐⇒ E[X t ] = (Id - p i=1 A i ) -1 E[ t ].
where matrix Id -p i=1 A i is invertible by conditions (3.2) and (3.3). The covariance matrix can be obtained using the companion CaR(1) form. More precisely, let us denote by V the (2p × 2p) covariance matrix of (Y t ) = (X t , X t-1 , ..., X t-p+1 ), i.e.:

V := E (Y t -E[Y t ])(Y t -E[Y t ]) =          Γ(0) Γ(1) • • • Γ(p -1) Γ (1) Γ(0) • • • Γ(p -2) • • • • • • • • • • • • Γ(p -1) Γ(p -2) • • • Γ(0)         
, where Γ(h) is the auto-covariance function of (X t ) at lag h. Since Y t has the VAR(1) representation:

E[Y t |Y t-1 ] =          A 1 A 2 • • • A p I 2 0 • • • • • • 0 I 2 • • • • • • • • • • • • • • • • • •          :=A Y t-1 +          E[ t ] 0 • • • • • •          ,
we have:

V = AV A +          γ + p i=1 E[X 1,t ]γ 12,i + X 2,t γ 34,i 0 • • • 0 0 0 0 0 • • • • • • • • • • • • 0 0 0 0          , (eq.a.3)
where γ 12,i , γ 34,i are defined in a similar way as in (2.7). Thus V = ∞ h=0 (A ) h V 0 A h , where V 0 is the second matrix term on the right hand side of (eq.a.3).

A.4. Proof of Proposition 7

We first reformulate the condition given in Proposition 7 in terms of the eigenvalues of A: Lemma 4. For any matrix A = (α i,j ) 1≤i,j≤2 with nonnegative entries only, the two following sets of conditions are equivalent:

1. α 11 < 1, α 22 < 1, and (1 -α 11 )(1 -α 22 ) > α 12 α 21 .
2. The eigenvalues of A are smaller than 1 in modulus.

The proof is straightforward and omitted [see e.g. [START_REF] Gouriéroux | Negative Binomial Autoregressive Process with Stochastic Intensity[END_REF]].

Let ( t ) be an i.i.d. sequence of innovations and (Z 1,t,i,j , Z 2,t,i,j ), (Z 3,t,i,j , Z 4,t,i,j ) i.i.d. bivariate Bernoulli variables for any t ∈ Z, i, j ∈ N. Then we define the doubly indexed sequences A.5.2. Strict and second-order stationarity of process (X t ) t By definition, for each n, process (X (n) t ) t is strictly stationary. Thus the limiting process (X t ) is also strictly stationary. By the L 1 (P) convergence, (X t ) is mean-stationary. To show the covariance stationarity, let us check that the sequence (X (n) 1t ) n is bounded in L 2 (P). We have:

V[X (n) 1t ] = V[ 1t ] + V n i=1 X (n-1) 1,t-i j=1 Z 1,t,i,j + n i=1 X (n-1) 2,t-i j=1 Z 3,t,i,j ≤ V[ 1t ] + E n i=1 X (n-i) 1,t-i j=1 Z 1,t,i,j + n i=1 X (n-i) 2,t-i j=1 Z 3,t,i,j 2 ≤ V[ 1t ] + E n i=1 X (n) 1,t j=1 Z 1,t,i,j + n i=1 X (n) 2,t j=1 Z 3,t,i,j 2 ≤ C 1 + V ( ∞ i=1 α 11,i )X (n) 1t + ( ∞ i=1 α 12,i )X (n) 2t ,
where constant C 1 is independent of t and n. Similar upper bounds can also be obtained for

V[X (n) 1t ] and V[X (n)
2t ] and in matrix form we have:

V n ≤ A ∞ V n A ∞ + C,
where V n is the covariance matrix of X (n) t , and C is a constant matrix. Thus we get:

V n ≤ ∞ i=0 A i ∞ C(A ∞ ) i ,
which is uniformly bounded. Thus by the dominated convergence theorem, (X (n) t ) n also converges to (X t ) in L 2 (P).

first briefly recall the GMM method, which has recently been suggested by [START_REF] Gouriéroux | Negative Binomial Autoregressive Process with Stochastic Intensity[END_REF] for CaR count processes.

A.7.1. The GMM method

We denote by U a grid of positive real numbers, then the conditional pgf provides a set of moment constraints:

E ψ u,v (X t-1 ) u X1,t v X2,t -a 1 (u, v) X1,t-1 a 2 (u, v) X2,t-1 b(u, v) =E[u X 1,t v X 2,t |Xt-1] = 0,
∀u, v ∈ U, (eq.a.6)

where ψ u,v (X t-1 ) can be any instrumental function depending only on X t-1 . We choose ψ to be:

ψ u,v (X t-1 ) = 1 V[u X1,t v X2,t |X t-1 ]
,

where

V[u X1,t v X2,t |X t-1 ] = E[u 2X1,t v 2X2,t |X t-1 ] -E[u X1,t v X2,t |X t-1 ] 2 = a 1 (u 2 , v 2 ) X1,t-1 a 2 (u 2 , v 2 ) X2,t-1 b(u 2 , v 2 ) -a 1 (u, v) 2X1,t-1 a 2 (u, v) 2X2,t-1 b 2 (u, v),
such that the integrand in equation (eq.a.6) has unitary unconditional variance.

Then a simple GMM estimator is obtained by minimizing the error of the empirical counterparts of these orthogonality conditions:

θT = arg min θ 1 T -1 T t=2 g(X t-1 , X t , θ) W 1 T -1 T t=2 g(X t-1 , X t , θ) ,
where g is a vector function of dimension Card(U) 2 given by:

g ψ,u,v (X t-1 , X t ) = ψ u,v (X t-1 ) u X1,t v X2,t -a 1 (u, v) X1,t-1 a 2 (u, v) X2,t-1 b(u, v) , ∀u, v ∈ U,
and W is a symmetric positive definite weighting matrix. To ensure the asymptotic consistency, the number of moment conditions [Card(U)] 2 has to be larger than the number of parameters.

Then, under standard regularity conditions, the GMM estimator is consistent and asymptotically We can remark that, first, since the sample size T is only moderately large, the distribution of the estimators are quite different from the limiting Gaussian distribution. In particular, due to the constraint that all components of the estimators have to be bounded from below (by zero), and most of them (including α 11 , .., α 22 , q 1 , q 2 are also from above by 1, their distribution can be either left skewed, or right skewed. Secondly, by comparing the above two figures, we can remark a slight better finite sample behavior for the MLE estimator, in the sense that for most components, the range of the MLE estimator is (slightly) more concentrated than for the GMM estimator. This is expected, since in most applications the MLE is asymptotically more efficient than a GMM estimator.

A.8. Additional figures

Figure 1 :

 1 Figure 1: Simulated trajectory of the BINAR(∞) process under constraint (3.12).

Figure 2 :

 2 Figure 2: ACF and CCF function of the same simulated BINAR(∞) process as in Figure 1. As is expected, both functions decay rather slowly.

  p = 4, X T = X t-1 = X T -2 = X T -3 = (1

Figure 3 :Figure 4 :

 34 Figure 3: Joint trajectory of the purchase count (red dashed line) and redemption count (black full line) during 100 trading days.

Figure 5 :

 5 Figure 5: ACF's and CCF of inflow count process (X 1,t ) and outflow process count (X 2t ).

Figure 6 :

 6 Figure 6: Joint evolution of X2,t+1 (in black full line) and the realized redemption count X 2,t+1 (in red dashed line).

  Figure A.7.1: Boxplots of the GMM estimators of the BINAR(1) model (2.1) with parameter values given by (A.7.1). The real value of each component of θ is plotted in red dotted line.

FigureFigure

  Figure A.8.1: ACF/CCF of the estimated process.

  or leave the queues after being served, with probabilities α 11 , α 21 and 1 -α 11 -α 21 (resp. α 22 , α 12 and 1 -α 22 -α 12 ).

	•. Moreover, this constrained
	model has also a similar queuing interpretation as the univariate INAR(1) model [see e.g. Schweer
	and Wichelhaus (2015)]. More precisely, we can think of X 1,t , X 2,t as the number of individuals
	in queues 1 and 2 at date t, respectively. Both queues have infinite capacity. At date t + 1, 1,t+1
	(resp. 2,t+1 ) new customers join queue 1 (resp. queue 2), whereas X 1,t (resp. X 2,t ) customers
	that were in queue 1 (resp. queue 2) at date t can either stay in the same queue, or go to the
	other queue,

Table 1 :

 1 Comparison of over-dispersion coefficients and correlation coefficients of the three models.

	Model 1	1.95	1.80	0.84
	Model 2	1.47	1.40	0.55
	Model 3	1	1	0

Table 4 :

 4 Summary empirical moments of the two processes.

	63	3.06	3.67	3.97	0.27

  The following table reports the maximum likelihood

		Parameter estimate
	α 11	0.174 (0.055)
	α 21	0.0355 (0.016)
	α 12	0.120 (0.052)
	α 22	0.285 (0.12)
	q 1	0.0219 (0.012)
	q 2	0.0173 (0.0073)
	λ 1	0.712 (0.21)
	λ 2	0.695 (0.16)
	λ 3	0.354 (0.18)
	d	1.83 (0.1)
	estimates:

(u, v) and b (h) (u, v) at (u, v) = (0, 0), where m is chosen such that the probability of either component processes taking values larger than m is negligible.Note that although

Acknowledgement

We thank an associate editor as well as two referees for insightful comments.

to extend our model into higher dimensions? From the fund manager's point of view, this can be of interest since different investors-large corporates and bank/insurance companies, say, can have different behaviour. If count series, two for each client category, becomes available, then a multivariate analysis allows to study cross sectional dynamic effects between different clients and may improve the quality of forecasts. The answer is (partially) affirmative, since equation (2.1) can be extended to the multivariate case, using the multivariate Bernoulli distribution [START_REF] Chaganty | Range of Correlation Matrices for Dependent Bernoulli Random Variables[END_REF]]. However such extensions are not without downsides. Firstly, this extended model is not closed under margins; for example the bivariate margins of the trivariate extension no longer have representation (2.1). Secondly and most importantly, as many multivariate models, when the dimension increases, both the number of parameters and the computational burden increases. These issues might be mitigated by introducing constrained specifications (see, e.g., Proposition 2), or by conducting pairwise analysis [see Pedeli and Karlis (2013a); [START_REF] Gouriéroux | Negative Binomial Autoregressive Process with Stochastic Intensity[END_REF]]These await future research.

Appendix

A.1. Proof of Proposition 4

Under the assumptions of the proposition, we have:

. By equation (2.12) we get:

which is the log p.g.f. of a bivariate independent Poisson distribution, with expectations:

Let us now check that (Y t ) = (X 1,t + X 2,t ) follows a Poisson INAR(1) process. The joint p.g.f. of (X 1,t-1 + X 2,t-1 , X 1t + X 2t ) is:

=λ1+λ2 by equation (2.20) distributed, the joint p.g.f. is:

) .

(eq.a.2)

By matching equations (eq.a.1) and (eq.a.

2), we conclude that

A.2. Proof of Proposition 5

The marginal expectation is obtained by taking expectation in equation (2.4). By the co-variance decomposition formula we get:

2,t ) t recursively by:

We will show that the n-indexed sequence (X

) n converges almost surely to a limit X t , and that the limiting process (X t ) satisfies the definition (3.7).

A.5.1. Almost sure and L 1 (P) convergence of X

(n) t

By induction (with respect to n) it is easily checked that:

Thus for each t, both n-indexed sequences (X

) n converge almost surely, to limits, say, X 1,t and X 2,t , respectively. Let us now consider the L 1 (P) convergence. We have:

where the inequalities hold componentwise. Thus we have:

Then we remark that the (2 × 2) matrix (Id -A ∞ ) is invertible and all the entries of its inverse are nonnegative. Thus we can multiply both sides by (Id -A ∞ ) -1 and deduce that:

is upper bounded when n increases. Hence (X (n) t ) converges also in L 1 (P) to X t for each t by monotonous convergence theorem.

A.5.3. Conditional distribution of X t given its past

It remains to check that the above limiting process (X t ) satisfies the representation (3.7). It suffices to show that, for a fixed t, the sequence

converges to X t almost surely (a.s.), or equivalently, since

converges to zero a.s. But this sequence is non-decreasing, thus it suffices to find an a.s. convergent subsequence. Its L 1 (P) norm is:

when n goes to infinity. Thus (r

and admits a subsequence that is a.s. convergent. Thus process (X t ) satisfies equation (3.7).

A.5.4. The necessary condition of stationarity

Thus under conditions (3.8), (3.10), (3.9), the BINAR(∞) process exists and is stationary. Let us now show that these conditions are also necessary. Taking expectation in (3.7), we get:

(eq.a.4) thus all the entries of A ∞ are finite, hence inequality (3.8). Finally, by iteration we have:

Thus the largest eigenvalue of A ∞ is smaller than one in modulus. Then by Lemma 4 we get conditions (3.9) and (3.10).

A.5. Proof of Lemma 2

Under Assumption 1, each δ i (X t-1 ) lies between 0 and 1, thus F (i|X t-1 ) is nondecreasing, and upper bounded by 1. Its lower limit is strictly positive under the stationarity condition since:

(eq.a.5)

As α 11,i + α 21,i → 0 and α 12,i + α 22,i → 0 when i goes to infinity, for large i, we have:

Then since

)X 2,t-i is finite, the right hand side of (eq.a.5) is positive.

A.6. Proof of Lemma 3

The proof is based on the fact that process (X t ) has the weak VAR(∞) representation:

where η t is a weak white noise. Then we revert this VAR(∞) representation into the Vector MA(∞) representation X t = (Id + ∞ i=1 B i L i )η t , where B j are matrices and L is the lag operator. By mimicking the proof of Theorem 2 in [START_REF] Zaffaroni | Stationarity and Memory of ARCH (∞) Models[END_REF], we can show thatB i = O(i d )D for some constant matrix D and the autocovariance function also decays at the same hyperbolic rate i d .

A.7. Comparison with the GMM and ML estimators

For illustrative purpose, in this section we focus on a BINAR(1) model with independent, Poisson innovations, and compare the small sample behavior of the GMM and ML estimators. Let us normally distributed:

where θ 0 is the true parameter value, G = E[∇g(X t-1 , X t , θ 0 )], with ∇ representing the differential with respect to argument θ, and Ω = E[g(X t-1 , X t , θ 0 )g(X t-1 , X t , θ 0 ) ]. Moreover, the optimal choice of the weighting matrix W is given by W = Ω -1 , and for this choice we have:

Thus the asymptotic behaviour of the GMM estimator will usually depend on the choice of the set U. As far as we know, this issue has yet to be considered for count data and thus is well beyond the scope of the present paper. In the simulation exercise below, we take: U = {0.11, 0.22, ..., 0.99, 1.1}, which corresponds to a total of [Card(U)] 2 = 121 moment conditions, which is well beyond the typical number of parameters of a BINAR(p) model.

A.7.2. A Monte-Carlo comparison of GMM and ML estimators

In this subsection we consider the following dependent BINAR(1) model with independent, Poisson P(λ 1 ) (resp. P(λ 2 )) distributed error terms 1,t and 2,t . The parameters of the model are set as:

Then we simulate m = 100 independent trajectories of the process, each with a sample size of T = 500. For each simulated dataset j = 1, ..., m, we estimate the model using both the GMM and ML methods to obtain two estimators θ1,j and θ2,j . Boxplots of the empirical distribution of the parameters ( θ1,j ) j and ( θ2,j ) j are given in the following figures.