A New Bivariate INAR(1) Model with Time-Dependent Innovation Vectors

Abstract

Recently, there has been a growing interest in integer-valued time series models, especially in multivariate models. Motivated by the diversity of the infinite-patch metapopulation models, we propose an extension to the popular bivariate INAR(1) model, whose innovation vector is assumed to be time-dependent in the sense that the mean of the innovation vector is linearly increased by the previous population size. We discuss the stationarity and ergodicity of the observed process and its subprocesses. We consider the conditional maximum likelihood estimate of the parameters of interest, and establish their large-sample properties. The finite sample performance of the estimator is assessed via simulations. Applications on crime data illustrate the model.

1. Introduction

Bivariate count data occur in many contexts, often as the counts of two events, objects or individuals during a certain period of time. For example, such counts occur in epidemiology when two kinds of related diseases are examined, in criminology when two kinds of crimes are committed, in business when the volume of sales of two correlated products are observed or in manufacturing when two similar products are produced.

In real application, the observed time series data are often discrete, over-dispersed (the empirical variance is greater than the empirical mean) and have other features such as time dependence. Many univariate models have been proposed to deal with integer-valued time series data based on the univariate binomial thinning operator $“\circ ”$, which is proposed by Steutel and van Harn [1]:

(1)$\begin{array}{c}\hfill \alpha \circ X={\sum }_{i=1}^X{W}_i,\end{array}$

where X is a non-negative integer-valued random variable and $P(W_i=1)=1-P(W_i=0)=\alpha$. The INAR(1) model [2], The BAR(1) model [3], the INAR(p) model [4], The PDINAR(1) model [5] and the BARIO model [6] are very popular in analyzing non-negative integer-valued time series; see Weiß [7], Scotto et al. [8] and Davis et al. [9] for recent reviews on this topic. Motivated by infinite-patch metapopulation models discussed in Buckley and Pollett [10], Weiß [11] proposed an extension to the popular Poisson INAR(1) model, which is characterized by time-dependent innovations, i.e., the mean of the innovation is linearly increased by the previous population size. An important advantage of this model is that it gives a reasonable interpretation for immigration, which becomes more attractive if the current population is large; see Weiß [11] for an application to iceberg order data.

Univariate models are extensively investigated in the literature, but relatively few multivariate models, especially for bivariate versions, have been studied in detail. Franke and Rao [12] proposed a multivariate INAR(1) model, which is generalized to the p-order case by Latour [13]. Pedeli and Karlis [14] discussed a tractable bivariate INAR(1) model, which can be used to deal with bivariate count data with equi-dispersion and over-dispersion, but with small flexibility. See Pedeli and Karlis [15] for the estimation of the BINAR model and Pedeli and Karlis [16] for a further discussion of the properties of the multivariate INAR(1) model. Based on hierarchical dynamic models, Ravishanker et al. [17] described a Bayesian framework for estimation and prediction for multivariate times series of counts. Popović [18] proposed a bivariate INAR(1) model with random coefficients based on different binomial thinning operators. The above models assumed the innovations of their marginal models are independent and identically distributed. Based on a finite range of counts, Scotto et al. [19] considered the density-dependent bivariate binomial autoregressive models by using their state-dependent thinning concept. Li et al. [20] proposed a bivariate random coefficient INAR(1) model with asymmetric Hermite innovations.

Inspired by Weiß [11], we aim at providing a bivariate INAR model to analyze bivariate time series with time dependence and cross-correlation. The first contribution is that this paper gives an available method to capture the time-dependence trend by imposing the past information in the distribution of the innovational vector, which in turn makes the cross-correlation between the two entries into an innovation vector. The second contribution is that the new model not only allows autocorrelation but also allows cross-correlation. The third contribution is that this paper illustrates the stationarity and ergodicity of the extended bivariate INAR process and its two subprocesses, which is important to derive the consistence and asymptotic normality of the CML estimation.

The remainder of the paper is organized as follows. In Section 2, we first give brief reviews of the bivariate Poisson distribution and the bivariate binomial thinning operator, based on which we give the definition of the new bivariate INAR(1) model. Conditional maximum likelihood (CML) estimates and the asymptotic properties of unknown parameters are discussed in Section 3. A simulation and two real data examples that show the effectiveness of the new model are given in Section 4 and Section 5, respectively. Conclusions are made in Section 6.

2. A New Bivariate INAR(1) Model

For readability, we first give a brief review of the bivariate Poisson distribution.

Definition 1.

If the joint probability mass function (pmf) of $(X,Y)$ satisfies

(2)$\begin{array}{cc}\hfill & P(X=x,Y=y)\hfill \\ \hfill & ={\mathrm{e}}^{-({\lambda }_1+{\lambda }_2-\varphi )}{\displaystyle \frac{{({\lambda }_1-\varphi )}^x}{x!}}{\displaystyle \frac{{({\lambda }_2-\varphi )}^y}{y!}}\sum _{i=0}^{min(x,y)}\left(\genfrac{}{}{0pt}{}{x}{i}\right)\left(\genfrac{}{}{0pt}{}{y}{i}\right)i!{\left[{\displaystyle \frac{\varphi }{({\lambda }_1-\varphi )({\lambda }_2-\varphi )}}\right]}^i,\hfill \end{array}$

where ${\lambda }_1,{\lambda }_2>0$ and $\varphi \in (0,min({\lambda }_1,{\lambda }_2))$, then $(X,Y)$ is called a bivariate Poisson random variable with parameters $({\lambda }_1,{\lambda }_2,{\lambda }_3)$, i.e.,BP$({\lambda }_1,{\lambda }_2,\varphi ).$

From Kocherlakota and Kocherlakota [21], we obtain the fact that if $(X,Y)$ follows BP$({\lambda }_1,{\lambda }_2,\varphi )$, there must exist three mutually independent random variables $Z_1$, $Z_2$, $Z_3$ such that $X=Z_1+Z_3$ and $Y=Z_2+Z_3$, where $Z_1,Z_2$ and $Z_3$ follow Poisson$({\lambda }_1-\varphi )$, Poisson$({\lambda }_2-\varphi )$ and Poisson$\left(\varphi \right)$, respectively. Then, we have the conclusion that $\mathrm{Cov}(X,Y)=\varphi$. In addition, $P(X=x,Y=y)$, given in (2), is continuous and differentiable. For convenience, we denote $f(x,y,{\lambda }_1,{\lambda }_2,\varphi )=P(X=x,Y=y)$. By using Lemma A3 in Li et al. [20], we obtain that

(3)$\begin{array}{cc}\hfill {\displaystyle \frac{\partial f(x,y,{\lambda }_1,{\lambda }_2,\varphi )}{\partial {\lambda }_1}}& =f(x-1,y,{\lambda }_1,{\lambda }_2,\varphi )-f(x,y,{\lambda }_1,{\lambda }_2,\varphi ),\hfill \end{array}$

(4)$\begin{array}{cc}\hfill {\displaystyle \frac{\partial f(x,y,{\lambda }_1,{\lambda }_2,\varphi )}{\partial {\lambda }_2}}& =f(x,y-1,{\lambda }_1,{\lambda }_2,\varphi )-f(x,y,{\lambda }_1,{\lambda }_2,\varphi ),\hfill \end{array}$

(5)$\begin{array}{cc}\hfill {\displaystyle \frac{\partial f(x,y,{\lambda }_1,{\lambda }_2,\varphi )}{\partial \varphi }}& =f(x,y,{\lambda }_1,{\lambda }_2,\varphi )-f(x-1,y,{\lambda }_1,{\lambda }_2,\varphi )-f(x,y-1,{\lambda }_1,{\lambda }_2,\varphi )\hfill \\ \hfill & +f(x-1,y-1,{\lambda }_1,{\lambda }_2,\varphi ).\hfill \end{array}$

Applying the univariate binomial thinning operator $“\circ ”$ given in (1) to the bivariate case with $\mathit{X}={\left(X_1,X_2\right)}^{\top }$ leads to the bivariate binomial thinning operator:

$\mathit{A}\circ \mathit{X}=\left(\begin{array}{c}{\alpha }_{11}\circ X_1+{\alpha }_{12}\circ X_2\\ {\alpha }_{21}\circ X_1+{\alpha }_{22}\circ X_2\end{array}\right)\mathrm{with}\mathit{A}={\left({\alpha }_{ij}\right)}_{2\times 2},$

where ${\alpha }_{ij}\in (0,1),i,j=1,2$, $X_1$ and $X_2$ are non-negative integer-valued random variables, and all the thinnings are performed independent of each other.

By calculation, $E(\mathit{A}\circ \mathit{X})=\mathit{A}E\left(\mathit{X}\right).$ Denoting $\mathit{V}$ as the $2\times 2$ variance matrix of the Bernoulli random variables ${\alpha }_{ij}\circ X_j$ with ${\left(\mathit{V}\right)}_{ij}={\alpha }_{ij}(1-{\alpha }_{ij})$, $i,j=1,2,$ we obtain that $E\left((\mathit{A}\circ \mathit{X}){(\mathit{A}\circ \mathit{X})}^{\top }\right)=\mathit{A}E\left(\mathit{X}{\mathit{X}}^{\top }\right){\mathit{A}}^{\top }+\mathrm{diag}\left(\mathit{V}E\left(\mathit{X}\right)\right).$ Furthermore, if all the counting series of $\mathit{A}\circ \mathit{X}$ and $\mathit{B}\circ \mathit{Y}$ are independent, $E\left((\mathit{A}\circ \mathit{X}){(\mathit{B}\circ \mathit{Y})}^{\top }\right)=\mathit{A}E\left(\mathit{X}{\mathit{Y}}^{\top }\right){\mathit{B}}^{\top }$.

In the following, we give the definition of the new bivariate INAR(1) model, which not only includes the property of the models defined by Pedeli and Karlis [14,16] but also allows the innovation vectors $\left\{{\mathit{ϵ}}_t\right\}$ to be time-dependent.

Definition 2.

Let ${\mathit{X}}_t={\left(X_{1t},X_{2t}\right)}^{\top }$ be non-negative integer-valued bivariate random vector. If the process $\left\{{\mathit{X}}_t\right\}$ satisfies

(6) ${\mathit{X}}_t=\mathit{A}\circ {\mathit{X}}_{t-1}+{\mathit{ϵ}}_t,t\in \mathbb{Z},$

then $\left\{{\mathit{X}}_t\right\}$ is said to follow the extended bivariate INAR(1) process, where $\mathit{A}={\left({\alpha }_{ij}\right)}_{2\times 2}$, $0<{\alpha }_{ij}<1$, for any $i,j=1,2,{\mathit{ϵ}}_t={\left({ϵ}_{1t},{ϵ}_{2t}\right)}^{\top }\sim \mathrm{BP}({\lambda }_{1t},{\lambda }_{2t},\varphi )$ with ${({\lambda }_{1t},{\lambda }_{2t})}^{\top }=\mathit{B}{\mathit{X}}_{t-1}+\mathit{C},\mathit{B}={\left(b_{ij}\right)}_{2\times 2}$, $\mathit{C}={(c_1,c_2)}^{\top }$, $0<b_{ij}<1,c_i>0,i,j=1,2$.

For simplicity, we denote the new model as the EBINAR(1) model. It is easy to see that the ith equation of model (6) is presented by:

(7)$X_{it}={\alpha }_{i1}\circ X_{1,t-1}+{\alpha }_{i2}\circ X_{2,t-1}+{ϵ}_{it},i=1,2.$

Notice that the model given by (7) is similar to the one discussed in Weiß [11], the main difference is that $X_{it}$ involves two paralleled survivors $X_{1,t-1}$ and $X_{2,t-1}$. It is known that the EBINAR(1) process $\{{\mathit{X}}_t,t\in \mathbb{Z}\}$ has two parts: the first part consists of the survivors of the elements of the system at the preceding time $t-1$, denoted by ${\mathit{X}}_{t-1}$; the other part is comprised by the time-dependent innovation vector ${\mathit{ϵ}}_t$, which implies that the mean of the innovation vector is linearly increased by the previous population size.

Remark 1. (1). If both $\mathit{A}$ and $\mathit{B}$ are diagonal matrices, the component equation given in (7) becomes to the one discussed by Weiß [11].

(2). If $\mathit{A}$ is diagonal and $\mathit{B}=\mathbf{0}$, model (6) becomes the one discussed in Pedeli and Karlis [14], but it is worth mentioning that the autoregression matrix in Pedeli and Karlis [14] is diagonal, which means that it causes no cross-correlation in the counts.

(3). If $\mathit{A}$ is non-diagonal and $\mathit{B}=\mathbf{0}$, model (6) becomes the one discussed in Pedeli and Karlis [16], which accounts for cross-correlation in the counts, but they still keep the innovations of their marginal models independent and identically distributed such that the time dependence can not to be captured.

To derive the pmf the EBINAR(1) process, we first denote $h(k,m_1,m_2,{\alpha }_1,{\alpha }_2):=P(X+Y=k)$ is the convolution of $X+Y$, $\forall k\ge 0$ with $X\sim \mathrm{Bin}(m_1,{\alpha }_1)$ and $Y\sim \mathrm{Bin}(m_2,{\alpha }_2)$. By calculation, we obtain that

(8)$\begin{array}{c}\hfill h(k,m_1,m_2,{\alpha }_1,{\alpha }_2)={\sum }_{j=0}^{s}P(X=j|m_1,{\alpha }_1)P(Y=k-j|m_2,{\alpha }_2).\end{array}$

Furthermore, by using Lemma A3 in Li et al. [20],

(9)$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial h(k,m_1,m_2,{\alpha }_1,{\alpha }_2)}{\partial {\alpha }_1}}=m_1\left(h(k-1,m_1-1,m_2,{\alpha }_1,{\alpha }_2)-h(k,m_1-1,m_2,{\alpha }_1,{\alpha }_2)\right),\hfill \end{array}$

