Bootstrap-based bias corrections for INAR count time series

ABSTRACT Integer-valued autoregressive (INAR) processes form a very useful class of processes suitable to model time series of counts. Several practically relevant estimators based on INAR data are known to be systematically biased away from their population values, e.g. sample autocovariances, sample autocorrelations, or the dispersion index. We propose to do bias correction for such estimators by using a recently proposed INAR-type bootstrap scheme that is tailor-made for INAR processes, and which has been proven to be asymptotically consistent under general conditions. This INAR bootstrap allows an implementation with and without parametrically specifying the innovations' distribution. To judge the potential of corresponding bias correction, we compare these bootstraps in simulations to several competitors that include the AR bootstrap and block bootstrap. Finally, we conclude with an illustrative data application.


Introduction
During the last decades, researchers paid more and more attention to discrete-valued time series, and especially count time series, which have the range N 0 = {0, 1, 2, . . .}, became a vivid research area [1]. For count time series with an autoregressive autocorrelation structure, i.e. where the autocorrelation function ρ(k) := Corr[X t , X t−k ] follows the well-known Yule-Walker equations for some autoregressive order p ∈ N = {1, 2, . . .} and corresponding autoregressive parameters α 1 , . . . , α p , a widely used approach is to consider the family of INAR(p) models according to [2]. These integer-valued autoregressive models, which contain the famous INAR(1) model by [3,4] as a special case for p = 1, are defined by a recursive scheme looking like the conventional AR(p) recursion, but with the multiplications being replaced by the binomial thinning operator of [5]: If X is a count random variable and if the thinning parameter satisfies α ∈ [0; 1], then the random variable α • X is said to arise from X by binomial thinning if α • X has a conditional binomial distribution given the value of X, i.e. α • X | X ∼ Bin(X, α). Using this probabilistic operator, Du and Li [2] define the INAR(p) model by the recursion where α • := p j=1 α j < 1. Here, the innovations ( t ) constitute an independent and identically distributed (i.i.d.) count process having some distribution G with range N 0 , with t being independent of (X s ) s<t and with the thinning operations being performed independently of each other and of ( t ). In particular, the thinnings (α 1 • t+1 X t , . . . , α p • t+p X t ) given X t are assumed to be conditionally independent.
For estimating the model parameters of an INAR(p) process, Yule-Walker (YW) estimation has been considered by [2,6]. Conditional least squares (CLS) estimation has been proposed by [2], and [7] investigated maximum likelihood (ML) estimation (which requires to specify the t 's distribution family). Besides these classical approaches, Drost et al. [8] proposed a semi-parametric maximum likelihood (SPML) estimator, which only relies on the INAR(p) model structure, but which, in contrast to the ML approach discussed by [7], does not require a distributional assumption with respect to the innovations t . In contrast to YW and CLS estimation, in turn, the SPML approach does not only estimate a couple of innovations' moments, but the full innovations' distribution (without assuming a parametric model). Such a distribution-free approach is very welcome in practice, since for model order p ≥ 2, it is rather difficult to identify an appropriate distribution for the INAR's innovations, because, e.g., in contrast to INAR(1), dispersion properties of the innovations do not carry over to the observations, see [9].
For the particular instance of a Poisson INAR(1) process, i.e. with autoregressive order p = 1 and with Poisson-distributed innovations ( t ), the performance of some estimators (or further statistics being computed from the given time series) has been investigated via simulations by several authors, e.g. [3,[10][11][12]; Bu et al. [7] considered the Poisson INAR(2) model in an analogous manner. Such estimators are known to be asymptotically consistent, but especially for small sample sizes, these authors reported about bias problems in estimation. For the specific Poisson INAR(1) model, it is possible to derive analytically asymptotic bias corrections for many relevant statistics [11,12]. However, the model complexity increases substantially if either the Poisson assumption or the limitation to model order p = 1 does not hold. Hence, it is not clear if such closed-form corrections could be derived also for differently distributed innovations or for higher model orders p, but this task clearly becomes a lot more involved. Therefore, in this work, we propose a general scheme for bias correction, which relies on the INAR bootstrap scheme recently developed by [13].
Bias correction for univariate AR processes has been put forward by a number of authors. Starting with [14], such approaches (with a focus on bias correction for the AR slope parameter) have been addressed by many authors during the last decades, see e.g. [15][16][17] for references. These include analytical as well as numerical or simulationbased approaches. In a multivariate context, analytical expressions for the finite-sample bias in estimated vector autoregressive (VAR) models have been developed by [18], see also [17] for references. The first to employ the (residual-based AR) bootstrap for bias correction for the VAR slope coefficient estimation was [15,16], see also [17,19] for discussions.
Whereas [18] derived closed-form expressions for the asymptotic first-order bias of the slope coefficients in the VAR model only without a deterministic time trend, the motivation of [15,16] to use the bootstrap for bias correction was to be able to deal easily with deterministic trends as well. A detailed simulation study for comparing different bias correction approaches in VAR models is provided by [17]. Similarly, in the context of INAR models, analytical expressions for the bias of most relevant statistics are available only for the special case of a Poisson INAR(1) model, so under order p = 1 and under a parametric assumption on the innovations' distribution, see Section 3.3 for further details. Hence, the idea to make use of distribution-free bootstrap methods appears to be rather natural to achieve robust bias correction. While most papers address only AR coefficient estimates, we look at different relevant estimators computed from INAR data.
In Section 2, we introduce the general INAR bootstrap framework considered in [13], before corresponding parametric and semi-parametric implementations are discussed. The latter approach crucially depends on the SPML estimator proposed by [8], which is addressed as well. In Section 3, we describe in detail how INAR bootstraps and competitors will be applied for bias correction. In Section 4, we investigate the finite sample performances of such bias corrections for several INAR DGPs. Two illustrative data examples are discussed in Section 5, and we conclude in Section 6. Theoretical details will be deferred to an appendix, and further simulation results to a supplement.

The INAR bootstrap
Recently, Jentsch and Weiß [13] proposed an INAR bootstrap scheme, which can be applied to statistics T n = T n (X), where X = (X 1 , . . . , X n ), belonging to the class of functions of generalized means, i.e. to statistics of the form where m ∈ N (fixed) and n m = n − m + 1, and where the functions g : R m → R d and f : R d → R have to satisfy some smoothness conditions. Note that the possible types of statistics include e.g. (functions of) the sample mean, sample autocovariances, and sample (partial) autocorrelations. The general INAR bootstrap scheme proposed by [13] applies to INAR(p) processes (1). Under some mild additional requirements, they were able to prove consistency of the INAR bootstrap. In our application of the INAR bootstrap scheme, we consider both a parametric implementation with G = G θ and corresponding estimatorĜ = Gθ as well as a semi-parametric implementation in the sense of [8]. The parametric implementation in general requires the proper choice and estimation of a parametric distribution (family), whereas the semi-parametric version comes without specifying the innovations distribution; see [13] for details. For the parametric case, the most prominent specification is a Poisson INAR model, i.e. G θ = Poi(θ ), which will be considered also below.
To implement the INAR bootstrap procedure, it is necessary to specify the estimation α 1 , . . . ,α p of the INAR coefficients as well as the estimationĜ of the innovations' distribution G, which can be done in multiple ways. One such implementation is the fully parametric approach according to [13]:
As an alternative, we also consider the semi-parametric approach in [13] to implement the INAR bootstrap procedure. It relies on the SPML estimator by [8] for INAR(p) models in (1), i.e. the semi-parametric joint estimation of the INAR coefficients α 1 , . . . , α p and the innovation distribution G: . . , n}. The conditional log-likelihood is computed by using the fact that the conditional distribution of X t given X t−1 , . . . , X t−p is a convolution of the binomial distributions Bin(X t−1 , α 1 ), . . . , Bin(X t−p , α p ) and of G; also see Appendix E.1 in [13] for computational issues.

INAR Bootstrap for bias correction
In [13], the finite-sample performance of both the parametric and the semi-parametric INAR bootstrap was investigated with respect to confidence intervals for several statistics T n (X) in the sense of (2). It turned out that both INAR bootstrap schemes (3) and (4) lead to satisfactory results, where, however, the performance of the parametric INAR bootstrap (3) deteriorated if the required parametric assumption (in this case, a Poisson assumption) was violated. The semi-parametric INAR bootstrap (4) was robust in these respects, and it also turned out to be superior to some further opponents: a naïve application of the conventional AR bootstrap (by just ignoring the discreteness of counts), a Markov bootstrap (which suffered severely from the huge number of parameters to be estimated), and a circular block bootstrap (the performance of which converged much more slowly with increasing sample size n). These positive results encourage us to apply the INAR bootstrap schemes for a completely different task, the bias correction of point estimators for some relevant statistics T n (X).
We provide the general approach to correct for a bias using the INAR bootstrap proposal, before we discuss the implementation for all statistics under consideration and, in particular, the corresponding centering schemes. Furthermore, we discuss several competitors based on analytical expressions and other resampling techniques.

Bootstrap-based bias correction
In general, let ϑ be a parameter of interest, which is estimated from the INAR(p) data X by the biased estimatorθ (computed according to an expression T n (X) in the sense of (2)). For the INAR-bootstrap bias correction, we apply either the parametric or the semi-parametric INAR bootstrap. Then, depending on the selected INAR bootstrap implementation, we compute an appropriate centering cent(θ * ), which we obtain from a functional relationship between the population quantity ϑ and the model parameters that are estimated for the bootstrap procedure itself (see details below). We execute the boot- Finally, we calculate the mean of the centred bootstrap estimates and define the bias-corrected estimate for ϑ bŷ
To be able to capture the finite sample bias ofθ, the use of an appropriate centering scheme is essential when computing the respective bias-corrected estimatesθ corr . These centering schemes are computed based on the respective model assumption, see Section 3.1. For example, the population mean of an INAR(p) model (1), which is estimated by the sample meanX, is known to satisfy the functional relationship This, in turn, leads to a centering forX * defined as cent(X * ) :=μ ; sp /(1 − α 1;sp − · · · −α p; sp ) (semi-parametric implementation); see details below. In the sequel, we provide expressions for centeringX * ,γ * (0) andγ * (1). Since the remaining estimators are functions thereof, we apply the corresponding function to cent(X * ), . . . to obtain an appropriate centering scheme, e.g. we centreρ * (1) : For the semi-parametric INAR bootstrap scheme (4), we assume an underlying INAR(p) process with specified model order p, but without further distributional assumptions. We obtain the estimated INAR coefficientsα 1;sp , . . . ,α p; sp and innovations' probabilitiesĝ sp (k) for k = 0, . . . , max X, and using the latter, we computeμ ; sp := max X k=0 k · g sp (k) andσ 2 ; sp := max X k=0 k 2 ·ĝ sp (k) −μ 2 ; sp . These are used to compute the following centerings: The corresponding functional relationships can be found, e.g. in [1,9]. Here and in Section 4, we restrict ourselves to model orders p ≤ 2, but centering for higher orders p could be computed by utilizing the respective Yule-Walker equations, see [2].
For the parametric INAR bootstrap scheme (3), we additionally assume a distribution family for the innovations ( t ). Here, we consider the Poisson INAR(p) model, i.e. we have t ∼ Poi(μ ) and thus equidispersed innovations, i.e. σ 2 /μ = 1. Note, however, that the observations become equidispersed (in fact even Poisson distributed) only for model order p = 1 [9]. Using the estimatorsX for the marginal mean μ andα 1;YW , . . . ,α p;YW (computed from the sample autocorrelation function according to the Yule-Walker equations of the INAR(p) model), we obtain the following centerings: • cent(X * ) :=X;

Possible competitors
In our simulation study presented in Section 4, we compare our above INAR-bootstrap bias corrections also to bias corrections obtained by (naïvely) applying the conventional AR bootstrap as well as the circular block bootstrap (using the corrected version of the automatic block-length selection procedure by [20,21]); for the corresponding centering schemes, we refer to Appendix E.2 in [13]. But in view of the bad experiences by the latter authors, we did not consider again the Markov bootstrap as a competitor. Furthermore, for model order p = 1 and assuming an underlying Poisson INAR(1) DGP, asymptotic approximations for the biases of the considered statisticsγ (0),Î disp ,γ (1), ρ(1) andX(1 −ρ(1)) are available in [11,12]. These are summarized in the centre column of Table 1. Plugging-in Yule-Walker estimates in these expressions instead of the true (but unknown) parameter values, a simple asymptotic bias correction is possible.
We also added the corresponding approximate bias corrections for an AR(1) DGP to Table 1 (right column; the derivation is sketched in Appendix), as this will help us in interpreting the simulation results below. Note that the distribution of the considered AR(1) model ist not specified, it is only assumed that relevant moments exist.

Performance of bias corrections
To evaluate the performance of the bias corrections presented in Section 3, a comprehensive simulation study was done, where mean and standard deviation of the (bias-corrected) estimates were computed. For each of the considered scenarios, 500 Monte-Carlo samples were simulated, and for each sample, in turn, 500 bootstrap replicates. For cases where the stationary distribution was not directly available, a pre-run of length 100 was used to approximately reach the steady state. Poisson innovations ('Poi-INAR') were parametrized as μ ∈ {1, 2.5}, and negative binomial innovations ('NB-INAR') as μ ∈ {1, 2.5} and σ 2 /μ ∈ {1.5, 2.0} (the parameters n, π of the NB(n, π)-distribution were computed from the relations σ 2 /μ = 1/π and μ = n((1 − π)/π)). The generated INAR (1) (1))) and 3 (concerning estimatorÎ disp ), any of the considered bias corrections leads to an improved mean, but with the CBB bootstrap being clearly worst. The latter is reasonable since the CBB bootstrap does not make use of any assumptions on the model structure. Since the DGP is Poi-INAR(1), and since both bias corrections 'pINAR' and 'asymp' rely on this Poisson assumption, it is clear that these corrections work very well in Tables 2 and 3. But notably, the spINAR-correction (which does not use the Poisson assumption) works equally well. Finally, let us look at the AR bootstrap. Although it has some corrective effect on estima-torÎ disp in Table 3, it is clearly worse than '(s)pINAR' and 'asymp', which is plausible as the DGP is not an AR process. In particular, we can see in Table 1 that the bias of the dispersion index based on AR data has a very different form in comparison to the bias when computed from INAR data. For the autocorrelation, the AR bootstrap is not consistent; see Appendix D in [13]. However, the bias expressions in Table 1 are very similar and differ only by a minor additional term for the INAR case, such that the AR bootstrap will work also well here. Furthermore, as the AR bootstrap will be generally consistent for the sample mean even for INAR data, it is not surprising that AR bootstrapping shows a good performance   in correcting the bias ofX (1 −ρ(1)) as well, see Table 2. Actually, the corresonding bias approximations in Table 1 are even identical. Besides the bias of an estimator, it is also important to look at the mean squared error (or the square-root thereof, hereafter abbreviated as 'RMSE'), because this measure is not only sensitive to the bias but also to the dispersion of the estimator. A reduction in bias often tends to be associated with an increase in the variance of an estimator. Hence, it will not be generally possible to reduce the (R)MSE as well. Looking again at Tables 2 and 3 (respective right part), we see that the RMSE is either also improved by using an INAR (or asymptotic) bias correction (Table 2), or it is only slightly increased as in Table 3. But compared to the strong improvement in bias, the small increase in RMSE appears to be acceptable for practice, which is also supported by the box plots shown in Figure 1: the uncorrected estimates ('no') look misplaced with respect to the true parameter value, but there is no visible difference in dispersion. The reason for a possible increase in RMSE is an increase in the standard deviation (SD) of the estimators, see Table 4.
Let us return to Table 3, where we concluded that the AR bootstrap performs worse than the INAR bootstraps because it corrects for a differing bias according to Table 1. But the picture is different in Table 5 about the variance estimatorγ (0). Note that the expression derived by [22] implies that asymptotically, all AR(1)-like DGPs lead to the same asymptotic bias forγ (0), namely −(1/n)((1 + α)/(1 − α))σ 2 (with σ 2 = μ in the Poisson case), also see Table 1. Hence, although we see some deviation for small sample size n = 100, INAR and AR bootstraps tend to be very similar for larger n in correcting the bias forγ (0).    But even more interesting is a comparison of both parts of Table 5: in the left part, the DGP is Poi-INAR(1) as assumed by 'pINAR' and 'asymp', whereas this is violated in the right part about an NB-INAR(1) DGP. As a result, these bias corrections perform clearly worse for the NB-INAR(1) DGP, while 'spINAR' still successfully corrects the bias and is clearly best among all approaches. Next, we turn our attention to autoregressive order p = 2, i.e. the DGP is an INAR(2) process. Nevertheless, it may happen that we falsely assume order p = 1 when applying our techniques for bias correction. The result of such a model underfitting becomes clear from Table 6, where the dispersion index is estimated by either falsely assuming p = 1, or correctly assuming p = 2. Note that the considered Poi-INAR(2) process is marginally overdispersed despite its Poisson innovations. Furthermore, if assuming model order p > 1, it makes no sense anymore to consider the asymptotic correction as this relies on order p = 1. Looking at Table 6, it becomes clear that both INAR bootstrap schemes suffer from an underfitted model (although they are still superior to not doing a correction), i.e. for practice, it is important to not choose the model order of the INAR(p) DGP too small. If the bootstrap schemes correctly assume p = 2, it can again be observed that the AR bootstrap is worse than the INAR bootstraps, for the reasons discussed above.   Finally, we consider the estimation of the first-order autocorrelation, both for the Poi-INAR(1) and Poi-INAR(2) DGP. Tables 7 and 8 (which refer to the medium autocorrelation level ρ(1) = 0.4) show that any of the bias corrections leads to an improvement compared to the case without correction, namely a reduced bias without increasing the RMSE. The INAR and AR bootstrap schemes perform similarly well, which is reasonable in view of our above discussion about the asymptotic bias formulae in Table 1. Comparing Tables 7 and 8, we see that the corrected estimates for ρ (1) show less bias and less RMSE for the lower model order p = 1, which is reasonable since we have one parameter less to be estimated. Analogous observations can be made for the highly correlated case ρ(1) = 0.8, which is illustrated by box plots in Figure 2.

