Integer-valued AR processes with Hermite innovations and time-varying parameters: An application to bovine fallen stock surveillance at a local scale

Abstract

In this article we present a new INteger-valued AutoRegressive (INAR) model with the aim of extracting baseline patterns of cattle fallen stock registered over an 5-year period at a local scale. We introduce HINAR as a generalization of the classical Poisson-based INAR models whose innovations follow a Hermite distribution. In order to assess trends and seasonality in these time series, we fit different models with time-dependent parameters by specifying proper functions. Using real world examples, we illustrate how to estimate parameters by maximum likelihood and validate the fitted models. We also show a detailed method to forecast. Our proposed model supposes a good solution for studying discrete time series when the counts have many zeros, low counts and moderate overdispersion. This model has been applied to the analysis of fallen cattle registered at a local scale as part of the development of a veterinary syndromic surveillance system.

Keywords

cattle fallen stock Hermite distribution INAR models Time dependent coefficients veterinary syndromic surveillance

1 Introduction

The well-known autoregressive moving average (ARMA) models are suitable for continuous or discrete time series which can be approximated conveniently by Gaussian distributions. However, in many situations in which discrete time series have low counts and/or many zeros, ARMA models can be inappropriate. Many models for discrete time series have been proposed in the literature. McKenzie (2003) summarizes some of them. One of the best-known proposals is the INteger-valued AutoRegressive models (abbreviated as INAR).

A classical INAR(p) process ${X_{t}; t \in ℤ}$ is defined by means of the recursion

X_{t} = ρ_{1} \circ X_{t - 1} + ρ_{2} \circ X_{t - 2} + \dots + ρ_{p} \circ X_{t - p} + Z_{t},

(1.1)

where

0 \leq ρ_{1}, ρ_{2}, \dots, ρ_{p} \leq 1

are fixed parameters and

Z_{t}

is a sequence of independent and identically distributed non-negative integer-valued random variables with mean ensuremath E(Z_t) = μ_z and finite variance Var(Z_t) = σ²_z, being also independent of

X_{t - j}

for all

j

In this process, the integer discreteness of the series is ensured by the binomial thinning operator or binomial subsampling $\circ$ , defined as

ρ_{j} \circ X_{t - j} = \sum_{i = 1}^{X_{t - j}} B_{i} (ρ_{j}),

(1.2)

where

{B_{i} (ρ_{j})}

are independent and identically distributed Bernoulli random variables with

P (B_{i} (ρ_{j}) = 1) = ρ_{j}

, which are also independent of

X_{t - j}

for all

j

. Steutel and van Harn (1979) introduced for the first time the binomial thinning operator and McKenzie (1985) proposed to use this operator for count time series. Comprehensive reviews on INAR models can be found in McKenzie (2003), Jung and Tremayne (2006) and Scotto et al. (2015).

INAR(1) models are those that are mostly used in practice, being extended and generalized in different ways. Some authors proposed other extensions of the classical Poisson INAR(1) model considering geometric or zero-inflated Poisson innovations (see Jazi et al., 2012a; Jazi et al., 2012b). Weiß (2008) reviewed INAR(1) models based on different kinds of thinning operators. Zheng et al. (2007) constructed INAR(1) processes with random coefficient thinning operators describing its properties and methods of estimation. Another extension is the RCINAR processes by Gomes and Canto e Castro (2009).

Some studies have been conducted concerning change-point detection in binary time series (see Fokianos et al., 2014; Hudecová, 2013). Moreover, some authors have taken into account not only change-points detections in binary time series but also in INAR models. In this sense, Hudecová et al. (2015a) and Hudecová et al. (2015b) have developed on-line procedures for detecting change-points in the parameters of INAR(1) models, proposing a statistical test based on the probability generating function of the innovations.

Different INAR(p) models ( $p > 1$ ) have been considered in the literature depending on the joint distribution of $(ρ_{1} \circ X_{t - 1}, \dots, ρ_{p} \circ X_{t - p})$ . Alzaid and Al-Osh (1990) assumed a conditional multinomial distribution, while Du and Li (1991) imposed conditional independence. On the other hand, Zheng et al. (2005) extended INAR(p) models considering random thinning coefficients.

It is important to remark that in the INAR(p) model of Du and Li (1991), $ρ_{j} \circ X_{t - j} | X_{t - j} = x_{t - j}$ are binomial random variables with parameters $ρ_{j}$ and $x_{t - j}$ , all of them being mutually independent. The model of Du and Li (1991) is more tractable in practice and, henceforth, it will be considered in this work.

Parameter $ρ_{j}$ in equation (1.1) represents the proportion of counts which emigrates from state $t - j$ to state $t$ , and $Z_{t}$ can be seen as the innovations or new events produced at time $t$ . In this way, the counts at time $t$ consist of survivors from the state $t - 1$ , from the state $t - 2, \dots$ , from the state $t - p$ , plus the innovations.

In the classical INAR(p) process, equation (1.1), the parameters $ρ_{1}, ρ_{2}, \dots, ρ_{p}$ are constants, but in many real world applications these parameters could vary over time taking into account, for instance, the trend and/or seasonality of the series. The distribution of the innovations $Z_{t}$ is usually assumed to be Poisson, although other alternatives have been considered in the literature (zero-inflated Poisson, weighted Poisson, negative binomial, $\dots$ ).

For the INAR(1) model, it is interesting to remark that only when $Z_{t}$ follows a rth-Hermite distribution, which includes the Poisson ( $r = 1$ ) and the Hermite distributions ( $r = 2$ ), will the marginal distribution of the process $X_{t}$ belong to the same family as the innovations $Z_{t}$ (see Weiß and Puig, 2015).

In this work, we propose a more general model that includes the classical INAR(p) of Du and Li (1991) with Poisson innovations. We will assume that the parameters can be time-dependent and the innovations will be Hermite distributed. In the next sections, we will see that our model is a flexible extension of the classical INAR(p) model, being applicable to a broader framework including the analysis of bovine fallen stock time series.

2 The HINAR(p) model

Let $X_{1}$ and $X_{2}$ be two independent Poisson variables with parameters $a_{1}$ and $a_{2}$ . Hence, the probability distribution of the variable $Y = X_{1} + 2 X_{2}$ is the Hermite distribution with parameters $a_{1}$ and $a_{2}$ (see Kemp and Kemp, 1965), and its probability mass function is expressed as

\begin{matrix} P (Y = k) = e^{- (a_{1} + a_{2})} \sum_{j = 0}^{[\frac{k}{2}]} \frac{a_{1}^{k - 2 j} a_{2}^{j}}{(k - 2 j)! j!}, \end{matrix}

(2.1)

for

k = 0, 1, 2, \dots, a_{1}, a_{2} \geq 0

, where

[\frac{k}{2}]

is the integer part of

\frac{k}{2}

Following is a summary of some important properties of the Hermite distribution whose probability generating function (pgf) is $G_{Y} (s) = e^{a_{1} (s - 1) + a_{2} (s^{2} - 1)}$ .

The moment generating function exists and is equal to $M_{Y} (t) = e^{a_{1} (e^{t} - 1) + a_{2} (e^{2 t} - 1)}$ , and its characteristic function is $ϕ_{Y} (t) = e^{a_{1} (e^{it} - 1) + a_{2} (e^{2 it} - 1)}$ .