(10)$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial h(k,m_1,m_2,{\alpha }_1,{\alpha }_2)}{\partial {\alpha }_2}}=m_2\left(h(k-1,m_1,m_2-1,{\alpha }_1,{\alpha }_2)-h(k,m_1,m_2-1,{\alpha }_1,{\alpha }_2)\right).\hfill \end{array}$

Second, we denote $\mathit{\varsigma }={\left({\varsigma }_1,{\varsigma }_2\right)}^{\top }$, $\mathit{\vartheta }={\left({\vartheta }_1,{\vartheta }_2\right)}^{\top }$, $\mathit{k}={(k_1,k_2)}^{\top }$ and let $x={\varsigma }_1-k_1$ and $y={\varsigma }_2-k_2$. Then, the conditional probability distribution of the EBINAR(1) process takes the following form:

(11)$\begin{array}{cc}\hfill & P\left(\mathit{\varsigma }\right|\mathit{\vartheta }):=P({\mathit{X}}_t=\mathit{\varsigma }|{\mathit{X}}_{t-1}=\mathit{\vartheta })={\sum }_{k_1=0}^{g_1}{\sum }_{k_2=0}^{g_2}P(\mathit{A}\circ {\mathit{X}}_{t-1}=\mathit{k}|{\mathit{ϵ}}_t=\mathit{\varsigma }-\mathit{k})P({\mathit{ϵ}}_t=\varsigma -\mathit{k})\hfill \\ \hfill & ={\sum }_{k_1=0}^{g_1}{\sum }_{k_2=0}^{g_2}(P({\alpha }_{11}\circ X_{1,t-1}+{\alpha }_{12}\circ X_{2,t-1}=k_1)\times P({\alpha }_{21}\circ X_{1,t-1}+{\alpha }_{22}\circ X_{2,t-1}=k_2).\hfill \\ \hfill & \times P({ϵ}_{1t}={\varsigma }_1-k_1,{ϵ}_{2t}={\varsigma }_2-k_2))\hfill \\ \hfill & =\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}h(k_1,{\vartheta }_1,{\vartheta }_2,{\alpha }_{11},{\alpha }_{12})h(k_2,{\vartheta }_1,{\vartheta }_2,{\alpha }_{21},{\alpha }_{22})f(x,y,{\lambda }_{1t},{\lambda }_{2t},\varphi ),\hfill \end{array}$

where $g_1=min({\varsigma }_1,{\vartheta }_1)$, $g_2=min({\varsigma }_2,{\vartheta }_2)$,

$\begin{array}{cc}\hfill & f(x,y,{\lambda }_{1t},{\lambda }_{2t},\varphi )=P({ϵ}_{1t}=x,{ϵ}_{2t}=y)\hfill \\ \hfill & \stackrel{\left(2\right)}{=}exp({\lambda }_{1t}+{\lambda }_{2t}-\varphi ){\displaystyle \frac{{({\lambda }_{1t}-\varphi )}^x}{x!}}{\displaystyle \frac{{({\lambda }_{2t}-\varphi )}^y}{y!}}\sum _{i=0}^{min(x,y)}\left(\genfrac{}{}{0pt}{}{x}{i}\right)\left(\genfrac{}{}{0pt}{}{y}{i}\right)i!{\left[{\displaystyle \frac{\varphi }{({\lambda }_{1t}-\varphi )({\lambda }_{2t}-\varphi )}}\right]}^i\hfill \end{array}$

with ${\lambda }_{1t}=b_{11}{X}_{1,t-1}+b_{12}{X}_{2,t-1}+c_1$ and ${\lambda }_{2t}=b_{21}{X}_{1,t-1}+b_{22}{X}_{2,t-1}+c_2$.

If the largest eigenvalue of non-negative matrix $\mathit{A}$ is less than 1, then the bivariate marginal distribution of model (6) can be expressed in terms of the bivariate innovation vectors:

(12)${\mathit{X}}_t\stackrel{d}{=}{\mathit{A}}^k\circ {\mathit{X}}_{t-k}+{\sum }_{j=0}^{k-1}{\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j}={\mathit{A}}^t\circ {\mathit{X}}_0+{\sum }_{j=0}^{t-1}{\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j},k=1,2,\cdots ,t,$

where ${\mathit{A}}^0$ is an identity matrix, and ${\mathit{X}}_0$ is the initial value of the process.

In what follows, we first discuss the stationarity and ergodicity of processes (6) and (7), respectively. Second, we obtain the first two-moment of $\left\{{\mathit{X}}_t\right\}$ and $\left\{{\mathit{ϵ}}_t\right\}$, respectively. Third, we give a necessary and sufficient condition for the existence of $E{\left(X_{1t}\right)}^k$ and $E{\left(X_{2t}\right)}^k$ for any fixed positive integer k. These properties are necessary to derive the asymptotic properties of the estimators.

Theorem 1.

Let $\{{\mathit{X}}_t={(X_{1t},X_{2t})}^{\top }\}$ follow (6), $\mathbf{\Gamma }=\mathit{A}+\mathit{B}={\left({\gamma }_{ij}\right)}_{i,j=1,2}$ with $0<{\gamma }_{ij}<1$. If the largest eigenvalue of Γ is less than 1, there exists a strictly stationary and ergodic process satisfying (7).

Proof. 

Let $W_{t,k}$, $V_{t,l}$ and ${\delta }_{1t}$ be independent of each other and each of them be independent and identically distributed, i.e., $W_{t,k}\sim \mathrm{Bin}(1,{\alpha }_{11})+\mathrm{Poi}\left(b_{11}\right)$, $V_{t,l}\sim \mathrm{Bin}(1,{\alpha }_{12})+\mathrm{Poi}\left(b_{12}\right)$ and ${\delta }_{1t}\sim \mathrm{Poi}\left(c_1\right)$, where $\mathrm{Bin}(1,{\alpha }_{1j})+\mathrm{Poi}\left(b_{1j}\right)$ means the convolution of the distributions $\mathrm{Bin}(1,{\alpha }_{1j})$ and $\mathrm{Poi}\left(b_{1j}\right)$, $k=1,2,\cdots ,X_{1,t-1}$ and $l=1,2,\cdots ,X_{2,t-1}$; see Weiß [11] for details. According to the concepts of bivariate binomial thinning and the additivity of binomial distribution and Poisson distribution, (7) can be rewritten as

(13)$\begin{array}{c}\hfill X_{1t}\stackrel{d}{=}{W}_{t,1}+\cdots +W_{t,X_{1,t-1}}+V_{t,1}+\cdots +V_{t,X_{2,t-1}}+{\delta }_{1t}.\end{array}$

Since $\gamma =max({\alpha }_{ij}+b_{ij})<1$, we have $E\left(W_{t,k}\right)={\alpha }_{11}+b_{11}<1$ and $E\left(V_{t,l}\right)={\alpha }_{12}+b_{12}<1$. Denote $H\left(n\right)={\sum }_{k=1}^n\frac{1}{k}$ and $H\left(0\right)=0$, then $E\left(H\left({\delta }_{1t}\right)\right)={\sum }_{k=1}^{+\infty }\frac{1}{k}P({\delta }_{1t}\ge k).$ In addition, that $H\left({\delta }_{1t}\right)\le {\delta }_{1t}$ and $E\left({\delta }_{1t}\right)=c_1<\infty$, thus, $EH\left({\delta }_{1t}\right)\le E\left({\delta }_{1t}\right)<\infty$. Therefore, the Theorem of Heathcote [22] holds. Hence, there exists a stationary marginal distribution of (7), i.e., there exists a strictly stationary process satisfying (7). Similarly, we also have a similar conclusion for $X_{2t}$. □

To prove the stationaity of the EBINAR(1) process, we first introduce a sequence of random variables $\left\{{\mathit{X}}_t^{\left(n\right)}\right\}$ that could be considered as approximations to $\left\{{\mathit{X}}_t\right\}$ with

