A Binomial Integer-Valued ARCH Model

Miroslav M. Ristić; Christian H. Weiß; Ana D. Janjić

doi:10.1515/ijb-2015-0051

Artikel Öffentlich zugänglich

A Binomial Integer-Valued ARCH Model

Miroslav M. Ristić , Christian H. Weiß und Ana D. Janjić

Veröffentlicht/Copyright: 5. Dezember 2015

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Manuskript einreichen Informationen für Autor*innen

Aus der Zeitschrift The International Journal of Biostatistics Band 12 Heft 2

Abstract

We present an integer-valued ARCH model which can be used for modeling time series of counts with under-, equi-, or overdispersion. The introduced model has a conditional binomial distribution, and it is shown to be strictly stationary and ergodic. The unknown parameters are estimated by three methods: conditional maximum likelihood, conditional least squares and maximum likelihood type penalty function estimation. The asymptotic distributions of the estimators are derived. A real application of the novel model to epidemic surveillance is briefly discussed. Finally, a generalization of the introduced model is considered by introducing an integer-valued GARCH model.

Keywords: binomial INARCH (p) model; binomial INGARCH (p,q) model; overdispersion; parameter estimation; epidemic surveillance

JEL Classification: 2000 MSC: 62M10

1 Introduction

In this paper, we introduce a model for integer-valued time series with finite range {0,1,…,n}, where n∈N denotes the (known) upper limit. In recent years, integer-valued time series with such a finite range have been widely reported in diverse real-life applications, such as the monitoring of computer pools with n workstations [1, 2], the number of transactions of n companies [3], the number of metapopulations with n patches [4], etc. Integer-valued time series with binomial marginals have been studied by many authors and different approaches have been used to construct them. The first approach is based on the binomial thinning operator “∘” as introduced by Steutel and van Harn [5]. McKenzie [6] defined the binomial AR(1) model as

Xt=α∘Xt−1+β∘(n−Xt−1),t≥1,

where X0 has the binomial Bin(n,p) distribution, ρ∈(max(−p/(1−p),−(1−p)/p),1), α=β+ρ, β=p(1−ρ) and p∈(0,1). All the counting series in “α∘” and “β∘” are mutually independent sequences of independent Bernoulli distributed random variables with parameters α and β, respectively, and the counting series at time t are independent of the random variables {Xs} for all s<t. The binomial AR(1) model is a stationary and ergodic Markov chain with binomial Bin(n,p) marginal distribution. Its autocorrelation function is of the same form as the autocorrelation function of the usual AR(1) process, and it is given by ρ(k)=ρk, k≥0. More properties of the binomial AR(1) model and some estimation issues can be found in Cui and Lund [7], Weiß and Pollett [4] and Weiß and Kim [2, 3]. Weiß [1] extended the binomial AR(1) model to the high-order binomial AR(p) model and defined it as

Xt=∑i=1pDt,i⋅[α∘t+iXt−i+β∘t+i(n−Xt−i)]≡∑i=1pDt,i⋅fi(Xt−i),

where {Dt=(Dt,1,…,Dt,p)} is a sequence of independent random vectors with the multinomial distribution MULT(1,ϕ1,…,ϕk), Dt are independent of Xs and fi(Xs) for all s<t and i=1,2,…,p, and the conditional probabilities P(f1(Xt)=i1,…,fp(Xt)=ip|Xt=xt,lht−1) and P(f1(Xt)=i1,…,fp(Xt)=ip|Xt=xt) are equal, where ht−1 is the process history of all random variables Xs and fj(Xs), j=1,2,…,p. A bivariate extension of the binomial AR(1) model can be found in Scotto et al. [8].

The second approach is based on the use of the hypergeometric thinning operator as introduced by Al-Osh and Alzaid [9]. Al-Osh and Alzaid [9] defined a binomial AR(1) model as

Xt=S(Xt−1)+εt,t≥1,

where {εt} is a sequence of independent and identically distributed random variables, independent of the initial state X0, and the random variable S(X) for given X=x has the hypergeometric distribution with parameters n, x and m, i.e.

PS(X)=k|X=x=xkn−xm−knm,max(0,m−n+x)≤k≤min(x,m).

From the definition of the model, it follows that the random variable S(X) has the binomial distribution with parameters m and p, and the random variable εt has the binomial distribution with parameters n−m and p. Al-Osh and Alzaid [9] showed that their binomial AR(1) model is a stationary Markov chain with binomial Bin(n,p) marginal distribution and autocorrelation function given by ρ(k)=(m/n)k, k≥0.

A third approach for count data time series with a finite range {0,1,…,n} is related to the so-called INGARCH models Ferland et al. [10], i.e. the integer-valued GARCH models. Weiß and Pollett [11] introduced the INARCH(1) model with binomial marginals as a boundary case of the binomial AR(1) processes with density dependent thinning. They considered an integer-valued time series model {Xt}t∈Z such that Xt|Ft−1:Bin(n,αt), t∈Z, where ft−1 is the σ-field generated by the random variables {Xt−k}k≥1, and αt is generated as αt=a0+a1nXt−1, t∈Z, where a0>0 and a1≥0. This model is referred to as the binomial INARCH(1) model. Weiß and Pollett [11] derived some properties of the binomial INARCH(1) model and compared it with the (infinite-range) Poisson INARCH(1) model introduced by Ferland et al. [10]. Such INGARCH models were further investigated and generalized by several authors including Zhu [12, 13]; Xu et al. [14]; Gonçalves et al. [15]. In contrast to the binomial INARCH(1) model, however, all these models are designed for processes with the infinite range N0.

In this paper, we follow the third approach and extend the binomial INARCH(1) model to the binomial INARCH(p) model. In Section 2, we introduce the binomial INARCH(p) model, and we prove strict stationarity and ergodicity of this model. In Section 3, we estimate the unknown parameters by three different estimation methods: conditional maximum likelihood, conditional least squares and maximum likelihood type penalty function estimation. The strong consistency and asymptotic normality of the obtained estimators are derived and discussed. In Section 4, we provide some simulation results to check the performance of the three estimation methods. Section 5 demonstrates that our new model is particularly beneficial in biostatistics. There, we discuss a possible application of the introduced model to a real data set from the field of epidemiology. It might be utilized for epidemic surveillance systems, where approaches based on (infinite-range) count data time series already have been successfully applied [16, 17]. Finally, in Section 6, we introduce the full binomial INGARCH(p,q) model as a generalization of the binomial INARCH(p) model.

2 The binomial INARCH model

In this section, we extend the binomial INARCH(1) model and consider an integer-valued time series model for {Xt}t∈Z given as

(1)Xt|ft−1: Bin(n,αt),∀t∈Z,

where ft−1 is the σ-field generated by the random variables {Xt−k}k≥1, n is a positive integer and αt is generated as

(2)αt=a0+1n∑i=1paiXt−i,∀t∈Z,

where a0>0, ai≥0, for i=1,2,…,p, and p∈{1,2,…}. We suppose that the parameters ai, i=0,1,…,p, satisfy the inequality a0+∑i=1pai<1, which implies that {αt} given by (2) is well-defined, i.e. that αt belongs to the interval (0,1) for all t∈Z. We will say that the time series {Xt}t∈Z given by (1) and (2) is the binomial integer-valued ARCH model and we will denote it as BINARCH(p).

Under the assumption (1), the conditional probability of the random variable Xt for given Ft−1 equals

P(Xt=x|ft−1)=nxαtx(1−αt)n−x,x=0,1,2,…,n.

Then the conditional mean and conditional variance are given as E(Xt|Ft−1)=nαt and Var(Xt|Ft−1)=nαt(1−αt). Thus, we have that Var(Xt|Ft−1)<E(Xt|ft−1). Let us now consider the unconditional mean and unconditional variance. The unconditional mean of the random variable Xt is given by E(Xt)=nE(αt). On the other hand, the unconditional variance of the random variable Xt is given by

