# Conditional-mean Multiplicative Operator Models for Count Time Series

Christian H. Weiß  Helmut Schmidt University, Department of Mathematics and
Statistics, Hamburg, Germany. E-Mail<EMAIL_ADDRESS>ORCID:
0000-0001-8739-6631. Fukang Zhu  School of Mathematics, Jilin University, 2699
Qianjin, Changchun 130012, China. E-Mail<EMAIL_ADDRESS>ORCID:
0000-0002-8808-8179.Corresponding author.

###### Abstract

Multiplicative error models (MEMs) are commonly used for real-valued time
series, but they cannot be applied to discrete-valued count time series as the
involved multiplication would not preserve the integer nature of the data.
Thus, the concept of a multiplicative operator for counts is proposed (as well
as several specific instances thereof), which are then used to develop a kind
of MEMs for count time series (CMEMs). If equipped with a linear conditional
mean, the resulting CMEMs are closely related to the class of so-called
integer-valued generalized autoregressive conditional heteroscedasticity
(INGARCH) models and might be used as a semi-parametric extension thereof.
Important stochastic properties of different types of INGARCH-CMEM as well as
relevant estimation approaches are derived, namely types of quasi-maximum
likelihood and weighted least squares estimation. The performance and
application are demonstrated with simulations as well as with two real-world
data examples.

Key words: Count time series; INGARCH models; multiplicative error model;
multiplicative operator; semi-parametric estimation.

###  1 Introduction

Since their introduction by Engle (2002), the _multiplicative error models_
(MEMs) received considerable research interest in the literature, see
Brownlees et al. (2012), Cipollini and Gallo (2023), and the references
therein. MEMs are defined for a positively real-valued process
$(Y_{t})_{t\in\mathbb{Z}=\\{\ldots,-1,0,1,\ldots\\}}$; let $\mathcal{F}_{s}$
denote the information being available up to time $s$. Then, $(Y_{t})$ follows
a MEM if it satisfies the equation

$Y_{t}\ =\ \mu_{t}\cdot\varepsilon_{t},$ (1.1)

where $\mu_{t}|\mathcal{F}_{t-1}$ is a deterministic, truly positive number,
and where $\varepsilon_{t}$ is a positively real-valued random variable that
is generated independently of $\mathcal{F}_{t-1}$ with (conditional) mean
$E[\varepsilon_{t}\ |\ \mathcal{F}_{t-1}]=E[\varepsilon_{t}]=1$ and variance
$\sigma_{t}^{2}:=V[\varepsilon_{t}\ |\
\mathcal{F}_{t-1}]=V[\varepsilon_{t}]\in(0,\infty)$. The observations
$(Y_{t})$ of the MEM satisfy

$E[Y_{t}\ |\ \mathcal{F}_{t-1}]\ =\ \mu_{t},\qquad V[Y_{t}\ |\
\mathcal{F}_{t-1}]=\sigma_{t}^{2}\,\mu_{t}^{2}.$ (1.2)

If $(\varepsilon_{t})$ are i. i. d. with mean 1 and variance $\sigma^{2}$, and
if the conditional means $(\mu_{t})$ satisfy the equation

$\mu_{t}\ =a_{0}+\sum_{i=1}^{p}a_{i}\,Y_{t-i}+\sum_{j=1}^{q}b_{j}\,\mu_{t-j},$
(1.3)

then the resulting MEM is said to be an autoregressive conditional duration
(ACD) model of order $(p,q)$, see Engle & Russell (1998). Here, the
constraints $a_{0}>0$, $a_{1},\ldots,a_{p},b_{1},\ldots,b_{q}\geq 0$, and
$\sum_{i=1}^{p}a_{i}+\sum_{j=1}^{q}b_{j}<1$ have to be satisfied to ensure the
existence of a truly positive conditional mean. The construction (1.1) is also
closely related to that of the generalized autoregressive conditional
heteroscedasticity (GARCH) models (Bollerslev, 1986).

The conditional-mean equation (1.3) of the ACD model gave rise to Heinen
(2003) to define an integer-valued counterpart for a count process $(X_{t})$,
i. e., where the $X_{t}$ are (unbounded) count random variables with range
$\mathbb{N}_{0}=\\{0,1,\ldots\\}$. With $M_{t}=E[X_{t}\ |\ \mathcal{F}_{t-1}]$
denoting the conditional mean (thus, now, $\mathcal{F}_{t-1}$ is the
$\sigma$-field generated by
$\\{(X_{t-1},M_{t-1}),(X_{t-2},M_{t-2}),\ldots\\}$), Heinen (2003) considered
the model equation

$M_{t}\ =\ a_{0}+\sum_{i=1}^{p}a_{i}\,X_{t-i}+\sum_{j=1}^{q}b_{j}\,M_{t-j},$
(1.4)

where $a_{0}>0,a_{1},\ldots,a_{p}\geq 0,b_{1},\ldots,b_{q}\geq 0$, and
$X_{t}|\mathcal{F}_{t-1}$ is generated by a conditional count distribution
having mean $M_{t}$, such as the Poisson (Poi) distribution
$\textup{Poi}(M_{t})$. Heinen (2003) referred to this model as the
autoregressive conditional Poisson (ACP) model, while Ferland et al. (2006)
used the term integer-valued GARCH (INGARCH) model. Although there is still a
debate about the appropriate terminology for models of the type (1.4), see
Remark 4.1.2 in Weiß (2018), it seems that the name INGARCH model is most
often used in the literature. Thus, in the sequel, we shall also refer to
count models satisfying (1.4) as INGARCH-type models.

In this article, our aim is to define MEM-type models for count time series,
but not starting at (1.3) as done by Heinen (2003), but one step earlier right
with the original definition (1.1) of a MEM. But since the multiplication
$\alpha\cdot\epsilon$ of a count random variable $\epsilon$ by a positive real
number $\alpha\in(0,\infty)$ does usually not lead to a count value anymore
(“multiplication problem”, see Weiß, 2018, p. 16), an integer-valued
substitute of the multiplication is required. This multiplication problem is
well known from the adaption of the autoregressive moving-average (ARMA) model
to the count-data case. The so-called integer-valued ARMA (INARMA) model
modifies the ordinary ARMA recursion by substituting its multiplications by
so-called “thinning operators” (plus several further assumptions, see Weiß
(2018) for details). During the last decades, a huge number of different types
of thinning operators have been proposed in the literature, see Weiß (2008)
and Scotto et al. (2015) for surveys. Generally, a thinning operator “thinn”
is a random operator from $\mathbb{N}_{0}$ to $\mathbb{N}_{0}$, whose
conditional mean satisfies $E[\textup{thinn}(\epsilon)\ |\
\epsilon]\leq\epsilon$. Many common thinning operators are defined in an
additive way and fall within the family of _generalized thinning_ operators
proposed by Latour (1998):

$\alpha\bullet_{\nu}\epsilon\ :=\ \sum_{j=1}^{\epsilon}\ Z_{j},$ (1.5)

where the $Z_{j}$ are i. i. d. count random variables (_counting series_ ,
generated independently of $\epsilon$) with mean $\alpha\in(0,1)$ and variance
$\nu>0$. As a result, $E[\alpha\bullet_{\nu}\epsilon\ |\
\epsilon]=\alpha\cdot\epsilon\ \leq\epsilon$, which justifies its use as an
integer-valued substitute of the multiplication.

###### 1.1 Example

Popular examples of generalized thinning are

  * (i)

_binomial thinning_ dating back to Steutel & van Harn (1979) (and commonly
denoted by “$\alpha\circ$”), where the $Z_{j}$ are Bernoulli random variables,
$Z_{j}\sim\textup{Bin}(1,\alpha)$ (thus $\nu=\nu(\alpha)=\alpha(1-\alpha)$),
such that $\alpha\circ\epsilon|\epsilon\ \sim\textup{Bin}(\epsilon,\alpha)$
with range $\\{0,\ldots,\epsilon\\}$;

  * (ii)

_negative-binomial (NB) thinning_ proposed by Ristić et al. (2009) (and
commonly denoted by “$\alpha*$”), where the $Z_{j}$ are geometrically
distributed with parameter $\pi:=1/(1+\alpha)$, $Z_{j}\sim\textup{NB}(1,\pi)$
(hence $E[Z_{j}]=\alpha$ and $\nu(\alpha)=\alpha(1+\alpha)$), such that
$\alpha*\epsilon|\epsilon\ \sim\textup{NB}\big{(}\epsilon,1/(1+\alpha)\big{)}$
with range $\mathbb{N}_{0}$;

  * (iii)

_Poisson thinning_ (Scotto et al., 2015) with $Z_{j}\sim\textup{Poi}(\alpha)$
(thus $\nu(\alpha)=\alpha$) such that $\alpha\bullet_{\nu}\epsilon|\epsilon\
\sim\textup{Poi}(\epsilon\cdot\alpha)$ with range $\mathbb{N}_{0}$.

In this article, we propose multiplicative operators for counts, which arise,
for example, by generalizing operator (1.5). Section 2 introduces new
multiplicative operators and discusses their basic properties. Using such
multiplicative operator, we adapt the MEM (1.1) to count time series, leading
to a _count MEM_ (CMEM). In particular, we propose CMEMs under the framework
of INGARCH model, referred to as INGARCH-CMEMs. Two specific cases are studied
in Section 3, namely CMEMs based on the compounding operator and the binomial
multiplicative operator, respectively. Stationarity and ergodicity of the new
models are established, and some probabilistic properties are also discussed.
Note that our novel INGARCH-CMEMs differ from the thinning-based MEM adaption
recently proposed by Aknouche and Scotto (2023), see Remark 3.1.1 for details.
In Section 4, two kinds of semi-parametric estimation methods are considered:
quasi-maximum likelihood estimation (QMLE) and two-stage weighted least
squares estimation (2SWLSE). The finite-sample performance of these methods is
investigated by simulations in Section 5. Section 6 then illustrates the
application of our novel CMEMs to two real-world data examples, namely to
count time series referring to infectious diseases and financial transactions,
respectively. Finally, Section 7 concludes the article and outlines topics for
future research.

###  2 Multiplicative Operators for Counts

For being able to adapt (1.1) in full generality, the restriction
$E[\textup{thinn}(\epsilon)\ |\ \epsilon]\leq\epsilon$ being inherent to any
thinning operator is too limiting, also the case “$>\epsilon$” (i. e., also a
“thickening”) should be permitted. Thus, any operator
“$\alpha\odot:\mathbb{N}_{0}\to\mathbb{N}_{0}$” satisfying

$\textstyle E[\alpha\odot\epsilon\ |\ \epsilon]\ =\
\alpha\cdot\epsilon\qquad\text{for all }\alpha>0$ (2.1)

is referred to as an integer-valued _multiplicative operator_. Note that the
multiplicative property refers to the conditional mean of
$\alpha\odot\epsilon$, but to keep it simple, we use the terminology
“multiplicative operator” for “$\alpha\odot$”. By (2.1), also
$E[\alpha\odot\epsilon]=\alpha\,E[\epsilon]$ is satisfied. To compute the
(conditional) variance of $\alpha\odot\epsilon$, further assumptions on
“$\alpha\odot$” are necessary. In any of the subsequent special cases, also
the conditional variance of $\alpha\odot\epsilon|\epsilon$ is proportional to
$\epsilon$, with the constant of proportionality given by some
$\nu=\nu(\alpha)>0$. In this special case, we immediately obtain

$V[\alpha\odot\epsilon\ |\ \epsilon]\ =\ \nu\cdot\epsilon,\qquad
V[\alpha\odot\epsilon]\ =\ \nu\,E[\epsilon]+\alpha^{2}\,V[\epsilon].$ (2.2)

While moment calculations can be done solely based on (2.1) and (2.2),
likelihood calculations or the simulation of “$\alpha\odot$” require to fully
specify the conditional distribution of $\alpha\odot\epsilon|\epsilon$. This
is done in the subsequent special cases.

The requirements imposed by (2.1) are easily achieved with an operator of the
form (1.5), namely by requiring that $\alpha$ can take any value in
$(0,\infty)$. To not confuse the terminology, we refer to

$\begin{array}[]{@{}l}\alpha\bullet_{\nu}\epsilon\ :=\ \sum_{j=1}^{\epsilon}\
Z_{j}\qquad\text{with i.\,i.\,d.\ counting series }(Z_{j}),\\\\[4.30554pt]
\text{which has arbitrary mean~{}$\alpha>0$ and
variance~{}$\nu>0$},\end{array}$ (2.3)

as a _compounding operator_. Important (conditional) moments of
“$\alpha\bullet_{\nu}$” are

$\begin{array}[]{@{}l@{\qquad}l}E[\alpha\bullet_{\nu}\epsilon\ |\ \epsilon]\
=\ \alpha\cdot\epsilon,&V[\alpha\bullet_{\nu}\epsilon\ |\ \epsilon]\ =\
\nu\cdot\epsilon,\qquad E\big{[}(\alpha\bullet_{\nu}\epsilon)^{2}\ \big{|}\
\epsilon\big{]}\ =\ \nu\,\epsilon+\alpha^{2}\,\epsilon^{2},\\\\[4.30554pt]
E[\alpha\bullet_{\nu}\epsilon]\ =\
\alpha\,E[\epsilon],&V[\alpha\bullet_{\nu}\epsilon]\ =\
\nu\,E[\epsilon]+\alpha^{2}\,V[\epsilon],\end{array}$ (2.4)

see (2.2). Also the probability generating function (pgf) of (2.3) is easily
expressed. Denoting the pgf of the counting series by
$\textup{pgf}_{Z}(u):=E[u^{Z}]$, and that of $\epsilon$ by
$\textup{pgf}_{\epsilon}(u)$, it holds that

$\textup{pgf}_{\alpha\bullet_{\nu}\epsilon}(u)\ =\
\textup{pgf}_{\epsilon}\big{(}\textup{pgf}_{Z}(u)\big{)}.$ (2.5)

Thus, the distribution of $\alpha\bullet_{\nu}\epsilon$ can be expressed as
the convolution

$P(\alpha\bullet_{\nu}\epsilon=k)\ =\ \sum_{l=0}^{\infty}P(\epsilon=l)\cdot
P(\alpha\bullet_{\nu}l=k).$ (2.6)

Obviously, a compounding operator with Poisson or geometric counting series
(recall Example 1.1) is still well-defined, with