${\mathit{X}}_t^{\left(n\right)}=\left\{\begin{array}{cc}\mathbf{0},\hfill & n<0,\hfill \\ {\mathit{R}}_t,\hfill & n=0,\hfill \\ \mathit{A}\circ {\mathit{X}}_{t-1}^{(n-1)}+\mathit{B}{\mathit{X}}_{t-1}^{(n-1)}+{\mathit{R}}_t,\hfill & n>0,\hfill \end{array}\right.$

where the largest eigenvalues of the non-negative matrices $\mathit{A}$, $\mathit{B}$ and $\mathbf{\Gamma }:=\mathit{A}+\mathit{B}$ are less than 1, all of the non-negative matrices $\mathit{A}$, $\mathit{B}$ and $I-\mathbf{\Gamma }$ are invertible, ${\mathit{R}}_t={(R_{1t},R_{2t})}^{\top }$, $R_{1t}$ is independent with $R_{2t}$ and $R_{it}$ follows a Poisson distribution with the parameter $c_i$, $i=1,2$.

Theorem 2.

If the conditions of Theorem 1 hold, there exists a strictly stationary process satisfying (6).

Proof. 

Because

$\left(\begin{array}{c}{\mathit{X}}_1^{\left(0\right)}\\ ⋮\\ {\mathit{X}}_k^{\left(0\right)}\end{array}\right)=\left(\begin{array}{c}{\mathit{R}}_1\\ ⋮\\ {\mathit{R}}_k\end{array}\right)\mathrm{and}\left(\begin{array}{c}{\mathit{X}}_{h+1}^{\left(0\right)}\\ ⋮\\ {\mathit{X}}_{h+k}^{\left(0\right)}\end{array}\right)=\left(\begin{array}{c}{\mathit{R}}_{h+1}\\ ⋮\\ {\mathit{R}}_{h+k}\end{array}\right)$

are identically distributed for ${({\mathit{R}}_1,\cdots ,{\mathit{R}}_k)}^{\top }$ and ${({\mathit{R}}_{h+1},\cdots ,{\mathit{R}}_{h+k})}^{\top }$ are identically distributed. Thus, $\left\{{\mathit{X}}_t^{\left(0\right)}\right\}$ is strictly stationary. Now, we suppose $\left\{{\mathit{X}}_t^{\left(n\right)}\right\}$ is strictly stationary. Then,

(14)$\begin{array}{c}\hfill \left(\begin{array}{c}{\mathit{X}}_1^{(n+1)}\\ ⋮\\ {\mathit{X}}_k^{(n+1)}\end{array}\right)=\left(\begin{array}{c}{\mathit{R}}_1\\ ⋮\\ {\mathit{R}}_k\end{array}\right)+\left(\begin{array}{ccc}\mathit{A}& \cdots & \mathbf{0}\\ ⋮& \ddots & ⋮\\ \mathbf{0}& \cdots & \mathit{A}\end{array}\right)\circ \left(\begin{array}{c}{\mathit{X}}_0^{\left(n\right)}\\ ⋮\\ {\mathit{X}}_{k-1}^{\left(n\right)}\end{array}\right)+\left(\begin{array}{ccc}\mathit{B}& \cdots & \mathbf{0}\\ ⋮& \ddots & ⋮\\ \mathbf{0}& \cdots & \mathit{B}\end{array}\right)\left(\begin{array}{c}{\mathit{X}}_0^{\left(n\right)}\\ ⋮\\ {\mathit{X}}_{k-1}^{\left(n\right)}\end{array}\right)\end{array}$

and

(15)$\begin{array}{c}\hfill \left(\begin{array}{c}{\mathit{X}}_{h+1}^{(n+1)}\\ ⋮\\ {\mathit{X}}_{h+k}^{(n+1)}\end{array}\right)=\left(\begin{array}{c}{\mathit{R}}_{h+1}\\ ⋮\\ {\mathit{R}}_{h+k}\end{array}\right)+\left(\begin{array}{ccc}\mathit{A}& \cdots & \mathbf{0}\\ ⋮& \ddots & ⋮\\ \mathbf{0}& \cdots & \mathit{A}\end{array}\right)\circ \left(\begin{array}{c}{\mathit{X}}_h^{\left(n\right)}\\ ⋮\\ {\mathit{X}}_{h+k-1}^{\left(n\right)}\end{array}\right)+\left(\begin{array}{ccc}\mathit{B}& \cdots & \mathbf{0}\\ ⋮& \ddots & ⋮\\ \mathbf{0}& \cdots & \mathit{B}\end{array}\right)\left(\begin{array}{c}{\mathit{X}}_h^{\left(n\right)}\\ ⋮\\ {\mathit{X}}_{h+k-1}^{\left(n\right)}\end{array}\right).\end{array}$

Because the joint distributions of the variables involved in the right-hand sides of (14) and (15) are identical, thus, ${\left({\mathit{X}}_1^{(n+1)},\cdots ,{\mathit{X}}_k^{(n+1)}\right)}^{\top }$ and ${\left({\mathit{X}}_{h+1}^{(n+1)},\cdots ,{\mathit{X}}_{h+k}^{(n+1)}\right)}^{\top }$ in the left-hand side of the two equations are identically distributed. Hence, the process $\left\{{\mathit{X}}_t^{\left(n\right)}\right\}$ is strictly stationary. Therefore, $\left\{{\mathit{X}}_t\right\}$ is a strictly stationary process. □

Theorem 3.

If the conditions of Theorem 1 hold, $\left\{{\mathit{X}}_t\right\}$ is a null recurrent and ergodic Markov chain.

Proof. 

First, we prove $\left\{{\mathit{X}}_t\right\}$ is null recurrent. Because ${\mathit{X}}_t\stackrel{d}{=}{\sum }_{j=0}^{\infty }{\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j},$ then $P_{\mathbf{0},\mathbf{0}}^t=P({\mathit{X}}_t=\mathbf{0}|{\mathit{X}}_0=\mathbf{0})={\prod }_{j=0}^{k-1}P({\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j}=\mathbf{0})$ with probability one. Let ${\mathit{A}}^j=\left(p_{ij}\right)$ and $\gamma =max\left(p_{ij}\right)$, $\forall i=1,2,\forall j=1,2$. Then, we obtain that

$\begin{array}{cc}\hfill & P({\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j}\ne \mathbf{0})=P(\{p_{11}\circ {ϵ}_{1,t-j}+p_{12}\circ {ϵ}_{2,t-j}\ge 1\}\cup \{p_{21}\circ {ϵ}_{1,t-j}+p_{22}\circ {ϵ}_{2,t-j}\ge 1\})\hfill \\ \hfill & \le P(p_{11}\circ {ϵ}_{1,t-j}+p_{12}\circ {ϵ}_{2,t-j}\ge 1)+P(p_{21}\circ {ϵ}_{1,t-j}+p_{22}\circ {ϵ}_{2,t-j}\ge 1)\hfill \\ \hfill & \le P(p_{11}\circ {ϵ}_{1,t-j}\ge 1)+P(p_{12}\circ {ϵ}_{2,t-j}\ge 1)+P(p_{21}\circ {ϵ}_{1,t-j}\ge 1)+P(p_{22}\circ {ϵ}_{2,t-j}\ge 1)\hfill \\ \hfill & \le 2\left[P({\gamma }^j\circ {ϵ}_{1,t-j}\ge 1)+P({\gamma }^j\circ {ϵ}_{2,t-j}\ge 1)\right]\hfill \\ \hfill & \le 2\left[E({\gamma }^j\circ {ϵ}_{1,t-j})+E({\gamma }^j\circ {ϵ}_{2,t-j})\right]=2{\gamma }^j({\mu }_{{ϵ}_1}+{\mu }_{{ϵ}_2}).\hfill \end{array}$

According to Theorem 2, there exists an $M>0$ such that ${\mu }_{{ϵ}_i}\le M/4,i=1,2.$ Then, we have $P({\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j}=\mathbf{0})\ge 1-M{\gamma }^j.$ Hence,

$\begin{array}{cc}\hfill P_{\mathbf{0},\mathbf{0}}^t\ge {\prod }_{j=0}^{t-1}P({\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j}=\mathbf{0})& \ge {\prod }_{j=0}^{t-1}(1-M{\gamma }^j)=exp\{{\sum }_{j=0}^{t-1}log(1-M{\gamma }^j)\}\hfill \\ \hfill & \ge exp\{log(1-M{\gamma }^{\tau })/\left(1-r\right)\}>0,\forall t>\tau .\hfill \end{array}$

Therefore, $\underset{t\to \infty }{lim}{P}_{\mathbf{0},\mathbf{0}}^t>0$. Thus, ${\sum }_{k=0}^{\infty }{P}_{\mathbf{0},\mathbf{0}}^t=\infty$, i.e., $\mathbf{0}$ is a recurrent state.

Second, we illustrate the ergodicity. For all states $\mathit{\varsigma }$, $\mathit{\vartheta }$, ${\mathit{\kappa }}_{t-2},{\mathit{\kappa }}_{t-3},\cdots$, we have

$P(\left.{\mathit{X}}_t=\mathit{\varsigma }\right|{\mathit{X}}_{t-1}=\mathit{\vartheta },{\mathit{X}}_{t-2}={\mathit{\kappa }}_{t-2},\cdots )=P(\left.{\mathit{X}}_t=\mathit{\varsigma }\right|{\mathit{X}}_{t-1}=\mathit{\vartheta })=P(\mathit{\vartheta },\mathit{\varsigma }),$

where $P(\mathit{\vartheta },\mathit{\varsigma })$ denotes the transition probability from state $\mathit{\vartheta }$ to state $\mathit{\varsigma }$. Thus, $\left\{{\mathit{X}}_t\right\}$ is a homogeneous Markov chain. Since ${\alpha }_{ij},b_{ij}\in (0,1)$, thus $P({ϵ}_{1t}={\varsigma }_1-k_1,{ϵ}_{2t}={\varsigma }_2-k_2)>0.$ Denote n-state transition probability from state $\mathit{\varsigma }$ to state $\mathit{\vartheta }$ with $P_{\mathit{\varsigma }\mathit{\vartheta }}^n$. For a given ${\mathit{X}}_{t-1}$, the conditional probability function of the random vector ${\mathit{X}}_t$ is derived by:

$\begin{array}{cc}\hfill P(X_{1t}& ={\varsigma }_1,X_{2t}={\varsigma }_2|X_{1,t-1}={\vartheta }_1,X_{2,t-1}={\vartheta }_2)\hfill \\ \hfill & =P(X_{1t}={\varsigma }_1|X_{1,t-1}={\vartheta }_1,X_{2,t-1}={\vartheta }_2)P(X_{2t}={\varsigma }_2|X_{1,t-1}={\vartheta }_1,X_{2,t-1}={\vartheta }_2,X_{1t}={\varsigma }_1),\hfill \end{array}$

then $P_{\mathit{\tau }\mathit{\upsilon }}^1>0$ for all $\mathit{\tau },\mathit{\upsilon }\in {\mathbb{N}}_0^2$. According to (12), every state can be reached from any other state with positive probability in a finite number of steps, analogously. Hence, $\left\{{\mathit{X}}_t\right\}$ is irreducible. By (12), k steps of conditional probability distribution $P_{\mathbf{0},\mathbf{0}}^k$ are obtained with:

$\begin{array}{cc}\hfill P_{\mathbf{0},\mathbf{0}}^k=P({\mathit{X}}_t=\mathbf{0}|{\mathit{X}}_{t-k}=\mathbf{0})& =P({\mathit{A}}^k\circ {\mathit{X}}_{t-k}+{\sum }_{j=0}^{k-1}{\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j}=\mathbf{0}|{\mathit{X}}_{t-k}=\mathbf{0})\hfill \\ \hfill & =\underset{\left(\mathrm{V}\right)}{\underbrace{P({\mathit{A}}^k\circ {\mathit{X}}_{t-k}=\mathbf{0}|{\mathit{X}}_{t-k}=\mathbf{0})}}\underset{\left(\mathrm{VI}\right)}{\underbrace{{\prod }_{j=0}^{k-1}P({\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j}=\mathbf{0})}}.\hfill \end{array}$

Note that the first multiplier (V) is positive, which can be obtained by a similar method to (11). Denoting ${\mathit{A}}^j={\left(p_{ij}\right)}_{i,j=1,2}$, then we have:

$\begin{array}{cc}\hfill & P({\mathit{A}}^j\circ {\mathit{ϵ}}_{t-j}=\mathbf{0})=P(p_{11}\circ {ϵ}_{1,t-j}+p_{12}\circ {ϵ}_{2,t-j}=0,p_{21}\circ {ϵ}_{1,t-j}+p_{22}\circ {ϵ}_{2,t-j}=0)\hfill \\ \hfill & ={\sum }_{k=0}^{\infty }{\sum }_{s=0}^{\infty }{(1-p_{11})}^k{(1-p_{12})}^s{(1-p_{21})}^k{(1-p_{22})}^{s}P({ϵ}_{1,t-j}=k,{ϵ}_{2,t-j}=s)>0,\hfill \end{array}$

thus, the second part, (VI), is also positive. Therefore, $P_{\mathbf{0},\mathbf{0}}^k>0$, with probability one, i.e., $\left\{{\mathit{X}}_t\right\}$, is aperiodic. Hence, $\left\{{\mathit{X}}_t\right\}$ is an ergodic Markov chain. □

Note that $E\left({\mathit{X}}_t^{\left(n\right)}\right)={(\mathit{I}-\mathit{A}-\mathit{B})}^{-1}\mathit{C}<\infty$ and

$\begin{array}{cc}\hfill & E\left({\mathit{X}}_t^{\left(n\right)}{{\mathit{X}}_t^{\left(n\right)}}^{\top }\right)=\mathbf{\Gamma }E\left({\mathit{X}}_t^{(n-1)}{{\mathit{X}}_t^{(n-1)}}^{\top }\right){\mathbf{\Gamma }}^{\top }+\mathbf{\Psi }\hfill \\ \hfill & =\cdots ={\mathbf{\Gamma }}^{n}E\left({\mathit{X}}_t^{\left(0\right)}{{\mathit{X}}_t^{\left(0\right)}}^{\top }\right){\left({\mathbf{\Gamma }}^n\right)}^{\top }+{\mathbf{\Gamma }}^{n-1}\mathbf{\Psi }{\left({\mathbf{\Gamma }}^{n-1}\right)}^{\top }+\cdots +\mathbf{\Gamma }\mathbf{\Psi }{\mathbf{\Gamma }}^{\top }+\mathbf{\Psi },\hfill \end{array}$

where $\mathbf{\Psi }$ involves the first moments of ${\mathit{X}}_t^{\left(n\right)}$ and ${\mathit{R}}_t$. Hence, the first two moments of ${\mathit{X}}_t^{\left(n\right)}$ are finite. Thus, $\left\{{\mathit{X}}_t^{\left(n\right)}\right\}$ is stationary and ergodic by Theorem 2, Theorem 3 and Shumway and Stoffer [23].

Theorem 4.

If the conditions of Theorem 1 hold, the first two moments and covariance matrix of $\left\{{\mathit{X}}_t\right\}$ exist and

(1). $E\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})=(\mathit{A}+\mathit{B}){\mathit{X}}_{t-1}+\mathit{C}$;

(2). $E\left({\mathit{X}}_t\right)={(\mathit{I}-\mathit{A}-\mathit{B})}^{-1}\mathit{C}$, if ${(\mathit{I}-\mathit{A}-\mathit{B})}^{-1}$ exists, where $\mathit{I}$ denotes the identity matrix;

(3). $R\left(k\right)=\mathrm{Cov}\left({\mathit{X}}_{t+k},{\mathit{X}}_t\right)={(\mathit{A}+\mathit{B})}^{k}R\left(0\right)$, $k=1,2,\cdots$.

In addition, if $k=0$, $R\left(0\right)=\mathit{A}R\left(0\right){\mathit{A}}^{\top }+{\mathit{H}}^{*}+\mathit{A}R\left(0\right){\mathit{B}}^{\top }+\mathit{B}R\left(0\right){\mathit{A}}^{\top }+\mathrm{\Sigma },$ where ${\mathit{H}}^{*}=\mathrm{diag}\left({\sum }_{j=1}^M{\mathit{V}}_{ij}E\left(X_{j,t-1}\right)\right)$, ${\mathit{V}}_{ij}={\alpha }_{ij}(1-{\alpha }_{ij})$ and $\mathrm{\Sigma }=\mathrm{Cov}({\mathit{ϵ}}_t,{\mathit{ϵ}}_t)$. Specifically, if $\mathit{A}$ and $\mathit{B}$ are diagonal matrices, $R\left(0\right)={(\mathit{I}-\mathit{A}{\mathit{A}}^{\top }-2\mathit{A}{\mathit{B}}^{\top })}^{-1}\mathrm{\Sigma }+{\mathit{H}}^{*}.$

Proof. 

(1) and (2) are easy to prove by the moment property of $\mathit{A}\circ$, and we omit them. Here, we only give the proof of (3):

$\begin{array}{cc}\hfill R\left(k\right)& =\mathrm{Cov}\left(\mathit{A}\circ {\mathit{X}}_{t+k-1}+{\mathit{ϵ}}_{t+k},{\mathit{X}}_t\right)=\mathrm{Cov}\left(\mathit{A}\circ {\mathit{X}}_{t+k-1},{\mathit{X}}_t\right)+\mathrm{Cov}\left({\mathit{ϵ}}_{t+k},{\mathit{X}}_t\right)\hfill \\ \hfill & =\mathit{A}\mathrm{Cov}\left({\mathit{X}}_{t+k-1},{\mathit{X}}_t\right)+\mathrm{Cov}\left(\mathit{B}{\mathit{X}}_{t+k-1}+\mathit{C},{\mathit{X}}_t\right)=(\mathit{A}+\mathit{B})\mathrm{Cov}\left({\mathit{X}}_{t+k-1},{\mathit{X}}_t\right)\hfill \\ \hfill & =(\mathit{A}+\mathit{B})\mathit{R}(k-1)=\cdots ={(\mathit{A}+\mathit{B})}^{k}R\left(0\right).\hfill \end{array}$

In fact, $\mathrm{Cov}\left({\mathit{X}}_{t-1},{\mathit{ϵ}}_t\right)=\mathrm{Cov}\left({\mathit{X}}_{t-1},\mathit{B}{\mathit{X}}_{t-1}+\mathit{C}\right)=\mathrm{Cov}\left({\mathit{X}}_{t-1},{\mathit{X}}_{t-1}\right){\mathit{B}}^{\top }$ and $\mathrm{Cov}\left({\mathit{ϵ}}_t,{\mathit{X}}_{t-1}\right)=\mathrm{Cov}\left(\mathit{B}{\mathit{X}}_{t-1}+\mathit{C},{\mathit{X}}_{t-1}\right)=\mathit{B}\mathrm{Cov}\left({\mathit{X}}_{t-1},{\mathit{X}}_{t-1}\right)$. Hence,

$\begin{array}{cc}\hfill R\left(0\right)& =\mathrm{Cov}\left({\mathit{X}}_t,{\mathit{X}}_t\right)=\mathrm{Cov}\left(\mathit{A}\circ {\mathit{X}}_{t-1}+{\mathit{ϵ}}_t,\mathit{A}\circ {\mathit{X}}_{t-1}+{\mathit{ϵ}}_t\right)\hfill \\ \hfill & =\mathit{A}\mathrm{Cov}\left({\mathit{X}}_{t-1},{\mathit{X}}_{t-1}\right){\mathit{A}}^{\top }+{\mathit{H}}^{*}+\mathit{A}\mathrm{Cov}\left({\mathit{X}}_{t-1},{\mathit{X}}_{t-1}\right){\mathit{B}}^{\top }+\mathit{B}\mathrm{Cov}\left({\mathit{X}}_{t-1},{\mathit{X}}_{t-1}\right){\mathit{A}}^{\top }+\mathrm{\Sigma }\hfill \\ \hfill & =\mathit{A}R\left(0\right){\mathit{A}}^{\top }+{\mathit{H}}^{*}+\mathit{A}R\left(0\right){\mathit{B}}^{\top }+\mathit{B}R\left(0\right){\mathit{A}}^{\top }+\mathrm{\Sigma },\hfill \end{array}$

