Parameters Estimation in Non-Negative Integer-Valued Time Series: Approach Based on Probability Generating Functions

Abstract

This manuscript deals with a parameter estimation of a non-negative integer-valued (NNIV) time series based on the so-called probability generating function (PGF) method. The theoretical background of the PGF estimation technique for a very general, stationary class of NNIV time series is described, as well as the asymptotic properties of the obtained estimates. After that, a particular emphasis is given to PGF estimators of independent identical distributed (IID) and integer-valued non-negative autoregressive (INAR) series. A Monte Carlo study of the thus obtained PGF estimates, based on a numerical integration of the appropriate objective function, is also presented. For this purpose, numerical quadrature formulas were computed using Gegenbauer orthogonal polynomials. Finally, the application of the PGF estimators in the dynamic analysis of some actual data is given.

1. Introduction

Non-negative integer-valued (NNIV) time series represent one of the most important stochastic models that have attracted the attention of numerous studies (for some more recent ones, see, e.g., [1,2,3,4,5,6,7]). Due to the specific structure of an NNIV series, the procedure for estimating their parameters may differ from that used for similar non-integer stochastic models. Thus, in the case of the most famous NNIV series, the so-called integer non-negative autoregressive (INAR) processes, various estimation techniques have been developed. For instance, in [8], conditional least square estimators were described, while in [9], estimators based on the conditional maximum likelihood method were discussed. On the other hand, the authors in [10,11] considered, among others, moment-based estimation techniques, i.e., the well-known Yule–Walker (YW) equations. Finally, some generalizations of these methods based on the high-order statistics can be found in [12,13].

However, in order to apply any of the mentioned methods, the estimation functions must usually be given in a finite, closed form and (most often) bounded within a certain parameter space. Unfortunately, sometimes, such conditions cannot be met, both in the case of independently identically distributed (IID) and stochastically dependent time series. Therefore, an estimation technique based on probability generating functions (PGFs) can be a useful alternative approach, since PGFs exist and are bounded for any NNIV time series, at least on the interval $[-1,1]$. Another advantage of PGF estimators is that the estimation procedure is close to the empirical characteristic function (ECF) method introduced in [14,15] and developed later in many studies, cf. [16,17,18]. Generally speaking, the idea behind the PGF method is based on the fact that PGF contains the same amount of information about the distribution of some random variable (RV) as the appropriate characteristic function (CF). Thus, it is expected that the minimization of ‘the distance’ between the PGF of some time series $\left(X_t\right)$ and its empirical PGF gives an efficient parameter estimator. In this way, the PGF method, similar to the ECF method for non-integer time series, could be applied to various kinds of stationary NNIV series. For these reasons, our main goal here is to obtain efficient PGF parameter estimates for the case of a rather general, stationary class of NNIV time series, with applications that will also be described in the following.

Let us point out that, to the best of our knowledge, there are only few studies where the PGF estimation procedure has been discussed so far. First, Esquivel [19] has shown the uniform laws of large numbers and the functional limit theorem for the empirical PGFs. Thereafter, Stojanović et al. [20] used a more general PGF estimation procedure, based on numerical integration techniques, in order to obtain parameter estimators of the so-called noise-indicator integer non-negative autoregressive (NIINAR) time series model. Our main motive here is to generalize these results, i.e., to give theoretical foundations of the PGF estimation procedure in the most general form. According to that, as it will be seen, the PGF estimates can be practically implemented in the estimation of some of the most important and often used NNIV times series, as well as in modelling the dynamics of an actual data set of the number of COVID-19 death cases.

For this purpose, in Section 2, we give some basic facts that are necessary to formulate the PGF estimation method. Additionally, under some regulatory conditions, we investigate the asymptotic properties and efficiency of the PGF estimators for an arbitrary stationary NNIV time series. After that, in Section 3, we analyze some specific PGF estimates for cases of independent identically distributed (IID) and integer-valued autoregressive (INAR) time series. For both types of time series, Monte Carlo simulations of the PGF estimates are computed and their asymptotic properties, such as efficiency and asymptotic normality, are examined too. Section 4 presents an application of the PGF method in parameters and distribution estimation of death cases in the period before and after the appearance of the SARS-CoV2 virus infection on the territory of the Republic of Serbia. Finally, Section 5 contains some concluding remarks.

2. PGF Estimates: Definition and Asymptotic Properties

Let $\left(X_t\right)$, $t\in \mathbb{Z}$ be a non-negative, integer-valued, and stationary time series with finite moments and probability mass distribution (PMF)

$p_X(x;\theta ):=P\left\{X_t=x\right\},x\in {\mathbb{Z}}^{+}=\left\{0,1,2,\cdots \right\},$

which depends on a vector of unknown parameters $\theta ={({\theta }_1,\cdots ,{\theta }_r)}^{\prime }\in {\mathbb{R}}^r$, $r\ge 1$. Additionally, let $\{X_1,X_2,\cdots ,X_T\}$ be some finite realization of the series $\left(X_t\right)$, of the length T. In order to obtain the PGF estimates of $\theta$, we introduce the following terms.

Let $\mathbf{u}={(u_1,\cdots ,u_r)}^{\prime }\in {\mathbb{R}}^r$ and ${\mathbf{X}}_t^{\left(r\right)}:={(X_t,\cdots ,X_{t+r-1})}^{\prime }$, $t\in \mathbb{Z}$ be the overlapping blocks of the process $\left(X_t\right)$, and

(1)${\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta ):=\mathrm{E}\left[u_1^{X_t}\cdots u_r^{X_{t+r-1}}\right]$

be the r-dimensional PGF of the random vector ${\mathbf{X}}_t^{\left(r\right)}$. Due to the stationarity of series $\left(X_t\right)$, it is obvious that the function in Equation (1) does not depend on the choice of overlapping blocks, i.e., moment of time t. The main aim of the PGF method is to minimize ‘the distance’ between the theoretical PGF, given by Equation (1), and the corresponding Empirical PGF (of order $r\ge 1)$:

${\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right):=\frac{1}{T-r+1}\sum _{t=1}^{T-r+1}{u}_1^{X_t}\cdots u_r^{X_{t+r-1}}.$

As the series $\left(X_t\right)$ is stationary, it follows that ${\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)$ is an unbiased estimator of ${\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)$, i.e.,

(2)$\mathrm{E}\left[{\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)\right]={\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0),$

where ${\theta }_0$ is the true (but unknown) value of parameter $\theta$. It is also well-known that theoretical PGF ${\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )$ is well defined at least for all values $\mathbf{u}\in {[-1,1]}^r$ (see, e.g., [20]). Thus, the objective function can be defined as

(3)$S_T^{\left(r\right)}\left(\theta \right):={\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right){\left({\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )-{\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)\right)}^2\mathrm{d}\mathbf{u},$

where $\mathrm{d}\mathbf{u}:=\mathrm{d}{u}_1\cdots \mathrm{d}{u}_r$ and $\omega :{\mathbb{R}}^r\to {\mathbb{R}}^{+}$ is a weight function, integrable on ${[-1,1]}^r$. The PGF estimates will be obtained by minimizing the objective function in Equation (3) with respect to $\theta$, i.e., they are solutions of the minimization equation

(4)${\widehat{\theta }}_T^{\left(r\right)}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\theta \in \mathsf{\Theta }}{min}{S}_T^{\left(r\right)}\left(\theta \right),$

where $\mathsf{\Theta }\subseteq {\mathbf{R}}^r$ is a parameter space. In order to solve this equation, we use numerical methods, i.e., procedures based on numerical integration, which is discussed in Section 3. In the sequel, we will examine strong consistency and asymptotic normality (AN) of the PGF estimators under some regularity conditions.

Theorem 1. 

Let ${\theta }_0$ be the exact value of the parameter θ, and let ${\widehat{\theta }}_T^{\left(r\right)}$, for an arbitrary $T=1,2,\dots$, be the solutions of Equation (4). Additionally, let us suppose that the following regularity conditions are fulfilled:

• (A₁) 

  Θ is a compact set, where ${\theta }_0\in \mathrm{int}\left(\mathsf{\Theta }\right)$, and ${\widehat{\theta }}_T^{\left(r\right)}\in \mathrm{int}\left(\mathsf{\Theta }\right)$ for T large enough.

• (A₂) 

  At the point $\theta ={\theta }_0$ function

  $S_0^{\left(r\right)}\left(\theta \right):={\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right){\left({\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )-{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)\right)}^2\mathrm{d}\mathbf{u}$

  has a unique minimum $S_0^{\left(r\right)}\left({\theta }_0\right)=0$.

• (A₃) 

  ${\displaystyle \frac{{\partial }^2{S}_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta \partial {\theta }^{\prime }}}$ is a regular matrix.

• (A₄) 

  ${\displaystyle \frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)}{\partial \theta }}{\displaystyle \frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)}{\partial {\theta }^{\prime }}}$ is a non-zero matrix, uniformly bounded by the strictly positive ω-integrable function $\mathcal{L}:{\mathbb{R}}^r\to {\mathbb{R}}^{+}$.

• (A₅) 

  The covariance function ${\gamma }_K(h;\theta ):=\mathrm{Cov}\left[K_t^{\left(r\right)}\left(\theta \right),K_{t+h}^{\left(r\right)}\left(\theta \right)\right]$ of the series

  (5)$K_t^{\left(r\right)}\left(\theta \right):={\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right)\left[{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )-u_1^{X_t}\cdots u_r^{X_{t+r-1}}\right]\frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )}{\partial \theta }\mathrm{d}\mathbf{u}.$

  is absolutely summable when $\theta ={\theta }_0$, i.e.,

  $\sum _{h=-\infty }^{\infty }\left|{\gamma }_K(h;{\theta }_0)\right|<+\infty .$

