Specifications tests for count time series models with covariates

We propose a goodness-of-fit test for a class of count time series models with covariates which includes the Poisson autoregressive model with covariates (PARX) as a special case. The test criteria are derived from a specific characterization for the conditional probability generating function and the test statistic is formulated as a $L_2$ weighting norm of the corresponding sample counterpart. The asymptotic properties of the proposed test statistic are provided under the null hypothesis as well as under specific alternatives. A bootstrap version of the test is explored in a Monte--Carlo study and illustrated on a real data set on road safety.


Introduction
Let tY t u be a time series of counts, i.e. a time series of variables taking values in N 0 , and let tX t u be a m-dimensional vector time series of exogenous covariates.Denote as I t the information set available up to time t and let tZ t u be a time series of d-dimensional random vectors such that Z t may contain past values of both endogenous as well as exogenous covariates.Our aim is to develop procedures for parametric specifications of the conditional distribution of Y t given I t´1 , i.e. we are interested in testing the null hypothesis H 0 : Y t |I t´1 " F λt , λ t " hpZ t ; θq, where F λ denotes a fixed count distribution with a mean λ, h is a specified function and θ is an unspecified vector of parameters.Notice that the manner in which λ t depends on the covariate vector Z t is fixed under H 0 apart from a finite-dimensional parameter θ which must be estimated from the data.We will restrict ourselves to a class of distributions tF λ u such that λ ą 0 is the mean of F λ and the family tF λ u satisfies λ 1 ď λ 2 ñ F λ 1 pyq ě F λ 2 pyq for all y P IR. (2) The condition (2) holds, for instance, for Poisson distribution with mean λ, for Negative Binomial distribution with mean λ and fixed ϕ, or for a zero-inflated Poisson distribution.Aknouche and Francq (2021) provide conditions for existence of a stationary and ergodic solution to (1) with a conditional distribution satisfying (2).If F λ denotes the Poisson distribution with mean λ and if no exogenous variables are considered then (1) corresponds to the INGARCH model of (Fokianos and Tjøstheim, 2009).An extension to this model, abbreviated sometimes as PARX, incorporating also exogenous variables was proposed and studied by Agosto et al. (2016), and the model has already found a number of important extensions and applications in various fields.In finance, Agosto et al. (2016) used PARX to model US corporate default counts, while Agosto and Raffinetti (2019) analyzed Italian corporate default counts data in the real estate sector.Other applications are in sports (Angelini and De Angelis, 2017), and disease modelling (Agosto and Giudici, 2020).Alternative specifications for F λ include the Negative Binomial (NB) distribution (Christou and Fokianos, 2014;Zhu, 2011;Liboschik et al., 2017), and the zero-inflated Poisson distribution (Aknouche and Francq, 2021), among others.
The problem of testing goodness-of-fit within INGARCH models has been considered by many authors, for instance (Fokianos and Neumann, 2013;Weiss and Schweer, 2015;Schweer, 2016;Hudecová et al., 2015Hudecová et al., , 2021) ) to mention some of them.We propose a method for testing H 0 for an arbitrary family of conditional distributions satisfying (2), a possibly non-linear function h and a setup with external covariates.This, of course, involves testing PARX and INGARCH specifications as special cases.Our test is based on a characterization of the conditional mean by Bierens (1982), which has ben successfully used in construction of various specification tests, e.g., see (Escanciano and Velasco, 2006;Hlávka et al., 2017) and references therein.In contrast to these works, we reformulate the characterization in terms of the conditional probability generating function (PGF) of Y t given Z t .PGF-based procedures are specifically tailored for count responses and are often convenient from the computational point of view.Earlier PGF-based methods, suggested by Hudecová et al. (2015), Schweer (2016) and Hudecová et al. (2021), are suitable for an analysis of time series without exogenous covariates.In contrast to this, the new tests suggested herein, although still PGF-based, adopt an alternative formulation which allows them to be readily applied to more complicated dynamic models.
The rest of the paper is outlined as follows.Section 2 introduces the considered model.The test statistic is proposed in Section 3, while its asymptotic properties are studied in Section 4. Section 5 describes the practical computation of the test and presents a bootstrap version.An application to the Poisson autoregression with covariates (PARX) is discussed in Section 6. Section 7 presents results of a Monte Carlo study.A real data application of road safety is included in Section 8. We finally conclude with discussion of our finding and further outlook in Section 9. Proofs are postponed to Appendix.
Notice that the random vector Z t implicitly depends on the value of the model parameter θ, but we do not stress this dependence if not necessary.If q " 0 then the vector Z t contains only past values of Y t and X t , and the conditional mean is specified using a (nonlinear) autoregressive model with exogenous covariates (ARX).
If the exogenous covariate process tX t u is stationary and ergodic, then there exists a stationary and ergodic solution tY t u to the model (3) with the link function h given by ( 5) provided that the function r satisfies |rpy 1 , λ 1 ; θq ´rpy 2 , λ 2 ; θq| ď for all y i " py 1i , . . ., y 1p q J P r0, 8q p , λ i " pλ i1 , . . ., λ iq q J P r0, 8q q , i " 1, 2, and for some α i , β j , i " 1, . . ., p, j " 1, . . ., q such that see (Aknouche and Francq, 2021).Model (5) includes as a special case a linear model, which can be for non-negative components of X t formulated as (3) for a link function hpy, l, x; θq " ω `αJ y `βJ l `πJ x (8) for θ " pω, α J , β J , π J q J with ω ą 0, α P r0, 1s p , β P r0, 1s q and π P p0, 8q m .The condition for existence of a stationary and ergodic solution then reduces to (7).In that case 3 Goodness-of-fit test Let tpY t , X t´1 qu T t"1 be observed data.We develop test for the null hypothesis 3) with a specified class of distributions F λ , and specified function hp¨; θ 0 q for some inner point θ 0 P Θ 0 , against a general alternative H 1 : not H 0 , where Θ 0 stands for a compact subset of the parametric space Θ.The model orders pp, qq are implicitly specified by the function h, so they are assumed to be known under the null hypothesis.
Our test utilizes the characterization of Bierens (1982) which for a random variable U and a d-dimensional random vector W may be formulated via the following equivalence relation where i " ?´1 and µ 0 p¨q denotes a specific regression function possibly involving an unknown parameter vector.
In the current context, we rephrase (10) for model (3) in terms of the conditional probability generating function (PGF).Recall that a PGF of a count variable Y is defined as g Y puq :" Epu Y q for u P IR for which the expectation is finite.The PGF is always well defined for u P r0, 1s and the distribution of Y is uniquely determined from the values of g Y from a neighborhood of the origin.Let g λ puq be the PGF corresponding to F λ .Under H 0 it holds that E ru Yt |I t´1 s " E ru Yt |Z t s " g λt puq " g hpZt;θ 0 q puq, u P r0, 1s.
Hence, the characterization featuring in (10) can be under H 0 reformulated as E rtu Yt ´gλt puque iv J Zt s " 0, @u P p0, 1q, v P IR d , λ t " hpZ t ; θ 0 q. (11) Our test procedure will be based on a weighted L 2 norm of an estimator of the expectation in (11 where The test statistic ∆ T,W in ( 13) is an L 2 -type distance statistic which, at least for large sample size T , should be small under the null hypothesis H 0 and large under alternatives.Therefore large values of the test statistic indicate that the null hypothesis is violated.The estimator p θ can be either a maximum likelihood estimator (MLE), or, more generally, a Poisson quasimaximum likelihood estimator (QMLE).The latter method has been studied in Ahmad and Francq (2016) and Aknouche and Francq (2021), where it is shown that the resulting estimator is consistent and asymptotically normal under mild regularity conditions.For practical computation of ∆ T,W , it is useful to notice that If the weight function can be decomposed as W pu, vq " wpuqωpvq, with wp¨q and ωp¨q satisfying some additional assumptions, then further simplification is possible; see Section 5 for more details on the practical aspects of computing the test statistic ∆ T,W and Section 6 for the special case of a linear PARX model.
Remark 3.1.Note that the integral in (13) is computed for u P r0, 1s, which is an approach consistent with a number of previous works.Alternatively, the integration can be carried over the interval r´1, 1s, or over an interval I Ă r´1, 1s containing the origin.
Remark 3.2.Our vector of covariates Z t defined as in (4) may also include an inter-cept.In such case, q Z t " p1, Y t´1 , . . ., Y t´q , λ t´1 , . . ., λ t´q , X J t´1 q J " p1, Z J t q J , where Z t is given (4) and d " p `q `m, then it follows from (14) and (13) that the resulting test statistic q ∆ T, | W based on a weight function | W : IR d`2 Ñ p0, 8q is equal to ∆ T,W computed from Z t and the weight function