where ${\mathit{H}}^{*}=\mathrm{diag}\left({\sum }_{j=1}^2{\mathit{V}}_{ij}E\left(X_{j,t-1}\right)\right).$ Let $\lambda$, ${\lambda }_1$ and ${\lambda }_2$ be the largest eigenvalues of ${\mathit{AA}}^{\top }+2{\mathit{AB}}^{\top }$, $\mathit{A}$ and $\mathit{B}$, respectively. If $\mathit{A}$ and $\mathit{B}$ are diagonal matrices,

$\left|\lambda \right|\le |{\lambda }_1^2+2{\lambda }_1{\lambda }_2|\le |{\lambda }_1({\lambda }_1+{\lambda }_2)+{\lambda }_1{\lambda }_2|\le {\lambda }_1|({\lambda }_1+{\lambda }_2)|+{\lambda }_2\left|{\lambda }_1\right|\le {\lambda }_1+{\lambda }_2<1,$

then $\mathit{I}-{\mathit{AA}}^{\top }-2{\mathit{AB}}^{\top }$ is a nonsingular matrix. Hence, $R\left(0\right)$ is obtained. □

Theorem 5.

If the conditions of Theorem 1 hold, the first two moments and covariance matrices of $\left\{{\mathit{ϵ}}_t\right\}$ exist and:

(1). $E\left({\mathit{ϵ}}_t\right|{\mathit{X}}_{t-1})=\mathit{B}{\mathit{X}}_{t-1}+\mathit{C}$;

(2). $E\left({\mathit{ϵ}}_t\right)=\left(\mathit{I}-\mathit{A}\right){(\mathit{I}-\mathit{A}-\mathit{B})}^{-1}\mathit{C},$ if ${(\mathit{I}-\mathit{A}-\mathit{B})}^{-1}$ exists;

(3). $R_{\mathit{ϵ}}\left(k\right)=\mathrm{Cov}\left({\mathit{ϵ}}_{t+k},{\mathit{ϵ}}_t\right)=\mathit{B}{\left(\mathit{A}+\mathit{B}\right)}^{k}R\left(0\right){\mathit{B}}^{\top }$, $k=0,1,2,\cdots$.

Proof. 

(1) is easy to prove by the distribution of ${\mathit{ϵ}}_t$. We only need to prove (2) and (3).

(2). $E\left({\mathit{ϵ}}_t\right)=\mathit{B}{\mathit{\mu }}_t+\mathit{C}=\mathit{B}{(\mathit{I}-\mathit{A}-\mathit{B})}^{-1}\mathit{C}+(\mathit{I}-\mathit{A}-\mathit{B}){(\mathit{I}-\mathit{A}-\mathit{B})}^{-1}\mathit{C}=(\mathit{I}-\mathit{A})E\left({\mathit{X}}_t\right)$ by the definition of ${\mathit{ϵ}}_t$. $E\left({\mathit{ϵ}}_t\right)$ can be obtained directly by (6).

(3). According to the construction of the EBINAR(1) model, we have:

(16)$\begin{array}{c}\mathrm{Cov}\left({\mathit{X}}_t,{\mathit{ϵ}}_t\right)=\mathrm{Cov}\left(\mathit{A}\circ {\mathit{X}}_{t-1},\mathit{B}{\mathit{X}}_{t-1}+\mathit{C}\right)+\mathrm{Cov}\left(\mathit{B}{\mathit{X}}_{t-1}+\mathit{C},\mathit{B}{\mathit{X}}_{t-1}+\mathit{C}\right)\hfill \\ =\mathit{A}\mathrm{Cov}\left({\mathit{X}}_{t-1},{\mathit{X}}_{t-1}\right){\mathit{B}}^{\top }+\mathit{B}\mathrm{Cov}\left({\mathit{X}}_{t-1},{\mathit{X}}_{t-1}\right){\mathit{B}}^{\top }=(\mathit{A}+\mathit{B})R\left(0\right){\mathit{B}}^{\top },\hfill \end{array}$

(17)$\begin{array}{c}{R}_{\mathit{ϵ}}\left(k\right)=\mathrm{Cov}\left(\mathit{B}{\mathit{X}}_{t+k-1}+\mathit{C},{\mathit{ϵ}}_t\right)=\mathit{B}\mathrm{Cov}\left({\mathit{X}}_{t+k-1},{\mathit{ϵ}}_t\right)\hfill \\ =\mathit{B}\mathrm{Cov}\left(\mathit{A}\circ {\mathit{X}}_{t+k-2}+{\mathit{ϵ}}_{t+k-1},{\mathit{ϵ}}_t\right)=\mathit{BA}\mathrm{Cov}\left({\mathit{X}}_{t+k-2},{\mathit{ϵ}}_t\right)+\mathit{B}\mathrm{Cov}\left({\mathit{ϵ}}_{t+k-1},{\mathit{ϵ}}_t\right)\hfill \\ =\mathit{BA}\mathrm{Cov}\left({\mathit{X}}_{t+k-2},{\mathit{ϵ}}_t\right)+\mathit{B}\mathrm{Cov}\left(\mathit{B}{\mathit{X}}_{t+k-2}+\mathit{C},{\mathit{ϵ}}_t\right)\hfill \\ =\mathit{BA}\mathrm{Cov}\left({\mathit{X}}_{t+k-2},{\mathit{ϵ}}_t\right)+{\mathit{B}}^2\mathrm{Cov}\left({\mathit{X}}_{t+k-2},{\mathit{ϵ}}_t\right)=\mathit{B}\left(\mathit{A}+\mathit{B}\right)\mathrm{Cov}\left({\mathit{X}}_{t+k-2},{\mathit{ϵ}}_t\right)\hfill \\ =\cdots =\mathit{B}{\left(\mathit{A}+\mathit{B}\right)}^{k-1}\mathrm{Cov}\left({\mathit{X}}_t,{\mathit{ϵ}}_t\right),\hfill \end{array}$

then $R_{\mathit{ϵ}}\left(k\right)$ is achieved by substituting (16) into (17), i.e., $R_{\mathit{ϵ}}\left(k\right)=\mathit{B}{(\mathit{A}+\mathit{B})}^{k}R\left(0\right){\mathit{B}}^{\top }$. Note that $\mathrm{Cov}\left({\mathit{ϵ}}_t,{\mathit{ϵ}}_t\right)=\mathrm{Cov}\left(\mathit{B}{\mathit{X}}_{t-1}+\mathit{C},\mathit{B}{\mathit{X}}_{t-1}+\mathit{C}\right)=\mathit{B}R\left(0\right){\mathit{B}}^{\top },$ i.e., the formula holds for $k=0$. □

Theorem 6.

For any fixed positive integer k, it is a necessary and sufficient condition that $E{\left(X_{it}\right)}^k<\infty$ is $\gamma <1$, $i=1,2$.

Proof. 

For convenience, let $\mathit{A}$ and $\mathit{B}$ be diagonal matrices.

Necessity. According to Lemma 2.1 of Silva and Oliveira [24], $E\left[{({\alpha }_{11}\circ X_{1,t-1})}^i{\left({ϵ}_{1t}\right)}^{k-i}\right]={\alpha }_{11}^i{b}_{11}^{k-i}E{\left(X_{1,t-1}\right)}^k+{\psi }_1$, where ${\psi }_1={\psi }_1\left(X_{1,t-1}\right)$ involves the moments of $X_{1,t-1}$ of order $\le (k-1)$ and $i=0,1,2,\cdots ,k$. Then,

(18)$\begin{array}{cc}\hfill & E{\left(X_{1t}\right)}^k=E{({\alpha }_{11}\circ X_{1,t-1}+{ϵ}_{1t})}^k={\sum }_{i=0}^k\left(\genfrac{}{}{0pt}{}{k}{i}\right)E\left[{({\alpha }_{11}\circ X_{1,t-1})}^i{\left({ϵ}_{1t}\right)}^{k-i}\right]\hfill \\ \hfill & ={\sum }_{i=0}^k\left(\genfrac{}{}{0pt}{}{k}{i}\right){\alpha }_{11}^i{b}_{11}^{k-i}E{\left(X_{1,t-1}\right)}^k+{\psi }_1={({\alpha }_{11}+b_{11})}^{k}E{\left(X_{1,t-1}\right)}^k+{\psi }_1.\hfill \end{array}$

Thus, $E{\left(X_{1t}\right)}^k={\displaystyle \frac{{\psi }_1}{1-{({\alpha }_{11}+b_{11})}^k}}$ by (18). Hence, $1-{({\alpha }_{11}+b_{11})}^k>0$ if $E{\left(X_{1t}\right)}^k<\infty$, i.e., ${\alpha }_{11}+b_{11}<1$. Similarly, ${\alpha }_{22}+b_{22}<1$ if $E{\left(X_{2t}\right)}^k<\infty$. Hence, $\gamma <1$ if $E{\left(X_{it}\right)}^k<\infty$, $i=1,2$.

Sufficiency. We know that $E{\left(X_{it}\right)}^k<\infty$ holds for $k=1,2$ by Theorems 4 and 5. The sufficient condition can be proved by induction with respect to k. Now suppose that $E{\left(X_{it}\right)}^{k-1}<\infty$, $k\ge 3$. According to (13), we define