At the same time, the integral in Equation (5) is non-singular and exists for each $\theta \in \mathsf{\Theta }$. Then, ${\widehat{\theta }}_T^{\left(r\right)}$ is a strictly consistent and AN estimator for θ.

Proof. 

Firstly, we check the sufficient consistency conditions of the so-called extremum estimators (for more details, see Newey and McFadden [21]). According to given assumptions, the series $\left(X_t\right)$ is ergodic, so that Equation (2) and the strong law of large numbers (SLLN) give

(6)${\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)\stackrel{\mathrm{as}}{⟶}{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0),T\to \infty ,$

where “as” denotes the almost sure convergence. Furthermore, according to assumption $\left(A_1\right)$, the function ${\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )$ is continuous on the compact ${[-1,1]}^r\times \mathsf{\Theta }$, and ${\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)$ is continuous on the compact $\mathsf{\Theta }$. Hence, there exist constants $C_1,C_2>0$ such that the inequalities

$\begin{array}{cc}\hfill \underset{(\mathbf{u};\theta )\in {[-1,1]}^r\times \mathsf{\Theta }}{max}\left|{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )\right|& \le C_1<+\infty ,\underset{\mathbf{u}\in {[-1,1]}^r}{max}\left|{\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)\right|\le C_2<+\infty \hfill \end{array}$

hold. According to these, one obtains

$\begin{array}{cc}\hfill |S_T^{\left(r\right)}\left(\theta \right)-S_0^{\left(r\right)}\left(\theta \right)|& =|{\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right){\left({\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )-{\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)\right)}^2\mathrm{d}\mathbf{u}\hfill \\ \hfill & -{\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right){\left({\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )-{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)\right)}^2\mathrm{d}\mathbf{u}|\hfill \\ \hfill & =|{\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right)\left[{\left({\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)\right)}^2-{\left({\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)\right)}^2\right]\mathrm{d}\mathbf{u}\hfill \\ \hfill & +2{\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right){\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )\left({\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)-{\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)\right)\mathrm{d}\mathbf{u}|\hfill \\ \hfill & \le (C_1+C_2){\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right)|{\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)-{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)|\mathrm{d}\mathbf{u}\hfill \\ \hfill & +2C_1{\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right)|{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)-{\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)|\mathrm{d}\mathbf{u}\hfill \\ \hfill & \le (3C_1+C_2){\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right)\left|{\tilde{\mathsf{\Psi }}}_T^{\left(r\right)}\left(\mathbf{u}\right)-{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)\right|\mathrm{d}\mathbf{u},\hfill \end{array}$

so the last inequality and Equation (6) imply

$\underset{\theta \in \mathsf{\Theta }}{sup}\left|S_T^{\left(r\right)}\left(\theta \right)-S_0^{\left(r\right)}\left({\theta }_0\right)\right|\stackrel{\mathrm{as}}{⟶}0,T\to +\infty .$

Thus, $S_T^{\left(r\right)}\left(\theta \right)$ uniformly converges almost surely toward $S_0^{\left(r\right)}\left(\theta \right)$. According to this, as well as the assumption $\left(A_2\right)$ and Theorem 1 in Newey and McFadden [21], it follows

${\widehat{\theta }}_T^{\left(r\right)}-{\theta }_0\stackrel{\mathrm{as}}{⟶}0,T\to +\infty ,$

i.e., ${\widehat{\theta }}_T^{\left(r\right)}$ is strictly consistent estimator for $\theta$.

Furthermore, let us notice that function $S_T^{\left(r\right)}\left(\theta \right)$ has continuous partial derivatives of the first and second orders, with respect to any component of the vector $\theta$. Hence, the Taylor expansion of $\partial S_T^{\left(r\right)}\left(\theta \right)/\partial \theta$ at $\theta ={\theta }_0$ gives

$\frac{\partial S_T^{\left(r\right)}\left(\theta \right)}{\partial \theta }=\frac{\partial S_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta }+\frac{{\partial }^2{S}_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta \partial {\theta }^{\prime }}·(\theta -{\theta }_0)+o(\theta -{\theta }_0).$

Substituting $\theta$ with ${\widehat{\theta }}_T^{\left(r\right)}$, for a T large enough, under assumption $\left(A_3\right)$ and the fact that $\partial S_T^{\left(r\right)}\left({\widehat{\theta }}_T^{\left(r\right)}\right)/\partial \theta =0$, we have

(7)${\widehat{\theta }}_T^{\left(r\right)}-{\theta }_0=-{\left[\frac{{\partial }^2{S}_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta \partial {\theta }^{\prime }}\right]}^{-1}\frac{\partial S_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta }+o\left({\widehat{\theta }}_T^{\left(r\right)}-{\theta }_0\right).$

Furthermore, the function $S_T^{\left(r\right)}\left(\theta \right)$ can be differentiated under the integral sign, i.e.,

(8)$\begin{array}{ccc}\hfill \frac{\partial S_T^{\left(r\right)}\left(\theta \right)}{\partial \theta }& =& 2{\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right)\left[{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )-{\tilde{\mathsf{\Psi }}}_T\left(\mathbf{u}\right)\right]\frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )}{\partial \theta }\mathrm{d}\mathbf{u},\hfill \\ \hfill \frac{{\partial }^2{S}_T^{\left(r\right)}\left(\theta \right)}{\partial \theta \partial {\theta }^{\prime }}& =& 2{\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right)\left\{\frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )}{\partial \theta }\frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )}{\partial {\theta }^{\prime }}\right.\hfill \end{array}$

(9)$\begin{array}{ccc}& & \left.+\left[{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )-{\tilde{\mathsf{\Psi }}}_T\left(\mathbf{u}\right)\right]\frac{{\partial }^2{\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta )}{\partial \theta \partial {\theta }^{\prime }}\right\}\mathrm{d}\mathbf{u}.\hfill \end{array}$

Now, Equations (8) and (9) give

(10)$E\left[\frac{\partial S_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta }\right]={\mathbf{0}}_{r\times 1},E\left[\frac{{\partial }^2{S}_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta \partial {\theta }^{\prime }}\right]=2\mathbf{V},$

where

$\mathbf{V}={\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right)\frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)}{\partial \theta }\frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)}{\partial {\theta }^{\prime }}\mathrm{d}\mathbf{u}.$

According to assumption $\left(A_4\right)$, there exists an $\omega$-integrable function $\mathcal{L}:{\mathbb{R}}^r\to {\mathbb{R}}^{+}$ such that

$0<∥\frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)}{\partial \theta }\frac{\partial {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};{\theta }_0)}{\partial {\theta }^{\prime }}∥\le \mathcal{L}\left(\mathbf{u}\right),\mathbf{u}\in {[-1,1]}^r,$

where $∥·∥$ is a matrix norm on ${\mathbb{R}}^r\times {\mathbb{R}}^r$. Hence, the inequalities

$0<∥\mathbf{V}∥\le {\int }_{-1}^1\cdots {\int }_{-1}^1\omega \left(\mathbf{u}\right)\mathcal{L}\left(\mathbf{u}\right)\mathrm{d}\mathbf{u}<+\infty$

hold, so Equation (10) and the SLLN give

(11)$\left(\frac{\partial S_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta },\frac{{\partial }^2{S}_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta \partial {\theta }^{\prime }}\right)\stackrel{as}{⟶}(0,2\mathbf{V}),T\to +\infty .$

To prove the AN property, note that using Equations (5) and (8), we can write the gradient of $S_T^{\left(r\right)}\left(\theta \right)$ as

(12)$\frac{\partial S_T^{\left(r\right)}\left(\theta \right)}{\partial \theta }=\frac{2}{T-r+1}\sum _{t=1}^{T-r+1}{K}_t^{\left(r\right)}\left(\theta \right).$

Therefore, according to assumption $\left(A_5\right)$, equality $\mathrm{E}\left[K_t\left({\theta }_0\right)\right]=0$ holds and it follows

$\begin{array}{ccc}\hfill {\mathbf{W}}^2& :=& \underset{T\to \infty }{lim}Var\left[\frac{1}{\sqrt{T-r+1}}\sum _{t=1}^{T-r+1}{K}_t^{\left(r\right)}\left({\theta }_0\right)\right]=\underset{T\to \infty }{lim}\frac{1}{T-r+1}\mathrm{E}{\left[\sum _{t=1}^{T-r+1}{K}_t^{\left(r\right)}\left({\theta }_0\right)\right]}^2\hfill \\ & =& \underset{T\to \infty }{lim}\frac{1}{T-r+1}\sum _{t=1}^{T-r+1}\sum _{s=1}^{T-r+1}E\left[K_t^{\left(r\right)}\left({\theta }_0\right)K_s^{\left(r\right)}\left({\theta }_0\right)\right]\hfill \\ & =& \underset{T\to \infty }{lim}\frac{1}{T-r+1}\sum _{h=-T+r-1}^{T-r+1}\left(T-r+1-h\right){\gamma }_K(h;{\theta }_0)\hfill \\ & =& \underset{T\to \infty }{lim}\sum _{h=-T}^T\left(1-\frac{h}{T}\right){\gamma }_K(h;{\theta }_0)=\sum _{h=-\infty }^{\infty }{\gamma }_K(h;{\theta }_0)<+\infty .\hfill \end{array}$

Note that the last equality follows from Kronecker’s lemma (see, c.f. Fuller [22]). Thus, the proven convergence and Equation (12) imply

