Next Article in Journal
New York FED Staff Nowcasts and Reality: What Can We Learn about the Future, the Present, and the Past?
Previous Article in Journal
Temperature Anomalies, Long Memory, and Aggregation
Previous Article in Special Issue
Maximum-Likelihood Estimation in a Special Integer Autoregressive Model
Article

Goodness–of–Fit Tests for Bivariate Time Series of Counts

1
Department of Probability and Mathematical Statistics, Faculty of Mathematics and Physics, Charles University, Sokolovská 83, 186 75 Prague 8, Czech Republic
2
Department of Economics, National and Kapodistrian University of Athens, 105 59 Athens, Greece
3
Unit for Pure and Applied Analytics, North–West University, Potchefstroom 2531, South Africa
*
Authors to whom correspondence should be addressed.
Received: 4 December 2020 / Revised: 25 February 2021 / Accepted: 26 February 2021 / Published: 4 March 2021
(This article belongs to the Special Issue Discrete-Valued Time Series: Modelling, Estimation and Forecasting)

Abstract

This article considers goodness-of-fit tests for bivariate INAR and bivariate Poisson autoregression models. The test statistics are based on an L2-type distance between two estimators of the probability generating function of the observations: one being entirely nonparametric and the second one being semiparametric computed under the corresponding null hypothesis. The asymptotic distribution of the proposed tests statistics both under the null hypotheses as well as under alternatives is derived and consistency is proved. The case of testing bivariate generalized Poisson autoregression and extension of the methods to dimension higher than two are also discussed. The finite-sample performance of a parametric bootstrap version of the tests is illustrated via a series of Monte Carlo experiments. The article concludes with applications on real data sets and discussion.
Keywords: goodness-of-fit; bivariate time series; probability generating function; time series of counts goodness-of-fit; bivariate time series; probability generating function; time series of counts

1. Introduction

Time series of counts enjoy numerous applications in such diverse fields as business, economics, engineering, and epidemiology, and corresponding inferential procedures have been intensively studied in recent years. The reader is referred to the earlier work of McKenzie (2003), as well as to the much updated full-book treatment of Davis et al. (2015) and Weiss (2018a) for an overview of the subject. Among the most popular models for count-time series are the integer autoregression (INAR) model and the Poisson autoregressive (PAR) model, also termed Poisson INARCH model. The INAR and PAR models, originally conceived for univariate counts, have been extended in order to accommodate bivariate, and more generally multivariate counts; see Latour (1997) and Liu (2012), respectively. One of the basic parametric elements of these generalizations involves moving from a univariate to a bivariate distribution, and eventually to a multivariate one.
In this connection, when confronted with a vectorial data-set of time series of counts, and since one is presented with a choice of several possible candidate models, it is extremely important to check model-adequacy via some goodness-of-fit (GOF) test. Otherwise a poorly fitted model might yield misleading inference. Various such procedures have been proposed in the literature for univariate series of counts. A brief overview of these approaches is provided in Section 2.
Inspired by the univariate test criteria, we propose GOF tests for bivariate INAR and PAR models. Our tests mainly target the (by far most popular) Poisson specification of these models, but also involve other structural aspects, such as order and functional specification of the underlying model. The suggested test statistic is constructed as a contrast between an estimator of the probability generating function (PGF) which is entirely model-free and a semiparametric counterpart that “respects” the model under test. Hence each new test may be viewed as a bivariate extension of the earlier PGF-based procedure suggested by Meintanis and Karlis (2014), Hudecová et al. (2015) and Schweer (2016).
The remainder of this work unfolds as follows. Section 2 provides a review of the models considered and of goodness-of-fit testing for univariate time series of counts. In Section 3 and Section 4 we introduce the bivariate models and the new test criteria. In Section 5 the asymptotic properties of the methods are investigated and resampling versions of the tests are proposed. An extension of one the tests to a more complicated model is discussed in Section 6. The finite-sample properties of the new criteria are studied by means of Monte Carlo methods in Section 7. Applications to real data sets are included in Section 8 while our final thoughts regarding the news methods and potential extensions thereof are summarized in Section 9. Proofs of asymptotics are provided in Appendix A.

2. Goodness-of-Fit Methods for Univariate Time Series of Counts

Many of the models for time series of low counts assume that conditionally on the past, the distribution of the current variable is fully specified by a family of laws indexed by a certain parameter. An important example is a class of integer autoregressive conditionally heteroscedastic (INGARCH) models, covering PAR models as special cases. Different models are based on the thinning operator see Steutel and van Harn (1979), with the INAR model being the most popular one. In the following, a count distribution refers to a discrete distribution on N 0 = N { 0 } and a count variable is a random variable with such distribution.
Let { Y t } = { Y t , t Z } be a univariate time series and let F t = σ { Y s , s t } is the information known up to time t. An integer GARCH model, abbreviated as nonlinear INGARCH( p , q ), is defined as
Y t | F t 1 F ( λ t ) , λ t = r ( Y t 1 , , Y t p , λ t 1 , , λ t q , η ) ,
where F ( λ ) is some count distribution with mean λ , and the regression function r : [ 0 , ) p + q × Θ belongs to some specific parametric family of functions R = { r ( · , η ) : η Θ R k } for some k 1 . If F is Poisson and R is a family of linear functions then this model is referred to as Poisson (linear) INGARCH, see Ferland et al. (2006). Some authors also use the name Poisson autoregression for cases where F is Poisson and r could be nonlinear, see, e.g., Fokianos and Tjøstheim (2009), Fokianos and Tjøstheim (2012), Fokianos (2012), for the case p = q = 1 . If q = 0 then the model has a purely autoregressive structure and is abbreviated as INARCH ( p ) . In the following the acronym PAR(p) is used for INARCH(p) with F Poisson and r linear. It has been shown that if p = q = 1 and r is linear, i.e., r ( x , η ) = η 1 + η 2 x 1 + η 3 x 2 , and if η i 0 and η 2 + η 3 < 1 and F belongs to the single-parameter exponential family of distributions (that includes the Poisson distribution as a special case), then there exists a strictly stationary and ergodic solution of (1), see Davis and Liu (2016). For an overview of conditions for strict stationarity and ergodicity for other choices of F see e.g., (Ahmad and Francq 2016, Section 3).
A different class of models consists of integer autoregressive moving average (INARMA) models. These models arise from the same structure as the classical ARMA time series models, but the multiplication sign is replaced by the Steutel and van Harn’s thinning operator ∘. If Y is a count random variable and α ( 0 , 1 ) then
α Y : = i = 1 Y U i ,
where { U i } are iid Bernoulli variables with α = 𝖯 ( U i = 1 ) = 1 𝖯 ( U i = 0 ) , which are independent of Y, with the convention that an empty sum (the case Y = 0 ) equals 0. Let { ε t } be a sequence of iid count random variables with distribution G with a finite variance, and let α i , β j ( 0 , 1 ) for i = 1 , , p , j = 1 , , q . The INARMA ( p , q ) model is defined as
Y t = i = 1 p α i Y t i + ε t + j = 1 q β j ε t j ,
where the Bernoulli variables involved in all the thinning operations are independent and independent of { ε t } . If p = 1 , q = 0 then the model (3) corresponds to the INAR(1) model introduced in McKenzie (1985) and Al-Osh and Alzaid (1987), taking the following form
Y t = α Y t 1 + ε t = i = 1 Y t 1 U t , i + ε t .
For α ( 0 , 1 ) there exists a strictly stationary solution of (4) and the law of the innovations { ε t } uniquely determines the marginal distribution of Y t as well as the conditional distribution Y t | F t 1 . In particular, if { ε t } are iid Poisson then Y t has Poisson distribution as well, and this special case has been considered in many applications. Model (3) with p > 1 , q = 0 , was introduced and studied by Du and Li (1991) and since then, many authors have considered variants of model (3), various extensions and modifications, see Scotto et al. (2015) for a comprehensive review.
Various GOF tests have been suggested in the literature for the aforementioned two classes of models. Neumann (2011) and Fokianos and Neumann (2013) considered GOF tests for the regression function r in a Poisson INGARCH(1,1); see model (1) with p = q = 1 and Poisson F. A slightly less formal assessment of model adequacy is explored in Aleksandrov and Weiss (2020) for a PAR(1) model as well as for a Poisson INAR(1). GOF tests based on the sample index of dispersion were considered in Schweer and Weiss (2014) and Weiss et al. (2019) for a Poisson INAR(1), and by Weiss and Schweer (2015) for a PAR(1). A test based on the classical Pearson’s χ 2 statistic for the marginal distribution specified by the null hypothesis is proposed by Weiss (2018b) for a Poisson INAR(p).
A different approach to GOF testing for time series of counts is based on the probability generating function (PGF). Recall that if Y is a nonnegative discrete random variable then its PGF is defined as g Y ( u ) = E u Y for all u for which the expectation is finite, which is always the case for u [ 0 , 1 ] . The distribution of Y can be easily obtained from g Y ( k ) ( 0 ) and thus the PGF uniquely determines the distribution of Y. If Y = ( Y 1 , Y 2 ) is a count random vector then its PGF is defined as g Y ( u ) = E u 1 Y 1 u 2 Y 2 and exists for all u = ( u 1 , u 2 ) [ 0 , 1 ] 2 . Test criteria involving the estimates of the PGF of ( Y t , Y t 1 ) have proved as useful not only for assessing the specification of F in (1) or G in (4), but also for determination the model itself, see Meintanis and Karlis (2014) for INAR(1) and Schweer (2016) for a more general setup involving INAR(1) and PAR(1) as special cases. In both mentioned articles the considered test statistic is an integrated L 2 distance between the empirical PGF of ( Y t , Y t 1 ) and the parametric estimate of g ( Y t , Y t 1 ) under the null hypothesis. Hudecová et al. (2015) also consider tests based on PGF, but instead of the estimator g ( Y t , Y t 1 ) , their criteria are constructed as an integrated L 2 distance between the nonparametric estimate g ^ n ( u ) = n 1 t = 1 n u Y t of the marginal PGF g Y and a semiparametric estimate of g Y , which is model-specific. In the following sections we extend this approach from the univariate to the bivariate setup.