$X_{1t}^{\left(n\right)}=\left\{\begin{array}{cc}0,\hfill & n<0;\hfill \\ {\delta }_{1t},\hfill & n=0;\hfill \\ {\delta }_{1t}+{\sum }_{j=1}^{X_{1,t-1}^{(n-1)}}{W}_{t,j},\hfill & n>0\hfill \end{array}\right.\mathrm{and}{X}_{2t}^{\left(n\right)}=\left\{\begin{array}{cc}0,\hfill & n<0;\hfill \\ {\delta }_{2t},\hfill & n=0;\hfill \\ {\delta }_{2t}+{\sum }_{s=1}^{X_{2,t-1}^{(n-1)}}{V}_{t,s},\hfill & n>0,\hfill \end{array}\right.$

where $W_{t,j}$, ${\delta }_{1t}$, $V_{t,s}$ and ${\delta }_{2t}$ are independent of each other and each of them is independent and identically distributed, i.e., $W_{t,j}\sim \mathrm{Bin}(1,{\alpha }_{11})+\mathrm{Poi}\left(b_{11}\right)$, ${\delta }_{1t}\sim \mathrm{Poi}\left(c_1\right)$, $V_{t,s}\sim \mathrm{Bin}(1,{\alpha }_{22})+\mathrm{Pois}\left(b_{22}\right)$ and ${\delta }_{2t}\sim \mathrm{Pois}\left(c_2\right)$. Using the univariate binomial thinning operator, $X_{1t}^{\left(n\right)}$ and $X_{2t}^{\left(n\right)}$ admit the representations:

(19)$\begin{array}{cc}\hfill & X_{1t}^{\left(n\right)}={\delta }_{1t}+({\alpha }_{11}\circ X_{1,t-1}^{(n-1)}+Z_{1t}),\hfill \end{array}$

(20)$\begin{array}{cc}\hfill & X_{2t}^{\left(n\right)}={\delta }_{2t}+({\alpha }_{22}\circ X_{2,t-1}^{(n-1)}+Z_{2t}),\hfill \end{array}$

where $Z_{1t}\sim \mathrm{Poi}\left(b_{11}{X}_{1,t-1}^{(n-1)}\right)$ and $Z_{2t}\sim \mathrm{Poisson}\left(b_{22}{X}_{2,t-1}^{(n-1)}\right)$. It is easy to see both ${\left\{X_{1t}^{\left(n\right)}\right\}}_{n\in \mathbb{N}}$ and ${\left\{X_{2t}^{\left(n\right)}\right\}}_{n\in \mathbb{N}}$ are non-decreasing. According to Lemma 2.1 of [24], we have:

$\begin{array}{cc}\hfill & E{({\alpha }_{11}\circ X_{1,t-1}^{(n-1)}+Z_{1t})}^k={({\alpha }_{11}+b_{11})}^{k}E{\left(X_{1,t-1}^{(n-1)}\right)}^k+{\psi }_2\le {({\alpha }_{11}+b_{11})}^{k}E{\left(X_{1,t-1}^{\left(n\right)}\right)}^k+{\psi }_4,\hfill \\ \hfill & E{({\alpha }_{22}\circ X_{2,t-1}^{(n-1)}+Z_{2t})}^k={({\alpha }_{22}+b_{22})}^{k}E{\left(X_{2,t-1}^{(n-1)}\right)}^k+{\psi }_3\le {({\alpha }_{22}+b_{22})}^{k}E{\left(X_{2,t-1}^{\left(n\right)}\right)}^k+{\psi }_5,\hfill \end{array}$

where ${\psi }_2={\psi }_2\left(X_{1,t-1}^{(n-1)}\right)$ and ${\psi }_3={\psi }_3\left(X_{2,t-1}^{(n-1)}\right)$ involve the moments of $X_{1,t-1}^{(n-1)}$ and $X_{2,t-1}^{(n-1)}$ of order $\le (k-1)$, and ${\psi }_4={\psi }_4\left(X_{1,t-1}^{\left(n\right)}\right)$ and ${\psi }_5={\psi }_5\left(X_{2,t-1}^{\left(n\right)}\right)$ involve the moments of $X_{1,t-1}^{\left(n\right)}$ and $X_{2,t-1}^{\left(n\right)}$ of order $\le (k-1)$, respectively. According to (19) and (20), we obtain:

(21)$\begin{array}{cc}\hfill E{\left(X_{1t}^{\left(n\right)}\right)}^k& =E{\left({\delta }_{1t}\right)}^k+{({\alpha }_{11}\circ X_{1,t-1}^{(n-1)}+Z_{1t})}^k+\sum _{j=1}^{k-1}\left(\genfrac{}{}{0pt}{}{k}{j}\right)E{\left({\delta }_{1t}\right)}^{k-j}E{({\alpha }_{11}\circ X_{1,t-1}^{(n-1)}+Z_{1t})}^j\hfill \\ \hfill & \le E{\left({\delta }_{1t}\right)}^k+{({\alpha }_{11}+b_{11})}^{k}E{\left(X_{1,t-1}^{\left(n\right)}\right)}^k+{\psi }_4+\sum _{j=1}^{k-1}\left(\genfrac{}{}{0pt}{}{k}{j}\right)E{\left({\delta }_{1t}\right)}^{k-j}E{({\alpha }_{11}\circ X_{1,t-1}^{\left(n\right)}+Z_{1t})}^j\hfill \\ \hfill & \le c_1^k+{\gamma }^{k}E{\left(X_{1,t-1}^{\left(n\right)}\right)}^k+{\psi }_4+\sum _{j=1}^{k-1}\left(\genfrac{}{}{0pt}{}{k}{j}\right)E{\left({\delta }_{1t}\right)}^{k-j}E{({\alpha }_{11}\circ X_{1,t-1}^{\left(n\right)}+Z_{1t})}^j,\hfill \end{array}$

(22)$\begin{array}{cc}\hfill E{\left(X_{2t}^{\left(n\right)}\right)}^k& =E{\left({\delta }_{2t}\right)}^k+{({\alpha }_{22}\circ X_{2,t-1}^{(n-1)}+Z_{2t})}^k+\sum _{j=1}^{k-1}\left(\genfrac{}{}{0pt}{}{k}{j}\right)E{\left({\delta }_{2t}\right)}^{k-j}E{({\alpha }_{22}\circ X_{2,t-1}^{(n-1)}+Z_{2t})}^j\hfill \\ \hfill & \le E{\left({\delta }_{2t}\right)}^k+{({\alpha }_{22}+b_{22})}^{k}E{\left(X_{2,t-1}^{\left(n\right)}\right)}^k+{\psi }_5+\sum _{j=1}^{k-1}\left(\genfrac{}{}{0pt}{}{k}{j}\right)E{\left({\delta }_{2t}\right)}^{k-j}E{({\alpha }_{22}\circ X_{2,t-1}^{\left(n\right)}+Z_{2t})}^j\hfill \\ \hfill & \le c_2^k+{\gamma }^{k}E{\left(X_{2,t-1}^{\left(n\right)}\right)}^k+{\psi }_5+\sum _{j=1}^{k-1}\left(\genfrac{}{}{0pt}{}{k}{j}\right)E{\left({\delta }_{2t}\right)}^{k-j}E{({\alpha }_{22}\circ X_{2,t-1}^{\left(n\right)}+Z_{2t})}^j.\hfill \end{array}$

Using (21) and (22),

(23)$\begin{array}{c}\hfill E{\left(X_{1t}^{\left(n\right)}\right)}^k+E{\left(X_{2t}^{\left(n\right)}\right)}^k\le {\displaystyle \frac{{\sum }_{i=1}^2\left[c_i^k+{\sum }_{j=1}^{k-1}\left(\genfrac{}{}{0pt}{}{k}{j}\right)E{\left({\delta }_{it}\right)}^{k-j}E{({\alpha }_{ii}\circ X_{i,t-1}^{\left(n\right)}+Z_{it})}^j\right]+{\psi }_6}{1-{\gamma }^k}},\end{array}$

where ${\psi }_6={\psi }_4+{\psi }_5$. Note that the numerator in (23) involves the moments of $X_{1,t-1}^{\left(n\right)}$ and $X_{2,t-1}^{\left(n\right)}$ of order $\le k-1$ and is finite; thus, $E{\left(X_{1t}^{\left(n\right)}\right)}^k+E{\left(X_{2t}^{\left(n\right)}\right)}^k$ is finite if $\gamma <1$. In addition that both $E\left(X_{1t}^{\left(n\right)}\right)$ and $E\left(X_{2t}^{\left(n\right)}\right)$ are non-negative; thus, $E\left(X_{1t}^{\left(n\right)}\right)$ and $E\left(X_{2t}^{\left(n\right)}\right)$ are finite. □

3. Parameter Estimation

In this section, we consider the conditional maximum likelihood estimation for model (6). Let $\mathit{\theta }={\left({\alpha }_{ij},b_{ij},c_i,\varphi \right)}^{\top }$, $i,j=1,2$. Suppose that ${\mathit{X}}_0,{\mathit{X}}_1,\cdots ,{\mathit{X}}_T$ are generated by the EBINAR(1) model with the true parameter value ${\mathit{\theta }}_0$.

By (11), the conditional log-likelihood function can be written as:

(24)$\ell \left(\mathit{\theta }\right)={\sum }_{t=1}^{T}ln{P}_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right),$

where

$\begin{array}{cc}\hfill & P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)=P(X_{1t}=X_{1t},X_{2t}=X_{2t}|X_{1,t-1}=X_{1,t-1},X_{2,t-1}=X_{2,t-1})=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}\hfill \\ \hfill & \left(h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})f(X_{1t}-k_1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi )\right)\hfill \end{array}$

with ${\lambda }_{1t}=b_{11}{X}_{1,t-1}+b_{12}{X}_{2,t-1}+c_1$, ${\lambda }_{2t}=b_{21}{X}_{1,t-1}+b_{22}{X}_{2,t-1}+c_2$, $g_1=min(X_{1t},X_{1,t-1})$, $g_2=min(X_{2t},X_{2,t-1})$, $f\left(\right)$ and $h\left(\right)$ are given in (2) and (8), respectively.

By using (3)–(5), and (9) and (10), we can derive the score equation:

(25)$\begin{array}{c}\hfill {\displaystyle \frac{\partial \ell \left({\mathit{\theta }}_0\right)}{\partial \mathit{\theta }}}=\sum _{t=1}^T{\displaystyle \frac{1}{P_{{\mathit{\theta }}_0}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}}{\displaystyle \frac{\partial P_{{\mathit{\theta }}_0}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial \mathit{\theta }}}=\mathbf{0},\end{array}$

where

$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial {\alpha }_{11}}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}{X}_{1,t-1}h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})f(X_{1t}-k_1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi )\hfill \\ \hfill & \times \left(h(k_1-1,X_{1,t-1}-1,X_{2,t-1},{\alpha }_{11},{\alpha }_{12})-h(k_1,X_{1,t-1}-1,X_{2,t-1},{\alpha }_{11},{\alpha }_{12})\right),\hfill \\ \hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial {\alpha }_{12}}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}{X}_{2,t-1}h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})f(X_{1t}-k_1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi )\hfill \\ \hfill & \times \left(h(k_1-1,X_{1,t-1},X_{2,t-1}-1,{\alpha }_{11},{\alpha }_{12})-h(k_1,X_{1,t-1},X_{2,t-1}-1,{\alpha }_{11},{\alpha }_{12})\right),\hfill \\ \hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial {\alpha }_{21}}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}{X}_{1,t-1}h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})f(X_{1t}-k_1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi )\hfill \\ \hfill & \times \left(h(k_2-1,X_{1,t-1}-1,X_{2,t-1},{\alpha }_{21},{\alpha }_{22})-h(k_2,X_{1,t-1}-1,X_{2,t-1},{\alpha }_{21},{\alpha }_{22})\right),\hfill \\ \hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial {\alpha }_{22}}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}{X}_{2,t-1}h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})f(X_{1t}-k_1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi )\hfill \\ \hfill & \times \left(h(k_2-1,X_{1,t-1},X_{2,t-1}-1,{\alpha }_{21},{\alpha }_{22})-h(k_2,X_{1,t-1},X_{2,t-1}-1,{\alpha }_{21},{\alpha }_{22})\right),\hfill \\ \hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial b_{11}}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})\hfill \\ \hfill & \times X_{1,t-1}f(X_{1t}-k_1-1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi ),\hfill \\ \hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial b_{12}}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})\hfill \\ \hfill & \times X_{2,t-1}f(X_{1t}-k_1-1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi ),\hfill \\ \hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial c_1}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})\hfill \\ \hfill & \times f(X_{1t}-k_1-1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi ),\hfill \\ \hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial b_{21}}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})\hfill \\ \hfill & \times X_{1,t-1}f(X_{1t}-k_1,X_{2t}-k_2-1,{\lambda }_{1t},{\lambda }_{2t},\varphi ),\hfill \\ \hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial b_{22}}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})\hfill \\ \hfill & \times X_{2,t-1}f(X_{1t}-k_1-1,X_{2t}-k_2-1,{\lambda }_{1t},{\lambda }_{2t},\varphi ),\hfill \\ \hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial c_2}}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})\hfill \\ \hfill & \times f(X_{1t}-k_1,X_{2t}-k_2-1,{\lambda }_{1t},{\lambda }_{2t},\varphi ),\hfill \end{array}$

$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial P_{\mathit{\theta }}\left(\left.{\mathit{X}}_t\right|{\mathit{X}}_{t-1}\right)}{\partial \varphi }}=\sum _{k_1=0}^{g_1}\sum _{k_2=0}^{g_2}h(k_1,X_{1,t-1},X_{2,t-1},{\alpha }_{11},{\alpha }_{12})h(k_2,X_{1,t-1},X_{2,t-1},{\alpha }_{21},{\alpha }_{22})\hfill \\ \hfill & \times (f(X_{1t}-k_1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi )-f(X_{1t}-k_1-1,X_{2t}-k_2,{\lambda }_{1t},{\lambda }_{2t},\varphi ).\hfill \\ \hfill & -f(X_{1t}-k_1,X_{2t}-k_2-1,{\lambda }_{1t},{\lambda }_{2t},\varphi )+f(X_{1t}-k_1-1,X_{2t}-k_2-1,{\lambda }_{1t},{\lambda }_{2t},\varphi )).\hfill \end{array}$

The maximizer ${\widehat{\mathit{\theta }}}_T$ of (24) is the CML estimate of ${\mathit{\theta }}_0$, which is obtained by numerically maximizing the log-likelihood (24) or by solving the score Equation (25). To study the asymptotic behaviour of the estimator, we make the following two Assumptions about the parameter space and the underlying process.

Assumption 1.

The parametric space Θ is compact with $\mathsf{\Theta }=\{\mathit{\theta },\mathit{\theta }={\{{\alpha }_{ij},b_{ij},c_1,c_2,\varphi \}}^{\top },i,j=1,2\}$, where $\underline{\delta }\le {\alpha }_{ij},b_{ij}\le \overline{\delta },\underline{c}\le c_i\le \overline{c},\underline{\varphi }\le \varphi \le \overline{\varphi },\gamma =max({\alpha }_{ij}+b_{ij})<1,$ $\underline{\delta }$, $\overline{\delta }$, $\underline{c}$, $\overline{c}$, $\underline{\varphi }$ and $\overline{\varphi }$ are finite positive constants, and ${\mathit{\theta }}_0$ is an interior point in Θ.

Assumption 2.

If there exists a $t\ge 1$, such that ${\mathit{X}}_t\left({\mathit{\theta }}_0\right)={\mathit{X}}_t\left(\mathit{\theta }\right)$, $P_{{\mathit{\theta }}_0}$ a.s., then $\mathit{\theta }={\mathit{\theta }}_0,$ where $P_{{\mathit{\theta }}_0}$ is the probability measure under the true parameter ${\mathit{\theta }}_0.$

To derive the identification of the EBINAR(1) model, we give the following two Lemmas.

Lemma 1.

Let $g_1(x,y,b_{11},b_{12},c_1)=b_{11}x+b_{12}y+c_1$, $b_{11},b_{12},c_1>0$ for $(x,y)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}.$ Then, if $g_1(x,y,b_{11},b_{12},c_1)=g_1(x,y,b_{11}^0,b_{12}^0,c_1^0)$, then $b_{11}=b_{11}^0,b_{12}=b_{12}^0,c_1=c_1^0$.

Proof. 

By the assumption:

$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial g_1(x,y,b_{11},b_{12},c_1)}{\partial x}}={\displaystyle \frac{\partial g_1(x,y,b_{11}^0,b_{12}^0,c_1^0)}{\partial x}},\hfill \\ \hfill & {\displaystyle \frac{\partial g_1(x,y,b_{11},b_{12},c_1)}{\partial y}}={\displaystyle \frac{\partial g_1(x,y,b_{11}^0,b_{12}^0,c_1^0)}{\partial y}},\hfill \\ \hfill & g_1(0,0,b_{11},b_{12},c_1)=g_1(0,0,b_{11}^0,b_{12}^0,c_1^0).\hfill \end{array}$

we obtain: $b_{11}=b_{11}^0,b_{12}=b_{12}^0,c_1=c_1^0$. □

Similarly, we denote $g_2(x,y,b_{21},b_{22},c_2)=b_{21}x+b_{22}y+c_2$. If $g_2(x,y,b_{21},b_{22},c_2)=g_2(x,y,b_{21}^0,b_{22}^0,c_2^0)$, then we have $b_{21}=b_{21}^0,b_{22}=b_{22}^0,c_2=c_2^0$ by the same method.

Lemma 2.

If $\left\{{\mathit{X}}_t\right\}$ is the strictly stationary and ergodic solution of model (6), Assumptions 1 and 2 hold, then model (6) is identifiable.

Proof. 