The mean and variance are $μ = a_{1} + 2 a_{2}$ and $σ^{2} = a_{1} + 4 a_{2}$ , respectively. The third central moment is $γ_{1} = \frac{a_{1} + 8 a_{2}}{{(a_{1} + 4 a_{2})}^{\frac{3}{2}}}$ and corresponds to the skewness. Whereas the fourth moment centred around the mean is $β_{2} = \frac{a_{1} + 16 a_{2}}{(a_{1} + 4 a_{2})^{2}} + 3$ and corresponds to the kurtosis.

Given two independent Hermite random variables such as $H_{1} \sim$ Hermite $(a_{1}, a_{2})$ and $H_{2} \sim$ Hermite $(b_{1}, b_{2})$ , the sum of both are also Hermite distributed, $M = H_{1} + H_{2} \sim$ Hermite $(a_{1} + b_{1}, a_{2} + b_{2})$ .

A consequence of property 2 is that the dispersion index $δ = σ^{2} / μ$ satisfies $1 \leq δ \leq 2$ . It can be noted that when $δ = 1$ , the Hermite distribution is just a Poisson distribution, and when $δ = 2$ , it is twice a Poisson distribution. Consequently, the Hermite distribution can be useful to model data which have a moderate overdispersion. In order to facilitate the computation of the Hermite probabilities and its quantiles, the hermite package in R could be a quality tool (see Moriña et al., 2015a, 2015b).

The third property implies that the Hermite distribution is closed under addition or convolution. This will be relevant for forecasting in Section 3. In fact, Puig and Valero (2007) proved that the rth-Hermite distributions are the unique count data distributions with $r$ parameters which are closed under binomial thinning and under addition. In this regard, Kemp and Kemp (1965) described more properties in further detail.

It is said that a distribution is Poisson zero-inflated if its proportion of zeros is greater than the proportion of zeros of a Poisson distribution with the same mean. Then, a good strategy to model time series, with low counts and many zeros, is to use Poisson zero-inflated distributions for the innovations. A direct solution is to consider the so-called zero-inflated Poisson distribution (Jazi et al., 2012a), but Hermite distribution is also Poisson zero-inflated because, according to equation (2.1), $P (Y = 0) = e^{- a_{1} - a_{2}} \geq e^{- μ}$ .

The novel extension of the conventional INAR(p) model proposed in this study considers that the innovations follow a Hermite distribution with time-dependent parameters $a_{1} (t)$ and $a_{2} (t)$ . This model will be called HINAR(p) and is expressed as

\begin{matrix} X_{t} = ρ_{1} (t) \circ X_{t - 1} + ρ_{2} (t) \circ X_{t - 2} + \dots + ρ_{p} (t) \circ X_{t - p} + Z_{t} (a_{1} (t), a_{2} (t)) . \end{matrix}

(2.2)

In addition, the former parameters of the model are functions of time

ρ_{i} (t) = f_{i} (β_{i}, t)

, provided that the outputs of these functions fall in the interval [0, 1]. Each function introduces a new set of parameters (

β_{i}, i = 1, \dots, p

) that have to be estimated. The functions acting as parameters of the Hermite distribution

a_{j} = f_{j} (γ_{j}, t), j = 1, 2

have to take positive values. In this sense, different choices of

ρ_{i}

and

a_{j}

have been considered for the purpose of exploring different patterns of trend and seasonality. For our examples, we have it found convenient to use a second-order trigonometric polynomial with a linear part, jointly with an appropriate link function. Specifically, for the Hermite parameters

a_{j}

, the choice has been,

\begin{matrix} a_{j} (t) & = & exp (γ_{j 0} + γ_{j 1} t + γ_{j 2} sin (\frac{2 π t}{T_{1}}) + γ_{j 3} cos (\frac{2 π t}{T_{1}}) + γ_{j 4} sin (\frac{2 π t}{T_{2}}) \\ + γ_{j 5} cos (\frac{2 π t}{T_{2}})) . \end{matrix}

(2.3)

Notice that $γ_{j 1}$ captures the possible trend of the series, while the other parameters capture the seasonalities of periods $T_{1}$ and $T_{2}$ respectively. Similarly, for $ρ_{1} (t), ρ_{2} (t), \dots, ρ_{p} (t)$ , we have taken a logit link function leading to the expression,

\begin{matrix} log (\frac{ρ_{j} (t)}{1 - ρ_{j} (t)}) & = β_{j 0} + β_{j 1} t + β_{j 2} sin (\frac{2 π t}{T_{1}}) + β_{j 3} cos (\frac{2 π t}{T_{1}}) + β_{j 4} sin (\frac{2 π t}{T_{2}}) \\ + β_{j 5} cos (\frac{2 π t}{T_{2}}) . \end{matrix}

(2.4)

In our examples, the time is measured in weeks. These are just epidemiological weeks, also called Epi or CDC weeks. Therefore, taking $T_{1} = 52$ and $T_{2} = 26$ , expressions in equations (2.3) and (2.4) would consider annual and biannual patterns of seasonality. It is quite common to find these patterns of seasonality in mortalities, in both beef and dairy cattle farms.

When coefficients $ρ_{i}$ and $a_{j}$ are constant, some results concerning stationarity will be provided in the next section.

2.1 HINAR(p) model with constant coefficients

Let $X_{t}$ be an INAR(p) process as in equation (1.1) with $ρ_{p} > 0$ . Du and Li (1991) showed that to ensure the stationarity of the process, all the roots of the characteristic polynomial $y^{p} - ρ_{1} y^{p - 1} - \dots - ρ_{p} = 0$ should be inside the unit circle which is equivalent to the more simple condition $\sum_{i = 1}^{p} ρ_{i} < 1$ . Du and Li (1991) also showed that the correlation structure of this process is the same as the classical AR(p) model. Accordingly, under the stationarity assumption, the marginal distribution represented by the random variable $X$ satisfies,

\begin{matrix} X \sim ρ_{1} \circ X + ρ_{2} \circ X + \dots + ρ_{p} \circ X + Z_{t} . \end{matrix}

(2.5)

Considering

s^{X}

, taking expectations and using the property of independence under conditioning of the INAR(p) model, we obtain the following expression for the pgf

ϕ (s)

of the marginal distribution in terms of the pgf of the innovations

ϕ_{Z} (s)

\begin{matrix} ϕ (s) = ϕ (P (s)) ϕ_{Z} (s), \end{matrix}

(2.6)

where

P (s) = \prod_{i = 1}^{p} ((1 - ρ_{i}) + ρ_{i} s)

. Notice that for

p = 1

, the expression in equation (2.6) is just the expression for discrete self-decomposability, and

P (s)

is just the pgf of a sum of

p

independent Bernoulli random variables with different probabilities

ρ_{i}

. The expectation E(X) = μ and variance Var(X) = σ² of the marginal distribution can also be derived from equation (2.6). Taking the first and second derivatives with respect to

s

in equation (2.6), evaluating them at

s = 1

and using the fact that