(13)$\sqrt{T-r+1}\frac{\partial S_T^{\left(r\right)}\left({\theta }_0\right)}{\partial \theta }\stackrel{d}{⟶}\mathcal{N}({\mathbf{0}}_{r\times 1},4{\mathbf{W}}^2),T\to +\infty ,$

where “d” denotes the convergence in distribution. Now, applying the central limit theorem for stationary processes, as well as Equations (7), (11) and (13), one obtains

$\sqrt{T-r+1}\left({\widehat{\theta }}_T^{\left(r\right)}-{\theta }_0\right)\stackrel{d}{⟶}\mathcal{N}({\mathbf{0}}_{r\times 1},{\mathbf{V}}^{-1}{\mathbf{W}}^2{\mathbf{V}}^{-1}),T\to +\infty ,$

and the statement of theorem is completely proven. □

At the end of this section, we point out that PGF estimators of the r-dimensional parameter $\theta$ are usually obtained using overlapping blocks of the same dimensions. This is also consistent with some of the basic assumptions of the ECF method (for more details, see [14]). For these reasons, below, we examine in detail the PGF estimation procedures for the parameter order (as well as the corresponding PGFs order) $r=1$ and $r=2$.

3. PGF Estimations of IID and INAR Time Series

In this section, we discuss the construction of PGF estimates for two important classes of stationary NNIV time series. First, using a very general form of the so-called power series (PS) distributions, the estimation procedure of independent identical distributed (IID) time series was examined. After that, using the obtained results, PGF estimators of integer-valued autoregressive (INAR) time series are considered.

3.1. Estimation of IID Time Series

The simplest case of the stationary NNIV time series are the IID series, which we denote as $\left({\epsilon }_t\right)$, $t\in \mathbb{Z}$. Similar to that in Stojanović et al. [20], and Bourguignon and Vasconcellos [23], we say that the series $\left({\epsilon }_t\right)$ has a PS distribution if its PMF is given as follows:

(14)$p_{\epsilon }(x;\theta ):=P\{{\epsilon }_t=x\}=\frac{a\left(x\right){\theta }^x}{f\left(\theta \right)},x\in \mathcal{S}.$

Here, $\mathcal{S}\subseteq {\mathbb{Z}}^{+}=\left\{0,1,2,\cdots \right\}$ is the so-called discrete support of the RV ${\epsilon }_t$ and

• $\left(i\right)$

  $a\left(x\right)\ge 0$ is a function that depends (only) on $x\in \mathcal{S}$;

• $\left(ii\right)$

  $\theta >0$ is a (one-dimensional) unknown parameter;

• $\left(iii\right)$

  $f\left(\theta \right):={\sum }_{x\in \mathcal{S}}a\left(x\right){\theta }^x$ is a function on $\theta$, so that $0<f\left(\theta \right)<+\infty$ when $\theta \in (0,R)$, $R>0$.

The expression given by Equation (14) represents the so-called PS distributions that, for some special choices of $a\left(x\right)$, $f\left(\theta \right)$ and $\theta$, allow for obtaining the most of the known distributions (see Table 1, below). In addition, according to assumption $\left(iii\right)$, it is obvious that, in fact, the power series $f\left(\theta \right)$ converges on the interval $(-R,R)$. However, a common assumption is $\theta >0$, so we say that power series $f\left(\theta \right)$ converges on $(0,R)$. Accordingly, the function $f\left(\theta \right)$ is positive and increasing at this interval. Moreover, for an arbitrary $u\in \left[-1,1\right]$, we denote by ${\mathsf{\Psi }}_{\epsilon }^{\left(1\right)}(u;\theta )={\mathsf{\Psi }}_{\epsilon }(u;\theta )$ the PGF (of the first order) of RVs $\left({\epsilon }_t\right)$. Thus, according to Equation (14), one obtains

(15)${\mathsf{\Psi }}_{\epsilon }(u;\theta ):=\mathrm{E}\left[u^{{\epsilon }_t}\right]=\frac{1}{f\left(\theta \right)}\sum _{x\in \mathcal{S}}a\left(x\right){\left(\theta u\right)}^x=\frac{f\left(\theta u\right)}{f\left(\theta \right)},$

where the sum above converges when $u<R/\theta$. Note that Equation (15) allows for a simple computation of first-order PGFs for some of the well-known NNIV distributions, which are also given in Table 1.

Furthermore, let ${\tilde{\mathsf{\Psi }}}_T\left(u\right):=\frac{1}{T}{\sum }_{t=1}^T{u}^{{\epsilon }_t}$ be the appropriate first-order empirical PGF of series $\left({\epsilon }_t\right)$, obtained by its realization ${\epsilon }_1,\cdots ,{\epsilon }_T$. Thus, the objective function can be written as

(16)$\begin{array}{cc}\hfill S_T\left(\theta \right)& :={\int }_{-1}^1\omega \left(u\right){\left({\mathsf{\Psi }}_{\epsilon }(u;\theta )-{\tilde{\mathsf{\Psi }}}_T\left(u\right)\right)}^2\mathrm{d}u=\frac{1}{T^2}{\int }_{-1}^1\omega \left(u\right){\left[\sum _{t=1}^T\left(\frac{f\left(\theta u\right)}{f\left(\theta \right)}-u^{{\epsilon }_t}\right)\right]}^2\mathrm{d}u\hfill \\ \hfill & =\frac{1}{T^2}\sum _{t=1}^T\sum _{s=1}^T{\int }_{-1}^1\omega \left(u\right)\left(\frac{f\left(\theta u\right)}{f\left(\theta \right)}-u^{{\epsilon }_t}\right)\left(\frac{f\left(\theta u\right)}{f\left(\theta \right)}-u^{{\epsilon }_s}\right)\mathrm{d}u,\hfill \end{array}$

where $\omega :\left[-1,1\right]\to {\mathbb{R}}^{+}$ is the weight function. Finally, in this case, the minimization Equation (5) reads as follows

(17)${\widehat{\theta }}_T=arg\underset{0<\theta <R}{min}{S}_T\left(\theta \right).$

As an illustration, Figure 1 shows plots of the PGFs of some typical NNIV IID time series, along with their empirical PGFs.

In the following, we assume that $R<+\infty$, so the set $\overline{(0,R)}=[0,R]$ will be a compact. Thus, assuming that other regularity conditions in Theorem 1 are fulfilled, we can apply the PGF method. For this purpose, we have taken two different types of distribution of the series $\left({\epsilon }_t\right)$. Firstly, we assume that $\left({\epsilon }_t\right)$ is the IID series with Poisson distribution and, later, that it has a geometric distribution. In the both cases, for estimating the parameter $\theta >0$, we generated $N=500$ independent Monte Carlo replications of the series $\left({\epsilon }_t\right)$, of the length $T=1000$. Note that this is close to the length of actual series of COVID-19 death cases, which will be considered in the following. To examine the convergence, i.e., the quality of the PGF estimators, the appropriate estimation errors will also be investigated.

The PGF estimates are computed according to a minimization of the integral in Equation (16). Note that it can be numerically approximated by using an N-point quadrature formula:

$I(\phi ;\omega ):={\int }_{-1}^1{\omega }_k\left(u\right)\phi \left(u\right)du\approx Q_N\left(h\right):=\sum _{j=1}^N{\omega }_k^{\left(j\right)}\phi \left(u_j\right),$

where $\left(u_j\right)$ are the quadrature nodes and ${\omega }_k^{\left(j\right)}$, $j=1,\cdots ,N$ are the corresponding weight coefficients. In our simulation study, quadrature formulas based on Gegenbauer orthogonal polynomials, with weights ${\omega }_k\left(u\right)={(1-u^2)}^{(k-1)/2}$ and $N=20$ nodes, were used. More specifically, we consider Gaussian quadratures based on three well-known special cases of Gegenbauer polynomials, i.e., the appropriate weights:

• $\left(i\right)$

  Chebyshev polynomials (of the first kind): ${\omega }_0\left(u\right)={(1-u^2)}^{-1/2}$;

• $\left(ii\right)$

  Legendre polynomials: ${\omega }_1\left(u\right)\equiv 1$;

• $\left(iii\right)$

  Chebyshev polynomials (of the second kind): ${\omega }_2\left(u\right)={(1-u^2)}^{1/2}$.

The construction of all these quadrature formulas was realized using the Wolfram Mathematica package “OrthogonalPolinomials”, authored by Cvetković and Milovanović [24]. After that, the objective function in Equation (17) is minimized by the “R” procedure “optimize”, based on the Brent minimization algorithm [25] and realized by the authors’ original codes in “R”.

Summary statistics of the PGF estimates ${\widehat{\theta }}_k$, $k=0,1,2$, obtained with the appropriate weights ${\omega }_k\left(u\right)$ and computed via aforementioned estimation procedure, are presented in Table 2. More precisely, it contains the mean values (Mean), minima (Min.), maxima (Max.), and the mean squared estimating errors (MSEE) of obtained PGF estimates of the IID series with Poisson and geometric distribution. (These two distributions are taken as an illustration of the application of the PGF method, that is, to describe the dynamics of actual data in Section 4, below). The values of objective functions ${\widehat{S}}_{T_k}=S_T\left({\widehat{\theta }}_k\right)$, $k=0,1,2$, as the reference estimation errors, are also shown. According to presented results, it is obvious that PGF estimates of both IID series are generally the efficient parameter estimators and have very similar properties. Slightly lower estimation errors are observed in the case of the Gauss–Chebyshev quadrature of the second type, i.e., the weight function ${\omega }_2\left(u\right)$. This is expected, since it is known that this weight forces the points around the coordinate origin and ignores those near the ends of the interval $[-1,1]$. On the other hand, the somewhat larger estimation error in the case of Gauss–Chebyshev quadrature of the first kind is a consequence of ‘forcing’ the ends of the interval $[-1,1]$, where the weight is infinite. In addition, the results of AN testing of the obtained PGF estimates are also presented in Table 2. For this purpose, the Anderson–Darling normality test was used and its test statistic (labelled as AD), along with the corresponding p-values, were calculated using the procedure from the R-package “nortest”, authored by Gross [26]. According to the values thus obtained, it is evident that the AN property is valid for all PGF estimators, which is in accordance with the previous theoretical results, i.e., the statement of Theorem 1.