According to Lemma 1, we conclude that if ${\lambda }_{1t}(b_{11},b_{12},c_1)={\lambda }_{1t}(b_{11}^0,b_{12}^0,c_1^0)$, then $b_{11}=b_{11}^0,b_{12}=b_{12}^0,c_1=c_1^0$. Similarly, if ${\lambda }_{2t}(b_{21},b_{22},c_2)={\lambda }_{2t}(b_{21}^0,b_{22}^0,c_2^0)$, then $b_{21}=b_{21}^0,b_{22}=b_{22}^0,c_2=c_2^0$. Thus, if ${\mathit{ϵ}}_t(b_{ij},c_i,\varphi )={\mathit{ϵ}}_t(b_{ij}^0,c_i^0,{\varphi }^0)$, $t\ge 1,i,j=1,2$, then $b_{ij}=b_{ij}^0$, $c_i=c_i^0$, $\varphi ={\varphi }^0.$ According to (7), we have ${ϵ}_{it}=X_{it}-{\alpha }_{i1}\circ X_{1,t-1}-{\alpha }_{i2}\circ X_{2,t-1}$, $i=1,2.$ If ${ϵ}_{it}\left(\mathit{\theta }\right)={ϵ}_{it}\left({\mathit{\theta }}_0\right)$, then we have ${\alpha }_{i1}={\alpha }_{i1}^0,{\alpha }_{i2}={\alpha }_{i2}^0$, otherwise

$0=E\left({ϵ}_{it}\left(\mathit{\theta }\right)\right)-E\left({ϵ}_{it}\left({\mathit{\theta }}_0\right)\right)=({\alpha }_{i1}-{\alpha }_{i1}^0)E\left(X_{1t}\right)+({\alpha }_{i2}-{\alpha }_{i2}^0)E\left(X_{2t}\right),$

then $E\left(X_{1t}\right)=0$ and $E\left(X_{2t}\right)=0$, which contradicts the fact that $E\left(X_{it}\right)>0$, $i=1,2$.

By Assumption 2, for given $X_{1,t-1}$ and $X_{2,t-1}$, we have

$\varphi =\mathrm{Cov}(X_{1t}\left(\mathit{\theta }\right),X_{2t}\left(\mathit{\theta }\right))=\mathrm{Cov}(X_{1t}\left({\mathit{\theta }}_0\right),X_{2t}\left({\mathit{\theta }}_0\right))={\varphi }^0.$

Thus, $\varphi ={\varphi }^0.$ Hence, model (6) is identifiable. □

Theorem 7.

Suppose that $\left\{{\mathit{X}}_t\right\}$ is the strictly stationary and ergodic solution of model (6) and Assumptions 1 and 2 hold. As $T\to \infty$, there exists an estimator ${\widehat{\mathit{\theta }}}_T$ such that ${\widehat{\mathit{\theta }}}_T\stackrel{a.s.}{⟶}{\mathit{\theta }}_0$.

Proof. 

To prove the strong consistence of ${\widehat{\mathit{\theta }}}_T$, we need to check all the assumptions given in Theorems 4.1.2 and 4.1.3 in Amemiya [25]. Let $W_t\left(\mathit{\theta }\right)=ln{P}_{\mathit{\theta }}\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})$, then $\ell \left(\mathit{\theta }\right)={\sum }_{t=1}^T{W}_t\left(\mathit{\theta }\right)$. We observe that $W_t\left(\mathit{\theta }\right)$ is a measurable function of ${\mathit{X}}_t$ for all $\mathit{\theta }\in \mathsf{\Theta }$, and is continuous in an open and convex neighborhood $N\left({\mathit{\theta }}_0\right)$ of ${\mathit{\theta }}_0$, then there at least exists a point $\overline{\mathit{\theta }}\in N\left({\mathit{\theta }}_0\right)$ such that $W_t\left(\mathit{\theta }\right)$ attains the maximum value at $\overline{\mathit{\theta }}$.

Thus,

$E\left(\underset{\mathit{\theta }\in N\left({\mathit{\theta }}_0\right)}{sup}{W}_t\left(\mathit{\theta }\right)\right)=E\left(ln{P}_{\overline{\mathit{\theta }}}\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})\right)\le lnE\left(P_{\overline{\mathit{\theta }}}\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})\right)<\infty .$

Note that $\left\{{\mathit{X}}_t\right\}$ is a stationary and ergodic time series, and in terms of Theorem 4.1.2 in Amemiya [25], $\frac{1}{T}{\sum }_{t=1}^T{W}_t\left(\mathit{\theta }\right)\to E{W}_t\left(\mathit{\theta }\right)$ in probability as $T\to \infty$. By Jensen’s inequality, we have:

(26)$E\left(W_t\left(\mathit{\theta }\right)\right)-E\left(W_t\left({\mathit{\theta }}_0\right)\right)=Eln\frac{P_{\mathit{\theta }}\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})}{P_{{\mathit{\theta }}_0}\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})}\le lnE\frac{P_{\mathit{\theta }}\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})}{P_{{\mathit{\theta }}_0}\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})}=0.$

Thus, $E{W}_t\left(\mathit{\theta }\right)$ attains a strict local maximum at ${\mathit{\theta }}_0$ by (26) and Lemma 2. Hence, the conditions of Theorem 4.1.2 in Amemiya [25] are fulfilled; thus, there exists an estimator ${\widehat{\mathit{\theta }}}_T$ such that ${\widehat{\mathit{\theta }}}_T\to {\mathit{\theta }}_0$, $T\to \infty$. □

Theorem 8.

If the conditions of Theorem 7 hold, as $T\to \infty$,

$\sqrt{T}({\widehat{\mathit{\theta }}}_T-{\mathit{\theta }}_0)\stackrel{d}{⟶}\mathcal{N}\left(\mathbf{0},{\left(\mathit{J}\left({\mathit{\theta }}_0\right)\right)}^{-1}\mathit{I}\left({\mathit{\theta }}_0\right){\left(\mathit{J}\left({\mathit{\theta }}_0\right)\right)}^{-1}\right),$

where $\mathit{I}\left({\mathit{\theta }}_0\right)=\underset{T\to \infty }{lim}{T}^{-1}E\left({\displaystyle \frac{\partial \ell \left({\mathit{\theta }}_0\right)}{\partial \mathit{\theta }}}{\left({\displaystyle \frac{\partial \ell \left({\mathit{\theta }}_0\right)}{\partial \mathit{\theta }}}\right)}^{\top }\right),\mathit{J}\left({\mathit{\theta }}_0\right)=\underset{T\to \infty }{lim}{T}^{-1}E\left({\displaystyle \frac{{\partial }^2\ell \left({\mathit{\theta }}_0\right)}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\top }}}\right).$

Proof. 

To prove the asymptotic normality of ${\widehat{\mathit{\theta }}}_T$, we need to verify the assumptions of Theorem 4.1.3 in Amemiya [25].

First, by Proposition 1 in Freeland and McCabe [26], it is easy to obtain all the partial derivatives in a similar way, i.e., ${\displaystyle \frac{\partial W_t\left(\mathit{\theta }\right)}{\partial {\theta }_i}}$ exist and are three times continuous differentiable in $\mathsf{\Theta }$; thus, ${\displaystyle \frac{{\partial }^2{W}_t\left(\mathit{\theta }\right)}{\partial {\theta }_i\partial {\theta }_j}}$ exists and is continuous in $N\left({\mathit{\theta }}_0\right)$, for any $i,j,k=1,2,\cdots ,11$. Thus, there at least exists a point $\tilde{\mathit{\theta }}\in N\left({\mathit{\theta }}_0\right)$ such that ${\displaystyle \frac{{\partial }^2{W}_t\left(\mathit{\theta }\right)}{\partial {\theta }_i\partial {\theta }_j}}$ attains the maximum value at $\tilde{\mathit{\theta }}$. Hence,

$E\left(\underset{\mathit{\theta }\in N\left({\mathit{\theta }}_0\right)}{sup}{\displaystyle \frac{{\partial }^2{W}_t\left(\mathit{\theta }\right)}{\partial {\theta }_i\partial {\theta }_j}}\right)=E\left({\displaystyle \frac{{\partial }^2{W}_t\left(\tilde{\mathit{\theta }}\right)}{\partial {\theta }_i\partial {\theta }_j}}\right)<\infty .$

For convenience, we denote: ${\displaystyle \frac{{\partial }^2\ell \left(\mathit{\theta }\right)}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\top }}}=\mathit{G}({\mathit{X}}_t,\mathit{\theta })=\left(g_{ij}({\mathit{X}}_t,\mathit{\theta })\right)$ and $E\left({\displaystyle \frac{{\partial }^2\ell \left(\mathit{\theta }\right)}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\top }}}\right)=\mathit{G}\left(\mathit{\theta }\right)=\left(g_{ij}({\mathit{X}}_t,\mathit{\theta })\right)$. We only need to prove $g_{ij}({\mathit{X}}_t,\mathit{\theta })$ converges to a finite and non-stochastic function $g_{ij}\left(\mathit{\theta }\right)=E\left(g_{ij}({\mathit{X}}_t,\mathit{\theta })\right)$. Let $h\left({\mathit{X}}_t,\mathit{\theta }\right)=g_{ij}({\mathit{X}}_t,\mathit{\theta })-E\left[g_{ij}({\mathit{X}}_t,\mathit{\theta })\right]$, then $Eh\left({\mathit{X}}_t,\mathit{\theta }\right)=0$. Hence, the conditions of Theorem 4.1.3 in [25] are fulfilled. Thus, $T^{-1}{\sum }_{t=1}^{T}h\left({\mathit{X}}_t,\mathit{\theta }\right)$ converges to 0 in probability uniformly in $\mathit{\theta }\in N\left({\mathit{\theta }}_0\right)$. Furthermore, $T^{-1}{\sum }_{t=1}^{T}h\left({\mathit{X}}_t,{\mathit{\theta }}_T^{\star }\right)$ converges to 0 in probability, when ${\mathit{\theta }}_T^{\star }\to {\mathit{\theta }}_0$, $T\to \infty$.

Second, it is easy to see $\mathrm{Cov}\left(\partial W_t\left({\mathit{\theta }}_0\right)/\partial \mathit{\theta }\right)=E\left[\partial W_t\left({\mathit{\theta }}_0\right)/\partial \mathit{\theta }\partial W_t\left({\mathit{\theta }}_0\right)/\partial {\mathit{\theta }}^{\top }\right]$ because $E\left(\partial W_t\left({\mathit{\theta }}_0\right)/\partial \mathit{\theta }\right)=\mathbf{0}$.

Using the ergodic theorem,

${\displaystyle \frac{1}{T}}{\displaystyle \frac{\partial \ell \left({\mathit{\theta }}_0\right)}{\partial \mathit{\theta }}}\stackrel{p}{⟶}E{\displaystyle \frac{1}{P_{{\mathit{\theta }}_0}\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})}}{\displaystyle \frac{\partial P_{{\mathit{\theta }}_0}\left({\mathit{X}}_t\right|{\mathit{X}}_{t-1})}{\partial \mathit{\theta }}}.$

Using the martingale central limit theorem and the Cramér–Wold device, it is direct to show that

$\frac{1}{\sqrt{T}}\partial \ell \left({\mathit{\theta }}_0\right)/\partial \mathit{\theta }\stackrel{d}{⟶}N(\mathbf{0},\mathit{I}\left({\mathit{\theta }}_0\right))\mathrm{with}\mathit{I}\left({\mathit{\theta }}_0\right)={lim}_{T\to \infty }{T}^{-1}E\left[\partial \ell \left({\mathit{\theta }}_0\right)/\partial \mathit{\theta }{\left(\partial \ell \left({\mathit{\theta }}_0\right)/\partial \mathit{\theta }\right)}^{\top }\right].$

Third, there exists an $H(X_{1t},X_{2t})$ such that $\left|{\displaystyle \frac{{\partial }^{3}ln\ell \left(\mathit{\theta }\right)}{\partial {\theta }_i\partial {\theta }_j\partial {\theta }_k}}\right|\le H(X_{1t},X_{2t})$ and $E\left[H(X_{1t},X_{2t})\right]$$<\infty$ by Theorem 4. By the Taylor expansion, we have

(27)$\frac{\partial \ell \left({\widehat{\mathit{\theta }}}_T\right)}{\partial \mathit{\theta }}=\frac{\partial \ell \left({\mathit{\theta }}_0\right)}{\partial \mathit{\theta }}+\frac{{\partial }^2\ell ({\mathit{\theta }}_T^{\star })}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\top }}({\widehat{\mathit{\theta }}}_T-{\mathit{\theta }}_0),$

where ${\mathit{\theta }}_T^{\star }$ lies in between ${\widehat{\mathit{\theta }}}_T$ and ${\mathit{\theta }}_0$. We observe that the ${\displaystyle \frac{\partial \ell \left({\mathit{\theta }}_0\right)}{\partial \mathit{\theta }}}=\mathbf{0}$ in (27) by (25), then

(28)$\sqrt{T}({\widehat{\mathit{\theta }}}_T-{\mathit{\theta }}_0)={\left[\frac{1}{T}{\displaystyle \frac{{\partial }^2\ell \left({\mathit{\theta }}_T^{\star }\right)}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\top }}}\right]}^{-1}\frac{1}{\sqrt{T}}{\displaystyle \frac{\partial \ell \left({\widehat{\mathit{\theta }}}_T\right)}{\partial \mathit{\theta }}}.$

Hence, the asymptotic normality of ${\widehat{\mathit{\theta }}}_T$ follows from (28). □

4. Simulation

In this section, we conduct a simulation study to illustrate the finite sample property of the CML estimate. The simulation is carried out in $\mathtt{R}$ by using the $\mathtt{optim}$ function for the optimization of the conditional log-likelihood function.

In the simulation, we generate data from the non-diagonal EBINAR(1) model and the diagonal EBINAR(1) model. The sizes of samples are chosen to be 50, 100, 200, 500 and 1000 to reflect relatively small, small, moderate, large and relatively large sample sizes, and we use 500 replications. For the simulated sample, performances of the estimators are evaluated by mean squared error (MSE) and mean absolute deviation error (MADE), where $\mathrm{MSE}=\frac{1}{m}{\sum }_{i=1}^m{({\widehat{\phi }}_i-\phi )}^2$, $\mathrm{MADE}=\frac{1}{m}{\sum }_{i=1}^m\left|{\widehat{\phi }}_i-\phi \right|$, where ${\widehat{\phi }}_i$ is the estimator of $\phi$ in the ith replication and m denotes replication times. The used parameter combinations of $\mathit{\theta }={({\alpha }_{11},{\alpha }_{12},{\alpha }_{21},{\alpha }_{22},b_{11},b_{12},b_{21},b_{22},c_1,c_2,\varphi )}^{\top }$ are listed as follows:

(1). For a non-diagonal model: $\mathit{\theta }={(0.3,0.1,0.1,0.1,0.2,0.1,0.1,0.3,0.6,0.6,0.5)}^{\top }$;

(2). For a diagonal model: $\mathit{\theta }=\mathrm{I}:{(0.2,0,0,0.3,0.3,0,0,0.2,0.6,0.6,0.5)}^{\top }$, $\mathrm{II}:(0.2,0,0,0.3,$${0.3,0,0,0.2,2,2,1)}^{\top },$ $\mathrm{III}:{(0.1,0,0,0.4,0.4,0,0,0.1,0.6,0.6,0.5)}^{\top },$ $\mathrm{IV}:(0.1,0,0,0.4,0.4,0,0,0.1,$${2,2,1)}^{\top }.$

Table 1, Table 2, Table 3, Table 4 and Table 5 show that the MSE and MADE decrease with the increase in T for diagonal and non-diagonal models, which implies that the estimators are consistent.

To illustrate the location and dispersion of the estimates, we present the boxplots of the estimates for the non-diagonal and I of diagonal parameter combinations in Figure 1 and Figure 2; the others are similar.

Figure 1 and Figure 2 illustrate the large sample properties of the estimators on a limited sample size. In general, the estimated medians are apparently closer to the real parameter values with the sample size increases. Regarding dispersion issues, both the interquartile ranges and the overall ranges of the produced values are narrower with the sample size increases.

5. Illustrative Examples