(3)Var(Xt)=EVar(Xt|ft−1)+VarE(Xt|ft−1)=nE(αt)−nE(αt2)+n2Var(αt)=n(n−1)Var(αt)+nE(αt)(1−E(αt)).

We shall further investigate the dispersion behaviour of the BINARCH(p) in Example 1.

In the next theorem, we derive some properties of the model given by (1) and (2).

Theorem 1

The BINARCH(p) process given by (1) and (2) is an ergodic, strict and second-order stationary process.

The proof of Theorem 1 is provided by Appendix A.1.

Remark 1. If {Xt}t∈Z is the BINARCH(p) process given by (1) and (2), then

μ≡E(Xt)=na01−∑i=1pai.

In the rest of this section, we will derive and discuss the autocovariance structure of the BINARCH(p) model given by (1) and (2).

Theorem 2

Let {Xt}t∈Z be the BINARCH(p) process given by (1) and (2). The autocovariance function γX(k)=Cov(Xt,Xt−k), k≥0, satisfies the equations

γX(0)=μ(1−μ/n)+(1−1/n)∑i=1paiγX(i),

γX(k)=∑i=1paiγX(|k−i|),k≥1.

The proof of Theorem 2 is provided in Appendix A.2.

Thus, we can see that the autocovariances of the BINARCH(p) model form equations similar to the Yule-Walker equations of the standard AR(p) model. Also, from the above theorem, we obtain that the autocorrelation function ρX(k)≡Corr(Xt,Xt−k) satisfies the equations ρX(k)=∑i=1paiρX(|k−i|), k≥1, and is independent of n and a0. In Figure 1, we present some examples of the autocorrelation function of the BINARCH(p) model, with different values of the parameters ai, i=1,2,…,p, and p∈{2,3,4}.

$Figure 1: Examples of autocorrelation functions of BINARCH(p)$(p)$ model for different values of parameters ai${a_i}$ and p.$

Figure 1:

Examples of autocorrelation functions of BINARCH(p) model for different values of parameters ai and p.

Example 1

Let us consider the special case p=1 in some more detail, i.e. the BINARCH(1) model as in Weiß and Pollett [11]. From Remark 1, we have that μ=na0/(1−a1), while Theorem 2 implies that the autocovariance structure of the time series {Xt}t∈Z is given by γX(k)=a1kγX(0) for k≥0 (exponentially decaying), where

γX(0)=n2a0(1−a0−a1)(1−a1)2[a12+n(1−a12)].

To further investigate the dispersion behaviour of the BINARCH(1) model, let us define the Poisson index of dispersion by IPois:=Var(Xt)/E(Xt). Then we say that we have underdispersion if IPois<1, equidispersion if IPois=1, and overdispersion if IPois>1 (all with respect to a Poisson distribution). For the BINARCH(1) model, we have that

IPois=n(1−a0−a1)(1−a1)[a12+n(1−a12)].

Thus, when a0=n−1na12(1−a1), we obtain equidispersion, when a0<n−1na12(1−a1), we obtain overdispersion, and when a0>n−1na12(1−a1), we have underdispersion. The attainable range of IPois is bounded by 0 (if a0+a1 is close to 1) and by

IPois<na12+n(1−a12)<n.

Instead of considering over- and underdispersion with respect to a Poisson distribution as before, we might also compare the unconditional variance-mean behavior with that of a binomial distribution with population parameter n. For this purpose, let us investigate the so-called binomial index of dispersion, defined by

IBin:=IBin(n,μ,σ2):=σ2μ(1−μ/n)

for a random variable X with range {0,…,n}, mean μ and variance σ2. For the case of the BINARCH(p) model with μ:=E(Xt) and σ2:=Var(Xt), it follows from eq. (3) that

IBin=n(n−1)Var(αt)+μ(1−μ/n)μ(1−μ/n)=1+n(n−1)Var(αt)μ(1−μ/n)=1+(1−1/n)∑i=1paiγX(i)μ(1−μ/n),

i.e. for any BINARCH(p) model, we have IBin>1. While the conditional distribution of Xt|Ft−1 satisfies IBin=1 (conditional binomial distribution), the unconditional distribution always exhibits extra-binomial variation (overdispersion with respect to a binomial distribution).

For the example of the BINARCH(1) model, also see Weiß and Pollett [11], we have

IBin=1+(n−1)a12n(1−a12)+a12.

3 Estimation of parameters

In this section, we consider the estimation of the unknown parameters θ=(a0,a1,…,ap)⊤ of the BINARCH(p) model, while we suppose that the parameter n of the conditional binomial distribution is known. We consider three estimation methods: conditional maximum likelihood estimation, conditional least squares estimation and maximum likelihood type penalty function estimation from Tjøstheim [18]. The second and the third estimation approach are based on the minimization of an objective function, so both approaches can be understood as penalty function approaches [18].The penalty function for the conditional least squares estimation as given by (5) accumulates squared deviations, while the penalty function (7) for the maximum likelihood type penalty function estimation is motivated by the Gaussian log-likelihood function. Also, these penalty functions are chosen in a way to provide some asymptotic properties of the resulting estimators.

We suppose that X1, X2, …, XN are the observations generated by the BINARCH(p) process {Xt}t∈Z, where N∈N represents the size of the sample. Here, the parameter n (upper limit of the range) is considered as a known quantity.

3.1 Conditional maximum likelihood estimation

From the definition of the BINARCH(p) model, we obtain that the conditional log-likelihood function is given by

(4)l(θ)=∑t=p+1NlognXt+Xtlogαt+(n−Xt)log(1−αt).

Since the BINARCH(p) process {Xt}t∈Z can be represented equivalently as the finite Markov chain {Xt}t∈Z (see Section 2), it is possible to apply the results in Billingsley [19] to investigate the properties of the conditional maximum likelihood (CML) estimators of the parameter vector θ:=(a0,a1,…,ap)⊤.

Theorem 3

There exists a consistent CML estimator of θ, maximizing (4), that is also asymptotically normally distributed.

The proof of Theorem 3 is provided by Appendix A.3.

3.2 Conditional least squares estimation

While the CML approach discussed in the previous section makes use of the complete conditional distribution, we shall now derive the semiparametric conditional least squares (CLS) estimators of the BINARCH(p) model {Xt}t∈Z and discuss their asymptotic properties. Let θ:=(a0,a1,…,ap)⊤ be the vector of the unknown parameters, and let Zt:=(1,Xt−1,…,Xt−p)⊤ and Wt:=(n,Xt−1,…,Xt−p)⊤ be the vectors of the observations. The CLS estimates of the vector θ:=(a0,a1,…,ap)⊤ are obtained by minimizing the function

(5)S(θ)=∑t=p+1NXt−W⊤θ2=∑t=p+1NXt−na0−∑i=1paiXt−i2

with respect to the vector θ, and they are given as

(6)θˆCLS=∑t=p+1NZtWt⊤−1∑t=p+1NXtZt.

Now we will derive the asymptotic properties of the CLS estimators θˆCLS. First, we start with the consistency of these estimators.

Theorem 4

The CLS estimators θˆCLS given by (6) are strongly consistent estimators of the unknown parameter θ.

The proof of Theorem 4 is provided by Appendix A.4.

The asymptotic distribution of the CLS estimators redθ^CLS given by (6) follows from the following theorem, the proof of which is provided by Appendix A.5.

Theorem 5

If the CLS estimators redθ^CLS are given by (6), then

N(redθ^CLS−θ)→dN(0,U−1RU−1), N→∞,

where U=E(WtWt⊤) and R=nEαt(1−αt)WtWt⊤.

Closed-forms expressions for the matrices U and R for higher order p are very cumbersome. Because of that, Appendix B provides closed-form expressions for the matrices U and R for the special case of the BINARCH(1) model, which may be used, in turn, to derive approximate standard errors of the CLS estimates.

3.3 Maximum likelihood type penalty function

This subsection is dedicated again to maximum likelihood estimators of the BINARCH(p) model. Instead of the standard approach to CML estimators, we shall use the maximum likelihood type penalty (MLTP) function (T1986) of the observed model given by

