- Email: [email protected]
First-order observation-driven integer-valued autoregressive processes
ARTICLE IN PRESS
Statistics & Probability Letters 78 (2008) 1–9 www.elsevier.com/locate/stapro
First-order observation-driven integer-valued autoregressive processes Haitao Zheng, Ishwar V. Basawa Department of Statistics, University of Georgia, Athens, GA 30602, USA Received 10 February 2006; received in revised form 22 February 2007; accepted 10 April 2007 Available online 8 May 2007
Abstract A first-order observation-driven integer-valued autoregressive model is introduced. Ergodicity of the process is established. Conditional least squares and maximum likelihood estimators of the model parameters are derived. The performances of these estimators are compared via simulation. The models are applied to a real data set. r 2007 Elsevier B.V. All rights reserved. Keywords: Models of count data; Thinning models; Observation-driven model; INAR models; Random coefficient models; Conditional least squares; Maximum likelihood
1. Introduction Time series models for count data have been studied by several authors in the recent literature. See Davis et al. (1999) for a review of state-space models. See also Fukasawa and Basawa (2002) for a class of observation-driven state-space models. MacDonald and Zucchini (1997) discuss a variety of models for count data including ‘‘thinning’’ models. The integer autoregressive (INAR) models belong to the class of thinning models. INAR models have been discussed extensively by McKenzie (1985a,b), Al-Osh and Alzaid (1987), Alzaid and Al-Osh (1988), among others. Recently, random coefficient INAR models were introduced by Zheng et al. (2006, 2007). In this paper we introduce an observation-driven state-space model for the first-order random coefficient INAR time series. See Zeger and Qaqish (1988) for some early work. A first-order autoregressive model with count (or integer-valued) data is defined through the ‘‘thinning’’ operator which is due to Steutal and Van Harn (1979). Let X be an integer-valued random variable and f 2 ½0; 1Þ, then the thinning operator ‘‘’’ is defined as fX ¼
X X
Bi ,
i¼1
Corresponding author. Tel.: +1 706 5423309; fax: +1 706 5423391.
E-mail address: [email protected] (I.V. Basawa). 0167-7152/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.spl.2007.04.017
(1)
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
2
where fBi g is an i.i.d. Bernoulli random sequence with PðBi ¼ 1Þ ¼ f that is independent of X . With this operator, the INAR(1) model is defined as X t ¼ f X t1 þ Z t ;
tX1,
(2)
where fZt g is a sequence of i.i.d. non-negative integer-valued random variables with mean l and variance s2Z and fZ t g is independent of X 0 . Throughout the rest of the paper, fZ t g will be referred to as innovations. It may be noted that model (2) is related to a branching process with immigration having a Bernoulli offspring distribution. See Guttorp (1991) for a review of inference for branching processes. If X t denotes the number of cancer patients in the tth month, f denotes the survival probability during ðt 1; tÞ, and Z t is the number of new cancer patients added in the tth month, the model in (2) gives a reasonable representation of X t . The survival probability f in (2) may vary with time and it may be random. We now replace f by a random sequence fft g, where ft depends on the past observation X t1 . We then get an observation-driven state-space model where ft represents the ‘‘state’’. Note that in the random coefficient models of Zheng et al. (2006, 2007), fft g is assumed to be a sequence of independent random variables, whereas in this paper the state process fft g is a sequence of dependent (on past observations) random variables. The paper is organized as follows. In Section 2, the observation-driven random coefficient model is described in detail and some statistical properties are established. In Section 3, we propose estimation methods for the model parameters. In Section 4, we present some simulation results for the estimation methods. In Section 5, we apply our model to a real data set. The paper ends with a concluding section. 2. The first-order observation-driven integer-valued autoregressive model A first-order observation-driven integer-valued autoregressive model is defined by the following recursive equation: X t ¼ ft X t1 þ Z t ;
tX1,
(3) 2
where {ft } is a random sequence on ½0; 1Þ and logitðft Þ ¼ a þ bX t1 þ t , t Nð0; s Þ and bp0; {Z t } is an i.i.d. non-negative integer-valued sequence with probability mass function (pmf) f z 40; X 0 , ft g and fZt g are independent. Let l ¼ EðZt Þ. It is easy to see that X t is a Markov chain on f0; 1; 2; . . .g with the following transition probabilities: ! Z 1 minði;jÞ X j ðft Þk1 ð1 ft Þjk1 Pij ¼ PðX t ¼ ijX t1 ¼ jÞ / f z ði kÞ k 0 k¼0 expfðlogitðft Þ a bjÞ2 =ð2s2 Þgdft .
ð4Þ
The Markovity follows from (3). The ergodicity property of the process is of independent interest. It will be very useful to derive the asymptotic properties of the parameter estimates. The moments and conditional moments will be useful in obtaining the appropriate estimating equations for parameter estimation. Proposition 2.1. Let at ¼ Eðft jX t1 Þ. For tX1, (i) EðX t jX t1 Þ ¼ at X t1 þ l; (ii) VarðX t jX t1 ; ft Þ ¼ ft ð1 ft ÞX t1 þ s2z ; (iii) VarðX t jX t1 Þ ¼ s2ft jX t1 X 2t1 þ ðat ð1 at Þ s2ft jX t1 ÞX t1 þ s2z , where Z 1 1 ðft at Þ2 2 expfðlogitðft Þ a bX t1 Þ2 =ð2s2 Þg dft ; sft jX t1 ¼ pffiffiffiffiffiffi 2ps 0 ft ð1 ft Þ (iv) CovðX t ; X tþ1 Þ ¼ Eðatþ1 X 2t Þ Eðatþ1 X t ÞEðX t Þ. The proof is omitted since it is elementary.
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
3
Proposition 2.2. The process fX t g is an ergodic Markov chain. Proof. According to Theorem 3.1 of Tweedie (1975) (see also, Meyn and Tweedie, 1993), the sufficient condition for fX t g to be ergodic is that there exists a set K and a non-negative measurable function g on state space w such that Z ¯ Pðx; dyÞgðyÞpgðxÞ 1; x 2 K, (5) w
and, for some fixed B Z Pðx; dyÞgðyÞ ¼ lðxÞpBo1;
x 2 K,
(6)
w
where Pðx; AÞ ¼ PðX 1 2 AjX 0 ¼ xÞ. We assume that logitðft Þ ¼ a þ bX t1 þ t and bp0. We have 1 1 1 Eðft jX t1 Þ ¼ E X t1 pE X t1 ¼ E . 1 þ eabX t1 t 1 þ eat 1 þ eat Choose a constant 0oco1, such that Eð1=ð1 þ eat ÞÞoc, we then have Eðft jX t1 Þoc for all t. Let gðxÞ ¼ x. We have Z gðyÞ dPðX 1 ¼ yjX 0 ¼ x0 Þ ¼ EðX 1 jX 0 ¼ x0 Þ ¼ x0 Eðf1 jX 0 ¼ x0 Þ þ locx0 þ l. Therefore, if we let N ¼ ½ðl þ 1Þ=ð1 cÞ þ 1, where [] denotes the integer part of a real number, then for x0 XN, we have cx0 þ lpx0 1 ¼ gðx0 Þ 1, and for 0px0 pN 1, Z gðyÞ dPðX 1 ¼ yjX 0 ¼ x0 Þ ¼ EðX 1 jX 0 ¼ x0 Þpx0 þ loN þ lo1. Let K¯ ¼ fx : N; N þ 1; . . . ; g and K ¼ fx : 0; 1; . . . ; N 1g. We have thus verified the sufficient conditions (5) and (6) for {X t } to be ergodic. & Models with covariates: In practice, we may have covariates related to the model. The model will then be non-stationary. In this situation, we can introduce the covariates into ft and/or the innovation in model (3). For instance, logitðft Þ ¼ a þ bX t1 þ b0 wt þ t , t Nð0; s2 Þ, bp0, wt is a non-random covariate vector and b is the corresponding coefficient vector; wt can contain the trend and seasonal terms. The parameters from Zt may depend on the covariates as well, for instance, if Z t is a Poisson sequence, then we can have logðlt Þ ¼ m þ s0 wt , where lt ¼ EðZt Þ. See Zeger and Qaqish (1988) and Davis et al. (1999) for related models with covariates. Models with covariates are used in Section 3. 3. Estimation methods In this paper, we consider two methods, namely, the conditional least squares (CLS) and the maximum likelihood (ML). An advantage of CLS is that it does not require specifying the exact family of distributions for fft g and the innovations. In this section, we focus on CLS and ML methods for the model in (3) without the covariates. 3.1. Conditional least squares P Let QðhÞ ¼ ni¼1 ðX i EðX i jX i1 ÞÞ2 and h be the parameter vector of the model contained in EðX t jX t1 Þ. Then the CLS estimator is given by h^ CLS ¼ argminh ðQðhÞÞ
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
4
i.e. h^ CLS is a solution of qQðhÞ=qh^ ¼ 0 or n X qEðX i jX i1 Þ
qh
i¼1
ðX i EðX i jX i1 ÞÞ ¼ 0.
(7)
For the model in (3), hT ¼ ða; b; s2 ; lÞ. Theorem 3.1. Under the regularity conditions (see Klimko and Nelson, 1978), the CLS estimator h^ CLS is consistent and has the following asymptotic distribution: pffiffiffi ^ L nðhCLS hÞ ! MVNð0p1 ; ðV T W 1 V Þ1 Þ, (8) where p is the dimension of h, W ¼ ðEðu21 ðhÞqEðX 1 jX 0 Þ=qyi qEðX 1 jX 0 Þ=qyj ÞÞ1pi;jpp , u1 ðhÞ ¼ X 1 EðX 1 jX 0 Þ and V ¼ ðEðqEðX 1 jX 0 Þ=qyi qEðX 1 jX 0 Þ=qyj ÞÞ1pi;jpp ðEðu1 ðhÞq2 EðX 1 jX 0 Þ=qyi qyj ÞÞ1pi;jpp . Proof. According to Klimko and Nelson (1978), if g ¼ EðX t jX t1 Þ, Theorem 3.1 holds if the following conditions are satisfied. (i) qg=qyi , q2 g=qyi qyj , q3 g=qyi qyj qyk , 1pi; j; kpp exist and continuous; (ii) for 1pi; jpp, EjðX t gÞqg=qyi jo1, EjðX t gÞq2 g=qyi qyj jo1, Ejqg=qyi qg=qyj jo1; (iii) for 1pi; j; kpp, there exist functions H ð0Þ ðX t1 ; . . . ; X 0 Þ;
H ð1Þ i ðX t1 ; . . . ; X 0 Þ;
H ð2Þ ij ðX t1 ; . . . ; X 0 Þ;
H ð3Þ ijk ðX t1 ; . . . ; X 0 Þ
such that jgjpH ð0Þ ;
jqg=qyi jpH ð1Þ i ;
jq2 g=qyi qyj jpH ð2Þ ij ;
jq3 g=qyi qyj qyk jpH ð3Þ ijk
and EjX t H ð3Þ ijk ðX t1 ; . . . ; X 0 Þjo1, H ð0Þ ðX t1 ; . . . ; X 0 Þ H ð3Þ ijk ðX t1 ; . . . ; X 0 Þo1, ð2Þ EjH ð1Þ i ðX t1 ; . . . ; X 0 Þ H ij ðX t1 ; . . . ; X 0 Þjo1.
(iv)
EðX t jX t1 ; . . . ; X 0 Þ ¼ EðX t jX t1 Þ
a:e:;
tX1,
Eðu2t ðhÞjqgðhÞ=qyi qgðhÞ=qyj jÞo1, where ut ðhÞ ¼ X t EðX t jX t1 Þ. R1 For our model in (3), we have EðX 1 jX 0 Þ ¼ 1 ð1=ð1 þ expða bX 0 1 ÞÞf ð1 Þ d1 þ l, where f ð1 Þ is the standard normal density function. Obviously, qEðX 1 jX 0 Þ=qyi , q2 EðX 1 jX 0 Þ=qyi qyj , q3 EðX 1 jX 0 Þ=qyi qyj qyk , 1pi; j; kpp exist and are continuous. We also have jqEðX 1 jX 0 Þ=qajp1; jqEðX 1 jX 0 Þ=qsjp2=s;
jqEðX 1 jX 0 Þ=qbjpX 0 , qEðX 1 jX 0 Þ=ql ¼ 1,
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
5
jq2 EðX 1 jX 0 Þ=qaqbjp3X 0 , jq2 EðX 1 jX 0 Þ=qaqsjp2=s, jq2 EðX 1 jX 0 Þ=qbqsjp2X 0 =s, jq3 EðX 1 jX 0 Þ=qaqbqsjp6X 0 =s. The second and third derivatives that involve l are 0. If EðX 40 Þo1, then all the conditions are satisfied and hence Theorem 3.1 holds. & 3.2. Maximum likelihood estimators (MLEs) It is easy to derive the likelihood function using the Markov property. From (3), we can get the likelihood function LðgÞ ¼ PðX 0 ¼ x0 Þ
n Y
PðX t ¼ xt jX t1 ¼ xt1 Þ,
t¼1
where g is the parameter vector and the parameters are from f z and Pf1 . The transition probabilities PðX t ¼ xt jX t1 ¼ xt1 Þ are given by (4). The MLE g^ ML is the value of g which maximizes the above likelihood function. Theorem 3.2. Under the regularity conditions (see Billingsley, 1961), the ML estimator g^ ML is consistent and has the following asymptotic distribution: pffiffiffi L nð^gML gÞ ! MVNð0q1 ; I 1 Þ,
(9)
where q is the dimension of g, I ¼ Eðq logðPðX 1 jX 0 ÞÞ=qZi q logðPðX 1 jX 0 ÞÞ=qZj Þ1pi;jpq , the Fisher information matrix. For our model in (3), gT ¼ ða; b; s2 ; l; s2z Þ. The proof is omitted here since it is straightforward to verify the regularity conditions. It may be noted that the ergodicity property proved in Proposition 2.2 plays a key role in the proofs of both Theorems 3.1 and 3.2. 4. Simulation studies Consider the model X t ¼ ft X t1 þ Z t , where logðft =ð1 ft ÞÞ ¼ a þ bX t1 þ t , and ft g is an i.i.d. sequence of normal random variables with mean 0 and variance s2 ; {Zt } is an i.i.d. Poisson sequence with mean l. In the simulation, we fixed X 0 at 1. Fig. 1 is a typical sample path from this model, for a sample size n ¼ 200, when a ¼ 1:4; b ¼ 0:67; s ¼ 3 and l ¼ 1. We use the above model to generate data, and then use CLS and MLE methods to estimate the parameters. The ML score equation under this assumed parametric model for the Z is given by R1 k PminðX t ;X t1 Þ X t1 X t k X t1 k =ðX t kÞ!Þ 0 ft ð1 ft Þ dPft ðl n ðq=qsÞ k¼0 X k ¼ 0, R1 k PminðX t ;X t1 Þ X t1 X t k X t1 k t¼1 ðl =ðX kÞ!Þ f ð1 f Þ dP t f t k¼0 t 0 t k
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
6
where s ¼ ða; b; l; sÞ and dPft ¼
1 ðlogðft =ð1 ft ÞÞ a bX t1 Þ2 pffiffiffiffiffiffi exp dft . 2s2 ft ð1 ft Þ 2ps
The (univariate) integrals in the above expressions were computed using Gaussian Quadrature routine in R. We chose various combinations of the parameters a, b and l for s ¼ 1 and 3. 6 5
y
4 3 2 1 0 0
50
100 x
150
200
Fig. 1. A typical sample path from a ODRCINAR(1) model. Table 1 Estimates for parameter a, b and l for s ¼ 1 a ¼ 1:4
b ¼ 0:67
l¼1
s¼1
CLS n ¼ 50 n ¼ 100 n ¼ 200
(1.35, 1.05) (1.38, 0.75) (1.43, 0.46)
(0.74, 0.39) (0.72, 0.24) (0.70, 0.13)
(1.03, 0.23) (1.01, 0.19) (1.00, 0.14)
(0.99, 1.04) (0.99, 0.82) (0.95, 0.53)
MLE n ¼ 50 n ¼ 100 n ¼ 200
(1.58, 0.95) (1.53, 0.77) (1.42, 0.39)
(0.76, 0.40) (0.73, 0.27) (0.69, 0.15)
(1.02, 0.22) (1.01, 0.17) (1.01, 0.12)
(0.90, 0.71) (1.00, 0.54) (0.97, 0.40)
a ¼ 1:4
b ¼ 0:67
l¼1
s¼3
CLS n ¼ 50 n ¼ 100 n ¼ 200
(1.69, 1.17) (1.58, 0.85) (1.51, 0.57)
(0.89, 0.47) (0.77, 0.31) (0.73, 0.20)
(1.00, 0.26) (0.99, 0.17) (1.00, 0.12)
(2.52, 1.44) (2.72, 1.07) (2.84, 0.71)
MLE n ¼ 50 n ¼ 100 n ¼ 200
(1.56, 1.14) (1.43, 0.78) (1.41, 0.49)
(0.78, 0.44) (0.69, 0.30) (0.69, 0.19)
(1.04, 0.22) (1.02, 0.15) (1.01, 0.10)
(2.39, 0.90) (2.65, 0.61) (2.80, 0.41)
Table 2 Estimates for parameter a, b and l for s ¼ 3
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
7
We computed the point estimates and the standard errors (SE) based on 500 replications for each parameter combination. These values are reported within parenthesis in Tables 1 and 2. In the following tables, the format (estimate, SE) is used; for example, ð0:70; 0:13Þ means that the estimate is 0:70 and standard error is 0.13. From the simulation results, we can see that these two methods are very close and the MLE is slightly better than the CLS. When s is small, CLS tends to underestimate a and b while MLE tends to overestimate a and underestimate b; when s is large, CLS and MLE overestimate a and underestimate b and s. As n increases, however, both the estimates seem reasonable with MLE being more efficient. 5. Polio counts data analysis The data set consists of the monthly number of cases of poliomyelitis in the US for the year 1970–1983 as reported by the Centers of Disease Control. These data were originally presented in Zeger and Qaqish (1988) and have become a standard data set in the field of longitudinal count data. Here we present a parametric analysis for this data, using the model we have introduced. See the sample path plot for the data in Fig. 2. The plot reveals some seasonality and the possibility of a slight decreasing trend. From the ACF and PACF plots in Figs. 3 and 4, the data may come from an AR(1) process. But the data are not stationary due to the presence of trend and seasonality. We can easily allow for trend and seasonality in our model as indicated below. The basic model is X t ¼ ft X t1 þ Z t , where fZ t g is a Poisson sequence with EðZ t Þ ¼ lt . With different model assumptions on ft and lt , we consider four such models below. Model I: ft ¼ f, Z t Poissonðlt Þ and logðlt Þ ¼ b0 þ b1 t=1000 þ b2 cosð2pt=12Þ þ b3 sinð2pt=12Þ þ b4 cosð2pt=6Þ þ b5 sinð2pt=6Þ. Model II: logitðft Þ ¼ a þ U t , Z t Poissonðlt Þ and logðlt Þ ¼ b0 þ b1 t=1000 þ b2 cosð2pt=12Þ þ b3 sinð2pt=12Þ þ b4 cosð2pt=6Þ þ b5 sinð2pt=6Þ.
14 12 10
y
8 6 4 2 0 0
50
100 x
Fig. 2. Polio counts plot.
105
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
8
Series y 1.0
0.8
ACF
0.6
0.4
0.2
0.0
0
5
10
15
20
Lag Fig. 3. Polio counts ACF plot.
Series y 0.3
Partial ACF
0.2
0.1
0.0
-0.1
5
10
15
20
Lag Fig. 4. Polio counts PACF plot.
Model III: Z t PoissonðlÞ and logitðft Þ ¼ a þ bX t1 þ
b1 t þ b2 cosð2pt=12Þ þ b3 sinð2pt=12Þ þ b4 cosð2pt=6Þ þ b5 sinð2pt=6Þ þ U t . 1000
Model IV: logitðft Þ ¼ a þ bX t1 þ U t , Z t Poissonðlt Þ, and log lt ¼ b0 þ b1 t=1000 þ b2 cosð2pt=12Þ þ b3 sinð2pt=12Þ þ b4 cosð2pt=6Þ þ b5 sinð2pt=6Þ.
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
9
Table 3
a b s l b0 b1 b2 b3 b4 b5 AIC
Model I
Model II
Model III
Model IV
2.01
5.02
3.65 0.35 4.50 1.03
2.70 0.33 5.00
5.16 0.43 4.61 0.37 0.40 0.28 0.35 551.82
0.29 3.61 0.35 0.39 0.22 0.28 530.12
2.42 2.56 3.29 2.26 3.09 550.77
0.26 3.82 0.36 0.38 0.23 0.28 534.59
Model I is a fixed coefficient INAR(1) process whereas Models II–IV are random coefficient models. Now we use the ML method to estimate the parameters in the above models. The results are summarized in Table 3. The AICs are given in the table for each model as well. Using the AIC criterion, the random coefficient Models II–IV are better than the fixed coefficient Model I. Model II is slightly better than Model IV for this data set. Model IV, however, is a reasonable alternative to model II with the additional advantage of incorporating the dependence of ft on the past observation X t1 . Models I and III have similar AIC. For this data set, we may choose Models II and IV in preference to Models I and III. 6. Summary and conclusions In this paper, we have introduced a new model by combining a thinning model with an observation-driven model. The autoregressive parameter is allowed to vary randomly over time and depends on the past observation. The ergodicity of the process is established. Conditional least squares (CLS) and maximum likelihood (ML) are used to estimate the parameters. In the simulation study, we compare CLS and ML estimators. The simulation results show that the two methods give similar estimates and the ML is slightly better. The model with covariates is applied to a real data set. It is shown that the random coefficient models are a better fit than the fixed coefficient model for this data set. References Al-Osh, M.A., Alzaid, A.A., 1987. First order integer-valued autoregressive (INAR(1)) processes. J. Time Series Anal. 8, 261–275. Alzaid, A.A., Al-Osh, M.A., 1988. First order integer-valued autoregressive (INAR(1)) processes: distributional and regression properties. Statist. Neerlandica 42, 53–61. Billingsley, P., 1961. Statistical Inference for Markov Processes. The University of Chicago Press, Chicago. Davis, R.A., Dunsmuir, W.T.M., Wang, Y., 1999. Modeling time series for count data. In: Ghosh, S. (Ed.), Asymptotics, Nonparametrics and Time Series. Marcel Dekker, New York, pp. 63–114. Fukasawa, T., Basawa, I.V., 2002. Estimation for a class of generalized state-space time series models. Statist. Probab. Lett. 60, 459–473. Guttorp, P., 1991. Statistical Inference for Branching Processes. Wiley, New York. Klimko, L.A., Nelson, P.I., 1978. On conditional least squares estimation for stochastic processes. Ann. Statist. 6, 629–642. MacDonald, I.L., Zucchini, W.Z., 1997. Hidden Markov and Other Models for Discrete-valued Time Series. Chapman & Hall, London. McKenzie, E., 1985a. Contribution to the discussion of Lawrence and Lewis. J. Roy. Statist. Soc. B 47, 187–188. McKenzie, E., 1985b. Some simple models for discrete variate time series. Water Resour. Bull. 21, 645–650. Meyn, S.P., Tweedie, R.L., 1993. Markov Chains and Stochastic Stability. Springer, Berlin. Steutal, F., Van Harn, K., 1979. Discrete analogues of self-decomposability and stability. Ann. Probab. 7, 893–899. Tweedie, R.L., 1975. Sufficient conditions for regularity, recurrence and ergocidicity of Markov processes. Stochastic Process Appl. 3, 385–403. Zeger, S.L., Qaqish, B., 1988. Markov regression models for time series: a quasi-likelihood approach. Biometrics 44, 1019–1031. Zheng, H.T., Basawa, I.V., Datta, S., 2006. Inference for pth order random coefficient integer-valued autoregressive processes. J. Time Series Anal. 27, 411–440. Zheng, H.T., Basawa, I.V., Datta, S., 2007. First order random coefficient integer-valued autoregressive processes. J. Statist. Plann. Inference 137 (1), 212–229.
Statistics & Probability Letters 78 (2008) 1–9 www.elsevier.com/locate/stapro
First-order observation-driven integer-valued autoregressive processes Haitao Zheng, Ishwar V. Basawa Department of Statistics, University of Georgia, Athens, GA 30602, USA Received 10 February 2006; received in revised form 22 February 2007; accepted 10 April 2007 Available online 8 May 2007
Abstract A first-order observation-driven integer-valued autoregressive model is introduced. Ergodicity of the process is established. Conditional least squares and maximum likelihood estimators of the model parameters are derived. The performances of these estimators are compared via simulation. The models are applied to a real data set. r 2007 Elsevier B.V. All rights reserved. Keywords: Models of count data; Thinning models; Observation-driven model; INAR models; Random coefficient models; Conditional least squares; Maximum likelihood
1. Introduction Time series models for count data have been studied by several authors in the recent literature. See Davis et al. (1999) for a review of state-space models. See also Fukasawa and Basawa (2002) for a class of observation-driven state-space models. MacDonald and Zucchini (1997) discuss a variety of models for count data including ‘‘thinning’’ models. The integer autoregressive (INAR) models belong to the class of thinning models. INAR models have been discussed extensively by McKenzie (1985a,b), Al-Osh and Alzaid (1987), Alzaid and Al-Osh (1988), among others. Recently, random coefficient INAR models were introduced by Zheng et al. (2006, 2007). In this paper we introduce an observation-driven state-space model for the first-order random coefficient INAR time series. See Zeger and Qaqish (1988) for some early work. A first-order autoregressive model with count (or integer-valued) data is defined through the ‘‘thinning’’ operator which is due to Steutal and Van Harn (1979). Let X be an integer-valued random variable and f 2 ½0; 1Þ, then the thinning operator ‘‘’’ is defined as fX ¼
X X
Bi ,
i¼1
Corresponding author. Tel.: +1 706 5423309; fax: +1 706 5423391.
E-mail address: [email protected] (I.V. Basawa). 0167-7152/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.spl.2007.04.017
(1)
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
2
where fBi g is an i.i.d. Bernoulli random sequence with PðBi ¼ 1Þ ¼ f that is independent of X . With this operator, the INAR(1) model is defined as X t ¼ f X t1 þ Z t ;
tX1,
(2)
where fZt g is a sequence of i.i.d. non-negative integer-valued random variables with mean l and variance s2Z and fZ t g is independent of X 0 . Throughout the rest of the paper, fZ t g will be referred to as innovations. It may be noted that model (2) is related to a branching process with immigration having a Bernoulli offspring distribution. See Guttorp (1991) for a review of inference for branching processes. If X t denotes the number of cancer patients in the tth month, f denotes the survival probability during ðt 1; tÞ, and Z t is the number of new cancer patients added in the tth month, the model in (2) gives a reasonable representation of X t . The survival probability f in (2) may vary with time and it may be random. We now replace f by a random sequence fft g, where ft depends on the past observation X t1 . We then get an observation-driven state-space model where ft represents the ‘‘state’’. Note that in the random coefficient models of Zheng et al. (2006, 2007), fft g is assumed to be a sequence of independent random variables, whereas in this paper the state process fft g is a sequence of dependent (on past observations) random variables. The paper is organized as follows. In Section 2, the observation-driven random coefficient model is described in detail and some statistical properties are established. In Section 3, we propose estimation methods for the model parameters. In Section 4, we present some simulation results for the estimation methods. In Section 5, we apply our model to a real data set. The paper ends with a concluding section. 2. The first-order observation-driven integer-valued autoregressive model A first-order observation-driven integer-valued autoregressive model is defined by the following recursive equation: X t ¼ ft X t1 þ Z t ;
tX1,
(3) 2
where {ft } is a random sequence on ½0; 1Þ and logitðft Þ ¼ a þ bX t1 þ t , t Nð0; s Þ and bp0; {Z t } is an i.i.d. non-negative integer-valued sequence with probability mass function (pmf) f z 40; X 0 , ft g and fZt g are independent. Let l ¼ EðZt Þ. It is easy to see that X t is a Markov chain on f0; 1; 2; . . .g with the following transition probabilities: ! Z 1 minði;jÞ X j ðft Þk1 ð1 ft Þjk1 Pij ¼ PðX t ¼ ijX t1 ¼ jÞ / f z ði kÞ k 0 k¼0 expfðlogitðft Þ a bjÞ2 =ð2s2 Þgdft .
ð4Þ
The Markovity follows from (3). The ergodicity property of the process is of independent interest. It will be very useful to derive the asymptotic properties of the parameter estimates. The moments and conditional moments will be useful in obtaining the appropriate estimating equations for parameter estimation. Proposition 2.1. Let at ¼ Eðft jX t1 Þ. For tX1, (i) EðX t jX t1 Þ ¼ at X t1 þ l; (ii) VarðX t jX t1 ; ft Þ ¼ ft ð1 ft ÞX t1 þ s2z ; (iii) VarðX t jX t1 Þ ¼ s2ft jX t1 X 2t1 þ ðat ð1 at Þ s2ft jX t1 ÞX t1 þ s2z , where Z 1 1 ðft at Þ2 2 expfðlogitðft Þ a bX t1 Þ2 =ð2s2 Þg dft ; sft jX t1 ¼ pffiffiffiffiffiffi 2ps 0 ft ð1 ft Þ (iv) CovðX t ; X tþ1 Þ ¼ Eðatþ1 X 2t Þ Eðatþ1 X t ÞEðX t Þ. The proof is omitted since it is elementary.
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
3
Proposition 2.2. The process fX t g is an ergodic Markov chain. Proof. According to Theorem 3.1 of Tweedie (1975) (see also, Meyn and Tweedie, 1993), the sufficient condition for fX t g to be ergodic is that there exists a set K and a non-negative measurable function g on state space w such that Z ¯ Pðx; dyÞgðyÞpgðxÞ 1; x 2 K, (5) w
and, for some fixed B Z Pðx; dyÞgðyÞ ¼ lðxÞpBo1;
x 2 K,
(6)
w
where Pðx; AÞ ¼ PðX 1 2 AjX 0 ¼ xÞ. We assume that logitðft Þ ¼ a þ bX t1 þ t and bp0. We have 1 1 1 Eðft jX t1 Þ ¼ E X t1 pE X t1 ¼ E . 1 þ eabX t1 t 1 þ eat 1 þ eat Choose a constant 0oco1, such that Eð1=ð1 þ eat ÞÞoc, we then have Eðft jX t1 Þoc for all t. Let gðxÞ ¼ x. We have Z gðyÞ dPðX 1 ¼ yjX 0 ¼ x0 Þ ¼ EðX 1 jX 0 ¼ x0 Þ ¼ x0 Eðf1 jX 0 ¼ x0 Þ þ locx0 þ l. Therefore, if we let N ¼ ½ðl þ 1Þ=ð1 cÞ þ 1, where [] denotes the integer part of a real number, then for x0 XN, we have cx0 þ lpx0 1 ¼ gðx0 Þ 1, and for 0px0 pN 1, Z gðyÞ dPðX 1 ¼ yjX 0 ¼ x0 Þ ¼ EðX 1 jX 0 ¼ x0 Þpx0 þ loN þ lo1. Let K¯ ¼ fx : N; N þ 1; . . . ; g and K ¼ fx : 0; 1; . . . ; N 1g. We have thus verified the sufficient conditions (5) and (6) for {X t } to be ergodic. & Models with covariates: In practice, we may have covariates related to the model. The model will then be non-stationary. In this situation, we can introduce the covariates into ft and/or the innovation in model (3). For instance, logitðft Þ ¼ a þ bX t1 þ b0 wt þ t , t Nð0; s2 Þ, bp0, wt is a non-random covariate vector and b is the corresponding coefficient vector; wt can contain the trend and seasonal terms. The parameters from Zt may depend on the covariates as well, for instance, if Z t is a Poisson sequence, then we can have logðlt Þ ¼ m þ s0 wt , where lt ¼ EðZt Þ. See Zeger and Qaqish (1988) and Davis et al. (1999) for related models with covariates. Models with covariates are used in Section 3. 3. Estimation methods In this paper, we consider two methods, namely, the conditional least squares (CLS) and the maximum likelihood (ML). An advantage of CLS is that it does not require specifying the exact family of distributions for fft g and the innovations. In this section, we focus on CLS and ML methods for the model in (3) without the covariates. 3.1. Conditional least squares P Let QðhÞ ¼ ni¼1 ðX i EðX i jX i1 ÞÞ2 and h be the parameter vector of the model contained in EðX t jX t1 Þ. Then the CLS estimator is given by h^ CLS ¼ argminh ðQðhÞÞ
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
4
i.e. h^ CLS is a solution of qQðhÞ=qh^ ¼ 0 or n X qEðX i jX i1 Þ
qh
i¼1
ðX i EðX i jX i1 ÞÞ ¼ 0.
(7)
For the model in (3), hT ¼ ða; b; s2 ; lÞ. Theorem 3.1. Under the regularity conditions (see Klimko and Nelson, 1978), the CLS estimator h^ CLS is consistent and has the following asymptotic distribution: pffiffiffi ^ L nðhCLS hÞ ! MVNð0p1 ; ðV T W 1 V Þ1 Þ, (8) where p is the dimension of h, W ¼ ðEðu21 ðhÞqEðX 1 jX 0 Þ=qyi qEðX 1 jX 0 Þ=qyj ÞÞ1pi;jpp , u1 ðhÞ ¼ X 1 EðX 1 jX 0 Þ and V ¼ ðEðqEðX 1 jX 0 Þ=qyi qEðX 1 jX 0 Þ=qyj ÞÞ1pi;jpp ðEðu1 ðhÞq2 EðX 1 jX 0 Þ=qyi qyj ÞÞ1pi;jpp . Proof. According to Klimko and Nelson (1978), if g ¼ EðX t jX t1 Þ, Theorem 3.1 holds if the following conditions are satisfied. (i) qg=qyi , q2 g=qyi qyj , q3 g=qyi qyj qyk , 1pi; j; kpp exist and continuous; (ii) for 1pi; jpp, EjðX t gÞqg=qyi jo1, EjðX t gÞq2 g=qyi qyj jo1, Ejqg=qyi qg=qyj jo1; (iii) for 1pi; j; kpp, there exist functions H ð0Þ ðX t1 ; . . . ; X 0 Þ;
H ð1Þ i ðX t1 ; . . . ; X 0 Þ;
H ð2Þ ij ðX t1 ; . . . ; X 0 Þ;
H ð3Þ ijk ðX t1 ; . . . ; X 0 Þ
such that jgjpH ð0Þ ;
jqg=qyi jpH ð1Þ i ;
jq2 g=qyi qyj jpH ð2Þ ij ;
jq3 g=qyi qyj qyk jpH ð3Þ ijk
and EjX t H ð3Þ ijk ðX t1 ; . . . ; X 0 Þjo1, H ð0Þ ðX t1 ; . . . ; X 0 Þ H ð3Þ ijk ðX t1 ; . . . ; X 0 Þo1, ð2Þ EjH ð1Þ i ðX t1 ; . . . ; X 0 Þ H ij ðX t1 ; . . . ; X 0 Þjo1.
(iv)
EðX t jX t1 ; . . . ; X 0 Þ ¼ EðX t jX t1 Þ
a:e:;
tX1,
Eðu2t ðhÞjqgðhÞ=qyi qgðhÞ=qyj jÞo1, where ut ðhÞ ¼ X t EðX t jX t1 Þ. R1 For our model in (3), we have EðX 1 jX 0 Þ ¼ 1 ð1=ð1 þ expða bX 0 1 ÞÞf ð1 Þ d1 þ l, where f ð1 Þ is the standard normal density function. Obviously, qEðX 1 jX 0 Þ=qyi , q2 EðX 1 jX 0 Þ=qyi qyj , q3 EðX 1 jX 0 Þ=qyi qyj qyk , 1pi; j; kpp exist and are continuous. We also have jqEðX 1 jX 0 Þ=qajp1; jqEðX 1 jX 0 Þ=qsjp2=s;
jqEðX 1 jX 0 Þ=qbjpX 0 , qEðX 1 jX 0 Þ=ql ¼ 1,
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
5
jq2 EðX 1 jX 0 Þ=qaqbjp3X 0 , jq2 EðX 1 jX 0 Þ=qaqsjp2=s, jq2 EðX 1 jX 0 Þ=qbqsjp2X 0 =s, jq3 EðX 1 jX 0 Þ=qaqbqsjp6X 0 =s. The second and third derivatives that involve l are 0. If EðX 40 Þo1, then all the conditions are satisfied and hence Theorem 3.1 holds. & 3.2. Maximum likelihood estimators (MLEs) It is easy to derive the likelihood function using the Markov property. From (3), we can get the likelihood function LðgÞ ¼ PðX 0 ¼ x0 Þ
n Y
PðX t ¼ xt jX t1 ¼ xt1 Þ,
t¼1
where g is the parameter vector and the parameters are from f z and Pf1 . The transition probabilities PðX t ¼ xt jX t1 ¼ xt1 Þ are given by (4). The MLE g^ ML is the value of g which maximizes the above likelihood function. Theorem 3.2. Under the regularity conditions (see Billingsley, 1961), the ML estimator g^ ML is consistent and has the following asymptotic distribution: pffiffiffi L nð^gML gÞ ! MVNð0q1 ; I 1 Þ,
(9)
where q is the dimension of g, I ¼ Eðq logðPðX 1 jX 0 ÞÞ=qZi q logðPðX 1 jX 0 ÞÞ=qZj Þ1pi;jpq , the Fisher information matrix. For our model in (3), gT ¼ ða; b; s2 ; l; s2z Þ. The proof is omitted here since it is straightforward to verify the regularity conditions. It may be noted that the ergodicity property proved in Proposition 2.2 plays a key role in the proofs of both Theorems 3.1 and 3.2. 4. Simulation studies Consider the model X t ¼ ft X t1 þ Z t , where logðft =ð1 ft ÞÞ ¼ a þ bX t1 þ t , and ft g is an i.i.d. sequence of normal random variables with mean 0 and variance s2 ; {Zt } is an i.i.d. Poisson sequence with mean l. In the simulation, we fixed X 0 at 1. Fig. 1 is a typical sample path from this model, for a sample size n ¼ 200, when a ¼ 1:4; b ¼ 0:67; s ¼ 3 and l ¼ 1. We use the above model to generate data, and then use CLS and MLE methods to estimate the parameters. The ML score equation under this assumed parametric model for the Z is given by R1 k PminðX t ;X t1 Þ X t1 X t k X t1 k =ðX t kÞ!Þ 0 ft ð1 ft Þ dPft ðl n ðq=qsÞ k¼0 X k ¼ 0, R1 k PminðX t ;X t1 Þ X t1 X t k X t1 k t¼1 ðl =ðX kÞ!Þ f ð1 f Þ dP t f t k¼0 t 0 t k
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
6
where s ¼ ða; b; l; sÞ and dPft ¼
1 ðlogðft =ð1 ft ÞÞ a bX t1 Þ2 pffiffiffiffiffiffi exp dft . 2s2 ft ð1 ft Þ 2ps
The (univariate) integrals in the above expressions were computed using Gaussian Quadrature routine in R. We chose various combinations of the parameters a, b and l for s ¼ 1 and 3. 6 5
y
4 3 2 1 0 0
50
100 x
150
200
Fig. 1. A typical sample path from a ODRCINAR(1) model. Table 1 Estimates for parameter a, b and l for s ¼ 1 a ¼ 1:4
b ¼ 0:67
l¼1
s¼1
CLS n ¼ 50 n ¼ 100 n ¼ 200
(1.35, 1.05) (1.38, 0.75) (1.43, 0.46)
(0.74, 0.39) (0.72, 0.24) (0.70, 0.13)
(1.03, 0.23) (1.01, 0.19) (1.00, 0.14)
(0.99, 1.04) (0.99, 0.82) (0.95, 0.53)
MLE n ¼ 50 n ¼ 100 n ¼ 200
(1.58, 0.95) (1.53, 0.77) (1.42, 0.39)
(0.76, 0.40) (0.73, 0.27) (0.69, 0.15)
(1.02, 0.22) (1.01, 0.17) (1.01, 0.12)
(0.90, 0.71) (1.00, 0.54) (0.97, 0.40)
a ¼ 1:4
b ¼ 0:67
l¼1
s¼3
CLS n ¼ 50 n ¼ 100 n ¼ 200
(1.69, 1.17) (1.58, 0.85) (1.51, 0.57)
(0.89, 0.47) (0.77, 0.31) (0.73, 0.20)
(1.00, 0.26) (0.99, 0.17) (1.00, 0.12)
(2.52, 1.44) (2.72, 1.07) (2.84, 0.71)
MLE n ¼ 50 n ¼ 100 n ¼ 200
(1.56, 1.14) (1.43, 0.78) (1.41, 0.49)
(0.78, 0.44) (0.69, 0.30) (0.69, 0.19)
(1.04, 0.22) (1.02, 0.15) (1.01, 0.10)
(2.39, 0.90) (2.65, 0.61) (2.80, 0.41)
Table 2 Estimates for parameter a, b and l for s ¼ 3
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
7
We computed the point estimates and the standard errors (SE) based on 500 replications for each parameter combination. These values are reported within parenthesis in Tables 1 and 2. In the following tables, the format (estimate, SE) is used; for example, ð0:70; 0:13Þ means that the estimate is 0:70 and standard error is 0.13. From the simulation results, we can see that these two methods are very close and the MLE is slightly better than the CLS. When s is small, CLS tends to underestimate a and b while MLE tends to overestimate a and underestimate b; when s is large, CLS and MLE overestimate a and underestimate b and s. As n increases, however, both the estimates seem reasonable with MLE being more efficient. 5. Polio counts data analysis The data set consists of the monthly number of cases of poliomyelitis in the US for the year 1970–1983 as reported by the Centers of Disease Control. These data were originally presented in Zeger and Qaqish (1988) and have become a standard data set in the field of longitudinal count data. Here we present a parametric analysis for this data, using the model we have introduced. See the sample path plot for the data in Fig. 2. The plot reveals some seasonality and the possibility of a slight decreasing trend. From the ACF and PACF plots in Figs. 3 and 4, the data may come from an AR(1) process. But the data are not stationary due to the presence of trend and seasonality. We can easily allow for trend and seasonality in our model as indicated below. The basic model is X t ¼ ft X t1 þ Z t , where fZ t g is a Poisson sequence with EðZ t Þ ¼ lt . With different model assumptions on ft and lt , we consider four such models below. Model I: ft ¼ f, Z t Poissonðlt Þ and logðlt Þ ¼ b0 þ b1 t=1000 þ b2 cosð2pt=12Þ þ b3 sinð2pt=12Þ þ b4 cosð2pt=6Þ þ b5 sinð2pt=6Þ. Model II: logitðft Þ ¼ a þ U t , Z t Poissonðlt Þ and logðlt Þ ¼ b0 þ b1 t=1000 þ b2 cosð2pt=12Þ þ b3 sinð2pt=12Þ þ b4 cosð2pt=6Þ þ b5 sinð2pt=6Þ.
14 12 10
y
8 6 4 2 0 0
50
100 x
Fig. 2. Polio counts plot.
105
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
8
Series y 1.0
0.8
ACF
0.6
0.4
0.2
0.0
0
5
10
15
20
Lag Fig. 3. Polio counts ACF plot.
Series y 0.3
Partial ACF
0.2
0.1
0.0
-0.1
5
10
15
20
Lag Fig. 4. Polio counts PACF plot.
Model III: Z t PoissonðlÞ and logitðft Þ ¼ a þ bX t1 þ
b1 t þ b2 cosð2pt=12Þ þ b3 sinð2pt=12Þ þ b4 cosð2pt=6Þ þ b5 sinð2pt=6Þ þ U t . 1000
Model IV: logitðft Þ ¼ a þ bX t1 þ U t , Z t Poissonðlt Þ, and log lt ¼ b0 þ b1 t=1000 þ b2 cosð2pt=12Þ þ b3 sinð2pt=12Þ þ b4 cosð2pt=6Þ þ b5 sinð2pt=6Þ.
ARTICLE IN PRESS H. Zheng, I.V. Basawa / Statistics & Probability Letters 78 (2008) 1–9
9
Table 3
a b s l b0 b1 b2 b3 b4 b5 AIC
Model I
Model II
Model III
Model IV
2.01
5.02
3.65 0.35 4.50 1.03
2.70 0.33 5.00
5.16 0.43 4.61 0.37 0.40 0.28 0.35 551.82
0.29 3.61 0.35 0.39 0.22 0.28 530.12
2.42 2.56 3.29 2.26 3.09 550.77
0.26 3.82 0.36 0.38 0.23 0.28 534.59
Model I is a fixed coefficient INAR(1) process whereas Models II–IV are random coefficient models. Now we use the ML method to estimate the parameters in the above models. The results are summarized in Table 3. The AICs are given in the table for each model as well. Using the AIC criterion, the random coefficient Models II–IV are better than the fixed coefficient Model I. Model II is slightly better than Model IV for this data set. Model IV, however, is a reasonable alternative to model II with the additional advantage of incorporating the dependence of ft on the past observation X t1 . Models I and III have similar AIC. For this data set, we may choose Models II and IV in preference to Models I and III. 6. Summary and conclusions In this paper, we have introduced a new model by combining a thinning model with an observation-driven model. The autoregressive parameter is allowed to vary randomly over time and depends on the past observation. The ergodicity of the process is established. Conditional least squares (CLS) and maximum likelihood (ML) are used to estimate the parameters. In the simulation study, we compare CLS and ML estimators. The simulation results show that the two methods give similar estimates and the ML is slightly better. The model with covariates is applied to a real data set. It is shown that the random coefficient models are a better fit than the fixed coefficient model for this data set. References Al-Osh, M.A., Alzaid, A.A., 1987. First order integer-valued autoregressive (INAR(1)) processes. J. Time Series Anal. 8, 261–275. Alzaid, A.A., Al-Osh, M.A., 1988. First order integer-valued autoregressive (INAR(1)) processes: distributional and regression properties. Statist. Neerlandica 42, 53–61. Billingsley, P., 1961. Statistical Inference for Markov Processes. The University of Chicago Press, Chicago. Davis, R.A., Dunsmuir, W.T.M., Wang, Y., 1999. Modeling time series for count data. In: Ghosh, S. (Ed.), Asymptotics, Nonparametrics and Time Series. Marcel Dekker, New York, pp. 63–114. Fukasawa, T., Basawa, I.V., 2002. Estimation for a class of generalized state-space time series models. Statist. Probab. Lett. 60, 459–473. Guttorp, P., 1991. Statistical Inference for Branching Processes. Wiley, New York. Klimko, L.A., Nelson, P.I., 1978. On conditional least squares estimation for stochastic processes. Ann. Statist. 6, 629–642. MacDonald, I.L., Zucchini, W.Z., 1997. Hidden Markov and Other Models for Discrete-valued Time Series. Chapman & Hall, London. McKenzie, E., 1985a. Contribution to the discussion of Lawrence and Lewis. J. Roy. Statist. Soc. B 47, 187–188. McKenzie, E., 1985b. Some simple models for discrete variate time series. Water Resour. Bull. 21, 645–650. Meyn, S.P., Tweedie, R.L., 1993. Markov Chains and Stochastic Stability. Springer, Berlin. Steutal, F., Van Harn, K., 1979. Discrete analogues of self-decomposability and stability. Ann. Probab. 7, 893–899. Tweedie, R.L., 1975. Sufficient conditions for regularity, recurrence and ergocidicity of Markov processes. Stochastic Process Appl. 3, 385–403. Zeger, S.L., Qaqish, B., 1988. Markov regression models for time series: a quasi-likelihood approach. Biometrics 44, 1019–1031. Zheng, H.T., Basawa, I.V., Datta, S., 2006. Inference for pth order random coefficient integer-valued autoregressive processes. J. Time Series Anal. 27, 411–440. Zheng, H.T., Basawa, I.V., Datta, S., 2007. First order random coefficient integer-valued autoregressive processes. J. Statist. Plann. Inference 137 (1), 212–229.