In this section, we apply the proposed model to two crime datasets coming from different number of car beats, which is the unique ID for the observation unit’s geography in Pittsburgh Police Department. The crime data is available online at “The Forecasting Principles” site (http://www.forecastingprinciples.com/index.php/crimedata) in the section about Crime data and download on 23 September 2016.

According to Cohen and Gorr [27], the occurrence of criminal mischief may be accompanied by burglary behavior, so does for the robbery. Hence, the monthly counts of burglary and CMIS (or those of burglary and robbery) may exhibit dependence. In this section, we take the monthly counts of burglary and CMIS in beat 11 and the monthly counts of burglary and robbery in beat 26 as examples.

5.1. Monthly Counts of Burglary and CMIS in Beat 11

In this part, we consider the monthly number of burglary and criminal mischief (CMIS) from January 1990 to December 2001 in the geographic ID = 11. Table 6 gives the statistics of the counts of burglaries and CMIS.

Table 6 shows that both the counts of burglaries and CMIS are over-dispersed because their variances are greater than their means. In contrast, this relationship can also be illustrated by the cross-correlation graph of the samples, which are given in Figure 3.

From Figure 3, the counts of burglaries are weakly dependent with those of CMIS. Their plots of sample path, autocorrelation function (ACF) and partial autocorrelation function (PACF) are given in Figure 4, which show that the analyzed data sets are bivariate integer-valued time series with some characteristics of mutual influence.

To give quantitative results about cross-correlation, we compare our model with the following models:

• Full BINAR-BP with ${\mathit{ϵ}}_t$ following $\mathrm{BP}({\lambda }_1,{\lambda }_2,\varphi )$ [16];

• Full BINAR-NB with ${\mathit{ϵ}}_t$ following bivariate negative binomial distribution with parameters $({\lambda }_1,{\lambda }_2,\beta )$; see [14,16] for detail.

As the goodness-of-fit criteria, we use the Akaike information criterion (AIC), the Bayesian information criterion (BIC) and the mean square error of the Pearson residuals (PRMS), which is equal to ${\sum }_{t=1}^n{Z}_t^2/(n-p^{*})$, where $p^{*}$ denotes the number of estimated parameters and $Z_t$ denotes standardized Pearson residuals.

The CML estimate and approximated standard error (SE) of the parameter, including the fitted values of PRMS, AIC, BIC and log-likelihood function (Log Lik), are summarized in Table 7, where the approximated standard error is computed by using the estimated version of the robust sandwich matrix ${\left(\mathit{J}\left({\mathit{\theta }}_0\right)\right)}^{-1}\mathit{I}\left({\mathit{\theta }}_0\right){\left(\mathit{J}\left({\mathit{\theta }}_0\right)\right)}^{-1}$; see Theorem 8 for details.

Table 7 shows that the EBINAR(1) model takes the highest Log Lik value and the lowest AIC, BIC and PRMS for the monthly number of burglaries and CMIS. Hence, the EBINAR(1) model is more suitable for the data sets.

5.2. Monthly Counts of Burglaries and Robberies in Beat 26

In this part, we consider the monthly number of burglaries and robberies from January 1990 to December 2001 in the geographic ID = 26; see Table 8 for some of their statistics.

Table 8 shows the monthly number of burglary and robbery are over-dispersed. In contrast, this relationship can also be illustrated by their cross-correlation graph given in Figure 5, which shows that the counts of burglary are significantly dependent on those of robbery.

To further illustrate the the monthly number of burglaries and robberies in beat 26, we present their sample path, ACF and PACF plots in Figure 6, from which we can conclude that the analyzed data set exhibits some characteristics of mutual influence.

To give quantitative result about cross-correlation, we compare our model with the Full BINAR-BP and Full BINAR-NB models. The CML estimate and SE, including the fitted PRMS, AIC, BIC and Log Lik, are summarized in Table 9.

Table 9 shows that the EBINAR(1) model takes the highest Log Lik value and the lowest AIC, BIC and PRMS for burglaries and robberies in beat 26. Hence, EBINAR(1) model is more suitable.

To sum up, our findings reveal that there are some connections for the burglary and CMIS in beat 11 and those for the burglary and robbery in beat 26, which agrees with the conclusion of Cohen and Gorr [27]. Of course, the counts of burglary may be affected by other crime activities, such as simple assault, vagrancy and trespassing, which will be studied in a further study.

Remark 2.

For the two real datasets, our EBINAR(1) model performs best, but it is not clear enough regarding predicting unseen data. To further illustrate the better performance of the new model in prediction, one available solution is to conduct a further experiment when dividing the considered data into a training set and test set. In addition such experiment will be considered in future study of the crime data.

6. Concluding Remarks

This paper proposes a more flexible model for bivariate integer-valued time series data, i.e., the EBINAR(1) model, whose innovation vector is time-dependent. It is a generalization of the EINAR(1) model [11] to the two-dimensional case as well as a generalization of the BINAR(1) model [14,16], but with more flexibility. We discuss some necessary properties of the model, the CML estimators of parameters and their large-sample properties. Simulation was conducted to examine the finite sample performance of estimators. Real data examples are provided to illustrate our model to be effective relative to existing models.

To make the bivariate INAR-type models more flexible with respect to real-data applications in some cases, it may be interesting to include explanatory covariates or periodicity in the model to account for dependence through thinning operations on several factors, which will be considered in another project: see Aknouche et al. [28] and Chen and Khamthong [29].

Author Contributions

Conceptualization, H.C. and F.Z.; methodology, H.C. and F.Z.; software, H.C., F.Z. and X.L.; validation, H.C., F.Z. and X.L.; formal analysis, H.C.; investigation, H.C.; resources, F.Z.; data curation, X.L.; writing—original draft preparation, H.C. and F.Z.; writing—review and editing, H.C. and F.Z.; visualization, H.C. and F.Z.; supervision, F.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Crime data are available online at The Forecasting Principles site and were downloaded on 23 September 2016 from http://www.forecastingprinciples.com/index.php/crimedata.

Acknowledgments

We are very grateful to anonymous referees for providing several exceptionally helpful comments which led to a significant improvement in the manuscript.

Conflicts of Interest

The authors declare no conflict of interests in the publication of this paper.

Abbreviations

The following abbreviations are used in this manuscript:

${\mathit{A}}^{-1}$	inverse of matrix $\mathit{A}$;
${\mathit{A}}^{\top }$	transpose of the matrix or vector $\mathit{A}$;
$\parallel ·\parallel$	Euclidean norm of a matrix or vector;
$|·|$	absolute value of a univariate variable;
$\stackrel{d}{⟶}$	convergence in distribution;
$\stackrel{p}{⟶}$	convergence in probability one;
pmf	probability mass function;
CML	conditional maximun likelihood;
AIC	Akaike information criterion;
BIC	Bayesian information criterion;
SE	standard error;
PRMS	mean square error of the Pearson residual;
Para.	parameter.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Boxplots of the CML estimates for non-diagonal EBINAR(1) model.

Figure 2. Boxplots of the CML estimates for diagonal EBINAR(1) model with parameter I.

Figure 3. Cross-correlation between the monthly number of burglaries and CMIS in beat 11.

Figure 4. Beat 11: (1) monthly number of burglary, (2) monthly number of CMIS, (3) ACF of burglary, (4) ACF of CMIS, (5) PACF of burglary, (6) PACF of CMIS.

Figure 5. Cross-correlation between the monthly number of burglaries and robberies in beat 26.

Figure 6. Beat 26: (1) monthly number of burglaries, (2) monthly number of robberies, (3) ACF of burglaries, (4) ACF of robberies, (5) PACF of burglaries, (6) PACF of robberies.

Table 1

Results for non-diagonal EBINAR(1) model.

Size		${\mathit{\alpha }}_{11}$	${\mathit{\alpha }}_{12}$	${\mathit{\alpha }}_{21}$	${\mathit{\alpha }}_{22}$	${\mathit{b}}_{11}$	${\mathit{b}}_{12}$	${\mathit{b}}_{21}$	${\mathit{b}}_{22}$	${\mathit{c}}_1$	${\mathit{c}}_2$	$\mathit{\varphi }$
50	MSE	0.0095	0.0122	0.0082	0.0177	0.0048	0.0187	0.0100	0.0380	0.0172	0.0246	0.0160
	MADE	0.0561	0.0537	0.0390	0.0460	0.0371	0.0487	0.0409	0.0750	0.0857	0.0991	0.0811
100	MSE	0.0069	0.0043	0.0033	0.0127	0.0063	0.0070	0.0099	0.0090	0.0128	0.0210	0.0150
	MADE	0.0473	0.0269	0.0287	0.0417	0.0399	0.0355	0.0378	0.0578	0.0763	0.0913	0.0736
200	MSE	0.0044	0.0019	0.0034	0.0108	0.0058	0.0063	0.0054	0.0063	0.0178	0.0185	0.0170
	MADE	0.0421	0.0281	0.0326	0.0291	0.0426	0.0363	0.0372	0.0532	0.0901	0.0876	0.0824
500	MSE	0.0033	0.0008	0.0011	0.0105	0.0044	0.0027	0.0008	0.0061	0.0029	0.0044	0.0063
	MADE	0.0317	0.0161	0.0213	0.0469	0.0379	0.0293	0.0225	0.0556	0.0446	0.0529	0.0465
1000	MSE	0.0002	0.0001	0.0005	0.0041	0.0014	0.0007	0.0002	0.0006	0.0006	0.0015	0.0048
	MADE	0.0116	0.0079	0.0167	0.0413	0.0280	0.0225	0.0114	0.0190	0.0199	0.0345	0.0379

Table 2

Results for diagonal EBINAR(1) model with parameter I.

Size		${\mathit{\alpha }}_{11}$	${\mathit{\alpha }}_{22}$	${\mathit{b}}_{11}$	${\mathit{b}}_{22}$	${\mathit{c}}_1$	${\mathit{c}}_2$	$\mathit{\varphi }$
50	MSE	0.0031	0.0146	0.0205	0.0030	1.2188	0.6134	0.3256
	MADE	0.0406	0.1021	0.1250	0.0398	0.7710	0.5880	0.5059
100	MSE	0.0020	0.0062	0.0134	0.0023	0.7887	0.4527	0.2778
	MADE	0.0323	0.0665	0.0976	0.0330	0.6125	0.4978	0.4843
200	MSE	0.0015	0.0045	0.0088	0.0010	0.4832	0.3250	0.2572
	MADE	0.0300	0.0524	0.0775	0.0217	0.5016	0.3995	0.4742
500	MSE	0.0007	0.0031	0.0043	0.0010	0.2240	0.1528	0.2198
	MADE	0.0202	0.0396	0.0495	0.0227	0.3655	0.2670	0.4312
1000	MSE	0.0005	0.0020	0.0022	0.0004	0.1965	0.1147	0.1789
	MADE	0.0156	0.0352	0.0377	0.0142	0.3100	0.2084	0.3954

Table 3

Results for diagonal EBINAR(1) model with parameter II.

Size		${\mathit{\alpha }}_{11}$	${\mathit{\alpha }}_{22}$	${\mathit{b}}_{11}$	${\mathit{b}}_{22}$	${\mathit{c}}_1$	${\mathit{c}}_2$	$\mathit{\varphi }$
50	MSE	0.0154	0.0183	0.0122	0.0191	0.7510	1.1007	0.3183
	MADE	0.0871	0.0988	0.0903	0.0949	0.6894	0.8269	0.5008
100	MSE	0.0059	0.0089	0.0059	0.0072	0.4333	0.6728	0.2336
	MADE	0.0470	0.0742	0.0599	0.0582	0.4957	0.5889	0.4442
200	MSE	0.0042	0.0044	0.0041	0.0053	0.2876	0.4939	0.1983
	MADE	0.0411	0.0499	0.0475	0.0470	0.3796	0.4866	0.4193
500	MSE	0.0027	0.0035	0.0036	0.0025	0.1038	0.3240	0.1899
	MADE	0.0344	0.0400	0.0414	0.0326	0.2636	0.4216	0.4107
1000	MSE	0.0013	0.0022	0.0017	0.0011	0.0730	0.0855	0.1352
	MADE	0.0238	0.0303	0.0307	0.0204	0.1978	0.2221	0.3512

Table 4

Results for diagonal EBINAR(1) model with parameter III.

Size		${\mathit{\alpha }}_{11}$	${\mathit{\alpha }}_{22}$	${\mathit{b}}_{11}$	${\mathit{b}}_{22}$	${\mathit{c}}_1$	${\mathit{c}}_2$	$\mathit{\varphi }$
50	MSE	0.0027	0.0048	0.0473	0.0083	0.0078	0.0083	0.0013
	MADE	0.0258	0.0428	0.1533	0.0546	0.0620	0.0586	0.0316
100	MSE	0.0036	0.0064	0.0429	0.0102	0.0089	0.0087	0.0017
	MADE	0.0359	0.0485	0.1486	0.0640	0.0680	0.0638	0.0330
200	MSE	0.0060	0.0059	0.0380	0.0059	0.0054	0.0047	0.0017
	MADE	0.0341	0.0469	0.1239	0.0469	0.0541	0.0507	0.0321
500	MSE	0.0018	0.0042	0.0082	0.0031	0.0046	0.0039	0.0016
	MADE	0.0312	0.0429	0.0638	0.0380	0.0477	0.0426	0.0287
1000	MSE	0.0011	0.0030	0.0037	0.0020	0.0020	0.0016	0.0005
	MADE	0.0253	0.0395	0.0380	0.0331	0.0384	0.0341	0.0182

Table 5

Results for diagonal EBINAR(1) model with parameter IV.

Size		${\mathit{\alpha }}_{11}$	${\mathit{\alpha }}_{22}$	${\mathit{b}}_{11}$	${\mathit{b}}_{22}$	${\mathit{c}}_1$	${\mathit{c}}_2$	$\mathit{\varphi }$
50	MSE	0.0031	0.0146	0.0205	0.0030	1.2188	0.6134	0.3256
	MADE	0.0406	0.1021	0.1250	0.0398	0.7710	0.5880	0.5059
100	MSE	0.0020	0.0062	0.0134	0.0023	0.7887	0.4527	0.2778
	MADE	0.0323	0.0665	0.0976	0.0330	0.6125	0.4978	0.4843
200	MSE	0.0015	0.0045	0.0088	0.0010	0.4832	0.3250	0.2572
	MADE	0.0300	0.0524	0.0775	0.0217	0.5016	0.3995	0.4742
500	MSE	0.0007	0.0031	0.0043	0.0010	0.2240	0.1528	0.2198
	MADE	0.0202	0.0396	0.0495	0.0227	0.3655	0.2670	0.4312
1000	MSE	0.0005	0.0020	0.0022	0.0004	0.1965	0.1147	0.1789
	MADE	0.0156	0.0352	0.0377	0.0142	0.3100	0.2084	0.3954

Table 6

Summary statistics for the monthly number of burglaries and CMIS in beat 11.

Data	Mean	Variance	Minimum	Median	Maximum
Burglary	2.8819	4.1188	0	3	10
CMIS	6.3819	10.0839	1	6	22

Table 7

Estimates for the monthly numbers of burglaries and those of CMIS in beat 11.

	EBINAR(1)			Full BINAR(1)-NB			Full BINAR(1)-BP
Para.	Estimate	SE	Para.	Estimate	SE	Para.	Estimate	SE
${\widehat{\alpha }}_{11}$	0.1689	0.1559	${\widehat{\alpha }}_{11}$	0.2784	0.0665	${\widehat{\alpha }}_{11}$	0.2993	0.0838
${\widehat{\alpha }}_{12}$	0.0179	0.0411	${\widehat{\alpha }}_{12}$	0.0217	0.0092	${\widehat{\alpha }}_{12}$	0.0217	0.0215
${\widehat{\alpha }}_{21}$	0.0390	0.1447	${\widehat{\alpha }}_{21}$	0.1060	0.0550	${\widehat{\alpha }}_{21}$	0.1060	0.0719
${\widehat{\alpha }}_{22}$	0.1131	0.1236	${\widehat{\alpha }}_{22}$	0.5010	0.0295	${\widehat{\alpha }}_{22}$	0.1934	0.0551
${\widehat{b}}_{11}$	0.0690	0.1460
${\widehat{b}}_{12}$	0.0093	0.0559
${\widehat{b}}_{21}$	0.1014	0.1809
${\widehat{b}}_{22}$	0.1354	0.1414
${\widehat{c}}_1$	1.0007	0.4372	${\widehat{\lambda }}_1$	1.9814	0.0186	${\widehat{\lambda }}_1$	1.5164	0.2561
${\widehat{c}}_2$	3.3190	0.5478	${\widehat{\lambda }}_2$	2.3137	0.0166	${\widehat{\lambda }}_2$	4.5258	0.4493
$\widehat{\varphi }$	0.5273	0.2628	$\widehat{\beta }$	0.1374	0.9759	$\widehat{\varphi }$	0.4044	0.2274
PRMS	0.0064			0.0245			0.0103
AIC	1315.4620			1387.8913			1350.9488
BIC	1348.1300			1408.6800			1371.7375
Log Lik	−646.7310			−686.9457			−668.4744

Table 8

Summary statistics for the monthly number of burglaries and robberies in beat 26.

Data	Mean	Variance	Minimum	Median	Maximum
Burglary	3.9306	9.7434	0	3	15
Robbery	3.0625	9.6394	0	2	17

Table 9

Estimates for the monthly number of burglaries and robberies in beat 26.

	EBINAR(1)			Full BINAR(1)-NB			Full BINAR(1)-BP
Para.	Estimate	SE	Para.	Estimate	SE	Para.	Estimate	SE
${\widehat{\alpha }}_{11}$	0.3117	0.0654	${\widehat{\alpha }}_{11}$	0.2314	0.3042	${\widehat{\alpha }}_{11}$	0.2765	0.0537
${\widehat{\alpha }}_{12}$	0.2086	0.0611	${\widehat{\alpha }}_{12}$	0.3172	0.2442	${\widehat{\alpha }}_{12}$	0.0927	0.0471
${\widehat{\alpha }}_{21}$	0.0900	0.0511	${\widehat{\alpha }}_{21}$	0.1099	0.2834	${\widehat{\alpha }}_{21}$	0.0001	0.0000
${\widehat{\alpha }}_{22}$	0.1906	0.1163	${\widehat{\alpha }}_{22}$	0.4361	0.2244	${\widehat{\alpha }}_{22}$	0.4249	0.0415
${\widehat{b}}_{11}$	0.0671	0.0706
${\widehat{b}}_{12}$	0.2280	0.0653
${\widehat{b}}_{21}$	0.1233	0.0511
${\widehat{b}}_{22}$	0.3358	0.1161
${\widehat{c}}_1$	0.2043	0.2048	${\widehat{\lambda }}_1$	2.2310	0.0026	${\widehat{\lambda }}_1$	1.7652	0.2048
${\widehat{c}}_2$	0.4139	0.1139	${\widehat{\lambda }}_2$	1.1708	0.0076	${\widehat{\lambda }}_2$	0.9604	0.1601
$\widehat{\varphi }$	0.5599	0.1187	$\widehat{\beta }$	0.4073	0.7189	$\widehat{\varphi }$	0.7778	0.1494
PRMS	0.0087			0.0748			0.0992
AIC	1320.8092			1344.6968			1357.7718
BIC	1353.4771			1365.4855			1378.5604
Log Lik	−649.4046			−665.3484			−671.8859

References

1. Steutel, F.W.; van Harn, K. Discrete analogues of self-decomposability and stability. Ann. Probab.; 1979; 7, pp. 893-899. [DOI: https://dx.doi.org/10.1214/aop/1176994950]

2. Al-Osh, M.A.; Alzaid, A.A. First-order integer-valued autoregressive process. J. Time Ser. Anal.; 1987; 8, pp. 261-275. [DOI: https://dx.doi.org/10.1111/j.1467-9892.1987.tb00438.x]

3. McKenzie, E. Some simple models for discrete variate time series. Water Resoure Bull.; 1985; 21, pp. 645-650. [DOI: https://dx.doi.org/10.1111/j.1752-1688.1985.tb05379.x]

4. Du, J.; Li, Y. The integer valued autoregressive INAR(p) model. J. Time Ser. Anal.; 1991; 12, pp. 129-142.

5. Alzaid, A.A.; Omair, M.A. Poisson difference integer valued autoressive model of order one. Bull. Malays. Math. Sci. Soc.; 2014; 37, pp. 465-485.

6. Chen, H.; Li, Q.; Zhu, F. Binomial AR(1) processes with innovational outliers. Commun. Stat. Theory Methods; 2021; 50, pp. 446-472. [DOI: https://dx.doi.org/10.1080/03610926.2019.1635704]

7. Weiß, C.H. Thinning operations for modeling time series of counts—A survey. Adv. Stat. Anal.; 2008; 92, pp. 319-341. [DOI: https://dx.doi.org/10.1007/s10182-008-0072-3]

8. Scotto, M.G.; Wei, C.H.; Gouveia, S. Thinning-based models in the analysis of integer-valued time series: A review. Stat. Model.; 2015; 15, pp. 590-618. [DOI: https://dx.doi.org/10.1177/1471082X15584701]

9. Davis, R.A.; Fokianos, K.; Holan, S.H.; Joe, H.; Livsey, J.; Lund, R.; Pipiras, V.; Ravishanker, N. Count time series: A methodological review. J. Am. Stat. Assoc.; 2021; 116, pp. 1533-1547. [DOI: https://dx.doi.org/10.1080/01621459.2021.1904957]

10. Buckley, F.M.; Pollett, P.K. Limit theorems for discrete-time metapopulation models. Probab. Surv.; 2010; 7, pp. 53-83. [DOI: https://dx.doi.org/10.1214/10-PS158]

11. Weiß, C.H. A Poisson INAR(1) model with serially dependent innovations. Metrika; 2015; 78, pp. 829-851. [DOI: https://dx.doi.org/10.1007/s00184-015-0529-9]

12. Franke, J.; Rao, T.S. Multivariate First-Order Integer Valued Autoregressions; Technical Report Department of Mathematics, UMIST: Manchester, UK, 1995.

13. Latour, A. The multivariate GINAR(p) process. Adv. Appl. Probab.; 1997; 29, pp. 228-248. [DOI: https://dx.doi.org/10.2307/1427868]

14. Pedeli, X.; Karlis, D. A bivariate INAR(1) processes with application. Stat. Model.; 2011; 11, pp. 325-349. [DOI: https://dx.doi.org/10.1177/1471082X1001100403]

15. Pedeli, X.; Karlis, D. On estimation of the bivariate Poisson INAR process. Commun. Stat. Simul. Comput.; 2013; 42, pp. 514-533. [DOI: https://dx.doi.org/10.1080/03610918.2011.639001]

16. Pedeli, X.; Karlis, D. Some properties of multivariate INAR(1) processes. Comput. Stat. Data Anal.; 2013; 67, pp. 213-225. [DOI: https://dx.doi.org/10.1016/j.csda.2013.05.019]

17. Ravishanker, N.; Serhiyenko, V.; Willig, M.R. Hierarchical dynamic models for multivariate times series of counts. Stat. Its Interface; 2014; 7, pp. 559-570. [DOI: https://dx.doi.org/10.4310/SII.2014.v7.n4.a11]

18. Popović, P.M. A bivariate INAR(1) model with different thinning parameters. Stat. Pap.; 2016; 57, pp. 517-538. [DOI: https://dx.doi.org/10.1007/s00362-015-0667-1]

19. Scotto, M.G.; Wei, C.H.; Silva, M.E.; Pereira, I. Bivariate binomial autoregressive models. J. Multivar. Anal.; 2014; 125, pp. 233-251. [DOI: https://dx.doi.org/10.1016/j.jmva.2013.12.014]

20. Li, Q.; Chen, H.; Liu, X. A new bivariate random coefficient INAR(1) model with applications. Symmetry; 2022; 14, 39. [DOI: https://dx.doi.org/10.3390/sym14010039]

21. Kocherlakota, S.; Kocherlakota, K. Bivariate Discrete Distributions; Marcel Dekker: New York, NY, USA, 1992; pp. 87-97.

22. Heathcote, C.R. Corrections and comments on the paper “A branching process allowing immigration”. J. R. Stat. Soc. Ser. B; 1966; 28, pp. 213-217. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1966.tb00634.x]

23. Shumway, R.H.; Stoffer, D.S. Time Series Analysis and Its Applications with R examples; 3rd ed. Springer: New York, NY, USA, 2011.

24. Silva, M.E.; Oliveira, V.L. Difference equations for the higher-order moments and cumulants of the INAR(1) model. J. Time Ser. Anal.; 2004; 25, pp. 317-333. [DOI: https://dx.doi.org/10.1111/j.1467-9892.2004.01685.x]

25. Amemiya, T. Advanced Econometrics; Harvard University Press: Cambridge, MA, USA, 1985; pp. 110-112.

26. Freeland, R.K.; McCabe, B.P.M. Analysis of low count time series data by poisson autoregression. J. Time Ser. Anal.; 2004; 25, pp. 701-722. [DOI: https://dx.doi.org/10.1111/j.1467-9892.2004.01885.x]

27. Cohen, J.; Gorr, W.L. Development of Crime Forecasting and Mapping Systems for Use by Police; Inter-University Consortium for Political and Social Research: New York, NY, USA, 2005; [DOI: https://dx.doi.org/10.3886/ICPSR04545.v1]

28. Aknouche, A.; Bentarzi, W.; Demouche, N. On periodic ergodicity of a general periodic mixed Poisson autoregression. Stat. Probab. Lett.; 2018; 134, pp. 15-21. [DOI: https://dx.doi.org/10.1016/j.spl.2017.10.014]

29. Chen, C.W.S.; Khamthong, K. Bayesian modelling of nonlinear negative binomial integer-valued GARCHX models. Stat. Model.; 2020; 20, pp. 537-561. [DOI: https://dx.doi.org/10.1177/1471082X19845541]