3. Bivariate Models for Time Series of Counts

The bivariate INAR and PAR models considered in this paper are based on the following bivariate distribution: We say that Y = ( Y 1 , Y 2 ) follows a bivariate Poisson distribution, see Kocherlakota and Kocherlakota (1992), denoted as BP ( λ 1 , λ 2 , ϕ ) , if its PGF is
g Y ( u ) = exp { λ 1 ( u 1 1 ) + λ 2 ( u 2 1 ) + ϕ ( u 1 1 ) ( u 2 1 ) } , u [ 0 , 1 ] 2 ,
where λ 1 ϕ > 0 , λ 2 ϕ > 0 , ϕ > 0 . This distribution arises via the trivariate reduction method, i.e., Y 1 = U 1 + U 2 , Y 2 = U 1 + U 3 , where U 1 , U 2 , U 3 are independent with Poisson distribution with mean ϕ , λ 1 ϕ , λ 2 ϕ respectively.
The models based on this kind of bivariate Poisson distribution seem to be the most popular in the literature despite some of their limitations. Please note that a different construction of a bivariate distribution with Poisson marginals is considered in Lakshminarayana et al. (1999).

3.1. Bivariate INAR Model

Suppose that we have a time series { Y t } composed by a pair of counts, i.e., Y t : = ( Y t , 1 , Y t , 2 ) . Following Latour (1997) we say that { Y t } follows a bivariate INAR model of the first order, in the following referred to as BINAR, if
Y t = A Y t 1 + ε t ,
where ε t = ( ε t , 1 , ε t , 2 ) are iid bivariate count random vectors with finite covariance matrix, A denotes a 2 × 2 matrix with elements a i j , i , j = 1 , 2 , and the operator ∘ is a multivariate generalization of (2). Namely the operator ∘ from (6) acts on a count random vector Y of dimension two by means of the equation
A Y = a 11 Y 1 + a 12 Y 2 a 21 Y 1 + a 22 Y 2 ,
where the univariate thinning operators was defined in (2).
Latour (1997) showed that if the spectral radius (the absolute value of the largest eigenvalue) of A is smaller than 1 and ε t has a finite covariance matrix then there exists a strictly stationary and ergodic process satisfying (6). Furthermore, the conditional least squares (CLS) estimate is shown to be consistent and asymptotically normal. Maximum likelihood estimation is considered in Pedeli and Karlis (2011) and Pedeli and Karlis (2013c), with a special emphasis on the Poisson specification for ε t . For estimation of BINAR with negative Binomial innovations the reader is referred to Mamode Khan et al. (2019) and references therein.

3.2. Bivariate PAR Model

Let F t = σ { Y s , s t } . We say that { Y t } follows a bivariate Poisson autoregression model of the first order, referred to as BPAR, model if,
Y t | F t 1 BP ( λ t , 1 , λ t , 2 , ϕ ) , λ t = λ t , 1 λ t , 2 = α 1 α 2 α + b 11 b 12 b 21 b 22 B Y t 1 , 1 Y t 1 , 2 ,
where α 1 , α 2 0 and the matrix B has non-negative entries and is of full-rank. This model is sometimes referred to as a bivariate Poisson INARCH(1) model, e.g., Lee et al. (2018). Liu (2012) proved that if 2 1 1 / p B p < 1 for some 1 p then { Y t , λ t } is strictly stationary and ergodic, with B p = sup x 0 , x C 2 B x p / x p . Strong consistency of the conditional maximum likelihood estimator (CMLE) of parameters was proved by Andreassen (2013) under the extra assumption ϕ < min { α 1 , α 2 } , while Lee et al. (2018) further showed that the CMLE has an asymptotic normal distribution.

4. Goodness-of-Fit Tests

Assume that we have observations Y 1 , , Y n , which come from a stationary bivariate time series of counts and we would like to perform a GOF test to a particular model for this data, with the model being fixed apart from finite-dimensional parameters. Let υ ( · , · ) be a non-negative weight function which will be further specified below. We propose as test statistic the weighted distance measure
S n , υ = n 0 1 0 1 [ g ^ n ( u 1 , u 2 ) g ^ n , 0 ( u 1 , u 2 ) ] 2 υ ( u 1 , u 2 ) d u 1 d u 2 ,
where g ^ n ( u 1 , u 2 ) is a non-parametric estimate of g Y ( u 1 , u 2 ) based on Y 1 , , Y n , given by g ^ n ( u 1 , u 2 ) = n 1 t = 1 n u 1 Y t , 1 u 2 Y t , 2 and where g ^ n , 0 ( u 1 , u 2 ) is a semi-parametric estimate of the PGF g Y ( u 1 , u 2 ) of Y , which is constructed specifically under the model being tested. The null hypothesis is rejected for large values of S n , υ .
By analogy to most, if not all, of the previously published work we consider u [ 0 , 1 ] 2 as our working interval, despite the fact that uniqueness of PGFs and corresponding consistency of the test might require working on a region of u containing the origin. Nevertheless we further investigate this aspect of our tests by simulations. See Esquível (2008) for a recent account of uniqueness of PGFs.
Although the idea of the proposed test statistics is analogous as for the univariate models mentioned in Section 2, the extension to a bivariate case is not straightforward, neither in terms of asymptotics nor on computational grounds as the multivariate nature of the data brings in some elevated technical difficulties. We state the necessary assumptions and provide a formal proof for the asymptotics of the suggested test statistic. In addition, extension of our GOF tests to more complex models is briefly discussed. In fact, and although here we only treat the bivariate case, extension to higher dimension also comes completely natural with our methods.

4.1. Tests for the BINAR Model

Let G = { g ε ( · ; θ ) : θ Θ } denote a parametric family of PGFs indexed by a parameter θ Θ , with Θ being a compact subset of R k , k N . We would like to test the null hypothesis
H 0 : { Y t } follows model ( 6 ) for some sequence ε t with PGF in G .
Assume that { Y t } is strictly stationary with PGF g Y . Then it follows from the properties of BINAR that
g Y ( u 1 , u 2 ) = E u 1 Y t , 1 u 2 Y t , 2 = E E u 1 Y t , 1 u 2 Y t , 2 | F t 1 = g ε ( u 1 , u 2 ; θ ) · g Y ( w 11 w 21 , w 12 w 22 ) ,
where
w i j = w i j ( u , a ) = 1 + a i j ( u i 1 ) ,
and a : = vec ( A ) = ( a 11 , a 21 , a 12 , a 22 ) denotes the vectorized version of A .
Here we used the fact that for a Binomial distribution with parameters ( m , p ) the PGF is ( 1 + p ( u 1 ) ) m , so that if additionally a ^ and θ ^ are suitable estimators of the unknown parameters, then (9) yields a semi-parametric estimate of g Y under model (6) as
g ^ n , 0 ( u 1 , u 2 ) = g ε ( u , θ ^ ) g ^ n ( w ^ 11 w ^ 21 , w ^ 12 w ^ 22 ) ,
where w ^ i j = w i j ( u , a ^ ) .
We will consider a GOF test for the null hypothesis H 0 with G being the family of bivariate Poisson distributions with unspecified parameters. Interest in testing H 0 lies in the fact that while this distribution is by far the most popular in the univariate as well as in the multivariate context, alternative specifications have also been employed such as models with innovations following a bivariate negative Binomial distribution; see for instance Pedeli and Karlis (2013a), Popović et al. (2018), and Kim and Lee (2017). Also BPAR models (see next section) may be considered to be alternative models of interest.

4.2. Tests for the BPAR Model

