2cm2.5cm2.5cm2cm

# Multivariate Count Time Series Modelling

Konstantinos Fokianos

###### Abstract

We review autoregressive models for the analysis of multivariate count time
series. In doing so, we discuss the choice of a suitable distribution for a
vectors of count random variables. This review focus on three main approaches
taken for multivariate count time series analysis: (a) integer autoregressive
processes, (b) parameter-driven models and (c) observation-driven models. The
aim of this work is to highlight some recent methodological developments and
propose some potentially useful research topics.

Keywords: auto-correlation, covariates, copula, estimation, multivariate count
distribution, prediction

##  1 Introduction

This work reviews three main approaches that have been put forward for
analysis and inference of multivariate count time series. By now there is an
extensive literature for modeling univariate count time series, see the recent
volume by Davis et al. (2016) and the review article by Davis et al. (2021),
for example. Theoretical and methodological development for multivariate count
time series is still on-going research area; see Paul et al. (2008) for a
medical application, Pedeli and Karlis (2013b) for a financial study and more
recently Ravishanker et al. (2015), for a marketing application, and Livsey et
al. (2018) for an environmental study. The interested reader is referred to
the review paper by Karlis (2016), for additional literature. The aim of this
work is to highlight some recent methodological developments and propose some
potentially useful research topics.

Following conventional theory, the standard venue for developing multivariate
count time series models requires specification of a joint conditional
distribution. Then, likelihood inference, for a given autoregressive model,
provides estimation, testing and all type of standard output. However,
choosing a joint count distribution is a challenging problem. There are
numerous proposals available in the literature generalizing univariate Poisson
probability mass function (p.m.f); some of these are reviewed in Sec. 2. The
main obstacle is that the p.m.f of a multivariate ”Poisson” discrete random
vector is often of complicated functional form and therefore maximum
likelihood inference is theoretically and numerically burdensome. The choice
of joint distribution for modeling multivariate count data is an interesting
topic and some comments will be made throughout this work.

Three main approaches are reviewed for modeling multivariate count time
series. The first approach is based on the theory of integer autoregressive
(INAR) models and was initiated by Franke and Rao (1995) and Latour (1997). It
was applied more recently by Pedeli and Karlis (2013a, b) and Scotto et al.
(2014) and Darolles et al. (2019). Estimation for INAR models is based on
least squares methodology and/or likelihood based methods. But, even for
univariate INAR models, likelihood theory is quite cumbersome, especially when
dealing with higher order autoregressive models. This methodology is reviewed
in Sec. 3.

The second model class reviewed is that of parameter driven models whose
dynamics–according to the broad categorization introduced by Cox (1981)– are
driven by an unobserved process. Such models are also called state space
models and have found numerous applications; see Zeger (1988), Harvey and
Fernandes (1989), Fahrmeir and Tutz (2001), West and Harrison (1997), Durbin
and Koopman (2000) and Früwirth-Schnatter and Wagner (2006), among others, for
contributions on univariate time series modeling. Multivariate state space
models were studied by Jørgensen et al. (1996) and Jung et al. (2011); see
also Ravishanker et al. (2014); Ravishanker et al. (2015), among others, for
more recent contributions. We review these models and we illustrate that, even
though their specification is simple, they still require extensive
computational efforts to be applied (using either frequentist or Bayesian
methods). Developments in this area are reviewed in Sec. 4.

Section 5 goes over observation-driven process. This is the third class of
models included in this work. Their main characteristic is that dynamics
evolve according to past values of the process plus some noise. For example,
ordinary autoregressive models belong to this class. Univariate
observation–driven models for count time series have been studied by Zeger and
Qaqish (1988), Fahrmeir and Tutz (2001), Rydberg and Shephard (2000), Kedem
and Fokianos (2002), Fokianos et al. (2009), Fokianos and Tjøstheim (2011),
Davis and Liu (2016), Ahmad and Franq (2016), Douc et al. (2017), among
others. There is a growing literature within the framework of multivariate
observation-driven count time series models; see Heinen and Rengifo (2007),
Liu (2012), Andreassen (2013), Ahmad (2016), Lee et al. (2018), Gouriéroux and
Lu (2019), Fokianos et al. (2020), Opschoor et al. (2020), Piancastelli et al.
(2020) for instance. Most of these studies are concerned with linear and/or
log-linear count time series models but other alternatives can be developed.

From a personal perspective, I think that this research area is still
underdeveloped and there is ample space for exciting new developments. Some
more recent works include that of Veraart (2019), who studies continuous-time
models for multivariate count time series whose marginal distribution is
infinitely divisible. This construction allows for separate modeling of serial
correlation and the cross-sectional dependence. Additional work by Hall et al.
(2019) considers high-dimensional count time series and studies the issues of
inference for autoregressive parameters and the corresponding network
structure by developing a sparsity-regularized maximum likelihood estimator.
Finally, works by Zhang and Wu (2017) and Dahlhaus (1997) are potentially
applicable to develop models for high-dimensional and non-stationary data.

The list of references focus heavily on works related to this research
direction and it is by no means complete. However, further information is
given by the references and interested readers should consult them for further
details.

##  2 Choosing a Multivariate Count Distributions

We outline some parametric multivariate count distributions for independent
data that have been found useful for regression analysis. Simple properties of
those models are discussed and their connection to time series data is
illustrated in Sec. 5. The goal is to show that some basic multivariate
distributions are directly applicable for fitting regression models and
develop inference. There are several alternative venues, for instance we can
rely on copulas and mixture models; see Johnson et al. (1997, Ch. 37), Joe
(1997, Sec 7.2) for numerous multivariate count distributions and Cameron and
Trivedi (2013, Ch.8) for an in-depth review of multivariate count regression
models. As a general remark, joint p. m. f. of a discrete random vector often
has complex functional form which is far from being useful to develop
likelihood based inference. Further recent work and good summaries of up-to
date parametric models can be found in Zhang et al. (2017), Inouye et al.
(2017), Koochemeshkian et al. (2020), among others.

A multinomial distribution, is traditionally employed for analysis of
multivariate count data by employing a multinomial logistic regression models.
The multinomial distribution and the case of Dirichlet-multinomial
distribution are not included in this work (see Zhang et al. (2017) for more).
Those distributions are defined by a conditioning argument which might not
extend to time series context. It is an elementary exercise to show that if
$Y_{i}$ is independent Poisson distributed with mean $\lambda_{i}$ for
$i=1,2,\ldots,d$, then the conditional distribution of $(Y_{1},\ldots Y_{d})$
given $\sum_{i}Y_{i}=y_{.}$ is multinomial with parameters $y_{.}$ and
$\lambda_{i}/\sum_{j}\lambda_{j}$, for $i=1,2,\ldots,d$. So, the multinomial
distribution applies to multivariate count modeling but subject to the
restriction that it is supported on $\\{y=(y_{1},\ldots
y_{d})^{T}:\sum_{i}y_{i}=y_{.}\\}$. This support constraint should be
considered cautiously in the context of dependent data. Additionally, this
approach provides a conditional likelihood function for estimating regression
parameters. The full likelihood function requires knowledge of the p.m.f. of
$\sum_{i}Y_{i}$. Similar remark holds for Dirichlet-multinomial regression.

In what follows, we denote by $Y=(Y_{1},\ldots,Y_{d})^{T}$ a $d$-dimensional
vector of counts whose components are not necessarily independent.

###  2.1 Multivariate Poisson Distributions

This class of distributions (Kocherlakota and Kocherlakota (1992) and Johnson
et al. (1997)), generalizes the univariate Poisson models. Put
$Y_{i}=W_{i}+W$, $i=1,2,\ldots,d$ where $W_{i}\sim\mbox{Poisson}(\lambda_{i})$
and $W\sim\mbox{Poisson}(\lambda_{0})$ and all $W$’s are independent. Then,
the joint p.m.f of $Y$ is given

$\displaystyle\mbox{P}[Y=y]=\exp\bigl{(}-\sum_{i=0}^{d}\lambda_{i}\bigr{)}\biggl{(}\prod_{i=1}^{d}\frac{\lambda_{i}^{y_{i}}}{y_{i}!}\biggr{)}\sum_{k=0}^{\min_{i}y_{i}}\biggl{(}\prod_{i=1}^{d}{y_{i}\choose
k}\biggr{)}k!\biggr{(}\frac{\lambda_{0}}{\prod_{i=1}^{d}\lambda_{i}}\biggl{)}^{k}.$
(1)

The marginals, $Y_{i}$, are Poisson with mean $\lambda_{i}+\lambda_{0}$, for
$i=1,2,\dots,d$ and it holds that $\mbox{Cov}(Y_{i},Y_{j})=\lambda_{0}$ which
is always positive. In addition, the parameter $\lambda_{0}$ determines all
possible pairwise correlations so the resulting model is of limited use. These
facts and the complicated form of (1) make this model suitable for relatively
low dimensional analysis where the marginals $Y_{i}$ are positively
correlated. An E-M type algorithm has been proposed by Karlis (2003) for
inference but this approach is still hard to implement when the dimension $d$
is large. Given a vector of regressors, say $X$, an appropriate regression
model $\lambda_{i}$ as a function of $X$ (in terms of linear of log-linear
link function) and $\lambda_{0}$ is taken as constant . Then using (1)
likelihood inference is straightforward, at least for low dimensions, for more
see Cameron and Trivedi (2013, Ch. 8.4.1).

###  2.2 Mixed Poisson Models

Mixed models provide a general class of multivariate count distributions, see
Marshall and Olkin (1988). Assume now that $Y_{i}$ are independent Poisson
distributed with $\mbox{E}[Y_{i}]=\lambda_{i}$, conditional on $\lambda_{i}$,
$i=1,2,\ldots,d$. Suppose that the vector
$\lambda=(\lambda_{1},\ldots,\lambda_{d})^{T}$ is distributed according to
some distribution $G(\lambda)$. The mixed Poisson distribution is defined as

$\displaystyle\mbox{P}[Y=y]=\int_{(R^{+})^{d}}\Bigl{[}\prod_{i=1}^{n}\frac{\exp(-\lambda_{i})\lambda_{i}^{y_{i}}}{y_{i}!}\Bigr{]}dG(\lambda).$
(2)

Several choices for the mixing distribution $G(\cdot)$ exist. But it is always
true (provided that appropriate moments exist) that

$\begin{split}\mbox{E}[Y]&=\mbox{E}[\lambda],\\\
\mbox{Var}[Y]&=\mbox{diag}(\mbox{E}(\lambda))+\mbox{Var}[\lambda]=\mbox{diag}(\mbox{E}(Y))+\mbox{Var}[\lambda],\end{split}$
(3)

because of the first equality and $\mbox{diag}(x)$ denotes a diagonal matrix
whose elements are given by a vector $x$. A trivial example of (2) is given
when $G(.)$ is the Dirac distribution placing its mass at $\lambda$. Then, $Y$
is just a vector which consists of independent Poisson random variables.
Finite mixtures of multivariate Poisson distributions, with application to
clustering, have been discussed by Karlis and Meligkotsidou (2007). Another
interesting case is when $\lambda$ follows the $d$-dimensional log-normal
distribution with parameters $\mu=(\mu_{1},\ldots,\mu_{d})^{T}$ and
$\Sigma=(\sigma_{ij})_{i,j=1,\ldots,d}$. Though no closed formula exists for
the p.m.f. of $Y$, eq. (3) implies that for $i=1,2,\ldots,d$

$\displaystyle\mbox{E}[Y_{i}]$ $\displaystyle=$
$\displaystyle\exp(\mu_{i}+0.5\sigma_{ii}),$ $\displaystyle\mbox{Var}[Y_{i}]$
$\displaystyle=$
$\displaystyle\mbox{E}[Y_{i}]+\mbox{E}^{2}[Y_{i}]\bigl{(}\exp(\sigma_{ii})-1\bigr{)},$
$\displaystyle\mbox{Cov}(Y_{i},Y_{j})$ $\displaystyle=$
$\displaystyle\mbox{E}[Y_{i}]\mbox{E}[Y_{j}]\bigl{(}\exp(\sigma_{ij})-1\bigr{)}.$

In general, denote by $y_{(k)}=y(y-1)....(y-k+1)$, the so called falling
factorial. Then (2) shows that multivariate factorial moments are computed by
using simple properties of Poisson distribution, i.e.