3.2. Estimation of INAR Time Series

In this part of our work, using the previous assumptions, we define the integer-valued autoregressive process (of the first order) or, simply, INAR(1) process $\left(X_t\right)$, $t\in \mathbb{Z}$ by the recurrence relation

(18)$X_t=\alpha \circ X_{t-1}+{\epsilon }_t.$

Here, $\left({\epsilon }_t\right)$ is the IID series of RVs commonly called innovations, $\alpha \in (0,1)$ is an unknown parameter, and

(19)$\alpha \circ X:=\sum _{j=1}^X{B}_j\left(\alpha \right)$

is the so-called binomial thinning operator (for more details, see [27]). More precisely, X is an arbitrary NNIV RV such that, for any $t\in \mathbb{Z}$, the RVs $B_j=B_j\left(\alpha \right)$ are mutually independent (also independent of X) with Bernoulli’s distribution

$P\{B_j=1\}=1-P\{B_j=0\}=\alpha .$

In that way, it is simply shown that the INAR(1) process $\left(X_t\right)$ is stationary for each $\alpha \in (0,1)$, with the autocorrelation function (ACF)

(20)${\rho }_X\left(k\right):=\mathrm{Corr}[X_t,X_{t+k}]={\alpha }^k,k=0,1,2,\cdots$

As an illustration, Figure 2 presents sample frequency distributions of the INAR(1) process $\left(X_t\right)$, along with their corresponding innovations $\left({\epsilon }_t\right)$, for Poisson and geometric distributions. Furthermore, it will be of interest to use the first-order PGF of the RV $\alpha \circ X$. For that purpose, we prove the following simple statement:

Lemma 1. 

For an arbitrary $\alpha \in (0,1)$ and NNIV RV X, the first-order PGF of the RV $\alpha \circ X$ reads as follows:

(21) ${\mathsf{\Psi }}_{\alpha \circ X}(u;\alpha )={\mathsf{\Psi }}_X\left(1+\alpha (u-1)\right).$

Proof. 

According to the definition of the first-order PGF, Equation (19) and the law of total expectation (see, e.g., Theorem 34.4 in [28]), we have

(22)$\begin{array}{cc}\hfill {\mathsf{\Psi }}_{\alpha \circ X}(u;\alpha )& :=E\left[u^{\alpha \circ X}\right]=E\left[u^{{\sum }_{j=1}^X{B}_j\left(\alpha \right)}\right]=E\left[E\left(u^{{\sum }_{j=1}^X{B}_j\left(\alpha \right)}|X\right)\right]\hfill \\ \hfill & =E\left[\prod _{j=1}^X{\mathsf{\Psi }}_j(u;\alpha )\right]=E\left[{\left({\mathsf{\Psi }}_j(u;\alpha )\right)}^X\right].\hfill \end{array}$

Here, ${\mathsf{\Psi }}_j(u;\alpha ):=E\left[u^{B_j\left(\alpha \right)}\right]=1+\alpha (u-1)$ is the first-order PGF of the RVs $B_j\left(\alpha \right)$. Hence, substituting this expression in Equation (22) obviously gives Equation (21). □

Let us also note that some of the important properties of the INAR(1) process can be found in many studies, cf. [29,30,31,32,33,34,35,36]. This kind of stochastic model was first introduced in the pioneer work of Al-Osh and Alzaid [37]. Here, among others, the first order PGF of the INAR(1) process was derived. In the case of PS-distributed innovations $\left({\epsilon }_t\right)$, we generalized this result in the following way:

Theorem 2. 

Let us assume that the innovations $\left({\epsilon }_t\right)$ of the INAR(1) process have the PS distribution given by Equation (14). Additionally, suppose that the function $g\left(\theta \right)=logf\left(\theta \right)$ has a finite derivative, uniformly bounded on a finite interval $\theta \in (0,R)$. Then, an arbitrary $\alpha \in (0,1)$ INAR (1) series $\left(X_t\right)$, given by Equation (19), has a PGF of the first order:

(23) $\begin{array}{c}\hfill {\mathsf{\Psi }}_X(u;\theta ,\alpha )=\prod _{k=0}^{\infty }\left(\frac{f\left(\left(1+{\alpha }^k(u-1)\right)\theta \right)}{f\left(\theta \right)}\right),\end{array}$

as well as PGFs of order $r\ge 2$:

(24) $\begin{array}{cc}\hfill {\mathsf{\Psi }}_X^{\left(r\right)}(\mathbf{u};\theta ,\alpha )& ={\mathsf{\Psi }}_X\left(\prod _{k=0}^{r-1}\left(1+{\alpha }^k(u_{k+1}-1)\right)\right)\prod _{\ell =2}^r{\mathsf{\Psi }}_{\epsilon }\left(\prod _{k=0}^{r-\ell }\left(1+{\alpha }^k(u_{k+\ell }-1)\right)\right).\hfill \end{array}$

Proof. 

According to a definition of the PS-distributed RVs $\left({\epsilon }_t\right)$, given by Equation (14), as well as their first-order PGF, given by Equation (15), the mathematical expectation of the series $\left({\epsilon }_t\right)$ can be obtained as follows:

${\mu }_{\epsilon }\left(\theta \right):=\mathrm{E}\left[{\epsilon }_t\right]=\frac{\partial {\mathsf{\Psi }}_{\epsilon }(u;\theta )}{\partial u}{|}_{u=1}=\theta \frac{f{}^{\prime }\left(\theta \right)}{f\left(\theta \right)}=\theta g^{\prime }\left(\theta \right).$

Based on this and the assumptions of the theorem, for some constant $M>0$ and any $\theta \in (0,R)$, we have

$\sum _{k=1}^{\infty }P\{{\epsilon }_t\ge k\}=\sum _{k=1}^{\infty }kP\left\{{\epsilon }_t=k\right\}={\mu }_{\epsilon }\left(\theta \right)\le M<+\infty ,$

where the sum above converges uniformly on $\theta \in (0,R)$. On the other hand, the sequence $1/k$, $k=1,2,\cdots$ is monotone and bounded, so Abel’s convergence criterion implies that, when uniformly on $\theta$, the following is valid:

(25)$\sum _{k=1}^{\infty }\frac{1}{k}P\{{\epsilon }_t\ge k\}<+\infty .$

As it is known from Alzaid and Al-Osh [38], the inequality (25) represents a necessary and sufficient condition for the infinite-order integer moving average (INMA) representation of the series $\left(X_t\right)$:

(26)$X_t\stackrel{d}{=}\sum _{k=0}^{\infty }{\alpha }^k\circ {\epsilon }_{t-k},$

where the sum above converges almost surely. According to Equation (26) and the previous Lemma, we have the following:

(27)${\mathsf{\Psi }}_X(u;\theta ,\alpha )=\prod _{k=0}^{\infty }{\mathsf{\Psi }}_{\epsilon }\left(1+{\alpha }^k(u-1);\theta \right),$

and the product above converges absolutely at least for all $u\in [-1,1]$. Hence, substituting Equation (15) in Equation (27), we have Equation (23). In order to prove the second part of the theorem, notice that, according to the definition (18) of the series $\left(X_t\right)$, for an arbitrary $k\in \mathbb{N}$, the following holds:

(28)$X_{t+k}\stackrel{d}{=}{\alpha }^k\circ X_t+\sum _{j=0}^{k-1}{\alpha }^j\circ {\epsilon }_{t+j}.$

Substituting $k=1,\cdots ,r-1$ into Equation (28) and after some computation similar to the previous one, an explicit expression of the PGF in Equation (24) can be easily obtained. □

Remark 1. 

Note that the distribution of the INAR(1) process depends on two unknown parameters $\theta ,\alpha$. Therefore, in order to estimate them, it is of interest to determine the two-dimensional PGF of $\left(X_t\right)$, which is simply obtained by replacing $r=2$ in Equation (24):

(29) ${\mathsf{\Psi }}_X^{\left(2\right)}(u_1,u_2;\theta ,\alpha )={\mathsf{\Psi }}_X\left(u_1(1+\alpha (u_2-1));\theta ,\alpha \right){\mathsf{\Psi }}_{\epsilon }\left(u_2;\theta \right).$

In the following, similar to the case of IID series, the procedure for estimating the parameters of the INAR(1) model with Poisson and geometric innovations $\left({\epsilon }_t\right)$ will be described. Using Poisson distribution for series $\left({\epsilon }_t\right)$, as well as Equations (23) and (29), the one- and two-dimensional PGFs of the INAR(1) process $\left(X_t\right)$, after some simple computations, can be obtained as follows:

$\begin{array}{cc}\hfill {\mathsf{\Psi }}_X(u;\theta ,\alpha )& =\prod _{k=0}^{\infty }\frac{exp\left[\theta \left(1+{\alpha }^k(u-1)\right)\right]}{exp\left(\theta \right)}=exp\left[\frac{\theta (u-1)}{1-\alpha }\right],\hfill \\ \hfill {\mathsf{\Psi }}_X^{\left(2\right)}(u_1,u_2;\theta ,\alpha )& =exp\left[\theta \left(u_1+u_2-1+\frac{\alpha u_1{u}_2-1}{1-\alpha }\right)\right].\hfill \end{array}$