The corresponding null hypothesis for the BPAR model is formulated as
H ˜ 0 : { Y t } follows model ( 7 ) .
If { Y t } is strictly stationary, then it follows from the properties of the model that the PGF of Y t is given by
g Y ( u 1 , u 2 ) = E u 1 Y t , 1 u 2 Y t , 2 = E E ( u 1 Y t , 1 u 2 Y t , 2 | F t 1 ) = E exp { λ t , 1 ( u 1 1 ) + λ t , 2 ( u 2 1 ) + ϕ ( u 1 1 ) ( u 2 1 ) } = e ϕ ( u 1 1 ) ( u 2 1 ) + α 1 ( u 1 1 ) + α 2 ( u 2 1 ) × E e Y t 1 , 1 [ b 11 ( u 1 1 ) + b 12 ( u 2 1 ) ] e Y t 1 , 2 [ b 21 ( u 1 1 ) + b 22 ( u 2 1 ) ]
= e ϕ ( u 1 1 ) ( u 2 1 ) + α 1 ( u 1 1 ) + α 2 ( u 2 1 ) g Y ( w 1 , w 2 ) ,
where
w i = e b 1 i ( u 1 1 ) + b 2 i ( u 2 1 ) , i = 1 , 2 .
Thus, if we have an appropriate estimator θ ^ of the parameter θ : = ( α , vec ( B ) , ϕ ) = ( α 1 , α 2 , b 11 , b 21 , b 12 , b 22 , ϕ ) , then in view of (11) we may define a semi-parametric estimate of g Y as
g ^ n , 0 ( u 1 , u 2 ) = e ϕ ^ ( u 1 1 ) ( u 2 1 ) + α ^ 1 ( u 1 1 ) + α ^ 2 ( u 2 1 ) g ^ n ( w ^ 1 , w ^ 2 ) : = c ( u , θ ^ ) g ^ n ( w ^ 1 , w ^ 2 ) ,
where
w ^ i = e b ^ 1 i ( u 1 1 ) + b ^ 2 i ( u 2 1 ) ,
and
c ( u , θ ) = e ϕ ( u 1 1 ) ( u 2 1 ) + α 1 ( u 1 1 ) + α 2 ( u 2 1 ) .
Deviations from the null hypothesis H ˜ 0 include non-Poisson conditionals (see Heinen and Rengifo (2007)) as well as model violations towards more general specifications.
Remark 1.
We should remark that the tests proposed in this article are not for BINAR or BPAR models per se. Specifically the test criterion for BINAR is for bivariate counts for which the PGF satisfies Equation (9), and analogously the test criterion for BPAR is for bivariate counts for which the PGF satisfies Equation (12). In the context of time series of counts however, BINAR and BPAR models are framed by Equations (9) and (12), respectively, to such an extend that for all practical purposes these equations may be regarded as characterizing equations for the models themselves. Thus, our tests could be viewed as being on an equal footing with universally consistent methods such as those suggested by Fokianos and Neumann (2013), Jiménez-Gamero et al. (2020), and Leucht et al. (2015).

4.3. Computations

From (8) and by means of (10) and (14), we have after straightforward algebra
S n , υ = 1 n i = 1 n j = 1 n Δ i , j ,
where Δ i , j = Δ ( Y i , Y j ) with
Δ i , j = 0 1 0 1 D i ( u 1 , u 2 ) D j ( u 1 , u 2 ) υ ( u 1 , u 2 ) d u 1 d u 2
where
D i ( u 1 , u 2 ) = u 1 Y i , 1 u 2 Y i , 2 g ε ( u 1 , u 2 ; θ ^ ) ( w ^ 11 w ^ 21 ) Y i , 1 ( w ^ 12 w ^ 22 ) Y i , 2 ,
for the BINAR model, and
D i ( u 1 , u 2 ) = u 1 Y i , 1 u 2 Y i , 2 c ( u , θ ^ ) w ^ 1 Y i , 1 w ^ 2 Y i , 2 ,
for the BPAR model.
The integral figuring in Equation (16) with weight function υ ( u 1 , u 2 ) = u 1 γ u 2 γ , γ 0 , may be expressed in more elementary terms. However the final fully reduced explicit forms, which are available from the authors upon request, are not convenient from the computational point of view. Thus, the tests were implemented by numerical evaluation of the corresponding integrals.

5. Asymptotics

In this section, we study the asymptotic distribution of the test statistic S n , υ under the null hypothesis of a BINAR model as well as the corresponding limit null distribution under a BPAR model. Results for the behavior under fixed alternatives are also provided and show that in both cases the test is consistent against certain fixed alternatives. Finally resampling versions are proposed that circumvent the problem of unknown components in the aforementioned limit null distributions. Please note that under the standing assumptions these results are valid for general innovation distribution and weight function as well as for arbitrary parameter-estimates.

5.1. Asymptotics of the Test Statistic: BINAR Case

