Semiparametric estimation of INAR models using roughness penalization

Popular models for time series of count data are integer-valued autoregressive (INAR) models, for which the literature mainly deals with parametric estimation. In this regard, a semiparametric estimation approach is a remarkable exception which allows for estimation of the INAR models without any parametric assumption on the innovation distribution. However, for small sample sizes, the estimation performance of this semiparametric estimation approach may be inferior. Therefore, to improve the estimation accuracy, we propose a penalized version of the semiparametric estimation approach, which exploits the fact that the innovation distribution is often considered to be smooth, i.e. two consecutive entries of the PMF differ only slightly from each other. This is the case, for example, in the frequently used INAR models with Poisson, negative binomially or geometrically distributed innovations. For the data-driven selection of the penalization parameter, we propose two algorithms and evaluate their performance. In Monte Carlo simulations, we illustrate the superiority of the proposed penalized estimation approach and argue that a combination of penalized and unpenalized estimation approaches results in overall best INAR model fits.


Introduction
According to Du and Li (1991), the INAR(p) model is defined by the recursion X t ¼ a 1 X tÀ1 þ a 2 X tÀ2 þ . . . þ a p X tÀp þ e t ; t 2 Z; with innovation process e t $ i.i.d G; where the distribution G has range N 0 ¼ f0; 1; 2; . . .g. Furthermore, let a ¼ ða 1 ; . . .; a p Þ 0 2 ð0; 1Þ p denote the vector of model coefficients with P p i¼1 a i \1 and where "" is the binomial thinning operator first introduced by Steutel and Van Harn (1979). Here, Z ðt;iÞ j ; j 2 N; t 2 Z ; i 2 1; . . .; p, are mutually independent Bernoulli-distributed random variables Z ðt;iÞ j $ Binð1; a i Þ with PðZ ðt;iÞ j ¼ 1Þ ¼ a i independent of ðe t ; t 2 ZÞ. The special case p ¼ 1 results in the INAR(1) model introduced by McKenzie (1985) and Al-Osh and Alzaid (1987). All the thinning operations "" are independent of each other and of e t ; t 2 Z. Furthermore, the thinning operation at time t and e t are independent of X s ; s\t.
Most researchers deal with parametric estimation of INAR models (see for example Franke and Seligmann (1993), Freeland and McCabe (2005), Brännäs and Hellström (2001) and Jung et al. (2005)), i.e. they assume G to lie in some parametric class of distributions ðG h j h 2 H & R q Þ for some finite q 2 N. In contrast, Drost et al. (2009) introduced a semiparametric estimator, which on the one hand keeps the parametric assumption of the binomial thinning operation, but on the other hand allows to estimate the innovation distribution nonparametrically. Using empirical process theory, they derive asymptotic theory in terms of consistency and asymptotic normality results and proved efficiency. Consequently, their estimation approach does not require any parametric assumption regarding the innovation distribution, and avoids the risk of a falsely specified parametric assumption and its undesirable consequences. The approach estimates the coefficients of INAR models and the innovation distribution simultaneously. The resulting semiparametric maximum likelihood estimator ðâ sp ;Ĝ sp Þ ¼ ðâ sp;1 ; . . .;â sp;p ;Ĝ sp ð0Þ;Ĝ sp ð1Þ;Ĝ sp ð2Þ; . . .Þ; whereâ sp ¼ ðâ sp;1 ; . . .;â sp;p Þ denotes the vector of the estimated INAR coefficients and fĜ sp ðkÞ; k 2 N 0 g are the estimated entries of the probability mass function (PMF) of G, maximizes the conditional log-likelihood function logðLða; GÞÞ, i.e. 8n 2 Z þ : ðâ sp ;Ĝ sp Þ 2 arg max ða;GÞ2½0;1 p ÂG Y n t¼0 P a;G ðX tÀ1 ;...;X tÀp Þ;X t ! : ð2Þ Here,G is the set of all probability measures on Z þ and P a;G ðX tÀ1 ;...;X tÀp Þ;X t are the transition probabilities under the true model parameters a and G, i.e. P a;G ðx tÀ1 ;...;x tÀp Þ;x t ¼ P a;G X p i¼1 a i X tÀi þ e t ¼ x t j X tÀ1 ¼ x tÀ1 ; . . .; X tÀp ¼ x tÀp ! ¼ ðBinðx tÀ1 ; a 1 Þ Ã . . . Ã Binðx tÀp ; a p Þ Ã GÞfx t g; with P the underlying probability measure and "Ã" denoting the convolution of distributions. In the special case of an INAR(1) model the transition probabilities are given by x tÀ1 j a j ð1 À aÞ x tÀ1 Àj P a;G ðe t ¼ x t À jÞ; where a is the coefficient of the INAR(1) model (McKenzie 1985;Al-Osh and Alzaid 1987). For k\minfX t À P p i¼1 X tÀi j t ¼ p þ 1; . . .; ng or k [ maxfX t j t ¼ 1; . . .; ng, the valuesĜ sp ðkÞ; k 2 N 0 ; are equal to 0. For further details, see Drost et al. (2009). In practice, discrete probability distributions such as the Poisson, the negative binomial or the geometric distribution are often used as innovation distribution G, see Weiß (2018), Yang (2019), Al-Osh and Alzaid (1987), Al-Osh and Alzaid (1990). The common feature of all these distributions is their smoothness in the sense that consecutive entries of their PMFs differ only slightly from each other. However, for a small sample size n, the semiparametric estimation approach of Drost et al. (2009) may lead to rather non-smooth innovation distributions with unnatural gaps in their PMF. For illustration, we consider a time series containing counts of transactions of structured products (factor long certificates with leverage) from on-market and offmarket trading per trading day between February 1, 2017 and July 31, 2018 (thus n ¼ 381). These data, which are plotted in Fig. 1, have first been presented by Homburg et al. (2021), who derived them from the Cascade-Turnoverdata of the Deutsche Börse Group. In the upper right corner, we see the estimated innovation distribution using the semiparametric procedure of Drost et al. (2009) which turns out to be smooth. In the second row, we consider only the first 100 observations of the time series, where the first plot shows indeed a bimodal estimated innovation distribution. In the third row, we only considered the first 20 observations. The lowerleft plot shows the resulting estimated PMF, which contains an unnatural gap witĥ G sp ð3Þ being estimated exactly equal to zero while its neighborsĜ sp ð2Þ andĜ sp ð4Þ are estimated positive. Hence, the resulting estimation is not smooth contrary to the estimated innovation distribution on the whole time series. In general, such nonsmooth innovation distributions are not common in practice and instead, smoothly estimated innovation distributions are often desired. In this paper, we want to use this prior knowledge and take advantage of a natural qualitative smoothness assumption on the innovation distribution by proposing a version of the semiparametric estimation approach, which penalizes the roughness of the innovation distribution. The resulting estimated PMFs of this approach are contained in the right plots in the second and third row, respectively. In comparison, the penalized estimation now leads to a smoother estimation of the PMF without any gaps. We will have a closer look at additional real data examples in Sect. 4. For long time series, the smoothing caused by penalization is not of such great importance, because the distribution estimated without penalization will be sufficiently smooth by itself. But for short time series, estimation without smoothing will commonly lead to jagged estimated innovation distributions although the true distribution behind the data might be smooth. So the need for smoothing is of particular importance for short time series.
The paper is organized as follows. In Sect. 2, we introduce a penalized estimation approach using roughness penalization and propose two algorithms for the datadriven selection of the penalization parameter. Section 3 examines our estimation approach in a comprehensive simulation study, where we compare the estimation performance of the penalized and the unpenalized approach for different settings. In a real data application in Sect. 4, we analyze the monthly demand of car spare parts to illustrate our method and its practical relevance. In the conclusion in Sect. 5, we summarize the results and give an outlook on further research questions. Penalized estimation of count data is a modern topic in current statistical research. Bui et al. (2021) consider parameter estimation in count data models using penalized likelihood methods. In a time series context, Nardi and Rinaldo (2011) studied LASSO penalization for fitting autoregressive time series models to get sparse solutions, i.e. where some autoregressive coefficients are estimated exactly as zero.
Fokianos (2010) proposed an alternative estimation scheme for the estimation of INAR models based on minimizing the least square criterion under ridge type of constraints. Wang (2020) proposed a variable selection procedure for INAR(1) models with Poisson distributed innovations including covariables by using penalized estimation and Wang et al. (2021) introduced an order selection procedure for INAR(p) and INARCH(p) models also by using penalized estimation. By contrast, in this paper, we propose a penalized estimation approach for INAR models which does not rely on a penalization of the INAR coefficients (towards zero), but on a penalization of the roughness of the innovation distribution (towards smoothness).

Penalized estimation approach using roughness penalty
The idea of our approach is to penalize the log-likelihood used in the semiparametric estimation of the INAR model according to Drost et al. (2009). Thus, we still do not assume a parametric class of distributions, we only use the assumed qualitative (i. e. nonparametric) property of smoothness. More precisely, this refers to a roughness penalization as introduced by Scott et al. (1980), which is e.g. used by Adam et al. (2019) for developing a nonparametric approach to fit hidden Markov models to time series of counts. We design the penalty term based on the idea of Tibshirani et al. (2005), where differences of successive parameters are penalized. In this regard, we allow for differences of order m 2 N. Applied to our setting, the estimation approach based on Drost et al. (2009) now maximizes the penalized log-likelihood (compare (2)) logðL pen ða; GÞÞ ¼ logðLða; GÞÞ À g Á d G;m ; where g [ 0 is the so-called smoothing or penalization parameter, d G;m denotes a suitable measure to quantify the roughness of G and m corresponds to the order of difference. According to Tibshirani et al. (2005), a first possible roughness measure for the penalization term is based on the L 1 distance (LASSO penalization), i.e.
The idea behind choosing this second roughness measure is that it does not shrink the differences of the successive entries of the PMF exactly to 0 (contrary to the first roughness measure), but the differences become close to 0, which is more in line with the idea of a smooth distribution (note the analogy of penalized regression, where the L 1 penalization is used for variable selection because of this property, see Fahrmeir et al. (2013)). The order of the differences m is a tuning parameter. For m ¼ 1, we penalize only the distance between two directly consecutive entries, for m ¼ 2 the smoothness is extended to a triple of values, etc.
Remark 1 A possible extension would be to allow for different penalization weights ðg i Þ for the individual (higher-order) differences of the entries of the PMF. For instance, in the case of L 1 penalization, the goal could be to maximize logðLða; GÞÞ À X maxðx 1 ;...;x n Þ i¼m g i j D m GðiÞ j; analogously for the case of L 2 penalization. Figure 2 shows a first exemplary result on a sample of an INAR(1) process with n ¼ 25 observations, order of difference m ¼ 1 and smoothing parameter g ¼ 1 roughly chosen by eye. In this example, the benefit of penalization already becomes clear. The penalized estimated innovation distributions are much closer to the true Poi(1) innovation distribution (which was truncated at value six for clarity) than the unpenalized estimated innovation distribution. Also, the difference between the L 1 and the L 2 penalization becomes visible. When using the L 2 penalization, the distances between the values of the PMF become small, when using the L 1 penalization they are shrinked to zero.

Selection of the penalization parameter
Now, we propose two approaches to determine for a fixed roughness measure the optimal smoothing/penalization parameter g, which is a trade-off between fit to the data and the smoothness assumption. For this purpose, we adapt as a first approach the cross-validation procedure described in Adam et al. (2019) to our setting. Therefore, we split the data set into s blocks F i ; i ¼ 1; . . .; s, of roughly equal size. In each fold i, F ðÀiÞ denotes the in-sample data (data without F i ) and F i the out-ofsample data. This replicates the correct dependence structure except for the "glue points", which only has a minor effect in practice when the data originate from an INAR model of small order. The greedy search algorithm is structured as follows: Algorithm 1 (1) Choose an initial g ð0Þ [ 0 and set z ¼ 0.
(2) For each fold i and for each value on a specified grid f. . .; g ðzÞ À 2c; g ðzÞ À c; g ðzÞ ; g ðzÞ þ c; g ðzÞ þ 2c; . . .g where c 2 R is a small constant, estimate the model with penalization on F ðÀiÞ and compute the penalized log-likelihood on F i .
(3) Average the resulting log-likelihood values across all folds i and choose g ðzþ1Þ as the penalization parameter on the grid that yields the maximum value. (4) Repeat steps 2) and 3) until g ðzþ1Þ ¼ g ðzÞ and define g opt :¼ g ðzþ1Þ .
Furthermore, to avoid a potentially non-optimal selection of the penalization parameter g caused by an inappropriate choice of the initial value g ð0Þ , we propose a second optimization algorithm. How we split the data in each fold j; j ¼ 1; . . .;s, in in-and out-of-sample data is specified later in Sect. 3.

Algorithm 2
(1) For each fold j and each value g on a specified grid f0;c; 2c; 3c; . . .; ug on the interval [0,u] for an appropriate upper bound u, estimate the model with penalization on the in-sample data and compute the penalized log-likelihood on the out-of-sample data.
(2) Average the resulting log-likelihood values across all folds j.
(3) Fit a polynomial of order r to the curve resulting from plotting the average outof-sample log-likelihood against the grid. (4) Choose g opt as the value on the grid, where the curve takes its maximum value.

Simulation study
We investigate the performance of the proposed procedure in a simulation study with K ¼ 500 Monte Carlo samples of size n 2 f20; 50; 100; 250; 500; 1000g generated from an INAR(1) process according to (1) for p ¼ 1 with different coefficients a 2 f0:2; 0:5; 0:8g and innovation distributions G 2 fPoið1Þ; NB 2; 2 3 À Á ; Geo 1 2 À Á ; ZIP 1 2 ; 2 À Á g, where ZIP denotes a zero-inflated Poisson distribution as in Jazi et al. (2012). The parameters of the negative binomial, geometric and zero-inflated Poisson distribution are hereby chosen to have the same expected value as the Poið1Þ distribution. But contrary to the Poið1Þ distribution which is equidispersed, i.e. the variance of the distribution equals its mean, they are overdispersed, i.e. their variances are larger than their mean values. Another difference between the considered innovation distributions is their (non-) smoothness, see also Fig. 12 in the appendix. The Poi(1), NB 2; 2 3 À Á and Geo 1 2 À Á distributions are rather smooth, but the ZIP 1 2 ; 2 À Á distribution, which shows a pronounced zero probability, is not. The effect of this property on the roughness penalization is investigated in Subsect. 3.5. Moreover, in Subsect. 3.2, we also provide a small simulation setting for higher-order INAR processes and consider the case of an INAR(2) model. The implementation is straightforward but is a lot more demanding such that we restrict the considered setting to a rather small extent. To ensure the stationarity of the time series, we actually generate n þ 100 observations and remove the first 100 observations. We consider first (m ¼ 1) and second (m ¼ 2) order differences in the penalization term (see Subsect. 3.4). As initialization for the smoothing parameter g ð0Þ , we set g ð0Þ ¼ 1 as in the example in Fig. 2 for the sample sizes n 2 f20; 50; 100; 250g and for computing time reasons g ð0Þ ¼ 0:5 for n 2 f500; 1000g. 1 For the considered grid around the smoothing parameter (see Algorithm 1) we choose c ¼ 0:05 resulting in fg ðzÞ À 0:1; g ðzÞ À 0:05; g ðzÞ ; g ðzÞ þ 0:05; g ðzÞ þ 0:1g. Unless stated otherwise, we use a ¼ 0:5 as true INAR(1) coefficient and Algorithm 1 with 10-fold cross validation (s ¼ 10) as optimization algorithm. For the realization of the simulation study, we use the statistical programming language R 4.1.2 (R Core Team 2021). for the different sample sizes and the respective estimation methods (unpenalized (up), L 1 penalization and L 2 penalization) for some large enough M. We use M ¼ 70 as upper bound for the observations x 1 ; . . .; x n since after this value the corresponding probabilities of occurrence are negligibly small. When the sample size n is small, the penalized estimation of the innovation distribution provides a large benefit compared to the unpenalized estimation: The L 2 distances of the penalized estimated to the true innovation distribution are much smaller than those of the unpenalized estimated to the true innovation distribution. Furthermore, the L 2 penalization performs better than the L 1 penalization. In Table 3 in the appendix, we also report the variance, the bias and the MSE of the first five estimated entries of the PMF resulting from the different procedures for the different sample sizes n. We see that the penalized estimation reduces both the variance, the absolute bias and consequently also the MSE of the estimated innovation distribution, especially for small n. Figures 17 and  18 and Tables 6 and 7 in the appendix show the analog results for a true NB 2; 2 3 À Á and Geo 1 2 À Á distribution, respectively. In general, regardless of the distribution and up to a sample size of n ¼ 100, we see a clear improvement concerning the estimation performance when using penalization. From a sample size of n ¼ 250 on, this improvement can only be seen marginally with the different methods essentially coinciding for large n. In Fig. 13 and Table 4 in the appendix, we show the results for INAR coefficient a ¼ 0:2 and Poi(1) innovation distribution and, correspondingly, in Fig. 15 and Table 5 for a ¼ 0:8. In the latter case, the benefit of the penalized estimation compared to the unpenalized estimation is even larger than in the case a ¼ 0:5. This is plausible because it is in general more difficult to estimate the innovation distribution for a larger value of a as this also leads to a larger observations mean with innovations mean remaining constant. Therefore, more entries of the PMF have to be estimated with the same amount of data. Contrary, for a ¼ 0:2, we have (with analog arguments) less entries of the PMF which have to be estimated with the same amount of data, which simplifies the estimation of the PMF in general and the benefit of penalization decreases. Altogether, we can conclude that the benefit of penalization is more pronounced with larger a, that is, with larger serial dependency.

Roughness penalty for smooth innovations distributions and first order differences
We get confirming conclusions, when we consider the values of the optimal smoothing parameter g, which approaches zero with increasing n, see  (1) innovation distribution of an INAR(1) process for different sample sizes n. We report results for unpenalized (up), L 1 and L 2 penalized estimation case of a true Poi(1) innovation distribution, Figs. 19 and 20 in the appendix for the cases of a true NB 2; 2 3 À Á and Geo 1 2 À Á innovation distribution and Figs. 14 and 16 in the appendix in case of a true Poi(1) innovation distribution with a ¼ 0:2 and a ¼ 0:8, respectively. Thus, for increasing n, the penalized and the unpenalized estimation coincide as intuitively expected: For large n, there are enough observations to learn the smoothness of the innovation distribution from the data even without imposing smoothness through penalization.

Higher-order INAR processes
To show that our proposed procedure is also applicable for higher-order INAR processes, we consider the case of a true INAR(2) process according to (1) for p ¼ 2 with coefficients a 1 ¼ 0:3, a 2 ¼ 0:2 and G ¼ Poið1Þ. Due to the high computing time for the semiparametric estimation, we only consider a small simulation setup with n ¼ 50 observations and K ¼ 100 Monte Carlo samples. We consider L 1 and L 2 penalization with first order differences and compare the performance with the case of estimation without penalization. In Fig. 21 in the appendix, we see that also for higher-order INAR models, penalized estimation of the innovation distribution provides a clear benefit compared to unpenalized estimation. With penalization we are closer to the true innovation distribution than without and we are able to reduce the variance, the absolute bias and consequently the MSE of our estimation, see Table 1. Again, L 2 penalization works best.

Alternative selection of the penalization parameter
To investigate whether the results depend on the chosen initial parameter, we now determine the optimal penalization parameter alternatively using Algorithm 2 with u ¼ 5;c ¼ 0:1 and r ¼ 5. In this context, we want to address a potential practical issue of Algorithm 1: the generation of the in-and out-of-sample data. For each of the 10 folds, 90% of the data becomes the in-sample data and the remaining 10% the out-sample data. For small n, 10% of the data is small. To avoid this, we now use an n-fold cross-validation (s ¼ n) for sample sizes n 2 f20; 50g with Algorithm 2, where starting from each observation the following 50% of the data is in-and the other 50% is out-of-sample. When reaching the end of the time series, we start again from its beginning. In Fig. 5, we see the results of this alternative procedure compared to the previous (iterative) procedure in Algorithm 1. It gives slightly better results than the iterative method, but overall the distances are very similar. The same can be concluded when considering Table 8. The alternative procedure leads to slightly lower MSE values, but altogether the values resemble each other. The 10-fold cross-validation also seems to be suitable and the resulting optimal parameters of the two procedures are close to each other (see Fig. 6). In conclusion, if we determine the optimal parameter from a sequence on a grid as in Algorithm 2, we tend to get slightly better results. However, the price to pay is a much higher computing time than with the iterative procedure. The iterative method needs a reasonably chosen starting value, but then it gives similarly good results in considerably less computing time. In addition, when using the alternative method, the question arises how to choose the upper limit of the interval adequately. In the following, we will continue to use the iterative method from Algorithm 1 but one should keep in mind that Algorithm 2 is also a practically useful procedure.

Higher-order differences in penalization term
So far we only considered first order differences (m ¼ 1). Now we want to see if penalizing higher-order differences (e.g. m ¼ 2) is able to improve the performance of our penalized estimation method. In Fig. 7 and Table 10 in the appendix it is visible for the case of a true Poi(1) innovation distribution and L 2 penalization that also the penalization of differences of higher order performs better than the unpenalized estimation in the cases of small sample sizes, and that it comes close to the penalization of first order differences. Similar results are in Fig. 22 and Table 9, both in the appendix, where we see the results of first and second order differences for the L 1 penalization. In case of L 1 penalization, we would prefer second order differences for small sample sizes. Overall, however, the L 2 penalization of first-order differences performs best.

Non-smooth innovation distribution
Finally, let us consider the case of ZIP 1 2 ; 2 À Á distributed innovations and consider Fig. 8. The results are as expected. Since the ZIP distribution is not smooth (see Fig. 12 in the appendix), the smoothness assumption and hence the penalization is not suitable. The boxplots reflect this: Except for sample size n ¼ 20, the penalized estimation procedure provides no benefit and for some n even leads to slightly higher L 2 distances from the true ZIP 1 2 ; 2 À Á distribution than the unpenalized procedure. As we can see in Table 11, the penalized estimation leads to a higher absolute bias when estimating the first (non-smooth) entry, G(0), of the PMF. As sample size n increases, the penalization has less impact, as there is enough data to detect the incorrect assumption such that the unpenalized and the penalized procedures coincide.
For comparison, let's take a look at the results with a true ZIP 1 2 ; 2 À Á distribution when we exclude G(0) from the penalization displayed in Fig. 9, i.e. when we considerd  (1) process for the different sample sizes n. We report results for unpenalized (up) and L 2 penalized estimation using either first order (diff1) or second order (diff2) differences expect: By excluding the "non-smooth entry" G(0) of the PMF of the innovation distribution from penalization, the penalized estimation works well again and provides a benefit for small n. In this case, the penalized estimation now results in a lower absolute bias of the estimated PMF's first entry compared to the unpenalized estimation (compare Table 12). However, this benefit is not as pronounced as in the cases of a true Poi(1), NB 2; 2 3 À Á and Geo 1 2 À Á innovation distribution. This can probably be explained by the fact that the ZIP 1 2 ; 2 À Á distribution has most of its mass in zero and the corresponding entry of the PMF, G(0), remains unaffected by the penalization. Consequently, the results from penalized and unpenalized estimation do not differ substantially from each other. In summary, if the smoothness assumption of the innovation distribution is correctly imposed, it provides a large benefit for small sample size n. This holds whether the true underlying distribution is equidispersed or overdispersed. The best results are obtained for L 2 penalization and first-order differences.

Estimation of the INAR coefficient
A drawback of the penalized estimation is that the estimation of the INAR coefficient a no longer works well for small sample size n, see Fig. 23 in the appendix. A strength of the semiparametric estimation approach of Drost et al. (2009) is the accurate joint estimation of the INAR coefficient and the innovation distribution. This joint estimation accuracy is not maintained when penalization is used for small n. The L 2 distances of the penalized estimated INAR coefficient a to the true value are higher than for the unpenalized estimated coefficient. For increasing n, the estimation of a improves, but since the benefit of the penalized estimation lies in the cases where n is small, this is no comfort.
Instead, we can solve this problem by taking only the estimator for the innovation distribution from the penalized approach and estimating the INAR coefficient with the unpenalized (efficient) estimation approach of Drost et al. (2009). In Fig. 23, we see that it is indeed preferable to combine the unpenalized estimation of the INAR coefficient a and the penalized estimation of the innovation distribution G. Also when looking at the MSE, it is clear that this combination outperforms all other estimation approaches under consideration.

Real data example
For modeling intermittent demand, Syntetos and Boylan (2021) consider the equidispersed Poisson distribution on the one hand, and, as the demand variability may be severe when demand is intermittent, overdispersed distributions from the Compound-Poisson family (such as the negative binomial distribution) on the other hand. All these parametric distributions are smooth. With our novel penalized semiparametric estimation approach, we get smooth distributions without parametric assumptions, and as we saw in our simulations, our penalization procedure works well for both equi-and overdispersed distributions. By contrast, if using an unpenalized non-parametric estimation approach such as the empirical distribution function (EDF), Syntetos and Boylan (2021) criticize that demand values not observed in the past are automatically assigned zero probabilities for the future. Furthermore, they state that an EDF provides a perfect fit to the historical data, but it does not ensure the goodness of fit to the demand over the forecast horizon, especially with respect to higher percentiles. Again, these drawbacks are omitted with our penalized estimation approach. Finally, historical demand time series are often rather short, see the demand count time series provided by Snyder (2002) as an example, such that smoothing approaches would be particularly welcome. For these reason, the forecasting of intermittent demand appears to be a promising application area for our proposed penalized semiparametric estimation procedure.
Therefore, we consider time series (n ¼ 51) of the monthly demand of different car spare parts offered by an Australian subsidiary of a Japanese car company from January 1998 to March 2002 (Snyder 2002). Figure 10 contains an exemplary time series of car part 2404. The observations vary between 0 and 5 and the up and down movements indicate a moderate autocorrelation level. After inspecting the corresponding (P)ACF also included in Fig. 10, we conclude that an AR(1)-like model might be appropriate for describing the serial dependence of the time series. Moreover, L 2 penalization with first order differences leads to an estimated innovation distribution without any unnatural gaps, i.e. zero values, in the PMF. Now consider the 1-step median prediction and the 90% quantile of the 1-step prediction of the demand for car spare part 2404. The latter serves here as a worstcase scenario for spare parts requirements. Therefore, we determine the median and the 90% quantile of the predictive distribution Pð. . . j yÞ, where y 2 0; . . .; 10. Based on the results of the simulation study in Subsect. 3.6, we use the penalized estimated innovation distribution and the unpenalized estimated INAR coefficient to determine the conditional predictive distribution. Table 2 shows that the penalized estimation tends to lead to higher predicted values (more conservative prediction). Consequently, without penalizing the innovation distribution, the predictions for the demand for spare parts may be too low, which can lead to a lack of spare parts. Moreover, the penalization of the innovation distribution (especially for such short In addition, we consider car spare part 1971. Figure 11 again suggests an AR(1)like model and a moderate autocorrelation level. The observations vary between 0 and 4 and there may be zero inflation in this time series. Therefore, in addition to the unpenalized and penalized estimates, we also consider the penalized estimate of the innovation distribution, where G(0) is not smoothed (see Subsect. 3.5). It becomes clear that this last estimation procedure yields more plausible results than when G(0) is smoothed. Again, the penalized estimation procedure yields a slightly smoother innovation distribution than the unpenalized estimation. In summary, if there is a reasonable suspicion of zero inflation, G(0) should not be smoothed.

Conclusion
Although semiparametric estimation yields a decent fit in INAR models, its performance is often not convincing for small sample sizes. Therefore, we proposed a penalization approach that exploits a qualitative smoothness assumption fulfilled by commonly used innovation distributions. A simulation study showed that our penalization approach provides a large benefit in estimating the innovation distribution, especially for small sample sizes. Additionally, we showed that the combination of unpenalized estimation of INAR coefficients and penalized estimation of the innovation distribution provided the best performance. Future research should investigate whether additional penalization of the INAR coefficients may result in further benefit. Furthermore, as the penalization approach proved to be beneficial for forecasting, one may also think of an application in statistical process control, i.e. for the design of control charts relying on a fitted INAR(1) model. Another interesting issue for future research is the application of our proposed method on integer-valued autoregressive models on Z, such as those proposed by Kim and Park (2008) or Liu et al. (2021).

Appendix
See Figures 12, 13 , 14, 15, 16, 17, 18, 19, 20, 21, 22 and 23. See Tables 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12.  (1) process for the different sample sizes n. We report results for unpenalized (up) and L 1 penalized estimation using either first order (diff1) or second order (diff2) differences     Table 9 Variance, bias and MSE of the first five estimated entries of the PMF for the different sample sizes n in case of a true Poi(1) innovation distribution of an INAR (1) process. We report results for unpenalized (up) and L 1 penalized estimation using either first order (diff1) or second order (diff2) Table 10 Variance, bias and MSE of the first five estimated entries of the PMF for the different sample sizes n in case of a true Poi(1) innovation distribution of an INAR (1) process. We report results for unpenalized (up) and L 2 penalized estimation using either first order (diff1) or second order (diff2) Table 11 Variance, bias and MSE of the first five estimated entries of the PMF for the different sample sizes n in case of a true ZIP( 1 2 ,2) innovation distribution of an INAR (1) process. We report results for unpenalized (up),  Table 12 Variance, bias and MSE of the first five estimated entries of the PMF for the different sample sizes n in case of a true ZIP( 1 2 ,2) innovation distribution of an INAR (1) process. We report results for unpenalized (up),