$\displaystyle\displaystyle\mbox{E}\Big{[}\prod_{i=1}^{d}Y_{i,(r_{i})}\Big{]}$
$\displaystyle=\mbox{E}\Big{\\{}\prod_{i=1}^{d}\mbox{E}\Big{[}Y_{i,(r_{i})}\mid\lambda_{i}\Big{]}\Big{\\}}=\mbox{E}\Big{[}\prod_{i=1}^{d}\lambda^{r_{i}}_{i}\Big{]},$

Furthermore, following Johnson et al. (1992) joint moments are given by

$\displaystyle\mbox{E}\Bigl{[}\prod_{i=1}^{d}Y^{r_{i}}_{i}\Bigr{]}$
$\displaystyle=\sum_{l_{1}}^{r_{1}}\ldots\sum_{l_{d}=0}^{r_{d}}\prod_{i=1}^{d}s(r_{i},l_{i})\mbox{E}\Bigl{[}\prod_{i=1}^{d}\lambda^{l_{i}}_{i}\Bigr{]},$

where $s(r,l)$ are the Stirling numbers of the second kind.

Next, it is shown that the multivariate negative-multinomial distribution is
recovered by means of (2). This is well-known in the univariate case. A
multivariate negative-multinomial distribution has p.m.f. which is given by

$\mbox{P}[Y=y]=\frac{(n+\sum_{i=1}^{d}y_{i}-1)!}{(\prod_{i=1}^{d}y_{i}!)(n-1)!}p_{0}^{n}\prod_{j=1}^{d}p^{y_{j}}_{j},$
(4)

where $y_{j}=0,1,2,\ldots$ for $j=1,\ldots,d$, $n>0$ and $0<p_{j}<1$,
$j=0,1,\ldots d$ satisfying $p_{0}=1-\sum_{j=1}^{d}p_{j}$. The parameters
$p_{j}$, $j=1,2,\ldots,d$ denote the probabilities of obtaining different
failures and thus $p_{0}$ is the probability of success in an experiment
terminating to $n$ failures; see Johnson et al. (1997, Ch. 36) and Joe (1997,
Ch. 7.2) for more. Note that the parameter $n$ might assume real values in
applications; it is short of ”dispersion” parameter though the concept of
dispersion–that is when the variance exceeds the mean–is quite vague in the
multivariate case we consider. In this case, it can be shown that all pairwise
correlations between the components of $Y$ are positive.

It is shown next that (4) is obtained as mixed Poisson model by using (2).
Suppose that $Y_{j}\sim\mbox{Poisson}(\lambda_{j}\theta)$ given $\theta$, and
suppose that $\theta\sim\mbox{Gamma}(\beta,\beta)$, for $j=1,2,\ldots,d$, so
that $\mbox{E}[\theta]=1$ and $\mbox{Var}[\theta]=1/\beta$. Then

$\displaystyle\mbox{P}[Y=y]$
$\displaystyle=\int_{0}^{\infty}\mbox{P}[Y=y\mid\theta]dG(\theta)$
$\displaystyle=\frac{\Gamma(\beta+\sum_{i=1}^{d}y_{i}-1)}{(\prod_{i=1}^{d}y_{i}!)\Gamma(\beta)}\Bigl{(}\frac{\beta}{\beta+\sum_{i=1}^{d}\lambda_{i}}\Bigr{)}^{\beta}\prod_{j=1}^{d}\Bigl{(}\frac{\lambda_{j}}{\beta+\sum_{i=1}^{d}\lambda_{i}}\Bigr{)}^{y_{j}},$
(5)

where $\Gamma(\cdot)$ denotes the Gamma function. Obviously (4) holds.
Clearly, this model shows that $\theta$ describes common unobserved
heterogeneity; see Munkin and Trivedi (1999) for simulated maximum likelihood
estimation for this particular class.

Recall (4). Then, given a covariate vector $X$, a multinomial logistic
regression model (see Agresti (2002)) is employed to link $X$ with the
probabilities $p_{j}$, $j=1,2,\ldots,(d+1))$. Furthermore, the model can be
extended to include a log-linear model for $n$ (which can be positive real in
general, see the previous case)–for more details see Zhang et al. (2017).

###  2.3 Copula approaches

Copula-based construction of multivariate count distributions is an active
topic of research; see Nikoloulopoulos (2013a) and Inouye et al. (2017) for
nice surveys. Copulas are useful because of Sklar’s theorem (Sklar (1959))
which shows that marginal distributions are combined to give a joint
distribution when applying a copula, i.e. a $d$-dimensional distribution
function all of whose marginals are standard uniforms; the book by Nelsen
(1999) gives a thorough introduction to copulas.

Even though copulas provide an appealing methodology for constructing joint
distribution functions they pose challenging issues when used for discrete
data analysis. First, the presence of ties in count data (several zeroes, for
example) makes them non-identifiable–Genest and Nešlehová (2007), in
particular pp. 507-508–illustrate the lack of identifiability. An additional
issue is that of the likelihood function’s computationally difficulty for
estimating unknown parameters. For a discrete random vector, whose cumulative
distribution function (c.d.f) is $F$, its p.m.f. involves $2^{d}$ finite
differences of $F$, i.e.

$\mbox{P}[Y=y]=\sum_{l_{1}=0,1}\ldots\sum_{l_{d}=0,1}(-1)^{l_{1}+\ldots+l_{d}}~{}\mbox{P}[Y_{1}\leq
y_{1}-l_{1},\cdots,Y_{d}\leq y_{d}-l_{d}].$

Computational methods using Bayesian data augmentation have been developed by
Smith and Khaled (2012), among others; the survey of Smith (2013) provides
references on Bayesian methodology for discrete data copula modeling.

In the rest of this section I describe a different approach for employing
copulas to model multivariate count data. The methodology is based on the work
by Fokianos et al. (2020) who advanced a particular data generating process
for multivariate count time series analysis; this topic is discussed in Sec.
5. Initially, the idea is illustrated for the case of i.i.d random variables.
The intent is to introduce a data generating process which keeps _all_
marginal distributions of the vector $Y$ to be Poisson distributed and, at the
same time, it allows for arbitrary dependence among them. This is accomplished
by appealing to elementary properties of Poisson process. An explicit account
of this construction is given by the following algorithm–recall that $Y_{i}$,
$i=1,2,\ldots,d$ is the $i$’th-component of the count vector $Y$ whose mean is
$\lambda_{i}$.

  1. 1.

Let ${U}^{l}=(U_{1,l},\ldots,U_{d,l})$ for $l=1,2,\ldots,K$, be a sample from
a $d$-dimensional copula $C(u_{1},\ldots,u_{d})$. Then $U_{i,l}$,
$l=1,2,\ldots,K$ follow marginally the uniform distribution on $(0,1)$, for
$i=1,2,\ldots,d$.

  2. 2.

Consider the transformation $X_{i,l}=-{\log
U_{i,l}}/{\lambda_{i}},~{}i=1,2,\ldots,d.$ Then, the marginal distribution of
$X_{i,l}$, $l=1,2,\ldots,K$ is exponential with parameter $\lambda_{i}$,
$i=1,2,\ldots,d$.

  3. 3.

Define now, (taking $K$ large enough), $Y_{i}=\max_{1\leq k\leq
K}\left\\{\sum_{l=1}^{k}X_{i,l}\leq 1\right\\},~{}i=1,2,\ldots,d.$ Then
${Y}=(Y_{1},\ldots,Y_{d})^{T}$ is marginally realization of a Poisson process
with parameter ${\lambda}$.

  4. 4.

Repeat steps 1-3 $n$ times to generate a sample $Y^{j}$, $j=1,2,\ldots,n$.

The algorithm generates i.i.d. random vectors whose dependence among their
components is introduced by a copula structure on the waiting times of the
Poisson process. In other words, the copula is imposed on the uniform random
variables generating the exponential waiting times. The end result gives a
sample of multivariate discrete random variables. The methodology can be
extended to other discrete marginal distributions provided that they can be
generated by continuous inter arrival times. For instance, suppose that
$Y_{i}$ is marginally mixed Poisson with mean $\theta\lambda_{i}$ where
$\theta$ satisfies $\mbox{E}[\theta]=1$. Many families of count distributions,
including the negative binomial, can be generated by this construction; see
Sec. 2.2. Then steps 1-4 of the above algorithm are still useful for
generating multivariate random vectors whose marginals are not necessarily
Poisson. Indeed, generating at the first step an additional random variable
$\theta$, define again at step 2 the waiting times by $X_{i,l}=-{\log
U_{i,l}}/{\theta\lambda_{i}},~{}~{}~{}i=1,2,\ldots,d.$ Then, the distribution
of $X_{i}$ is mixed exponential and therefore steps 3-4 deliver a realization
of a count vector whose marginal distribution is mixed Poisson.

Visualizing the joint p.m.f., based on the above construction, is illustrated
by means of Fig. 1. Plots (a), (b) and (c) show the case of independence,
positive and negative correlation, respectively when applying a Gaussian
copula. The algorithm delivers desired marginals and, in addition, shows that
all type of different correlations can be achieved. Fig. 1(d) shows that joint
p.m.f. of the negative multinomial vector, see (5) with the same $\lambda$
parameters values used for the previous algorithm and $\beta=10$. The positive
correlation between the vector components is obvious.

The algorithm, when implemented, requires clear distinction between the
resulting copula to the vector of counts and the copula imposed on waiting
times. The transformation from waiting times to counts is stochastic, and
while the copula as such is invariant to one-to-one deterministic
transformations, we do not have such a transformation in this case. Hence,
instantaneous correlation among the components of the count vector is not
equal to the correlation induced by the copula imposed to the vector of
waiting times. _Interpretation of the instantaneous correlation found in data
is related to the correlation of the vector of waiting times and should be
done with care._ Recalling Fig. 1 and calculating the sample correlation
coefficient, it is found that $\hat{\rho}=0.34(-0.34)$ when the data are
generated with $\rho>0$ ($\rho<0$), respectively.

Figure 1: Joint p.m.f of a bivariate count distribution using the copula
construction as outlined by steps 1–4. Results are based on a Gaussian copula
with correlation coefficient $\rho$. (a) $\rho=0$ (independence) (b)
$\rho=0.8$ (positive correlation) (c) $\rho=-0.8$ (negative correlation).
Plots are based on 10000 independent observations where the marginals are
Poisson with $\lambda_{1}=3$ and $\lambda_{2}=10$. (d) Joint p.m.f of negative
multinomial distribution (5). Results are based on 10000 independent
observations with $\beta=10$ and $\lambda_{1}=3$ and $\lambda_{2}=10$.

This approach is different from the methodology advocated by Denuit and
Lambert (2005). These authors employ the continued extension, i.e. they add
noise of the form $U-1$, where $U$ is standard uniform, to counts so that
those are transformed into continuous random variables. In doing so, the
copula identifiability problem of is bypassed. An analogous approach, based on
the distributional transform which adds a random jump to the c.d.f. of the
discrete variable, has been studied by Rüschendorf (2013). An interesting
decomposition of the joint p.m.f. of a discrete random vector has been
discussed by Panagiotelis et al. (2012) using the idea of pair-copula
construction (Czado (2010)) by utilizing the concept D-vine copulas (Bedford
and Cooke (2001, 2002)).

Further work on copulas, in the context of generalized linear models, is given
by Song (2000) and Song et al. (2009). The latter reference employs Gaussian
copulas, for multivariate regression analysis of continuous, discrete, and
mixed correlated outcomes under the generalized linear models (GLM) framework
(Nelder and Wedderburn (1972) and McCullagh and Nelder (1989)). More recently,
Yang et al. (2020) consider discrete regression models and copula estimation
arguing that inclusion of continuous covariates implies consistent estimation
of the unknown copula.

###  2.4 Additional models

There are additional approaches for defining a multivariate count
distribution; Joe (1997, Ch.9) and Inouye et al. (2017) review construction of
multidimensional Poisson p.m.f by appealing to full conditional distributions
and Markov random fields; see Besag (1974) for the so called auto-Poisson
model. Additional models include the Sarmanov and bivariate Hurdle
distributions among others; see Cameron and Trivedi (2013, Ch. 8). These
models are mentioned for completeness of presentation but their properties
have not been fully explored in the literature, to the best of my knowledge.

##  3 Integer AR models