ϕ^{'} (1) = μ

ϕ^{″} (1) = σ^{2} + μ^{2} - μ

and

P^{'} (1) = \sum_{i = 1}^{p} ρ_{i}

P^{″} (1) = (\sum_{i = 1}^{p} ρ_{i})^{2} - \sum_{i = 1}^{p} ρ_{i}^{2}

, we find,

\begin{matrix} μ = & \frac{μ_{z}}{1 - \sum_{i = 1}^{p} ρ_{i}} \\ σ^{2} = & \frac{μ_{z} (\sum_{i = 1}^{p} ρ_{i} - \sum_{i = 1}^{p} ρ_{i}^{2}) + σ_{z}^{2} (1 - \sum_{i = 1}^{p} ρ_{i})}{(1 - (\sum_{i = 1}^{p} ρ_{i})^{2}) (1 - \sum_{i = 1}^{p} ρ_{i})} . \end{matrix}

(2.7)

δ_{z} = σ_{z}^{2} / μ_{z}

indicates the dispersion index of the innovations, the dispersion index of the marginal

δ = σ^{2} / μ

remains,

\begin{matrix} δ = \frac{1 - \sum_{i = 1}^{p} ρ_{i}^{2}}{1 - (\sum_{i = 1}^{p} ρ_{i})^{2}} + \frac{δ_{z} - 1}{1 + \sum_{i = 1}^{p} ρ_{i}} . \end{matrix}

(2.8)

An immediate consequence is that $δ \geq 1 + \frac{δ_{z} - 1}{1 + \sum_{i = 1}^{p} ρ_{i}}$ . Consequently, for Poisson innovations and $p > 1$ , the marginal distribution is not Poisson distributed because it is overdispersed. Hermite innovations $Z_{t} (a_{1}, a_{2})$ have mean and variance equal to $μ_{z} = a_{1} + 2 a_{2}$ and $σ_{z}^{2} = a_{1} + 4 a_{2}$ , respectively. Because $1 \leq δ_{z} \leq 2$ , Hermite innovations produce marginal overdispersed distributions with dispersion index in the interval,

\begin{matrix} δ \in [\frac{1 - \sum_{i = 1}^{p} ρ_{i}^{2}}{1 - (\sum_{i = 1}^{p} ρ_{i})^{2}}, \frac{1 - \sum_{i = 1}^{p} ρ_{i}^{2}}{1 - (\sum_{i = 1}^{p} ρ_{i})^{2}} + \frac{1}{1 + \sum_{i = 1}^{p} ρ_{i}}] . \end{matrix}

(2.9)

It is interesting to point out that overdispersed innovations can produce marginals having the same dispersion index. For instance, an HINAR(2) model with

ρ_{1} = ρ_{2} = 1 / 4

and Hermite innovations

Z_{t} (2, 1)

has a marginal distribution with dispersion index

δ = δ_{z} = 3 / 2

As we have commented before, for $p = 1$ the marginal distribution is Poisson (Hermite), if and only if the innovations $Z_{t}$ are Poisson (Hermite) distributed. However, for $p > 1$ , it is impossible to obtain a marginal Poisson (Hermite) distributed. To show this, we first need the following lemma:

Lemma 1. (Milne and Wescott (1993), p. 374). Let $log (f (s))$ be a polynomial of degree m, that is,

f (s) = exp (a_{0} + a_{1} s + a_{2} s^{2} + \dots + a_{m - 1} s^{m - 1} + a_{m} s^{m}),

where at least

a_{1}

a_{m - 1}

and

a_{m}

are non-zero. Then,

a_{1} > 0

a_{m - 1} > 0

and

a_{m} > 0

are necessary conditions for

f (s)

to be a pgf.

Proposition Consider an INAR(p) process, $p > 1$ , like in (1.1). It does not exist any sequence of independent and identically distributed count innovations $Z_{t}$ giving a stationary marginal Poisson (Hermite) distribution.

Proof. Suppose that the marginal distribution is Poisson (Hermite). Then $log (ϕ (s))$ in equation (2.6) is a polynomial of degree 1 (2) and $log (ϕ (P (s)))$ is a polynomial of degree p (2p). Then, $log (ϕ_{Z} (s))$ is a polynomial of degree p (2p) with a negative coefficient of the leading term. Consequently, Lemma 1 shows that $ϕ_{Z} (s)$ is not a pgf arriving to a contradiction.

The statement of proposition can also be extended to any distribution having a polynomic log-pgf (r-Hermite distributions).

In general, for $p > 1$ the pgf of the marginal distribution can not be obtained in a closed form, but a MA-type representation can be derived from equation (2.6). When evaluating $ϕ (s)$ at $P (s)$ , $ϕ (P (s)) = ϕ (P (P (s))) ϕ_{Z} (P (s))$ and $ϕ (s) = ϕ (P (P (s))) ϕ_{Z} (s) ϕ_{Z} (P (s))$ . By iterating this scheme we obtain, $ϕ (s) = ϕ (P^{(n)} (s)) ϕ_{Z} (s) ϕ_{Z} (P (s)) ϕ_{Z} (P^{(2)} (s)) \dots ϕ_{Z} (P^{(n - 1)} (s))$ , where $P^{(k + 1)} (s) = P (P^{(k)} (s))$ , $k = 1, 2, . . .$ . As $n$ tends to infinity, $P^{(n)} (s)$ tends to its fixed point $s = 1$ , being an attractor, because $P^{'} (1) = \sum_{i = 1}^{p} ρ_{i} < 1$ . Therefore, we obtain the representation

log (ϕ (s)) = log (ϕ_{Z} (s)) + \sum_{i = 1}^{\infty} log (ϕ_{Z} (P^{(i)} (s))) .

(2.10)

For Hermite innovations

Z (a_{1}, a_{2})

log (ϕ_{Z} (s)) = a_{1} (s - 1) + a_{2} (s^{2} - 1)

and the log-pgf of the marginal distribution can be numerically computed from equation (2.10). In this way, the probability of zero can be directly computed. Numerical derivatives at

s = 0

are required to calculate the other probabilities.

2.2 Parameter estimation

Parameters $Θ = (β_{1}, β_{2}, \dots, β_{l}, γ_{1}, γ_{2}, \dots, γ_{m})$ are estimated by means of the maximum likelihood method. According to the Total Probability Theorem, the likelihood function can be expressed as the product of probabilities conditioned to a set of past observations. Moreover, because of the existing dependence of order $p$ , the likelihood function can be expressed as

\begin{matrix} L (X, Θ) & = P (x_{1}, Θ) P (x_{2} | x_{1}, Θ) \dots P (x_{k} | x_{k - 1}, x_{k - 2}, \dots, x_{k - p + 1}, Θ) \\ \dots P (x_{n} | x_{n - 1}, x_{n - 2}, \dots, x_{n - p + 1}, Θ) . \end{matrix}

(2.11)

Considering that the following probabilities

P (x_{1}, Θ)

P (x_{2} | x_{1}, Θ)

P (x_{3} | x_{2}, x_{1}, Θ), \dots,

P (x_{p} | x_{p - 1}, \dots, x_{1} Θ)

are equal to

1

, the likelihood function conditioned to the

p

initial observations remains,

L (x_{t}; Θ) = \prod_{t = p + 1}^{n} P (x_{t} | x_{t - 1}, \dots, x_{t - p}; Θ) .

(2.12)

In order to evaluate this likelihood function we need the following proposition.

Proposition Given a HINAR(p) process like in (2.2), then

\begin{matrix} P (x_{t} | x_{t - 1}, \dots, x_{t - p}; Θ) = & e^{- (a_{1} (t) + a_{2} (t))} \sum_{i_{1} = 0}^{m_{1}} x_{t - 1} choose i_{1} ρ_{1}^{i_{1}} (t) {(1 - ρ_{1} (t))}^{x_{t - 1} - i_{1}} \\ \cdot & \sum_{i_{2} = 0}^{m_{2}} x_{t - 2} choose i_{2} ρ_{2}^{i_{2}} (t) {(1 - ρ_{2} (t))}^{x_{t - 2} - i_{2}} . \dots \\ \cdot & \sum_{i_{p} = 0}^{m_{p}} x_{t - p} choose i_{p} ρ_{p}^{i_{p}} (t) {(1 - ρ_{p} (t))}^{x_{t - p} - i_{p}} \sum_{j = 0}^{[\frac{M}{2}]} \frac{a_{1} (t)^{M - 2 j} a_{2} (t)^{j}}{(M - 2 j)! j!}, \end{matrix}

(2.13)

where

M = x_{t} - i_{1} - i_{2} - \dots - i_{p}

m_{1} = min (x_{t}, x_{t - 1})

m_{2} = min (x_{t - 2}, x_{t} - i_{1})

\dots

m_{p} = min (x_{t - p}, x_{t} - (i_{1} + \dots + i_{p - 1}))

Proof Let X and Y be two independent count random variables, where $Z$ is the sum of X and Y. Then, according to the Total Probability Theorem, it can be showed that

\begin{matrix} P (Z = z) & = P (X + Y = z) = P (Y = z - X) = \sum_{j = 0}^{z} P (Y = z - X | X = j) P (X = j) \\ = \sum_{j = 0}^{z} P (Y = z - j) P (X = j) . \end{matrix}

On the other hand,

X_{t}

conditioned to the last

p

observation, is the sum of

p

independent binomial distributions plus a Hermite distribution. That is,

\begin{matrix} X_{t} = x_{t} | X_{t - 1} = x_{t - 1}, X_{t - 2} = x_{t - 2}, \dots, X_{t - k} = x_{t - k} \\ \sim B (x_{t - 1}, ρ_{1} (t)) + B (x_{t - 2}, ρ_{2} (t)) + \dots + B (x_{t - k}, ρ_{k} (t)) + Z_{t} \end{matrix}

(2.14)

By induction, using the expression of the binomial probability function and the expression in equation (2.1), the proof holds.

Due to the complexity of the likelihood function of the model in equation (2.12), there are no closed expressions for the maximum likelihood estimation (MLE) of the parameters. Nevertheless, an appropriate numerical method can be used to maximize the likelihood function and obtain the MLE. In this article, a Newton-type algorithm provided by the function nlm in em R was used.

The code of the R-program built to estimate the parameters is included in the supplementary material. Notice that the number of parameters could be very large, so it is important to find parsimonious models. In this work, the model selection is based on both statistical significance of the parameters and Akaike information criterion (AIC) (or Bayesian information criterion, BIC). That is, the best model among all possible models is the one in which all parameters are statistically significant and which also shows a small value of AIC (or BIC).

3 Forecasting

Two kinds of predictions will be considered: the average behaviour of the series and the k-time-ahead distribution that will be studied for $p = 1$ , due to the nature of our examples. McCabe et al. (2011) presented a non-parametric approach to estimate the forecasting distribution of $X_{t + k}$ for $k \geq 1$ , when the process $X_{t}$ follows a classical INAR(p) model with arbitrary innovations. They estimate the forecasting distributions using non-parametric likelihood maximum estimators (NPMLE) which are basically obtained by maximizing the likelihood function of an INAR(p) model, but considering that the innovations are distributed following some discrete parameter-free distribution called $G$ , which is also estimated. As a result, the NPMLE estimated vector $\hat{θ} = ({\hat{ρ}}_{1}, \dots, {\hat{ρ}}_{p}, {{\hat{g}}_{r}})$ is obtained, where sequence ${{\hat{g}}_{r}}$ corresponds to the estimate of the parameter-free discrete distribution $G$ . It is worth mentioning that this method does not allow to consider seasonality and/or time-dependent parameters, while the method proposed in this work allows to handle both seasonality and time-dependent parameters. However, for comparative purposes, the method of McCabe et al. (2011) will be used in the examples of Section 4.

3.1 The average behaviour of the time series

In order to compute point predictions for the average behaviour of the series, we apply an extension of the methodology described by Moriña et al. (2011). Accordingly, taking expectations in equation (2.2) and calling E(X_t)= _t, we obtain the following recurrence relation,

μ_{t} = ρ_{1} (t) μ_{t - 1} + ρ_{2} (t) μ_{t - 2} + \dots + ρ_{p} (t) μ_{t - p} + a_{1} (t) + 2 a_{2} (t) .

(3.1)

For the first

p

initial values of

μ_{t}

, we can consider arbitrary values. In our case, we have used the overall mean of the series.

The variance of $μ_{t}$ is calculated using the delta method. First, for each $t$ , the expression of $μ_{t}$ in equation (3.1) can be understood as a function of the parameters $Θ$ ,

μ_{t} = μ (β_{1}, β_{2}, \dots, β_{l}, γ_{1}, γ_{2}, \dots, γ_{m}) .

(3.2)

Next, the gradient of

μ_{t}

is computed in order to apply the delta method technique. That is,

\nabla μ_{t} = (\frac{\partial μ_{t}}{\partial β_{1}}, \frac{\partial μ_{t}}{\partial β_{2}}, \dots, \frac{\partial μ_{t}}{\partial γ_{1}}, \frac{\partial μ_{t}}{\partial γ_{2}}, \dots, \frac{\partial μ_{t}}{\partial γ_{m}}), t = 1, \dots, n .

(3.3)

Finally, the variance of

μ_{t}

is computed by using the following expression which depends on the variance–covariance matrix of the MLE of

Θ

(Σ)

and the above estimated gradient,

̂ Var (\hat{μ_{t}}) = {\hat{\nabla μ_{t}}}^{t} Σ \hat{\nabla μ_{t}}, t = 1, \dots, n .

(3.4)

Furthermore, as the closed form of the gradient, equation (3.3), is intractable, it is numerically calculated by a first-order approximation of the derivatives, taking successive values of

h

until the difference between two consecutive

h

becomes insignificant

\nabla {\hat{μ}}_{t} = (\frac{μ_{t} (β_{1} + h, \dots, β_{l}, γ_{1}, \dots, γ_{m}) - {\hat{μ}}_{t}}{h}, \dots, \frac{μ_{t} (β_{1}, \dots, β_{l}, γ_{1}, \dots, γ_{m} + h) - {\hat{μ}}_{t}}{h}) .

(3.5)

Computationally, we have built a R-program (see supplementary material) in which the gradient is calculated numerically using the numDeriv package (see Gilbert and Varadham, 2012). This method improves this first-order approximation by applying Richardson's extrapolation.

As usual, the $(1 - α) (100) %$ confidence intervals of these point predictions are obtained by means of this commonly expression

μ_{t} = {\hat{μ}}_{t} \pm z_{1 - \frac{α}{2}} \sqrt{} Var (\hat{μ_{t}}) .

(3.6)

3.2 The

k

time-ahead forecasting distribution for

p = 1

Given $n$ observations $x_{1}, \dots, x_{n}$ , to make predictions with the HINAR( $1$ ) model, it is required to know the distribution of $X_{n + k}$ . To do this, we first need the following lemma.

Lemma Let $Z (a_{1}, a_{2})$ Hermite distributed according to equation (2.1), and $ρ \in (0, 1)$ . Then,

ρ \circ Z (a_{1}, a_{2}) \sim Z (ρ a_{1} + 2 a_{2} ρ (1 - ρ), ρ^{2} a_{2}) .

Proof The pgf of a $ρ$ -thinning of a count random variable $X$ has the form $ϕ_{ρ \circ X} (s) = ϕ_{X} (1 - ρ + ρ s)$ (see Puig and Valero, 2007). Because a $ρ$ -thinning of a Hermite distribution is also Hermite, and the log-pgf is a second order polynomial, equating its coefficients holds the proof.

In order to simplify the notation used in equation (2.2), we will denote $ρ_{1} (t) = ρ_{t}$ , $a_{1} (t) = a_{t}$ and $a_{2} (t) = b_{t}$ . Assuming first that $k = 1$ , notice that $X_{n + 1}$ conditioned to $X_{n} = x_{n}$ is the sum of an independent binomial distribution with parameters $(x_{n}, ρ_{n + 1})$ , plus a Hermite distribution with parameters $a_{n + 1}$ and $b_{n + 1}$ . Consider now the case $k = 2$ . Then, the distribution of $X_{n + 2}$ can be determined through the last observed value of the series in the following way,

\begin{matrix} X_{n + 2} & = ρ_{n + 2} \circ X_{n + 1} + Z (a_{n + 2}, b_{n + 2}) = \\ = ρ_{n + 2} \circ (ρ_{n + 1} \circ X_{n} + Z (a_{n + 1}, b_{n + 1})) + Z (a_{n + 2}, b_{n + 2}) = \\ = ρ_{n + 2} ρ_{n + 1} \circ X_{n} + Z (c_{n + 2}, d_{n + 2}), \end{matrix}

where

c_{n + 2} = ρ_{n + 2} (a_{n + 1} + 2 (1 - ρ_{n + 2}) b_{n + 1}) + a_{n + 2}

and

d_{n + 2} = ρ_{n + 2}^{2} b_{n + 1} + b_{n + 2}

Therefore, the distribution of $X_{n + 2}$ conditioned to $X_{n} = x_{n}$ is the sum of independent binomial random variables with parameters $(x_{n}, ρ_{n + 2} ρ_{n + 1})$ , plus a Hermite distribution with parameters $c_{n + 2}$ and $d_{n + 2}$ . This recursive procedure can be extended to $X_{n + k}, k > 2$ , leading to the following proposition:

Proposition 3.

\begin{matrix} X_{n + k} = f_{n + k} \circ X_{n} + Z (c_{n + k}, d_{n + k}), k = 1, 2, . . . \end{matrix}

(3.7)

where

f_{n + k} = \prod_{i = 1}^{k} ρ_{n + i}

(3.8)

and

\begin{matrix} c_{n + k} & = ρ_{n + k} (c_{n + k - 1} + 2 d_{n + k - 1} (1 - ρ_{n + k})) + a_{n + k} \end{matrix}

(3.9)

d_{n + k} = ρ_{n + k}^{2} d_{n + k - 1} + b_{n + k}, k = 2, 3, . . .,

(3.10)

where

c_{n + 1} = a_{n + 1}

and

d_{n + 1} = b_{n + 1}

. Note that the property of closure under convolution of the Hermite distribution mentioned in Section 2 is essential to obtain this result. Consequently, the distribution of

X_{n + k}

is the sum of a binomial with parameters

(x_{n}, f_{n + k})

and a Hermite(

c_{n + k}, d_{n + k}

) both independent. When

f_{n + k}

tends to zero as

k

increases, the long-term forecasting distribution of

X_{n + k}

(

k

large) follows an approximate Hermite distribution with parameters

c_{n + k}

and

d_{n + k}

. The case in which the parameters are time constant (

ρ_{n} = ρ

a_{n} = a

and

b_{n} = b

) has special practical interest. In this situation,

f_{n + k} = ρ^{k}

tends to zero as

k

tends to infinity and

X_{n + k} \sim

Hermite

(a^{*}, b^{*})

, where the mean and variance are those of the marginal distribution given in equation (2.7), that is,

a^{*} + 2 b^{*} = (a + 2 b) / (1 - ρ)

and

a^{*} + 4 b^{*} = (a (ρ + 1) + 2 b (2 + ρ)) / (1 - ρ^{2})

By replacing parameters in equations (3.9), (3.10) and (3.10) by their MLE, we can estimate regions of prediction of size $1 - α$ , finding the lower and upper limits $r_{1}$ and $r_{2}$ that satisfy

\begin{array}{l} \sum_{i = 0}^{r 1} P (X_{n + k} = i) \approx \frac{a}{2} \\ \sum_{i = 0}^{r 2} P (X_{n + k} = i) \approx 1 - a . \end{array}

(3.11)

Moriña et al. (2011) applied a similar methodology to predict the number of beds needed by month at the hospital level in the emergency service units. They used an INAR(2) model incorporating seasonality only in the innovations.

4 Examples of application

Our model has been used to extract the baseline patterns of mortality recorded at a farm level in two small populations of beef and dairy cattle. The continuous collection and analysis of automated data of fallen stock of cattle has shown a good potential as a health indicator (see Alba et al., 2015). If abnormal values of fallen cattle are out of the expected ranges according to the baseline patterns previously studied in a specific population, these peaks could indicate unusual events related to animal health issues. These signs could also pre-date the diagnosis and predict possible disease outbreaks. Several studies have been focused on the analysis of fallen stock in large populations that resulted in observations over time. However, many series should be studied at a local scale because important events may occur at this geographical level. That is, a region in Spain is organized at different geographical levels; province, county and municipality. Data recorded at a region or province is usually suitable to be studied by using ARMA models, while data recorded at a county and municipality has sometimes low values and/or many zeros. Also, in order to build a robust syndromic surveillance system based on cattle mortality, it would be quite necessary to explore and analyze fallen stock patterns at these different geographical levels. These patterns of fallen stock can be monitored by observing the number of carcasses collected at a farm level over time, although in small animal populations the series has many zeros and/or low counts. Therefore, our HINAR model could be a suitable tool to solve this problem.

4.1 Example 1. Modelling of fallen stock from a small population of beef cattle

The data of this study are provided by the Ministry of Agriculture, Food and Environment of Spain under strict agreement of confidentiality. These data record the fallen stock collected between January 2007 and December 2011 from 16 beef cattle farms with a total census of 1140 bovines placed in a small region of Catalonia. This dataset registers the number of bovine carcasses collected by the disposal services, the date of collection and the identification of each farm. For the purpose of studying the baseline patterns of fallen bovine of this region, the number of fallen cattle was aggregated by a week.

The range of carcasses collected by week is between 0 and 9 (see Figure 1). From 2007 to 2011 the median number of fallen cattle is 1, with a weekly average of 1.5. The standard deviation of the series is 1.7, a slight overdispersion is detected.

Figure 1:

Evolution of fallen stock in a small beef cattle population from 2007 to 2011

It seems to be a change-point at $m = 104$ (the beginning of 2009), according to the behaviour of the series shown in Figure 1. Change points in the parameters of a HINAR model can be tested using a likelihood ratio test (LRT). Under the null hypothesis, parameters $ρ$ , $a_{1}$ and $a_{2}$ are equal throughout the periods $1$ to $m$ and $m + 1$ to $T$ , that is

H_{0} : ρ_{{m \leq 104}} = ρ_{{m > 104}}; {a_{i}}_{{m \leq 104}} = {a_{i}}_{{m > 104}}, i = 1, 2 .

The estimates of the parameters under the null hypothesis are

\hat{ρ} = 0.110

{\hat{a}}_{1} = 0.572

and

{\hat{a}}_{2} = 0.334

, while the estimates without constraints are

{\hat{ρ}}_{{m \leq 104}} = 0.000

{\hat{a_{1}}}_{{m \leq 104}} = 0.817

{\hat{a_{2}}}_{{m \leq 104}} = 0.651

{\hat{ρ}}_{{m > 104}} = 0.038

{\hat{a_{1}}}_{{m > 104}} = 0.611

and

{\hat{a_{2}}}_{{m > 104}} = 0.135

. It leads to a LRT statistic equal to

38.744

corresponding to a p-value equal to

0.000

. Therefore, a model with constant parameters does not seem a good choice.

In order to fit the series properly, different hermite-based integer-valued autoregressive models (HINAR) are tested. That is, we have fitted several HINAR models with different temporal lags and patterns of trend and seasonality. The model we propose here is the best, based on our criteria for model selection (i.e., lowest AIC and BIC, and models with statistically significant coefficients). Accordingly, the estimated proposed model is a HINAR(1) with $ρ_{1} (t) = ρ_{1}$ and a decreasing trend included in the Hermite parameter $a_{1}$ ,

\begin{matrix} {\hat{X}}_{t} = ρ_{1} \circ X_{t - 1} + Z_{t} (e^{γ_{1} t}, a_{2}) . \end{matrix}

(4.1)

The MLE of the parameters of the model are presented in Table 1. The low value of the standard error of the Hermite parameter $a_{2}$ implies that $a_{2} > 0$ , and consequently the innovations of the model are not Poisson distributed. The AIC = 841.47 and BIC = 852.15 of this model are lower than the AIC = 915.13 and BIC = 922.25 of the equivalent Poisson-based integer-valued autoregressive model (PoINAR), leading to a LRT statistic equal to 75.65. Under the null hypothesis, $H_{0} : a_{2} = 0$ , the LRT does not have an asymptotic $χ_{1}^{2}$ distribution as expected, because $a_{2} = 0$ belongs to the boundary of the domain of parameters. It can be shown that the asymptotic distribution is a 50:50 mixture of the constant 0 and the $χ_{1}^{2}$ distribution, leading in this case to a p-value equal to $p = 0.0000$ .

Table 1:

Maximum likelihood estimates of the model proposed for the series of beef cattle

Parameter	ML estimate	s.e.
$ρ_{1}$	0.104	0.044
$γ_{1}$	$-$ 0.004	0.001
$a_{2}$	0.306	0.046

Figure 2:

ACF and PACF of the residuals of model, in equation (4.1)

Figure 3:

Predictions made in static cross-validation and forecasting of the HINAR(1) model for the beef cattle. The average behaviour of the time series (top) and the k time-ahead forecasting distribution (bottom). Confidence limits at 95%

The auto-correlation function (ACF) and partial auto-correlation function (PACF) plots (Figure 2) allow us to validate the residuals of the model. It can be noted that their behaviour is close to white noise. The predictive performance of the proposed model is also assessed using two techniques of cross-validation. On the one hand, we use a static cross-validation in order to explore the predictive capability of the model. Using a static form, the dataset is divided into training and testing sets: the first is used to estimate the model; whereas the second is used for the model validation. In our case, the training set includes the first 80% of the observations registered over the period of study (up to 2011), while the testing sample includes the last 20% of the observations (2011). The model is fitted with the training set to predict the expected values over the period corresponding to the testing set. Then, the training set predictions are compared with the testing set observations. The average prediction of the series and the $k$ time-ahead forecasting distribution based on the estimated parameters according to equation (3.1) is shown in Figure 3. The estimates of these parameters based on the training sample (up to 2011) are similar ( $\hat{ρ_{1}} = 0.108$ , $\hat{γ_{1}} = - 0.006$ and $\hat{a_{2}} = 0.328$ ) to those obtained based on the whole sample (see Table 1). In addition, the 95% confidence limits calculated from the $k$ time-ahead forecasting distribution show that the number of carcasses observed after 2011 agree with the proposed model. The approach proposed by McCabe et al. (2011) to estimate the $k$ time-ahead forecasting distributions is also performed and the results are shown in Figure 9. In this case, the NPMLE of the parameters (and their standard errors) are ${\hat{ρ}}_{1}^{NP} = 0.137$ ( $0.048$ ), with an estimated distribution of the innovations equal to ${\hat{g}}_{0} = 0.445$ ( $0.044$ ), ${\hat{g}}_{1} = 0.266$ ( $0.036$ ), ${\hat{g}}_{2} = 0.132$ ( $0.029$ ), ${\hat{g}}_{3} = 0.075$ ( $0.019$ ), ${\hat{g}}_{4} = 0.028$ ( $0.013$ ), ${\hat{g}}_{5} = 0.026$ ( $0.012$ ), ${\hat{g}}_{6} = 0.008$ ( $0.008$ ), ${\hat{g}}_{7} = 0.010$ ( $0.008$ ), ${\hat{p}}_{8} = 0.003$ ( $0.005$ ) and ${\hat{g}}_{9} = 0.007$ ( $0.006$ ). Figure 9 shows that the non-parametric-based 95% confidence limits (approximately $[0, 6]$ ) are larger than those computed using our method (approximately $[0, 4]$ ), summarized in Figure 3.

Finally, because this HINAR(1) model has a constant $ρ_{1}$ , the approximated long-term forecasting distribution can also be considered. However, this distribution is almost indistinguishable from the $k$ time-ahead forecasting distribution because ${\hat{ρ}}_{1}$ has a very low value. On the other hand, we use a dynamic cross-validation which considers that the training and testing samples are changing over time. That is, an initial sample is used for estimating the model and making predictions for a near future time period. Once the predicted values are available, these are compared with the corresponding observed values which will be incorporated into the series over time. In that case, it is required to identify whether some abnormal peak is caused by a disease or a catastrophe, or whether it is only a regular random value of the series. If the peak is really anomalous, it will be replaced by its predicted value in order to keep the capability of the fitted model to detect abnormalities of this magnitude. This procedure is repeated in such a way that all the last available information is re-used to estimate the model, and it is also used to re-make predictions for a period without information. When such information is available, then predictions are re-validated. Hence, dynamic cross-validation can be performed at the same time the information is updated.

The dynamic cross-validation of the model in equation (4.1) is shown in Figure 4. This validation is built in the following way. First, an initial sample based on the information of 2007 is considered. Second, this initial sample is used to predict a week ahead. Next, this prediction is compared with the corresponding observed value; if this value is anomalous, it is studied to determine whether it is actually abnormal or not. Finally, if it is abnormal, it will be replaced by its predicted median value. This process is repeated as new information is incorporated in the sample.

Figure 4:

Dynamic cross-validation of the HINAR(1) model for the beef cattle

Notice that the proposed model forecasts one week ahead. This elapsed time is considered appropriate to provide continuous information and gives responses avoiding the distortion that exists at a daily level due to the lack of carcasses collected over the weekend.

In this example, the values registered in weeks 58, 66, 75 and 81 could be anomalous according to the forecasting estimates. In our example, based on the available information, they have been considered as regular random values.

4.2 Example 2. Modelling of fallen stock from a small population of dairy cattle

In this example, the fallen stock data are also collected between January 2007 and December 2011 from three dairy cattle farms with a total census of 647 bovines placed in the same region as the previous example. A range between 0 and 7 is collected (see Figure 5), and from 2007 to 2011 the median number of collection recorded is 0.5, with a weekly average of 1 carcass. A slight overdispersion in the series with a dispersion index of 1.33 is also detected.

Figure 5:

Evolution of fallen stock in a small dairy cattle population from 2007 to 2011

Table 2:

Maximum likelihood estimates of the model proposed for the series of dairy cattle

Parameter	ML estimate	s.e.
$β_{0}$	$-$ 3.894	1.404
$β_{2}$	$-$ 3.304	1.594
$a_{1}$	0.388	0.065
$γ_{0}$	$-$ 1.540	0.197
$γ_{3}$	0.698	0.240

As in the previous example, different HINAR models are tested assessing trends and seasonality based on equations (2.3) and (2.4). Taking into account our criteria for model selection, the best model among all possible models is a HINAR(1) leading to the process,

\begin{matrix} {\hat{X}}_{t} = (\frac{e^{β_{0} + β_{2} sin (\frac{2 π t}{52})}}{1 + e^{β_{0} + β_{2} sin (\frac{2 π t}{52})}}) \circ X_{t - 1} + Z_{t} (a_{1}, e^{γ_{0} + γ_{3} cos (\frac{2 π t}{52})}), \end{matrix}

(4.2)

where annual seasonality has been included in both parameters

ρ_{1}

and

a_{2}

. A summary of the MLE of the parameters is shown in Table 2. The AIC = 694.96 and BIC = 712.76 of this model are also lower than the AIC = 745.25 and BIC = 752.37 of the equivalent PoINAR. In this case, the asymptotic distribution of the LRT statistic is unknown since the null hypothesis takes the form of

H_{0} : e^{γ_{0} + γ_{3} \cos (\frac{2 π t}{52})} = 0

. However, the p-value of the LRT could be calculated using a parametric bootstrap method.

The plots of the ACF and PACF are shown in Figure 6. Notice that their behaviour is also close to white noise.

Figure 6:

ACF and PACF of the residuals of model (4.2)

Figure 7:

Predictions made in static cross-validation and forecasting of the HINAR(1) model for the dairy cattle. The average behaviour of the time series (top) and the k time-ahead forecasting distribution (bottom). Confidence limits at 95%

Finally, in order to validate our model not only the checking of its residuals is taken into account but also its capability as a predictive tool. As in the example 1, the static and dynamic cross-validation have been conducted. Figure 7 shows the average prediction of the series as well as the $k$ time-ahead forecasting distribution based on the estimated parameters according to equation (3.1). The estimates of the parameters based on the training sample (up to 2011) have also been similar ( $\hat{β_{0}} = - 4.045$ , $\hat{β_{2}} = - 3.311$ , $\hat{a_{1}} = 0.439$ , $\hat{γ_{0}} = - 1.433$ and $\hat{γ_{3}} = 0.694$ ) to those shown in Table 2. The 95% confidence limits calculated from the $k$ time-ahead forecasting distribution also agree with the proposed model. As in the first example, the approach of McCabe et al. (2011) is applied and the corresponding 95% confidence limits are shown in Figure 9. Here, the non-parametric estimates of the parameters are ${\hat{ρ}}_{1}^{NP} = 0.093$ ( $0.049$ ), and the distribution of the innovations are ${\hat{g}}_{0} = 0.555$ ( $0.040$ ), ${\hat{g}}_{1} = 0.190$ ( $0.032$ ), ${\hat{g}}_{2} = 0.142$ ( $0.024$ ), ${\hat{g}}_{3} = 0.082$ ( $0.019$ ), ${\hat{g}}_{4} = 0.082$ ( $0.011$ ) and ${\hat{g}}_{\geq 5} = 0.008$ ( $0.006$ ). Probabilities $g_{6}$ and $g_{7}$ are aggregated to $g_{5}$ ( $g_{\geq 5}$ ) to avoid convergence problems. Comparing with the results obtained using our method (Figure 7), it is clear that the non-parametric-based approach of McCabe et al. (2011) completely ignores the seasonality pattern of the series, losing an important information.

Figure 8:

Dynamic cross-validation of the HINAR(1) model for fallen dairy cattle

Figure 8 shows the dynamic cross-validation made for the model in equation (4.2). This validation has been conducted in the same way as the validation of example 1. The observed values agree with the model except for the value recorded in the week 153. However, further study at a field level is required to determine the causes associated with this event.

5 Discussion

Over the last few years the classical Poisson INAR model has been studied in detail and applied in different situations. The HINAR model is a reasonable generalization of the classical model allowing to cope with many time series in small scale studies having low values with many zeros. The HINAR model presented in this work allow to consider time-dependent coefficients, which are useful to include the trend and/or seasonality of the series. These reasons make the HINAR model more flexible than the classical Poisson INAR model. The main advantage of the INAR(1) model with Hermite innovations with respect to other innovations proposed in the literature is that the marginal distribution of the process (with constant coefficients) is also Hermite distributed. For time-dependent models and also for $p > 1$ , the property of the Hermite distribution of being closed under independent addition simplifies notoriously the k time-ahead forecasting distribution.

Our model has been motivated by the study of mortalities in beef and dairy cattle farms at a small level. Specifically, these mortalities recorded at farm level have the potential of being used as an early warning system for assessing the impact and magnitude of abnormal events occurred over time. Moreover, these events can also have different impact at various geographical and temporal scales. For instance, a disease outbreak that causes mortality can occur in a wide geographical area or can be localized in a small area. In the first case, the analysis at a large scale is able to detect the outbreak, although in the second case it is required to analyze the data at a small scale in order to detect the anomaly. Accordingly, ARMA models are suitable for describing patterns of series at a large scale, but they are not as appropriate as the HINAR models for describing the profiles of the series at a small scale.

Figure 9:

The 95% confidence limits of the forecasting distributions from the 1-week ahead to 48-week ahead based on the approach proposed by McCabe et al. (2011)

The series of beef cattle has been fitted using a HINAR(1) model with a decreasing trend included in the innovations, while the series of dairy cattle has been fitted using a HINAR(1) with an annual seasonality included in the $ρ$ -parameter and also in the innovations. These models have been selected based on different kind of validations. In both examples, the model validation was conducted analyzing the crude residuals ( $r_{t} = x_{t} - {\hat{x}}_{t}$ ) to identify if they look like white noise (see Figures 2 and 5), and then perform a cross-validation (see Figures 3 and 6). The static cross-validation is useful to validate the baseline patterns of the model, while the dynamic cross-validation is a better choice, if the aim of the study is to forecast future values. The new models appear to be more appropriate than the classical Poisson INAR models because the AIC and BIC are lower. We have explored several patterns of the trend of seasonality using linear and trigonometric functions with appropriate links. The use of non-parametric approaches like P-splines (Eilers et al., 2015) is a topic of further research.

Why are two very different models needed to model beef versus dairy cattle farms? The answer comes from the fact that beef cattle farms are characteristic of being extensive farms. These kind of farms are mainly kept in fields and may be housed for a part of the year. Furthermore, data are derived from a few number of farms which are quite heterogeneous and some of them were closed during the period of study, agreeing with the decreasing trend detected by the model. This trend can be explained as a consequence of the impact of the major economic crisis in Spain after 2008, where many farms decided to close leading to a decreasing number of beef cattle in the regions and a decreasing number of carcasses recorded. On the other hand, the profile of the dairy cattle farm is more sophisticated, appearing to be an annual seasonality. In dairy cattle, the farms are intensive in such a way that cattle may be housed throughout their lives. Furthermore, dairy cattle are quite sensitive to temperature changes. This fact agrees with the seasonality introduced in equation (4.2) and described in the Figure 7.

To sum up, the model proposed in this article can be an appropriate and less restrictive method for a fitting discrete time series with low counts and many zeros. Specifically, it allows us to fit series assuming non-Poisson innovations and time-dependent parameters.

Footnotes

Acknowledgments

This work was funded by grants MTM2012-31118 and MTM2015-69493-R from the Ministry of Economy and Competitiveness of Spain. We thank to the Ministry of Agriculture, Food and Environment of Spain for providing the data analyzed in the article. We also give special thanks to David Moriña who helped us with classical INAR models and their implementation in the R-software and Sandy Beitsch who carefully reviewed the article for English language. The authors thank the two referees and the associate editor for their valuable comments and suggestions.

References

Alba

Dórea

Arinero

Sánchez

Cordón

Puig

Crawford

(2015) Exploring the surveillance potential of mortality data: Nine years of bovine fallen stock data collected in Catalonia (Spain). PLoS ONE 10(4), e0122547, doi:10.1371/journal.pone.0122547.

Alzaid

Al-Osh

(1990) An integer-valued pth-order autoregressive structure (INAR(p)) process. Journal of Applied Probability 27, 314–24.

(1991) The integer-valued autoregressive (INAR(p)) model. Journal of Time Series Analysis , 12, 129–42.

Eilers

PHC

Marx

Durban

(2015) Twenty years of P-splines. Statistics and Operations Research Transactions , 39, 149–86.

Fokianos

Gombay

Hussein

(2014) Retrospective change detection for binary time series models. Journal of Statistical Planning and Inference , 145, 102–12.

Gilbert

Varadham

(2012) numDeriv: Accurate numerical derivatives (R package) . URL http://CRAN.R-project.org/package=numDeriv (last accessed 28 June 2016).

Gomes

Canto e Castro

(2009) Generalized integer-valued random coefficient for a first order structure autoregressive (RCINAR) process. Journal of Statistical Planning and Inference , 139, 4088–97.

Hudecová

(2013) Structural changes in autoregressive models for binary time series. Journal of Statistical Planning and Inference , 143, 1744–52.

Hudecová

Meintanis

(2015a) Tests for time series of counts based on the probability generating function. Statistics , 49, 316–37.

10.

Hudecová

Meintanis

(2015b) Detection of changes in INAR models. In Steland

., Stochastic models, statistics and their applications, Springer Proceedings in Mathematics and Statistics. Vol. 122, pages 11–18.

11.

Jazi

Jones

Lai

(2012a) First-order integer valued AR processes with zero inflated poisson innovations. Journal of Time Series Analysis , 33, 954–63.

12.

Jazi

Jones

Lai

(2012b) Integer valued AR(1) with geometric innovations. Journal of the Iranian Statistical Society , 11, 173–90.

13.

Jung

Tremayne

(2006) Binomial thinning models for integer time series. Statistical Modelling , 6, 81–96.

14.

Kemp

(1965) Some properties of the Hermite distribution. Biometrika Trust , 52, 381–94.

15.

McCabe

BPM

Martin

Harris

(2011) Efficient probabilistic forecasts for counts. Journal of the Royal Statistical Society: Series B , 73, 253–72.

16.

McKenzie

(1985a) Some simple models for discrete variate time series. Water Resources Bull , 21, 645–50.

17.

McKenzie

(2003) Discrete variate time series. In Handbook of Statistics . pages 573–606. Amsterdam: Elsevier.

18.

Milne

Westcott

(1993) Generalized multivariate Hermite distributions and related point processes. Annals of the Institute of Statistical Mathematics , 45, 367–81.

19.

Moriña

Puig

Rȷos

Vilella

Trilla

(2011) A statistical model for hospital admissions caused by seasonal diseases. Statistics in Medicine , 30, 3125–36.

20.

Moriña

Higueras

Puig

Oliveira

(2015a) Hermite: Generalized hermite distribution (R package) . URL http://CRAN.Rproject. org/package=hermite

21.

Moriña

Higueras

Puig

Oliveira

(2015b) Generalized hermite distribution modelling with the R package hermite. The R Journal , 7, 263–74.

22.

Puig

Valero

(2007) Characterization of count data distributions involving additivity and binomial subsampling. Bernoulli , 13, 544–55.

23.

Scotto

Weiß

Gouveia

(2015) Thinning-based models in the analysis of integer-valued time series: A review. Statistical Modelling , 15, 590–618.

24.

Steutel

van Harn

(1979) Discrete analogues of self-decomposability and stability. The Annals of Probability , 7, 893–9.

25.

Weiß

(2008) Thinning operations for modeling time series of counts---a survey. Advances in Statistical Analysis , 92, 319–41.

26.

Weiß

Puig

(2015) The marginal distribution of compound Poisson INAR(1) processes. Stochastic Models, Statistics and Their Applications , 122, 351–9.

27.

Zheng

Basawa

Datta

(2005) Inference for pth-order random coefficient integer-valued autoregressive processes. Journal of Time Series Analysis , 27, 411–40.

28.

Zheng

Basawa

Datta

(2007) First-order random coefficient integervalued autoregressive processes. Journal of Statistical Planning and Inference , 173, 212–29.