AA new INARMA(1, 1) model with Poisson marginals
Johannes BracherOctober 17, 2019
Epidemiology, Biostatistics and Prevention Institute, University of Zurich,Hirschengraben 84, 8001 Zurich, Switzerland [email protected]
Bibliographic note:
This is a pre-print (submitted version before peer review) of a contributionin Steland, A., Rafajlowicz, E., Okhrin, O. (Eds.):
Stochastic Models, Statistics and Their Ap-plications , p. 323–333, published by Springer Nature Switzerland, 2019. The final authenticatedversion is available at https://doi.org/10.1007/978-3-030-28665-1_24 . Abstract
We suggest an INARMA(1, 1) model with Poisson marginals which extends the INAR(1) in asimilar way as the INGARCH(1, 1) does for the INARCH(1) model. The new model is equivalentto a binomially thinned INAR(1) process. This allows us to obtain some of its stochastic propertiesand use inference methods for hidden Markov models. The model is compared to various othermodels in two case studies.
Time series of counts are encountered in a broad variety of contexts. Two popular modellingapproaches are the INAR (integer-valued autoregressive [1]) and INGARCH (integer-valued gener-alized autoregressive conditional heteroscedasticity [5]) classes. In this article an extension of thePoisson INAR(1) model is proposed which parallels the generalization of the INARCH(1) to the IN-GARCH(1, 1) model. We give some properties of the new model, which we refer to as INARMA(1,1), and point out how to do inference via methods for hidden Markov processes. The performanceof the model is compared to various INAR and INGARCH-type models in two case studies.
The Poisson INAR(1) model { X t , t ∈ Z } with parameters ν > < α < X t = I t + α ◦ X t − . (1)The operator ◦ denotes binomial thinning, i.e. α ◦ Y = (cid:80) Yi =1 Z i with Z i iid ∼ Bernoulli( α ), implying α ◦ Y | Y ∼ Bin(
Y, α ). The sequence { I t , t ∈ Z } consists of independent Poisson random variables withrate ν . All thinning operations are performed independently of each other and of { I t } . Moreover,at each time t , α ◦ X t and I t are independent of all X u , u < t . Marginally the X t are then Poissondistributed with rate ν/ (1 − α ); the autocorrelation function is ρ ( h ) = α h . a r X i v : . [ s t a t . M E ] O c t he Poisson INARCH(1) model { X t , t ∈ Z } is usually defined as [10] X t | X t − , X t − , . . . ∼ Pois( λ t ); λ t = ν + αX t − (2)with ν > , α ≥
0, but can also be formulated as (compare [11]) X t = I t + α ∗ X t − (3)where { I t } is again a sequence of independent Poisson random variables with rate ν . We definethe operator ∗ as α ∗ Y | Y ∼ Pois( αY ) where for αY = 0 we include the degenerate Poissondistribution Pr( α ∗ Y = 0) = 1. Note that while X t − in (3) is integer-valued, α ∗ Y is also definedfor real-valued Y ≥
0. If α <
1, the process { X t } is stationary with E ( X t ) = ν/ (1 − α ) andVar( X t ) = E ( X t ) / (1 − α ), i.e. X t is over-dispersed for α >
0. The autocorrelation function is again ρ ( h ) = α h . The INARCH(1) model can be extended to the Poisson INGARCH(1, 1) model [5] X t | X t − , X t − , . . . ∼ Pois( λ t ) (4) λ t = ν + αX t − + βλ t − (5)with ν > α, β ≥
0. In the following we assume { X t } to be stationary, which is the case if α + β <
1. We can then express it using the operator ∗ from (3). Consider S t = λ t − ν − β − β which after some simple algebra leads to λ t = (1 − β ) S t + ν − βS t = βS t − + α − β · X t − . Note that the recursive definition (5) of λ t implies λ t ≥ (cid:80) ∞ d =0 νβ d = ν/ (1 − β ) so that non-negativityof S t is ensured. An alternative display of (4)–(5) is then X t = φ ∗ S t + I t (6) S t = (1 − φ ) S t − + κX t − (7)with I t iid ∼ Pois( τ ) and τ = ν − β ; φ = 1 − β ; κ = α − β . Stationarity of { X t } implies 0 ≤ κ < E ( X t ) = τ − κ ; Var( X t ) = 1 − ξ + κ φ − ξ · E ( X t ) ρ ( h ) = 1 − ξ + κ φ + κφ (1 − φ )1 − ξ + κ φ · κφξ h − with ξ = 1 − φ (1 − κ ), i.e. the second-order properties of an ARMA(1, 1) process. e now suggest a similar generalization of the Poisson INAR(1) model which we call PoissonINARMA(1, 1). It is defined as { X t , t ∈ Z } with X t = φ ◦ S t + I t (8) S t = S t − − ( X t − − I t − ) + κ ◦ X t − (9)where I t iid ∼ Pois( τ ) , τ > < φ ≤ , < κ <
1. Again, all thinning operations are independentof each other and of { I t } . At each t , φ ◦ S t , κ ◦ X t and I t are independent of all X u , S u , u < t and, given X t , κ ◦ X t is independent of S t . This formulation parallels (6)–(7) as, using X t − − I t − = φ ◦ S t − ,it is easily seen that S t d = (1 − φ ) ◦ S t − + κ ◦ X t − . However, (9) implies a dependence between the two thinnings φ ◦ S t and (1 − φ ) ◦ S t , entering into X t and S t +1 , respectively, as they are forced to sum up to S t . Unlike in the INGARCH model (6)–(7), S t is discrete-valued here (it can be shown to be an INAR(1) process with S t = J t + ξ ◦ S t − ; J t ∼ Pois( κτ )). This is necessary to ensure well-definedness of φ ◦ S t and achieved by replacing themultiplications from (7) by binomial thinnings.As in an INAR(1) model, the X t are marginally Poisson distributed under model (8)–(9), therate being E ( X t ) = Var( X t ) = τ / (1 − κ ) . (10)The autocorrelation function is ρ ( h ) = φκξ h − (11)where again ξ = 1 − φ (1 − κ ). These properties are easy to show using the representation of { X t } asa binomially thinned INAR(1) process, see next section. Thus the new model, too, has the second-order properties of an ARMA(1, 1) process, justifying the name INARMA(1, 1). Note, however,that the formulation differs from other models referred to as INARMA in the literature (e.g. [7]).The INAR(1) model corresponds to the boundary case φ = 1 of the new class. In comparison tothe INGARCH(1, 1) model with the same parameters the new model has lower dispersion and itsautocorrelation function is damped if φ < The INARMA(1, 1) model can be interpreted as follows: X t is the number of fertile females ina population and S t is the (unobserved) number of juvenile, i.e. not yet fertile females. I t is thenumber of fertile female immigrants (there is no immigration of juveniles). Females do not die beforereaching fertility and at each time of their juvenile period have a probability of φ to transition tothe fertile state. They stay fertile for exactly one time period and can have at most one femaleoffspring, the probability of which is κ . A graphical display of such a system can be found in Figure1.The time from a female’s birth to its fertile period obviously follows a geometric distribution withparameter φ . We can use this to express the model as an INAR( ∞ ) model X t = ∞ (cid:88) i =1 α i ◦ X t − i + I t ; I t iid ∼ Pois( τ ) (12) The INGARCH(1, 1) model, too, can be expressed with a discrete-valued process { S t } , just set S t = S t − − ( X t − − I t − ) + κ ∗ X t − in (8)–(9). Details are omitted due to space constraints. As mentioned in [5], Lemma 2, the fact that the autocovariance structure of { X t } coincides with that of a stationaryARMA(1, 1) process is sufficient for { X t } to be an ARMA(1, 1) process itself. t X t I t S t + X t + I t + . . .. . . enter fertile periodwith probability φ immigrationat rate τ remain juvenile withprobability 1 − φ femaleoffspring withprobability κ Figure 1:
Interpretation of the INARMA(1, 1) process in the form of a flow diagram.with α i = κφ (1 − φ ) i − , i = 1 , , . . . and dependent thinning operations given by B t | X t ∼ Bin( X t , κ ) (13) A ( j ) t iid ∼ Geom( φ ) , j = 1 , . . . B t (14) α i ◦ X t = B t (cid:88) j =1 I ( A ( j ) t = i ) , i = 1 , , . . . (15)Here, B t is the number of female offspring born in t and A ( j ) t is the waiting time until fertility forthe j th of the females born at time t . The definition (13)–(15) of the dependent thinnings impliesthat α i ◦ X t | X t ∼ Bin( X t , α i ) for i = 1 , , . . . under the constraint (cid:80) ∞ i =1 α i ◦ X t ≤ X t .The representation (12)–(15) nicely illustrates the relationship of the INARMA(1, 1) model toother common models. Replacing the geometric waiting time distribution in (14) by a one-pointdistribution with Pr( A ( j ) t = 1) = 1 yields the INAR(1) model while a categorical distribution withsupport 1 , . . . , p gives the INAR( p ) model by Alzaid and Al-Osh [2] (see [3] for details). Replacingthe binomial offspring distribution in (13) by a Poisson distribution, i.e. setting B t | X t ∼ Pois( κX t )yields the INGARCH(1, 1) model. Due to space restrictions, we do not detail on this, but it isstraightforward to show using the INARCH( ∞ ) representation of the INGARCH(1, 1) model ([12],p.76) and some basic properties of the Poisson distribution. The INARCH(1) and INARCH( p )models can be obtained by using again one-point and categorical waiting time distributions in (14).We recently encountered INAR( ∞ ) models of type (12)–(15) in [3] where we extended workby Fern´andez-Fontelo et al [6] on underreported INAR models. We showed that { X t , t ∈ Z } isequivalent to a binomially thinned INAR(1) model { ˜ Y t , t ∈ Z } given by Y t = J t + ξ ◦ Y t − (16)˜ Y t | Y t ∼ Bin( Y t , φκ/ξ ) (17)with J t iid ∼ Pois( τ ξ/κ ) and, as before, ξ = 1 − φ (1 − κ ). This represents an interesting parallel tothe Gaussian ARMA(1, 1) model which, as shown for instance in [9], can be obtained by addinghomoscedastic measurement error to a Gaussian AR(1) process. This third representation of theprocess makes the derivation of equations (10)–(11) easy (see [6], Section 2 and Appendix A; ourmodel corresponds to the special case ω = 1 of the class discussed there). Also, it implies that manyproperties of the INAR(1) process translate to the INARMA(1, 1) model, e.g. the marginal Poissondistribution and time-reversibility [8]. Inference for higher-order INAR models with dependent thinning operations is challenging as thelikelihood is generally intractable [2]. For our model, however, we can exploit the representation ω = 1 of the class treated in [6]). As the state space of the latent process { Y t } is infinite,truncation at some reasonably large value Y max is required. The maximum of the log-likelihood isthen obtained by numerical optimization. We apply the four models from Sections 2 and 3 to two data sets. The first example consists ofthe gold particle counts from Westgren [13], a data set which is often used in the literature. Forinstance Weiß ([12], p.48) applies Poisson INAR(1), INAR(2) and CINAR(2) models to these data.To make our results comparable to these analyses we fit all models to observations 501–870 of Series(C). As a second example we use weekly counts of mumps in Bavaria, Germany, from week 1/2014to week 52/2017 (downloaded from on 8 Oct 2018). Mumps, a viral disease,used to be a common childhood disease. Since the introduction of a vaccine in the 1950s it hasbecome rare in Germany, but remains under routine surveillance. The data are displayed in Figure6. Both time series show slowly decaying autocorrelation functions, indicating that relaxing theAR(1) assumption may be beneficial. While the particle counts are approximately equidispersed(mean 1.55, variance 1.65) the mumps data show some overdispersion (mean 2.49, variance 3.93).
Westgren's gold particle data p a r ti c l e c oun t s A C F mean = 1.55 var = 1.65 Mumps in Bavaria (weekly) ca s e c oun t s A C F mean = 2.49 var = 3.93 Figure 2:
Case studies: gold particle counts and weekly numbers of reported mumps cases in Bavaria.Table 1 shows parameter estimates and AIC values for the gold particle data. For comparison weadded the results Weiß [12] obtains for INAR(2) and CINAR(2) models (see there for details). Toassess the out-of-sample predictive performance we computed mean log scores (logS, [4]) of one-step-ahead forecasts for the second half of the time series. For each of these forecasts the models werere-fitted to the data which were already available at the respective time point (“rolling forecast”).Note that the log score is negatively oriented, i.e. smaller values are better.The INARMA(1, 1) model has the best in-sample and out-of-sample performance. Interestingly italso outperforms the two AR(2) models from [12], indicating that observations more than two timesteps back still contain additional information.The corresponding results for the mumps data can be found in Table 2. While the INARMA(1,1) model again represents a considerable improvement compared to the INAR(1), the INGARCH(1,1) model performs best. This is not surprising given the overdispersion in the data.
Parameter estimates, AIC values and mean log scores for the gold particle data. Mean log scores forINAR(2) and CINAR(2) were computed using code from the supplement of [12].
Model Parameter AIC logSˆ ν ˆ α ˆ α ˆ λ Poisson INAR(1) 0.73 0.53 1040 1.642Poisson INAR(2) 0.54 0.47 0.18 1027 1.610Poisson CINAR(2) 0.60 0.41 0.19 1027 1.611Poisson INARCH(1) 0.75 0.52 0.00 1057 1.624ˆ τ ˆ φ ˆ κ ˆ S AICPoisson INARMA(1, 1) 0.31 0.67 0.80 1014 1.577Poisson INGARCH(1, 1) 0.47 0.54 0.70 1.85 1047 1.592
Table 2:
Parameter estimates, AIC values and mean log scores for the mumps data.
Model Parameter AIC logSˆ ν ˆ α ˆ λ Poisson INAR(1) 2.21 0.11 842 2.017Poisson INARCH(1) 2.08 0.15 9.01 832 2.010ˆ τ ˆ φ ˆ κ ˆ S AICPoisson INARMA(1, 1) 1.21 0.25 0.52 827 1.963Poisson INGARCH(1, 1) 1.12 0.26 0.52 18.97 816 1.955
We suggested an INARMA(1, 1) model with Poisson marginals which draws conceptually fromthe INGARCH(1, 1) model [5] and the INAR( p ) model by Al-Osh and Alzaid [2]. We provided analternative representation in terms of an offspring distribution and a waiting time distribution whichenlightens the close relation to several existing models from the literature. A third representationas a binomially thinned INAR(1) process turned out to be useful for inference and to obtain somestochastic properties. In our case studies the model performed favourably for equidispersed data.For overdispersed data it outperformed the INAR(1) model, but achieved lower performance thanthe INGARCH(1, 1) model. This raises the question how overdispersion could be accommodatedin the model. Alternative immigration distributions could be considered, but the model would thenno longer have a representation as a binomially thinned version of a Markov process. Obtaining itsstochastic properties and evaluating the likelihood would get more involved. Another open questionis how higher-order INARMA( p, q ) models could be defined. Data and code
Data and R code are available at github.com/jbracher/inarma. Acknowledgements
The author thanks Christian H. Weiß for helpful feedback.
References [1] M.A. Al-Osh and A.A. Alzaid. First-order integer-valued autoregressive (INAR(1)) process.
JTime Ser Anal , 8(3):261–275, 1987.
2] A.A. Alzaid and M.A. Al-Osh. An integer-valued pth-order autoregressive structure (INAR(p))process.
J Appl Probab , 27(2):314–324, 1990.[3] J. Bracher. Comment on “Under-reported data analysis with INAR-hidden Markov chains”.
Stat Med , 38(5):893–898, 2019.[4] C. Czado, T. Gneiting, and L. Held. Predictive model assessment for count data.
Biometrics ,65(4):1254–1261, 2009.[5] R. Ferland, R. Latour, and D. Oraichi. Integer-valued GARCH process.
J Time Ser Anal ,27(6):923–942, 2006.[6] A. Fern´andez-Fontelo, A. Caba˜na, P. Puig, and D. Mori˜na. Under-reported data analysis withINAR-hidden Markov chains.
Stat Med , 35(26):4875–4890, 2016.[7] P. Neal and T. Subba Rao. MCMC for integer-valued ARMA processes.
J Time Ser Anal ,28(1):92–110, 2006.[8] S. Schweer. On the time-reversibility of integer-valued autoregressive processes of general order.In A. Steland, E. Rafaj(cid:32)lowicz, and K. Szajowski, editors,
Stochastic Models, Statistics and TheirApplications , pages 169–177, Cham, 2015. Springer.[9] J. Staudenmayer and J.P. Buonaccorsi. Measurement error in linear autoregressive models.
JAm Stat Assoc , 100(471):841–852, 2005.[10] C.H. Weiß. The INARCH(1) model for overdispersed time series of counts.
Commun StatA-Theor , 39(6):1269–1291, 2010.[11] C.H. Weiß. A Poisson INAR(1) model with serially dependent innovations.
Metrika , 78(7):829–851, 2015.[12] C.H. Weiß.
An Introduction to Discrete-Valued Time Series . Wiley, Hoboken, NJ, 2018.[13] A. Westgren. Die Ver¨anderungsgeschwindigkeit der lokalen Teilchenkonzentration in kolloidenSystemen.
Arkiv f¨or Matematik, Astronomi och Fysik , 11(14):1–24, 1916.[14] W. Zucchini and I. MacDonald.
Hidden Markov Models for Time Series . Chapman andHall/CRC, New York, NY, 2009.. Chapman andHall/CRC, New York, NY, 2009.