Integer Autoregressive (INAR) models deserve special consideration due to the
thinning operation. The calculus of thinning operators provides useful insight
into the probabilistic properties of those processes by employing the simple
device of summing up a _random number_ of integer-valued random variables (see
Steutel and van Harn (1979)). The case that has attracted more attention is
when the summands consist of an independent and identically distributed (iid)
sequence of Bernoulli random variables.

###  3.1 The thinning operator

Define the generalized Steutel and van Harn operator (see Latour (1997, Def.
1.1)) as follows:

###### Definition 3.1

Suppose that $X$ is a non–negative integer random variable. The generalized
thinning operator, denoted by $\circ$, is defined as

$\alpha\circ
X=\left\\{\begin{array}[]{ll}\sum_{k=1}^{X}I_{k},&\hbox{$X>0$;}\\\
0,&\hbox{$X=0$.}\end{array}\right.$

where $\left\\{I_{k},k\in\mathbb{N}\right\\}$ is a sequence of iid integer
random variables–independent of $X$–with mean $\alpha$ and variance $\beta$.

The sequence $\left\\{I_{k},k\in\mathbb{N}\right\\}$ is called counting
series. If $\left\\{I_{k}\right\\}$ is an iid sequence of Bernoulli random
variables, then $\alpha\circ X$ counts the number of successes in a random
number of Bernoulli trials where the probability of success $\alpha$ remains
constant throughout the experiment so that given $X$, $\alpha\circ X$ is a
binomial random variable with parameters $X$ and $\alpha$. Numerous properties
can be proved for the thinning operator, for instance it can be shown that
$\mbox{E}\left[\alpha\circ X\mid X\right]=\alpha X$,
$\mbox{E}\left[\alpha\circ X\right]=\alpha\mbox{E}\left[X\right]$,
$\mbox{Var}\left[\alpha\circ X\mid X\right]=\beta X$ and
$\mbox{Var}\left[\alpha\circ
X\right]=\alpha^{2}\mbox{Var}\left[X\right]+\beta\mbox{E}\left[X\right]$,
provided that appropriate moments of $X$ exist.

Definition 3.1 can be extended to integer $d$-dimensional random vectors. For
$i,j=1,2,\ldots,d$, define $\left\\{I_{ij;k},k\in\mathbb{N}\right\\}$ an array
of counting series such that $\mbox{E}[I_{ij}]=\alpha_{ij}$ and
$\mbox{Var}[I_{ij}]=\beta_{ij}$. Let $A=(\alpha_{ij})$, $B=(\beta_{ij})$ be
the corresponding $d\times d$ matrices. Then, Definition 3.1 can be extended
as follows:

###### Definition 3.2

Suppose that $X=(X_{1},X_{2},\ldots,X_{d})^{T}$ is $d$-dimensional integer-
valued random vector with all components being non-negative and denote by
$A\circ=(\alpha_{ij}\circ)$ a $d\times d$ matrix of thinning operators whose
each element is given by Def. 3.1 with corresponding array of counting series
$\left\\{I_{ij;k},k\in\mathbb{N}\right\\}$. Then the _multivariate thinning_
operator is defined as

$\displaystyle A\circ X$
$\displaystyle=\left(\begin{array}[]{c}\sum_{j=1}^{d}\alpha_{1j}\circ X_{j}\\\
\vdots\\\ \sum_{j=1}^{d}\alpha_{dj}\circ X_{j}\\\ \end{array}\right).$

Similar to the univariate case, it can be shown that $\mbox{E}[A\circ
X]=A\mbox{E}[X]$ and $\mbox{Var}[A\circ
X]=\mbox{diag}(BE[X])+A\mbox{Cov}[X]A^{T}$, where $\mbox{diag}(.)$ denotes a
diagonal matrix of appropriate order.

###  3.2 The multivariate INAR model

The multivariate thinning operator serves as basic tool to develop
multidimensional integer AR models (abbreviated by MINAR) of order $p$. A
$d$-dimensional time series $\\{Y_{t},t\in\mathbb{Z}\\}$ is called
multivariate INAR($p$) process if it satisfies

$\displaystyle Y_{t}$ $\displaystyle=\sum_{i=1}^{p}A_{i}\circ
Y_{t-i}+\epsilon_{t},$ (6)

where $\\{\epsilon_{t},t\in\mathbb{Z}\\}$ is a sequence of iid integer-valued
random vectors with mean $\mbox{E}[\epsilon_{t}]=\mu_{\epsilon}$ and
$\mbox{Var}[\epsilon_{t}]=\Sigma_{\epsilon}$ which is independent of all
thinning operators $A_{i}$, $i=1,2,\ldots,p$ and $A_{p}\neq 0$. Denote by
$I_{d}$ the $d$-dimensional identity matrix. Provided that the roots of
polynomial $\det(I_{d}-A_{1}z-\cdots A_{p}z^{p})$ are all located outside the
unit circle, then Latour (1997, Prop. 3.1) shows that there exists an almost
surely unique integer-valued strictly stationary process that satisfies (6)
and such that $\epsilon_{t}$ is independent of $Y_{s},~{}s<t$. For the case of
$d=1$, this condition is equivalent to $\sum_{i=1}^{p}A_{i}<1$, see Du and Li
(1991). The univariate INAR($p$) processes have been introduced by Al-Osh and
Alzaid (1987), Alzaid and Al-Osh (1990). In the same vein, multivariate moving
average (MINARMA() models are defined but this will not be discussed any
further.

We assume the stability condition for MINAR($p$) processes holds true. Recall
(6) and consider the special case of $p=1$. By taking expectations in both
sides of (6), $\mbox{E}[Y_{t}]=(I_{d}-A_{1})^{-1}\mu_{\epsilon}$. Multiplying
both sides of (6) by $Y_{t+h}$ and taking expectations, it follows that
$\mbox{Cov}(Y_{t},Y_{t+h})=A_{1}^{h}\mbox{Var}(Y_{t})$, where
$\mbox{Var}(Y_{t})=A_{1}\mbox{Var}(Y_{t})A_{1}^{T}+\mbox{diag}(B\mbox{E}[Y_{t}])+\Sigma_{\epsilon}$
for $h=0,\pm 1,\pm 2\cdots$. As a final remark, Latour (1997, Prop. 4.1) shows
that a MINAR($p$), which satisfies the stability condition discussed earlier,
has identical second order properties with an ordinary vector AR($p$) (VAR)
model (see Tsay (2014)). Consider again the case $p=1$ for (6). Then
$\\{Y_{t}\\}$ is represented by a VAR(1) process of the form

$\displaystyle Y_{t}=\omega+A_{1}Y_{t-1}+\zeta_{t}$ (7)

where $\\{\zeta_{t}\\}$ is a white noise process with covariance matrix
$\Sigma_{\zeta}=\mbox{diag}(B\mbox{E}[Y_{t}])+\Sigma_{\epsilon}$ and
$\omega=\mbox{E}[Y_{t}]$. This fact has important consequences on estimation.
For instance least squares estimators (LSE) directly applies to this class of
models but subject to restriction that all unknown coefficients are positive.

###  3.3 Estimation

Besides LSE, likelihood estimation has been also developed for estimating the
unknown parameters of the model (6). Both methods are discussed next by
assuming that $Y_{1},\ldots,Y_{n}$ is a sample from a MINAR(1) model–this is
done mostly for convenience.

LSE are computed and studied by using (7). Let

$\underbrace{Y}_{(n-1)\times
d}=\underbrace{X^{\star}}_{(n-1)\times(d+1)}~{}\underbrace{\beta}_{(d+1)\times
d}+\underbrace{Z}_{(n-1)\times d}$

where the $i$’th row of $Y$, $X^{\star}$ and $Z$ is given by $Y_{i+1}^{T}$,
$(1,Y^{T}_{i})$ and $\zeta^{T}_{i}$, $i=1,2,\ldots,(n-1)$, respectively. The
regression matrix parameter is denoted by $\beta$, i.e.
$\beta=(\omega,A_{1})^{T}$. Then the LSE of $\beta$ is denoted by
$\widetilde{\beta}$ and is equal to

$\displaystyle\begin{aligned} \widetilde{\beta}&=\biggl{(}X^{\star
T}X^{\star}\biggr{)}^{-1}X^{\star
T}Y=\biggl{(}\sum_{t=2}^{n}\begin{bmatrix}1&Y^{T}_{t}\\\
Y_{t}&Y_{t}Y^{T}_{t}\end{bmatrix}\biggr{)}^{-1}\biggl{(}\sum_{t=2}^{n}\begin{bmatrix}Y^{T}_{t}\\\
Y_{t-1}Y^{T}_{t}\end{bmatrix}\biggr{)}.\end{aligned}$ (8)

Based on this define the residual matrix
$\widetilde{Z}=Y-X^{\star}\widetilde{\beta}$ to obtain an estimator of
$\Sigma_{\zeta}$ by

$\displaystyle\widetilde{\Sigma}_{\zeta}=\frac{1}{n-d-2}\widetilde{Z}\widetilde{Z}^{T},$

where the numerator $n-d-2=(n-1)-(d+1)$ is equal to the effective degrees of
freedom minus the number of parameters estimated for each component series.
More precisely, it can be shown that $\mbox{E}[\widetilde{\beta}]=\beta$ and
$\mbox{E}[\widetilde{\Sigma}_{\zeta}]=\Sigma_{\zeta}$. Moreover, as
$n\rightarrow\infty$, and assuming suitable regularity conditions are
fulfilled

$\displaystyle\sqrt{n}\operatorname{vec}\bigl{(}\widetilde{\beta}-\beta\bigr{)}\stackrel{{\scriptstyle
D}}{{\Longrightarrow}}N(0,\Sigma_{\zeta}\otimes H^{-1}),$

where $\operatorname{vec}(.)$ denotes the vec operator, $\otimes$ is the
Kronecker product and the $(d+1)\times(d+1)$ matrix $H$ is the limit (in
probability) of $X^{\star T}X^{\star}/n$. A rigorous statement is developed
along the lines of Latour (1997) and Lütkepohl (2005, Lemma 3.1). This result
holds when the true parameters belong to the interior of parameter space
employing unconstrained optimization. The asymptotic distribution of LSE,
under the constraint that all elements of the matrix $A_{1}$ are positive, is
an open problem.

Besides LSE, conditional likelihood estimation is developed (recall again (6)
with $p=1$) by maximizing the likelihood function

$\displaystyle L(\beta)=\prod_{t=1}^{n}\mbox{P}_{\beta}[Y_{t}=y_{t}\mid
Y_{t-1}=y_{t-1}],$

when imposing a multivariate distribution on the error term $\epsilon_{t}$. In
general the conditional transition is given by the $d$-dimensional convolution

$\displaystyle\mbox{P}_{\beta}[Y_{t}=y_{t}\mid
Y_{t-1}=y_{t-1}]=\sum_{k=0}^{y_{t}}P_{\beta}[A_{1}\circ
Y_{t-1}=y_{t}-k]P[\epsilon_{t}=k].$

Both the above equation implies that the log-likelihood function is given by

$\displaystyle l(\beta)\equiv\log
L(\beta)=\sum_{t=1}^{n}\log\Bigl{(}\sum_{k=0}^{y_{t}}P_{\beta}[A_{1}\circ
Y_{t-1}=y_{t}-k]P[\epsilon_{t}=k]\Bigr{)}.$

This short discussion shows that the task of computing the log-likelihood
function is daunting even in the simple case $p=1$. Some simplifications occur
when we assume that the matrix $A_{1}$ is diagonal and by applying pairwise
likelihood methodology as in Pedeli and Karlis (2013a). Such an approach can
deliver some insights for data analysis but it will be difficult to be
justified for multivariate dynamic systems. Generally, likelihood methods are
not suitable for this class of models because of their complicated structure.
The problem complexity increases when both $p$ and $d$ grow but sparsity ideas
(see Hastie et al. (2015)) will be helpful.

###  3.4 Prediction

Consider again the case of model (6) for $p=1$. We briefly discuss prediction
for this particular case to overcome cumbersome notation. Similar to the case
of VAR($p$) models, a MINAR($p$) model is written as a ”big” MINAR(1) so these
results suffice to develop a general point of view.

The one-step ahead predictor of the MINAR(1) process is easily calculated by

$\displaystyle\mbox{E}[Y_{t+1}\mid Y_{t}]=\mbox{E}\left[A_{1}\circ
Y_{t}+\epsilon_{t}\mid Y_{t}\right]=A_{1}Y_{t}+\mu_{\epsilon}.$

Therefore, by recursion

$\displaystyle\mbox{E}[Y_{t+h}\mid
Y_{t}]=A_{1}^{h}Y_{t}+\bigl{(}I_{d}+A_{1}+\cdots+A_{1}^{h-1}\bigr{)}\mu_{\epsilon}.$

Define $V(h)=\mbox{Var}[Y_{t+h}\mid Y_{t}]$. Then

$\displaystyle V(1)=\mbox{Var}\bigl{(}A_{1}\circ Y_{t}+\epsilon_{t}\mid
Y_{t}\bigr{)}=\mbox{diag}(B\mbox{E}[Y_{t}])+\Sigma_{\epsilon}.$

The law of total variance shows that

$\displaystyle V(h)$
$\displaystyle=\mbox{E}\Bigl{[}\mbox{Var}\bigl{(}Y_{t+h}\mid
Y_{t+1}\bigr{)}\mid
Y_{t}\Bigr{]}+\mbox{Var}\Bigl{[}\mbox{E}\bigl{(}Y_{t+h}\mid
Y_{t+1}\bigr{)}\mid Y_{t}\Bigr{]}$
$\displaystyle=\mbox{E}[V(h-1)]+A^{h-1}_{1}V(1)(A^{h-1}_{1})^{T}$

Provided that $A_{1}$ is diagonal, some simplifications of the previous
formulas have been proved by Pedeli and Karlis (2013b). The complexity of
prediction problem increases as both dimension and model order increase. In
addition, care should be taken when unknown parameters are replaced by their
estimators.

Closing this section, we mention that properties of multivariate INAR models
are well understood for low-dimensional data and under tangible assumptions.
It is interesting to consider properties of the multivariate thinning operator
in high-dimensions and investigate the problem of estimation and prediction.

##  4 Parameter-Driven Models

A parameter driven model, according to classification introduced by Cox
(1981), is a time series driven by an unobserved process (as opposed to past
process values; see Sec 5). For multivariate count series, state-space models
were studied in Jørgensen et al. (1999) and Jung et al. (2011); see
Ravishanker et al. (2014), and Ravishanker et al. (2015), among others, for
more recent contributions. The approaches that have been taken for estimation
are based either on likelihood or full Bayesian methods. We review some of
these works.

The model proposed by Jørgensen et al. (1999), is closely related to the
theory of mixed Poisson distributions, see Sec. 2.2. It assumes that the
conditional distribution of the $i$’th component of the multivariate count
series at time $t$, $Y_{i,t}$, given an unobserved univariate time-varying
process 111For this part of the text, I replace the notation $\lambda_{t}$ by
$\theta_{t}$ because it denotes a univariate time-varying mean process.
$\theta_{t}$, is Poisson distributed with mean $a_{i,t}\theta_{t}$ such that
$a_{i,t}=\exp({c}_{t}^{T}\mathbf{\alpha}_{i})$, where ${c}_{t}$ are
$k$-dimensional time-varying covariate vectors and ${\alpha}_{i}$ are
$k$-dimensional regression parameters, $i=1,\ldots,d$. Assume that
$\theta_{0}=1$ and the conditional distribution of $\theta_{t}$ given
$\theta_{t-1}$ is Gamma with mean $b_{t}\theta_{t-1}$ and a squared
coefficient of variation of form $\sigma^{2}/\theta_{t-1}$. The parameters
$b_{t}$ depend on the so called long–term covariates ${c}_{t}$ through
$b_{t}=\exp(\Delta{c}_{t}^{T}\beta)$, where $\Delta{c}_{t}={c}_{t}-{c}_{t-1}$
and ${c}_{0}=0$ and $\sigma^{2}$ denotes a dispersion parameter. It can be
shown that

$\displaystyle\mbox{E}[\theta_{t}]$ $\displaystyle=b_{1}\ldots b_{t}$

which implies that $\log(\mbox{E}[\theta_{t}])={c}_{t}^{\prime}\beta$. In
addition, for $h\geq 0$,

$\displaystyle\mbox{Var}(\theta_{t})=\phi_{t}\mbox{E}[\theta_{t}]\sigma^{2},~{}~{}\mbox{Cov}(\theta_{t},\theta_{t+h})$
$\displaystyle=\phi_{t}\mbox{E}[\theta_{t+h}]\sigma^{2},$

where $\phi_{t}=b_{t}+b_{t}b_{t-1}+b_{t}b_{t-1}\ldots b_{1}$. Set
${a}_{t}=(a_{1,t},\ldots,a_{d,t})^{T}$ and
$\Theta_{t}=\mbox{Diag}(a_{1,t},\ldots,a_{d,t})$. Then (compare with (3) when
$a_{i,t}=1$ for all $i$ and $t$)

$\displaystyle\mbox{E}[{Y}_{t}]={a}_{t}\mbox{
E}[\theta_{t}],~{}~{}\mbox{Var}({Y}_{t})={\Theta}_{t}\mbox{E}[\theta_{t}]+{a}_{t}{a}_{t}^{T}\phi_{t}\sigma^{2}\mbox{E}[\theta_{t}].$

This last result shows that the variance matrix of ${Y}_{t}$ consists of two
components; (a) a Poisson variance and (b) a type of ”overdispersion”
component. The authors discuss Kalman prediction and filtering, for the log-
linear model
$\mbox{E}[Y_{i,t}]=\exp({c}_{t}^{T}{\alpha}_{i}+{z}_{t}^{T}\beta)$, relying on
previous calculations.

The above model is closely related to the model of Jung et al. (2011) which in
turn generalizes that of Wedel et al. (2003) who developed a comprehensive
class of factor models for multivariate truncated count data. The former
authors assume that $Y_{i,t}$, conditionally on $\lambda_{i,t}$, are
independent Poisson distributed random variables with mean $\lambda_{i,t}$,
$i=1,2,\ldots,d$ and for all $t$. By considering the $d$-dimensional time-
varying vector process $\lambda_{t}=(\lambda_{1,t},\ldots,\lambda_{d,t})^{T}$,
it is assumed that

$\displaystyle\log\lambda_{t}=\omega+\Gamma f_{t},$

where $\omega$ is $d$-dimensional vector of parameters, $\Gamma$ is a $d\times
s$ matrix of factor loadings, $f_{t}$ is an $s$-dimensional vector of latent
random factors and the $\log(.)$ function is taken componentwise. Further, the
components of $f_{t}$ are decomposed to similar subsets which are assumed to
follow independently Gaussian AR(1) model. Jung et al. (2011) develop
estimation under Poisson and Negative binomial distribution by employing
efficient importance sampling and apply this methodology to numbers of trades,
in 5-min intervals, for five New York Stock Exchange stocks from two
industrial sectors.

Similarly, Wang and Wang (2018) assume that
$\mbox{E}[\lambda_{i,t}]=\mbox{E}[\lambda_{i,t}]\varepsilon_{i,t}$, where
$\varepsilon_{i,t}$ is the $i$’th component of a $d$-dimensional hidden
process $\varepsilon_{t}$ such that $\mbox{E}[\varepsilon_{i,t}]=1$ (see also
Zeger (1988), Davis and Wu (2009) and Christou and Fokianos (2014)). This
assumption implies that

$\displaystyle\mbox{Cov}[Y_{i,t},Y_{j,u}]=\mbox{E}[\lambda_{i,t}]\mbox{E}[\lambda_{j,u}]\mbox{Cov}[\varepsilon_{i,t},\varepsilon_{j,u}],$

for any $t,u$ and $i\neq j$ or $t\neq u$ and all $i,j$. Effectively the
autocovariance function of the hidden process is identical to the
autocovariance function of the standardized process
$Y_{i;t}/\mbox{E}[\lambda_{i,t}]$. Wang and Wang (2018) assume further that
$\mbox{E}[a_{i,t}]=\exp({c}_{t}^{T}{\alpha}_{i})$ using the previous notation.
To reduce the dimensionality of hidden process $\varepsilon_{t}$, it is
assumed to satisfy $\varepsilon_{t}=\Gamma f_{t}$ as in Jung et al. (2011) but
with $s=\dim(f_{t})$ unknown. Correlation is taken into account by this
construction since the dynamics of $f_{t}$ drive the time-evolution of
$\varepsilon_{t}$. Inference proceeds in two steps: (a) pseudo-maximum-
likelihood estimation for regression coefficients and (b) identification of
common factor(s) utilizing eigenanalysis on a positive definite matrix.

In a related article, Zhang, Chen, and Li (2017) discuss a model which is
based on the the multivariate lognormal mixture Poisson distribution (see Sec.
2.2) and allows for serial correlations by assuming that the Poisson mean
vector is a latent process driven by a nonlinear autoregressive model. The
authors employ Monte Carlo Expectation Maximization algorithm together with
particle filtering and smoothing methods to develop inference. Similarly, Al-
Wahsh and Hussein (2020) motivated by an application concerning asthma related
visits to emergency rooms, consider a hidden autoregressive process which
drives the dynamics of a positively correlated bivariate time series of counts
whose conditional distribution is assumed to be the multivariate Poisson
distribution (1). The authors use a Bayesian data cloning approach to compute
maximum likelihood estimators and their standard errors.

From a fully Bayesian point of view, Aktekin et al. (2018) (see also Gamerman
et al. (2013)) assume that $Y_{i;t}$ are independent, conditionally on
univariate parameters $\alpha_{i}$ and a process $\varepsilon_{t}$, Poisson
distributed with mean $\alpha_{i}\theta_{t}$. The parameter $\alpha_{i}$ are
individual specific rates and $\varepsilon_{t}$ is a common process that
drives the dynamics of the observed process and satisfies
$\varepsilon_{t}=(\varepsilon_{t-1}/\gamma)d_{t}$, where $\gamma\in(0,1)$ and
$d_{t}$ are independent Beta random variables with suitable parameters. The
authors study, in addition, a negative binomial model and they implement
inference by particle learning algorithm. Another fully Bayesian approach is
that of Berry and West (2020) who introduce models, within the framework of
dynamic GLM (see West and Harrison (1997)), that allow use of time-varying
covariates for binary and Poisson conditionally distributed time series. The
recent review article by West (2020) gives further insight for Bayesian
modeling of multivariate count time series. Other works, along these lines
include Serhiyenko (2015), Ravishanker et al. (2014), Ravishanker et al.
(2015). The previous articles and the recent work of Davis et al. (2021) give
further references and list other approaches.

##  5 Observation-Driven Models

In this section, we discuss observation-driven models, that is processes whose
dynamics are driven by past observations plus noise. A convenient model would
postulate a multivariate conditional count distribution to the observed
process such that likelihood inference is feasible. But the discussion so far
has shown the obstacles of choosing appropriate count distribution. In this
section we will be studying the GLM approach. It will be argued that this
framework generalizes the traditional ARMA methodology (see Shumway and
Stoffer (2011) for example) to the count time series framework. Model fitting
is based, in general, on quasi-likelihood inference; Godambe (1991) Heyde
(1997); therefore testing, diagnostics and all type of likelihood arguments
are directly applicable for this case.

###  5.1 Linear Models

To initiate the discussion, consider the standard VAR(1) model, but in the
context of a multivariate Poisson autoregression as it is discussed next.
Denote by ${\cal F}_{t}$ the $\sigma$–field generated by all past values of
the process $\\{Y_{s},~{}s\leq t\\}$. Let
$\\{\lambda_{t}=(\lambda_{i,t}),~{}i=1,2,\ldots,d,t\in\mathbb{Z}\\}$ be the
corresponding $d$-dimensional intensity process, vis. $\mbox{E}[Y_{t}\mid{\cal
F}_{t-1}]=\lambda_{t}$. The univariate linear autoregressive model discussed
by Rydberg and Shephard (2000), Heinen (2003), Ferland et al. (2006) and
Fokianos et al. (2009), among others, serves as basic building block to
construct a multivariate Poisson linear VAR(1) process by defining

$Y_{i,t}\mid{\cal F}_{t-1}~{}~{}\sim~{}~{}\mbox{marginally
Poisson}(\lambda_{i,t}),~{}~{}~{}~{}~{}\lambda_{t}=\omega+B_{1}Y_{t-1},$ (9)

where ${\omega}$ is a $d$-dimensional vector and $B_{1}$ is a $d\times d$
unknown matrices. The elements of $\omega$ and $B_{1}$ are assumed to be
positive such that $\lambda_{i,t}>0$. for all $i$ and $t$. It is instructive
to consider (9) in more detail. For the simple case $d=2$, it implies that

$\displaystyle\lambda_{1t}$ $\displaystyle=$
$\displaystyle\omega_{1}+b_{1}^{11}Y_{1,(t-1)}+b_{1}^{12}Y_{2,(t-1)},$
$\displaystyle\lambda_{2t}$ $\displaystyle=$
$\displaystyle\omega_{2}+b_{1}^{21}Y_{1,(t-1)}+b_{1}^{22}Y_{2,(t-1)},$

where $\omega_{i}$ is the $i$’th element of $\omega$ and $b_{1}^{ij}$ is the
$(i,j)$th element of $B_{1}$. Then setting $b_{1}^{12}=0$ implies that the
past values of $Y_{2,t}$ do not affect the evolution of $Y_{1,t}$. Similarly,
$b_{1}^{21}=0$ shows that past values of $Y_{1,t}$ do not affect the evolution
of $Y_{2,t}$. These arguments extend naturally to the case $d>2$.

Generating data using (9) is accomplished by imposing (1), for instance; see
Liu (2012) and Pedeli and Karlis (2013b) for some examples. But the discussion
in Sec. 2 shows the challenges of fitting (9) to data under a bivariate
Poisson (and more generally multivariate Poisson) distribution. To overcome
this challenge we appeal to copulas–but other suitable p.m.f are applicable–by
introducing a joint distribution constructed by utilizing the data generating
process described in Sec. 2.3 but taking into account (9). The algorithm is
repeated for completeness of presentation. Let ${\lambda}_{0}$ be a starting
value and assume that $\omega$, $B_{1}$ are given. Then

  1. 1.

Let ${U}^{l}=(U_{1,l},\ldots,U_{d,l})$ for $l=1,2,\ldots,K$, be a sample from
a $d$-dimensional copula $C(u_{d},\ldots,u_{d})$. Then $U_{i,l}$,
$l=1,2,\ldots,K$ follow marginally the uniform distribution on $(0,1)$, for
$i=1,2,\ldots,d$.

  2. 2.

Consider the transformation $X_{i,l}=-{\log
U_{i,l}}/{\lambda_{i,0}},~{}~{}~{}i=1,2,\ldots,d.$ Then, the marginal
distribution of $X_{i,l}$, $l=1,2,\ldots,K$ is exponential with parameter
$\lambda_{i,0}$, $i=1,2,\ldots,d$.

  3. 3.

Define now (taking $K$ large enough) $Y_{i,0}=\max_{1\leq k\leq
K}\left\\{\sum_{l=1}^{k}X_{i,l}\leq 1\right\\},~{}i=1,2,\ldots,d.$ Then
${Y}_{0}=(Y_{1,0},\ldots,Y_{d,0})^{T}$ is marginally a realization of a
Poisson process with parameter ${\lambda}_{0}$.

  4. 4.

Use (9) to obtain ${\lambda}_{1}$.

  5. 5.

Return back to step 1 to obtain ${Y}_{1}$, and so on.

Figure 2: Lag 1 autocorrelation matrix of a bivariate count time series as a
function of the Gaussian copula parameter $\rho$. Results are based on 5000
data points.

Figure 2 shows plots of the sample autocorrelation matrix
$\mbox{Cov}(Y_{t},Y_{t+1})$ where $Y_{t}$ is a bivariate count time series
generated using the above algorithm and model (9) with a Gaussian copula
having parameter $\rho$. The plot reinforces the point raised in Sec. 2.3
where it is noted that interpretation of the instantaneous correlation found
in data is related to the correlation of the vector of waiting times and
should be done with care. This simple example shows that, as $\rho$ varies
between -1 and 1, all the lag 1 correlation functions do not exceed 0.45. But
the plot also shows that this approach takes into account properly the
correlation found in data.

It is worth pointing out that this approach is different when compared to the
work of Heinen and Rengifo (2007). These authors replace the original counts
by employing the continued extension method of Denuit and Lambert (2005), as
it was discussed in Sec. 2.3. The continued extension method of Heinen and
Rengifo (2007) has been investigated in a simulation study by Nikoloulopoulos
(2013b). Other copula-based models for multivariate count distributions with
serial dependence are discussed in Joe (1997, Ch. 8).

From now on, define (9) as vector Integer Autoregressive Conditional
Heteroscedastic model of order 1, i.e. V-INARCH(1) model. The reason for
choosing this terminology will be explained below. Following identical
arguments as those of Fokianos (2012, Sec. 3.1), model (9) is rewritten as

$\displaystyle
Y_{t}=\lambda_{t}+(Y_{t}-\lambda_{t})=\omega+B_{1}Y_{t-1}+\zeta_{t},$ (10)

which shows that the values of $Y_{t}$ depend on $Y_{t-1}$ plus the white
noise sequence $\\{\zeta_{t}\\}$. Indeed, if $\\{Y_{t}\\}$ is assumed to be
stationary, then it is easily shown that $\mbox{E}[\zeta_{t}]=0$,
$\mbox{Var}[\zeta_{t}]=\mbox{E}[\Sigma_{t}]$, where
$\Sigma_{t}=\mbox{Var}[Y_{t}\mid{\cal F}_{t-1}]$ and
$\mbox{Cov}(\zeta_{t},\zeta_{t+k})=0$ for $k\in\mathbb{Z}$. The matrix
$\Sigma_{t}$ is not determined as in the univariate case whereby the
conditional variance of $Y_{t}$ is $\lambda_{t}$ (with some abuse of
notation). Model (9) implies that the diagonal elements of $\Sigma_{t}$ are
$\lambda_{i,t}$ but the off-diagonal entries depend on the copula employed to
generate data but are unknown because of the contemporaneous correlation
between $Y_{i,t}$ and $Y_{j,t}$, $i\neq j$.

Because of the assumed stationarity, (9) shows that
$\mbox{E}[Y_{t}]=\omega+B_{1}\mbox{E}[Y_{t-1}]$. Then
$\mbox{E}[Y_{t}]=(I_{d}-B_{1})^{-1}\omega$, provided that $\rho(B_{1})<1$,
where $\rho(.)$ denotes the spectral radius of a matrix. Furthermore, (see
Lütkepohl (2005, Ch.2))

$\displaystyle Y_{t}$ $\displaystyle=$
$\displaystyle\omega+B_{1}Y_{t-1}+\zeta_{t}$ (11) $\displaystyle=$
$\displaystyle\omega+B_{1}(\omega+B_{1}Y_{t-2}+\zeta_{t-1})+\zeta_{t}$
$\displaystyle=$
$\displaystyle(I_{d}+B_{1})\omega+B_{1}^{2}Y_{t-2}+B_{1}\zeta_{t-1}+\zeta_{t}$
$\displaystyle=$ $\displaystyle\cdots\cdots\cdots\cdots\cdots\cdots$
$\displaystyle=$ $\displaystyle(1+B_{1}+B_{1}^{2}+\cdots
B_{1}^{t})\omega+\sum_{i=0}^{t}B_{1}^{i}\zeta_{t-i}.$

Therefore, as in the case of ordinary VAR(1) model, assuming that
$\rho(B_{1})<1$, we obtain (in mean square sense) by (11) and for large $t$,
the useful one-sided infinite order moving average representation]

$\displaystyle
Y_{t}=(I_{d}-B_{1})^{-1}\omega+\sum_{i=0}^{\infty}B_{1}^{i}\zeta_{t-i}.$ (12)

In addition

$\displaystyle\Gamma_{Y}(h)\equiv\mbox{Cov}(Y_{t},Y_{t+h})=\sum_{i=0}^{\infty}B_{1}^{i+h}\mbox{E}[\Sigma_{t}](B^{i})^{T},~{}~{}~{}~{}h\geq
0,$

Several other results are readily available because of (10), see Lütkepohl
(2005) or Tsay (2014).

Generalizations of (9), such as the V-INARCH($p$) model

$\displaystyle\lambda_{t}=\omega+\sum_{i=1}^{p}B_{i}Y_{t-i},$

or the vector V-INGARCH($p$, $q$) (where ”G” stands for Generalized)

$\displaystyle Y_{i,t}\mid{\cal F}_{t-1}~{}~{}\mbox{is
marginally}~{}~{}\mbox{Poisson}(\lambda_{i,t}),~{}~{}{\lambda}_{t}=\omega+\sum_{i=1}^{p}B_{i}Y_{t-i}+\sum_{j=1}^{q}A_{j}\lambda_{t-j},$
(13)

where $({A_{j}})_{j=1}^{q}$, $({B}_{i})_{i=1}^{p}$ are $d\times d$ unknown
matrices and all the elements of $\omega$, $({A_{j}})_{j=1}^{q}$,
$({B}_{i})_{i=1}^{p}$ are positive such that $\lambda_{i,t}>0$. for all $i$
and $t$, are studied along the previous arguments. The abbreviation ”INGARCH”
for model (13) just indicates its structural connection to ordinary GARCH
model, Bollerslev (1986), because each component of the vector process $Y_{t}$
is distributed as a Poisson random variable. But the mean of a Poisson random
variable equals its variance; therefore the structure of (13) bears some
resemblance to that of a multivariate GARCH model, see Francq and Zakoïan
(2010). Though the term ”V-INGARCH” does not reflect accurately the true data
generating process it will be used as a generalization of the terminology
introduced for univariate models by Ferland et al. (2006).

It is proved that (13) is a VARMA($\max(p,q)$, $q$) process. Recall that
$\zeta_{t}={Y}_{t}-{\lambda}_{t}$. Then assuming first order stationarity of
$Y_{t}$ and taking expectations on both sides of (13), we obtain that
$\mbox{E}[Y_{t}]=(I_{d}-\sum_{i=1}^{p}A_{i}-\sum_{j=1}^{q}B_{j})^{-1}\omega$,
provided that $\rho(\sum_{i=1}^{p}A_{i}+\sum_{j=1}^{q}B_{j})<1$. Then, by
manipulating (10) but for model (13), it is easily shown that

$\displaystyle\biggl{(}{Y}_{t}-\mbox{E}[Y_{t}]\biggr{)}=\sum_{i=1}^{\max(p,q)}\biggl{(}{A}_{i}+{B}_{i}\biggl{)}\biggr{(}{Y}_{t-i}-\mbox{E}[Y_{t}]\biggr{)}+\zeta_{t}-\sum_{j=1}^{q}{A}_{j}\zeta_{t-j},$
(14)

where we set $A_{i}=0_{d}$ if $q<p$ for $i=q+1,\ldots,p$ or $B_{i}=0_{d}$ if
$q>p$ for $i=p+1,\ldots,q$. In the case of $p=q=1$ the one-sided MA($\infty$)
is given by
${Y}_{t}=\mbox{\boldmath$\bf\mu$}+\sum_{j=0}^{\infty}\Phi_{j}\zeta_{t-j}$ with
${\Phi}_{0}={I}_{d}$ and ${\Phi}_{j}=(A_{1}+B_{1})^{j-1}B_{1}$, for $j\geq 1$,
a fact that shows for any $h>0$

$\displaystyle\Gamma_{Y}(h)=\sum_{j=0}^{\infty}\bigl{(}{A}_{1}+{B}_{1}\bigr{)}^{j+h-1}{B}_{1}\mbox{E}[{\Sigma}_{t}]{B}_{1}^{T}\bigl{(}{A}_{1}^{T}+{B}_{1}^{T}\bigr{)}^{j+h-1},$

by using properties of the linear multivariate processes.

###  5.2 Log-linear Models

The log-linear model we consider is the multivariate analogue of the
univariate log-linear model proposed by Fokianos and Tjøstheim (2011). In a
more general form, assume that for each $i=1,2,\ldots,d$

$\displaystyle Y_{i,t}\mid{\cal F}_{t-1}~{}~{}\mbox{is
marginally}~{}~{}\mbox{Poisson}(\lambda_{i,t}),~{}~{}{\nu}_{t}=\omega+\sum_{i=1}^{p}B_{i}\log(Y_{t-i}+1_{d})+\sum_{j=1}^{q}A_{j}\nu_{t-j}$
(15)

where ${\nu}_{t}\equiv\log\lambda_{t}$ is defined component wise (i.e.
$\nu_{i,t}=\log\lambda_{i,t}$) and ${1}_{d}$ denotes the $d$–dimensional
vector which consists of ones. For this model, there is no need to impose any
constraints on the matrix coefficients. Additionally, the log-linear model
accommodates covariates much easier than the linear model which requires that
any such inclusion has to satisfy $\lambda_{i,t}>0$. In the case of (15)
though, if ${Z}_{t}$ is a covariate vector of dimension $d$, then the second
equation of (15) becomes
${\nu}_{t}=\omega+\sum_{i=1}^{p}B_{i}\log(Y_{t-i}+1_{d})+\sum_{j=1}^{q}A_{j}\nu_{t-j}+CX_{t}$
for a $d\times d$ matrix ${C}$. Interpretation of model parameters for the
log-linear model (15) is identical to the case of the linear model but in
terms of the vector process $\nu_{t}$.

It is more challenging to derive formulas for the mean and autocovariances of
model (15). However, some approximations are possible by considering the
process ${W}_{t}\equiv\log\bigl{(}{Y}_{t}+{1}_{d}\bigr{)}$. Indeed, define now
${\zeta_{t}}={W}_{t}-{\nu}_{t}$ and use the results of Fokianos and Tjøstheim
(2011); Fokianos et al. (2020) to see that ${W}_{t}$ is _approximated_ by a
VARMA model of the form

$\displaystyle\biggl{(}{W}_{t}-\mbox{E}[W_{t}]\biggr{)}=\sum_{i=1}^{\max(p,q)}\biggl{(}{A}_{i}+{B}_{i}\biggl{)}\biggr{(}{W}_{t-i}-\mbox{E}[W_{t}]\biggr{)}+\zeta_{t}-\sum_{j=1}^{q}{A}_{j}\zeta_{t-j},$

similar to (14). An approximate formula for the sequence of autocovariance
matrices for ${W}_{t}$ (but not for $Y_{t}$) is then derived but with suitable
adjustments.

###### Remark 5.1

Stability conditions for model (13) have been developed by Liu (2012) (under
the framework of multivariate Poisson distribution (1)) and Fokianos et al.
(2020) under the copula construction as outlined in Sec. 2.3 for the case
$p=q=1$. Recently Debaly and Truquet (2019) have improved these conditions
again considering the copula-based data generating process. Without
introducing any further notation, we note that the condition
$\rho(\sum_{i}A_{i}+\sum_{j}B_{j})<1$ guarantees stability of the process. For
the log-linear model (15) the desired conditions are more complicated; see
Fokianos et al. (2020) who consider the case $p=q=1$ and prove that either
$\||A_{1}\||_{2}+\||B_{1}\||_{2}<1$ or
$\||{A}_{1}\||_{1}+\||{B}_{1}\||_{1}<1$, where
$\||{A}\||_{d}=\max_{\|{x}\|_{d}=1}\|A{x}\|_{d}$, guarantee ergodicity of the
process. Related stability conditions are discussed in Debaly and Truquet
(2019).

The main notions used to derive such conditions are those of Markov chain
theory (Meyn and Tweedie (1993)), weak dependence (Doukhan and Louhichi
(1999), Dedecker et al. (2007)) and convergence of backward iterations of
random maps (Wu and Shao (2004)). Following the discussions of Neumann (2011)
and Tjøstheim (2012, 2015), the main difficulty is that the process itself
consists of integer valued random variables; however the mean process takes
values on the positive real line and therefore it is quite challenging to
prove stability of the joint process (see also Andrews (1984)). The study of
theoretical properties of univariate models was initiated by the perturbation
method suggested in Fokianos et al. (2009) and was further developed in
Neumann (2011) (using the notion of $\beta$-mixing), Doukhan et al. (2012)
(weak dependence approach), Woodard et al. (2011) and Douc et al. (2013)
(Markov chain theory without irreducibility assumptions) and Wang et al.
(2014) (based on the theory of $e$-chains).

###  5.3 Quasi-Likelihood Inference

Suppose that $\\{Y_{t},t=1,2,\ldots,n\\}$ is an available sample from a count
time series and for the sake of presentation assume model (13) for $p=q=1$.
Inference is analogously developed to the case of log–linear model and for
$p,q>1$. Denote by
${\theta}=({d}^{T},\operatorname{vec}^{T}(A_{1}),\operatorname{vec}^{T}(B_{1}))$,
$\dim(\theta)\equiv\kappa=d(1+2d)$. Following Fokianos et al. (2020), the
estimation problem is approached by employing the theory of estimating
functions. Consider the following conditional quasi–likelihood function, given
a starting value $\lambda_{0}$, for the parameter vector ${\theta}$,

$\displaystyle
L({\theta})=\prod_{t=1}^{n}\prod_{i=1}^{d}\Bigl{\\{}\frac{\exp(-\lambda_{i,t}({\theta}))\lambda_{i,t}^{y_{i,t}}({\theta})}{y_{i,t}!}\Bigr{\\}},$

which is identical to consider (13) assuming independence among time series.
This is strong assumption yet it simplifies computation of estimators and
their respective standard errors. At the same time, it guarantees consistency
and asymptotic normality of the maximizer. Furthermore, the dependence
structure in (13) and (15) is taken into account through because of the
dependence of the likelihood function on the matrices ${A_{1}}$ and ${B_{1}}$;
see Fokianos et al. (2020) for more. The quasi log-likelihood function is
equal to

$\displaystyle
l({\theta})=\sum_{t=1}^{n}\sum_{i=1}^{d}\Bigl{(}y_{i,t}\log\lambda_{i,t}({\theta})-\lambda_{i,t}({\theta})\Bigr{)}.$

We denote by $\widehat{{\theta}}\equiv\arg\max_{{\theta}}l({\theta}),$ the
QMLE of ${\theta}$. The score function is given by

$\displaystyle
S_{n}(\theta)=\sum_{t=1}^{n}\sum_{i=1}^{d}\Bigl{(}\frac{y_{i,t}}{\lambda_{i,t}({\theta})}-1\Bigr{)}\frac{\partial\lambda_{i,t}(\theta)}{\partial{\theta}}=\displaystyle\sum_{t=1}^{n}\frac{\partial{{\lambda}}^{T}_{t}({\theta})}{\partial{{\theta}}}{D}^{-1}_{t}({\theta})\Bigl{(}{Y}_{t}-{{\lambda}}_{t}({{\theta}})\Bigr{)}$
(16)

where ${\partial{{\lambda}}_{t}}/{\partial{{\theta}^{T}}}$ is a
$d\times\kappa$ matrix and ${D}_{t}$ is the $d\times d$ diagonal matrix with
the $i$’th diagonal element equal to $\lambda_{i,t}({\bf\theta})$,
$i=1,2,\ldots,p$. Furthermore

$\displaystyle\frac{\partial{{\lambda}}_{t}}{\partial{{d}^{T}}}$
$\displaystyle=$
$\displaystyle{I}_{d}+{A_{1}}\frac{\partial{{\lambda}}_{t-1}}{\partial{{d}}^{T}},$
$\displaystyle\frac{\partial{{\lambda}}_{t}}{\partial\operatorname{vec}^{T}({A_{1}})}$
$\displaystyle=$
$\displaystyle({\lambda}_{t-1}\otimes{I}_{d})^{T}+{A_{1}}\frac{\partial{{\lambda}}_{t-1}}{\partial\operatorname{vec}^{T}({A_{1}})},$
$\displaystyle\frac{\partial{{\lambda}}_{t}}{\partial\operatorname{vec}^{T}({B_{1}})}$
$\displaystyle=$
$\displaystyle({Y}_{t-1}\otimes{I}_{d})^{T}+{A_{1}}\frac{\partial{{\lambda}}_{t-1}}{\partial\operatorname{vec}^{T}({B_{1}})}.$

The Hessian matrix is given by

$\displaystyle{H}_{n}({\theta})$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\sum_{i=1}^{p}\frac{y_{i;t}}{\lambda_{i,t}^{2}({\theta})}\frac{\partial\lambda_{i,t}({\theta})}{\partial{\theta}}\frac{\partial\lambda_{i,t}({\theta})}{\partial{\theta}^{T}}-\sum_{t=1}^{n}\sum_{i=1}^{p}\Bigl{(}\frac{y_{i,t}}{\lambda_{i,t}({\theta})}-1\Bigr{)}\frac{\partial^{2}\lambda_{i,t}(\theta)}{\partial{\theta}\partial{\theta}^{T}}.$
(17)

Therefore, the conditional information matrix is equal to

$\displaystyle{G}_{n}({\theta})$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\frac{\partial{{\lambda}}^{T}_{t}({\theta})}{\partial{{\theta}}}{D}^{-1}_{t}({\theta}){\Sigma}_{t}({\theta}){D}^{-1}_{t}({\theta})\frac{\partial{{\lambda}}_{t}({\theta})}{\partial{{\theta}^{T}}},$
(18)

where the matrix ${\Sigma}_{t}(\cdot)$ denotes the _true_ covariance matrix of
the vector ${Y}_{t}$. In case that the process $\\{{Y}_{t}\\}$ consists of
uncorrelated components then ${\Sigma}_{t}({\theta})={D}_{t}({\theta})$. Under
suitable conditions, Fokianos et al. (2020) show that

$\sqrt{n}(\widehat{{\theta}}-{\theta}_{0})\xrightarrow{\text{d}}N(0,{H}^{-1}{G}{H}^{-1})$

where the matrices ${H}$ and ${G}$ are defined by the limits (in probability)
of (17) and (18), respectively. The same result is true for the log-linear
model (15); details are omitted. To estimate the copula parameter, it is
desirable to compare the conditional distribution of
${Y}_{t}\mid{\lambda}_{t}$ to that of ${Y}^{*}_{t}\mid{\lambda}_{t}$, where
${Y}^{*}_{t}$ is a count time series generated by a suitable choice of a
copula. There are several ways of comparing such distributions and this topic
is still under investigation. In Fokianos et al. (2020) an initial approach,
based on the newly developed concept of local Gaussian correlation (see
Berentsen et al. (2014)) was shown to be satisfactory. But the problem of
estimating the copula parameter remains unexplored.

###### Remark 5.2

Equation (16) motivates a more general framework that can be applied to the
analysis of multivariate count time series modes. A natural generalization, is
to consider the following estimating functions

$\displaystyle
S_{{v}}({\theta})=\sum_{t=1}^{n}\frac{\partial{{\lambda}}^{T}_{t}({\theta})}{\partial{{\theta}}}{V}^{-1}_{t}(\rho,{\lambda}_{t}({\theta}))\Bigl{(}{Y}_{t}-{{\lambda}}_{t}(\theta)\Bigr{)},$
(19)

where the notation is completely analogous to (16) and
${V}_{t}(\rho,{\lambda}_{t}({\theta}))$ is a $d\times d$ ”working” conditional
covariance matrix which depend upon the process $\\{{\lambda}_{t}\\}$ and
possibly some other parameters $\rho$. Several choices for the working
conditional covariance matrix are available in the literature; we list some
possibilities. If ${V}={I}_{d}$, then (19) corresponds to a least squares
minimization problem for estimating ${\theta}$. If $V=D_{t}$ then we obtain
(16). More generally, the choice

${V}(\rho,{{\lambda}})=\left(\begin{array}[]{cccc}\lambda_{1,t}&\rho_{12}\sqrt{\lambda_{1,t}}\sqrt{\lambda_{2,t}}&\cdots&\rho_{1d}\sqrt{\lambda_{1,t}}\sqrt{\lambda_{d,t}}\\\
\rho_{12}\sqrt{\lambda_{1,t}}\sqrt{\lambda_{2,t}}&\lambda_{2,t}&\cdots&\rho_{2d}\sqrt{\lambda_{2,t}}\sqrt{\lambda_{d,t}}\\\
\cdots&\cdots&\cdots&\cdots\\\
\rho_{1p}\sqrt{\lambda_{1,t}}\sqrt{\lambda_{p,t}}&\rho_{2p}\sqrt{\lambda_{2,t}}\sqrt{\lambda{p,t}}&\cdots&\lambda_{d,t}\\\
\end{array}\right)$

yields to a constant conditional correlation type of model for multivariate
count time series, see Teräsvirta et al. (2010), among others. This topic
deserved further research; a possible method might rely on the Francq and
Zakoïan (2016) who consider estimation of multivariate volatility models
equation by equation.

##  6 Some remarks

This article communicates some recent progress in the area of multivariate
count time series and advocates the point of view that a successful approach
towards to the resolution of the aforementioned problems is via the theory of
GLM. As it was mentioned in the introduction, it is envisaged that this review
will motivate further research on modeling and inference for multivariate
count time series. For instance, simple questions like predictions,
diagnostics, testing or development of other more suitable models with/without
covariates will require further studies and theoretical developments.

## Acknowledgments

Many thanks to M. Armillotta who provided several comments that improved the
presentation. This work has been funded by the European Regional development
Fund and the Republic of Cyprus through the Research and innovation
Foundation, under the project INFRASTRUCTURES/1216/0017 (IRIDA).

## References

  * Agresti (2002) Agresti, A. (2002).  Categorical Data Analysis (Second ed.).  New York: John Wiley & Sons.
  * Ahmad (2016) Ahmad, A. (2016).  Contributions à l’éconemétrie des séries temporelles à valeurs entières.  Ph. D. thesis, University Charles De Gaulle-Lille III, France.
  * Ahmad and Franq (2016) Ahmad, A. and C. Franq (2016).  Poisson QMLE of count time series models.  Journal of Time Series Analysis 37, 291–314.
  * Aktekin et al. (2018) Aktekin, T., N. Polson, and R. Soyer (2018).  Sequential bayesian analysis of multivariate count data.  Bayesian Analysis 13, 385 – 409.
  * Al-Osh and Alzaid (1987) Al-Osh, M. A. and A. A. Alzaid (1987).  First-order integer-valued autoregressive (INAR($1$)) process.  Journal of Time Series Analysis 8, 261–275.
  * Al-Wahsh and Hussein (2020) Al-Wahsh, H. and A. Hussein (2020).  A bivariate autoregressive poisson model and its application to asthma-related emergency room visits.  Statistics in Medicine 39, 3184–3194.
  * Alzaid and Al-Osh (1990) Alzaid, A. A. and M. Al-Osh (1990).  An integer-valued $p$th-order autoregressive structure (INAR($p$)) process.  Journal of Applied Probability 27, 314–324.
  * Andreassen (2013) Andreassen, C. M. (2013).  Models and inference for correlated count data.  Ph. D. thesis, Aaarhus University, Denmark.
  * Andrews (1984) Andrews, D. (1984).  Non-strong mixing autoregressive processes.  Journal of Applied Probability 21, 930–934.
  * Bedford and Cooke (2001) Bedford, T. and R. M. Cooke (2001).  Probability density decomposition for conditionally dependent random variables modeled by vines.  Annals of Mathematics and Artificial Intelligence 32, 245–268.
  * Bedford and Cooke (2002) Bedford, T. and R. M. Cooke (2002).  Vines–a new graphical model for dependent random variables.  Annals of Statistics 30, 1031–1068.
  * Berentsen et al. (2014) Berentsen, G. D., B. Støve, D. Tjøstheim, and T. Nordbø (2014).  Recognizing and visualizing copulas: an approach using local Gaussian approximation.  Insurance: Mathematics & Economics 57, 90–103.
  * Berry and West (2020) Berry, L. R. and M. West (2020).  Bayesian forecasting of many count-valued time series.  Journal of Business & Economic Statistics 38, 872–887.
  * Besag (1974) Besag, J. (1974).  Spatial interaction and the statistical analysis of lattice systems.  Journal of the Royal Statistical Society, Series B 36, 192–236.  with discussion.
  * Bollerslev (1986) Bollerslev, T. (1986).  Generalized autoregressive conditional heteroskedasticity.  Journal of Econometrics 31, 307–327.
  * Cameron and Trivedi (2013) Cameron, A. C. and P. K. Trivedi (2013).  Regression analysis of count data (Second ed.), Volume 53 of Econometric Society Monographs.  Cambridge University Press, Cambridge.
  * Christou and Fokianos (2014) Christou, V. and K. Fokianos (2014).  Quasi-likelihood inference for negative binomial time series models.  Journal of Time Series Analysis 35, 55–78.
  * Cox (1981) Cox, D. R. (1981).  Statistical analysis of time series: Some recent developments.  Scandinavian Journal of Statistics 8, 93–115.
  * Czado (2010) Czado, C. (2010).  Pair-copula constructions of multivariate copulas.  In P. Jaworski, F. Durante, W. K. Härdle, and T. Rychlik (Eds.), Copula Theory and Its Applications, pp. 93–109. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Dahlhaus (1997) Dahlhaus, R. (1997).  Fitting time series models to nonstationary processes.  Ann. Statist. 25, 1–37.
  * Darolles et al. (2019) Darolles, S., G. L. Fol, Y. Lu, and R. Sun (2019).  Bivariate integer-autoregressive process with an application to mutual fund flows.  Journal of Multivariate Analysis 173, 181–203.
  * Davis and Wu (2009) Davis, R. and R. Wu (2009).  A negative binomial model for time series of counts.  Biometrika 96, 735–749.
  * Davis et al. (2021) Davis, R. A., K. Fokianos, S. H. Holan, H. Joe, J. Livsey, R. Lund, V. Pipiras, and N. Ravishanker (2021).  Count time series: A methodological review.  Journal of the American Statistical Association.  to appear.
  * Davis et al. (2016) Davis, R. A., S. H. Holan, R. Lund, and N. Ravishanker (Eds.) (2016).  Handbook of Discrete-Valued Time Series.  Handbooks of Modern Statistical Methods. London: Chapman & Hall/CRC.
  * Davis and Liu (2016) Davis, R. A. and H. Liu (2016).  Theory and inference for a class of observation-driven models with application to time series of counts.  Statistica Sinica 26, 1673–1707.
  * Debaly and Truquet (2019) Debaly, Z. M. and L. Truquet (2019).  Stationarity and moment properties of some multivariate count autoregressions.  arXiv preprint arXiv:1909.11392.
  * Dedecker et al. (2007) Dedecker, J., P. Doukhan, G. Lang, J. R. León R., S. Louhichi, and C. Prieur (2007).  Weak dependence: with examples and applications, Volume 190 of Lecture Notes in Statistics.  New York: Springer.
  * Denuit and Lambert (2005) Denuit, M. and P. Lambert (2005).  Constraints on concordance measures in bivariate discrete data.  Journal of Multivariate Analysis 93, 40–57.
  * Douc et al. (2013) Douc, R., P. Doukhan, and E. Moulines (2013).  Ergodicity of observation-driven time series models and consistency of the maximum likelihood estimator.  Stochastic Processes and their Applications 123, 2620–2647.
  * Douc et al. (2017) Douc, R., K. Fokianos, and E. Moulines (2017).  Asymptotic properties of quasi-maximum likelihood estimators in observation-driven time series models.  Electron. J. Stat. 11, 2707–2740.
  * Doukhan et al. (2012) Doukhan, P., K. Fokianos, and D. Tjøstheim (2012).  On weak dependence conditions for Poisson autoregressions.  Statistics & Probability Letters 82, 942–948.  with a correction in Vol. 83, pp. 1926-1927.
  * Doukhan and Louhichi (1999) Doukhan, P. and S. Louhichi (1999).  A new weak dependence condition and applications to moment inequalities.  Stochastic Processes and their Applications 84, 313–342.
  * Du and Li (1991) Du, J. G. and Y. Li (1991).  The integer-valued autoregressive INAR(p) model.  Journal of Time Series Analysis 12, 129–142.
  * Durbin and Koopman (2000) Durbin, J. and S. J. Koopman (2000).  Time series analysis of non-Gaussian observations based on state space models from both classical and Bayesian perspective.  Journal of the Royal Statistical Society, Series B 62, 3–56.  with discussion.
  * Fahrmeir and Tutz (2001) Fahrmeir, L. and G. Tutz (2001).  Multivariate Statistical Modelling Based on Generalized Linear Models (2nd ed.).  New York: Springer.
  * Ferland et al. (2006) Ferland, R., A. Latour, and D. Oraichi (2006).  Integer–valued GARCH processes.  Journal of Time Series Analysis Analysis 27, 923–942.
  * Fokianos (2012) Fokianos, K. (2012).  Count time series models.  In T. S. Rao, S. S. Rao, and C. R. Rao (Eds.), Handbook of Statistics: Time Series Analysis–Methods and Applications, Volume 30, pp. 315–347. Amsterdam: Elsevier B. V.
  * Fokianos et al. (2009) Fokianos, K., A. Rahbek, and D. Tjøstheim (2009).  Poisson autoregression.  Journal of the American Statistical Association 104, 1430–1439.
  * Fokianos et al. (2020) Fokianos, K., B. Støve, D. Tjøstheim, and P. Doukhan (2020).  Multivariate count autoregressions.  Bernoulli 26, 471–499.
  * Fokianos and Tjøstheim (2011) Fokianos, K. and D. Tjøstheim (2011).  Log–linear Poisson autoregression.  Journal of Multivariate Analysis 102, 563–578.
  * Francq and Zakoïan (2010) Francq, C. and J.-M. Zakoïan (2010).  GARCH models: Stracture, Statistical Inference and Financial Applications.  United Kingdom: Wiley.
  * Francq and Zakoïan (2016) Francq, C. and J.-M. Zakoïan (2016).  Estimating multivariate volatility models equation by equation.  J. R. Stat. Soc. Ser. B. Stat. Methodol. 78, 613–635.
  * Franke and Rao (1995) Franke, J. and T. S. Rao (1995).  Multivariate first-order integer values autoregressions.  Technical report, Department of Mathematics, UMIST.
  * Früwirth-Schnatter and Wagner (2006) Früwirth-Schnatter, S. and H. Wagner (2006).  Auxiliary mixture sampling for parameter-driven models of time series of counts with applications to state space modelling.  Biometrika 93, 827–841.
  * Gamerman et al. (2013) Gamerman, D., T. R. dos Santos, and G. C. Franco (2013).  A non-gaussian family of state-space models with exact marginal likelihood.  Journal of Time Series Analysis 34, 625–645.
  * Genest and Nešlehová (2007) Genest, C. and J. Nešlehová (2007).  A primer on copulas for count data.  Astin Bull. 37, 475–515.
  * Godambe (1991) Godambe, V. P. (1991).  Estimating Functions.  Oxford: Oxford Science Publications.
  * Gouriéroux and Lu (2019) Gouriéroux, C. and Y. Lu (2019).  Negative binomial autoregressive process with stochastic intensity.  Journal of Time Series Analysis 40, 225–247.
  * Hall et al. (2019) Hall, E. C., G. Raskutti, and R. M. Willett (2019).  Learning high-dimensional generalized linear autoregressive models.  IEEE Transactions on Information Theory 65, 2401–2422.
  * Harvey and Fernandes (1989) Harvey, A. C. and C. Fernandes (1989).  Time series models for count or qualitative observations.  Journal of Business & Economic Statistics 7, 407–422.  with discussion.
  * Hastie et al. (2015) Hastie, T., R. Tibshirani, and M. Wainwright (2015).  Statistical Learning with Sparsity:the Lasso and Generalizations.  Boca Raton: CRC Press.
  * Heinen (2003) Heinen, A. (2003).  Modelling time series count data: An autoregressive conditional poisson model.  Technical Report MPRA Paper 8113, University Library of Munich, Germany.  availabel at http://mpra.ub.uni-muenchen.de/8113/.
  * Heinen and Rengifo (2007) Heinen, A. and E. Rengifo (2007).  Multivariate autoregressive modeling of time series count data using copulas.  Journal of Empirical Finance 14, 564 – 583.
  * Heyde (1997) Heyde, C. C. (1997).  Quasi-Likelihood and its Applications: A General Approach to Optimal Parameter Estimation.  New York: Springer.
  * Inouye et al. (2017) Inouye, D. I., E. Yang, G. I. Allen, and P. Ravikumar (2017).  A review of multivariate distributions for count data derived from the Poisson distribution.  Wiley Interdisciplinary Reviews. Computational Statistics (WIREs) 9(3), e1398, 25.
  * Joe (1997) Joe, H. (1997).  Multivariate Models and Dependence Concepts.  London: Chapman & Hall.
  * Johnson et al. (1997) Johnson, N. L., S. Kotz, and N. Balakrishnan (1997).  Discrete multivariate distributions.  John Wiley, New York.
  * Johnson et al. (1992) Johnson, N. L., S. Kotz, and A. W. Kemp (1992).  Univariate Discrete Distributions (second ed.).  New York: Wiley.
  * Jørgensen et al. (1996) Jørgensen, B., S. Lundbye-Christensen, P. X.-K. Song, and L. Sun (1996).  State-space models for multivariate longitudinal data of mixed types.  The Canadian Journal of Statistics 24, 385–402.
  * Jørgensen et al. (1999) Jørgensen, B., S. Lundbye-Christensen, P. X.-K. Song, and L. Sun (1999).  A state space model for multivariate longitudinal count data.  Biometrika 86, 169–181.
  * Jung et al. (2011) Jung, R., R. Liesenfeld, and J.-F. Richard (2011).  Dynamic factor models for multivariate count data: an application to stock–market trading activity.  Journal of Business & Economic Statistics 29, 73–85.
  * Karlis (2003) Karlis, D. (2003).  An EM algorithm for multivariate Poisson distribution and related models.  J. Appl. Stat. 30(1), 63–77.
  * Karlis (2016) Karlis, D. (2016).  Modelling multivariate times series for counts.  In R. Davis, S. Holan, R. Lund, and N. Ravishanker (Eds.), Handbook of Discrete-Valued Time Series, Handbooks of Modern Statistical Methods, pp. 407–424. London: Chapman & Hall.
  * Karlis and Meligkotsidou (2007) Karlis, D. and L. Meligkotsidou (2007).  Finite mixtures of multivariate Poisson distributions with application.  J. Statist. Plann. Inference 137, 1942–1960.
  * Kedem and Fokianos (2002) Kedem, B. and K. Fokianos (2002).  Regression Models for Time Series Analysis.  Hoboken, NJ: Wiley.
  * Kocherlakota and Kocherlakota (1992) Kocherlakota, S. and K. Kocherlakota (1992).  Bivariate discrete distributions.  Marcel Dekker, Inc., New York.
  * Koochemeshkian et al. (2020) Koochemeshkian, P., N. Zamzami, and N. Bouguila (2020).  Flexible distribution-based regression models for count data: Application to medical diagnosis.  Cybernetics and Systems 51, 442–466.
  * Latour (1997) Latour, A. (1997).  The multivariate GINAR(p) process.  Advances in Applied Probability 29, 228–248.
  * Lee et al. (2018) Lee, Y., S. Lee, and D. Tjøstheim (2018).  Asymptotic normality and parameter change test for bivariate Poisson INGARCH models.  TEST 27, 52–69.
  * Liu (2012) Liu, H. (2012).  Some models for time series of counts.  Ph. D. thesis, Columbia University, USA.
  * Livsey et al. (2018) Livsey, J., R. Lund, S. Kechagias, and V. Pipiras (2018).  Multivariate integere-valued time series with flexible autocovariances and their application to major hurricane counts.  Annals of Applied Statistics 12, 408–431.
  * Lütkepohl (2005) Lütkepohl, H. (2005).  New Introduction to Multiple Time Series Analysis.  Berlin: Springer-Verlag.
  * Marshall and Olkin (1988) Marshall, A. W. and I. Olkin (1988).  Families of multivariate distributions.  Journal of the American Statistical Association 83, 834–841.
  * McCullagh and Nelder (1989) McCullagh, P. and J. A. Nelder (1989).  Generalized Linear Models (2nd ed.).  London: Chapman & Hall.
  * Meyn and Tweedie (1993) Meyn, S. P. and R. L. Tweedie (1993).  Markov Chains and Stochastic Stability.  London: Springer.
  * Munkin and Trivedi (1999) Munkin, M. K. and P. K. Trivedi (1999).  Simulated maximum likelihood estimation of multivariate mixed-poisson regression models, with application.  The Econometrics Journal 2, 29–48.
  * Nelder and Wedderburn (1972) Nelder, J. A. and R. W. M. Wedderburn (1972).  Generalized linear models.  Journal of the Royal Statistical Society, Series A 135, 370–384.
  * Nelsen (1999) Nelsen, R. B. (1999).  An Introduction to Copulas, Volume 139 of Lecture Notes in Statistics.  New York: Springer-Verlag.
  * Neumann (2011) Neumann, M. (2011).  Absolute regularity and ergodicity of Poisson count processes.  Bernoulli 17, 1268–1284.
  * Nikoloulopoulos (2013a) Nikoloulopoulos, A. K. (2013a).  On the estimation of normal copula discrete regression models using the continuous extension and simulated likelihood.  Journal of Statistical Planning and Inference 143, 1923–1937.
  * Nikoloulopoulos (2013b) Nikoloulopoulos, A. K. (2013b).  On the estimation of normal copula discrete regression models using the continuous extension and simulated likelihood.  Journal of Statistical Planning and Inference 143, 1923–1937.
  * Opschoor et al. (2020) Opschoor, A., A. Lucas, I. Barra, and D. van Dijk (2020).  Closed-form multi-factor copula models with observation-driven dynamic factor loadings.  Journal of Business & Economic Statistics.  to appear.
  * Panagiotelis et al. (2012) Panagiotelis, A., C. Czado, and H. Joe (2012).  Pair copula constructions for multivariate discrete data.  Journal of the American Statistical Association 107, 1063–1072.
  * Paul et al. (2008) Paul, M., L. Held, and A. M. Toschke (2008).  Multivariate modelling of infectious disease surveillance data.  Statistics in Medicine 27, 6250–6267.
  * Pedeli and Karlis (2013a) Pedeli, X. and D. Karlis (2013a).  On composite likelihood estimation of a multivariate INAR(1) model.  Journal of Time Series Analysis 34, 206–220.
  * Pedeli and Karlis (2013b) Pedeli, X. and D. Karlis (2013b).  Some properties of multivariate INAR(1) processes.  Computational Statistics & Data Analysis 67, 213 – 225\.
  * Piancastelli et al. (2020) Piancastelli, L. S. C., W. Barreto-Souza, and H. Ombao (2020).  Flexible bivariate INGARCH process with a broad range of contemporaneous correlation.  Technical report.  available at https://arxiv.org/abs/2011.08799.
  * Ravishanker et al. (2014) Ravishanker, N., V. Serhiyenko, and M. R. Willig (2014).  Hierarchical dynamic models for multivariate times series of counts.  Statistics and its Interface 7, 559–570.
  * Ravishanker et al. (2015) Ravishanker, N., R. Venkatesan, and S. Hu (2015).  Dynamic models for time series of counts with a marketing application.  In R. Davis, S. Holan, R. Lund, and N. Ravishanker (Eds.), Handbook of Discrete-Valued Time Series, Handbooks of Modern Statistical Methods, pp. 425–446. London: Chapman & Hall.
  * Rüschendorf (2013) Rüschendorf, L. (2013).  Copulas, sklar’s theorem, and distributional transform.  In Mathematical Risk Analysis: Dependence, Risk Bounds, Optimal Allocations and Portfolios, pp. 3–34. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Rydberg and Shephard (2000) Rydberg, T. H. and N. Shephard (2000).  A modeling framework for the prices and times of trades on the New York stock exchange.  In W. J. Fitzgerlad, R. L. Smith, A. T. Walden, and P. C. Young (Eds.), Nonlinear and Nonstationary Signal Processing, pp. 217–246. Cambridge: Isaac Newton Institute and Cambridge University Press.
  * Scotto et al. (2014) Scotto, M. G., C. H. Weiß, M. E. Silva, and I. Pereira (2014).  Bivariate binomial autoregressive models.  Journal of Multivariate Analysis 125, 233–251.
  * Serhiyenko (2015) Serhiyenko, V. (2015).  Dynamic Modeling of Multivariate Counts - Fitting, Diagnostics, and Applications.  Ph. D. thesis, University of Connecticut, USA.
  * Shumway and Stoffer (2011) Shumway, R. H. and D. S. Stoffer (2011).  Time Series Analysis and its Applications (Third ed.).  New York: Springer.  With R examples.
  * Sklar (1959) Sklar, A. (1959).  Fonctions de répartition à $n$ dimensions et leurs marges.  Publ. Inst. Statist. Univ. Paris 8, 229–231.
  * Smith (2013) Smith, M. S. (2013).  Bayesian approaches to copula modelling.  In P. Damien, P. Dellaportas, N. G. Polson, and D. A. Stephens (Eds.), Bayesian Theory and Applications, pp. 336–358. USA: Oxford University Press.
  * Smith and Khaled (2012) Smith, M. S. and M. A. Khaled (2012).  Estimation of copula models with discrete margins via Bayesian data augmentation.  Journal of the American Statistical Association 107, 290–303.
  * Song (2000) Song, P. X.-K. (2000).  Multivariate dispersion models generated from Gaussian copula.  Scandinavian Journal of Statistics 27, 305–320.
  * Song et al. (2009) Song, P. X.-K., M. Li, and Y. Yuan (2009).  Joint regression analysis of correlated data using Gaussian copulas.  Biometrics 65, 60–68.
  * Steutel and van Harn (1979) Steutel, F. W. and K. van Harn (1979).  Discrete analogues of self-decomposability and stability.  The Annals of Probability 7, 893–899.
  * Teräsvirta et al. (2010) Teräsvirta, T., D. Tjøstheim, and C. W. J. Granger (2010).  Modelling Nonlinear Economic Time Series.  Oxford: Oxford University Press.
  * Tjøstheim (2012) Tjøstheim, D. (2012).  Some recent theory for autoregressive count time series.  TEST 21, 413–438.
  * Tjøstheim (2015) Tjøstheim, D. (2015).  Count Time Series with Observation-Driven Autoregressive Parameter Dynamics.  In R. Davis, S. Holan, R. Lund, and N. Ravishanker (Eds.), Handbook of Discrete-Valued Time Series, Handbooks of Modern Statistical Methods, pp. 77–100. London: Chapman & Hall.
  * Tsay (2014) Tsay, R. S. (2014).  Multivariate Time Series Analysis.  John Wiley & Sons, Inc., Hoboken, NJ.
  * Veraart (2019) Veraart, A. E. (2019).  Modeling, simulation and inference for multivariate time series of counts using trawl processes.  Journal of Multivariate Analysis 169, 110–129.
  * Wang et al. (2014) Wang, C., H. Liu, J.-F. Yao, R. A. Davis, and W. K. Li (2014).  Self-excited threshold Poisson autoregression.  Journal of the American Statistical Association 109, 777–787.
  * Wang and Wang (2018) Wang, F. and H. Wang (2018).  Modelling non-stationary multivariate time series of counts via common factors.  Journal of the Royal Statistical Society: Series B 80, 769–791.
  * Wedel et al. (2003) Wedel, M., U. Böckenholt, and W. A. Kamakura (2003).  Factor models for multivariate count data.  Journal of Multivariate Analysis 87, 356–369.
  * West (2020) West, M. (2020).  Bayesian forecasting of multivariate time series: scalability, structure uncertainty and decisions.  Annals of the Institute of Statistical Mathematics 72, 1–31.
  * West and Harrison (1997) West, M. and P. Harrison (1997).  Bayesian Forecasting and Dynamic Models (second ed.).  New York: Springer.
  * Woodard et al. (2011) Woodard, D. W., D. S. Matteson, and S. G. Henderson (2011).  Stationarity of count-valued and nonlinear time series models.  Electronic Journal of Statistics 5, 800–828.
  * Wu and Shao (2004) Wu, W. B. and X. Shao (2004).  Limit theorems for iterated random functions.  Journal of Applied Probability 41, 425?436.
  * Yang et al. (2020) Yang, L., E. W. Frees, and Z. Zhang (2020).  Nonparametric estimation of copula regression models with discrete outcomes.  Journal of the American Statistical Association 115, 707–720.
  * Zeger (1988) Zeger, S. L. (1988).  A regression model for time series of counts.  Biometrika 75, 621–629.
  * Zeger and Qaqish (1988) Zeger, S. L. and B. Qaqish (1988).  Markov regression models for time series: a quasi-likelihood approach.  Biometrics 44, 1019–1031.
  * Zhang et al. (2017) Zhang, C., N. Chen, and Z. Li (2017).  State space modeling of autocorrelated multivariate poisson counts.  IISE Transactions 49, 518–531.
  * Zhang and Wu (2017) Zhang, D. and W. B. Wu (2017).  Gaussian approximation for high dimensional time series.  The Annals of Statistics 45, 1895–1919.
  * Zhang et al. (2017) Zhang, Y., H. Zhou, J. Zhou, and W. Sun (2017).  Regression models for multivariate count data.  J. Comput. Graph. Statist. 26, 1–13.