(7)L(θ)=∑t=p+1Nlog(n)+log(αt)+log(1−αt)+(Xt−nαt)2nαt(1−αt)=∑t=p+1Nϕt.

Tjøstheim [18] gave two motivations for using a penalty term in the conditional log-likelihood function. First, in the case of a conditional Gaussian process, the maximum likelihood penalty function L(θ) coincides with the conditional log-likelihood function of this process except a multiplicative constant. The second motivation is that ϕt has the marginal property, which can be used to derive asymptotic properties of the considered estimators. Differentiating the function L(θ) with respect to the parameters ai, i=0,1,…,p, we obtain the estimators of the unknown parameters as the solutions of the nonlinear system of the equations

∑t=p+1N∂ϕt∂ai=∑t=p+1N−Xt2+nαt+2αtXt2−3nαt2+n2αt2−2nαt2Xt+2nαt3nαt2(1−αt)2⋅∂αt∂ai=0,

for i=0,1,2,…,p. Now we will focus on deriving consistency and asymptotic properties of the estimators obtained by using above function. The notation introduced in previous section is retained.

Theorem 6

The MLTP estimators obtained by minimizing (7) are strongly consistent estimators of unknown parameter θ.

The proof of Theorem 6 is provided by Appendix A.6.

Finally, the asymptotic normality of the estimators obtained by minimizing (7) is established by the following theorem, the proof of which is provided by Appendix A.7.

Theorem 7

If {θˆMLTP} are the MLTP estimators obtained by minimizing (7), then

N(θˆMLTP−θ)→dN(0,V−1+V−1SV−1),N→∞,

where

V=12n2E1+2(n−2)αt−2(n−2)αt2αt2(1−αt)2WtWt⊤

and

S=14n4E−Xt2+nαt+2αtXt2−3nαt2+n2αt2−2nαt2Xt+2nαt3αt2(1−αt)22WtWt⊤−V.

4 Simulation study

In this section, we provide some results from a simulation study to check the finite-sample performance of the three estimation methods considered in the previous section. We simulated samples of size 500, and the number of replications is m=10,000. We perform the estimation for the subsamples consisting of the first 50, 100, 200 and 500 elements, thus considering subsamples of four different sizes N=50, N=100, N=200, and N=500. We provide the estimation of the parameters of BINARCH(p) models for p=1 and p=2. For the case p=1, we conduct simulations for the following cases: (1) (n,a0,a1)=(5,0.8,0.1); (2) (n,a0,a1)=(5,0.6,0.1); (3) (n,a0,a1)=(5,0.3,0.6); (4) (n,a0,a1)=(5,0.1,0.8); (5) (n,a0,a1)=(5,0.05,0.9). We also consider four cases for p=2: (6) (n,a0,a1,a2)=(5,0.6,0.2,0.1); (7) (n,a0,a1,a2)=(5,0.3,0.5,0.1); (8) (n,a0,a1,a2)=(5,0.1,0.1,0.75); and (9) (n,a0,a1,a2)=(5,0.1,0.75,0.1). For each choice of the true values of the parameters, we compute the estimates according to the three considered estimation methods for each generated subsample. The mean of the estimates and the mean absolute deviation errors provided within parentheses are given in Tables 1 and 2. From these tables, we can conclude that all three estimation methods give good estimates that quickly approach their true values as the size of the sample increases. Also, we can see that the mean absolute deviations are small and they decrease as the size of the sample increases.

Table 1:

Mean of estimates and mean absolute deviation errors in parentheses for BINARCH(1) model.

Model	N	a0CLS	a1CLS	a0CML	a1CML	a0MLTP	a1MLTP
(1)	50	0.8263(0.1076)	0.0710(0.1178)	0.8009(0.0851)	0.1016(0.0903)	0.7935(0.0868)	0.1067(0.0972)
	100	0.8137(0.0758)	0.0848(0.0833)	0.8039(0.0669)	0.0973(0.0719)	0.7995(0.0695)	0.1002(0.0777)
	200	0.8073(0.0541)	0.0921(0.0591)	0.8044(0.0514)	0.0959(0.0556)	0.8025(0.0545)	0.0973(0.0602)
	500	0.8026(0.0343)	0.0972(0.0376)	0.8023(0.0340)	0.0976(0.0372)	0.8022(0.0369)	0.0979(0.0405)
(2)	50	0.6176(0.0807)	0.0741(0.1147)	0.5965(0.0598)	0.1014(0.0878)	0.5983(0.0648)	0.1047(0.0911)
	100	0.6089(0.0560)	0.0867(0.0797)	0.5996(0.0468)	0.0973(0.0692)	0.6004(0.0490)	0.0989(0.0708)
	200	0.6043(0.0401)	0.0934(0.0569)	0.6012(0.0371)	0.0965(0.0538)	0.6012(0.0378)	0.0972(0.0546)
	500	0.6018(0.0251)	0.0974(0.0356)	0.6015(0.0249)	0.0976(0.0353)	0.6014(0.0252)	0.0979(0.0359)
(3)	50	0.3463(0.0898)	0.5412(0.1076)	0.3432(0.0867)	0.5461(0.1031)	0.3370(0.0891)	0.5570(0.1060)
	100	0.3230(0.0594)	0.5711(0.0703)	0.3213(0.0576)	0.5739(0.0674)	0.3178(0.0595)	0.5799(0.0703)
	200	0.3111(0.0410)	0.5861(0.0485)	0.3103(0.0400)	0.5875(0.0467)	0.3085(0.0415)	0.5904(0.0491)
	500	0.3046(0.0255)	0.5941(0.0304)	0.3043(0.0247)	0.5947(0.0290)	0.3036(0.0256)	0.5959(0.0304)
(4)	50	0.1303(0.0494)	0.7396(0.0849)	0.1302(0.0485)	0.7401(0.0825)	0.1229(0.0486)	0.7553(0.0797)
	100	0.1135(0.0293)	0.7731(0.0490)	0.1131(0.0281)	0.7737(0.0468)	0.1083(0.0291)	0.7833(0.0473)
	200	0.1064(0.0193)	0.7877(0.0313)	0.1062(0.0183)	0.7880(0.0296)	0.1035(0.0195)	0.7932(0.0310)
	500	0.1025(0.0116)	0.7952(0.0184)	0.1025(0.0109)	0.7954(0.0173)	0.1013(0.0119)	0.7976(0.0187)
(5)	50	0.0922(0.0540)	0.8151(0.0976)	0.0907(0.0524)	0.8191(0.0923)	0.0886(0.0531)	0.8265(0.0884)
	100	0.0670(0.0272)	0.8663(0.0454)	0.0653(0.0249)	0.8697(0.0416)	0.0621(0.0249)	0.8765(0.0398)
	200	0.0571(0.0156)	0.8860(0.0243)	0.0556(0.0134)	0.8888(0.0212)	0.0535(0.0136)	0.8928(0.0211)
	500	0.0526(0.0089)	0.8950(0.0131)	0.0519(0.0074)	0.8964(0.0112)	0.0511(0.0080)	0.8979(0.0117)

Table 2:

Mean of estimates and mean absolute deviation errors in parentheses for BINARCH(2) model.

Model	N	a0CLS	a1CLS	a2CLS	a0CML	a1CML	a2CML	a0MLTP	a1MLTP	a2MLTP
(6)	50	0.6700(0.1456)	0.1706(0.1290)	0.0490(0.1246)	0.5905(0.0899)	0.2092(0.1151)	0.1022(0.0956)	0.5843(0.1007)	0.2114(0.1216)	0.1077(0.1016)
	100	0.6368(0.0998)	0.1845(0.0886)	0.0733(0.0862)	0.6007(0.0744)	0.2055(0.0883)	0.0933(0.0742)	0.5972(0.0820)	0.2059(0.0931)	0.0977(0.0793)
	200	0.6199(0.0686)	0.1918(0.0620)	0.0854(0.0601)	0.6055(0.0576)	0.2024(0.0640)	0.0911(0.0554)	0.6036(0.0628)	0.2027(0.0676)	0.0932(0.0594)
	500	0.6082(0.0428)	0.1970(0.0387)	0.0937(0.0380)	0.6053(0.0399)	0.1996(0.0392)	0.0943(0.0373)	0.6046(0.0432)	0.2003(0.0420)	0.0946(0.0406)