$\alpha\bullet_{\nu}l\ \sim\
\left\\{\begin{array}[]{ll}\textup{Poi}(l\,\alpha)&\text{if
}Z_{j}\sim\textup{Poi}(\alpha),\\\\[2.15277pt]
\textup{NB}\big{(}l,1/(1+\alpha)\big{)}&\text{if
}Z_{j}\sim\textup{NB}\big{(}1,1/(1+\alpha)\big{)}.\\\ \end{array}\right.$

However, a Bernoulli counting series with $\alpha>1$ does not exist, and the
extension of the binomial thinning operator “$\alpha\circ$” to a full
multiplicative operator is not obvious. In what follows, we propose a
generalized definition of the “$\alpha\circ$” operator, which contains both
thinning (if $\alpha<1$) and thickening (if $\alpha>1$):

$\alpha\otimes\epsilon=\ \lfloor\alpha\rfloor\cdot\epsilon\ +\
\textup{Bin}\big{(}\epsilon,\ \alpha-\lfloor\alpha\rfloor\big{)}.$ (2.7)

Here, $\lfloor\alpha\rfloor$ denotes the greatest integer less than or equal
to $\alpha$, $\textup{Bin}\big{(}\epsilon,\
\alpha-\lfloor\alpha\rfloor\big{)}$ expresses a count random variable being
generated according to this (conditional) distribution, and the two terms in
the above addition are conditionally independent. The operator $\otimes$ in
(2.7) is referred to as the _binomial multiplicative operator_ , and it
differs from the (purely deterministic and non-multiplicative) rounding
operator of Kachour and Yao (2009). Note that if $\alpha>1$, it is guaranteed
that $\alpha\otimes\epsilon$ gives a value $\geq\epsilon$, which differs from
the compounding operator $\alpha\bullet_{\nu}$ with $\alpha>1$. $\alpha=1$
corresponds to a (non-random) preservation of $\epsilon$. Further stochastic
properties of the binomial multiplicative operator are as follows.

###### 2.1 Proposition

Operator (2.7) has the conditional range
$\big{\\{}\lfloor\alpha\rfloor\,\epsilon,\ldots,(\lfloor\alpha\rfloor+1)\,\epsilon\big{\\}}$,
the conditional moments

$E[\alpha\otimes\epsilon\ |\ \epsilon]\ =\ \alpha\cdot\epsilon,\qquad
V[\alpha\otimes\epsilon\ |\ \epsilon]\ =\
(\alpha-\lfloor\alpha\rfloor)(1-\alpha+\lfloor\alpha\rfloor)\,\epsilon\
=:\nu(\alpha)\,\epsilon,$

the unconditional moments

$E[\alpha\otimes\epsilon]\ =\ \alpha\cdot E[\epsilon],\qquad
V[\alpha\otimes\epsilon]\ =\ \alpha^{2}\cdot V[\epsilon]\ +\
(\alpha-\lfloor\alpha\rfloor)(1-\alpha+\lfloor\alpha\rfloor)\cdot
E[\epsilon],$

and the pgf

$\textup{pgf}_{\alpha\otimes\epsilon}(u)\ =\
\textup{pgf}_{\epsilon}\Big{(}u^{\lfloor\alpha\rfloor}\,\big{(}1-(\alpha-\lfloor\alpha\rfloor)+(\alpha-\lfloor\alpha\rfloor)\,u\big{)}\Big{)}.$

###### Proof.

The conditional mean of the operator is

$\displaystyle E[\alpha\otimes\epsilon\ |\ \epsilon]$ $\displaystyle=$
$\displaystyle E\big{[}\lfloor\alpha\rfloor\,\epsilon\ |\ \epsilon]\ +\
E\big{[}\textup{Bin}\big{(}\epsilon,\ \alpha-\lfloor\alpha\rfloor\big{)}\ |\
\epsilon\big{]}$ $\displaystyle=$
$\displaystyle\lfloor\alpha\rfloor\,\epsilon\ +\
(\alpha-\lfloor\alpha\rfloor)\,\epsilon\ =\ \alpha\cdot\epsilon.$

The conditional variance is

$\displaystyle V[\alpha\otimes\epsilon\ |\ \epsilon]$ $\displaystyle=$
$\displaystyle
V\big{[}\underbrace{\lfloor\alpha\rfloor\,\epsilon}_{\text{(*)}}\ +\
\textup{Bin}\big{(}\epsilon,\ \alpha-\lfloor\alpha\rfloor\big{)}\ |\
\epsilon\big{]}$ $\displaystyle=$ $\displaystyle
V\big{[}\textup{Bin}\big{(}\epsilon,\ \alpha-\lfloor\alpha\rfloor\big{)}\ |\
\epsilon\big{]}\ =\
(\alpha-\lfloor\alpha\rfloor)(1-\alpha+\lfloor\alpha\rfloor)\,\epsilon,$

where at (*), we used that $\lfloor\alpha\rfloor\,\epsilon\big{|}\epsilon$ is
conditionally constant. So the unconditional moments immediately follow from
(2.2). The pgf is

$\displaystyle\textup{pgf}_{\alpha\otimes\epsilon}(u)$ $\displaystyle=$
$\displaystyle E\big{[}u^{\alpha\otimes\epsilon}\big{]}\ =\
E\Big{[}E\big{[}u^{\lfloor\alpha\rfloor\cdot\epsilon\ +\
\textup{Bin}(\epsilon,\ \alpha-\lfloor\alpha\rfloor)}\ \big{|}\
\epsilon\big{]}\Big{]}$ $\displaystyle=$ $\displaystyle
E\Big{[}u^{\lfloor\alpha\rfloor\cdot\epsilon}\cdot
E\big{[}u^{\textup{Bin}(\epsilon,\ \alpha-\lfloor\alpha\rfloor)}\ \big{|}\
\epsilon\big{]}\Big{]}$ $\displaystyle=$ $\displaystyle
E\Big{[}u^{\lfloor\alpha\rfloor\cdot\epsilon}\cdot\big{(}1-(\alpha-\lfloor\alpha\rfloor)+(\alpha-\lfloor\alpha\rfloor)\,u\big{)}^{\epsilon}\Big{]}$
$\displaystyle=$
$\displaystyle\textup{pgf}_{\epsilon}\Big{(}u^{\lfloor\alpha\rfloor}\,\big{(}1-(\alpha-\lfloor\alpha\rfloor)+(\alpha-\lfloor\alpha\rfloor)\,u\big{)}\Big{)},$

which completes the proof of Proposition 2.1. ∎

While our main interest is in CMEMs based on the compounding operator (2.3)
and on the binomial multiplicative operator (2.7), respectively, let us
conclude Section 2 by pointing out that multiplicative operators do not need
to be of such an additive structure. Actually, any count distribution with
closed-form mean might give rise to a multiplicative operator. For
illustration, let us take the zero-inflated Poisson (ZIP) distribution
$\textup{ZIP}(\lambda,\omega)$ with $\lambda>0$ and $\omega\in(0,1)$ as the
starting point. Its mean $\mu$ and variance $\sigma^{2}$ are given by the
relations $\mu=(1-\omega)\,\lambda$ and $\sigma^{2}/\mu=\lambda-\mu+1$. Thus,
a ZIP-multiplicative operator “$\alpha\varoast_{\kappa}$” with $\alpha>0$ and
$\kappa>1$ might be defined by requiring that
$\alpha\varoast_{\kappa}\epsilon|\epsilon$ is conditionally ZIP-distributed
with mean $E[\alpha\varoast_{\kappa}\epsilon\,|\,\epsilon]=\alpha\,\epsilon$
and variance $V[\alpha\varoast_{\kappa}\epsilon\,|\,\epsilon]=\kappa\cdot
E[\alpha\varoast_{\kappa}\epsilon\,|\,\epsilon]$, i. e., with conditional ZIP-
parameters $\lambda=\alpha\,\epsilon+\kappa-1$ and
$\omega=(\kappa-1)/(\alpha\,\epsilon+\kappa-1)$. Note that (2.2) also applies
to “$\alpha\varoast_{\kappa}$”.

###  3 MEMs for Count Time Series

Inspired by Equation (1.1), a CMEM for a process $(X_{t})$ consisting of
unbounded counts is defined by

$X_{t}\ =\ M_{t}\odot\epsilon_{t},$ (3.1)

where $M_{t}|\mathcal{F}_{t-1}$ is a deterministic truly positive number.
$(\epsilon_{t})$ is i. i. d., and $\epsilon_{t}$ is a count random variable
that is generated independently of $\mathcal{F}_{t-1}$ with mean
$E[\epsilon_{t}]=1$ and variance $V[\epsilon_{t}]=\sigma^{2}$. The
multiplicative operator “$\odot$” defined in (2.1) is executed independently
of $\mathcal{F}_{t-1}$ and $\epsilon_{t}$.

With these assumptions, the conditional mean of (3.1) satisfies

$\begin{array}[]{@{}rl}E[X_{t}\ |\ \mathcal{F}_{t-1}]\
=&E\big{[}E[M_{t}\odot\epsilon_{t}\ |\ \epsilon_{t},\mathcal{F}_{t-1}]\
\big{|}\ \mathcal{F}_{t-1}\big{]}\\\\[4.30554pt] =&E[M_{t}\cdot\epsilon_{t}\
|\ \mathcal{F}_{t-1}]\ =\ M_{t}\cdot 1\ =\ M_{t}.\end{array}$ (3.2)

Thus, if the conditional mean $M_{t}$ follows a recursive scheme of the form

$M_{t}\ =\ f(X_{t-1},\ldots,X_{t-p},\ M_{t-1},\ldots,M_{t-q})\quad\text{with
}p\geq 1,\ q\geq 0,$ (3.3)

then the CMEM becomes an INGARCH-type model of order $(p,q)$, denoted as
_INGARCH $(p,q)$-CMEM_. The most well-known special case of (3.3) is the
linear specification (1.4). But also non-linear mean specifications are
possible for $M_{t}$, e. g., by adapting the log-linear structure discussed by
Fokianos & Tjøstheim (2011) or the approximately linear softplus structure
proposed by Weiß et al. (2022), both allowing for negative model parameters
(also see Section 5).

###### 3.1 Remark

One of the formal advantages of the proposed model is that it is defined
through a stochastic equation and not by a conditional distribution. In
particular, the standard Poi-INGARCH model of Ferland et al. (2006) is nested
in this model, see Section 3.1 for details, and is thus defined by an equation
instead of a conditional distribution.

###### 3.2 Remark

We have already mentioned some specific operators before. An important class
of operators defined by (2.1) is given by a hierarchical mixture of INGARCHs.
We already indicated that the negative binomial INGARCH is a specific case of
our model. Indeed, the mixed Poisson autoregression
$\\{Y_{t},t\in\mathbb{Z}\\}$ given by

$Y_{t}=N_{t}(\varepsilon_{t}M_{t}),~{}~{}M_{t}=h(Y_{t-1},M_{t-1}),~{}~{}t\in\mathbb{Z},$

is a specific case of (3.1), where $(N_{t})$ is an independent sequence of
homogeneous Poisson processes with unit intensity. More generally, any
hierarchical mixture process given by

$\displaystyle Y_{t}|\mathcal{F}_{t-1},\varepsilon_{t}$ $\displaystyle\sim
F_{M_{t},\varepsilon_{t}},$ $\displaystyle\varepsilon_{t}$ $\displaystyle\sim
F_{\varepsilon_{t}},$

($F_{M_{t},\varepsilon_{t}}$ and $F_{\varepsilon_{t}}$ being given
distributions) is a particular instance of model (2.1). So

$f_{Y_{t}|\mathcal{F}_{t-1}}(y)=\sum_{\varepsilon}f_{M_{t},\varepsilon_{t}}(y,\varepsilon)f_{\varepsilon_{t}}(\varepsilon)$

is a finite or infinite mixture of distributions over $\varepsilon_{t}$
depending on whether the range of $\varepsilon_{t}$ is finite or infinite. In
this sense, the proposed model (3.1) is deeply related to a mixture of INGARCH
representations, see Zhu et al. (2010), Mao et al. (2020), Doukhan et al.
(2021), among others.

In the special case, where the multiplicative operator “$\odot$” also
satisfies $V[\alpha\odot\epsilon\,|\,\epsilon]=\nu(\alpha)\cdot\epsilon$, we
conclude from (2.2) that

$\begin{array}[]{@{}rl}V[X_{t}\ |\ \mathcal{F}_{t-1}]\
=&E\big{[}V[M_{t}\odot\epsilon_{t}\ |\ \epsilon_{t},\mathcal{F}_{t-1}]\
\big{|}\ \mathcal{F}_{t-1}\big{]}+V\big{[}E[M_{t}\odot\epsilon_{t}\ |\
\epsilon_{t},\mathcal{F}_{t-1}]\ \big{|}\
\mathcal{F}_{t-1}\big{]}\\\\[4.30554pt] =&E[\nu(M_{t})\cdot\epsilon_{t}\ |\
\mathcal{F}_{t-1}]+V[M_{t}\cdot\epsilon_{t}\ |\ \mathcal{F}_{t-1}]\ =\
\nu(M_{t})+\sigma^{2}M_{t}^{2}.\end{array}$ (3.4)

So the conditional variance is caused by two sources: the variance of the
operator “$\odot$” and the one of the innovation term $\epsilon_{t}$.

To establish the existence of an INGARCH$(p,q)$-CMEM, one can benefit from the
existing INGARCH literature. To apply Theorem 2.1 of Doukhan and Neumann
(2019), for example, one has to ensure that on the one hand, assumptions
(A$1^{\prime}$) and (A2) of Doukhan and Neumann (2019) hold for the given
INGARCH specification (3.3), and on the other hand, certain distributional
requirements for the multiplicative operator “$\odot$” have to be satisfied.
More precisely, if also the conditional distribution of
$X_{t}|\mathcal{F}_{t-1}=M_{t}\odot\epsilon_{t}|M_{t}=:Q(M_{t})$ satisfies
(A3) of Doukhan and Neumann (2019), then their Theorem 2.1 can be used to
conclude on the existence of a stationary solution to (3.1) and (3.3), as well
as on mixing properties. This, certainly, requires to first fix a specific
choice of the multiplicative operator “$\odot$”. In what follows, we discuss
two kinds of CMEMs (3.1) with the linear specification (1.4), namely based on
the compounding operator (2.3) and on the binomial multiplicative operator
(2.7), respectively.

####  3.1 CMEMs Using a Compounding Operator

In this section, let us consider the special case where the CMEM (3.1) uses
the compounding operator (2.3) $\alpha\bullet_{\nu}$ with arbitrary
$\alpha,\nu>0$ as the multiplicative operator, i. e.,

$X_{t}\ =\ M_{t}\bullet_{\nu}\epsilon_{t},$ (3.5)

where $M_{t}$ is given in (1.4) and $\epsilon_{t}$ is the same one defined in
(3.1). Note that $X_{t}$ can be written as $\sum_{j=1}^{\epsilon_{t}}Z_{tj}$,
where $(Z_{tj})$ are count random variables with $E[Z_{tj}\ |\
\mathcal{F}_{t-1}]=M_{t}$, and $(\epsilon_{t})$ are independent of $(Z_{tj})$.
Given $\mathcal{F}_{t-1},Z_{sj}$ are conditionally independent for $s\geq t$
and different $j$. Furthermore, if $t$ is fixed, then the $Z_{tj}$ are
conditionally i. i. d. with some $Z_{t}$ for different $j$, and more
discussions about $Z_{t}$ will be given later in the context of Theorem 3.1.6.

The formulae (2.4)–(2.5) are now used to derive properties of the compounding
CMEM (3.5). While still $E[X_{t}\ |\ \mathcal{F}_{t-1}]=M_{t}$ holds, we get

$\begin{array}[]{@{}rl}E[X_{t}^{2}\ |\ \mathcal{F}_{t-1}]\
=&E\big{[}E[(M_{t}\bullet_{\nu}\epsilon_{t})^{2}\ |\
\epsilon_{t},\mathcal{F}_{t-1}]\ \big{|}\
\mathcal{F}_{t-1}\big{]}\\\\[4.30554pt]
\overset{\eqref{CompoundingCMom}}{=}&E\big{[}\nu(M_{t})\,\epsilon_{t}+M_{t}^{2}\,\epsilon_{t}^{2}\
\big{|}\ \mathcal{F}_{t-1}\big{]}\ =\
\nu(M_{t})+(\sigma^{2}+1)M_{t}^{2},\end{array}$

where $\nu=\nu(M_{t})$. Thus,

$V[X_{t}\ |\ \mathcal{F}_{t-1}]\ =\ \nu(M_{t})+\sigma^{2}M_{t}^{2},$ (3.6)

also see (3.4). In the special case of a Poi-counting series, we have
$\nu(M_{t})=M_{t}$, whereas a NB-counting series leads to
$\nu(M_{t})=M_{t}(1+M_{t})$, recall Example 1.1. Compared with an ordinary
INGARCH model with conditional mean $M_{t}$, the main difference in
conditional mean and variance between the CMEM (3.5) and the INGARCH model is
as follows: while having the same conditional mean $M_{t}$, there is an
additional term $\sigma^{2}M_{t}^{2}$ in the conditional variance of the CMEM,
which makes the CMEM more flexible for applications in practice. Only in the
boundary case where $\epsilon_{t}=1$ deterministically (thus $\sigma^{2}=0$),
the CMEM specification reduces to the ordinary INGARCH approach. In
particular, the special case $\epsilon_{t}\equiv 1$ and Poi-counting series
just agrees with the ordinary Poi-INGARCH model of Ferland et al. (2006).

###### 3.1.1 Remark

Aknouche and Scotto (2023) proposed a multiplicative thinning-based INGARCH
model defined by

$X_{t}=\lambda_{t}\cdot\epsilon_{t},~{}~{}~{}\lambda_{t}=1+\omega\circ
m+\sum_{i=1}^{q}\alpha_{i}\circ
Y_{t-i}+\sum_{j=1}^{p}\beta_{j}\circ\lambda_{t-j},$ (3.7)

where $\epsilon_{t}$ is given like in (3.1), and where $\circ$ is the binomial
thinning operator discussed in Example 1.1 (i). Model (3.7) has a similar mean
structure as model (3.5), and the inclusion of “$\omega\circ$” makes model
(3.7) analogous to (3.5). But there are also substantial differences between
both approaches. First, the positive integer $m$ in (3.7) is not treated as a
parameter (for practice, Aknouche and Scotto (2023) suggest to choose $m$ as
the integer part of the sample mean). Second, the conditional variance of
model (3.7) is a function of two recursive items, see equation (2.11) in
Aknouche and Scotto (2023), while the conditional variance of our model (3.1)
is just a function of $M_{t}$, see (3.6). The complexity of the conditional
variance also affects the estimation method discussed in Section 4.

The conditional pgf of the compounding CMEM (3.5) satisfies

$\textup{pgf}_{X_{t}|\mathcal{F}_{t-1}}(u)\ =\
E[u^{M_{t}\bullet_{\nu}\epsilon_{t}}\ |\ \mathcal{F}_{t-1}]\
\overset{\eqref{CompoundingPgf}}{=}\
\textup{pgf}_{\epsilon}\big{(}\textup{pgf}_{Z|\mathcal{F}_{t-1}}(u)\big{)},$
(3.8)

where $Z|\mathcal{F}_{t-1}$ refers to the counting series of the operator
“$M_{t}\bullet_{\nu}$” given $\mathcal{F}_{t-1}$, recall (3.5).

###### 3.1.2 Example

In the special case of Poisson-distributed innovations,
$\epsilon_{t}\sim\textup{Poi}(1)$, (3.8) implies a conditional compound-
Poisson (CP) distribution, namely
$\textup{pgf}_{X_{t}|\mathcal{F}_{t-1}}(u)=\exp\big{(}\textup{pgf}_{Z|\mathcal{F}_{t-1}}(u)-1\big{)}$.
Thus, if $M_{t}$ also satisfies (1.4), we get a type of CP-INGARCH model as
discussed by Gonçalves et al. (2015a). If, in addition, we have a Poisson-
counting series $Z|\mathcal{F}_{t-1}\sim\textup{Poi}(M_{t})$ with
$\textup{pgf}_{X_{t}|\mathcal{F}_{t-1}}(u)=\exp\big{(}M_{t}\,(u-1)\big{)}$,
(3.8) implies
$\textup{pgf}_{X_{t}|\mathcal{F}_{t-1}}(u)=\exp\big{(}\exp\big{(}M_{t}\,(u-1)\big{)}-1\big{)}$,
i. e., we get a Neyman type-A (NTA) INGARCH model as discussed by Gonçalves et
al. (2015b).

Another special case is
$\epsilon_{t}\sim\textup{ZIP}(\frac{1}{1-\omega},\omega)$ with
$\omega\in(0,1)$, i. e., a ZIP distribution with mean 1 and
$\textup{pgf}_{\epsilon}(u)=\omega+(1-\omega)\,\exp(\frac{u-1}{1-\omega})$,
which leads to a zero-inflated CP-INGARCH model as discussed by Gonçalves et
al. (2016). Generally, there are two possible reasons for the compounding
operator $\alpha\bullet_{\nu}\epsilon=\sum_{j=1}^{\epsilon}Z_{j}$ from (2.3)
to become zero: either $\epsilon=0$, or $\epsilon>0$ but all
$Z_{1}=\ldots=Z_{\epsilon}=0$, i. e.,

$\textstyle P(\alpha\bullet_{\nu}\epsilon=0)\ =\ P(\epsilon=0)\ +\
\sum_{i=1}^{\infty}P(\epsilon=i)\,P(Z=0)^{i}.$

Thus, if $\epsilon_{t}\sim\textup{ZIP}(\frac{1}{1-\omega},\omega)$, we get

$\begin{array}[]{@{}rl}P(X_{t}=0\ |\ \mathcal{F}_{t-1})\ =&\omega\ +\
(1-\omega)\,{\rm
e}^{-1/(1-\omega)}\,\sum\limits_{i=0}^{\infty}\frac{(1-\omega)^{-i}}{i!}\,\textup{pgf}_{Z|\mathcal{F}_{t-1}}(0)^{i}\\\\[3.87498pt]
\ =&\omega\ +\
(1-\omega)\,\exp\Big{(}\frac{\textup{pgf}_{Z|\mathcal{F}_{t-1}}(0)-1}{1-\omega}\Big{)},\end{array}$

in accordance with the expression for
$\textup{pgf}_{X_{t}|\mathcal{F}_{t-1}}(0)$.

The following proposition gives a necessary and sufficient condition for the
first-order stationarity of the model introduced in (3.5), which is analogous
to the corresponding property of Ferland et al. (2006). The proof of
Proposition 3.1.3 is done with the same arguments used in Proposition 1 of
Ferland et al. (2006), so we omit the details.

###### 3.1.3 Proposition

The process $(X_{t})$ defined in (3.5) is first-order stationary if and only
if $\sum_{i=1}^{p}a_{i}+\sum_{j=1}^{q}b_{j}<1$.

As the next step, we discuss the second-order stationarity. We make the
following assumption

$\frac{V[X_{t}\ |\ \mathcal{F}_{t-1}]}{E[X_{t}\ |\
\mathcal{F}_{t-1}]}=v_{0}+v_{1}M_{t},$ (3.9)

with $v_{0}\geq 0$ and $v_{1}\geq 0$ not being simultaneously zero. It is
obvious that the Poisson- and NB-counting series satisfy the above assumption,
recall (3.6), leading to $v_{1}=\sigma^{2}$ and $v_{1}=1+\sigma^{2}$,
respectively. Then, the results of Section 3.2 in Gonçalves et al. (2015a)
hold for the model (3.5). In particular, we have the following proposition.

###### 3.1.4 Proposition

Consider a first-order stationary INGARCH$(1,1)$-CMEM defined according to
(3.5), which also satisfies (3.9). A necessary and sufficient condition for
second-order stationarity is $(a_{1}+b_{1})^{2}+v_{1}a_{1}^{2}<1$.

For second-order stationarity of the INGARCH$(p,q)$-CMEM model with general
order $(p,q)$, we refer to Theorem 4 in Gonçalves et al. (2015a), which will
be not repeated here.

If $\sum_{i=1}^{p}a_{i}+\sum_{j=1}^{q}b_{j}<1$, we obtain the unconditional
mean of $X_{t}$ as

$\mu:=E[X_{t}]=\frac{a_{0}}{1-\sum_{i=1}^{p}a_{i}-\sum_{j=1}^{q}b_{j}}.$
(3.10)

Next, we derive the unconditional variance of $X_{t}$. From Theorem 1 in Weiß
(2009), we know that the autocovariances $\gamma_{X}(k):={\rm
Cov}[X_{t},X_{t-k}]$ and $\gamma_{M}(k):={\rm Cov}[M_{t},M_{t-k}]$ satisfy the
linear equations

$\begin{array}[]{@{}l}\displaystyle\gamma_{X}(k)=\sum_{i=1}^{p}a_{i}\gamma_{X}(|k-i|)+\sum_{j=1}^{\min(k-1,q)}b_{j}\gamma_{X}(k-j)+\sum_{j=k}^{q}b_{j}\gamma_{M}(j-k),\\\
\\\\[-8.61108pt]
\displaystyle\gamma_{M}(k)=\sum_{i=1}^{\min(k,p)}a_{i}\gamma_{M}(|k-i|)+\sum_{i=k+1}^{p}a_{i}\gamma_{X}(i-k)+\sum_{j=1}^{q}b_{j}\gamma_{M}(|k-j|),\end{array}$
(3.11)

for $k\geq 1$ and $k\geq 0$, respectively. Then, using (3.6), the variance can
be obtained from

$V[X_{t}]=E[\nu(M_{t})+\sigma^{2}M_{t}^{2}]+V[M_{t}]\ =\
E[\nu(M_{t})]+\mu^{2}\sigma^{2}+(\sigma^{2}+1)\,V[M_{t}],$ (3.12)

where
$V[M_{t}]=\sum_{i=1}^{p}a_{i}\gamma_{X}(i)+\sum_{j=1}^{q}b_{j}\gamma_{M}(j)$.
For the Poisson-counting case, we have $\nu(M_{t})=M_{t}$, thus

$V[X_{t}]=\mu+\mu^{2}\sigma^{2}+(\sigma^{2}+1)\,V[M_{t}],$ (3.13)

while for the NB-counting case with $\nu(M_{t})=M_{t}(1+M_{t})$, we have

$V[X_{t}]=E[M_{t}+M_{t}^{2}]+\mu^{2}\sigma^{2}+(\sigma^{2}+1)\,V[M_{t}]=\mu+\mu^{2}(\sigma^{2}+1)+(\sigma^{2}+2)\,V[M_{t}].$
(3.14)

Let us illustrate the required computations for the case $p=q=1$, which is
most common in applications.

###### 3.1.5 Example

If $p=q=1$, we get from (3.11), after few calculations, that

$\big{(}1-(a_{1}+b_{1})^{2}+a_{1}^{2}\big{)}\,V[M_{t}]=a_{1}^{2}\,V[X_{t}].$

This can be used to obtain $V[X_{t}]$ after substituting $V[M_{t}]$ into
(3.13) or (3.14), respectively. The autocorrelation function (ACF)
$\rho(k)=Corr[X_{t},X_{t-k}]$ takes the usual form (see Weiß, 2009)

$\rho(k)=(a_{1}+b_{1})^{k-1}\,\frac{a_{1}\,(1-b_{1}(a_{1}+b_{1}))}{1-(a_{1}+b_{1})^{2}+a_{1}^{2}}.$

Together with (3.10), it can be utilized for deriving method-of-moment
estimators.

In practice, equations (3.11) can be used together with the sample ACF for the
identification of the model order $(p,q)$ as well as for getting initial
parameter estimates (method of moments). This is discussed in some more
details in Sections 4.3 and 6.

Before discussing the stationary properties of $(X_{t})$, we first recall the
definition of the CP-INGARCH model given in Gonçalves et al. (2015a). The
count process $(Z_{t})$ is said to satisfy a CP-INGARCH$(p,q)$ model if the
conditional characteristic function of $Z_{t}\ |\ \mathcal{F}_{t-1}$ is
$\Phi_{Z_{t}|\mathcal{F}_{t-1}}(u)=\exp\\{{\rm
i}M_{t}^{*}[\varphi_{t}(u)-1]/\varphi^{\prime}_{t}(0)\\}$ with
$u\in\mathbb{R}$ and $t\in\mathbb{Z}$, and if $M_{t}^{*}$ is defined like in
(1.4), i. e.,
$M_{t}^{*}=a_{0}+\sum_{i=1}^{p}a_{i}Z_{t-i}+\sum_{j=1}^{q}b_{j}M_{t-j}^{*}$.
Here, $(\varphi_{t})$ is a family of characteristic functions on $\mathbb{R}$,
and “i” denotes the imaginary unit.

The following theorem establishes the stationarity and ergodicity of the
compounding CMEM $(X_{t})$ if the counting series arises from the CP-family.
Recall from (3.5) that $X_{t}=\sum_{j=1}^{\epsilon_{t}}Z_{tj}$, where
$(Z_{tj})$ is supposed to be identically distributed with some $(Z_{t})$ that
satisfies $E[Z_{t}\ |\ \mathcal{F}_{t-1}]=E[X_{t}\ |\
\mathcal{F}_{t-1}]=M_{t}$.

###### 3.1.6 Theorem

Assume that $(Z_{t})$ satisfies a CP-INGARCH process. If $\varphi_{t}$ is
deterministic and independent of $t$, and if
$\sum_{i=1}^{p}a_{i}+\sum_{j=1}^{q}b_{j}<1$ holds, then the process $(X_{t})$
defined according to the compounding CMEM (3.5) is stationary and ergodic.

###### Proof.

Using arguments in Gonçalves et al. (2015a), we know that $M_{t}$ can be
written as
$M_{t}=a_{0}/\sum_{j=1}^{q}b_{j}+\sum_{k=1}^{\infty}\psi_{k}X_{t-k}$, where
the coefficient $\psi_{k}$ is given by the power expansion of a rational
function in the neighbourhood of 0. Based on the construction method for
proving stationarity and ergodicity in Ferland et al. (2006) and Gonçalves et
al. (2015a), let
$M_{t}^{(n)}=a_{0}/\sum_{j=1}^{q}b_{j}+\sum_{k=1}^{n}\psi_{k}X_{t-k}$,
$M_{t}^{\ast(n)}=a_{0}/\sum_{j=1}^{q}b_{j}+\sum_{k=1}^{n}\psi_{k}Z_{t-k}$,
$r^{(n)}=E\big{|}(M_{t}^{(n)}-M_{t}^{\ast(n)})-(M_{t}^{(n-1)}-M_{t}^{\ast(n-1)})\big{|}$.
Note that $0\leq
r^{(n)}=|\psi_{n}|\,E|X_{t-n}-Z_{t-n}|=|\psi_{n}|\,\big{(}E[(X_{t-n}-Z_{t-n})\,\mathds{1}(A_{n})]+E[(Z_{t-n}-X_{t-n})\,\mathds{1}(A_{n}^{c})]\big{)}\leq|\psi_{n}|\,\big{(}E(X_{t-n}-Z_{t-n})+E(Z_{t-n}-X_{t-n})\big{)}=0$,
where $A_{n}=\\{X_{t-n}-Z_{t-n}>0\\}$ and $\mathds{1}(\cdot)$ denotes the
indicator function. Thus, $r^{(n)}=0$. Therefore,
$(M_{t}^{(n)}-M_{t}^{\ast(n)})$ is a Cauchy sequence and thus
$(M_{t}^{(n)}-M_{t}^{\ast(n)})$ converges in $L^{1}$ and a.s., i. e.,
$\lim_{n\rightarrow\infty}(M_{t}^{(n)}-M_{t}^{\ast(n)})=0$ a.s. Using Theorem
5 in Gonçalves et al. (2015a), we know that $(Z_{t})$ is stationary and
ergodic.

Notice that $X_{t}=\sum_{j=1}^{\epsilon_{t}}Z_{tj}$, and recall that
$(Z_{tj})$ is supposed to be identically distributed with $(Z_{t})$. Then,
there exists a real function $\xi$, which does not depend on $t$, such that
$X_{t}=\xi\big{(}(Z_{tj},\epsilon_{t}),t\in\mathbb{Z},j\in\mathbb{N}\big{)}$.
Note that $\epsilon_{t}$ is i. i. d. and independent of $\\{Z_{tj}\\}$, so the
result follows by applying Theorem 36.4 in Billingsley (1995). ∎

We conclude Section 3.1 by pointing out that if the distribution of
$\epsilon_{t}$ is not further specified, then the resulting INGARCH-CMEM can
be understood as a kind of semi-parametric INGARCH model, in analogy to the
semi-parametric INAR models proposed by Drost et al. (2009) and Liu et al.
(2021). Then, the conditional distribution of $X_{t}|\mathcal{F}_{t-1}$ is a
convolution between the distribution of counting series and the one of
$\epsilon_{t}$, recall (2.6). Numerically exact computations are possible if
the range of $\epsilon_{t}$ is assumed to be bounded, because then, (2.6)
becomes a finite sum. Generally, denoting $p_{\epsilon,l}=P(\epsilon_{t}=l)$,
the distribution of $\epsilon_{t}$ is completely determined by
$p_{\epsilon,2},p_{\epsilon,3},\ldots$ with the constraint $\sum_{l\geq
2}l\cdot p_{\epsilon,l}<1$, where $p_{\epsilon,0}=\sum_{l\geq 2}(l-1)\cdot
p_{\epsilon,l}$ and $p_{\epsilon,1}=1-p_{\epsilon,0}-\sum_{l\geq
2}p_{\epsilon,l}$ because of the requirement $E[\epsilon_{t}]=1$.

####  3.2 CMEMs Using a Binomial Multiplicative Operator

In this section, we consider the special case where the CMEM (3.1) uses the
binomial multiplicative operator $\otimes$ defined in (2.7) as the
multiplicative operator, i. e.,

$X_{t}\ =\ M_{t}\otimes\epsilon_{t},$ (3.15)

where $M_{t}$ is given in (1.4), and where $\epsilon_{t}$ is defined like in
(3.1). Note that $X_{t}$ can be expressed as $\lfloor
M_{t}\rfloor\cdot\epsilon_{t}\ +\ \textup{Bin}\big{(}\epsilon_{t},\
M_{t}-\lfloor M_{t}\rfloor\big{)}$. As pointed out in Section 2, the two terms
in the above addition are independent.

Let us derive properties of the CMEM (3.15). First, we note that the first-
order stationarity given in Proposition 3.1.3 still holds here. Note that
$E[X_{t}\ |\ \mathcal{F}_{t-1}]=M_{t}$, and that the variance function
$\nu(\alpha)=(\alpha-\lfloor\alpha\rfloor)(1-\alpha+\lfloor\alpha\rfloor)$
according to Proposition 2.1. Hence, we have

$\displaystyle E[X_{t}^{2}\ |\ \mathcal{F}_{t-1}]\ =$ $\displaystyle\
E\Big{[}E\big{[}\big{(}\lfloor M_{t}\rfloor\cdot\epsilon_{t}\ +\
\textup{Bin}(\epsilon_{t},\ M_{t}-\lfloor M_{t}\rfloor)\big{)}^{2}\ |\
\epsilon_{t},\mathcal{F}_{t-1}\big{]}\ \big{|}\ \mathcal{F}_{t-1}\Big{]}$
$\displaystyle=$ $\displaystyle\ E\Big{[}(M_{t}-\lfloor M_{t}\rfloor)(1-\
M_{t}+\lfloor M_{t}\rfloor)\epsilon_{t}+M_{t}^{2}\,\epsilon_{t}^{2}\ \big{|}\
\mathcal{F}_{t-1}\Big{]}$ $\displaystyle\ =$ $\displaystyle\ (M_{t}-\lfloor
M_{t}\rfloor)(1-\ M_{t}+\lfloor M_{t}\rfloor)+(\sigma^{2}+1)M_{t}^{2},$
$\displaystyle V[X_{t}\ |\ \mathcal{F}_{t-1}]=$ $\displaystyle\ (M_{t}-\lfloor
M_{t}\rfloor)(1-\ M_{t}+\lfloor M_{t}\rfloor)+\sigma^{2}M_{t}^{2},$ (3.16)

also see (3.4). It is easy to see that $\nu(\alpha)\in[0,0.25]$, so
$\sigma^{2}M_{t}^{2}\leq V[X_{t}\ |\ \mathcal{F}_{t-1}]\leq
0.25+\sigma^{2}M_{t}^{2}$ holds. Similar to arguments in Section 3.1, if
$\sum_{i=1}^{p}a_{i}+\sum_{j=1}^{q}b_{j}<1$, we obtain the unconditional mean
of $X_{t}$ as

$\mu:=E[X_{t}]=\frac{a_{0}}{1-\sum_{i=1}^{p}a_{i}-\sum_{j=1}^{q}b_{j}}.$

Then, the variance is given by

$\displaystyle V[X_{t}]$ $\displaystyle=E\Big{[}(M_{t}-\lfloor
M_{t}\rfloor)(1-\ M_{t}+\lfloor
M_{t}\rfloor)+\sigma^{2}M_{t}^{2}\Big{]}+V[M_{t}]$
$\displaystyle=E\big{[}(M_{t}-\lfloor M_{t}\rfloor)(1-\ M_{t}+\lfloor
M_{t}\rfloor)\big{]}+\mu^{2}\sigma^{2}+(\sigma^{2}+1)\,V[M_{t}],$

where
$V[M_{t}]=\sum_{i=1}^{p}a_{i}\gamma_{X}(i)+\sum_{j=1}^{q}b_{j}\gamma_{M}(j)$.
If $p=q=1$, for example, then like in Section 3.1, we can insert
$(1-(a_{1}+b_{1})^{2}+a_{1}^{2})\,V[M_{t}]=a_{1}^{2}\,V[X_{t}]$. However, as
there is a floor operation in $V[X_{t}]$, we cannot obtain a closed-form
expression for $V[X_{t}]$ this time. But we have
$\mu^{2}\sigma^{2}+(\sigma^{2}+1)V[M_{t}]\leq V[X_{t}]\leq
0.25+\mu^{2}\sigma^{2}+(\sigma^{2}+1)V[M_{t}]$, so we can approximate
$V[X_{t}]$ by using these bounds.

The following theorem establishes the stationarity and ergodicity of
$(X_{t})$, which can be proven using the technique similar to that in Agosto
et al. (2016). Note that because of $E[X_{t}\ |\ \mathcal{F}_{t-1}]=M_{t}$,
also $E[X_{t}]=E[M_{t}]$ holds.

###### 3.2.1 Theorem

If $\sum_{i=1}^{p}a_{i}+\sum_{j=1}^{q}b_{j}<1$ holds, there exists a weakly
dependent stationary and ergodic solution to (3.15).

###### Proof.

We first verify that the process $(X_{t})$ satisfies the conditions of Theorem
3.1 in Doukhan and Wintenberger (2008) from which the theorem will follow. To
simplify the notation, assume without loss of generality $p=q$ in the
following. With $\beta(z)=1-\sum_{i=1}^{q}b_{i}z^{i},z\in\mathbb{C}$, the
assumption $\sum_{i=1}^{p}a_{i}+\sum_{j=1}^{q}b_{j}<1$ implies that
$\beta(z)^{-1}=\psi(z):=\sum_{i=0}^{\infty}\psi_{i}z^{i}$ is well-defined for
$|z|\leq 1+\eta$ with some $\eta>0$, where $\psi$ is exponentially decreasing
and defined recursively by $\psi_{0}=1$ and
$\psi_{n}=\sum_{i=1}^{n}b_{i}\psi_{n-i}\geq 0$ for $n\geq 1$. Next, with
$\alpha(z)=\sum_{i=1}^{q}a_{i}z^{i-1}$, note that $M_{t}$ defined in (1.4) can
be rewritten in terms of the backshift operator $B$ as
$\beta(B)M_{t}=a_{0}+\alpha(B)X_{t-1}$. Then,

$M_{t}=a_{0}/\beta(1)+\psi(B)\alpha(B)X_{t-1}=a_{0}/\beta(1)+\theta(B)X_{t-1},$

where $\theta(z):=\psi(z)\alpha(z)=\sum_{i=0}^{\infty}\theta_{i}z^{i}$ with
$\theta_{i}\geq 0$. Thus, $X_{t}$ satisfies
$X_{t}=F(X_{t-1},X_{t-2},\ldots;\eta_{t})$, where

$\displaystyle F(X_{t-1},X_{t-2},\ldots;\eta_{t})$ $\displaystyle=\lfloor
a_{0}/\beta(1)+\theta(B)X_{t-1}\rfloor\cdot\epsilon_{t}\ $
$\displaystyle~{}~{}~{}+\ \textup{Bin}_{t}\big{(}\epsilon_{t},\
a_{0}/\beta(1)+\theta(B)X_{t-1}-\lfloor
a_{0}/\beta(1)+\theta(B)X_{t-1}\rfloor\big{)}.$

Here, $\eta_{t}=(\epsilon_{t},\textup{Bin}_{t})^{\prime}$ is an i. i. d.
sequence with $\textup{Bin}_{t}$ denoting a binomial random variable being
executed at time $t$.

From Klenke and Mattner (2010), we know that the stochastic ordering
$\textup{Bin}(n,p_{1})\leq_{\rm st}\textup{Bin}(n,p_{2})$ holds if and only if
$p_{1}\leq p_{2}$. Denote $F_{m}$ as the cumulative distribution function
(cdf) of $\lfloor m\rfloor\cdot\epsilon\ +\ \textup{Bin}\big{(}\epsilon,\
m-\lfloor m\rfloor\big{)}$, with $m$ being its expectation. Let
$F_{m}^{-1}(u):=\inf\\{t\geq 0,F_{m}(t)\geq u\\}$ for $u\in[0,1]$. Let $U$ be
a uniform random variable on $(0,1)$, and define
$Z_{i}=F_{m_{i}}^{-1}(U),i=1,2$. Then, using arguments similar to those in
Lemma 4 of Chen et al. (2022), we have $E|Z_{1}-Z_{2}|=|m_{1}-m_{2}|$.

For any two deterministic sequences $(x_{t-i})_{i\geq 1}$ and
$(\tilde{x}_{t-i})_{i\geq 1}$, we can define $m_{t}$ and $\tilde{m}_{t}$
accordingly. Then, we obtain

$\displaystyle
E\big{|}F(x_{t-1},x_{t-2},\ldots;\eta_{t})-F(\tilde{x}_{t-1},\tilde{x}_{t-2},\ldots;\eta_{t})\big{|}$
$\displaystyle=E\Big{|}\lfloor m_{t}\rfloor\cdot\epsilon_{t}\ +\
\textup{Bin}(\epsilon_{t},\ m_{t}-\lfloor
m_{t}\rfloor)-\big{(}\lfloor\tilde{m}_{t}\rfloor\cdot\epsilon_{t}\ +\
\textup{Bin}(\epsilon_{t},\
\tilde{m}_{t}-\lfloor\tilde{m}_{t}\rfloor)\big{)}\Big{|}$
$\displaystyle=|m_{t}-\tilde{m}_{t}|\leq\sum_{i=0}^{\infty}\theta_{i}|x_{t-1-i}-\tilde{x}_{t-1-i}|=\sum_{i=1}^{\infty}\theta_{i-1}|x_{t-i}-\tilde{x}_{t-i}|,$

where we have used that $\theta_{i}\geq 0$. Observe that by the identity
$\alpha(z)=\beta(z)\theta(z)$, it follows that $\theta(1)<1$ if and only if
$\sum_{i=1}^{p}a_{i}+\sum_{j=1}^{q}b_{j}<1$. This verifies the conditions of
Theorem 3.1 in Doukhan and Wintenberger (2008), which, in turn, implies that
$(X_{t})$ is weakly dependent. ∎

###  4 Parameter Estimation and Model Diagnostics

In this section, we consider two kinds of estimation methods for the INGARCH-
CMEMs discussed in Section 3: QMLE and 2SWLSE, where the former has been
discussed by Christou and Fokianos (2014), Ahmad and Francq (2016) and
Aknouche et al. (2018), and the latter by Aknouche and Scotto (2023). To
implement these estimation approaches, the INGARCH structure as well as the
multiplicative operator have to be specified, but not the distribution of
$\epsilon_{t}$. Thus, QMLE and 2SWLSE approaches can be understood as semi-
parametric estimation approaches.

Throughout this section, let $X_{1},\ldots,X_{n}$ be observations from models
(3.5) or (3.15), respectively, with true parameters
$\mbox{\boldmath$\theta$}_{0}:=(a_{00},a_{01},\ldots,a_{0p},b_{01},\ldots,b_{0q})^{\prime}$
and $\sigma_{0}^{2}$. We assume that
$\mbox{\boldmath$\theta$}_{0}\in\mbox{\boldmath$\Theta$}\subset(0,\infty)\times[0,1)^{p+q}$
and $\sigma_{0}^{2}\in\Delta\subset(0,\infty)$, where the latter condition
ensures the non-degeneracy of $\epsilon_{t}$.

####  4.1 Quasi-Maximum Likelihood Estimation

Because models (3.5) and (3.15) have the same mean, the QMLE approach will
yield the same estimate of the regression parameters, although observations
come from different models.

#####  4.1.1 QMLE of Regression Parameters

Christou and Fokianos (2014) considered Poisson QMLE for negative-binomial
INGARCH models, and Ahmad and Francq (2016) studied the asymptotic
distribution of the estimator when the parameter belongs to the interior of
the parameter space and when it lies at the boundary.

For any generic
$\mbox{\boldmath$\theta$}:=(a_{0},a_{1},\ldots,a_{p},b_{1},\ldots,b_{q})^{\prime}\in\mbox{\boldmath$\Theta$}$,
define

$M_{t}(\mbox{\boldmath$\theta$})\ =\
a_{0}+\sum_{i=1}^{p}a_{i}\,X_{t-i}+\sum_{j=1}^{q}b_{j}\,M_{t-j}(\mbox{\boldmath$\theta$}),~{}~{}~{}t\in\mathbb{Z},$
(4.1)

where $M_{t}(\mbox{\boldmath$\theta$}_{0})=M_{t}$. Let
$l_{t}(\mbox{\boldmath$\theta$})=X_{t}\log
M_{t}(\mbox{\boldmath$\theta$})-M_{t}(\mbox{\boldmath$\theta$})$, then the
Poisson QMLE (PQ) of $\theta_{0}$ is defined as

$\hat{}\mbox{\boldmath$\theta$}_{\textup{PQ}}=\arg\max_{\mbox{\scriptsize\boldmath$\theta$}\in\mbox{\scriptsize\boldmath$\Theta$}}\sum_{t=1}^{n}l_{t}(\mbox{\boldmath$\theta$}).$
(4.2)

Under some regularity conditions, using techniques similar to those in Theorem
2 of Christou and Fokianos (2014),
$\hat{}\mbox{\boldmath$\theta$}_{\textup{PQ}}$ can be proven to be consistent
and asymptotically normal, i. e.,

$\begin{array}[]{@{}l}\sqrt{n}(\hat{}\mbox{\boldmath$\theta$}_{\textup{PQ}}-\mbox{\boldmath$\theta$}_{0})\xrightarrow{d}\textup{N}\big{(}\mbox{\boldmath$0$},\
\mbox{{{G}}}^{-1}(\mbox{\boldmath$\theta$}_{0})\mbox{{{G}}}_{1}(\mbox{\boldmath$\theta$}_{0},\sigma_{0}^{2})\mbox{{{G}}}^{-1}(\mbox{\boldmath$\theta$}_{0})\big{)},\quad\text{where}\\\\[8.61108pt]
\displaystyle\mbox{{{G}}}(\mbox{\boldmath$\theta$})=-E\left[\frac{\partial^{2}l_{t}(\mbox{\boldmath$\theta$})}{\partial\mbox{\boldmath$\theta$}\partial\mbox{\boldmath$\theta$}^{\prime}}\right]=E\left[\frac{1}{M_{t}(\mbox{\boldmath$\theta$})}\,\frac{\partial
M_{t}(\mbox{\boldmath$\theta$})}{\partial\mbox{\boldmath$\theta$}}\,\frac{\partial
M_{t}(\mbox{\boldmath$\theta$})}{\partial\mbox{\boldmath$\theta$}^{\prime}}\right],\\\
\\\\[-8.61108pt]
\displaystyle\mbox{{{G}}}_{1}(\mbox{\boldmath$\theta$},\sigma^{2})=E\left[V\left[\frac{\partial
l_{t}(\mbox{\boldmath$\theta$})}{\partial\mbox{\boldmath$\theta$}}\bigg{|}\mathcal{F}_{t-1}\right]\right]=E\left[\frac{v_{t}(\mbox{\boldmath$\theta$},\sigma^{2})}{M_{t}^{2}(\mbox{\boldmath$\theta$})}\,\frac{\partial
M_{t}(\mbox{\boldmath$\theta$})}{\partial\mbox{\boldmath$\theta$}}\,\frac{\partial
M_{t}(\mbox{\boldmath$\theta$})}{\partial\mbox{\boldmath$\theta$}^{\prime}}\right],\end{array}$

where $v_{t}(\mbox{\boldmath$\theta$},\sigma^{2}):=V[X_{t}\ |\
\mathcal{F}_{t-1}]$. In practice, maximization of (4.2) is implemented by
assuming some (fixed or random) positive starting values
$\tilde{M}_{0}(\mbox{\boldmath$\theta$}),\ldots,\tilde{M}_{1-q}(\mbox{\boldmath$\theta$}),X_{0},\ldots,X_{1-p}$.
Let $\tilde{M}_{t}(\mbox{\boldmath$\theta$})$ be an observable counterpart of
$M_{t}(\mbox{\boldmath$\theta$})$, given by the recursion

$\tilde{M}_{t}(\mbox{\boldmath$\theta$})=a_{0}+\sum_{i=1}^{p}a_{i}X_{t-i}+\sum_{j=1}^{q}b_{j}\tilde{M}_{t-j}(\mbox{\boldmath$\theta$}),~{}~{}~{}t\geq
1.$ (4.3)

Then, we have an observable counterpart of $l_{t}(\mbox{\boldmath$\theta$})$
as
$\tilde{l}_{t}(\mbox{\boldmath$\theta$})=X_{t}\log\tilde{M}_{t}(\mbox{\boldmath$\theta$})-\tilde{M}_{t}(\mbox{\boldmath$\theta$})$,
and we can approximate $l_{t}(\mbox{\boldmath$\theta$})$ by
$\tilde{l}_{t}(\mbox{\boldmath$\theta$})$.

Besides the Poisson QMLE, there are other versions of QMLE as well. Aknouche
et al. (2018) and Aknouche and Francq (2021) discussed the negative-binomial
QMLE and the exponential QMLE, respectively. As noted by Aknouche and Scotto
(2023), the asymptotic efficiency of the QMLE depends on the model’s variance-
to-mean relationship $V[X_{t}\ |\ \mathcal{F}_{t-1}]=g(M_{t})$, where $g$ is a
positive real-valued function. Thus, each one of the Poisson, negative-
binomial, and exponential QMLE is asymptotically optimal whenever
$g(x)=c_{1}x$, $g(x)=c_{2}x(1+c_{3}x)$, and $g(x)=c_{4}x^{2}$, respectively,
for any $c_{1},c_{2},c_{3},c_{4}>0$. From (3.6) and (3.2), we may conclude
that negative-binomial and exponential QMLEs are asymptotically efficient for
models (3.5) and (3.15), respectively. We should point out that they are not
universally (asymptotically) efficient but only within the class of all QMLEs
in the linear exponential family (see also Wooldridge (1999) and the
references therein).

For model (3.5), the negative-binomial QMLE (NQ) is defined as

$\hat{}\mbox{\boldmath$\theta$}_{\textup{NQ}}=\arg\max_{\mbox{\scriptsize\boldmath$\theta$}\in\mbox{\scriptsize\boldmath$\Theta$}}\sum_{t=1}^{n}\Big{\\{}X_{t}\log\tilde{M}_{t}(\mbox{\boldmath$\theta$})-(r+X_{t})\log(r+\tilde{M}_{t}(\mbox{\boldmath$\theta$}))\Big{\\}}$

for some fixed $r>0$ (in their simulations, Aknouche and Scotto (2023) used
$r=1$ throughout). Note that Aknouche et al. (2018) also proposed another
estimation approach, which consists of four stages estimating both conditional
mean and dispersion parameters.

For model (3.15), the exponential QMLE (EQ) is defined as

$\hat{}\mbox{\boldmath$\theta$}_{\textup{EQ}}=\arg\max_{\mbox{\scriptsize\boldmath$\theta$}\in\mbox{\scriptsize\boldmath$\Theta$}}\sum_{t=1}^{n}\Big{\\{}-\log\tilde{M}_{t}(\mbox{\boldmath$\theta$})-X_{t}/\tilde{M}_{t}(\mbox{\boldmath$\theta$})\Big{\\}}.$

Large-sample properties of $\hat{}\mbox{\boldmath$\theta$}_{\textup{NQ}}$ and
$\hat{}\mbox{\boldmath$\theta$}_{\textup{EQ}}$ are discussed by Aknouche et
al. (2018) and Aknouche and Francq (2021), respectively, which will not be
repeated here. From a practical point of view, the crucial difference are
modified expressions for the matrices $\mbox{{{G}}}(\mbox{\boldmath$\theta$})$
and $\mbox{{{G}}}_{1}(\mbox{\boldmath$\theta$},\sigma^{2})$ involved in the
asymptotic normal distribution, namely

$\textstyle\mbox{{{G}}}(\mbox{\boldmath$\theta$})=E\left[\frac{1}{M_{t}(\mbox{\scriptsize\boldmath$\theta$})(r+M_{t}(\mbox{\scriptsize\boldmath$\theta$}))}\,\frac{\partial
M_{t}(\mbox{\scriptsize\boldmath$\theta$})}{\partial\mbox{\scriptsize\boldmath$\theta$}}\,\frac{\partial
M_{t}(\mbox{\scriptsize\boldmath$\theta$})}{\partial\mbox{\scriptsize\boldmath$\theta$}^{\prime}}\right],\
\mbox{{{G}}}_{1}(\mbox{\boldmath$\theta$},\sigma^{2})=E\left[\frac{v_{t}(\mbox{\scriptsize\boldmath$\theta$},\sigma^{2})}{M_{t}^{2}(\mbox{\scriptsize\boldmath$\theta$})(r+M_{t}(\mbox{\scriptsize\boldmath$\theta$}))^{2}}\,\frac{\partial
M_{t}(\mbox{\scriptsize\boldmath$\theta$})}{\partial\mbox{\scriptsize\boldmath$\theta$}}\,\frac{\partial
M_{t}(\mbox{\scriptsize\boldmath$\theta$})}{\partial\mbox{\scriptsize\boldmath$\theta$}^{\prime}}\right]\
\text{for }\hat{}\mbox{\boldmath$\theta$}_{\textup{NQ}},$

and

$\textstyle\mbox{{{G}}}(\mbox{\boldmath$\theta$})=E\left[\frac{1}{M_{t}^{2}(\mbox{\scriptsize\boldmath$\theta$})}\,\frac{\partial
M_{t}(\mbox{\scriptsize\boldmath$\theta$})}{\partial\mbox{\scriptsize\boldmath$\theta$}}\,\frac{\partial
M_{t}(\mbox{\scriptsize\boldmath$\theta$})}{\partial\mbox{\scriptsize\boldmath$\theta$}^{\prime}}\right],\
\mbox{{{G}}}_{1}(\mbox{\boldmath$\theta$},\sigma^{2})=E\left[\frac{v_{t}(\mbox{\scriptsize\boldmath$\theta$},\sigma^{2})}{M_{t}^{4}(\mbox{\scriptsize\boldmath$\theta$})}\,\frac{\partial
M_{t}(\mbox{\scriptsize\boldmath$\theta$})}{\partial\mbox{\scriptsize\boldmath$\theta$}}\,\frac{\partial
M_{t}(\mbox{\scriptsize\boldmath$\theta$})}{\partial\mbox{\scriptsize\boldmath$\theta$}^{\prime}}\right]\
\text{for }\hat{}\mbox{\boldmath$\theta$}_{\textup{EQ}}.$

For any of the three QMLE approaches, the estimates are determined by
numerical optimization. To initialize the optimization routine, a moment
estimate of $\mbox{\boldmath$\theta$}_{0}$ might be used, which is easily
obtained from sample mean and ACF by solving the moment properties provided by
Section 3.

#####  4.1.2 Estimation of Parameter $\sigma^{2}$

It remains to estimate the parameter value $\sigma_{0}^{2}$, which corresponds
to the variance of $\epsilon_{t}$. For model (3.5) with Poisson-counting
series, the conditional variance is given by $V[X_{t}\ |\
\mathcal{F}_{t-1}]=M_{t}+\sigma^{2}M_{t}^{2}$, recall (3.6). Thus, let
$e_{t}=(X_{t}-M_{t})^{2}-(M_{t}+\sigma^{2}\,M_{t}^{2})$. Then, we have
$(X_{t}-M_{t})^{2}=M_{t}+\sigma^{2}\,M_{t}^{2}+e_{t}$, or equivalently,

$\frac{(X_{t}-M_{t}(\mbox{\boldmath$\theta$}))^{2}-M_{t}(\mbox{\boldmath$\theta$})}{M_{t}^{2}(\mbox{\boldmath$\theta$})}=\sigma^{2}+\tau_{t},$
(4.4)

where $\tau_{t}=e_{t}/M_{t}^{2}(\mbox{\boldmath$\theta$})$ is a term of
martingale difference. Replacing $M_{t}(\mbox{\boldmath$\theta$})$ by
$\hat{M}_{t}:=\tilde{M}_{t}(\hat{}\mbox{\boldmath$\theta$}_{\textup{PQ}})$ or
$\tilde{M}_{t}(\hat{}\mbox{\boldmath$\theta$}_{\textup{NQ}})$ or
$\tilde{M}_{t}(\hat{}\mbox{\boldmath$\theta$}_{\textup{EQ}})$, respectively
(one might also use the moment estimate of $\mbox{\boldmath$\theta$}_{0}$), we
obtain the least-squares estimate of $\sigma_{0}^{2}$ as

$\hat{\sigma}_{\textup{Poi}}^{2}=\frac{1}{n}\sum_{t=1}^{n}\frac{(X_{t}-\hat{M}_{t})^{2}-\hat{M}_{t}}{\hat{M}_{t}^{2}}.$

Under some regularity conditions, using techniques similar to those in Theorem
3.2 of Aknouche and Scotto (2023), $\hat{\sigma}_{\textup{Poi}}^{2}$ can be
proven to be consistent and asymptotically normal, i. e.,

$\sqrt{n}(\hat{\sigma}_{\textup{Poi}}^{2}-\sigma_{0}^{2})\xrightarrow{d}\textup{N}(0,\Lambda),\quad\text{where}\quad\Lambda=E\left[\left(\frac{(X_{t}-M_{t}(\mbox{\boldmath$\theta$}_{0}))^{2}-(M_{t}(\mbox{\boldmath$\theta$}_{0})+\sigma^{2}_{0}\,M_{t}^{2}(\mbox{\boldmath$\theta$}_{0}))}{M_{t}^{2}(\mbox{\boldmath$\theta$}_{0})}\right)^{2}\right].$

In the case of an NB-counting series, we have $V[X_{t}\ |\
\mathcal{F}_{t-1}]=M_{t}(1+M_{t})+\sigma^{2}M_{t}^{2}$ according to (3.6), so
we would define the least-squares estimate of $\sigma_{0}^{2}$ as

$\hat{\sigma}_{\textup{NB}}^{2}=\frac{1}{n}\sum_{t=1}^{n}\frac{(X_{t}-\hat{M}_{t})^{2}-\hat{M}_{t}(1+\hat{M}_{t})}{\hat{M}_{t}^{2}}\
=\ \hat{\sigma}_{\textup{Poi}}^{2}-1.$

But as can be seen by the last equality, we just end up with the estimator of
the Poi-case reduced by one. This seems plausible as now a larger part of the
variation is explained by the counting series.

Finally, for model (3.15), the conditional variance is given by $V[X_{t}\ |\
\mathcal{F}_{t-1}]=(M_{t}-\lfloor M_{t}\rfloor)(1-\ M_{t}+\lfloor
M_{t}\rfloor)+\sigma^{2}M_{t}^{2}$, recall (3.2). Thus, let
$e_{t}=(X_{t}-M_{t})^{2}-\big{(}(M_{t}-\lfloor M_{t}\rfloor)(1-\ M_{t}+\lfloor
M_{t}\rfloor)+\sigma^{2}\,M_{t}^{2}\big{)}$. Then, we have

$\hat{\sigma}_{\textup{Bin}}^{2}=\frac{1}{n}\sum_{t=1}^{n}\frac{(X_{t}-\hat{M}_{t})^{2}-(\hat{M}_{t}-\lfloor\hat{M}_{t}\rfloor)(1-\hat{M}_{t}+\lfloor\hat{M}_{t}\rfloor)}{\hat{M}_{t}^{2}},$

and
$\sqrt{n}(\hat{\sigma}_{\textup{Bin}}^{2}-\sigma_{0}^{2})\xrightarrow{d}\textup{N}(0,\Lambda)$
with

$\Lambda=E\left[\left(\frac{(X_{t}-M_{t}(\mbox{\boldmath$\theta$}_{0}))^{2}-\big{(}(M_{t}(\mbox{\boldmath$\theta$}_{0})-\lfloor
M_{t}(\mbox{\boldmath$\theta$}_{0})\rfloor)(1-\
M_{t}(\mbox{\boldmath$\theta$}_{0})+\lfloor
M_{t}(\mbox{\boldmath$\theta$}_{0})\rfloor)+\sigma^{2}_{0}\,M_{t}^{2}(\mbox{\boldmath$\theta$}_{0})\big{)}}{M_{t}^{2}(\mbox{\boldmath$\theta$}_{0})}\right)^{2}\right].$

####  4.2 Two-Stage Weighted Least-Squares Estimation

As an alternative to the above QMLEs, Aknouche and Scotto (2023) proposed a
2SWLSE approach. The 2SWLSE is not asymptotically efficient in the sense that
its asymptotic variance reaches the Cramer–Rao bound, but is only “never
asymptotically less efficient” than any QMLE in the linear exponential family.
When the conditional distribution of the model is known and is a member of the
exponential family, then the 2SWLS estimate is asymptotically efficient (in
the sense of the Cramer–Rao bound). See Aknouche (2012) and Aknouche and
Scotto (2023) for more details and other models. We shall describe the
estimation procedure as follows. For convenience, we again denote
$v_{t}(\mbox{\boldmath$\theta$},\sigma^{2}):=V[X_{t}\ |\ \mathcal{F}_{t-1}]$.
For an arbitrary fixed weighting point, say
$(\mbox{\boldmath$\theta$}^{\prime}_{\ast},\sigma^{2}_{\ast})^{\prime}\in\mbox{\boldmath$\Theta$}\times\Delta$
(in practice, the aforementioned moment estimates of
$(\mbox{\boldmath$\theta$}_{0},\sigma^{2}_{0})$ might be used), a first-stage
weighted least-squares estimator (1WLSE, briefly “1W”) is defined as

$\hat{}\mbox{\boldmath$\theta$}_{\textup{1W}}=\arg\min_{\mbox{\scriptsize\boldmath$\theta$}\in\mbox{\scriptsize\boldmath$\Theta$}}\sum_{t=1}^{n}\frac{(X_{t}-\tilde{M}_{t}(\mbox{\boldmath$\theta$}))^{2}}{v_{t}(\mbox{\boldmath$\theta$}_{\ast},\sigma^{2}_{\ast})}.$
(4.5)

We next estimate $\sigma^{2}$ like in Section 4.1.2, and we denote the
estimator as $\hat{\sigma}_{n}^{2}$. But here, we replace
$M_{t}(\mbox{\boldmath$\theta$})$ by
$\hat{M}_{t}=\tilde{M}_{t}(\hat{}\mbox{\boldmath$\theta$}_{\textup{1W}})$.
After having consistent estimates of $\mbox{\boldmath$\theta$}_{0}$ and
$\sigma_{0}^{2}$, the optimal weight is estimated by
$v_{t}(\hat{}\mbox{\boldmath$\theta$}_{\textup{1W}},\hat{\sigma}^{2}_{n})$,
and a second-stage weighted least-squares estimator (2WLSE, briefly “2W”) of
$\mbox{\boldmath$\theta$}_{0}$ has the form

$\hat{}\mbox{\boldmath$\theta$}_{\textup{2W}}=\arg\min_{\mbox{\scriptsize\boldmath$\theta$}\in\mbox{\scriptsize\boldmath$\Theta$}}\sum_{t=1}^{n}\frac{(X_{t}-\tilde{M}_{t}(\mbox{\boldmath$\theta$}))^{2}}{v_{t}(\hat{}\mbox{\boldmath$\theta$}_{\textup{1W}},\hat{\sigma}^{2}_{n})}.$
(4.6)

Under some regularity conditions, using techniques similar to those in Theorem
3.3 in Aknouche and Scotto (2023),
$\hat{}\mbox{\boldmath$\theta$}_{\textup{2W}}$ can be proven to be strongly
consistent and asymptotically normal, i. e.,
$\sqrt{n}(\hat{}\mbox{\boldmath$\theta$}_{\textup{2W}}-\mbox{\boldmath$\theta$}_{0})\xrightarrow{d}\textup{N}(\mbox{\boldmath$0$},\mbox{{{J}}}^{-1}(\mbox{\boldmath$\theta$}_{0},\sigma_{0}^{2})),$
where

$\mbox{{{J}}}(\mbox{\boldmath$\theta$}_{0},\sigma_{0}^{2})=E\left[\frac{1}{v_{t}(\mbox{\boldmath$\theta$}_{0},\sigma^{2}_{0})}\,\frac{\partial
M_{t}(\mbox{\boldmath$\theta$}_{0})}{\partial\mbox{\boldmath$\theta$}}\,\frac{\partial
M_{t}(\mbox{\boldmath$\theta$}_{0})}{\partial\mbox{\boldmath$\theta$}^{\prime}}\right].$

Aknouche and Scotto (2023) stated that $\hat{\theta}_{\textup{2W}}$ is
asymptotically more efficient than
$\hat{}\mbox{\boldmath$\theta$}_{\textup{1W}}$ and the QMLEs. Note that if
$b_{0j}=0$ for $1\leq j\leq q$, then
$\hat{}\mbox{\boldmath$\theta$}_{\textup{1W}}$ and
$\hat{}\mbox{\boldmath$\theta$}_{\textup{2W}}$ have a closed form, which is
omitted here.

####  4.3 Approaches for Model Diagnostics

To check the adequacy of the fitted INGARCH-CMEMs (or to compare the
performance of different candidate models), several approaches are possible,
see Section 2.4 in Weiß (2018) for a general survey. It should be noted,
however, that we solely considered semi-parametric estimation approaches, i.
e., we did not specify the full distribution of the considered CMEM. Thus,
mainly moment-related diagnostic checks are relevant here.

The first option is to compare the fitted model’s mean, variance, and ACF
(which can be computed from the results derived in Section 3) to the
corresponding sample moments.

Second, to check the adequacy of conditional moments, it is common to consider
the standardized Pearson residuals $(R_{t})$, which are given by

$R_{t}(\hat{}\mbox{\boldmath$\theta$},\hat{\sigma}^{2})\ =\
\frac{X_{t}-E[X_{t}\ |\ \mathcal{F}_{t-1};\
\hat{}\mbox{\boldmath$\theta$},\hat{\sigma}^{2}]}{\sqrt{V[X_{t}\ |\
\mathcal{F}_{t-1};\
\hat{}\mbox{\boldmath$\theta$},\hat{\sigma}^{2}]}}\overset{\eqref{CountMEMcvar}}{=}\frac{X_{t}-M_{t}(\hat{}\mbox{\boldmath$\theta$})}{\sqrt{\nu\big{(}M_{t}(\hat{}\mbox{\boldmath$\theta$})\big{)}+M_{t}^{2}(\hat{}\mbox{\boldmath$\theta$})\cdot\hat{\sigma}^{2}}}$
(4.7)

for the models considered here. For an adequate model, the Pearson residuals
have mean zero, variance one, and are uncorrelated. In particular, Aknouche
and Scotto (2023) used the mean squared Pearson residual,
$\textup{MSPR}=\frac{1}{n}\sum_{t=1}^{n}R_{t}^{2}(\hat{}\mbox{\boldmath$\theta$},\hat{\sigma}^{2})$,
for checking the adequacy of the dispersion structure, also see Aleksandrov
and Weiß (2020) for a more general analysis of this statistic. For an
andequate model, the MSPR should be close to one.

Third, for INGARCH-type models, also another type of residuals is commonly
used, where the $X_{t}$ (or powers thereof) are scaled by $M_{t}$, see Weiß et
al. (2017) for a general discussion. For INGARCH-type and multiplicative
models such as the one proposed by Zhu and Wang (2010) and Aknouche and Scotto
(2023) and the CMEMs considered here, especially the scaled residuals
$S_{t}(\hat{}\mbox{\boldmath$\theta$})=X_{t}/M_{t}(\hat{}\mbox{\boldmath$\theta$})$
are reasonable, and Aknouche and Scotto (2023) proposed to compute the sample
mean and variance thereof (MSR and VSR, respectively). Note that if $(X_{t})$
stems from the CMEM (3.5) satisfying (3.6), then

$\begin{array}[]{@{}rl}E[S_{t}(\mbox{\boldmath$\theta$})]\
=&E\big{[}E[X_{t}/M_{t}\ |\ \mathcal{F}_{t-1}\big{]}\ =\ 1,\\\\[4.30554pt]
V[S_{t}(\mbox{\boldmath$\theta$})]\ =&0\ +\ E\big{[}V[X_{t}/M_{t}\ |\
\mathcal{F}_{t-1}\big{]}\ =\
\sigma^{2}+E\big{[}\nu(M_{t})/M_{t}^{2}\big{]}.\end{array}$

The last term, $E\big{[}\nu(M_{t})/M_{t}^{2}\big{]}$, depends on the specific
type of CMEM. For the compounding CMEM with Poisson counting series, we have
$V[S_{t}(\mbox{\boldmath$\theta$})]=\sigma^{2}+E[1/M_{t}]$, while
$V[S_{t}(\mbox{\boldmath$\theta$})]=\sigma^{2}+1+E[1/M_{t}]$ for the NB-
counting series. The latter result also shows that an NB-counting series is
possible only if the scaled residuals’ variance is larger than one. Generally,
if $\nu$ is differentiable, we can approximate
$E\big{[}\nu(M_{t})/M_{t}^{2}\big{]}$ by a second-order Taylor expansion
around the marginal mean $\mu=E[X_{t}]$ as

$E\left[\frac{\nu(M_{t})}{M_{t}^{2}}\right]\ \approx\
\frac{\nu(\mu)}{\mu^{2}}+\frac{\mu^{2}\,\nu^{\prime\prime}(\mu)-4\,\mu\,\nu^{\prime}(\mu)+6\,\nu(\mu)}{2\,\mu^{4}}\,V[X_{t}],\quad
E\left[\frac{1}{M_{t}}\right]\ \approx\
\frac{1}{\mu}+\frac{1}{\mu^{3}}\,V[X_{t}].$

For the binomial CMEM (3.15), the function
$\nu(\alpha)=(\alpha-\lfloor\alpha\rfloor)(1-\alpha+\lfloor\alpha\rfloor)$ is
not differentiable, but it is bounded by $\nu(\alpha)\in[0,0.25]$. Thus,
$V[S_{t}(\mbox{\boldmath$\theta$})]\leq\sigma^{2}+E[0.25/M_{t}^{2}]$, i. e.,
$V[S_{t}(\mbox{\boldmath$\theta$})]$ often takes a value close to
$\sigma^{2}$.

As a fourth type of diagnostic check, Aknouche and Scotto (2023) proposed to
compute the mean absolute residual (MAR) defined as
$\textup{MAR}=\frac{1}{n}\sum_{t=1}^{n}|X_{t}-M_{t}(\hat{}\mbox{\boldmath$\theta$})|$,
where smaller values indicate a better model fit.

###  5 Simulation Study

To analyze the finite-sample performance of the estimation approaches
discussed in Section 4, a simulation study was done. For different types of
data-generating process (DGP) and sample size $n\in\\{300,600,1000\\}$, we
simulated 10,000 time series and computed estimates according to the methods
PQ, NQ, EQ, and 2W (by using moment estimates for initialization, recall
Example 3.1.5). As the DGP, we considered twelve types of INGARCH$(1,1)$-CMEM,
namely:

  * •

$(a_{0},a_{1},b_{1})=(2.8,0.4,0.2)$ and $\sigma^{2}=1$ or $\sigma^{2}=0.4$,
similar to the simulation scenario in Section 4 of Aknouche and Scotto (2023);

  * •

$(a_{0},a_{1},b_{1})=(3,0.35,0.5)$ and $\sigma^{2}=1$ or $\sigma^{2}=0.1$,
similar to the first data example in Section 6;

  * •

$(a_{0},a_{1},b_{1})=(1,0.25,0.65)$ and $\sigma^{2}=1$ or $\sigma^{2}=0.4$,
similar to the second data example in Section 6;

and these parametrizations are combined with either the compounding CMEM (3.5)
with Poisson-counting series, or the binomial CMEM (3.15). The model order
$(1,1)$ is relevant for many INGARCH applications in practice, including our
subsequent data examples.

The considered values of $\sigma^{2}$ are realized by an appropriate choice of
the innovations’ distribution. Recall that the standard INGARCH model
corresponds to a compounding INGARCH-CMEM with deterministic $\epsilon_{t}=1$
(thus $\sigma^{2}=0$). In applications, the variance
$\sigma^{2}=V[\epsilon_{t}]$ of a fitted INGARCH-CMEM is often rather small,
see Section 6 for such examples. Furthermore, an upper bound for $\sigma^{2}$
is also implied by Proposition 3.1.4 if requiring for second-order
stationarity. In the case of a Poisson-counting series, for example,
$\sigma^{2}<(1-(a_{1}+b_{1})^{2})/a_{1}^{2}$ has to hold. If using
$\sigma^{2}=1$ in the above DGP specifications, this value is achieved by
$\textup{Poi}(1)$-distributed $(\epsilon_{t})$. Our model for achieving
$\sigma^{2}\in(0,1)$, by contrast, is a three-point distribution with range
$\\{0,1,2\\}$. From the discussion at the end of Section 3.1, we know that we
can choose $p_{\epsilon,2}\in(0,0.5)$ arbitrarily. Then,
$p_{\epsilon,0}=p_{\epsilon,2}$ and $p_{\epsilon,1}=1-2\,p_{\epsilon,2}$
follow such that $\sigma^{2}=E[\epsilon_{t}^{2}]-1=2\,p_{\epsilon,2}\in(0,1)$.

Estimation of $a_{0}=2.8$ |  Estimation of $a_{1}=0.4$
---|---
|
Estimation of $b_{1}=0.2$ |  Estimation of $\sigma^{2}=1$
|
Figure 1: INGARCH$(1,1)$-CMEM with Poi-counting series: four groups of
boxplots for different estimation approaches, resulting from simulated count
time series of lengths $n\in\\{300,600,1000\\}$ per group.

Detailed simulation results are summarized in Tables S.1–S.6 in the
Supplement. Furthermore, illustrative boxplots for one example scenario
(corresponding to the left part of Table S.1) are presented in Figure 1.
Before discussing these results, a few notes on Tables S.1–S.6. The standard
errors (SEs) of the considered estimators are approximated based on the
asymptotic normal distributions discussed in Section 4. For this purpose, the
means $E[\cdot]$ involved in the expressions for
$\mbox{{{G}}}(\mbox{\boldmath$\theta$}_{0})$,
$\mbox{{{G}}}_{1}(\mbox{\boldmath$\theta$}_{0},\sigma_{0}^{2})$, etc. are
replaced by sample means, and the respective estimates are plugged-in instead
of the true parameter values. The required partial derivatives $\frac{\partial
M_{t}(\mbox{\scriptsize\boldmath$\theta$})}{\partial\mbox{\scriptsize\boldmath$\theta$}}$
can be computed recursively based on (1.4). For model order $p=q=1$, for
example, we get the recursive scheme

$\frac{\partial
M_{t}(\mbox{\boldmath$\theta$})}{\partial\mbox{\boldmath$\theta$}}=\big{(}1,X_{t-1},M_{t-1}(\mbox{\boldmath$\theta$})\big{)}\
+\ b_{1}\,\frac{\partial
M_{t-1}(\mbox{\boldmath$\theta$})}{\partial\mbox{\boldmath$\theta$}},$

also see Ferland et al. (2006, pp. 936–937). In Tables S.1–S.6, two types of
SE are presented: the “true” SE, computed as the (trimmed) sample standard
deviation among the simulated estimates and abbreviated as “SSE”, and the mean
among the aforementioned approximate SEs, abbreviated as “ASE”.

The results in Figure 1 show that bias and SE decrease with increasing sample
size $n$ (also see Tables S.1–S.6), confirming the consistency of the
estimators. But Figure 1 also demonstrates one of the general patterns in
Tables S.1–S.6, namely that the estimates of the AR-parameter $a_{1}$ and the
variance $\sigma^{2}$ are only mildly affected by bias, whereas the estimates
of the feedback parameter $b_{1}$ and especially of the intercept parameter
$a_{0}$ might be biased also for larger sample sizes. From Tables S.1–S.6, it
gets clear that this bias is mainly caused by the extent of $b_{1}$: the
larger $b_{1}$, the stronger the bias (especially of $a_{0}$). The problems in
estimating the intercept term of INGARCH-type models are well known in the
literature, see Fokianos et al. (2009) for example.

Further details can be recognized from Tables S.1–S.6. If $\sigma^{2}=1$ (so a
large $\sigma^{2}$, always see the left part of Tables S.1–S.6), the PQ-method
often causes larger bias and SE than NQ, EQ (especially for $a_{0}$), whereas
we often see the opposite ranking if $\sigma^{2}<1$ (right part of Tables
S.1–S.6). This goes along with the discussion of asymptotic optimality in
Aknouche and Scotto (2023), recall Section 4.1.1, because for small
$\sigma^{2}$, the conditional variance of (3.5) is close to a linear function
in $M_{t}$. The 2W-method often has bias and SE between the ones of PQ and NQ,
EQ. Comparing Tables S.1 vs. S.2, S.3 vs. S.4, and S.5 vs. S.6, i. e., the
case of a Poi-counting series to a binomial one, we do not see big differences
in performance, i. e., the estimation approaches generally perform similar for
both classes of DGP. An exception is the bias of $a_{0}$, which is somewhat
larger for the Bin-counting series if $b_{1}\geq 0.5$. Finally, the values of
SSE and ASE are often close to each other. If there are visible differences
between SSE and ASE, these mainly happen for $n=300$, and the strongest
absolute discrepancies are observed for the intercept term $a_{0}$. The
relative deviation, however, is sometimes stronger for $\sigma^{2}$. In
particular, the 2W-method sometimes leads to visibly larger SSE than ASE (e.
g., lower left block in Table S.6), which is caused by a few extreme values of
the estimates. The QMLE have less tendency to produce extreme estimates.
Overall, the approximate SE seems to give reasonable insight into the true SE
of the estimator in most cases.

We conclude this section with some simulation experiments on the effects of a
possible model misspecification. Here, we focus on the DGPs with
parametrization $(a_{0},a_{1},b_{1})=(2.8,0.4,0.2)$ and $\sigma^{2}=1$ or
$\sigma^{2}=0.4$ (see Tables S.1 and S.2). In our first experiment, we fitted
an INGARCH$(1,1)$-CMEM to the generated time series by assuming the wrong
multiplicative operator:

  * •

for the DGP with Poi-counting series from Table S.1, the misspecified model
assumes a Bin-counting series, leading to the simulation results of Table S.7;

  * •

for the DGP with Bin-counting series from Table S.2, the misspecified model
assumes a Poi-counting series, leading to the simulation results of Table S.8.

As the estimation approaches PQ, NQ, and EQ for $(a_{0},a_{1},b_{1})$ do not
depend on the specific structure of the multiplicative operator, the
corresponding simulated means and SEs in Tables S.7 and S.8 are identical to
those in Tables S.1 and S.2, respectively. The 2W-method, by contrast, depends
on the choice of the multiplicative operator, but as can be seen by comparing
Table S.7 to S.1 and Table S.8 to S.2, the effect of the misspecified operator
is negligible. So altogether, the estimation of $(a_{0},a_{1},b_{1})$ is
robust against such model misspecification. The approximate SEs, in turn, are
always affected by the misspecified operator, but the difference between
correct and misspecified ASE is again small for $(a_{0},a_{1},b_{1})$.

A clear difference is only observed for the estimation of $\sigma^{2}$ (and
its corresponding ASE), which is plausible in view of Section 4.1.2, where
tailor-made estimation methods are presented for each type of counting series.
However, in view of George Box’ famous words “all models are wrong but some
are useful” (i. e., real-world data will never follow exactly one of the
considered models, but some of these models might provide a good approximation
to the true DGP), the practitioner is probably more interested in the
performance of the fitted models. In Table S.9, the mean MAR values are
provided if estimation is done based on the correct vs. the misspecified
model. For the estimation approaches PQ, NQ, and EQ, the MAR values are
identical as the estimates for $(a_{0},a_{1},b_{1})$ do not depend on the
specific multiplicative operator. Also for the 2W-method, there is an at most
slight difference in the mean MARs, so the MAR performance is very robust
against such misspecification. Analogous results but concerning the MSPR
performance are shown in Table S.10. If the Poi-counting series is
misspecified as a Bin-counting series (upper part of Table S.10), there is
again an at most tiny difference in the resulting MSPRs. A somewhat stronger
effect (although still quite small) is observed in the opposite scenario,
where the Bin-counting series is misspecified as a Poi-counting series. So it
seems that the CMEM with Bin-counting series can better adapt to the given DGP
than the Poi-one.

As a final experiment on the effects of model misspecification, we introduce
non-linearity into the DGP while model fitting still assumes the conditional
linearity in (1.4). More precisely, we modify (1.4) by adding the softplus
response function $s_{c}(x)=c\,\ln\big{(}1+\exp(x/c)\big{)}$ with adjustment
parameter $c>0$, in analogy to the class of softplus INGARCH models proposed
by Weiß et al. (2022):

$\textstyle M_{t}\ =\
s_{c}\Big{(}a_{0}+\sum_{i=1}^{p}a_{i}\,X_{t-i}+\sum_{j=1}^{q}b_{j}\,M_{t-j}\Big{)}.$
(5.1)

Increasing $c$ increases the extent of non-linearity, whereas $\lim_{c\to
0}s_{c}(x)=\max\\{0,x\\}$ is piecewise linear. We focus on the illustrative
choice $c=2$ (for the default choice $c=1$, we hardly observed any effect on
estimates and diagnostics), where the results are summarized in Tables
S.11–S.14. Comparing Table S.11 to S.1 and Table S.12 to S.2, we note a slight
effect of the model misspecification on the estimates for all parameters, and
with the strongest effect again on $\sigma^{2}$. For practice, however, the
main question is the achieved model performance. In Tables S.13 and S.14, the
MAR and MSPR performance, respectively, is analyzed in comparison to that of
the truly linear DGPs. For the MSPR values in Table S.14, there is an at most
tiny effect if DGPs with Poi-counting series are considered, whereas the
increase in MSPR is slightly more pronounced in the Bin-case. The main
difference, however, is observed in the MAR performance, see Table S.13,
namely clearly increased MAR values under neglected non-linearity. This is
plausible as the softplus response function directly affects the conditional
mean $M_{t}$, which, in turn, is the crucial quantity when computing the MAR,
see Section 4.3.

###  6 Real-World Data Examples

(a)$t$ (b) $k$

Figure 2: Ecoli counts from Section 6.1: (a) time series plot and (b) sample
ACF against lag $k$.

####  6.1 Disease Counts

As our first data example, we use one of the data sets analyzed by Aknouche
and Scotto (2023), as these authors proposed another multiplicative model for
counts, recall Remark 3.1.1. More precisely, we consider the weekly counts of
disease cases caused by Escherichia coli (Ecoli), as reported for North Rhine-
Westphalia (Germany) from January 2001 to May 2013. These data are provided
through the command ecoli of the R-package tscount, see Liboschik et al.
(2017), and they were originally taken from SurvStat@RKI 2.0 at
https://survstat.rki.de/. Plots of this time series, $x_{1},\ldots,x_{n}$ of
length $n=646$, as well as of its sample ACF are shown in Figure 2. We
recognize a slow decay of the ACF, and also the partial ACF (not displayed)
exhibits a gradual fading. This supports the recommendation of Aknouche and
Scotto (2023) to use an INGARCH$(1,1)$-type model for these data. Initial
model fitting is done by the method of moments (MM), in analogy to Section 5.
From sample mean and ACF at lags 1–2, recall Example 3.1.5, we get the
estimates
$(\hat{a}_{0,\textup{MM}},\hat{a}_{1,\textup{MM}},\hat{b}_{1,\textup{MM}})\approx(2.465,0.431,0.448)$.
These can be used to compute the scaled residuals
$S_{t}(\hat{}\mbox{\boldmath$\theta$})=X_{t}/M_{t}(\hat{}\mbox{\boldmath$\theta$})$,
recall Section 4.3, the sample variance of which is given by $\approx 0.121$.
As this value is much smaller than one, we conclude that an
INGARCH$(1,1)$-CMEM with NB-counting series is not suitable for the Ecoli
counts, so we restrict the remaining model fitting to the cases of a Poi- and
Bin-counting series. Using the approach of Section 4.1.2, the corresponding
initial estimates of $\sigma^{2}$ are $0.068$ (Poi-case) and $0.120$ (Bin-
case), respectively.

Table 1: Ecoli counts from Section 6.1: Estimation results from different methods, approximate SEs in parentheses. |  Poi-counting series (3.5) |  Bin-counting series (3.15)
---|---|---
| $a_{0}$ | $a_{1}$ | $b_{1}$ |  $\sigma^{2}$ | $a_{0}$ | $a_{1}$ | $b_{1}$ | $\sigma^{2}$
PQ | 2.887 | 0.378 | 0.481 | 0.063 | 2.887 | 0.378 | 0.481 | 0.115
| (0.620) | (0.040) | (0.055) | (0.012) | (0.649) | (0.043) | (0.057) | (0.012)
NQ | 3.054 | 0.337 | 0.512 | 0.063 | 3.054 | 0.337 | 0.512 | 0.115
| (0.616) | (0.038) | (0.055) | (0.012) | (0.577) | (0.037) | (0.052) | (0.012)
WQ | 3.081 | 0.336 | 0.511 | 0.063 | 3.081 | 0.336 | 0.511 | 0.114
| (0.626) | (0.038) | (0.055) | (0.012) | (0.580) | (0.037) | (0.053) | (0.012)
2W | 2.938 | 0.351 | 0.505 | 0.063 | 3.084 | 0.339 | 0.508 | 0.114
| (0.590) | (0.038) | (0.053) | (0.012) | (0.581) | (0.037) | (0.053) | (0.012)

Using these initial estimates, the final estimates according to the methods
PQ, NQ, EQ, and 2W are computed, see Table 1 for an overview. Note that the
methods PQ, NQ, and EQ lead to the same estimates for $a_{0},a_{1},b_{1}$ for
both models (3.5) and (3.15), whereas the estimates of $\sigma^{2}$ and the
approximate SEs differ. The method 2W, by contrast, leads to different
estimates throughout. According to Table 1, the NQ- and EQ-estimates are quite
similar to each other for a given type of model (and for (3.15), they are also
close to the 2W-estimates), but they differ visibly from the PQ-estimates.
This goes along with our findings from Section 5, where we noted the same
(dis)agreements in Tables S.3 and S.4 (right blocks, as $\sigma^{2}$ is
estimated to about 0.1). Also from these tables, recalling that $n=646$ for
the Ecoli counts, we expect the estimates of $a_{0},a_{1},b_{1}$ to be biased.
For $\sigma^{2}$, by contrast, we do not expect a notable bias. Its estimate
is larger for model (3.15) than for (3.5), which is plausible as the Bin-
counting series causes less variation than the Poi-one, so more variation has
to be caused by the innovations $\epsilon_{t}$.

Table 2: Ecoli counts from Section 6.1: Mean, variance, and ACF at lags 1–5:
sample properties vs. properties of fitted models (with Poi- or Bin-counting
series).

| Mean | Var | $\rho(1)$ | $\rho(2)$ | $\rho(3)$ | $\rho(4)$ | $\rho(5)$
---|---|---|---|---|---|---|---
Sample | 20.334 | 88.753 | 0.632 | 0.555 | 0.471 | 0.426 | 0.343
Poi, PQ | 20.462 | 75.261 | 0.547 | 0.470 | 0.404 | 0.347 | 0.298
NQ | 20.223 | 66.333 | 0.485 | 0.412 | 0.350 | 0.297 | 0.252
EQ | 20.214 | 65.992 | 0.483 | 0.409 | 0.347 | 0.294 | 0.249
2W | 20.306 | 69.539 | 0.509 | 0.435 | 0.372 | 0.319 | 0.272
Bin, PQ | 20.462 | $[79.56,79.97]$ | 0.547 | 0.470 | 0.404 | 0.347 | 0.298
NQ | 20.223 | $[69.13,69.50]$ | 0.485 | 0.412 | 0.350 | 0.297 | 0.252
EQ | 20.214 | $[68.74,69.10]$ | 0.483 | 0.409 | 0.347 | 0.294 | 0.249
2W | 20.223 | $[69.20,69.57]$ | 0.487 | 0.413 | 0.350 | 0.296 | 0.251

Table 3: Ecoli counts from Section 6.1: model diagnostics according to Section
4.3, applied to CMEMs and some competing models.

|  Poi-counting series (3.5) |  Bin-counting series (3.15) |  Model (3.7)
---|---|---|---
| MAR | MSR | VSR | MSPR | MAR | MSR | VSR | MSPR | MAR | MSR | VSR | MSPR
PQ | 5.154 | 1.000 | 0.116 | 0.989 | 5.154 | 1.000 | 0.116 | 1.000 | 5.166 | 1.051 | 0.215 | 0.999
NQ | 5.143 | 1.000 | 0.115 | 0.995 | 5.143 | 1.000 | 0.115 | 1.000 | 5.150 | 1.054 | 0.218 | 1.004
EQ | 5.143 | 1.000 | 0.115 | 0.995 | 5.143 | 1.000 | 0.115 | 1.000 | 5.150 | 1.053 | 0.217 | 1.004
2W | 5.145 | 1.000 | 0.115 | 0.992 | 5.144 | 1.000 | 0.115 | 1.000 | 5.154 | 1.053 | 0.216 | 1.002
|  Poi-INGARCH$(1,1)$ |  NB-INGARCH$(1,1)$ |  |  |  |
CML | 5.154 | 1.000 | 0.116 | 2.267 | 5.144 | 1.000 | 0.116 | 1.035 |  |  |  |
Note: Italic values taken from Aknouche and Scotto (2023, Tables III–IV).

To check the adequacy of the fitted models, we consider the approaches
discussed in Section 4.3. In Table 3, marginal mean and variance as well as
the ACF at lags 1–5 of the fitted models are compared with the respective
sample values. Here, it should be recalled that for the Bin-model (3.15), we
do not have a closed-form expression, but only a lower and upper bound, see
Section 3.2. But the resulting intervals in the lower part of Table 3 are
quite narrow and, thus, provide reasonable insight into the true variance of
these models. While the means agree very well, the fitted models’ variances
are somewhat lower than the sample variance, where a closer agreement is
reached for the PQ-fitted model. The fitted models’ ACF is slowly decaying
like the sample ACF, but its actual values are too low (again, PQ gives a
better fit). The statistics in Table 3 have also been considered by Aknouche
and Scotto (2023, Tables III–IV) for their type of multiplicative count model
(3.7), recall Remark 3.1.1. Therefore, we provide the corresponding values for
MAR, MSR, VSR, and MSPR from Aknouche and Scotto (2023, Tables III–IV) within
Table 3 for comparison. As further competitors, we fit the (fully
parametrized) Poi-INGARCH$(1,1)$ model of Ferland et al. (2006) and NB-
INGARCH$(1,1)$ model of Zhu (2011) by conditional maximum likelihood (CML)
estimation, see the lower part of Table 3. The MARs in Table 3 are lower than
the MARs for the respective fit of model (3.7), and we get lower MARs for NQ,
EQ, 2W than for PQ. The MSPRs for model (3.5) are again closer to the target
value of 1 if using NQ-, EQ- or 2W-estimates, but a nearly perfect agreement
is achieved for the Bin-model (3.15). Again, we observe a better MSPR
performance than Aknouche and Scotto (2023) did for their model fits of (3.7).
Regarding the scaled residuals, the MSR is extremely close to the target value
1 without exception. The VSR values do not differ much between the different
model fits and are close to the estimate of $\sigma^{2}$ within the Bin-model
(3.15). For the fully parametrized INGARCH$(1,1)$ models, the Poi-
INGARCH$(1,1)$ model clearly falls short in terms of MSPR, because it is not
able to capture the (conditional) variance of the Ecoli counts. The NB-
INGARCH$(1,1)$ model is more powerful in this respect, but its MSPR value
still shows stronger deviation from the target value of 1. So altogether, the
CMEMs show the best overall performance.

To sum up, while the unconditional moments are better met if using PQ-
estimates, the conditional properties are better captured by NQ, EQ, and 2W
for the Ecoli data. Overall, the model (3.15) with Bin-counting series
provides a better fit, but with visible deviations regarding marginal variance
and ACF. A reason might be the period $t=544,\ldots,550$, where exceptionally
large count values $x_{t}$ are observed (larger than 50). This segment might
affect the obtained estimates (as well as the printed sample properties)
beyond the “natural” bias that is to be expected for an INGARCH$(1,1)$-CMEM
time series of length $n=646$.

Table 4: Ecoli counts from Section 6.1: predictive performance for CMEMs and some competing models. |  Poi-counting series (3.5) |  Bin-counting series (3.15)
---|---|---
| MAR | MSR | VSR | MSPR | MAR | MSR | VSR | MSPR
PQ | 6.249 | 1.066 | 0.124 | 1.136 | 6.249 | 1.066 | 0.124 | 1.092
NQ | 6.323 | 1.075 | 0.124 | 1.167 | 6.323 | 1.075 | 0.124 | 1.115
EQ | 6.329 | 1.075 | 0.124 | 1.169 | 6.329 | 1.075 | 0.124 | 1.117
2W | 6.278 | 1.071 | 0.124 | 1.152 | 6.323 | 1.075 | 0.124 | 1.116
|  Poi-INGARCH$(1,1)$ |  NB-INGARCH$(1,1)$
CML | 6.249 | 1.066 | 0.124 | 2.812 | 6.272 | 1.070 | 0.123 | 1.201

Finally, let us analyze the forecast performance of the fitted CMEMs in
comparison to the competing Poi- and NB-INGARCH$(1,1)$ model. In Table 3, we
applied the MAR, MSR, VSR, and MSPR to the same data that was also used for
model fitting, i. e., they provide insight into the goodness of fit. Next, we
apply these criteria to 84 new Ecoli counts, taken from SurvStat@RKI 2.0 at
https://survstat.rki.de/ for the remainder of 2013 as well as for the year
2014, to get insight into the predictive abilities of the fitted models.
Results are summarized in Table 4. While all models perform similar in terms
of MSR and VSR, the best MAR performance is achieved for the PQ-fit of both
types of CMEM as well as for the Poi-INGARCH$(1,1)$ model. The latter,
however, again falls short in terms of the MSPR. The NB-INGARCH$(1,1)$ model
has a MSPR value much closer to 1, but worse than for any of the CMEMs. Among
the CMEMs, the models with Bin-counting series have a better MSPR performance
than their Poi-counterparts, and the MSPR closest to 1 is achieved again for
the PQ-fit. Hence, the PQ-fitted CMEM (3.15) with Bin-counting series shows
the best predictive performance regarding the 84 Ecoli counts from June 2013
until December 2014.

(a)$t$ (b) $k$

Figure 3: WPP counts from Section 6.2: (a) time series plot and (b) sample ACF
against lag $k$.

####  6.2 Transaction Counts

Besides the surveillance of infectious diseases, another common application
area for count time series is financial transactions. Here, we consider one of
the time series provided by Aknouche et al. (2022), namely the number of stock
transactions concerning the Wausau Paper Corporation (WPP), measured in 5-min
intervals between 9:45 AM and 4:00 PM for the period from January 3 to
February 18 in 2005. More precisely, in analogy to Aknouche et al. (2022,
Section 6.2), we leave out the last 100 observations for a later forecast
experiment. The resulting count time series $x_{1},\ldots,x_{n}$ of length
$n=2825$ is plotted in Figure 3 together with its sample ACF. The ACF is
slowly decaying, which indicates that an INGARCH$(1,1)$-type model with a
large value for the feedback parameter $b_{1}$ could be appropriate for these
data, also see Aknouche et al. (2022). In analogy to Section 6.1, we start
with an initial analysis of the WPP counts. The moment estimates regarding the
conditional mean (1.4) are
$(\hat{a}_{0,\textup{MM}},\hat{a}_{1,\textup{MM}},\hat{b}_{1,\textup{MM}})\approx(1.874,0.361,0.409)$,
which are used to compute the scaled residuals and the sample variance
thereof. The obtained value $\approx 0.437$ is again much smaller than one, so
an NB-counting series is not suitable and we restrict the following model
fitting to INGARCH$(1,1)$-CMEMs with Poi- and Bin-counting series. The
corresponding initial estimates of $\sigma^{2}$ are $0.299$ (Poi-case) and
$0.433$ (Bin-case), respectively.

Table 5: WPP counts from Section 6.2: Estimation results from different methods, approximate SEs in parentheses. |  Poi-counting series (3.5) |  Bin-counting series (3.15)
---|---|---
| $a_{0}$ | $a_{1}$ | $b_{1}$ |  $\sigma^{2}$ | $a_{0}$ | $a_{1}$ | $b_{1}$ | $\sigma^{2}$
PQ | 0.792 | 0.268 | 0.634 | 0.301 | 0.792 | 0.268 | 0.634 | 0.439
| (0.108) | (0.018) | (0.026) | (0.022) | (0.110) | (0.019) | (0.027) | (0.022)
NQ | 0.768 | 0.266 | 0.640 | 0.302 | 0.768 | 0.266 | 0.640 | 0.439
| (0.103) | (0.017) | (0.025) | (0.022) | (0.100) | (0.017) | (0.025) | (0.022)
WQ | 0.769 | 0.265 | 0.640 | 0.301 | 0.769 | 0.265 | 0.640 | 0.439
| (0.104) | (0.017) | (0.025) | (0.022) | (0.100) | (0.017) | (0.025) | (0.022)
2W | 0.769 | 0.266 | 0.640 | 0.302 | 0.769 | 0.265 | 0.640 | 0.439
| (0.103) | (0.017) | (0.025) | (0.022) | (0.100) | (0.017) | (0.025) | (0.022)

The final estimates together with their approximate SEs are summarized in
Table 5. Because of the very large sample size, namely $n=2825$, we do not
expect any bias problems in estimation, recall our simulation study in Section
5. In fact, all estimation methods lead to very similar results. The only
visible deviations are observed for the intercept parameter $a_{0}$ on the one
hand, where PQ leads to a slightly larger estimate than the other methods, and
for $\sigma^{2}$ on the other hand, which has larger estimates for the Bin-
than for the Poi-counting series (for the same reason as explained in Section
6.1). It is worth mentioning that our final estimates are close to the ones
obtained by Aknouche et al. (2022, Supplement 1.2), although these authors
used fully parametric INGARCH$(1,1)$ models. In particular, we note a very
close agreement to their “NB2-INGARCH model” in Table 3, which is the NB-
INGARCH$(1,1)$ model of Zhu (2011) and constitutes their recommended model
choice.

Table 6: WPP counts from Section 6.2: Mean, variance, and ACF at lags 1–5:
sample properties vs. properties of fitted models (with Poi- or Bin-counting
series).

| Mean | Var | $\rho(1)$ | $\rho(2)$ | $\rho(3)$ | $\rho(4)$ | $\rho(5)$
---|---|---|---|---|---|---|---
Sample | 8.144 | 35.976 | 0.460 | 0.354 | 0.327 | 0.309 | 0.273
Poi, PQ | 8.150 | 44.375 | 0.446 | 0.403 | 0.364 | 0.329 | 0.297
NQ | 8.164 | 44.697 | 0.446 | 0.404 | 0.366 | 0.332 | 0.301
EQ | 8.163 | 44.613 | 0.446 | 0.404 | 0.366 | 0.331 | 0.300
2W | 8.166 | 44.786 | 0.447 | 0.405 | 0.367 | 0.332 | 0.301
Bin, PQ | 8.150 | $[48.92,49.34]$ | 0.446 | 0.403 | 0.364 | 0.329 | 0.297
NQ | 8.164 | $[49.34,49.76]$ | 0.446 | 0.404 | 0.366 | 0.332 | 0.301
EQ | 8.163 | $[49.23,49.65]$ | 0.446 | 0.404 | 0.366 | 0.331 | 0.300
2W | 8.161 | $[49.24,49.66]$ | 0.446 | 0.404 | 0.366 | 0.331 | 0.300

Table 7: WPP counts from Section 6.2: model diagnostics according to Section
4.3, applied to CMEMs and some competing models.

|  Poi-counting series (3.5) |  Bin-counting series (3.15)
---|---|---
| MAR | MSR | VSR | MSPR | MAR | MSR | VSR | MSPR
PQ | 3.916 | 1.000 | 0.443 | 0.975 | 3.916 | 1.000 | 0.443 | 0.998
NQ | 3.916 | 1.000 | 0.443 | 0.975 | 3.916 | 1.000 | 0.443 | 0.998
EQ | 3.916 | 1.000 | 0.443 | 0.975 | 3.916 | 1.000 | 0.443 | 0.998
2W | 3.916 | 1.000 | 0.443 | 0.975 | 3.916 | 1.000 | 0.443 | 0.998
|  Poi-INGARCH$(1,1)$ |  NB-INGARCH$(1,1)$
CML | 3.916 | 1.000 | 0.443 | 3.181 | 3.916 | 1.000 | 0.444 | 1.036

Thus, in analogy to Section 6.1, we also fit the (fully parametrized) Poi-
INGARCH$(1,1)$ model of Ferland et al. (2006) and NB-INGARCH$(1,1)$ model of
Zhu (2011) as further competitors. Looking at the adequacy checks in Tables
7–7, we generally recognize a very good performance of the CMEMs with only
minor deviations between the different estimation methods. Table 7 shows that
the fitted CMEMs’ mean and ACF closely agree with the corresponding sample
properties. Only for the marginal variance, we observe an exceedance for the
fitted models, which is less pronounced in case of the Poi-counting series
(3.5). The opposite ranking is implied by the MSPR values in Table 7, i. e.,
the conditional variance is better captured by the Bin-counting series (3.15).
But in any case, the differences between the different models are rather
small. Compared with the CML-fitted Poi- and NB-INGARCH$(1,1)$ model, we note
that the Poi-INGARCH$(1,1)$ model fails again in terms of the MSPR, which is
much larger than 1. The NB-INGARCH$(1,1)$ model does clearly better, but still
shows a stronger deviation of the MSPR from 1 than any of the CMEMs.
Altogether, any of the fitted CMEMs performs reasonably well, where the Poi-
counting series (3.5) is advantageous regarding the marginal variance, and the
Bin-counting series (3.15) regarding the conditional variance.

Table 8: WPP counts from Section 6.2: predictive performance for CMEMs and some competing models. |  Poi-counting series (3.5) |  Bin-counting series (3.15)
---|---|---
| MAR | MSR | VSR | MSPR | MAR | MSR | VSR | MSPR
PQ | 3.613 | 1.026 | 0.552 | 1.164 | 3.613 | 1.026 | 0.552 | 1.232
NQ | 3.616 | 1.026 | 0.553 | 1.165 | 3.616 | 1.026 | 0.553 | 1.233
EQ | 3.616 | 1.026 | 0.553 | 1.165 | 3.616 | 1.026 | 0.553 | 1.233
2W | 3.616 | 1.026 | 0.553 | 1.166 | 3.616 | 1.026 | 0.553 | 1.233
|  Poi-INGARCH$(1,1)$ |  NB-INGARCH$(1,1)$
CML | 3.610 | 1.025 | 0.551 | 3.400 | 3.618 | 1.027 | 0.554 | 1.236

Finally, like in Section 6.1, we do a forecast experiment based on the WPP’s
last 100 observations that were left out during model fitting, see Table 8 for
the results. The CMEMs have rather similar MAR, MSR, and VSR values, but the
MSPR of the CMEMs with Poi-counting series is closer to 1, with the best
predictive performance for the PQ-fitted CMEM. Regarding the competitors, the
Poi-INGARCH$(1,1)$ model fails again in terms of the MSPR, while the NB-
INGARCH$(1,1)$ model performs slightly worse than the preferred CMEM with
respect to all four criteria. Hence, the PQ-fitted CMEM (3.5) with Poi-
counting series shows the best predictive performance regarding the last 100
WPP counts.

###  7 Conclusions and Future Research

In this article, we proposed a novel framework for transferring the well-
established MEMs for real-valued time series to the case of count time series.
The proposed CMEMs make use of an appropriate multiplicative operator for
counts. After discussing such operators and the corresponding CMEMs in
general, we focused on two special cases, namely a compounding operator using,
e. g., a Poi-counting series, and a binomial multiplicative operator. We
established the existence of stationary and ergodic solutions to such INGARCH-
CMEMs and derived relevant moment properties. Afterwards, we developed two
approaches for semi-parametric model estimation, namely QMLE and 2WLSE. The
finite-sample performance of these estimation approaches was analyzed in a
simulation study, and the applicability of our novel CMEMs was demonstrated
with two real-world data examples.

There are several directions for future research. First, CMEMs using another
multiplicative operator (such as the ZIP-multiplicative operator mentioned in
the end of Section 2) or a non-linear mean specification might be analyzed.
Second, the semi-parametric INGARCH model sketched in the end of Section 3.1
deserves further investigation, e. g., with respect to the question if a semi-
parametric maximum likelihood (ML) approach, analogous to Drost et al. (2009)
and Liu et al. (2021), is possible for INGARCH models. Such an ML-fitted CMEM
would be attractive for applications as it would allow, for example, for
density forecasts like in Aknouche et al. (2022). Third, related to the
previous research direction, it could also be studied if a combination of
semi-parametric ML estimation with the smoothing method discussed by
Faymonville et al. (2023) could be beneficial for practice. Finally,
$\mathbb{Z}$-valued GARCH models also get attention, see Cui et al. (2021) and
Xu and Zhu (2022) among others, so how to generalize the CMEM framework to the
$\mathbb{Z}$-valued case is interesting.

##### Acknowledgements

The authors thank the associate editor and the two referees for their useful
and constructive comments on an earlier draft of this article. Zhu’s work is
supported by National Natural Science Foundation of China (No. 12271206) and
Natural Science Foundation of Jilin Province (No. 20210101143JC), and the
Science and Technology Research Planning Project of Jilin Provincial
Department of Education (No. JJKH20231122KJ).

### References

  * Agosto et al. (2016) Agosto, A., Cavaliere, G., Kristensen, D., Rahbek, A., 2016. Modeling corporate defaults: Poisson autoregressions with exogenous covariates (PARX). Journal of Empirical Finance 38, 640–663.
  * Ahmad and Francq (2016) Ahmad, A., Francq, C., 2016. Poisson QMLE of count time series models. Journal of Time Series Analysis 37, 291–314.
  * Aknouche (2012) Aknouche, A., 2012. Multistage weighted least squares estimation of ARCH processes in the stable and unstable cases. Statistical Inference for Stochastic Processes 15, 241–256.
  * Aknouche et al. (2022) Aknouche, A., Almohaimeed, B.S., Dimitrakopoulos, S., 2022. Forecasting transaction counts with integer-valued GARCH models. Studies in Nonlinear Dynamics & Econometrics 26, 529–539.
  * Aknouche et al. (2018) Aknouche, A., Bendjeddou, S., Touche, N., 2018. Negative binomial quasi-likelihood inference for general integer-valued time series models. Journal of Time Series Analysis 39, 192–211.
  * Aknouche and Francq (2021) Aknouche, A., Francq, C., 2021. Count and duration time series with equal conditional stochastic and mean orders. Econometric Theory 37, 248–280.
  * Aknouche and Scotto (2023) Aknouche, A., Scotto, M.G., 2023. A multiplicative thinning-based integer-valued GARCH model. Journal of Time Series Analysis, forthcoming, https://doi.org/10.1111/jtsa.12682
  * Aleksandrov and Weiß (2020) Aleksandrov, B., Weiß, C.H., 2020. Testing the dispersion structure of count time series using Pearson residuals. AStA Advances in Statistical Analysis 104, 325–361.
  * Billingsley (1995) Billingsley, P., 1995. Probability and Measure (3rd ed.) John Wiley & Sons, New York.
  * Bollerslev (1986) Bollerslev, T., 1986. Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics 31, 307–327.
  * Brownlees et al. (2012) Brownlees, C.T., Cipollini, F., Gallo, G.M., 2012. Multiplicative error models. Chapter 9 in Bauwens et al. (eds.): Handbook of Volatility Models and Their Applications, John Wiley & Sons, Inc., 223–247.
  * Chen et al. (2022) Chen, H., Li, Q., Zhu, F., 2022. A new class of integer-valued GARCH models for time series of bounded counts with extra-binomial variation. AStA Advances in Statistical Analysis 106, 243–270.
  * Christou and Fokianos (2014) Christou, V., Fokianos, K., 2014. Quasi-likelihood inference for negative binomial time series models. Journal of Time Series Analysis 35, 55–78.
  * Cipollini and Gallo (2023) Cipollini, F., Gallo, G.M., 2023. Multiplicative error models: 20 years on. Econometrics and Statistics, forthcoming, https://doi.org/10.1016/j.ecosta.2022.05.005.
  * Cui et al. (2021) Cui, Y., Li, Q., Zhu, F., 2021. Modeling $\mathbb{Z}$-valued time series based on new versions of the Skellam INGARCH model. Brazilian Journal of Probability and Statistics 35, 293–314.
  * Doukhan et al. (2021) Doukhan, P., Fokianos, K. and Rynkiewicz, J., 2021. Mixtures of Nonlinear Poisson Autoregressions. Journal of Time Series Analysis 42, 107–135.
  * Doukhan and Neumann (2019) Doukhan, P., Neumann, M.H., 2019. Absolute regularity of semi-contractive GARCH-type processes. Journal of Applied Probability 56, 91–115.
  * Doukhan and Wintenberger (2008) Doukhan, P., Wintenberger, O., 2008. Weakly dependent chains with infinite memory. Stochastic Processes and their Applications 118, 1997–2013.
  * Drost et al. (2009) Drost, F.C., van den Akker, R., Werker, B.J.M., 2009. Efficient estimation of auto-regression parameters and innovation distributions for semiparametric integer-valued AR(p) models. Journal of the Royal Statistical Society, Series B 71, 467–485.
  * Engle (2002) Engle, R.F., 2002. New frontiers for ARCH models. Journal of Applied Econometrics 17, 425–446.
  * Engle & Russell (1998) Engle, R.F., Russell, J.R., 1998. Autoregressive conditional duration: a new model for irregularly spaced transaction data. Econometrica 66, 1127–1162.
  * Faymonville et al. (2023) Faymonville, M., Jentsch, C., Weiß, C.H., Aleksandrov, B., 2023. Semiparametric estimation of INAR models using roughness penalization. Statistical Methods and Applications 32, 365–400.
  * Ferland et al. (2006) Ferland, R., Latour, A., Oraichi, D., 2006. Integer-valued GARCH processes. Journal of Time Series Analysis 27, 923–942.
  * Fokianos et al. (2009) Fokianos, K., Rahbek, A., Tjøstheim, D., 2009. Poisson autoregression. Journal of the American Statistical Association 104, 1430–1439.
  * Fokianos & Tjøstheim (2011) Fokianos, K., Tjøstheim, D., 2011. Log-linear Poisson autoregression. Journal of Multivariate Analysis 102, 563–578.
  * Gonçalves et al. (2015a) Gonçalves, E., Mendes-Lopes, N., Silva, F., 2015a. Infinitely divisible distributions in integer-valued GARCH models. Journal of Time Series Analysis 36, 503–527.
  * Gonçalves et al. (2015b) Gonçalves, E., Mendes-Lopes, N., Silva, F., 2015b. A new approach to integer-valued time series modeling: the Neyman type-A INGARCH model. Lithuanian Mathematical Journal 55, 231–242.
  * Gonçalves et al. (2016) Gonçalves, E., Mendes-Lopes, N., Silva, F., 2016. Zero-inflated compound Poisson distributions in integer-valued GARCH models. Statistics 50, 558–578.
  * Heinen (2003) Heinen, A., 2003. Modelling time series count data: an autoregressive conditional Poisson model. CORE Discussion Paper 2003/62, University of Louvain, Belgium.
  * Kachour and Yao (2009) Kachour, M., Yao, J.F., 2009. First-order rounded integer-valued autoregressive (RINAR(1)) process. Journal of Time Series Analysis 30, 417–448.
  * Klenke and Mattner (2010) Klenke, A., Mattner, L., 2010. Stochastic ordering of classical discrete distributions. Advances in Applied Probability 42, 392–410.
  * Latour (1998) Latour, A., 1998. Existence and stochastic structure of a non-negative integer-valued autoregressive process. Journal of Time Series Analysis 19, 439–455.
  * Liboschik et al. (2017) Liboschik, T., Fokianos, K., Fried, R., 2017. tscount: An R package for analysis of count time series following generalized linear models. Journal of Statistical Software 82, 1–51.
  * Liu et al. (2021) Liu, Z., Li, Q., Zhu, F., 2021. Semiparametric integer-valued autoregressive models on $\mathbb{Z}$. Canadian Journal of Statistics 49, 1317–1337.
  * Mao et al. (2020) Mao, H., Zhu, F., Cui, Y., 2020. A generalized mixture integer-valued GARCH model. Statistical Methods and Applications 29, 527–552.
  * Ristić et al. (2009) Ristić, M.M., Bakouch, H.S., Nastić, A.S., 2009. A new geometric first-order integer-valued autoregressive (NGINAR$(1)$) process. Journal of Statistical Planning and Inference 139, 2218–2226.
  * Scotto et al. (2015) Scotto, M.G., Weiß, C.H., Gouveia, S., 2015. Thinning-based models in the analysis of integer-valued time series: a review. Statistical Modelling 15, 590–618.
  * Steutel & van Harn (1979) Steutel, F.W., van Harn, K., 1979. Discrete analogues of self-decomposability and stability. Annals of Probability 7, 893–899.
  * Weiß (2008) Weiß, C.H., 2008. Thinning operations for modelling time series of counts — a survey. AStA Advances in Statistical Analysis 92, 319–341.
  * Weiß (2009) Weiß, C.H., 2009. Modelling time series of counts with overdispersion. Statistical Methods and Applications 18, 507–519.
  * Weiß (2018) Weiß, C.H., 2018. An Introduction to Discrete-Valued Time Series. John Wiley & Sons, Chichester.
  * Weiß et al. (2017) Weiß, C.H., Gonçalves, E., Mendes Lopes, N., 2017. Testing the compounding structure of the CP-INARCH model. Metrika 80, 571–603.
  * Weiß et al. (2022) Weiß, C.H., Zhu, F., Hoshiyar, A., 2022. Softplus INGARCH models. Statistica Sinica 32, 1099–1120.
  * Wooldridge (1999) Wooldridge, J.M., 1999. Quasi-likelihood methods for count data. Chapter 8 in Pesaran H, Schmidt P. (eds.): Handbook of Applied Econometrics, Volume 2, Microeconomics, Blackwell, 321–368.
  * Xu and Zhu (2022) Xu, Y., Zhu, F., 2022. A new GJR-GARCH model for $\mathbb{Z}$-valued time series. Journal of Time Series Analysis 43, 490–500.
  * Zhu (2011) Zhu, F., 2011. A negative binomial integer-valued GARCH model. Journal of Time Series Analysis 32, 54–67.
  * Zhu et al. (2010) Zhu, F., Li, Q., Wang, D., 2010. A mixture integer-valued ARCH model. Journal of Statistical Planning and Inference 140, 2025–2036.
  * Zhu and Wang (2010) Zhu, F., Wang, D., 2010. Diagnostic checking integer-valued ARCH$(p)$ models using conditional residual autocorrelations. Computational Statistics and Data Analysis 54, 496–508.

## Supplementary Material:

### Conditional-mean Multiplicative Operator Models
for Count Time Series

###  Appendix S Detailed Simulation Results

In what follows, we present tables with detailed simulation results, which
correspond to the simulation study presented in Section 5 of the main
manuscript. The results rely on 10,000 replications per scenario. To achieve
robustness with respect to possible outliers, we used the trimmed versions of
sample mean and standard deviation, where the smallest and largest 0.1 % of
simulated values are omitted during calculation. We observed such outliers
only in very few simulation runs, but these would have a strongly distorting
influence on ordinary sample mean and standard deviation.

##### Acronyms

Bin = binomial

NB = negative binomial

Poi = Poisson

DGP = data-generating process

CMEM = count multiplicative error model

INGARCH = integer-valued generalized autoregressive conditional
heteroscedasticity

QMLE = quasi-maximum likelihood estimation

PQ / NQ / EQ = Poi- / NB- / exponential QMLE

2W, 2SWLSE = two-stage weighted least squares estimation

Table S.1: Different estimation approaches applied to simulated count time
series from INGARCH$(1,1)$-CMEM with Poi-counting series: Mean, simulated
standard error (SSE), and mean approximated standard error (ASE) for different
sample sizes $n$.

$n$ |  | Param. | PQ | NQ | EQ | 2W | Param. | PQ | NQ | EQ | 2W
---|---|---|---|---|---|---|---|---|---|---|---
300 | Mean | $a_{0}=2.8$ | 2.997 | 2.922 | 2.915 | 2.935 | $a_{0}=2.8$ | 2.965 | 2.936 | 2.932 | 2.942
| SSE |  | 0.722 | 0.686 | 0.688 | 0.701 |  | 0.681 | 0.668 | 0.671 | 0.666
| ASE |  | 0.780 | 0.694 | 0.692 | 0.699 |  | 0.720 | 0.686 | 0.687 | 0.687
600 | Mean |  | 2.904 | 2.858 | 2.854 | 2.866 |  | 2.889 | 2.875 | 2.874 | 2.877
| SSE |  | 0.529 | 0.492 | 0.493 | 0.500 |  | 0.488 | 0.474 | 0.476 | 0.473
| ASE |  | 0.536 | 0.480 | 0.479 | 0.481 |  | 0.498 | 0.477 | 0.478 | 0.477
1000 | Mean |  | 2.871 | 2.839 | 2.837 | 2.843 |  | 2.855 | 2.844 | 2.843 | 2.845
| SSE |  | 0.406 | 0.375 | 0.376 | 0.378 |  | 0.380 | 0.366 | 0.367 | 0.365
| ASE |  | 0.412 | 0.369 | 0.369 | 0.370 |  | 0.382 | 0.366 | 0.367 | 0.366
300 | Mean | $a_{1}=0.4$ | 0.384 | 0.393 | 0.394 | 0.392 | $a_{1}=0.4$ | 0.392 | 0.395 | 0.396 | 0.395
| SSE |  | 0.081 | 0.077 | 0.077 | 0.079 |  | 0.062 | 0.060 | 0.060 | 0.060
| ASE |  | 0.082 | 0.074 | 0.074 | 0.074 |  | 0.061 | 0.059 | 0.059 | 0.059
600 | Mean |  | 0.392 | 0.397 | 0.397 | 0.397 |  | 0.397 | 0.398 | 0.399 | 0.398
| SSE |  | 0.057 | 0.053 | 0.053 | 0.053 |  | 0.044 | 0.042 | 0.043 | 0.042
| ASE |  | 0.059 | 0.052 | 0.052 | 0.052 |  | 0.043 | 0.042 | 0.042 | 0.042
1000 | Mean |  | 0.395 | 0.399 | 0.399 | 0.398 |  | 0.398 | 0.399 | 0.399 | 0.399
| SSE |  | 0.045 | 0.041 | 0.041 | 0.041 |  | 0.034 | 0.033 | 0.033 | 0.033
| ASE |  | 0.046 | 0.041 | 0.041 | 0.041 |  | 0.034 | 0.032 | 0.032 | 0.032
300 | Mean | $b_{1}=0.2$ | 0.180 | 0.184 | 0.185 | 0.183 | $b_{1}=0.2$ | 0.182 | 0.183 | 0.184 | 0.183
| SSE |  | 0.121 | 0.116 | 0.117 | 0.120 |  | 0.112 | 0.110 | 0.111 | 0.110
| ASE |  | 0.128 | 0.116 | 0.115 | 0.116 |  | 0.119 | 0.114 | 0.114 | 0.114
600 | Mean |  | 0.189 | 0.192 | 0.193 | 0.191 |  | 0.190 | 0.190 | 0.190 | 0.190
| SSE |  | 0.089 | 0.083 | 0.084 | 0.085 |  | 0.080 | 0.079 | 0.079 | 0.079
| ASE |  | 0.088 | 0.080 | 0.080 | 0.081 |  | 0.082 | 0.079 | 0.080 | 0.079
1000 | Mean |  | 0.192 | 0.194 | 0.195 | 0.194 |  | 0.194 | 0.194 | 0.194 | 0.194
| SSE |  | 0.068 | 0.064 | 0.064 | 0.064 |  | 0.063 | 0.061 | 0.062 | 0.061
| ASE |  | 0.068 | 0.062 | 0.062 | 0.062 |  | 0.063 | 0.061 | 0.061 | 0.061
300 | Mean | $\sigma^{2}=1$ | 1.001 | 0.998 | 0.998 | 1.002 | $\sigma^{2}=0.4$ | 0.398 | 0.397 | 0.397 | 0.397
| SSE |  | 0.114 | 0.113 | 0.113 | 0.130 |  | 0.050 | 0.050 | 0.050 | 0.050
| ASE |  | 0.128 | 0.127 | 0.127 | 0.128 |  | 0.046 | 0.046 | 0.046 | 0.046
600 | Mean |  | 1.002 | 0.999 | 0.999 | 1.000 |  | 0.400 | 0.399 | 0.399 | 0.399
| SSE |  | 0.080 | 0.080 | 0.080 | 0.081 |  | 0.035 | 0.035 | 0.035 | 0.035
| ASE |  | 0.092 | 0.091 | 0.091 | 0.091 |  | 0.033 | 0.033 | 0.033 | 0.033
1000 | Mean |  | 1.000 | 0.999 | 0.999 | 0.999 |  | 0.399 | 0.399 | 0.399 | 0.399
| SSE |  | 0.061 | 0.061 | 0.061 | 0.062 |  | 0.027 | 0.027 | 0.027 | 0.027
| ASE |  | 0.071 | 0.071 | 0.071 | 0.071 |  | 0.025 | 0.025 | 0.025 | 0.025

Table S.2: Different estimation approaches applied to simulated count time
series from INGARCH$(1,1)$-CMEM with Bin-counting series: Mean, simulated
standard error (SSE), and mean approximated standard error (ASE) for different
sample sizes $n$.

$n$ |  | Param. | PQ | NQ | EQ | 2W | Param. | PQ | NQ | EQ | 2W
---|---|---|---|---|---|---|---|---|---|---|---
300 | Mean | $a_{0}=2.8$ | 2.995 | 2.926 | 2.920 | 2.938 | $a_{0}=2.8$ | 2.967 | 2.942 | 2.939 | 2.942
| SSE |  | 0.713 | 0.667 | 0.667 | 0.676 |  | 0.655 | 0.629 | 0.629 | 0.629
| ASE |  | 0.768 | 0.673 | 0.670 | 0.673 |  | 0.693 | 0.650 | 0.648 | 0.649
600 | Mean |  | 2.913 | 2.869 | 2.865 | 2.875 |  | 2.892 | 2.879 | 2.878 | 2.878
| SSE |  | 0.520 | 0.477 | 0.477 | 0.483 |  | 0.479 | 0.454 | 0.453 | 0.453
| ASE |  | 0.531 | 0.467 | 0.466 | 0.466 |  | 0.480 | 0.451 | 0.451 | 0.451
1000 | Mean |  | 2.875 | 2.844 | 2.842 | 2.847 |  | 2.865 | 2.856 | 2.855 | 2.856
| SSE |  | 0.406 | 0.369 | 0.368 | 0.371 |  | 0.366 | 0.347 | 0.347 | 0.346
| ASE |  | 0.407 | 0.359 | 0.358 | 0.358 |  | 0.369 | 0.347 | 0.347 | 0.347
300 | Mean | $a_{1}=0.4$ | 0.384 | 0.393 | 0.394 | 0.393 | $a_{1}=0.4$ | 0.392 | 0.394 | 0.395 | 0.395
| SSE |  | 0.077 | 0.072 | 0.072 | 0.072 |  | 0.056 | 0.054 | 0.054 | 0.054
| ASE |  | 0.079 | 0.070 | 0.070 | 0.070 |  | 0.056 | 0.053 | 0.053 | 0.053
600 | Mean |  | 0.391 | 0.396 | 0.396 | 0.396 |  | 0.396 | 0.398 | 0.398 | 0.398
| SSE |  | 0.056 | 0.050 | 0.050 | 0.050 |  | 0.040 | 0.038 | 0.038 | 0.038
| ASE |  | 0.057 | 0.050 | 0.050 | 0.050 |  | 0.040 | 0.037 | 0.037 | 0.037
1000 | Mean |  | 0.395 | 0.398 | 0.398 | 0.398 |  | 0.398 | 0.399 | 0.399 | 0.399
| SSE |  | 0.043 | 0.038 | 0.038 | 0.038 |  | 0.030 | 0.029 | 0.029 | 0.029
| ASE |  | 0.044 | 0.039 | 0.039 | 0.039 |  | 0.031 | 0.029 | 0.029 | 0.029
300 | Mean | $b_{1}=0.2$ | 0.181 | 0.185 | 0.185 | 0.183 | $b_{1}=0.2$ | 0.181 | 0.183 | 0.183 | 0.182
| SSE |  | 0.118 | 0.112 | 0.112 | 0.114 |  | 0.106 | 0.103 | 0.103 | 0.103
| ASE |  | 0.124 | 0.112 | 0.111 | 0.111 |  | 0.112 | 0.106 | 0.106 | 0.106
600 | Mean |  | 0.189 | 0.192 | 0.192 | 0.191 |  | 0.189 | 0.190 | 0.190 | 0.190
| SSE |  | 0.086 | 0.080 | 0.080 | 0.081 |  | 0.078 | 0.075 | 0.075 | 0.075
| ASE |  | 0.086 | 0.078 | 0.078 | 0.078 |  | 0.078 | 0.074 | 0.074 | 0.074
1000 | Mean |  | 0.192 | 0.194 | 0.194 | 0.193 |  | 0.192 | 0.193 | 0.193 | 0.193
| SSE |  | 0.067 | 0.062 | 0.062 | 0.062 |  | 0.059 | 0.057 | 0.057 | 0.057
| ASE |  | 0.066 | 0.060 | 0.060 | 0.060 |  | 0.060 | 0.057 | 0.057 | 0.057
300 | Mean | $\sigma^{2}=1$ | 1.004 | 1.001 | 1.001 | 1.002 | $\sigma^{2}=0.4$ | 0.402 | 0.401 | 0.401 | 0.401
| SSE |  | 0.101 | 0.100 | 0.100 | 0.101 |  | 0.042 | 0.041 | 0.041 | 0.041
| ASE |  | 0.100 | 0.099 | 0.099 | 0.100 |  | 0.029 | 0.029 | 0.029 | 0.029
600 | Mean |  | 1.003 | 1.001 | 1.001 | 1.002 |  | 0.401 | 0.400 | 0.400 | 0.400
| SSE |  | 0.071 | 0.071 | 0.071 | 0.071 |  | 0.029 | 0.029 | 0.029 | 0.029
| ASE |  | 0.071 | 0.071 | 0.071 | 0.071 |  | 0.020 | 0.020 | 0.020 | 0.020
1000 | Mean |  | 1.003 | 1.001 | 1.001 | 1.001 |  | 0.400 | 0.400 | 0.400 | 0.400
| SSE |  | 0.054 | 0.053 | 0.053 | 0.054 |  | 0.022 | 0.022 | 0.022 | 0.022
| ASE |  | 0.055 | 0.055 | 0.055 | 0.055 |  | 0.016 | 0.016 | 0.016 | 0.016

Table S.3: Different estimation approaches applied to simulated count time
series from INGARCH$(1,1)$-CMEM with Poi-counting series: Mean, simulated
standard error (SSE), and mean approximated standard error (ASE) for different
sample sizes $n$.

$n$ |  | Param. | PQ | NQ | EQ | 2W | Param. | PQ | NQ | EQ | 2W
---|---|---|---|---|---|---|---|---|---|---|---
300 | Mean | $a_{0}=3$ | 3.723 | 3.551 | 3.553 | 3.579 | $a_{0}=3$ | 3.981 | 4.092 | 4.103 | 3.997
| SSE |  | 1.398 | 1.243 | 1.247 | 1.356 |  | 1.217 | 1.222 | 1.226 | 1.199
| ASE |  | 1.304 | 1.075 | 1.075 | 1.083 |  | 1.132 | 1.135 | 1.139 | 1.112
600 | Mean |  | 3.363 | 3.256 | 3.258 | 3.265 |  | 3.500 | 3.579 | 3.587 | 3.529
| SSE |  | 0.876 | 0.737 | 0.738 | 0.786 |  | 0.778 | 0.783 | 0.787 | 0.767
| ASE |  | 0.877 | 0.711 | 0.712 | 0.713 |  | 0.749 | 0.747 | 0.749 | 0.736
1000 | Mean |  | 3.238 | 3.166 | 3.168 | 3.171 |  | 3.290 | 3.354 | 3.361 | 3.318
| SSE |  | 0.647 | 0.553 | 0.555 | 0.567 |  | 0.574 | 0.580 | 0.582 | 0.566
| ASE |  | 0.671 | 0.539 | 0.539 | 0.539 |  | 0.563 | 0.560 | 0.562 | 0.553
300 | Mean | $a_{1}=0.35$ | 0.336 | 0.349 | 0.349 | 0.347 | $a_{1}=0.35$ | 0.361 | 0.365 | 0.366 | 0.363
| SSE |  | 0.073 | 0.066 | 0.066 | 0.068 |  | 0.050 | 0.050 | 0.050 | 0.050
| ASE |  | 0.079 | 0.066 | 0.066 | 0.066 |  | 0.051 | 0.050 | 0.051 | 0.050
600 | Mean |  | 0.342 | 0.349 | 0.350 | 0.349 |  | 0.358 | 0.361 | 0.361 | 0.360
| SSE |  | 0.054 | 0.046 | 0.046 | 0.047 |  | 0.035 | 0.034 | 0.034 | 0.034
| ASE |  | 0.058 | 0.046 | 0.046 | 0.046 |  | 0.036 | 0.036 | 0.036 | 0.035
1000 | Mean |  | 0.345 | 0.350 | 0.351 | 0.350 |  | 0.354 | 0.357 | 0.357 | 0.356
| SSE |  | 0.042 | 0.036 | 0.036 | 0.036 |  | 0.028 | 0.027 | 0.027 | 0.027
| ASE |  | 0.045 | 0.036 | 0.036 | 0.036 |  | 0.028 | 0.028 | 0.028 | 0.027
300 | Mean | $b_{1}=0.5$ | 0.468 | 0.468 | 0.467 | 0.468 | $b_{1}=0.5$ | 0.439 | 0.430 | 0.429 | 0.437
| SSE |  | 0.112 | 0.097 | 0.097 | 0.105 |  | 0.081 | 0.079 | 0.079 | 0.080
| ASE |  | 0.106 | 0.090 | 0.090 | 0.090 |  | 0.083 | 0.083 | 0.083 | 0.082
600 | Mean |  | 0.485 | 0.485 | 0.485 | 0.485 |  | 0.467 | 0.460 | 0.460 | 0.464
| SSE |  | 0.074 | 0.061 | 0.061 | 0.065 |  | 0.056 | 0.055 | 0.055 | 0.055
| ASE |  | 0.073 | 0.061 | 0.061 | 0.061 |  | 0.058 | 0.057 | 0.057 | 0.057
1000 | Mean |  | 0.490 | 0.490 | 0.490 | 0.490 |  | 0.481 | 0.475 | 0.475 | 0.478
| SSE |  | 0.055 | 0.046 | 0.046 | 0.048 |  | 0.043 | 0.043 | 0.043 | 0.042
| ASE |  | 0.056 | 0.047 | 0.047 | 0.047 |  | 0.044 | 0.044 | 0.044 | 0.043
300 | Mean | $\sigma^{2}=1$ | 1.018 | 1.008 | 1.008 | 1.019 | $\sigma^{2}=0.1$ | 0.101 | 0.100 | 0.100 | 0.101
| SSE |  | 0.121 | 0.113 | 0.112 | 0.159 |  | 0.020 | 0.020 | 0.020 | 0.020
| ASE |  | 0.118 | 0.114 | 0.114 | 0.120 |  | 0.020 | 0.020 | 0.020 | 0.020
600 | Mean |  | 1.011 | 1.005 | 1.005 | 1.009 |  | 0.101 | 0.100 | 0.100 | 0.100
| SSE |  | 0.082 | 0.078 | 0.078 | 0.094 |  | 0.014 | 0.014 | 0.014 | 0.014
| ASE |  | 0.083 | 0.081 | 0.081 | 0.083 |  | 0.014 | 0.014 | 0.014 | 0.014
1000 | Mean |  | 1.007 | 1.003 | 1.003 | 1.004 |  | 0.100 | 0.100 | 0.100 | 0.100
| SSE |  | 0.063 | 0.060 | 0.060 | 0.065 |  | 0.011 | 0.011 | 0.011 | 0.011
| ASE |  | 0.064 | 0.063 | 0.063 | 0.064 |  | 0.011 | 0.011 | 0.011 | 0.011

Table S.4: Different estimation approaches applied to simulated count time
series from INGARCH$(1,1)$-CMEM with Bin-counting series: Mean, simulated
standard error (SSE), and mean approximated standard error (ASE) for different
sample sizes $n$.

$n$ |  | Param. | PQ | NQ | EQ | 2W | Param. | PQ | NQ | EQ | 2W
---|---|---|---|---|---|---|---|---|---|---|---
300 | Mean | $a_{0}=3$ | 3.750 | 3.574 | 3.575 | 3.601 | $a_{0}=3$ | 4.131 | 4.256 | 4.267 | 4.169
| SSE |  | 1.412 | 1.234 | 1.235 | 1.369 |  | 1.217 | 1.212 | 1.213 | 1.181
| ASE |  | 1.310 | 1.070 | 1.069 | 1.077 |  | 1.100 | 1.080 | 1.082 | 1.064
600 | Mean |  | 3.376 | 3.272 | 3.275 | 3.279 |  | 3.587 | 3.687 | 3.696 | 3.654
| SSE |  | 0.887 | 0.748 | 0.749 | 0.789 |  | 0.767 | 0.768 | 0.770 | 0.746
| ASE |  | 0.880 | 0.705 | 0.704 | 0.706 |  | 0.733 | 0.715 | 0.716 | 0.710
1000 | Mean |  | 3.234 | 3.165 | 3.167 | 3.165 |  | 3.359 | 3.440 | 3.448 | 3.427
| SSE |  | 0.654 | 0.542 | 0.543 | 0.553 |  | 0.562 | 0.571 | 0.573 | 0.549
| ASE |  | 0.671 | 0.532 | 0.532 | 0.532 |  | 0.554 | 0.538 | 0.539 | 0.536
300 | Mean | $a_{1}=0.35$ | 0.335 | 0.348 | 0.349 | 0.346 | $a_{1}=0.35$ | 0.366 | 0.369 | 0.369 | 0.367
| SSE |  | 0.072 | 0.065 | 0.065 | 0.067 |  | 0.047 | 0.045 | 0.045 | 0.046
| ASE |  | 0.079 | 0.065 | 0.065 | 0.065 |  | 0.048 | 0.045 | 0.045 | 0.045
600 | Mean |  | 0.342 | 0.350 | 0.350 | 0.349 |  | 0.361 | 0.363 | 0.363 | 0.362
| SSE |  | 0.054 | 0.045 | 0.045 | 0.046 |  | 0.033 | 0.031 | 0.031 | 0.031
| ASE |  | 0.057 | 0.046 | 0.046 | 0.046 |  | 0.034 | 0.032 | 0.032 | 0.032
1000 | Mean |  | 0.345 | 0.350 | 0.350 | 0.350 |  | 0.357 | 0.359 | 0.359 | 0.359
| SSE |  | 0.042 | 0.035 | 0.035 | 0.035 |  | 0.026 | 0.024 | 0.024 | 0.024
| ASE |  | 0.045 | 0.035 | 0.035 | 0.035 |  | 0.026 | 0.025 | 0.025 | 0.025
300 | Mean | $b_{1}=0.5$ | 0.467 | 0.467 | 0.466 | 0.467 | $b_{1}=0.5$ | 0.427 | 0.419 | 0.418 | 0.425
| SSE |  | 0.110 | 0.095 | 0.095 | 0.105 |  | 0.077 | 0.074 | 0.074 | 0.075
| ASE |  | 0.106 | 0.089 | 0.089 | 0.089 |  | 0.078 | 0.076 | 0.076 | 0.075
600 | Mean |  | 0.484 | 0.484 | 0.484 | 0.484 |  | 0.460 | 0.453 | 0.452 | 0.455
| SSE |  | 0.073 | 0.061 | 0.061 | 0.064 |  | 0.053 | 0.051 | 0.051 | 0.050
| ASE |  | 0.073 | 0.060 | 0.060 | 0.060 |  | 0.054 | 0.053 | 0.053 | 0.052
1000 | Mean |  | 0.490 | 0.490 | 0.490 | 0.490 |  | 0.475 | 0.469 | 0.469 | 0.470
| SSE |  | 0.055 | 0.046 | 0.046 | 0.047 |  | 0.041 | 0.040 | 0.040 | 0.039
| ASE |  | 0.056 | 0.046 | 0.046 | 0.046 |  | 0.042 | 0.040 | 0.040 | 0.040
300 | Mean | $\sigma^{2}=1$ | 1.020 | 1.011 | 1.010 | 1.022 | $\sigma^{2}=0.1$ | 0.101 | 0.101 | 0.101 | 0.101
| SSE |  | 0.112 | 0.105 | 0.105 | 0.158 |  | 0.018 | 0.018 | 0.018 | 0.018
| ASE |  | 0.107 | 0.103 | 0.103 | 0.109 |  | 0.017 | 0.017 | 0.017 | 0.017
600 | Mean |  | 1.013 | 1.007 | 1.007 | 1.010 |  | 0.101 | 0.101 | 0.101 | 0.101
| SSE |  | 0.081 | 0.075 | 0.075 | 0.087 |  | 0.013 | 0.013 | 0.013 | 0.013
| ASE |  | 0.076 | 0.073 | 0.073 | 0.075 |  | 0.012 | 0.012 | 0.012 | 0.012
1000 | Mean |  | 1.007 | 1.004 | 1.003 | 1.005 |  | 0.100 | 0.100 | 0.100 | 0.100
| SSE |  | 0.058 | 0.057 | 0.056 | 0.060 |  | 0.010 | 0.010 | 0.010 | 0.010
| ASE |  | 0.057 | 0.056 | 0.056 | 0.057 |  | 0.010 | 0.010 | 0.009 | 0.010

Table S.5: Different estimation approaches applied to simulated count time
series from INGARCH$(1,1)$-CMEM with Poi-counting series: Mean, simulated
standard error (SSE), and mean approximated standard error (ASE) for different
sample sizes $n$.

$n$ |  | Param. | PQ | NQ | EQ | 2W | Param. | PQ | NQ | EQ | 2W
---|---|---|---|---|---|---|---|---|---|---|---
300 | Mean | $a_{0}=1$ | 1.452 | 1.435 | 1.446 | 1.474 | $a_{0}=1$ | 1.474 | 1.527 | 1.544 | 1.481
| SSE |  | 0.738 | 0.695 | 0.699 | 0.857 |  | 0.640 | 0.632 | 0.638 | 0.660
| ASE |  | 0.612 | 0.541 | 0.544 | 0.556 |  | 0.549 | 0.535 | 0.541 | 0.523
600 | Mean |  | 1.208 | 1.206 | 1.217 | 1.211 |  | 1.225 | 1.268 | 1.283 | 1.238
| SSE |  | 0.388 | 0.361 | 0.367 | 0.416 |  | 0.355 | 0.356 | 0.364 | 0.350
| ASE |  | 0.382 | 0.333 | 0.335 | 0.334 |  | 0.342 | 0.330 | 0.334 | 0.325
1000 | Mean |  | 1.127 | 1.126 | 1.135 | 1.123 |  | 1.130 | 1.163 | 1.175 | 1.145
| SSE |  | 0.279 | 0.257 | 0.263 | 0.264 |  | 0.256 | 0.255 | 0.261 | 0.249
| ASE |  | 0.283 | 0.244 | 0.246 | 0.244 |  | 0.252 | 0.241 | 0.244 | 0.239
300 | Mean | $a_{1}=0.25$ | 0.246 | 0.257 | 0.258 | 0.255 | $a_{1}=0.25$ | 0.256 | 0.263 | 0.264 | 0.261
| SSE |  | 0.061 | 0.058 | 0.058 | 0.061 |  | 0.051 | 0.049 | 0.049 | 0.050
| ASE |  | 0.067 | 0.059 | 0.059 | 0.059 |  | 0.053 | 0.051 | 0.051 | 0.050
600 | Mean |  | 0.248 | 0.255 | 0.256 | 0.255 |  | 0.254 | 0.260 | 0.261 | 0.258
| SSE |  | 0.044 | 0.040 | 0.040 | 0.041 |  | 0.036 | 0.034 | 0.034 | 0.034
| ASE |  | 0.048 | 0.041 | 0.041 | 0.041 |  | 0.037 | 0.035 | 0.036 | 0.035
1000 | Mean |  | 0.248 | 0.253 | 0.254 | 0.253 |  | 0.252 | 0.256 | 0.257 | 0.255
| SSE |  | 0.035 | 0.031 | 0.031 | 0.031 |  | 0.027 | 0.026 | 0.026 | 0.026
| ASE |  | 0.037 | 0.032 | 0.032 | 0.032 |  | 0.029 | 0.027 | 0.027 | 0.027
300 | Mean | $b_{1}=0.65$ | 0.600 | 0.593 | 0.591 | 0.590 | $b_{1}=0.65$ | 0.594 | 0.582 | 0.580 | 0.589
| SSE |  | 0.114 | 0.104 | 0.105 | 0.123 |  | 0.090 | 0.085 | 0.086 | 0.091
| ASE |  | 0.102 | 0.091 | 0.092 | 0.093 |  | 0.087 | 0.084 | 0.085 | 0.083
600 | Mean |  | 0.628 | 0.622 | 0.620 | 0.622 |  | 0.622 | 0.613 | 0.611 | 0.618
| SSE |  | 0.066 | 0.059 | 0.059 | 0.067 |  | 0.055 | 0.053 | 0.054 | 0.053
| ASE |  | 0.067 | 0.059 | 0.059 | 0.059 |  | 0.057 | 0.055 | 0.055 | 0.054
1000 | Mean |  | 0.637 | 0.633 | 0.631 | 0.633 |  | 0.634 | 0.627 | 0.625 | 0.630
| SSE |  | 0.049 | 0.044 | 0.044 | 0.045 |  | 0.042 | 0.040 | 0.041 | 0.040
| ASE |  | 0.051 | 0.044 | 0.044 | 0.044 |  | 0.043 | 0.041 | 0.042 | 0.041
300 | Mean | $\sigma^{2}=1$ | 1.039 | 1.020 | 1.018 | 1.051 | $\sigma^{2}=0.4$ | 0.414 | 0.408 | 0.407 | 0.415
| SSE |  | 0.155 | 0.125 | 0.123 | 0.258 |  | 0.057 | 0.050 | 0.050 | 0.072
| ASE |  | 0.143 | 0.131 | 0.130 | 0.153 |  | 0.048 | 0.044 | 0.044 | 0.050
600 | Mean |  | 1.026 | 1.015 | 1.013 | 1.025 |  | 0.410 | 0.406 | 0.405 | 0.409
| SSE |  | 0.111 | 0.090 | 0.088 | 0.147 |  | 0.039 | 0.035 | 0.035 | 0.041
| ASE |  | 0.100 | 0.094 | 0.093 | 0.101 |  | 0.034 | 0.032 | 0.031 | 0.034
1000 | Mean |  | 1.017 | 1.011 | 1.010 | 1.014 |  | 0.407 | 0.405 | 0.404 | 0.406
| SSE |  | 0.081 | 0.070 | 0.068 | 0.090 |  | 0.029 | 0.027 | 0.027 | 0.029
| ASE |  | 0.076 | 0.072 | 0.072 | 0.075 |  | 0.026 | 0.024 | 0.024 | 0.025

Table S.6: Different estimation approaches applied to simulated count time
series from INGARCH$(1,1)$-CMEM with Bin-counting series: Mean, simulated
standard error (SSE), and mean approximated standard error (ASE) for different
sample sizes $n$.

$n$ |  | Param. | PQ | NQ | EQ | 2W | Param. | PQ | NQ | EQ | 2W
---|---|---|---|---|---|---|---|---|---|---|---
300 | Mean | $a_{0}=1$ | 1.477 | 1.460 | 1.470 | 1.496 | $a_{0}=1$ | 1.517 | 1.581 | 1.600 | 1.531
| SSE |  | 0.771 | 0.712 | 0.712 | 0.884 |  | 0.642 | 0.637 | 0.644 | 0.674
| ASE |  | 0.619 | 0.540 | 0.541 | 0.554 |  | 0.546 | 0.524 | 0.528 | 0.512
600 | Mean |  | 1.215 | 1.214 | 1.226 | 1.230 |  | 1.247 | 1.300 | 1.316 | 1.274
| SSE |  | 0.396 | 0.361 | 0.367 | 0.447 |  | 0.362 | 0.365 | 0.373 | 0.369
| ASE |  | 0.383 | 0.328 | 0.330 | 0.331 |  | 0.340 | 0.323 | 0.326 | 0.319
1000 | Mean |  | 1.133 | 1.133 | 1.142 | 1.135 |  | 1.142 | 1.180 | 1.193 | 1.169
| SSE |  | 0.283 | 0.255 | 0.260 | 0.274 |  | 0.255 | 0.253 | 0.259 | 0.253
| ASE |  | 0.285 | 0.241 | 0.242 | 0.242 |  | 0.250 | 0.235 | 0.236 | 0.233
300 | Mean | $a_{1}=0.25$ | 0.246 | 0.257 | 0.258 | 0.255 | $a_{1}=0.25$ | 0.257 | 0.264 | 0.265 | 0.262
| SSE |  | 0.060 | 0.055 | 0.056 | 0.059 |  | 0.049 | 0.047 | 0.047 | 0.048
| ASE |  | 0.066 | 0.058 | 0.058 | 0.058 |  | 0.051 | 0.048 | 0.048 | 0.048
600 | Mean |  | 0.248 | 0.255 | 0.256 | 0.255 |  | 0.255 | 0.260 | 0.261 | 0.259
| SSE |  | 0.044 | 0.039 | 0.039 | 0.040 |  | 0.034 | 0.031 | 0.031 | 0.032
| ASE |  | 0.047 | 0.040 | 0.040 | 0.040 |  | 0.036 | 0.033 | 0.034 | 0.033
1000 | Mean |  | 0.248 | 0.253 | 0.254 | 0.254 |  | 0.254 | 0.258 | 0.258 | 0.257
| SSE |  | 0.035 | 0.030 | 0.030 | 0.031 |  | 0.027 | 0.025 | 0.025 | 0.025
| ASE |  | 0.037 | 0.031 | 0.031 | 0.031 |  | 0.028 | 0.026 | 0.026 | 0.026
300 | Mean | $b_{1}=0.65$ | 0.598 | 0.591 | 0.589 | 0.588 | $b_{1}=0.65$ | 0.589 | 0.576 | 0.574 | 0.583
| SSE |  | 0.114 | 0.102 | 0.102 | 0.123 |  | 0.087 | 0.083 | 0.083 | 0.090
| ASE |  | 0.102 | 0.090 | 0.090 | 0.092 |  | 0.084 | 0.081 | 0.081 | 0.079
600 | Mean |  | 0.628 | 0.622 | 0.620 | 0.619 |  | 0.619 | 0.609 | 0.607 | 0.613
| SSE |  | 0.066 | 0.058 | 0.059 | 0.069 |  | 0.054 | 0.052 | 0.052 | 0.054
| ASE |  | 0.066 | 0.058 | 0.058 | 0.058 |  | 0.055 | 0.053 | 0.053 | 0.052
1000 | Mean |  | 0.636 | 0.632 | 0.631 | 0.632 |  | 0.632 | 0.624 | 0.623 | 0.626
| SSE |  | 0.049 | 0.043 | 0.043 | 0.046 |  | 0.040 | 0.039 | 0.039 | 0.039
| ASE |  | 0.050 | 0.043 | 0.043 | 0.043 |  | 0.041 | 0.039 | 0.040 | 0.039
300 | Mean | $\sigma^{2}=1$ | 1.040 | 1.022 | 1.020 | 1.067 | $\sigma^{2}=0.4$ | 0.413 | 0.409 | 0.408 | 0.415
| SSE |  | 0.144 | 0.115 | 0.113 | 0.525 |  | 0.047 | 0.044 | 0.043 | 0.059
| ASE |  | 0.121 | 0.110 | 0.108 | 0.137 |  | 0.035 | 0.032 | 0.032 | 0.036
600 | Mean |  | 1.024 | 1.014 | 1.012 | 1.025 |  | 0.409 | 0.406 | 0.406 | 0.409
| SSE |  | 0.096 | 0.080 | 0.078 | 0.136 |  | 0.034 | 0.031 | 0.030 | 0.037
| ASE |  | 0.083 | 0.077 | 0.076 | 0.083 |  | 0.024 | 0.023 | 0.022 | 0.024
1000 | Mean |  | 1.016 | 1.010 | 1.009 | 1.014 |  | 0.406 | 0.404 | 0.404 | 0.405
| SSE |  | 0.072 | 0.061 | 0.060 | 0.083 |  | 0.026 | 0.024 | 0.023 | 0.026
| ASE |  | 0.063 | 0.060 | 0.059 | 0.062 |  | 0.018 | 0.017 | 0.017 | 0.018

Table S.7: Different estimation approaches applied to simulated count time
series from INGARCH$(1,1)$-CMEM with Poi-counting series, like in Table S.1,
but by assuming a Bin-counting series (model misspecification): Mean,
simulated standard error (SSE), and mean approximated standard error (ASE) for
different sample sizes $n$.

$n$ |  | Param. | PQ | NQ | EQ | 2W | Param. | PQ | NQ | EQ | 2W
---|---|---|---|---|---|---|---|---|---|---|---
300 | Mean | $a_{0}=2.8$ | 2.997 | 2.922 | 2.915 | 2.936 | $a_{0}=2.8$ | 2.965 | 2.936 | 2.932 | 2.937
| SSE |  | 0.722 | 0.686 | 0.688 | 0.696 |  | 0.681 | 0.668 | 0.671 | 0.669
| ASE |  | 0.790 | 0.688 | 0.685 | 0.689 |  | 0.727 | 0.673 | 0.671 | 0.672
600 | Mean |  | 2.904 | 2.858 | 2.854 | 2.865 |  | 2.889 | 2.875 | 2.874 | 2.875
| SSE |  | 0.529 | 0.492 | 0.493 | 0.498 |  | 0.488 | 0.474 | 0.476 | 0.475
| ASE |  | 0.544 | 0.477 | 0.475 | 0.476 |  | 0.503 | 0.468 | 0.467 | 0.467
1000 | Mean |  | 2.871 | 2.839 | 2.837 | 2.841 |  | 2.855 | 2.844 | 2.843 | 2.843
| SSE |  | 0.406 | 0.375 | 0.376 | 0.379 |  | 0.380 | 0.366 | 0.367 | 0.366
| ASE |  | 0.419 | 0.367 | 0.366 | 0.367 |  | 0.387 | 0.360 | 0.359 | 0.359
300 | Mean | $a_{1}=0.4$ | 0.384 | 0.393 | 0.394 | 0.393 | $a_{1}=0.4$ | 0.392 | 0.395 | 0.396 | 0.396
| SSE |  | 0.081 | 0.077 | 0.077 | 0.077 |  | 0.062 | 0.060 | 0.060 | 0.060
| ASE |  | 0.085 | 0.075 | 0.075 | 0.075 |  | 0.064 | 0.059 | 0.059 | 0.059
600 | Mean |  | 0.392 | 0.397 | 0.397 | 0.397 |  | 0.397 | 0.398 | 0.399 | 0.398
| SSE |  | 0.057 | 0.053 | 0.053 | 0.053 |  | 0.044 | 0.042 | 0.043 | 0.042
| ASE |  | 0.061 | 0.053 | 0.053 | 0.053 |  | 0.045 | 0.042 | 0.042 | 0.042
1000 | Mean |  | 0.395 | 0.399 | 0.399 | 0.399 |  | 0.398 | 0.399 | 0.399 | 0.399
| SSE |  | 0.045 | 0.041 | 0.041 | 0.041 |  | 0.034 | 0.033 | 0.033 | 0.033
| ASE |  | 0.048 | 0.042 | 0.041 | 0.041 |  | 0.035 | 0.033 | 0.033 | 0.033
300 | Mean | $b_{1}=0.2$ | 0.180 | 0.184 | 0.185 | 0.182 | $b_{1}=0.2$ | 0.182 | 0.183 | 0.184 | 0.183
| SSE |  | 0.121 | 0.116 | 0.117 | 0.118 |  | 0.112 | 0.110 | 0.111 | 0.110
| ASE |  | 0.130 | 0.116 | 0.115 | 0.116 |  | 0.120 | 0.112 | 0.112 | 0.112
600 | Mean |  | 0.189 | 0.192 | 0.193 | 0.191 |  | 0.190 | 0.190 | 0.190 | 0.190
| SSE |  | 0.089 | 0.083 | 0.084 | 0.084 |  | 0.080 | 0.079 | 0.079 | 0.079
| ASE |  | 0.090 | 0.081 | 0.080 | 0.080 |  | 0.083 | 0.078 | 0.078 | 0.078
1000 | Mean |  | 0.192 | 0.194 | 0.195 | 0.194 |  | 0.194 | 0.194 | 0.194 | 0.194