Assume the test statistic S n , υ (8) is constructed for testing H 0 with specified G = { g ε ( · ; θ ) : θ Θ } , for a compact Θ , i.e., the semi-parametric estimate g ^ n , 0 is constructed by (10). Consider the following assumptions:
(A.1)
Let υ ( · , · ) be a non-negative function satisfying 0 < 0 1 0 1 υ ( u 1 , u 2 ) d u 1 d u 2 < .
(A.2)
Let { Y t } follow model (6), where the spectral radius of A is smaller than 1.
(A.3)
Let G correspond to distributions with finite second moment. Furthermore, let the second partial derivatives of g ε with respect to θ exist and be continuous in θ and suppose that
sup u [ 0 , 1 ] 2 , θ Θ g ε θ i ( u , θ ) < , sup u [ 0 , 1 ] 2 , θ Θ u i g ε θ j ( u , θ ) < , [ 0 , 1 ] 2 2 g ε θ i θ j ( u , θ ) υ ( u ) d u <
hold true for all θ Θ and i , j { 1 , , k } , where (recall) that k is the dimension of θ .
(A.4)
Let ζ ^ : = ( θ ^ , a ^ ) be an estimator of ζ : = ( θ , a ) such that for some q 1
n ( ζ ^ ζ ) = 1 n t = q + 1 n l ( Y t , Y t 1 , , Y t q , ζ ) + o P ( 1 ) ,
where l t = l ( Y t , Y t 1 , , Y t q , ζ ) = ( l t , 1 , l t , 2 ) form a strictly stationary and ergodic sequence of martingale differences with finite variance. Here, l t , 1 is of the same dimension as θ and l t , 2 has the dimension of vec ( A ) , i.e., four.
Under (A.2) and (A.3), { Y t } is strictly stationary and ergodic with finite second order moments, see Franke and Rao (1995). Regarding possible estimators ζ ^ in (A.4), Franke and Rao (1995) considered the CMLE and proved its asymptotic normality under a set of regularity conditions. These conditions involve the finiteness of E ε t 3 and some further assumptions on the distribution of ε t .
Theorem 1.
Under (A.1)–(A.4) and as n , the limit distribution of S n , υ under the null hypothesis H 0 is the same as the distribution of
[ 0 , 1 ] 2 Q 2 ( u , ζ ) υ ( u ) d u ,
where { Q ( u , ζ ) , u [ 0 , 1 ] 2 } is a Gaussian random field with zero mean and covariance structure
E ( Q ( u , ζ ) Q ( v , ζ ) ) = E ( [ u 1 Y q + 1 , 1 u 2 Y q + 1 , 2 g ε ( u , θ ) W q ( u , a ) l q + 1 , 1 ( ζ ) h 2 ( u , ζ ) l q + 1 , 2 ( ζ ) h 3 ( u , ζ ) ] × [ v 1 Y q + 1 , 1 v 2 Y q + 1 , 2 g ε ( v , θ ) W q ( v , a ) l q + 1 , 1 ( ζ ) h 2 ( v , ζ ) l q + 1 , 2 ( ζ ) h 3 ( v , ζ ) ] ) ,
where l q + 1 , 1 and l q + 1 , 2 are from assumption (A.4) and
W t ( u , a ) = i , j = 1 2 w i j ( u , a ) Y t , j , h 2 ( u , ζ ) = θ g ε ( u , θ ) g Y ( w 11 w 21 , w 12 w 22 ) , h 3 ( u , ζ ) = g ε ( u , θ ) E a W 1 ( u , a ) .
The proof of the assertion is postponed to the Appendix A.
The asymptotic distribution of the test statistic S n , υ depends on several unknown quantities. One possibility is to generate the Gaussian random field figuring in Theorem 1 with the theoretical quantities replaced by some consistent estimators and then compute the critical values. Another possibility is the parametric bootstrap, which is quite natural here. The justification of the bootstrap approximation under the null hypothesis proceeds along similar lines as the proof of Theorem 1, conducted conditionally on the observed data Y 1 , , Y n and with the help of assertions for triangular arrays and sums of martingale difference arrays.
Write S n : = S n ( Y 1 , , Y n ; ζ ^ ) for the test statistic based on the original observations Y 1 , , Y n , and the resulting parameter estimate ζ ^ . (Here for simplicity we suppress the dependence of S n , υ on the weight function υ ( · , · ) ).
  • Generate ε t , t = 1 , n , where ε t are iid and follow the distribution with PGF g ε ( u , θ ^ ) .
  • Compute pseudo-observations Y t , t = 1 , n , using Equation (6) with ε t replaced by ε t , t = 1 , n , and A replaced by A ^ .
  • Fit the model (6) using Y 1 , , Y n , and compute the bootstrap estimator ζ ^ of ζ .
  • Compute the corresponding test statistic S n : = S n ( Y 1 , , Y n ; ζ ^ ) .
  • Repeat steps 1–4 several times, say B, and obtain the sequence of test statistics, S 1 , n , , S B , n .
  • Compute the p-value as ( # b : S b , n S n ) / ( B + 1 ) .
Next we shortly discuss the limit behavior of the test statistic under alternatives of type g ε G . We assume that model (6) holds true but g ε in the null hypothesis H 0 does not belong to G . Moreover suppose that the estimators ( θ ^ , a ^ ) satisfy
( θ ^ , a ^ ) 𝖯 ( θ A , a A ) ,
for some θ A Θ and some a A such that the largest eigenvalue of the respective matrix A A is in absolute value smaller then 1.
Theorem 2.
Let { Y t } be strictly stationary with PGF g A ( u 1 , u 2 ) that is continuous in ( u 1 , u 2 ) [ 0 , 1 ] 2 and let g ε ( u , θ ) be continuous in θ. Let (A.1) and (17) be satisfied. Then, as n ,
S n , v n 𝖯 [ 0 , 1 ] 2 g A ( u ) g ε ( u , θ A ) g A ( w 11 , A w 21 , A , w 12 , A w 22 , A ) 2 v ( u ) d u ,
where w i j , A = w i j ( u , a A ) , i , j = 1 , 2 .
The proof is omitted since it suffices to follow the line of proof of Theorem 1 and use stationarity and ergodicity of { Y t } .
The right-hand side of (18) is strictly positive unless the true PGF g ε ( · ) coincides with the PGF g ε ( · ; θ a ) from the null hypothesis H 0 . This fact and the uniqueness of the PGF implies the consistency of the test which rejects the null hypothesis H 0 for large values of the test statistic S n , v under such fixed alternatives. The test is also consistent for other types of fixed alternatives, e.g., against model violation. This feature of the test is further illustrated by Monte Carlo simulations in Section 7.
Please note that the test even has (non-negligible) power for some local alternatives, i.e., when the difference g ε ( u ) g ε ( u ; θ a ) tends to 0 not too fast as n and g ε ( u ) depends on n. However, a rigorous proof of this result is quite technical, and therefore it is not discussed here any further.

5.2. Asymptotics of the Test Statistic: BPAR Case

This section considers the problem of testing the BPAR model. Recall that θ = ( α 1 , α 2 , b 11 , , b 22 , ϕ ) stands for the vector of the parameters of the model (7), and suppose that:
(B.1)
The series { Y t } is strictly stationary solution of (7) with parameters θ Θ such that ϕ < min { α 1 , α 2 } and Θ is compact.
(B.2)
The estimator θ ^ of the parameter θ is such that
n ( θ ^ θ ) = 1 n t = q + 1 n l ( Y t , Y t 1 , , Y t q , θ ) + o P ( 1 ) ,
where l t = l ( Y t , Y t 1 , , Y t q , θ ) = ( l t , 1 , , l t , 7 ) form a strictly stationary and ergodic sequence of martingale differences with finite variances.
An estimator θ ^ that satisfies (B.2) is for instance the CMLE; see Lee et al. (2018).
Theorem 3.
Under (A.1),(B.1)–(B.2) the limit distribution of S n , υ as n is the same as the distribution of
[ 0 , 1 ] 2 Q 2 ( u , θ ) υ ( u ) d u ,
where { Q ( u , θ ) , u [ 0 , 1 ] 2 } is a Gaussian random field with zero mean and covariance structure
E Q ( u , θ ) Q ( v , θ ) = E ( [ u 1 Y t , 1 u 2 Y t , 2 c ( u , θ ) w 1 ( u ) Y t , 1 w 2 ( u ) Y t , 2 + l t d ( θ , u ) ] × [ v 1 Y t , 1 v 2 Y t , 2 c ( v , θ ) w 1 ( v ) Y t , 1 w 2 ( v ) Y t , 2 + l t d ( θ , v ) ] ) ,
with l t defined in assumption (B.2), c ( u , θ ) defined in (15),
w i ( u ) = e β i 1 ( u 1 1 ) + β i 2 ( u 2 1 ) , i = 1 , 2 ,
d ( θ , u ) = c ( u , θ ) E D 1 ( u ) w 1 Y 1 , 1 w 2 Y 1 , 2
and
D t ( u ) = ( u 1 1 , u 2 1 , ( u 1 1 ) Y t , 1 , ( u 1 1 ) Y t , 2 , ( u 2 1 ) Y t , 1 , ( u 2 1 ) Y t , 2 , ( u 1 1 ) ( u 2 1 ) ) .
Similarly as for the BINAR model, the parametric bootstrap can be carried out in a very natural way and is recommended for practical use. The justification of the bootstrap approximation would again proceed along similar lines as the proof of Theorem 3.
The following theorem is analogous to Theorem 2 and describes the behavior of the test statistic under some alternatives.
Theorem 4.
Let { Y t } be strictly stationary with PGF, g ˜ A ( u 1 , u 2 ) that is continuous in ( u 1 , u 2 ) [ 0 , 1 ] 2 and let the estimator θ ^ of θ = ( α , vec ( B ) , ϕ ) satisfy θ ^ P θ A for some θ A Θ . Let (A.1) be satisfied. Then, as n ,
S n , v n [ 0 , 1 ] 2 { g ˜ A ( u ) c ( u , θ A ) g ˜ A ( w 1 , A , w 2 , A ) 2 υ ( u ) d u ,
where w 1 , A , w 2 , A are defined by (13) with b i j replaced by b i j , A , and c ( u , θ ) is defined in (15).
The proof is omitted since it suffices to follow the line of the proof of Theorem 3 and to use stationarity and ergodicity of { Y t } . Possible comments for this theorem are quite parallel to Theorem 1 and therefore are omitted.

6. Extension to Generalized BPAR Model

Extension of the proposed test to arbitrary higher order or dimension will be discussed in Section 9. In this section we aim at a new PGF-based GOF test for a relatively mild, yet very important, generalization of the BPAR model. Specifically consider the following model for { Y t } :
Y t | F t 1 BP ( λ t , 1 , λ t , 2 , ϕ ) , λ t = α + A λ t 1 + B Y t 1 ,
where α = ( α 1 , α 2 ) with α 1 , α 2 0 , A , B are matrices of non-negative entries and B is of full rank. (Please note that if A is replaced by the zero matrix, then model (19) reduces to the BPAR model defined in (7)). Model (19) was studied, e.g., in Andreassen (2013), Liu (2012), and Lee et al. (2018), with the acronym INGARCH(1,1). If A p + 2 1 1 / p B p < 1 for some 1 p , then { Y t , λ t } is strictly stationary and ergodic, see Liu (2012), and if ϕ < min { a 1 , a 2 } , where ( a 1 , a 2 ) = ( I A ) 1 α , then the CMLE is consistent and asymptotically normal, see Andreassen (2013), Lee et al. (2018).
Assume that { Y t } follows a stationary model (19). Following analogous arguments as in Section 4.2 we have that the PGF of Y under this model is given by,
g Y ( u 1 , u 2 ) = E u 1 Y t , 1 u 2 Y t , 2 = E E u 1 Y t , 1 u 2 Y t , 2 | F t 1 = E exp { λ t , 1 ( u 1 1 ) + λ t , 2 ( u 2 1 ) + ϕ ( u 1 1 ) ( u 2 1 ) } = e α 1 ( u 1 1 ) + α 2 ( u 2 1 ) + ϕ ( u 1 1 ) ( u 2 1 ) E w 1 , b Y t 1 , 1 w 2 , b Y t 1 , 2 w 1 , a λ t 1 , 1 w 2 , a λ t 1 , 2 = e α 1 ( u 1 1 ) + α 2 ( u 2 1 ) + ϕ ( u 1 1 ) ( u 2 1 ) h Y , λ ( w 1 , b , w 2 , b , w 1 , a , w 2 , a ) ,
where
w i , a = e a 1 i ( u 1 1 ) + a 2 i ( u 2 1 ) , w i , b = e b 1 i ( u 1 1 ) + b 2 i ( u 2 1 ) , i = 1 , 2 ,
with the joint PGF of the vector ( Y , λ ) defined by
h Y , λ ( u 1 , u 2 , v 1 , v 2 ) : = E u 1 Y t , 1 u 2 Y t , 2 v 1 λ t , 1 v 2 λ t , 2 .
It is straightforward to device a test for the model (19) by using Equation (20) and proceeding analogously as with Equation (8). However the asymptotics of such a test as well as its actual implementation require a separate investigation. In this connection preliminary Monte Carlo results showed some promise but there were also problems, and therefore we decide not to pursue this extension any further here.

7. Simulations

The finite sample behavior of the suggested bootstrap test is investigated in the following simulation study. We consider the null hypotheses H 0 of Poisson BINAR and H 0 ˜ of BPAR model and investigate the size of the test under the null hypothesis and the power for various alternatives.
The unknown parameters of the BINAR model are estimated using the CLS method, and ϕ is estimated using the moment method, see Pedeli and Karlis (2013b). The parameters of BPAR model are estimated by the CMLE method. For simplicity the weight function is set to υ ( · , · ) 1 , i.e., we take γ = 0 in υ ( u 1 , u 2 ) = u 1 γ u 2 γ . The simulations were conducted in the R-computing environment R Core Team (2019) and by employing the warp-speed bootstrap of Giacomini et al. (2013) for M = 1000 repetitions. When using this method, B=1 bootstrap samples are generated for each Monte Carlo repetition and the resulting p-value is computed from the overall bootstrap sample of M replicas.
Our results are for sample size n = 100 , 200 and n = 500 at level of significance α = 0.01 , 0.05 and 0.1 for both H 0 and H 0 ˜ . For the BINAR model, a reasonable alternative might be BPAR model or model (6) with innovations ε t following a distribution other than the bivariate Poisson. Such a popular alternative is a bivariate distribution with negative Binomial marginals. There are several possibilities as to how to generate such variables. Here we consider the bivariate negative Binomial distribution BNB ( λ 1 , λ 2 , r ) of Dunn (1967), whereby ( U 1 , U 2 ) BNB ( λ 1 , λ 2 , r ) , with U i , i = 1 , 2 , being marginally negative Binomial with mean λ i and variance λ i ( 1 + λ i / r ) and cov ( U 1 , U 2 ) = λ 1 λ 2 / r . This bivariate negative Binomial distribution is also used in an alternative considered to a BPAR model. Namely a model of form (7) with conditional distribution BNB ( λ t 1 , λ t 2 , r ) instead of bivariate Poisson is considered, with the dependence of λ t i on F t 1 the same as specified in Equation (7). We will refer to this model as the negative Binomial BINARCH. Finally, we also explored the power of the test when testing H 0 ˜ and the data follow BINAR model and vice versa.
The results for the size of the two tests are summarized in Table 1 and for the BINAR model were obtained by using either A = A 1 or A = A 2 , where
A 1 = 0.3 0.1 0.2 0.2 , A 2 = 0.6 0.0 0.0 0.7 ,
and ε t following a BP ( 6 , 4 , 3 ) distribution. For BPAR, we set α = ( 6 , 4 ) , ϕ = 3 , and again B equal either to A 1 or A 2 . As it may be seen from Table 1, the tests are conservative if the matrix A 1 is used in the data generating process. On the other hand, the size slightly exceeds the nominal size α if A 2 is used in the data generating process Please note that a similar phenomenon was already observed for the univariate series in Hudecová et al. (2015): The test may be slightly conservative for certain parametric settings of the model at hand, while it can be slightly anticonservative for other settings. Recall also that the unknown parameters are estimated by a CLS method for BINAR models whereas CMLE is used for BPAR models, which might also partly explain the different behaviour observed in Table 1. In either case however, the observed size approaches the prescribed significance level as the sample size increases.
In this connection, the rather poor small sample behavior observed in the left part of Table 1 for the BINAR model parameterised by the matrix A 1 , can be improved by considering a modified test statistic for which integration in Equation (8) is carried over [ 1 , 1 ] 2 rather than over [ 0 , 1 ] 2 ; see the corresponding discussion at the end of Section 4. Specifically this approach (performed on the same simulated data) and for nominal size α = 0.05 leads to empirical size equal to 0.033, 0.039 and 0.052, for n = 100 , 200 , and 500, respectively. For all other settings however, the results for the modified test statistic (available from the authors upon a request) are very similar and hence we only present here results for the original test statistic S n , υ defined by Equation (8).
The power of the test for H 0 with data from a BINAR model with negative Binomial innovations is provided in Table 2. The considered bivariate negative Binomial distribution is BNB ( 6 , 4 , r ) for r = 5 , 10 . Please note that for larger r, the BNB distribution is closer to a bivariate Poisson distribution, and this fact is also reflected in the power of the test, which is lower for r = 10 . For example, a sample size n = 100 seems to be insufficient for distinguishing between Poisson and negative binomial INAR for this larger r. However, as the sample size n growths, the test performs very well even for r = 10 irrespective of the matrix A being used to simulate the model.
A similar observation also holds for the results in Table 3 which correspond to the null hypothesis H 0 ˜ with data from a negative Binomial BINARCH. These results were obtained with the same values of α and B used for the null hypothesis, and r = 30 , 50 .
On the other hand, Table 4 reports the power of the test for H 0 with data following a BPAR model and the power of the test for H ˜ 0 with data from a Poisson BINAR model. The results in Table 4 show that if the matrix B in the BPAR model equals A 1 then the test for H 0 fails to distinguish between the two models for sample sizes up to n = 500 . In contrast, if the matrix B is equal to A 2 in BPAR and we test for Poisson BINAR then the power is satisfactory even for n = 200 . The same observation holds for the opposite situation when one tests for BPAR and the data come from a Poisson BINAR.

8. Real-Data Application

Joint modelling of count observations finds important applications in the insurance industry, see for instance Partrat (1994), and Vernic (1997). In this connection, it is a common practice for insurance companies to split the reported claims into several types. Typically, it is reasonable to expect that aggregate amounts (daily or monthly) of these different types of claims to be dependent, see, e.g., Shi et al. (2016). If the mean size of claims counts is high, then classical models for continuous variables could be applied. On the other hand, if the observed counts of claims are formed by small integers it is appropriate to treat the data as genuine counts, and, consequently, engage models and methods specifically tailored for count time series.
We illustrate this kind of application on real data sets on the monthly number of claims of short-term disability benefits made by injured workers to the British Columbia Workers Compensation Board. The time period is from January 1985 to December 1994. The original data set from Freeland (1998) contains five time series corresponding to five different injury categories: burn injuries, soft tissue injuries, cuts, dermatitis and dislocations. These five time series have been previously analyzed by several authors, and separate univariate models were fitted. It has been found that the Poisson INAR is appropriate for all five series, except for series # 3 (cuts), for which this model is not appropriate, see e.g., Freeland (1998), Zhu and Joe (2006), Hudecová et al. (2015). In particular, Freeland and McCabe (2004) and Zhu and Joe (2006) suggest to model the time series of cuts claims using an extension of INAR(1) model with a seasonal component. On the other hand, Biswas and Song (2009) argue that the ACF and PACF do not indicate a significant seasonal effect.
We first consider a bivariate series of soft tissue injuries claims and dermatitis claims; see Figure 1 (left panel). Possible dependence among these two series may be due to the fact that a major accident causes several injuries, often of different types. Previous analyses in Freeland (1998) or Hudecová et al. (2015) reveal that the marginal INAR(1) models might be appropriate for these two series. Hence, we consider a Poisson BINAR model with a diagonal matrix A and estimate the parameters using the conditional least squares method (note that if a general, i.e., non-diagonal, matrix A is considered, then the estimates of the off-diagonal entries are very close to zero). Estimators of parameters were obtained by the methods used in the simulations of the previous section. The matrix A is estimated as A ^ = diag ( 0.452 , 0.293 ) , and the parameters of the bivariate Poisson distribution of the innovations ε t are estimated as λ ^ = ( 5.376 , 0.202 ) , ϕ ^ = 0.109 . The resulting GOF test applied with υ 1 and B = 999 bootstrap samples leads to S n , υ = 8.618 · 10 4 with p-value 0.327 . Thus, one may conclude that the Poisson BINAR model with a diagonal matrix A seems to fit the data well.
On the other hand, if one considers a bivariate series of soft tissue injuries claims and cuts claims, see Figure 1 (right panel), and tests the GOF for a Poisson BINAR model with a diagonal matrix A (which corresponds to univariate INAR models for the two series), the null hypothesis is rejected with p-value 0.0370 . A Poisson BINAR with a general matrix A is rejected as well (p-value 0.002 ). This is in accordance with the findings of previous papers mentioned above. However, the hypothesis of BPAR is not rejected with p-value 0.277 . The MLEs in this model are α ^ = ( 5.144 , 2.944 ) , ϕ ^ = 0.721 and B ^ is such that vec ( B ^ ) = ( 0.455 , 0.004 , 0.043 , 0.521 ) . One could further scrutinize this data set by postulating an extended (nonstationary) BPAR model with a seasonal component. This might improve the models’ predictions compared to the fitted stationary model, but such a model lies outside the scope of this paper.

9. Concluding Remarks

We suggest consistent goodness-of-fit tests for bivariate INAR and bivariate Poisson autoregressive models, estimated by least squares and maximum likelihood, respectively, with well defined limit null distributions. Since these limit distributions are complicated we suggest parametric bootstrap resampling which engages distributional assumptions featuring in the null hypothesis in order to actually carry out the tests. Monte Carlo results show that this bootstrap version of the new tests is generally reasonably sized and has good power against certain popular alternative configurations. Our real-data applications are in the direction of further scrutiny and better understanding of the mechanisms generating the data at hand.
At this point, we wish to discuss potential extension of the tests to models of higher order or dimension, such as the multivariate INAR type processes studied by Franke and Rao (1995), Latour (1997), Pedeli and Karlis (2011), Pedeli and Karlis (2013c), and Pedeli and Karlis (2013a), and corresponding extensions of PAR type models considered by Liu (2012), Andreassen (2013), Lee et al. (2018), Ciu and Zhu (2018) and Ciu et al. (2020). In this connection, and while it is conceptually straightforward to extend the tests for INAR or PAR models of higher order, see also Section 6 of Hudecová et al. (2015), we wish to emphasize the fact that a certain order and/or dimension increase brings about serious challenges in estimation as well as on the actual implementation of the methods due to the potentially great number of new parameters introduced.
Finally, before we close, we wish to point out that while this paper deals solely with stationary time series models, certain real world phenomena, including time series of counts, often exhibit deterministic trends or seasonal patterns, and thus GOF tests for non-stationary models would be of great practical interest. However, despite the fact that univariate INAR and PAR models with deterministic components have been extensively studied in the literature, multivariate extensions of such nonstationary models are rather scarce; we refer to Santos et al. (2019) for some recent works.

Author Contributions

Conceptualization, Š.H., M.H., S.G.M.; methodology, Š.H., M.H., S.G.M.; software Š.H.; writing—original draft preparation, Š.H., S.G.M.; writing—review and editing, Š.H., M.H., S.G.M. All authors have read and agreed to the published version of the manuscript.

Funding

The work of Šárka Hudecová and Marie Hušková is funded by the Czech Science Foundation project GA18-08888S.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in Freeland (1998) or on request from the corresponding author.

Acknowledgments

The authors are grateful to the Associated Editor and the anonymous reviewers for their careful inspection of the paper and valuable comments, which led to an improved manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Proofs

Appendix A.1. Proof of Theorem 1

Proof. 
Denote as A [ 0 , 1 ] 4 such that if a A then A has the largest eigenvalue smaller than 1. The test statistic S n , υ is of the form
S n , υ = [ 0 , 1 ] 2 1 n t = 1 n u 1 Y t , 1 u 2 Y t , 2 g ε ( u , θ ^ ) W t ( u , a ^ ) 2 υ ( u ) d u = [ 0 , 1 ] 2 1 n t = 2 n u 1 Y t , 1 u 2 Y t , 2 g ε ( u , θ ^ ) W t 1 ( u , a ^ ) 2 υ ( u ) d u + o P ( 1 ) .
Let A t = A t ( θ , a ) = g ε ( u , θ ) W t ( u , a ) . In the following, we omit for a moment the argument u in all the considered functions. Namely we consider W t = W t ( a ) as a function of the argument a alone, i.e., W t is a differential of W t with respect to a , W t = a W t = ( W / a 11 , W / a 12 , W / a 21 , W / a 22 ) . Likewise, g ε = θ g ε is a differential of g ε with respect to θ . Then a Taylor expansion gives:
A t ( θ ^ , a ^ ) = A t ( θ , a ) + θ ^ θ a ^ a W t ( a ) g ε ( θ ) g ε ( θ ) W t ( a ) + R t ,
where
R t = 1 2 θ ^ θ a ^ a W t ( a ) H g ε ( θ ) g ε ( θ ) W t ( a ) W t ( a ) g ε ( θ ) g ε ( θ ) H W t ( a ) θ ^ θ a ^ a ,
with a = a ^ + c ( a a ^ ) and θ = θ ^ + d ( θ θ ^ ) for some c , d ( 0 , 1 ) and where H g ε = 2 g ε / ( θ k θ l ) k , l is the k × k matrix of second partial derivatives of g ε with respect to θ . Similarly, H W t is the 4 × 4 matrix of second partial derivatives of W t with respect to a . For W t we have
W t a i j = W t w i j Y t , j Y t 1 , j w i j Y t , j 1 ( u i 1 ) = W t Y t , j w i j ( u i 1 ) ,
and clearly W t / a i j Y t j for all u [ 0 , 1 ] 2 and a A . Furthermore,
2 W t a i j a k l = W t Y t 1 , j Y t 1 , l w i j w k l ( u i 1 ) ( u k 1 ) , ( i , j ) ( k , l ) , W t Y t 1 , j ( Y t 1 , j 1 ) w i j 2 ( u i 1 ) 2 ( i , j ) = ( k , l ) ,
and thus, 2 W t / ( a i j a k l ) Y t 1 , j Y t 1 , l for all u [ 0 , 1 ] 2 and all a A .
Due to the finite second order moments of { Y t } and Assumption (A.2) we have
1 n t = 2 n R t K 1 n θ ^ θ 2 + n a ^ a 2 K 2 + K 3 n ( θ ^ θ ) 1 k × 4 ( a ^ a ) ,
where K i are constants and 1 k × 4 is a matrix of 1’s of dimension k × 4 . Hence, it follows from the Cauchy–Schwartz inequality that [ 0 , 1 ] 2 1 n t = 2 n R t 2 υ ( u ) d u 𝖯 0 , as n , and thus the asymptotic distribution of S n , υ is the same as the asymptotic distribution of
[ 0 , 1 ] 2 J 1 n ( θ , a , u ) + J 2 n ( θ , a , u ) + J 3 n ( θ , a , u ) 2 υ ( u ) d u ,
where
J 1 n ( θ , a , u ) = 1 n t = 2 n u 1 Y t , 1 u 2 Y t , 2 g ε ( u , ε ) W t 1 ( u , a ) , J 2 n ( θ , a , u ) = n ( θ ^ θ ) g ε ( u , θ ) 1 n t = 2 n W t 1 ( u , a ) , J 3 n ( θ , a , u ) = g ε ( u , θ ) n ( a ^ a ) 1 n t = 2 n a W t 1 ( u , a ) .
Regarding the behavior of J 2 n , notice that | W t ( a ) | 1 and thus, for a given a and as n
1 n t = 2 n W t 1 ( a ) E W 1 ( u , a ) = E i , j = 1 2 w i j Y 1 , j = g Y ( w 11 w 21 , w 12 w 22 ) ,
uniformly in u [ 0 , 1 ] 2 due to the uniform ergodicity theorem. Now define
J 2 n 0 ( θ , a , u ) = 1 n t = q + 1 n l t 1 h 2 ( θ , a , u ) ,
where h 2 ( θ , a , u ) = θ g ε ( θ , u ) g Y ( w 11 w 21 , w 12 w 22 ) . Due to (A.2), it follows that [ 0 , 1 ] 2 [ J 2 n ( θ , a , u ) J 2 n 0 ( θ , a , u ) ] 2 υ ( u ) d u = o P ( 1 ) as n . Similarly, in view of the form of W t 1 ( a ) given above and the finiteness of E Y 1 j , we get from the uniform ergodicity theorem that as n
g ε ( θ , a ) 1 n t = 2 n W t 1 ( u , a ) 𝖯 g ε ( u , θ ) E W 1 ( u , a ) = : h 3 ( θ , a , u )
uniformly in u [ 0 , 1 ] 2 . Define
J 3 n 0 ( θ , a , u ) = g ε ( θ , u ) 1 n t = q + 1 n l t 2 h 3 ( θ , a , u ) .
As | h 3 ( θ , a , u ) | max j = 1 , 2 E Y 1 , j then [ 0 , 1 ] 2 [ J 3 n ( θ , a , u ) J 3 n 0 ( θ , a , u ) ] 2 υ ( u ) d u = o P ( 1 ) . Hence, it suffices to study the asymptotic behavior of [ 0 , 1 ] 2 Q n 2 ( u , θ , a ) υ ( u ) d u , where
Q n ( u , θ , a ) = J 1 n 0 ( θ , a , u ) + J 2 n 0 ( θ , a , u ) + J 3 n 0 ( θ , a , u )
with J 1 n 0 ( θ , a , u ) = 1 n t = q + 1 n u 1 Y t , 1 u 2 Y t , 2 g ε ( u , ε ) W t 1 ( u , a ) . We will make use Theorem 22 from Ibragimov and Chasminskij (1981) and the subsequent remark in order to show that the integral [ 0 , 1 ] 2 Q n 2 ( u , θ , a ) υ ( u ) d u converges in distribution to [ 0 , 1 ] 2 Q 2 ( u , θ , a ) υ ( u ) d u , where Q is a Gaussian random field specified in the statement of the theorem. To this end, we need to verify that
  • sup n , u , θ , a E Q n 2 ( u , θ , a ) < .
  • There exist constants α > 0 and H > 0 such that
    sup n , θ , a E | Q n ( u , θ , a ) Q n ( v , θ , a ) | 2 H u v α .
  • The marginal distributions of Q n converge to the marginal distributions of Q uniformly in θ and a .
First, notice that Q n is of a form 1 n t = q + 1 n L t , where { L t } is a martingale difference sequence L t = u 1 Y t 1 u 2 Y t 2 g ε ( θ , u ) W t 1 ( u , a ) l t 1 ( θ , a ) h 2 ( θ , a , u ) l t 2 ( θ , a ) h 3 ( θ , a , u ) . Thus E Q n 2 1 n t = q + 1 n E | L t | 2 4 ( 1 + g ε 2 + h 2 2 E l 11 2 + h 3 2 E l 12 2 ) . Since we assume finite variances of l t , j and thus h j 2 are bounded (due to finiteness of E Y t , j and assumption A.2), condition I. directly follows. Condition III. follows from the central limit theorem for martingale difference sequences on L t .
In order to prove condition II., we will show that E | J j n 0 ( θ , a , u ) J j n 0 ( θ , a , v ) | 2 D u v α for j = 1 , 2 , 3 , for some D > 0 and α > 0 . In the following K is a generic constant. For j = 1 we have
E | J 1 n 0 ( θ , a , u ) J 1 n 0 ( θ , a , v ) | 2 = 1 n E t = q + 1 n u 1 Y t , 1 u 2 Y t , 2 v 1 Y t , 1 v 2 Y t , 2 g ε ( u ) W t 1 ( u ) + g ε ( v , θ ) W t 1 ( v , a ) 2 K { E [ u 1 Y t , 1 u 2 Y t , 2 v 1 Y t , 1 v 2 Y t , 2 ] 2 + g ε ( u , θ ) 2 E [ W 1 ( u , a ) W 1 ( v , a ) ] 2 + ( g ε ( u , θ ) g ε ( v , θ ) ) 2 E | W 1 ( u , a ) | 2 } .
By the mean value theorem applied on the function f ( x , y ) = x a y b and due to the fact that E | Y t , j | is finite, we get that E [ u 1 Y t , 1 u 2 Y t , 2 v 1 Y t , 1 v 2 Y t , 2 ] 2 K u v 2 . Next, | g ε | 1 and the partial derivatives | W 1 / u j | are bounded by | Y 1 , 1 + Y 1 , 2 | , which has finite expectation. This implies that g ε ( u , θ ) 2 E [ W 1 ( u , a ) W 1 ( v , a ) ] 2 K u v 2 . Finally, it follows from the definition of g ε as a PGF that it is Lipschitz. Furthermore, E | W 1 ( u , a ) | 2 1 , which implies ( g ε ( u , θ ) g ε ( v , θ ) ) 2 E | W 1 ( u , a ) | 2 K u v 2 . Together, we get E | J 1 n 0 ( θ , a , u ) J 1 n 0 ( θ , a , v ) | 2 K u v 2 . Similar arguments are used to show that the condition holds also for j = 3 . If we use the assumption that l t has finite variances, it remains to show that h j ( θ , a , u ) h j ( θ , a , v ) 2 K u v α . First, notice that the partial derivatives of h 3 with respect to u are continuous functions for u [ 0 , 1 ] 2 . Thus, they are bounded and we can apply the mean value theorem on the components of h 3 , which implies that h 3 ( θ , a , u ) h 3 ( θ , a , v ) 2 K u v 2 . Similarly, if the assumption (A.2) holds then the partial derivatives of h 2 are bounded and this implies that h 2 ( θ , a , u ) h 2 ( θ , a , v ) 2 K u v 2 , which completes the proof of condition II., and thus the assertion of the theorem follows. □

Appendix A.2. Proof of Theorem 3

Proof. 
The proof proceeds along the same lines as the proof of Theorem 1. A Taylor expansion is used to show that
S n , υ = 0 1 0 1 [ J 0 ( θ , u ) J 1 ( θ , u ) ] 2 υ ( u ) d u + o P ( 1 ) ,
where D t ( u ) are defined in the statement of the theorem,
J 0 ( θ , u ) = 1 n t = 2 n u 1 Y t , 1 u 2 Y t , 2 f t 1 ( θ , u ) , J 1 ( θ , u ) = n 1 / 2 t = 2 n f t 1 ( θ , u ) D t 1 ( u ) ( θ ^ θ ) , f t 1 ( θ , u ) = c ( u , θ ) w 1 ( u ) Y t 1 , 1 w 2 ( u ) Y t 1 , 2 .
The term J 1 ( θ , u ) can be replaced by J 10 ( θ , u ) = E D 1 ( u ) f 1 ( θ , u ) n 1 / 2 t = q n l t . Finally, the conditions I.–III. from Theorem 22 of Ibragimov and Chasminskij (1981) can be verified similarly as it is done for BINAR using the fact that the fourth moments of Y t are finite and the parametric space Θ is compact. □

References

  1. Ahmad, Ali, and Christian Francq. 2016. Poisson QMLE of count time series models. Journal of Time Series Analysis 37: 291–314. [Google Scholar] [CrossRef]
  2. Al-Osh, Mohamed A., and Abdulhamid A. Alzaid. 1987. First–order integer–valued autoregressive (INAR(1)) process. Journal of Time Series Analysis 8: 261–75. [Google Scholar] [CrossRef]
  3. Aleksandrov, Boris, and Christian H. Weiss. 2020. Testing the dispersion structure of count time series using Pearson residuals. AStA Advances in Statistical Analysis 104: 325–61. [Google Scholar] [CrossRef]
  4. Andreassen, Camilla Mondrup. 2013. Models and Inference for Correlated Count Data. Ph.D. thesis, Aarhus University, Aarhus, Denmark. [Google Scholar]
  5. Biswas, Atanu, and Peter Xue Kun Song. 2009. Discrete-valued ARMA processes. Statistics and Probability Letters 79: 1884–89. [Google Scholar] [CrossRef]
  6. Ciu, Yan, Qi Li, and Fukang Zhu. 2020. Flexible bivariate Poisson integer-valued GARCH model. Annals of the Institute of Statistical Mathematics 72: 1449–77. [Google Scholar]
  7. Ciu, Yan, and Fukang Zhu. 2018. A new bivariate integer-valued GARCH model allowing for negative cross-correlation. TEST 27: 428–52. [Google Scholar]
  8. Davis, Richard A., Scott H. Holan, Robert Lund, and Nalini Ravishanker. 2015. Handbook of Discrete-Valued Time Series. New York: Chapman and Hall/CRC. [Google Scholar]
  9. Davis, Richard A., and Heng Liu. 2016. Theory and inference for a class of nonlinear models with application to time series of counts. Statistica Sinica 26: 1673–707. [Google Scholar] [CrossRef]
  10. Du, Jin-Guan, and Yuan Li. 1991. The integer-valued autoregressive (INAR(p)) model. Journal of Time Series Analysis 12: 129–42. [Google Scholar]
  11. Dunn, James E. 1967. Characterization of the bivariate negative binomial distribution. Journal of the Arkansas Academy of Science 21: 77–86. [Google Scholar]
  12. Esquível, Manuel L. 2008. Probability generating functions for discrete real–valued random variables. Theory of Probability and Its Applications 52: 40–57. [Google Scholar] [CrossRef]
  13. Ferland, René, Alain Latour, and Driss Oraichi. 2006. Integer-valued GARCH processes. Journal of Time Series Analysis 27: 923–42. [Google Scholar] [CrossRef]
  14. Fokianos, Konstantinos. 2012. Count time series. In Handbook of Statistics 30: Time Series—Methods and Applications. Edited by Tata Subba Rao, Suhasini Subba Rao and Calyampudi Radhakrishna Rao. Amsterdam: Elsevier, pp. 315–47. [Google Scholar]
  15. Fokianos, Konstantinos, and Michael H. Neumann. 2013. A goodness–of–fit tests for Poisson count processes. Electronic Journal of Statistics 7: 793–819. [Google Scholar] [CrossRef]
  16. Fokianos, Konstantinos, Andres Rahbek, and Dag Tjøstheim. 2009. Poisson autoregression. Journal of the American Statistical Association 104: 1430–39. [Google Scholar] [CrossRef]
  17. Fokianos, Konstantinos, and Dag Tjøstheim. 2012. Nonlinear Poisson autoregression. Annals of the Institute of Statistical Mathematics 64: 1205–25. [Google Scholar] [CrossRef]
  18. Franke, Jürgen, and T. Subba Rao. 1995. Multivariate First-Order Integer-Valued Autoregressions. Technical Report. Kaiserslautern: Universität Kaiserslautern. [Google Scholar]
  19. Freeland, R. Keith. 1998. Statistical Analysis of Discrete Time series with Application to the Analysis of Workers’ Compensation Claims Data. Ph.D. thesis, Management Science Division, Faculty of Commerce and Business Administration, University of British Columbia, Vancouver, BC, Canada. Available online: https://open.library.ubc.ca/cIRcle/collections/ubctheses (accessed on 25 February 2021).
  20. Freeland, R. Keith, and Brendan P. M. McCabe. 2004. Analysis of low count time series data by Poisson autoregression. Journal of Time Series Analysis 25: 701–22. [Google Scholar] [CrossRef]
  21. Giacomini, Raffaella, Dmitris Politis, and Halbert White. 2013. A warp-speed method for conducting Monte Carlo experiments involving bootstrap estimators. Econometric Theory 29: 567–89. [Google Scholar] [CrossRef]
  22. Heinen, Andréas, and Erick Rengifo. 2007. Multivariate autoregressive modeling of time series count data using copulas. Journal of Empirical Finance 14: 564–583. [Google Scholar] [CrossRef]
  23. Hudecová, Šárka, Marie Hušková, and Simos G. Meintanis. 2015. Tests for time series of counts based on the probability generating function. Statistics 49: 316–37. [Google Scholar] [CrossRef]
  24. Ibragimov, Il’dar Abdulovich, and Rafail Zalmonovich Chasminskij. 1981. Statistical Estimation, Asymptotic Theory. New York: Springer. [Google Scholar]
  25. Jiménez-Gamero, M. Dolores, Sangyeol Lee, and Simos G. Meintanis. 2020. Goodness-of-fit tests for parametric specifications of conditionally heteroscedastic models. TEST 29: 682–703. [Google Scholar] [CrossRef]
  26. Kim, Hanwool, and Sangyeol Lee. 2017. On first order integer-valued autoregressive process with Katz family innovations. Journal of Statistical Computation and Simulation 87: 546–62. [Google Scholar] [CrossRef]
  27. Kocherlakota, Subrahmaniam, and Kathleen Kocherlakota. 1992. Bivariate Discrete Distributions. New York: Marcel Dekker Inc. [Google Scholar]
  28. Lakshminarayana, J., S. N. Narahari Pandit, and K. Srinivasa Rao. 1999. On a bivariate Poisson distribution. Communications in Statistics—Theory and Methods 28: 267–76. [Google Scholar] [CrossRef]
  29. Latour, Alain. 1997. The multivariate GINAR(p) process. Advances in Applied Probability 29: 228–48. [Google Scholar] [CrossRef]
  30. Lee, Youngmi, Sangyeol Lee, and Dag Tjøstheim. 2018. Asymptotic normality and parameter change test for bivariate Poisson INGARCH models. TEST 27: 52–69. [Google Scholar] [CrossRef]
  31. Leucht, Anne, Jens-Peter Kreiss, and Michael H. Neumann. 2015. A model specification test for GARCH(1,1) processes. Scandinavian Journal of Statistics 42: 1167–93. [Google Scholar] [CrossRef]
  32. Liu, Heng. 2012. Some Models for Time Series of Counts. Ph.D. thesis, Columbia University, New York, NY, USA. [Google Scholar]
  33. Mamode Khan, Naushad Ali, Yuvraj Sunecher, Vandna Jowaheer, Miroslav M. Ristić, and Maleika Heenaye–Mamode Khan. 2019. Investigating GQL–based inferential approaches for non-stationary BINAR(1) model under different quantum of over-dispersion with application. Computational Statistics 34: 1275–313. [Google Scholar] [CrossRef]
  34. McKenzie, Eddie. 1985. Some simple models for discrete variate time series. Water Resources Bulletin 21: 645–50. [Google Scholar] [CrossRef]
  35. McKenzie, Eddie. 2003. Discrete variate time series. Stochastic Processes: Modelling and Simulation. In Handbook of Statistics. Amsterdam: Elsevier, vol. 21, pp. 573–606. [Google Scholar]
  36. Meintanis, Simos G., and Dimitris Karlis. 2014. Validation tests for the innovation distribution in INAR time series models. Computational Statistics 29: 1221–41. [Google Scholar] [CrossRef]
  37. Neumann, Michael H. 2011. Absolute regularity and ergodicity of Poisson count processes. Bernoulli 17: 1258–84. [Google Scholar] [CrossRef]
  38. Partrat, Christian. 1994. Compound model for two dependent kinds of claim. Insurance: Mathematics and Economics 15: 219–31. [Google Scholar] [CrossRef]
  39. Pedeli, Xanthi, and Dimitris Karlis. 2011. A bivariate INAR(1) process with application. Statistical Modelling 11: 325–49. [Google Scholar] [CrossRef]
  40. Pedeli, Xanthi, and Dimitris Karlis. 2013a. On composite likelihood estimation of a multivariate INAR(1) model. Journal of Time Series Analysis 34: 206–20. [Google Scholar] [CrossRef]
  41. Pedeli, Xanthi, and Dimitris Karlis. 2013b. On estimation of the bivariate Poisson INAR process. Communications in Statistics—Simulation and Computation 42: 514–33. [Google Scholar] [CrossRef]
  42. Pedeli, Xanthi, and Dimitris Karlis. 2013c. Some properties of multivariate INAR(1) processes. Computational Statistics & Data Analysis 67: 213–25. [Google Scholar]
  43. Popović, Predrag M., Aleksandar S. Nastić, and Miroslav M. Ristić. 2018. Residual analysis with bivariate INAR models. REVSTAT 16: 349–64. [Google Scholar]
  44. R Core Team. 2019. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. [Google Scholar]
  45. Santos, Cláudia, Isabel Pereira, and Manuel G. Scotto. 2019. On the theory of periodic multivariate INAR processes. Statistical Papers. [Google Scholar] [CrossRef]
  46. Schweer, Sebastian. 2016. A goodness–of–fit test for integer valued autoregressive processes. Journal of Time Series Analysis 37: 77–98. [Google Scholar] [CrossRef]
  47. Schweer, Sebastian, and Christian H. Weiss. 2014. Compound Poisson INAR(1) processes: Stochastic properties and testing for overdispersion. Computational Statistics & Data Analysis 77: 267–84. [Google Scholar]
  48. Scotto, Manuel G., Christian H. Weiss, and Sónia Gouveia. 2015. Thinning-based models in the analysis of integer-valued time series: A review. Statistical Modelling 15: 590–618. [Google Scholar] [CrossRef]
  49. Shi, Peng, Xiaoping Feng, and Jean-Philippe Boucher. 2016. Multilevel modeling of insurance claims using copulas. The Annals of Applied Statistics 10: 834–36. [Google Scholar] [CrossRef]
  50. Steutel, Fred W., and Klaas van Harn. 1979. Discrete analogues of self-decomposability and stability. Annals of Probability 7: 893–99. [Google Scholar] [CrossRef]
  51. Vernic, Raluca. 1997. On the bivariate generalized poisson distribution. ASTIN Bulletin 27: 22–32. [Google Scholar] [CrossRef]
  52. Weiss, Christian H. 2018a. Discrete-Valued Time Series. Hoboken: John Wiley & Sons. [Google Scholar]
  53. Weiss, Christian H. 2018b. Goodness–of–fit testing of a count series’ marginal distribution. Metrika 81: 619–51. [Google Scholar] [CrossRef]
  54. Weiss, Christian H., Annika Homburg, and Pedro Puig. 2019. Testing for zero inflation and overdispersion in INAR(1) models. Statistical Papers 60: 473–98. [Google Scholar] [CrossRef]
  55. Weiss, Christian H., and Sebastian Schweer. 2015. Detecting overdispersion in INARCH(1) processes. Statistica Neerlandica 69: 281–97. [Google Scholar] [CrossRef]
  56. Zhu, Rong, and Harry Joe. 2006. Modelling count data time series with Markov processes based on binomial thinning. Journal of Time Series Analysis 27: 725–38. [Google Scholar] [CrossRef]
Figure 1. Monthly number of claims. (a) Soft tissue injuries claims (black dots) and dermatitis claims (gray triangles). (b) Cuts claims.
Figure 1. Monthly number of claims. (a) Soft tissue injuries claims (black dots) and dermatitis claims (gray triangles). (b) Cuts claims.
Econometrics 09 00010 g001
Table 1. Size of the test for testing Poisson BINAR in H 0 and BPAR in H 0 ˜ .
Table 1. Size of the test for testing Poisson BINAR in H 0 and BPAR in H 0 ˜ .
Poisson BINARBPAR
A = A 1 A = A 2 B = A 1 B = A 2
n α 0.010.050.100.010.050.100.010.050.100.010.050.10
1000.0000.0100.0300.0110.0390.1040.0020.0290.0750.0130.0770.145
2000.0010.0100.0480.0170.0650.1140.0060.0400.0760.0280.0940.168
5000.0020.0250.0690.0050.0490.1050.0050.0500.0940.0140.0530.127
Table 2. Power of the test for H 0 for data from BINAR model with negative Binomial marginals.
Table 2. Power of the test for H 0 for data from BINAR model with negative Binomial marginals.
A = A 1 A = A 2
r = 5 r = 10 r = 5 r = 10
n α 0.01 0.05 0.10 0.01 0.05 0.10 0.01 0.05 0.10 0.01 0.05 0.10
1000.0050.1390.3280.0040.0650.1350.3570.5810.6940.0800.1960.296
2000.1200.6390.9020.0090.1190.3310.7150.8970.9550.1670.4320.540
5000.9511.0001.0000.1320.8840.9841.0001.0001.0000.6370.8610.920
Table 3. Power of the test for H 0 ˜ with data generated by a negative Binomial BINARCH model.
Table 3. Power of the test for H 0 ˜ with data generated by a negative Binomial BINARCH model.
B = A 1 B = A 2
r = 30 r = 50 r = 30 r = 50
n α 0.01 0.05 0.10 0.01 0.05 0.10 0.01 0.05 0.10 0.01 0.05 0.10
1000.0850.3570.4690.0600.1690.2650.3470.6720.8210.1240.3500.526
2000.4460.7400.8330.1430.3450.4600.4480.8710.9560.2030.5680.746
5000.9420.9890.9940.5020.7010.7970.8100.9961.0000.5240.8740.957
Table 4. Power of the test of Poisson BINAR in H 0 for data following BPAR and vice versa.
Table 4. Power of the test of Poisson BINAR in H 0 for data following BPAR and vice versa.
Test H 0 for Data from BPARTest H 0 ˜ for Data from BINAR(1)
B = A 1 B = A 2 A = A 1 A = A 2
n α 0.01 0.05 0.10 0.01 0.05 0.10 0.01 0.05 0.10 0.01 0.05 0.10
1000.0030.0220.0610.6130.7910.8520.0030.0240.0820.0110.2440.626
2000.0060.0270.1020.9550.9880.9940.0040.0580.1680.1090.8100.942
5000.0090.1360.2781.0001.0001.0000.0300.1510.3460.7240.9800.995
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Back to TopTop