(7)	50	0.3676(0.1102)	0.4643(0.1267)	0.0489(0.1230)	0.3554(0.1001)	0.4449(0.1142)	0.0856(0.0864)	0.3476(0.1039)	0.4551(0.1199)	0.0901(0.0923)
	100	0.3323(0.0709)	0.4846(0.0877)	0.0739(0.0862)	0.3295(0.0695)	0.4755(0.0813)	0.0887(0.0703)	0.3249(0.0718)	0.4805(0.0856)	0.0912(0.0749)
	200	0.3167(0.0483)	0.4917(0.0612)	0.0869(0.0598)	0.3161(0.0477)	0.4890(0.0581)	0.0914(0.0540)	0.3135(0.0495)	0.4920(0.0611)	0.0925(0.0574)
	500	0.3061(0.0296)	0.4972(0.0376)	0.0949(0.0380)	0.3060(0.0289)	0.4971(0.0365)	0.0954(0.0364)	0.3052(0.0301)	0.4986(0.0384)	0.0953(0.0389)
(8)	50	0.2015(0.1160)	0.0549(0.1096)	0.6396(0.1270)	0.2101(0.1218)	0.0826(0.0805)	0.6087(0.1524)	0.2108(0.1260)	0.0840(0.0859)	0.6163(0.1500)
	100	0.1495(0.0652)	0.0824(0.0672)	0.6957(0.0721)	0.1600(0.0735)	0.0877(0.0576)	0.6832(0.0819)	0.1562(0.0745)	0.0875(0.0626)	0.6914(0.0813)
	200	0.1240(0.0387)	0.0923(0.0440)	0.7237(0.0440)	0.1279(0.0412)	0.0929(0.0398)	0.7211(0.0456)	0.1250(0.0425)	0.0926(0.0437)	0.7269(0.0465)
	500	0.1087(0.0210)	0.0976(0.0262)	0.7402(0.0244)	0.1087(0.0199)	0.0977(0.0238)	0.7409(0.0234)	0.1073(0.0212)	0.0975(0.0263)	0.7436(0.0248)
(9)	50	0.1691(0.0883)	0.7036(0.1402)	0.0493(0.1310)	0.1649(0.0833)	0.6724(0.1228)	0.0895(0.0906)	0.1575(0.0830)	0.6854(0.1260)	0.0892(0.0963)
	100	0.1316(0.0494)	0.7312(0.0910)	0.0754(0.0907)	0.1304(0.0474)	0.7206(0.0776)	0.0892(0.0711)	0.1248(0.0482)	0.7313(0.0826)	0.0881(0.0776)
	200	0.1145(0.0294)	0.7417(0.0617)	0.0884(0.0630)	0.1147(0.0286)	0.7388(0.0545)	0.0921(0.0551)	0.1115(0.0302)	0.7454(0.0599)	0.0908(0.0611)
	500	0.1058(0.0169)	0.7466(0.0382)	0.0955(0.0398)	0.1058(0.0161)	0.7466(0.0355)	0.0959(0.0364)	0.1043(0.0173)	0.7499(0.0395)	0.0949(0.0412)

In the case of the BINARCH(1) model, the best results are usually obtained by CML estimation. If a1 is large and a0 is small (cases 4 and 5), then MLTP estimation provides the best results for a1 for samples of small size, while otherwise, the best results are provided by CML estimation.

In the case of the BINARCH(2) model, CML and CLS estimation give the best results, while MLTP estimation was superior only in one case. The CLS estimation is very sensitive to small true values of the parameters. In this case, if some CLS estimates are smaller than 0, we set these estimates to take the value 0. For this method of estimation, it is also possible that the sum of the estimates is equal or greater than 1. On the other hand, many statistical packages allow to minimize and maximize functions with boundary conditions, which imply that the other two estimation methods give estimates that never take values outside (0,1), and that their sum is always less than 1.

The simplest method for use is the CLS estimation method, and the corresponding estimates are used as the starting values for other two methods. An interesting conclusion from our simulations is that, as expected, the efficiency of the CML estimates with respect to the CLS estimates is high.

5 Real-data example

In this section, we present a possible application of the novel BINARCH(p) model in the field of biostatistics. As a real-data set, we consider the infection counts as previously discussed by Weiß and Pollett [11]. This data set, taken from the “SurvStat” data base of the Robert-Koch-Institut [20], contains the number of districts with new cases of hantavirus infections per week in the year 2011 (N=52 counts) reported in n=38 Germany’s districts, i.e. the counts express the regional spread of the hantavirus infections. More than 200,000 infections by hantavirus are reported from all over the world per year. For most European countries, rising numbers of cases are reported, and in Germany, hantaviral infections meanwhile became the most common endemic rodent-borne human illness [21, 22]. The hantavirus mainly causes two diseases: The hantavirus cardiopulmonary syndrome (HCPS) with case-fatality rates >35% is mainly reported from North, Central, and South America, while in Europe (and Asia), the less severe haemorrhagic fever with renal syndrome (HFRS; case-fatality rates >10%) is usually observed [21, 23]. Among the viral species causing HFRS in Central Europe are the Puumala virus (carried by the red bank vole) and the Dobrava virus (carried by the striped field mouse), see Schilling et al. [24]; Heyman et al. [22] for further details.

Let us return to the particular infection counts for 2011 in Germany. The minimum number of districts with new cases is 0 and the maximum number is 11. The sample mean is 4.173 and the sample variance is 7.793. As a result, the binomial index of dispersion is 2.098, which indicates overdispersion with respect to a binomial distribution (extra-binomial variation), also see Example 1. Therefore, it is plausible that Weiß and Pollett [11] found the binomial INARCH(1) model to be superior to the binomial AR(1) model. The empirical autocorrelation and partial autocorrelation function (ACF and PACF, respectively) are given in Figure 2. Inspecting the plot of the PACF in more detail, however, an autoregressive model of order >1 appears to be reasonable. Hence, we shall now apply the novel BINARCH(p) model to the data, and we shall also compare it to the binomial AR(p) model (BINAR(p)) as introduced in Weiß [1], see Section 1 for further details.

Figure 2:

Infection counts: (a) autocorrelation function; (b) partial autocorrelation function.

First, we need to determine the order p. The plot of the PACF indicates that p may take values ≤3. According to this, we consider models up to order 4 in our study. Next, we estimate the unknown parameters of both types of models by a maximum likelihood method. Finally, we compare the performance of these models. Often, one of the information criteria AIC or BIC is applied in this context. But since we have to use a conditional ML approach, where the number of terms involved in the log-likelihood function (4), namely N−p, decreases with increasing p, these information criteria are misleading in the present situation. Therefore, we shall evaluate the performance by comparing the root mean square errors of the models. Let aˆ0, aˆ1, …, aˆp represent the CML estimates of the parameters of the BINARCH(p) model. Then the root mean square error of the BINARCH(p) model is given by

RMS=1N−p∑t=p+1NXt−naˆ0−∑i=1paˆiXt−i2.

On the other hand, let πˆ, ρˆ, ϕˆ1, …, ϕˆp represent the CML estimates of the parameters of the BINAR(p) model. The estimates ϕˆi, i=1,2,…,p satisfy the condition ∑i=1pϕˆi=1. Because of that, we only estimate the parameters ϕ1, ϕ2, …, ϕp−1, while the estimate of the last parameter ϕp is obtained as ϕˆp=1−∑i=1p−1ϕˆi. This implies that the estimated standard error for ϕˆp will always be equal to 0. Now, the root mean square error of the BINAR(p) model is given by

RMS=1N−p∑t=p+1NXt−nπˆ(1−ρˆ)−ρˆ∑i=1pϕˆiXt−i2.