Illustrative data examples
We illustrate the considered bias corrections with two data examples taken from the book by [1]. The first time series (length n = 267) gives daily counts of downloads of a TEX-editor and was shown in Section 2.5 of [1] to be of INAR(1)-type with a strong degree of overdispersion. The second time series (length n = 370) consisting of counts of gold particles over time, in contrast, appears to be of INAR(2)-type with a nearly equidispersed marginal distribution, see Example 3.1.6 in [1]. This can also be seen from Table 9, which shows some results from SPML estimation according to [8]: for the downloads counts, the estimate for the second autoregressive parameter deviates only slightly from 0 (thus also supporting a first-order autoregression), and the innovations' variance is much larger than the mean (within the fitted model), whereas the gold particles counts have nearly equidispersed innovations (as it would be the case for a Poisson INAR model) and a clearly non-zero value forα sp,2 . Certainly, we do not know the true data-generating mechanism behind these time series, but these observations will help us in interpreting the subsequent findings.
Let us concentrate on the estimation of the varianceγ (0), the dispersion indexÎ disp , and the first-order autocorrelationρ(1) in the sequel, which are meaningful and relevant quantities for both model orders p = 1 and p = 2. For bias correction, we applied the bootstrap schemes presented in Section 3 (now we managed even 10 000 bootstrap replications), where the optimal block length for the block bootstrap [20,21] was calculated as 6 (downloads) and 33 (gold particles), respectively. For model order p = 1, we also applied the asymptotic bias corrections described in Section 3.3.
Results for the downloads counts are shown in Table 10. It becomes clear that all methods for bias correction have a similar effect on the autocorrelationρ(1), which is not surprising in view of the discussions in Section 4. It is also plausible that there is nearly no difference between model orders p = 1 and 2, since the first-order model is included as a boundary case within the second-order one (overfitting increases computational costs, whereas underfitting detoriates the performance of bias correction, remember the discussion of Table 6).  Looking at the varianceγ (0) and dispersion indexÎ disp , it is apparent that the methods 'pINAR' and 'asymp' have a smaller corrective effect than 'spINAR' and 'CBB'. This is plausible since both rely on a parametric Poisson assumption, which appears to be severly violated by the data.
Let us now discuss the results for the gold particles counts, see Table 11. For these data, there is strong evidence for a second-order autoregressive structure, also see Example 3.1.6 in [1]. Hence, it is not surprising that the INAR bootstrap schemes show visible deviations if being based on model order p = 1 or 2. The difference between 'spINAR' and 'pINAR', however, is rather small, which can be explained by the fact that the gold particles counts are close to equidispersion, i.e. a Poisson INAR(2) model appears to be a reasonable choice for the data.

Conclusions
We propose to use INAR-type bootstrap procedures for finite sample bias correction of several practically relevant statistics computed from INAR data. In general, we found that bias correction works quite well in simulations. We compare the potential of INAR bootstrap bias correction to other competitors including bias correction exploiting asymptotic approximations as well as other resampling techniques as AR bootstrap and block bootstrap. The INAR bootstraps tend to outperform their competitors in simulations. In particular, the block bootstrap converges more slowly and performs worse for small sample sizes. The AR bootstrap shows good performance for some statistics, but is clearly worse for others. We were able to explain these phenomena by making use of asymptotic approximations of the bias in special cases. Furthermore, we found that the parametric INAR bootstrap and the asymptotic bias correction, which both rely on Poisson assumptions, performed worse whenever these assumptions were violated. In such cases, the semiparametric INAR bootstrap were robust to non-Poisson innovations and still led to good results.

Disclosure statement
No potential conflict of interest was reported by the authors.