In a similar way, for the appropriate PGFs of the INAR(1) process with geometric innovations, one obtains

$\begin{array}{cc}\hfill {\mathsf{\Psi }}_X(u;\theta ,\alpha )& =\prod _{k=0}^{\infty }\frac{1-\theta }{1-\theta \left(1+{\alpha }^k(u-1)\right)}=exp\left[-\sum _{k=0}^{\infty }ln\left(1-\frac{\theta (u-1)}{1-\theta }{\alpha }^k\right)\right]\hfill \\ \hfill & =exp\left[\sum _{j=1}^{\infty }\frac{{(-1)}^{j-1}}{j(1-{\alpha }^j)}{\left(\frac{\theta (u-1)}{1-\theta }\right)}^j\right],\hfill \\ \hfill {\mathsf{\Psi }}_X^{\left(2\right)}(u_1,u_2;\theta ,\alpha )& =exp\left[\sum _{j=1}^{\infty }\frac{{(-1)}^{j-1}}{j(1-{\alpha }^j)}{\left(\frac{\theta (u_1(1+\alpha (u_2-1))-1)}{1-\theta }\right)}^j\right]·\frac{1-\theta }{1-\theta u_2}.\hfill \end{array}$

Three-dimensional plots of the theoretical and empirical PGFs of INAR(1) processes with Poisson and geometric innovations, respectively, are shown in Figure 3.

Therefore, notice that the PGFs of the INAR(1) process with geometric innovations are not in closed form. However, they depend on alternating sums and can be easily approximated by finite k-term sums. Thus, for instance, in the case of a one-dimensional PGF, the approximation error amounts to

$R_{k+1}(u;\theta ,\alpha )\le \frac{{\theta }^{k+1}{(u-1)}^{k+1}}{(k+1)(1-{\alpha }^{k+1}){(1-\theta )}^{k+1}}{|}_{u=\frac{1}{\theta }}\le \frac{1}{(k+1)(1-{\alpha }^{k+1})}.$

Similar to that in the previous estimation procedure, the PGF estimates $(\widehat{\theta },\widehat{\alpha })$ of unknown parameters $(\theta ,\alpha )$ are computed according to a minimization of the double integral

(30)$S_T^{\left(2\right)}(\theta ,\alpha )=\underset{{[-1,1]}^2}{\int \int }\omega (u_1,u_2){\left({\mathsf{\Psi }}_X^{\left(2\right)}(u_1,u_2;\theta ,\alpha )-{\tilde{\mathsf{\Psi }}}_T(u_1,u_2)\right)}^2\mathrm{d}{u}_1\mathrm{d}{u}_2,$

with respect to some weight function $\omega :{[-1,1]}^2\to {\mathbb{R}}^{+}$. Similar to earlier, the integral in (30) can be numerically approximated by using some N-point cubature formula:

$I(\phi ;\omega ):=\underset{{[-1,1]}^2}{\int \int }\omega (u_1,u_2)\phi (u_1,u_2)d{u}_{1}d{u}_2\approx \sum _{j=1}^N{\omega }_j\phi (u_{1j},u_{2j}),$

where $(u_{1j},u_{2j})$ are the cubature nodes, and ${\omega }_j$ are the corresponding weight coefficients. In our simulations study, we used the same kind of (two-dimensional) weights as in the IID case, i.e., the Gaussian qubature formulas based on Gegenbauer orthogonal polynomials, with weights

${\omega }_k(u_1,u_2)={\left((1-u_1^2)(1-u_2^2)\right)}^{(k-1)/2}$

and $N=36$ nodes. After that, the objective function (30) is minimized by “R” procedure for linearly constrained minimization “constrOptim”, based on the Nelder–Mead optimization method [39].

As the initial values of this estimation procedure, we used the Yule–Walker (YW) estimates, obtained by solving the equations

${\mu }_X={\overline{X}}_T,{\rho }_X\left(1\right)=\widehat{\rho }\left(1\right).$

Here,

$\begin{array}{c}\hfill {\mu }_X:=E\left[X_t\right]=\frac{{\mu }_{\epsilon }}{1-\alpha },{\rho }_X\left(1\right):=\mathrm{Corr}[X_t,X_{t+1}]=\alpha \end{array}$

are, respectively, the mean and first autocorrelation of the INAR(1) series $\left(X_t\right)$, while

${\overline{X}}_T:=\frac{1}{T}\sum _{t=1}^T{X}_t,\widehat{\rho }\left(1\right):=\frac{T^{-1}{\sum }_{t=1}^T{X}_t{X}_{t-1}-{\overline{X}}_T^2}{T^{-1}{\sum }_{t=1}^T{X}_t^2-{\overline{X}}_T^2}$

are their sample estimators. Thus, in the case of Poisson distributed innovations $\left({\epsilon }_t\right)$, the YW estimators are obtained as follows:

$\tilde{\alpha }=\widehat{\rho }\left(1\right),\tilde{\theta }=(1-\tilde{\alpha }){\overline{X}}_T.$

Similarly, for geometric distributed innovations $\left({\epsilon }_t\right)$, the YW estimates are

$\tilde{\alpha }=\widehat{\rho }\left(1\right),\tilde{\theta }=\frac{(1-\tilde{\alpha }){\overline{X}}_T}{1+(1-\tilde{\alpha }){\overline{X}}_T}.$

Summary statistics of the PGF estimates $({\widehat{\theta }}_k$, ${\widehat{\alpha }}_k)$, obtained using the aforementioned numerical estimation procedure with the appropriate weights ${\omega }_k(u_1,u_2)$, $k=0,1,2$, are presented in the following Table 3. Similar to in the IID case, two different types of innovations $\left({\epsilon }_t\right)$, with Poisson and geometric distributions, were used. Thus, in both cases, the unknown parameters $(\theta ,\alpha )$ were estimated through the set of $N=500$ independent Monte Carlo simulations of the innovation series $\left({\epsilon }_t\right)$ and INAR(1) process $\left(X_t\right)$ of the length $T=1000$.

In the same way as before, Table 3 contains the mean values (Mean), minima (Min.), maxima (Max.), as well as the mean squared estimating errors (MSEE) and the objective function $S_T^{\left(2\right)}$ for all obtained PGF estimates. Notice that the thus obtained PGF estimators have very similar properties as in the one-dimensional case, with somewhat higher estimation errors, which is expected. Namely, they are computed using a two-step procedure, using the YW estimated values as the initial ones.

Similar to earlier, the summarized values of $AD$ test statistics, along with the appropriate p-values, are also shown in Table 3. It can be seen that the AN property is confirmed in almost all cases, that is, for most of the PGF estimators of the parameters $(\theta ,\alpha )$, as well as the corresponding weight functions ${\omega }_k(u_1,u_2)$. We assume, in the same way as it is shown in Stojanović et al. [40], that the variations that occur in some simulations, at the significant level of $0.01<p<0.05$, are also reasonable as a consequence of this specific, two-stage estimation procedure.

4. Application of the PGF Estimates

In this section, we indicate some possibilities of the practical application of PGF estimators in stochastic modeling of real data. Additionally, our motivation is also to point out the different dynamics of time series of death cases in the time period before and after the emergence of the SARS-CoV2 virus. We emphasize that much attention has been paid to the COVID-19 disease, from the beginning of the pandemic until today. For these reasons, various theoretical models were proposed for research of this still current issue. Rigorous mathematical ones, usually based on solving coupled partial equations, have been proposed, e.g., in [41,42,43]. Some other authors [44,45,46] used both deterministic and stochastic approaches to predict the various effects caused by the COVID-19 pandemic. A particular interesting approach is given in [47,48] where the NNIV time series were used for modelling the dynamic of some the COVID-19 count data. In this regard, here, we explore some more possibilities in modeling the dynamics of such data.

More precisely, we observed two real data sets, one of which refers to the period before and the other after the appearance of the SARS-CoV2 virus infection on the territory of the Republic of Serbia. The first time series, labelled as Series A, contains the number of daily mortality in Niš, the third largest city in the Republic of Serbia, in the period from 1 January 2017 until 31 December 2019. The second one (Series B) contains the daily number of deaths per day in the Republic of Serbia, from 20 March 2020 (the date of the first death caused by the disease COVID-19) to 1 December 2022. Dynamics of the observed time series, lengths $T=1095$ and $T=990$, respectively, are shown on the left plots in Figure 4. On the other hand, the right plots in Figure 4 show the autocorrelation functions (ACFs) of these time series. In the case of Series A, it is obvious that there is a very weak correlation, so it can be seen as an IID time series. On the contrary, with Series B, it is evident that there is a positive and significant but exponentially decreasing correlation, depending on the lag $k=0,1,2,\cdots$. Therefore, in accordance with the previous theoretical results, primarily Equation (20), it could be expected that the INAR(1) process is an adequate stochastic model for describing the dynamics of the number of deaths caused by the COVID-19 disease.

Some confirmation of these observations can also be seen in Table 4, where its upper part shows the basic descriptive statistical analysis of the both time series. In the case of Series A, the existence of an equal-dispersion property can be recognized through ’close’ values of its sample mean and variance. Thus, this series can be modeled with a Poisson distribution. On the other hand, with Series B, there is a wide range and variability of observed data. Additionally, it can be seen that frequencies of ‘smaller’ numerical values are emphasized (its mode equals 1). Therefore, it could be assumed that, in the INAR(1) model, innovations with a Poisson distribution (and a ‘small’ parameter $\theta$) can be taken, as well as with a geometric distribution. In accordance with the mentioned facts, we modeled the dynamics of the Series B by using the INAR(1) process $\left(X_t\right)$ with geometric distributed innovations $\left({\epsilon }_t\right)$. Finally, the lower part of Table 4 shows some of the estimated values for the ACFs of both observed time series. Obviously, the obtained autocorrelation values of Series A are not statistically significant, so it can be considered as an IID time series. Conversely, the estimated ACF values of Series B indicate the possibility of modeling the dynamics of this time series using the INAR(1) process.