In Table 3, we give the CML estimates of the unknown parameters with the estimated standard errors in parentheses and the respective root mean square errors. Also, we give the negative log-likelihood functions for each fitted model. It becomes clear that the BINARCH models are always superior to their BINAR counterpart, which is plausible in view of the extra-binomial variation observed in the data. With increasing autoregressive order p (and hence increasing number of parameters), some of the parameter estimates are not significant anymore. This is mainly a problem of the sample size, which equals only N=52 for the analyzed data set. The smallest root mean square error is obtained for the BINARCH(3) model, which indicates that our novel model with p=3 is most appropriate for the infection counts. Recalling that the counts are sampled on a weekly base, the order 3 is indeed plausible, since both for the Puumala virus and the Dobrava virus, the incubation period is known to strongly vary, namely between 2 and 4 weeks [22]. Finally, in Figure 3, we present plots of the observed and expected values for each considered model, which highlight that the higher-order BINARCH models are best suited to adapt to the higher level observed in the second half of 2011.

Table 3:

Maximum likelihood estimates of the parameters of the BINARCH(p) and BINAR(p) models with the corresponding root mean square errors and negative log-likelihood functions.

Model			Estimates				-log L	RMS
BINAR(1)	πˆ	ρˆ	ϕˆ1
	0.1151	0.5353	1.0000				109.3891	2.1226
	(.01337)	(0.0707)	(0.0000)
BINARCH(1)	aˆ0	aˆ1
	0.0303	0.7476					103.6976	0.0981
	(.0111)	(0.1085)
BINAR(2)	πˆ	ρˆ	ϕˆ1	ϕˆ2
	0.1176	0.6789	0.4283	0.5717			101.3702	1.9140
	(.0193)	(0.0688)	(0.1690)	(0.0000)
BINARCH(2)	aˆ0	aˆ1	aˆ2
	0.0142	0.4685	0.4333				95.9150	1.8653
	(.0109)	(0.1396)	(0.1405)
BINAR(3)	πˆ	ρˆ	ϕˆ1	ϕˆ2	ϕˆ3
	0.1226	0.7556	0.2841	0.4478	0.2681		96.1284	1.8208
	(.0258)	(0.0639)	(0.1644)	(0.1634)	(0.0000)
BINARCH(3)	aˆ0	aˆ1	aˆ2	aˆ3
	0.0110	0.3500	0.3216	0.2791			92.3401	1.7697
	(0.0108)	(0.1591)	(0.1506)	(0.1534)
BINAR(4)	πˆ	ρˆ	ϕˆ1	ϕˆ2	ϕˆ3	ϕˆ4
	0.1235	0.7622	0.2909	0.4159	0.2594	0.0338	94.4509	1.8359
	(0.0269)	(0.0657)	(0.1672)	(0.1896)	(0.1552)	(0.0000)
BINARCH(4)	aˆ0	aˆ1	aˆ2	aˆ3	aˆ4
	0.0105	0.3599	0.3045	0.2941	0.000002		90.6773	1.7861
	(.0110)	(0.1644)	(0.1660)	(0.1620)	(0.1692)

Figure 3:

Plots of observed and expected values for each considered model.

To further check the adequacy of the BINARCH(3) model, we use the parametric bootstrap based on the fitted model introduced by Tsay [25], which was also considered in Jung and Tremayne [26]; Weiß [27]. For parameter values a0=0.0110, a1=0.3500, a2=0.3216 and a3=0.2791, we simulate 10,000 samples of size 52 from the BINARCH(3) model. For each simulated sample, we compute the sample ACF, and for each fixed lag, we derive the 2.5% and 97.5% quantiles. By using these quantiles, we draw the bootstrap confidence intervals in Figure 4. From this graph, we can conclude that the BINARCH(3) model adequately describes the autocorrelation structure of the infection counts. Certainly, the bootstrap confidence intervals are rather wide, but this was to be expected since the sample size is small in our data example.

Figure 4:

ACF for infection counts with 95% bootstrap confidence intervals.

If we consider the other two estimation methods, the CLS method and the MLTP method, we obtain the following results for the BINARCH(3) model. The CLS estimates are aˆ0=0.0094, aˆ1=0.2713, aˆ2=0.3379, and aˆ3=0.3277, with the RMS=1.7634. We can see that this method gives little smaller RMS than the CML method. On the other hand, the MLTP estimates are aˆ0=0.0032, aˆ1=0.4298, aˆ2=0.2995, and aˆ3=0.2596, with RMS=1.7917. This method gives the largest RMS.

6 The binomial INGARCH model

In this section, we generalize the binomial INARCH(p) model and consider an integer-valued time series model for {Xt}t∈Z following (1) with ft−1 and n defined as in Section 2, while αt is generated as

(8)αt=a0+1n∑i=1paiXt−i+∑j=1qbjαt−j,∀t∈Z,

where a0>0, a1≥0, …, ap≥0, b1≥0, …, bq≥0, p∈{1,2,…} and q∈{0,1,2,…}. Similarly as in the case of the BINARCH(p) model, we suppose that the parameters ai and bj, i=0,1,…,p, and j=1,2,…,q, satisfy the inequality a0+∑i=1pai+∑j=1qbj<1, which implies that {αt} given by (8) is well-defined. We will say that the time series {Xt}t∈Z given by (1) and (8) is the binomial integer-valued GARCH model and we will denote it as BINGARCH(p,q).

Obviously, for q=0, we obtain BINARCH(p) model from before. Another important special case is p=q=1, the BINGARCH(1,1) model, which is also an instance of the observation-driven models as introduced by Davis and Liu [28], see Example 3 therein.

In the next theorem, we show that the model given by (1) and (8) is a first-order stationary process.

Theorem 8

The BINGARCH(p,q) process given by (1) and (8) is a first-order stationary process, where

μ≡E(Xt)=na01−∑i=1pai−∑j=1qbj.

The proof of Theorem 8 is provided in Appendix A.8. For the particular case of a BINGARCH(1,1) process, Davis and Liu [28] even established strict stationarity. Furthermore, also consistency and asymptotic normality of the maximum likelihood estimators are proven in that work.

Now, we derive and discuss the autocovariance structure of the BINGARCH(p,q) model given by (1) and (8).

Theorem 9

Let {Xt}t∈Z be a stationary BINGARCH(p,q) process given by (1) and (8). The autocovariance functions γX(k)=Cov(Xt,Xt−k) and γα(k)=Cov(αt,αt−k) satisfy the following equations:

γX(k)=∑i=1paiγX(|k−i|)+∑j=1min(k−1,q)bjγX(k−j)+n2∑j=kqbjγα(j−k),k≥1,

γα(k)=∑i=1min(k,p)aiγα(k−i)+1n2∑i=k+1paiγX(i−k)+∑j=1qbjγα(|k−j|),k≥0.

The proof of Theorem 9 is provided in Appendix A.9.

Example 2

For the stationary BINGARCH(1,1) model [28], we derive the autocovariance structure explicitly. From Theorem 9, we have that the autocovariance structure of the process {Xt}t∈Z is given by

(9)γX(1)=a1γX(0)+n2b1γα(0),

γX(k)=(a1+b1)k−1γX(1),k≥2.

On the other hand, the autocovariance structure of the process {αt}t∈Z is given by

(10)γα(0)=a1n2γX(1)+b1γα(1),

(11)γα(k)=(a1+b1)kγα(0),k≥1.

Substituting eqs (9) and (11) for k=1 into eq. (10), we obtain that

a12n2γX(0)+(2a1+b1)b1−1γα(0)=0.

Also from eq. (3), we have that

γX(0)−n(n−1)γα(0)=na0(1−a0−a1−b1)(1−a1−b1)2.

Solving the last two equations with respect to γX(0) and γα(0), we obtain that the variances of the random variables Xt and αt are given as, respectively,

Var(Xt)=n2a0(1−a0−a1−b1)(1−2a1b1−b12)(1−a1−b1)2[a12+n(1−(a1+b1)2)],