Asymptotic distribution
This section describes the asymptotic behavior of the test statistic ∆ T,W under the null hypothesis (Section 4.1) as well as under some alternatives (Section 4.3).Since some of the assumptions can be substantially weakened for an autoregressive type of model, where q " 0, the behavior under the null for these models is treated in more detail in Section 4.2.

Behavior under the null hypothesis
We start with some assumptions on the link function and the weight function.
Assumption A2.Let tX t u be a strictly stationary and ergodic sequence of m-dimensional random vectors.Let tF λ u satisfy (2) and let tY t u be the stationary and ergodic solution of (3) with λ t " hpZ t ; θ 0 q, where Z t is defined by (4), for a true parameter θ 0 P intpΘ 0 q, where Θ 0 is a compact subset of Θ Ă IR K .Assume further that E }Z t } 2 ă 8.
(Here intpAq stands for the interior of a set A and }x} is the Euclidean norm of a vector x).
The existence of a stationary ergodic solution from Assumption A2 is proved for h satisfying Assumption A1 in (Aknouche and Francq, 2021, Theorem 3.3).Under the same conditions, the process pY t , X J t q J is stationary and ergodic and the same holds for the multivariate process tZ t u.Sufficient conditions for the existence of moments for the linear case with p " q " 1 are provided in (Aknouche and Francq, 2021, Theorem 3.2).Also notice that we assume θ 0 lies in the interior of the parameter space Θ 0 ; for estimation and specification testing which allow the parameter to lie on the boundary of the parameter space we refer to Francq and Zakoïan (2007) and Cavaliere et al. (2023).
Assumption A3.Let r and π be the functions from Assumption A1.Let rpy, ¨, ¨q be three times continuously differentiable on p0, 8q q ˆΘ for all y P p0, 8q p and denote r λ py, λ; θq " Brpy, λ; θq Bλ , r θ py, λ; θq " Brpy, λ; θq Bθ for λ P IR q and θ P Θ. Suppose r λ py, λ; θq " dpθq is constant with respect to y and λ such that where ρp¨q is a spectral radius and I k is a k ˆk identity matrix.
where Y t´1:t´p " pY t´1 , . . ., Y t´p q J and Λ t´1:t´q " pλ t´1 , . . ., λ t´q q J .Assumption A3 is related to the smoothness of the mean function and the required conditions ensure that the effect of the initial values is asymptotically negligible.A similar set of assumptions is invoked by Aknouche and Francq (2021) to prove the asymptotic normality of the QMLE.Notice that Assumption A3 is always satisfied by a linear function (8).
Assumption A4.Assume that the PGF g λ puq of the distribution F λ , considered as a function g z puq of pu, zq, is twice continuously differentiable and that
Assumption A4 ensures that the PGF g λ as a function λ is sufficiently smooth.For a Poisson distribution, one can take M 1 " M 2 " 1 and Hpzq " z `1.
Assumption A5.Let pY t , λ t pθqq be the stationary solution of (3) for parameter θ P Θ.
Let the second derivative B 2 λtpθq BθBθ J exist and be continuous in θ and let there exist a neighborhood Vpθ 0 q of θ 0 such that where θ 0 P Θ 0 is the true value of θ, ts t pθ 0 , Y t , I t´1 qu is a stationary martingale difference sequence with respect to the filtration tI t u with a finite variance such that T ´1{2 ř T t"1 s t pθ 0 , Y t , I t´1 q D Ñ Np0, Σq for some K ˆK positive semidefinite matrix Σ.
Assumption A6 is satisfied under some additional regularity conditions by the MLE, see Agosto et al. (2016) for a linear Poisson model.More generally, the estimator p θ T may be the Poisson QMLE, see Ahmad and Francq (2016); Aknouche and Francq (2021).
Define for v 1 , v 2 P R d and u P r0, 1s Note that βpu, v 1 q is finite for all u P r0, 1s and v P IR d due to Assumptions A4 and A5.The following theorem describes the asymptotic distribution of the proposed test statistic under the null hypothesis H 0 .
Theorem 4.1.Let Assumptions A1-A7 be satisfied.Then the test statistic ∆ T,W converges in distribution for T Ñ 8 to a random variable where Z " tZpu, vq, s P r0, 1s, v P IR d u is a centered Gaussian process with a covariance function The assertion provides an approximation of the distribution of ∆ T,W under the null hypothesis.It follows from Theorem 4.1 that ∆ T,W is asymptotically distributed as an infinite weighted sum of χ 2 1 distributed random variable, where, however, the weights depend on the unknown parameters via a highly non-trivial way.Hence, it is extremely complicated to use this asymptotic distribution in order to determine critical points and actually carry out the test.Instead, in Section 5 we recommend the use of a bootstrap test which is straightforward to apply.

Behavior under the null for an ARX type of model
An important special case of the model (3) is obtained for q " 0, when the conditional mean λ t depends only on the past values of Y t and on X t´1 via some kind of an ARX structure.Some of the assumptions specified in the previous section can be weakened for these models, so we treat this situation in more detail.If q " 0, then hpy, x; θq " rpy; θq `πpx; θq, y P r0, 8q p , x P IR m , θ P IR K for some functions r and π, and Z t " pY t´1 , . . ., Y t´q , X J t´1 q J depends solely on past values of Y t and on X t´1 .This fact enables us to reduce and simplify the assumptions required in A3, A5 and A7, into the following weaker versions: Assumption A3 ˚.Let r : r0, 8q p ˆIR K Ñ p0, 8q and π : IR K`m Ñ p0, 8q be functions such that rpy; ¨9 q and πpx; ¨q are twice continuously differentiable for all y P r0, 8q p and x P IR m .Assume that there exists a neighborhood Vpθ 0 q of θ 0 such that for all i, j " 1, . . ., K, where Y t´1:t´p is defined below (15) and Y t , X t are the stationary variables from Assumption A2.
where the function Hp¨q is from Assumption A4 and pY t , λ t q are the stationary variables from Assumption A2.

Behavior under some alternatives
The limit behavior of the test statistic is provided also for some types of alternatives.To this end analogously to the notation below (15), we write Y k:l for pY k , Y k´1 , . . ., Y k´pk´lq q J and Λ k:l for pλ k , λ k´1 , . . ., λ k´pk´lq q J , for any integers pk, lq with k ě l.
Let p 0 ą 0, q 0 ě 0, K 0 ě 1 and write h 0 : r0, 8q p 0 `q0 ˆIR m`K 0 Ñ p0, 8q for a particular link function.Consider the null hypothesis H 0 whereby the link function is specified by h 0 and a fixed family of distributions tF λ , λ ą 0u, that is H 0 : there exists θ P intpΘ 0 q such that Y t |I t´1 " F λt λ t " h 0 pY t´1:t´p 0 , Λ t´1:t´q 0 , X t´1 ; θq, for Θ 0 Ă Θ a compact set.A general alternative admits an incorrect specification of the conditional mean function as well as the conditional distribution.Specifically suppose that the actual model is given by for some family of distributions tF A λ , λ ą 0u, a link function h A : r0, 8q p A `qA ÎR m`K A Ñ p0, 8q, with p A ą 0, q A ě 0, K A ě 1, and for some θ A P IR K A .Hence if the true pair pF A λ , h A q differs from pF λ , h 0 q then the model is not correctly specified under H 0 .Let ∆ T,W be the test statistic defined by ( 13) for testing H 0 , based on some estimator p θ T computed under H 0 , with the sequence t p λ t,0 u defined recursively using the link function h 0 , and a d-dimensional sequence t p Z t,0 u where d " p 0 `q0 `m.
Assume that the PGF g λ puq of F λ satisfies |Bg λ puq{Bλ| ă M 1 for all u P r0, 1s and all λ P p0, 8q for some M 1 ą 0 and that there exists θ 0 P intpΘ 0 q such that Consider first an alternative whereby the mean function h is correctly specified, but the conditional distribution is different from the hypothesized one, i.e.F λ ‰ F A λ , and write g A λ for the PGF corresponding to F A λ .Note that in this case the true parameter θ 0 can still be consistently estimated using, for instance, the Poisson QMLE, and accordingly we may assume that (18) holds with θ 0 and θ A being identical.
Consider now a situation whereby the conditional distribution is correctly specified under the null hypothesis but the link function is misspecified, i.e.F A λ " F λ but h A ‰ h 0 .We provide the asymptotic behavior of ∆ T,W in two practically important situations, when either h 0 corresponds to a linear model, or when h 0 is non-linear with q 0 " 0, i.e. the tested link h 0 corresponds to an ARX type of a model.For both situations we show that the test based on ∆ T,W is consistent.
Let us start with the latter (non-linear) case, whereby h 0 py, x; θq " r 0 py; θq π0 px; θq for some r : r0, 8q p 0 ˆIR K Ñ p0, 8q and π : IR K 0 `m Ñ p0, 8q and p Z t,0 " Z t,0 " pY J t´1:t´p 0 , X J t´1 q J for t ě p 0 `1.Notice that nothing is assumed about the true orders pp A , q A q, so these can be completely arbitrary.
Theorem 4.4.Let h 0 : r0, 8q p 0 ˆIR K 0 `m Ñ p0, 8q be a function satisfying A1 on Θ Ă IR K 0 for some p 0 ą 0, K 0 ě 1.Let tpY t , λ t,A qu and tX t u satisfy A2 with a link function h A and distribution F λ .Assume that (17), ( 18) and A7 ˚hold, and there exists a neighborhood Vpθ 0 q of θ 0 such that Consider now a situation, where h 0 corresponds to a linear model with orders p " q " 1, that is h 0 py, λ, x, ; θq " ω `αy `βλ `γJ x (19) for θ " pω, α, β, γ J q J .Let Θ be a subset of IR 3`m such that for any θ P Θ, θ " pω, α, β, γ J q J , it holds that ω, α, β ą 0 and α `β ă 1.Again, nothing is assumed about the true orders pp A , q A q and only the standard assumptions are posed for the true link function h A .
In order to further discuss the conditions of consistency, write ∆ ζ 8,W for each probability limit of ∆ T,W {T figuring in Theorems 4.3, 4.4 and 4.5.Since ζ is a continuous function on p0, 1q ˆIR d , ∆ ζ 8,W " 0 only if ζpu, vq is equal to zero identically in pu, vq P p0, 1q ˆIR d .On the other hand if ζpu, vq ‰ 0 for some pu, vq P p0, 1q ˆIR d , then it follows from assumption A7 that ∆ ζ 8,W ą 0, and consequently the corresponding test that rejects the null hypothesis H 0 for large values of ∆ T,W is consistent, in the sense that the rejection probability is equal to one asymptotically under the violations of H 0 considered herein.
In this connection note that in Theorem 4.3 g A λ ‰ g λ , and therefore g A λ 1 puq ´gλ 1 puq is generally a non-zero random variable.If in addition g A λ puq ą g λ puq for all u P p0, 1q (which is the case for F λ Poisson and F A λ a Negative Binomial), then g A λ 1 puq´g λ 1 puq ą 0 for all u P p0, 1q.In the case of Theorems 4.4 and 4.5, if the condition (2) holds sharply in the sense that tλ 1 ă λ 2 ùñ g λ 1 puq ą g λ 2 puq for all u P p0, 1qu (as, e.g., for the Poisson distribution), then λ 1,A ´λ1,0 ‰ 0 a.s.implies that g λ 1,A puq ´gλ 1,0 puq ‰ 0 a.s.for all u P p0, 1q.Also observe that αpv 1 , v 2 q " 0 only if v J 1 v 2 " p3π{4q `kπ, k " 0, ˘1, ˘2, ... .Hence in each of the three cases of alternatives considered by Theorems 4.3-4.5, the random function within the corresponding expectation in the definition of ζpu, vq is generally non-zero.This of course does not rule out the case of ζp¨, ¨q being identically equal to zero, but alludes to the very special circumstances under which the test based on ∆ T,W will miss any of the alternatives considered herein, at least for large sample size T .In this connection, the small sample behavior of the test based on ∆ T,W under various realizations of the alternatives considered in Theorems 4.3-4.5 is investigated by means of a Monte Carlo simulation study in Section 7 and shows that the new test has non-trivial power against several instances of violations of the null hypothesis.
5 Practical computation and bootstrap test

Analytic computations
Let the weight function W be of the form W pu, vq " wpuqωpvq.Due to ( 14) the test statistic ∆ T,W takes the form with and r K ω,ts :" r K ω p p Z t ´p Z s q, where r K ω pzq " The integral r K ω p¨q can be easily analytically computed if the weight function ωp¨q corresponds to the density of a spherical distribution.Then we immediately have r K ω pzq " Ξp}z}q with the function Ξp¨q being the generator of the characteristic function of the chosen spherical distribution.Hereafter we use Ξ :" Ξ γ,η , with Ξ γ,η pξq " e ´γξ η , ξ, γ ą 0, η P p0, 2s, corresponding to the spherical stable characteristic function.Notice that Assumption A7 ˚is satisfied for all γ ą 0 and η P p0, 2s,, while only η " 2 fulfills the requirements of Assumption A7.
The integral K w,ts can be further simplified depending on the PGF g λ and the choice of wpuq.Commonly used function is w ρ puq " u ρ for ρ ě 0. Another possibility is to take a weight function proportional to a density of a beta distribution w ρ,κ puq " u ρ p1 ´uq κ for ρ, κ ě 0, which enables to regulate the weight given to neighborhoods of the origin and point u " 1.The case of linear Poisson PARX model and weight function w ρ is considered in more detail in Section 6.

Bootstrap version of the test
Theorem 4.1 provides the asymptotic distribution of the test statistics ∆ T,W under the null hypothesis H 0 , but this distribution depends on several unknown quantities.To get an approximation for the desired critical value one can generate the Gaussian process described in Theorem 4.1 with the unknown quantities replaced by their estimators and then calculate the integral and its quantiles.However, here we propose a parametric bootstrap.
(iii) Generate the bootstrap observations Y 1 , . . ., Y T recursively from the distribution F r λ t p p θ T q with mean p λ t p p θ T q defined as p λ t pθq with θ replaced by p θ T and Y t´i replaced by Y t ´i and X t replaced by X t .
The simulation study in next section uses overlapping block bootstrap with fixed block length l with l « n 1{3 .The initial values for computation of p λ t and p λ t are chosen as p λ 0 " 0, Y 0 " Y 1 and X 0 " X 1 .
Using similar tools as in Neumann ( 2021) while examining step-by-step the proof of Theorem 4.1, one will come to the conclusion that the approximation of the critical value by the above parametric bootstrap is asymptotically correct if the original data follow the null hypothesis.

Application to PARX model
The linear Poisson autoregression with exogenous covariates (PARX) assumes that F λ is Poisson, so the corresponding PGF takes form g λ puq " e λpu´1q , and the conditional mean is modelled as (9) for some p ě 1 and q ě 0 and an unspecified model parameter θ.In the following this particular hypothesis is referred to as H P 0 .Then where p λ t are computed recursively as described in Section 3. The integral K w,ts from (20) can be computed as K w,ts :" K w pp ε t , p ε s q for where λ t and y t , t " 1, 2, are used as generic notations for pairs of parameters and associated responses, respectively.To proceed further hereafter we adopt the weight function wpuq " w ρ puq " u ρ , ρ ě 0, and write K pρq pε 1 , ε 2 q :" K pρq py 1 , y 2 ; λ 1 , λ 2 q, for the resulting integral which after some straightforward algebra is rendered as where Ipa, θq " ż 1 0 u a e θu du " p´θq ´a θ pΓp1 `a, ´θq ´Γp1 `aqq .
Then the test statistic figuring in (20) may be written as with p K pρq ts " K pρq py t , y s ; p λ t , p λ s q and p Z t defined recursively as described in Section 3.

Simulations
This simulation study explores the finite sample behavior of the suggested bootstrap test based on ∆ T,W for the null hypothesis H P 0 , where the conditional distribution F λ is Poisson.The conditional mean is specified either as an ARX(1) kind of model or as a GARCHX(1,1) type of model The exogenous series tX t u was generated from a stationary AR(1) model X t " ρ X X t´1 ὲt for ρ X " 0.5 and tε t u iid centred normal with variance p1 ´ρ2 X q ´1.The model parameters were estimated by the maximum likelihood method.The performance of the testing procedure was investigated under the null hypothesis.Moreover the alternative hypothesis settings considered are the cases whereby we have: (a) Incorrect model specification.Observations were generated from a model of an incorrect order, while the function π from (5) and the conditional distribution F λ were specified correctly.Namely, the test for (S1) was conducted for observations from the GARCHX(1,1) model (S2) and from a ARX(2) kind of model which was used as an alternative for testing (S2) as well.
(b) Incorrect conditional distribution.Observations were generated from a model with a correctly specified conditional mean λ t but the conditional distribution was the Negative Binomial NBpλ t , rq with a dispersion parameter r.Recall that if a random variable follows NBpλ, rq then its mean is λ and variance is λp1 `λ{rq.
For model (S1) all these settings satisfy the assumptions required by the weight function, but for (S2) we have theoretical justification only for the last two choices of pγ, ηq due to the stronger requirements in Assumption A7.Despite this fact, we explore the behavior of the bootstrap test for all the choices of pγ, ηq listed above, so we can investigate how sensitive the procedure is to violation of Assumption A7.
In addition, the behavior of the test statistic ∆ T,W was compared to tests based on the following two test statistics: The significance of ∆ T,W and ∆ piq T , i " 0, 1, was evaluated using the parametric bootstrap described in Section 5 with resampling size B " 499.The simulations were conducted in the R-computing environment for M " 1 000 repetitions.
The results in Table 3 are quantitatively similar with those of Table 2 for the NB distributional alternative and the "abs" misspecification of the π function, and qualitatively similar for the other two types of alternatives, but in this case it is harder to identify a ARX(2) alternative while the power against a "sin" misspecification of the π function is somewhat lower for all tests.Comparison-wise, the conclusions are also analogous to those drawn from the results of Table 2, overall favouring the tests based on ∆ T,W , with γ " η " 1{4 and γ " η " 1{2.It is interesting to see that the procedure leads to reasonable results even if η ă 2. The corresponding bootstrap test keeps the prescribed level, and this choice seems to be even preferable against some alternatives.

Real data analysis
Road crashes have not only a devastating impact on victims and their families, but they also bear a heavy economic cost, see, e.g., (Peden et al., 2004, Table 2.9).The total ∆ T,W with ρ " 0 and pγ, ηq  S3), data from the Negative Binomial distribution with r " 3 (NB), and data from incorrectly specified π, where the data are generated using πpx; cq " cpsin x `1q (sin) or πpx; cq " cp2 ´|x|qIr|x| ă 2s (abs), while we test πpx; cq " cpcos x `1q.cost of motor vehicle collisions involve lost productivity, medical costs, legal and court costs, emergency services, insurance administration, travel delay, property damage and workplace losses, (Blincoe et al., 2015).Among all accidents, injuries to pedestrians and bicyclists caused 7 % of the total economic costs and 10 % of the total societal harm in US in 2010 (Blincoe et al., 2015).
In this application, we analyze daily numbers of accidents involving at least one pedestrian in Prague, the capital of the Czech Republic, in years 2018-2019.The data are available at the website of the Ministry of Interior of the Czech Republic, https://www.mvcr.cz.The considered time series has 730 observations and it consists of rather small counts, with minimum 0, maximum 8, and mean 1.70 collisions per day, see also Figure 1.Zero accidents are observed in 22 % of days.The plot of the sample autocorrelation function (ACF) in Figure 1 reveals an obvious weekly pattern.The 1-lag autocorrelation is rather small in magnitude, but it should not be ignored, as already shown in several previous analyses of traffic accident counts, see, e.g., Brijs et al. (2008), Pedeli and Karlis (2011).
It is well known that weather affects both the driver capabilities and road conditions via visibility impairments, precipitation, and temperature extremes.Hence, the following explanatory variables were considered in a model for number of collisions: • average daily temperature (AT) in Celsius degrees, • daily precipitation (P) in mm, • deterministic dummy covariate indicating a working day (W).
The corresponding weather variables were downloaded from the Czech Hydrometeorological Institute, https://www.chmi.cz/.The AT ranges from ´8.40 to 31.00 with mean 12.83.The daily precipitation takes values between 0.00 and 42.30 with mean 1.02, with 67% of values being 0.
During the model building the average daily temperature was transformed either to a categorical variable with 3 categories (AT in p5, 10s , p10, 25s, otherwise), considered linearly via a variable 35´AT, or different linear dependencies were considered for ATă 5 and ATě 5, modelled via variables IpATă 5q¨pAT`10q and IpATě 5q¨p35´AT).The transformations (shifts of AT) are performed, because only positive association can be captured by a PARX model and the covariates are assumed to take positive values.Furthermore, the precipitation was considered either linearly, or via an indicator IpPą 0q.Competitive models were fitted by MLE and compared using AIC and BIC, similarly as in Agosto et al. (2016).We considered AR(1), AR(2) and GARCH(1,1) structure, but the simplest AR(1) type model of structure `cJ X t seems to be the most appropriate, with X t including the external regressors.Note that in this application we model λ t as a function of the current value of X t .The bootstrap goodness-of-fit model was performed for B " 1000 with various tuning parameters pρ, γ, ηq.Following the results of the simulation study, we set ρ " 0, γ, η P t1{4, 1{2, 1u.
A block bootstrap was applied to the stochastic variables with block length l " 10 and l " 30, while the deterministic variables are held fixed.
The results of the goodness-of-fit tests for a Poisson model with λ t specified by M 1 and M 2 are presented in Table 4 together with the corresponding information criteria.The obtained p-values are non-significant for all choices of the tuning parameters pγ, ηq and both considered bootstrap block lengths l.AIC slightly prefers M 2 , while the opposite conclusion is obtained for BIC.As the p-values for the test are generally larger for M 1 and due to an easier interpretation, the simpler model M 1 is chosen as the final.The estimated parameters suggest that the difference in number of pedestrian collisions in working days and during weekends is in average equal to 1.As expected, precipitation (and consequently worse visibility) increases a risk of pedestrian collision, but the marginal effect is rather small in magnitude, and the same can be concluded for the marginal effect of temperature.

Discussion and conclusion
A goodness-of-fit test for count time series with exogenous covariates is proposed.The new test statistic is formulated in terms of Bieren's characterization and may be applied to general models of this type provided that the conditional PGF is specified.The limit distribution of this statistic is obtained under the null hypothesis of fixed, but otherwise arbitrary, order and rather general form of the link function.Asymptotic results for non-null situations are also presented, including consistency of the test against distributional misspecification as well as against certain types of violations of the hypothesized functional form of the conditional mean.The finite-sample performance of a bootstrap version of the suggested test is investigated by a Monte Carlo study against competing procedures, and corroborates the consistency features of our test.Finally the proposed methodology is used in order to perform explanatory analysis of the main factors affecting road accidents in the city of Prague.Weiss, C. H. and Schweer, S. (2015).Detecting overdispersion in INARCH(1) processes.

A Proofs
We start with some auxiliary technical lemmas.
Lemma A.1.Let tw t u be a sequence of random variables such that E |w t | κ ă 8 for some κ ą 0 and all t ě 1.Let ρ P p0, 1q.Then ρ t w t a.s.
Proof.The first statement follows from the Cantelli lemma, because by the Markov inequality and the assumption E |w t | κ ă 8.The second claim then follows from the property of the Cesàro limit.

■
In the following we often make use of the following simple relations, which hold for a, b, c, r a, r b, r c P IR: ab ´r a r b " apb ´r bq `r bpa ´r aq, abc ´r a r br c " abpc ´r cq `r cpab ´r a r bq " abpc ´r cq `ar cpb ´r bq `r br cpa ´r aq.
Proof.The assertions are proved in (Aknouche and Francq, 2021, Lemma A.1., Lemma A.2, Lemma A.3).Using the notation from ( 24), it follows that under our Assumption A3 it holds that r λθ " r θλ " Bdpθq{Bθ, which is a continuously differentiable function on Θ. Hence the supremum of its norm over the compact set Θ 0 is finite.In addition, r λλ " 0, so assumption A7 of Aknouche and Francq ( 2021) is satisfied.

Denote
which is the desired formula.
Similarly, let r Z t pθq be defined as where we now stress the dependence on θ.Then Recall the notation αpv 1 , v 2 q " cos `vJ 1 v 2 ˘`sin `vJ 1 v 2 ˘, and define where Also notice that assumption A7 implies W pu, vq " wpuqωpvq, with ż I R d sinpv J xqωpvqdv " 0 for all x P R d .Hence, we can get an equivalent expression for ∆ T,W as ∆ T,W " Thus, it suffices to investigate the asymptotic behavior of for some K 4 ă 8. Hence, the sequence tℓ 2,t u can be treated by the same arguments and it follows that 1 T Ñ E ℓ 0 2,t pu, v, θq " 0.
Therefore, the term A 2 behaves asymptotically as ?T p p θ ´θ0 q J βpu, vq, so in view of Assumption A6 we have, s t pθ 0 , Y t , I t´1 q J βpu, vq `oP p1q.

ˇˇˇ,
where as T Ñ 8, the first term converges to a finite number due to (16) in Assumption A5.Now let V η pθ 0 q be a neighborhood of θ 0 of radius 1{η.Then for η and T large enough V η pθ 0 q Ă Vpθ 0 q and the second term can be bounded by which is a sum of stationary and ergodic variables, so as T Ñ 8 this sum converges to the expectation Bλ t pθq Bθ i ´Bλ t pθ 0 q Bθ j Bλ t pθ 0 q Bθ i ˇˇw hich converges to 0 as η Ñ 8.The "coefficients" of T i p¨, ¨, ¨q , i " 3, 4, 6 in (31) can be treated in a completely analogous manner, while for the "coefficients" of T i p¨, ¨, ¨q i " 2, 5, in (31), we use Lemma A.5, otherwise the approximations are the same.
To sum up, we have that wpuqωpvqdudv converges to a finite distribution.In order to show this, let Z be the Gaussian process from the statement of the Theorem 4.1.It is clear that all the marginal distributions of Q " tQpu, vq, u P r0, 1s, v P IR d u converge to the marginal distributions of Z.We will show that for any compact F Ă IR d .To this end, and in view of Theorem 22 from Ibragimov and Chasminskij (1981) it suffices to verify that sup and that there exist constants C P p0, 8q and κ ą 0 such that sup for some finite constant H 2 .The condition (33) then follows from the properties of the weight function W . Condition (34): The stationarity of tδ t pu, vqu, and the the property of martingale differences imply that E " Q T pu 1 , v 1 q ´QT pu 2 , v 2 q ‰ 2 " E rδ 1 pu 1 , v 1 q ´δ1 pu 2 , v 2 qs 2 .
In the following let s t :" s t pθ 0 , Y where M 1 is from A4, and for some finite constant H 4 due to the Cauchy-Schwartz inequality and assumptions A5 and A2.Hence, the mean value theorem implies that To prove (34) we use again the Cauchy-Schwartz inequality to see that ď C 2 }pu 1 , v J 1 q J ´pu 2 , v J 2 q J }.
for some finite constant C 2 .

Figure 1 :
Figure 1: Daily numbers of pedestrian collisions (left panel) and the sample ACF (right panel).
zFε |Zpu, vq| 2 wpuqωpvqdudv ă ε.This finishes the proof.■ Proof of Theorem 4.2 p θq| 2 wpuqωpvqdvdu `oP p1q.Hence, the Taylor expansion is used directly on Bpu, v, θq and a careful inspection of the proof of Theorem 4.1 reveals under the assumptions of Theorem 4.2 it holds that Bpu, v, p θq " Qpu, vq `RT pu, vq, where R T pu, vq " o P p1q uniformly in u P r0, 1s and v P IR d .The rest of the proof then proceeds along the same lines.■ Proof of Theorem 4.3.Recall that we use gpu, λq for the PGF g λ puq of F λ and similarly, let g A pu, λq " g A λ puq stand for the PGF of F A λ .It is shown in the proof of Theorem 4Cpu, v, θq " T ´1{2 r Bpu, v, θq for r B given in (29).The Taylor expansion of r Cpu, v, p θq in θ 0 gives r Cpu, v, p θq " r Cpu, v, θ 0 q `B r Cpu, v, θq Bθ | θ"θ ˚pp θ ´θ0 q,where θ ˚is a point between θ 0 and p θ.A careful inspection of the proof of Theorem 4 Z t being an estimated version of Z t , computed from some initial values for Y t , p λ t , and X t for t ď 0, i.e. p Z t " pY t´1 , . . ., Y t´p , p λ t´1 , . . ., p λ t´q , X J t´1 q J for t ě p `1.Let W : IR d`1 Ñ p0, 8q be a nonnegative weight function.Define the test statistic ).Let p θ :" p θ T be a suitable estimator of the model parameter θ and define p ε t puq " u Yt ´gp λt puq, t " 1, ..., T, (12) where p λ t are for t " 1, . . ., T , defined recursively as p λ t " hp p Z t ; p θq, with p Assume that for the stationary variables Y t , λ t , X t there exists κ ą 0 such that

Table 4 :
p-values for the bootstrap goodness-of-fit test for Poisson models M 1 and M 2 for daily series of pedestrian collisions.