In the upper part of Table 5, the estimated values of parameter $\theta >0$ for Series A, as well as estimates $(\theta ,\alpha )$ of the INAR(1) process for Series B, are shown. As in the previous simulation study, PGF estimates for IID Series A were obtained using an unconstrained Brent’s optimization procedure. On the other hand, for Series B, PGF estimates were computed using the two-step estimation procedure described in the previous section, that is, by linearly constrained optimization and YW estimates as initial parameter values. In both cases, as before, weight functions ${\omega }_k(u_1,u_2)$, $k=0,1,2$ based on Gegenbauer orthogonal polynomials were used.

It can be easily seen that in the case of Series A, estimated parameter $\theta$ is ’between’ to sample mean and variance, which is expected. For Series B, the estimated values of the parameter $\theta$ are relatively ‘small’. This is obviously a consequence of the emphasized frequencies of smaller integer values of the observed time series (which can also be seen in the following Figure 5). Conversely, the estimated values of the parameter $\alpha$ are close to 1. This can also be explained by the previously described properties of this time series, which has a significantly high estimated value of the first correlation.

Furthermore, we analyzed the fitting efficiency for both time series, when all weights are applied. For this purpose, we generated 500 independent simulations for the IID Series A, as well as the INAR(1) model for Series B, for all estimated parameters’ values. In order to check the efficiency of proposed models, we computed two typical goodness-of-fit statistics: the mean squared estimation error (MSEE) and the Akaike information criterion (AIC). The mean values of these statistics are shown in the middle part of Table 5. It is noticeable that, in all the mentioned cases, the PGF estimates have similar efficiency because the MSEE and AIC statistics have relatively close and small estimated values.

Finally, an analysis of the forecast accuracy of models obtained with PGF estimates was also checked. For this purpose, we used the one-tailed Diebold–Mariano test [49] for prediction accuracy. The null hypothesis is that the fitted time series have the same accuracy as the actual data, while the alternative is that the series obtained by the PGF estimation procedure are less accurate. The test statistic (denoted as DM), as well as the corresponding p-values, were calculated using the R-package “forecast”, authored by Hyndman [50]. In the case of Series A, the data set from 1 January 2020 to 20 March 2020 was used, as a forecast horizon of length $h=80$. For Series B, the forecast horizon from 1 December 2022 to 9 January 2023 was used, the length of which is equal to $h=40$. As a result, it is noticeable that the realized values of the DM statistics, shown in the lower part of Table 5, indicate that all models fitted by PGF estimators have the same forecast accuracy as empirical data, i.e., the null hypothesis is valid. Some of these aforementioned facts can be also seen in Figure 5, where the empirical frequencies distribution, along with the appropriate distributions of IID and INAR(1) processes, are shown.

5. Conclusions

In this paper, an estimation method based on probability generating functions (PGFs) is proposed. For this purpose, the stochastic properties of this method were analyzed, as well as its application in estimating the parameters of some typical stationary time series. Notice that the PGF method, with certain modifications, may be suitable for estimating and fitting time series of a similar, or some other kind to those considered here. These can be, for instance, INAR processes with negative binomial thinning [51,52], integer generalized autoregressive (IN-GARCH) processes [53], or the so-called self-exciting processes [54]. Moreover, the application of the PGF estimates is also particularly interesting in the case of some non-linear time series, when the usual estimation methods cannot be used. We point out again that, in Stojanović et al. [20], the PGF method was applied in parameter estimation of the so-called noise-indicator INAR process. Consequently, the application of the PGF estimators presented here can be a starting point for some future research related to these estimators. Therefore, notice that, in the case of Series B, the estimated values of the parameter $\alpha$ indicate that the hypothesis $H_0:\alpha =1$ is potentially valid. This could indicate a possible non-stationarity in the dynamics of the observed time series and could be the subject of some further research.

Author Contributions