Var(αt)=a0a12(1−a0−a1−b1)(1−a1−b1)2[a12+n(1−(a1+b1)2)].

Note that these variances are well-defined. Now the autocovariance and the autocorrelation functions of the process {Xt}t∈Z are given as, respectively,

γX(k)=n2a0a1(1−a0−a1−b1)(1−a1b1−b12)(1−a1−b1)2[a12+n(1−(a1+b1)2)](a1+b1)k−1,k≥1,

ρX(k)=a1(1−a1b1−b12)1−2a1b1−b12(a1+b1)k−1,k≥1.

From the last two equations, we can conclude that the autocovariance and autocorrelation functions are positive and exponentially decaying in time.

Remark 2. Using the results from the Example 2, we can show that the BINGARCH(1,1) model can handle underdispersion, equidispersion or overdispersion with respect to a Poisson distribution. We have that

IPois=n(1−a0−a1−b1)(1−2a1b1−b12)(1−a1−b1)[a12+n(1−(a1+b1)2)].

When the sum a0+a1+b1 is close to 1, we have that IPois is close to 0. Since 1−a0−a1−b1<1−a1−b1, we obtain that

IPois<n(1−2a1b1−b12)a12+n(1−(a1+b1)2)<n,

since 1−2a1b1−b12<a12+n(1−(a1+b1)2). So the range of IPois is bounded.

For example, when (a0,a1,b1,n)=(0.1,0.493686,0.1,10), we obtain the equidispersion, when (a0,a1,b1,n)=(0.52,0.42,0.05,10), we obtain that IPois=0.0226 (underdispersion), and when (a0,a1,b1,n)=(0.001,0.99,0.001,10) we obtain that IPois=7.6523 (overdispersion).

Expressing the dispersion behaviour with respect to a binomial distribution, in contrast, we get

IBin=1+(n−1)Var(αt)μ/n(1−μ/n)=1+(n−1)a12n(1−(a1+b1)2)+a12,

i.e. we always have extra-binomial variation.

7 Conclusions

The BINARCH approach offers a way of modelling time series of counts with a finite range which exhibit extra-binomial variation. We introduced an extension of the basic BINARCH(1) model to arbitrary orders p, with a generalization to a full BINGARCH(p,q) model. After having discussed stochastic properties of these models, we analyzed three approaches for parameter estimation, and we successfully applied our novel model to a time series of counts of infections by the hantavirus.

Acknowledgments

The authors are very grateful to the Referees and to Professor Dimitris Karlis for their valuable suggestions and comments, which greatly improved this manuscript.

Conflict of Interest: The authors have declared no conflict of interest.

Appendix A Appendix A Proofs

Appendix A.1 A.1 Proof of Theorem 1

The sequence of random variables {Xt}t∈Z is a pth order Markov process and, hence, the vector-valued process

Xt:=(Xt,…,Xt−p+1)⊤for t∈Z

constitutes a finite 1st order Markov process, i.e. a finite Markov chain. Its 1-step-ahead transition probabilities are

pk|l(a):=P(Xt=k|Xt−1=l)=δk2l1⋯δkplp−1⋅P(Xt=k1|Xt−1=l)

(12)=δk2l1⋯δkplp−1⋅nk1α0k1(1−α0)n−k1with α0=a0+1n∑i=1paili.

Let ki=(ki,ki+1,…,ki+p−1)⊤. Then the p-step-ahead transition probabilities are

P(Xt=k0|Xt−p=kp)=∏j=0p−1P(Xt−j=kj|Xt−j−1=kj+1,…,Xt−j−p=kj+p)

=∏j=0p−1nkjαjkj(1−αj)n−kjwith αj=a0+1n∑i=1pailj+i,

which are truly larger than 0 because of a0>0. Hence, the finite Markov chain {Xt}t∈Z is primitive, implying that it is also irreducible and aperiodic, and therefore ergodic with a unique stationary distribution [29]. Since the range of Xt is finite, any moments exist, and the strict stationarity of {Xt}t∈Z also implies its second-order stationarity. □

Appendix A.2 A.2 Proof of Theorem 2

Let k≥1. From the definition of the BINARCH(p) process {Xt}t∈Z, we have that

(13)γX(k)=E(Xt−μ)(Xt−k−μ)=E(Xt−k−μ)E(Xt|Ft−1)−μ=Cov(nαt,Xt−k)=nCov(αt,Xt−k)=nCova0+1n∑i=1paiXt−i,Xt−k=∑i=1paiγX(|k−i|).

Let us consider now the variance of the random variable Xt. According to (3) and the fact that E(αt)=μ/n, we need the variance of the random variable αt for derivation of Var(Xt). Thus, we have from (2) and (13) that

Var(αt)=Vara0+1n∑i=1paiXt−i=1n2∑i=1pai∑j=1pajγX(|i−j|)=1n2∑i=1paiγX(i).

Replacing the last expression in (3), we obtain the expression for γX(0) □.

Appendix A.3 A.3 Proof of Theorem 3

To prove Theorem 3, we have to check if Condition 5.1 in Billingsley [19] is satisfied. First, the transition probabilities pk|l(θ) from (12) are truly positive for each (k,l)∈D, where

D={(k,l)∈{0,…,n}p×2|k2=l1,…,kp=lp−1}.

Obviously, the set D is independent of θ. Secondly, the non-zero probabilities pk|l(θ) are polynomials in θ and, thus, continuously differentiable in θ up to any order. Also, see above, the transition matrix is primitive. Defining d:=|D|=(n+1)p+1, it suffices to check if the d×(p+1)-matrix

(∂pk|l(θ)/∂au)with (k,l)∈D and u=0,…,p

has rank p+1 for all θ. For this purpose, consider the submatrix corresponding to

k=(n,l1,…,lp−1)⊤, and l=0 or l=(n,0,…,0)⊤ or …or l=(0,…,0,n)⊤,

i.e. the submatrix corresponding to the transition probabilities

a0n,(a0+a1)n,…,(a0+ap)n.

This submatrix equals

na0n−10⋯0n(a0+a1)n−1n(a0+a1)n−10⋮⋱0n(a0+a1)n−10n(a0+ap)n−1,

the determinant of which is ≠0 and, thus, has full rank p+1.

So Condition 5.1 in Billingsley [19] is satisfied, and there exists a consistent CML estimator of θ that is also asymptotically normally distributed [19, Theorems 2.1 and 2.2].

Appendix A.4 A.4 Proof of Theorem 4

We will prove this theorem by using Theorem 3.1 in Tjøstheim [18]. Let us show that all the conditions C1–C3 of this theorem are satisfied. According to the results of the previous section, the BINARCH(p) model {Xt}t∈Z is strictly stationary ergodic process with finite moments. Next, the function gt(θ)=E(Xt|Ft−1)=na0+∑i=1paiXt−i is almost surely three times continuously differentiable with respect to the vector θ in an open set, which contains true vector θ0. The condition C1 is trivially satisfied since

E∂gt(θ)∂a02=n2;E∂gt(θ)∂ai2=E(Xt−i2)<∞,i=1,2,…,p;E∂2gt(θ)∂ai∂aj2=0,

for i,j∈{0,1,…,}. In a similar way, we can show that the condition C3 is satisfied. Now we will prove that the condition C2 is satisfied, i.e. that all the components of the vector ∂gt(θ)/∂θ=(n,Xt−1,…, Xt−p)⊤ are linearly independent. This means that we need to prove that if d0, d1, …, dp are arbitrary real numbers such that

(14)E∑i=0pdi∂gt(θ)∂ai2=End0+∑i=1pdiXt−i2=0,

then d0=d1=…=dp=0. We will prove that the condition C2 is satisfied following the proof of Lemma 6 [30]. Condition (14) implies that nd0+∑i=1pdiXt−i=a.s.0. Now, for simplicity, let us suppose that d1≠0. Then Xt−1=−nd0d1−∑i=2pdid1Xt−i and

P(Xt−1=xt−1|Xt−i=xt−i, i≥2)=1(nd0+∑i=1pdixt−i=0)≠(nxt−1)αt−1xt−1(1−αt−1)n−xt−1,

since Xt−1 for given Xt−i, i≥2, has the binomial distribution with parameters n and αt−1. Thus, we obtain a contradiction which implies that d1=0. In a similar way, we can show that all di=0, i=0,1,…,p.

Thus all the conditions of Theorem 3.1 in Tjøstheim [18] are satisfied, which implies that the CLS estimators θˆCLS given by (6) are strongly consistent estimators of the unknown parameter θ□.

Appendix A.5 A.5 Proof of Theorem 5

We only need to prove that all the elements of the matrix R are finite. Then the condition D1 of Theorem 3.2 [18] will be satisfied, which implies the asymptotic normality of the conditional least squares estimators θˆCLS. Let us denote the joint moments of the random variables Xt−i, Xt−j, Xt−k and Xt−l, i<j<k<l, by μi,j,k,l(r1,r2,r3,r4)=EXt−ir1Xt−jr2Xt−kr3Xt−lr4, r1≥0, r2≥0, r3≥0, r4≥0. Then the elements of the matrix R are the joint moments μi,j,k,l(r1,r2,r3,r4), where 1≤r1+r2+r3+r4≤4. By applying the Lemma 2 [30] when ri≥1 with si=∑j=1irj, j∈{2,3,4}, we obtain that

μi,j(r1,r2)≤μi(s2)r1s2μj(s2)r1s2<∞,

μi,j,k(r1,r2,r3)≤μi(s3)r1s3μj(s3)r2s3μk(s3)r3s3<∞,