Conceptualization, V.S. and E.L.; methodology, V.S. and M.T.; software, V.S. and E.L.; validation, E.L.; formal analysis, V.S.; data curation, V.S. and M.T.; writing—original draft preparation, V.S., E.L. and M.T.; writing—review and editing, E.L.; visualization, V.S. and M.T.; supervision, V.S.; project administration, E.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Official data on deaths from the disease COVID-19 can be found on the website of the Office for IT and e-Government of the Republic of Serbia (http://covid19.data.gov.rs/?locale=en, accessed on 6 December 2022).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

ACF	Auto Correlation Function
AD test	Anderson–Darling test
AIC	Akaike Information Criterion
AN	Asymptotic Normality
CF	Characteristic Function
ECF	Empirical Characteristic Function
IID series	Independent Identical Distributed series
INAR process	Integer-Valued Autoregressive process
INMA	Integer Moving Average
MSEE	Mean Squared Estimating Error
NIINAR process	Noise-Indicator Integer-Valued Autoregressive process
NNIV series	Non-Negative Integer-Valued series
PGF	Probability Generating Function
PMF	Probability Mass Distribution
PS	Power Series
RV	Random Variable
SLLN	Strong Law of Large Numbers
YW	Yule–Walker

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Graphs of PGFs (lines) and empirical PGFs (dots) of the IID series: (a) Poisson distribution; (b) geometric distribution. (True parameter value is [Forumla omitted. See PDF.]).

Figure 2. Plots of the sample frequency distributions of INAR(1) processes (a) and their innovations (b) with Poisson distribution (plots above) and geometric distribution (plots below). True parameters’ values are [Forumla omitted. See PDF.].

Figure 3. Three-dimensional plots of the two-dimensional PGFs (a) and empirical PGFs (b) of the INAR(1) series with Poisson innovations (panels above) and geometric innovations (panels below). True parameters’ values are [Forumla omitted. See PDF.].

Figure 4. (a) Dynamics of the number of deaths before and after the onset of the disease COVID-19 in Serbia. (b) The corresponding ACFs of the observed time series.

Figure 5. Frequency distribution of the actual and fitted data, obtained by (a) IID Poisson process; (b) INAR(1) process with geometric innovations.

Table 1

NNIV distributions, along with their PGFs, interpreted in the form of PS distributions.

Distributions	$\mathit{S}$	$\mathit{a}\left(\mathit{x}\right)$	$(0,\mathit{R})$	$\mathit{f}\left(\mathit{\theta }\right)$	${\mathbf{\Psi }}_{\mathit{\epsilon }}(\mathit{u};\mathit{\theta })$	$\mathit{R}/\mathit{\theta }$
1. Bernoulli	$\left\{0,1\right\}$	1	$(0,\infty )$	$1+\theta$	${\displaystyle \frac{1+\theta u}{1+\theta }}$	∞
2. Binomial	$\left\{0,\cdots ,n\right\}$	${\displaystyle \left(\genfrac{}{}{0pt}{}{n}{x}\right)}$	$(0,\infty )$	${(1+\theta )}^n$	${\displaystyle {\left(\frac{1+\theta u}{1+\theta }\right)}^n}$	∞
3. Poisson	$\left\{0,\cdots ,\infty \right\}$	${\displaystyle \frac{1}{x!}}$	$(0,\infty )$	$exp\left(\theta \right)$	$exp\left(\theta (u-1)\right)$	∞
4. Geometric	$\left\{0,\cdots ,\infty \right\}$	1	$(0,1)$	${\displaystyle \frac{1}{1-\theta }}$	${\displaystyle \frac{1-\theta }{1-\theta u}}$	$1/\theta$
5. Neg. Binomial	$\left\{0,\cdots ,\infty \right\}$	${\displaystyle \frac{\Gamma (x+n)}{x!\Gamma \left(n\right)}}$	$(0,1)$	${\displaystyle \frac{1}{{(1-\theta )}^n}}$	${\displaystyle {\left(\frac{1-\theta }{1-\theta u}\right)}^n}$	$1/\theta$
6. Pascal	$\left\{n,\cdots ,\infty \right\}$	${\displaystyle \left(\genfrac{}{}{0pt}{}{x-1}{n-1}\right)}$	$(0,1)$	${\displaystyle {\left(\frac{\theta }{1-\theta }\right)}^n}$	${\displaystyle {\left(\frac{(1-\theta )u}{1-\theta u}\right)}^n}$	$1/\theta$

Table 2

Estimated parameters’ values of IID series $\left({\epsilon }_t\right)$ series with Poisson and geometric distribution. (True parameter value is $\theta =0.5$.)

Statistics		PGF Estimates
		${\widehat{\mathit{\theta }}}_0$	${\widehat{\mathit{S}}}_{{\mathit{T}}_0}$		${\widehat{\mathit{\theta }}}_1$	${\widehat{\mathit{S}}}_{{\mathit{T}}_1}$		${\widehat{\mathit{\theta }}}_2$	${\widehat{\mathit{S}}}_{{\mathit{T}}_2}$
Poisson	Min.	0.4222	$4.63\times {10}^{-8}$		0.4293	$3.20\times {10}^{-9}$		0.4267	$3.83\times {10}^{-8}$
	Mean	0.5006	$9.19\times {10}^{-5}$		0.4994	$4.47\times {10}^{-5}$		0.5005	$2.80\times {10}^{-5}$
	Max.	0.5967	$7.17\times {10}^{-4}$		0.6098	$4.00\times {10}^{-4}$		0.5938	$2.38\times {10}^{-4}$
	MSEE	0.0220	–		0.0200	–		0.0196	–
	$AD$	0.3849	–		0.1965	–		0.2911	–
	p-value	0.3919	–		0.8889	–		0.6076	–
Geometric	Min.	0.4491	$1.09\times {10}^{-6}$		0.4520	$3.27\times {10}^{-7}$		0.4467	$2.40\times {10}^{-7}$
	Mean	0.5004	$2.30\times {10}^{-4}$		0.5001	$1.06\times {10}^{-4}$		0.4996	$6.17\times {10}^{-5}$
	Max.	0.5535	$2.71\times {10}^{-3}$		0.5353	$3.62\times {10}^{-3}$		0.5580	$5.60\times {10}^{-4}$
	MSEE	0.0129	–		0.0121	–		0.0117	–
	$AD$	0.4583	–		0.2214	–		0.2490	–
	p-value	0.2627	–		0.8307	–		0.7467	–

Table 3

Estimated parameter values of the INAR(1) process $\left(X_t\right)$ with Poisson and geometric innovations. (True parameters’ values are $\theta =\alpha =1/2$).

Statistics		PGF Estimates
		${\widehat{\mathit{\theta }}}_0$	${\widehat{\mathit{\alpha }}}_0$	${\mathit{S}}_{{\mathit{T}}_0}^{\left(2\right)}$		${\widehat{\mathit{\theta }}}_1$	${\widehat{\mathit{\alpha }}}_1$	${\mathit{S}}_{{\mathit{T}}_1}^{\left(2\right)}$		${\widehat{\mathit{\theta }}}_2$	${\widehat{\mathit{\alpha }}}_2$	${\mathit{S}}_{{\mathit{T}}_2}^{\left(2\right)}$
Poisson	Min.	0.3845	0.3606	$2.12\times {10}^{-5}$		0.3653	0.3352	$5.48\times {10}^{-6}$		0.3710	0.3365	$1.65\times {10}^{-6}$
	Mean	0.5003	0.4986	$8.08\times {10}^{-4}$		0.5002	0.5006	$1.73\times {10}^{-4}$		0.5015	0.4951	$7.36\times {10}^{-5}$
	Max.	0.6711	0.5971	$7.12\times {10}^{-3}$		0.6948	0.6292	$1.32\times {10}^{-3}$		0.7049	0.6251	$5.13\times {10}^{-4}$
	MSEE	0.0329	0.0350	–		0.0364	0.0337	–		0.0381	0.0321	–
	$AD$	0.7468	0.3089	–		0.7272	0.4153	–		0.8609	0.4724	–
	p-value	0.0515	0.5574	–		0.0576	0.3326	–		0.0269 *	0.2427	–
Geometric	Min.	0.3962	0.3592	$3.96\times {10}^{-5}$		0.3673	0.3129	$7.98\times {10}^{-6}$		0.3965	0.3199	$2.75\times {10}^{-6}$
	Mean	0.5002	0.4974	$7.31\times {10}^{-4}$		0.4981	0.4977	$2.23\times {10}^{-4}$		0.4986	0.5005	$8.98\times {10}^{-5}$
	Max.	0.6071	0.6646	$9.02\times {10}^{-3}$		0.5960	0.6037	$1.57\times {10}^{-3}$		0.6033	0.6747	$7.39\times {10}^{-4}$
	MSEE	0.0291	0.0521	–		0.0265	0.0446	–		0.0271	0.0425	–
	$AD$	0.6155	0.8198	–		0.3150	0.5174	–		0.4126	0.7898	–
	p-value	0.1088	0.0340 *	–		0.5426	0.1882	–		0.3375	0.0403 *	–

* p < 0.05.

Table 4

Summary statistics of the total number of deaths per day before and after the appearance of the COVID-19 disease in the Republic of Serbia.

Statistics		Series A		Series B
Sample size		1095		990
Minimum		0		0
Mode		6		1
1st Quartile		4		3
Median		6		9
Mean		5.758		17.58
3rd Quartile		7		27
Maximum		24		79
St. deviation		2.428		18.48
Variance		5.895		341.6
Skewness		1.797		1.159
Kurtosis		10.803		3.161
ACF(1)		0.085		0.987
⋯		⋯		⋯
ACF(10)		0.020		0.911
⋯		⋯		⋯
ACF(20)		$-0.024$		0.754
⋯		⋯		⋯
ACF(50)		0.049		0.254
⋯		⋯		⋯
ACF(100)		$-0.020$		0.078

Table 5

Estimated values of parameters, along with error and predictive statistics of the actual data.

Parameters/Statistics		Series A			Series B
		${\mathbf{\omega }}_{\mathbf{0}}$	${\mathbf{\omega }}_{\mathbf{1}}$	${\mathbf{\omega }}_{\mathbf{2}}$	${\mathbf{\omega }}_{\mathbf{0}}$	${\mathbf{\omega }}_{\mathbf{1}}$	${\mathbf{\omega }}_{\mathbf{2}}$
$\theta$		5.8069	5.8218	5.8176	0.1812	0.1808	0.1810
$\alpha$		–	–	–	0.9872	0.9865	0.9871
$S_T^{\left(2\right)}$		$1.06\times {10}^{-3}$	$1.06\times {10}^{-3}$	$7.16\times {10}^{-5}$	$6.48\times {10}^{-3}$	$3.46\times {10}^{-3}$	$1.92\times {10}^{-3}$
MSEE		0.0726	0.0789	0.0771	0.0221	0.0288	0.0233
AIC		4.5992	4.6023	4.6016	5.3035	5.3041	5.3042
$DM$		−0.3641	0.1039	0.0974	0.5855	0.7797	0.6456
p-value		0.3579	0.5414	0.5388	0.7208	0.7821	0.7406

References

1. Johnson, N.L.; Kemp, A.W.; Kotz, S. Univariate Discrete Distributions; 3rd ed. Wiley Series in Probability and Statistics; Wiley: Washington, DC, USA, 2005.

2. Xu, H.-Y.; Xie, M.; Goh, T.N.; Fu, X. A Model for Integer–Valued Time Series With Conditional Overdispersion. Comput. Stat. Data Anal.; 2012; 56, pp. 4229-4242. [DOI: https://dx.doi.org/10.1016/j.csda.2012.04.011]

3. Jazi, M.A.; Jones, G.; Lai, C.-D. First-Order Integer Valued AR Processes With Zero Inflated Poisson Innovations. J. Time Ser. Anal.; 2012; 33, pp. 954-963. [DOI: https://dx.doi.org/10.1111/j.1467-9892.2012.00809.x]

4. Weiß, C.H.; Pollet, P.K. Binomial Autoregressive Processes With Density-Dependent Thinning. J. Time Ser. Anal.; 2014; 35, pp. 115-132. [DOI: https://dx.doi.org/10.1002/jtsa.12054]

5. Graziadei, H.; Lijoi, A.; Lopes, H.F.; Marques, P.C.; Prünster, I. Prior Sensitivity Analysis in a Semi-Parametric Integer-Valued Time Series Model. Entropy; 2020; 22, 69. [DOI: https://dx.doi.org/10.3390/e22010069] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33285844]

6. Khoo, W.C.; Ong, S.H.; Atanu, B. Coherent Forecasting for a Mixed Integer-Valued Time Series Model. Mathematics; 2022; 10, 2961. [DOI: https://dx.doi.org/10.3390/math10162961]

7. El-Morshedy, M. A Discrete Linear-Exponential Model: Synthesis and Analysis with Inference to Model Extreme Count Data. Axioms; 2022; 11, 531. [DOI: https://dx.doi.org/10.3390/axioms11100531]

8. Du, J.-G.; Li, Y. The Integer-Valued Autoregressive (INAR(p)) Model. J. Time Ser. Anal.; 1991; 12, pp. 129-142.

9. Franke, J.; Seligmann, T.H. Conditional Maximum Likelihood Estimates for INAR(1) Processes and Their Application to Modelling Epileptic Seizure Counts. Developments in Time Series; Chapman & Hall: London, UK, 1993; pp. 310-330.

10. Latour, A. Existence and Stochastic Structure of a Non-negative Integer-valued Autoregressive Process. J. Time Ser. Anal.; 1998; 19, pp. 439-455. [DOI: https://dx.doi.org/10.1111/1467-9892.00102]

11. Silva, I.; Silva, M.E. Asymptotic Distribution of the Yule-Walker Estimator for INAR(p) Processes. Stat. Prob. Lett.; 2006; 76, pp. 1655-1663. [DOI: https://dx.doi.org/10.1016/j.spl.2006.04.008]

12. Silva, I.; Silva, M.E. Parameter Estimation for INAR Processes Based on High-Order Statistics. Revstat; 2009; 7, pp. 105-117.

13. Martin, V.L.; Tremayne, A.R.; Jung, R.C. Efficient Method of Moments Estimators for Integer Time Series Models. J. Time Ser. Anal.; 2014; 35, pp. 491-516. [DOI: https://dx.doi.org/10.1111/jtsa.12078]

14. Knight, J.L.; Yu, J. Empirical Characteristic Function in Time Series Estimation. Econom. Theory; 2002; 18, pp. 691-721. [DOI: https://dx.doi.org/10.1017/S026646660218306X]

15. Yu, J. Empirical Characteristic Function Estimation and Its Applications. Econom. Rev.; 2004; 23, pp. 93-123. [DOI: https://dx.doi.org/10.1081/ETC-120039605]

16. Kotchoni, R. Applications of the Characteristic Function-Based Continuum GMM in Finance. Comput. Stat. Data Anal.; 2012; 56, pp. 3599-3622. [DOI: https://dx.doi.org/10.1016/j.csda.2010.08.011]

17. Meintanis, S.G.; Swanepoel, J.; Allison, J. The probability weighted characteristic function and goodness-of-fit testing. J. Stat. Plan. Infer.; 2014; 146, pp. 122-132. [DOI: https://dx.doi.org/10.1016/j.jspi.2013.09.011]

18. Stojanović, V.; Milovanović, G.V.; Jelić, G. Distributional Properties and Parameters Estimation of GSB Process: An Approach Based on Characteristic Functions. ALEA Lat. Am. J. Probab. Math. Stat.; 2016; 13, pp. 835-861. [DOI: https://dx.doi.org/10.30757/ALEA.v13-33]

19. Esquìvel, M.L. Some Applications of Probability Generating Function Based Methods to Statistical Estimation. Discuss. Math.; 2009; 29, pp. 131-153. [DOI: https://dx.doi.org/10.7151/dmps.1111]

20. Stojanović, V.; Randjelović, D.; Kuk, K. Noise-Indicator Nonnegative Integer-Valued Autoregressive Time Series of the First Order. Braz. J. Probab. Stat.; 2018; 32, pp. 147-171. [DOI: https://dx.doi.org/10.1214/16-BJPS335]

21. Newey, W.K.; McFadden, D. Large Sample Estimation and Hypothesis Testing. Handbook of Econometrics; Elsevier: Amsterdam, The Netherlands, 1994; Volume 4.

22. Fuller, W.A. Introduction to Statistical Time Series; John Wiley & Sons: New York, NY, USA, 1996.

23. Bourguignon, M.; Vasconcellos, K.L.P. First Order Non-Negative Integer Valued Autoregressive Processes With Power Series Innovations. Braz. J. Probab. Stat.; 2015; 29, pp. 71-93. [DOI: https://dx.doi.org/10.1214/13-BJPS229]

24. Cvetković, A.S.; Milovanović, G.V. The Mathematica Package “Orthogonal Polynomials”. Facta Univ. Ser. Math. Inform.; 2004; 19, pp. 17-36.

25. Brent, R. Algorithms for Minimization Without Derivatives; Prentice-Hall: Englewood Cliffs, NJ, USA, 1973.

26. Gross, L. Tests for Normality. R Package Version 1.0-2. 2013; Available online: http://CRAN.R-project.org/package=nortest (accessed on 6 December 2022).

27. Weiß, C.H. Thinning Operations for Modelling Time Series of Counts—A Survey. Adv. Statist. Anal.; 2008; 92, pp. 319-341. [DOI: https://dx.doi.org/10.1007/s10182-008-0072-3]

28. Billingsley, P. Probability and Measure; John Wiley & Sons: New York, NY, USA, 1995.

29. Silva, M.E.; Oliveira, V.L. Difference Equations for the Higher-Order Moments and Cumulants of the INAR(1) Model. J. Time Ser. Anal.; 2004; 25, pp. 317-333. [DOI: https://dx.doi.org/10.1111/j.1467-9892.2004.01685.x]

30. Jung, R.C.; Ronning, G.; Tremayne, A.R. Estimation in Conditional First Order Autoregression with Discrete Support. Stat. Pap.; 2005; 46, pp. 195-224. [DOI: https://dx.doi.org/10.1007/BF02762968]

31. Ristić, M.M.; Nastić, A.S.; Miletić, A.V. A Geometric Time Series Model with Dependent Bernoulli Counting Series. J. Time Ser. Anal.; 2013; 34, pp. 466-476. [DOI: https://dx.doi.org/10.1111/jtsa.12023]

32. Schweer, S.; Weiß, C.H. Compound Poisson INAR(1) Processes: Stochastic Properties and Testing for Over-Dispersion. Comput. Stat. Data Anal.; 2014; 77, pp. 267-284. [DOI: https://dx.doi.org/10.1016/j.csda.2014.03.005]

33. Bermúdez, L.; Karlis, D. Multivariate INAR(1) Regression Models Based on the Sarmanov Distribution. Mathematics; 2021; 9, 505. [DOI: https://dx.doi.org/10.3390/math9050505]

34. Li, Q.; Chen, H.; Liu, X. A New Bivariate Random Coefficient INAR(1) Model with Applications. Symmetry; 2022; 14, 39. [DOI: https://dx.doi.org/10.3390/sym14010039]

35. Maya, R.; Chesneau, C.; Krishna, A.; Irshad, M.R. Poisson Extended Exponential Distribution with Associated INAR(1) Process and Applications. Stats; 2022; 5, pp. 755-772. [DOI: https://dx.doi.org/10.3390/stats5030044]

36. Maya, R.; Irshad, M.R.; Chesneau, C.; Nitin, S.L.; Shibu, D.S. On Discrete Poisson–Mirra Distribution: Regression, INAR(1) Process and Applications. Axioms; 2022; 11, 193. [DOI: https://dx.doi.org/10.3390/axioms11050193]

37. Al-Osh, M.A.; Alzaid, A.A. First-Order Integer-Valued Autoregressive (INAR(1)) Process. J. Time Ser. Anal.; 1987; 8, pp. 261-275. [DOI: https://dx.doi.org/10.1111/j.1467-9892.1987.tb00438.x]

38. Alzaid, A.A.; Al-Osh, M.A. An Integer-Valued p^th-order Autoregressive Structure (INAR(p)) Process. J. App. Prob.; 1990; 27, pp. 314-324. [DOI: https://dx.doi.org/10.2307/3214650]

39. Lange, K. Numerical Analysis for Statisticians. Book Series: Statistics and Computing; Springer: New York, NY, USA, 2001.

40. Stojanović, V.; Popović, B.Č.; Milovanović, G.V. The Split-SV model. Comput. Statist. Data Anal.; 2016; 100, pp. 560-581. [DOI: https://dx.doi.org/10.1016/j.csda.2014.08.010]

41. Vaz, S.; Torres, D.F.M. A Discrete-Time Compartmental Epidemiological Model for COVID-19 with a Case Study for Portugal. Axioms; 2021; 10, 314. [DOI: https://dx.doi.org/10.3390/axioms10040314]

42. Ghosh, S.; Volpert, V.; Banerjee, M. An Epidemic Model with Time Delay Determined by the Disease Duration. Mathematics; 2022; 10, 2561. [DOI: https://dx.doi.org/10.3390/math10152561]

43. Sivakumar, B.; Deepthi, B. Complexity of COVID-19 Dynamics. Entropy; 2022; 24, 50. [DOI: https://dx.doi.org/10.3390/e24010050] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35052076]

44. Hassan, S.M.; Riveros Gavilanes, J.M. First to React Is the Last to Forgive: Evidence from the Stock Market Impact of COVID 19. J. Risk Financ. Manag.; 2021; 14, 26. [DOI: https://dx.doi.org/10.3390/jrfm14010026]

45. Zakharov, V.; Balykina, Y.; Ilin, I.; Tick, A. Forecasting a New Type of Virus Spread: A Case Study of COVID-19 with Stochastic Parameters. Mathematics; 2022; 10, 3725. [DOI: https://dx.doi.org/10.3390/math10203725]

46. Jovanović, M.; Stojanović, V.; Kuk, K.; Popović, B.; Čisar, P. Asymptotic Properties and Application of GSB Process: A Case Study of the COVID-19 Dynamics in Serbia. Mathematics; 2022; 10, 3849. [DOI: https://dx.doi.org/10.3390/math10203849]

47. Stapper, M. Count Data Time Series Modelling in Julia—The CountTimeSeries.jl Package and Applications. Entropy; 2021; 23, 666. [DOI: https://dx.doi.org/10.3390/e23060666]

48. Eliwa, M.S.; Tyagi, A.; Almohaimeed, B.; El-Morshedy, M. Modelling Coronavirus and Larvae Pyrausta Data: A Discrete Binomial Exponential II Distribution with Properties, Classical and Bayesian Estimation. Axioms; 2022; 11, 646. [DOI: https://dx.doi.org/10.3390/axioms11110646]

49. Diebold, F.X.; Mariano, R.S. Comparing Predictive Accuracy. J. Bus. Econ. Stat.; 1995; 13, pp. 253-263.

50. Hyndman, R. Forecasting Functions for Time Series and Linear Models. R Package Version 7.1. 2016; Available online: http://CRAN.R-project.org/package=forecast (accessed on 9 January 2023).

51. Nastić, A.S.; Ristić, M.M.; Bakouch, H.S. A Combined Geometric INAR(p) Model Based on Negative Binomial Thinning. Math. Comput. Model.; 2012; 55, pp. 1665-1672. [DOI: https://dx.doi.org/10.1016/j.mcm.2011.10.080]

52. Ristić, M.M.; Nastić, A.S.; Bakouch, H.S. Estimation in an Integer-Valued Autoregressive Process With Negative Binomial Marginals (NBINAR(1)). Comm. Stat. Theory Methods; 2012; 41, pp. 606-618. [DOI: https://dx.doi.org/10.1080/03610926.2010.529528]

53. Weiß, C. INGARCH Models for Count Time Series. An Introduction to Discrete-Valued Time Series; Weiß, C. Elsevier: Amsterdam, The Netherlands, 2018; [DOI: https://dx.doi.org/10.1002/9781119097013.ch4]

54. Reinhart, A. A Review of Self-Exciting Spatio-Temporal Point Processes and Their Applications. Stat. Sci.; 2018; 33, pp. 299-318.