μi,j,k,l(r1,r2,r3,r4)≤(μi(s4))r1/s4(μj(s4))r2/s4(μk(s4)r3/s4(μl(s4))r4/s4<∞,

since the rth order moments of the random variable Xt are finite. Thus, we have proven that all elements of the matrix R are finite.□

Appendix A.6 A.6 Proof of Theorem 6

To prove this theorem, we will follow the proof of Theorem 5.1 in Tjøstheim [18], i.e. we will show that the conditions E1–E3 of the mentioned theorem are satisfied. As it was mentioned earlier in the proof of Theorem 4, the BINARCH(p) model {Xt}t∈Z is strictly stationary and ergodic process with finite moments of any order. Also, functions gt(θ)=E(Xt|Ft−1)=na0+∑i=1paiXt−i and Var(Xt|ft−1)=nαt(1−αt) are almost surely three times continuously differentiable with respect to the vector θ in an open set which contains true vector θ0. Having in mind that ∂αt∂a0=1,∂αt∂ai=1nXt−i,i=1,2…p, and the results given in Theorem 5 (recall Lemma 2 in Zhu and Wang [30]), it is easily seen that

E|∂φt∂ai|≤E(Xt2+nαt+2αtXt2+3nαt2+n2αt2+2nαt2Xt+2nαt3nαt2(1−αt)2)=E(E(Xt2+nαt+2αtXt2+3nαt2+n2αt2+2nαt2Xt+2nαt3nαt2(1−αt)2|Ft−1))=E((1+2αt)[nαt(1−αt)+n2αt2]+nαt+3nαt2+n2αt2+2n2αt3+2nαt3nαt2(1−αt)2)=2E((1+2αt)(1+nαt)αt(1−αt)2)<∞,

since 0<a0≤αt≤∑i=0pai<1. Similarly, based on the fact that ∂2αt∂ai∂aj=0,i,j∈{0,1,…,p}, we find

E∂2ϕt∂ai∂aj<∞,i,j∈{0,1,…,p}.

We conclude that condition E1 is satisfied. In a similar way, it can be shown that condition E3 is also satisfied. The only thing left to prove is that condition E2 is fulfilled, i.e. if we assume that for arbitrary real numbers d0, d1, …, dp, it holds that

(15)E((d0n+∑i=1pdiXt−i)2nαt(1−αt))+E((1−2αt)(d0n+∑i=1pdiXt−i)nαt(1−αt))2=0,

then d0=d1=⋯=dp=0. Equality (15) implies

d0n+∑i=1pdiXt−i2αt(1−αt)(n+4)+1nαt(1−αt)2=a.s.0.

It is easy to conclude that αt(1−αt)(n+4)+1≠0, and further implies that nd0+∑i=1pdiXt−i=a.s.0. The rest of the proof is equivalent to the proof of Theorem 4.□

Appendix A.7 A.7 Proof of Theorem 7

Following the proof of Theorem 5 and using the fact that 0<a0≤αt≤∑i=0pai<1, we obtain that all the elements of the matrix S are finite. This implies that condition F1 of Theorem 5.2 in Tjøstheim [18] is satisfied, which proves the asymptotic normality of the estimators {θˆMLTP} obtained by minimizing (7).□

Appendix A.8 A.8 Proof of Theorem 8

Let μt=E(Xt). Using the fact that E(αt)=μt/n, we obtain from the eq. (8) that

(16)μt=na0+∑i=1paiμt−i+∑j=1qbjμt−j.

The eq. (16) represents a non-homogeneous difference equation. According to Goldberg [31], this equation has a finite stable solution which is independent of t if all the roots of the equation

1−∑i=1paiz−i−∑j=1qbjz−j=0

lie inside the unit circle. Since the parameters ai, i=0,1,…,p, and bj, j=1,2,…,q, are non-negative, it follows that the roots lie inside the unit circle if ∑i=1pai+∑j=1qbj<1. Finally, this condition is satisfied since a0+∑i=1pai+∑j=1qbj<1.□

Appendix A.9 A.9 Proof of Theorem 9

The proof of Theorem 9 is an immediate consequence of the following Lemma 1, which provides the dependence between the random variables Xt and αt−k for k∈Z.

Lemma 1

Let {Xt}t∈Z be a BINGARCH(p,q) model given by (1) and (8). Then the covariance function between the random variables Xt and αt−k, k∈Z, is given by

Cov(Xt,αt−k)={nCov(αt,αt−k),k≥0,1nCov(Xt,Xt−k),k<0.

Proof. Let us first derive the covariance function of the random variables Xt and Xt−k, k>0. Since {Xt}t∈Z is a first-order stationary process, we have that μ=E(Xt), for all t∈Z. Since E(Xt|ft−1)=nαt and Xt−k is a Ft−1-measurable function, we obtain that

Cov(Xt,Xt−k)=E(Xt−μ)(Xt−k−μ)=E(Xt−k−μ)E(Xt|Ft−1)−μ=nCov(αt,Xt−k),

which implies that Cov(Xt,αt−k)=1nCov(Xt,Xt−k), for k<0.

Let us now derive the covariance Cov(Xt,αt−k) for k≥0. Let At represents the σ-field generated by the random variables {αt}t≥0. Since E(Xt|At)=E(E(Xt|At,Ft−1)|At)=nαt, the random variable αt−k is an At-measurable function and E(αt)=μ/n, we obtain that

Cov(Xt,nαt−k)=E(Xt−μ)(nαt−k−μ)=E(nαt−k−μ)E(E(Xt|At,Ft−1)|At)−μ

=Cov(nαt,nαt−k),

which implies that Cov(Xt,αt−k)=nCov(αt,αt−k) for k≥0.□

Appendix B Appendix B CLS estimators of the BINARCH(1) model

Proposition 1 below provides closed-form expressions for the matrices U and R for BINARCH(1) model. Before, we provide the following lemma with the expressions for the first four moments, which will be used to derive the elements of the matrices U and R.

Lemma 2

Let {Xt}t∈Z be the BINARCH(1) process given by (1) and (2). Then μ1≡E(Xt)=nA1, A1≡E(αt)=a0/(1−a1), and the second, third and fourth moments μi and Ai, i=2,3,4, of the random variables Xt and αt, respectively, are given as solutions of the following equations

μ2A2=1−n(n−1)a12n2−1−1nA1−a02−2a0a1A1

μ3A3=1−n(n−1)(n−2)a13n3−1−1nA1+3n(n−1)A2−a03−3a02a1A1−3a0a12n2μ2

μ4A4=1−24C4na14n4−1−1nA1+7n(n−1)A2+6n(n−1)(n−2)A3−a04−4a03a1A1−6a02a12n2μ2−4a0a13n3μ3,

where C4n=n(n−1)(n−2)(n−3)/24.

The proof of Lemma 2 follows after standard calculations, which are based on the fact that the random variable Xt for given Xt−1 has the binomial distribution with parameters n and αt, and expressions for its moments about 0.

Proposition 1

Let {Xt}t∈Z be the BINARCH(1) process given by (1) and (2). Then the matrix U is given by

U=n2nμ1nμ1μ2,

where μ1 and μ2 are given as in the previous Lemma. The matrix R follows as

R=nn2(A1−A2)na0(1−a0)μ1+a1(1−2a0)μ2−a12nμ3na0(1−a0)μ1+a1(1−2a0)μ2−a12nμ3a0(1−a0)μ2+a1n(1−2a0)μ3−a12n2μ4.

To derive the matrix U, we need the first- and second-order moments of the random variable Xt−1. On the other hand, for derivation of the matrix R, we need up to the fourth-order moment of the random variable Xt−1, since

E(αt(1−αt)Xt−1i)=a0(1−a0)E(Xt−1i)+a1n(1−2a0)E(Xt−1i+1)−a12n2E(Xt−1i+2),i=0,1,2.

All these moments can be obtained by using the results of the previous Lemma 2.

References

1. Weiß CH. A new class of autoregressive models for time series of binomial counts. Commun Stat Theor Method 2009;38:447–60.10.1080/03610920802233937Suche in Google Scholar

2. Weiß CH, Kim HY. Binomial AR (1) processes: moments, cumulants, and estimation. Statistics 2013;47:494–510.10.1080/02331888.2011.605893Suche in Google Scholar

3. Weiß CH, Kim HY. Parameter estimation for binomial AR (1) models with applications in finance and industry. Stat Papers 2013;54:563–90.10.1007/s00362-012-0449-ySuche in Google Scholar

4. Weiß CH, Pollett PK. Chain binomial models and binomial autoregressive processes. Biometrics 2012;68:815–24.10.1111/j.1541-0420.2011.01716.xSuche in Google Scholar PubMed

5. Steutel FW, van Harn K. Discrete analogues of self-decomposability and stability. Ann Prob 1979;7:893–9.10.1214/aop/1176994950Suche in Google Scholar

6. McKenzie E. Some simple models for discrete variate time series. Water Resour Bull 1985;21:645–50.10.1111/j.1752-1688.1985.tb05379.xSuche in Google Scholar

7. Cui Y, Lund R. Inference in binomial AR (1) models. Stat Prob Lett 2010;80:1985–90.10.1016/j.spl.2010.09.003Suche in Google Scholar

8. Scotto MG, Weiß CH, Silva ME, Pereira I. Bivariate binomial autoregressive models. J Multivar Anal 2014;125:233–51.10.1016/j.jmva.2013.12.014Suche in Google Scholar

9. Al-Osh MA, Alzaid AA. Binomial autoregressive moving average models. Stochastic Models 1991;7:261–82.10.1080/15326349108807188Suche in Google Scholar

10. Ferland R, Latour A, Oraichi D. Integer-valued GARCH process. J Time Ser Anal 2006;27:923–42.10.1111/j.1467-9892.2006.00496.xSuche in Google Scholar

11. Weiß CH, Pollett PK. Binomial autoregressive processes with density-dependent thinning. J Time Ser Anal 2014;35:115–32.10.1002/jtsa.12054Suche in Google Scholar

12. Zhu F. A negative binomial integer-valued GARCH model. J Time Ser Anal 2011;32:54–67.10.1111/j.1467-9892.2010.00684.xSuche in Google Scholar

13. Zhu F. Modeling overdispersed or underdispersed count data with generalized Poisson integer-valued GARCH models. J Math Anal Appl 2012;389:58–71.10.1016/j.jmaa.2011.11.042Suche in Google Scholar

14. Xu H-Y, Xie M, Goh TN, Fu X. A model for integer-valued time series with conditional overdispersion. Comput Stat Data Anal 2012;56:4229–42.10.1016/j.csda.2012.04.011Suche in Google Scholar

15. Gonçalves E, Mendes-Lopes N, Silva F. Infinitely divisible distributions in integer-valued GARCH models. J Time Ser Anal 2015;36:503–27.10.1111/jtsa.12112Suche in Google Scholar

16. Held L, Hofmann M, Höhle M, Schmid V. A two-component model for counts of infectious diseases. Biostatistics 2006;7:422–37.10.1093/biostatistics/kxj016Suche in Google Scholar

17. Höhle M, Paul M. Count data regression charts for the monitoring of surveillance time series. Comput Stat Data Anal 2008;52:4357–68.10.1016/j.csda.2008.02.015Suche in Google Scholar

18. Tjøstheim D. Estimation in nonlinear time series models. Stochastic Processes Appl 1986;21:251–73.10.1016/0304-4149(86)90099-2Suche in Google Scholar

19. Billingsley P. Statistical inference for Markov processes. Chicago: University of Chicago Press, 1961.10.2307/1401956Suche in Google Scholar

20. Robert-Koch-Institut. SurvStat@RKI http://www3.rki.de/SurvStat, data status: June 12, 2013.Suche in Google Scholar

21. Faulde M, Sobe D, Kimmig P, Scharninghausen J. Renal failure and hantavirus infection in Europe. Nephrol Dialysis Transpl 2000;15:751–3.10.1093/ndt/15.6.751Suche in Google Scholar PubMed

22. Heyman P, Vaheri A, Lundkvist Å, Avsic-Zupanc T. Hantavirus infections in Europe: from virus carriers to a major public-health problem. Expert Rev Anti-Infect Therapy 2009;7:205–17.10.1586/14787210.7.2.205Suche in Google Scholar PubMed

23. MacNeil A, Comer JA, Ksiazeka TG, Rollin PE. Sin Nombre virus-specific immunoglobulin M and G kinetics in hantavirus pulmonary syndrome and the role played by serologic responses in predicting disease outcome. J Infect Dis 2010;202:242–6.10.1086/653482Suche in Google Scholar PubMed

24. Schilling S, Emmerich P, Klempa B, Auste B, Schnaith E, Schmitz H, et al. Hantavirus disease outbreak in Germany: limitations of routine serological diagnostics and clustering of virus sequences of human and rodent origin. J Clin Microbiol 2007;45:3008–14.10.1128/JCM.02573-06Suche in Google Scholar PubMed PubMed Central

25. Tsay RS. Model checking via parametric bootstraps in time series analysis. Appl Stat 1992;41:1–15.10.2307/2347612Suche in Google Scholar

26. Jung RC, Tremayne AR. Convolution-closed models for count time series with applications. J Time Ser Anal 2011;32:268–80.10.1111/j.1467-9892.2010.00697.xSuche in Google Scholar

27. Weiß CH. A Poisson INAR (1) model with serially dependent innovations. Metrika 2015;78:829–51.10.1007/s00184-015-0529-9Suche in Google Scholar

28. Davis RA, Liu H. Theory and inference for a class of observation-driven models with application to time series of counts. eprint arXiv:1204.3915. 2012.Suche in Google Scholar

29. Seneta E. Non-negative matrices and Markov chains, 2nd ed. New York: Springer Verlag, 1983.Suche in Google Scholar

30. Zhu F, Wang D. Estimation and testing for a Poisson autoregressive model. Metrika 2011;73:211–30.10.1007/s00184-009-0274-zSuche in Google Scholar

31. Goldberg S. Difference equations. New York: Wiley, 1958.Suche in Google Scholar

Published Online: 2015-12-5

Published in Print: 2016-11-1

Artikel in diesem Heft

https://doi.org/10.1515/ijb-2015-0051

Schlagwörter für diesen Artikel

binomial INARCH (p) model; binomial INGARCH (p,q) model; overdispersion; parameter estimation; epidemic surveillance