First–order integer–valued autoregressive processes with Generalized Katz innovations

The probability mass function (pmf), $P(X=x)=p_{x}$, of the Generalized Lagrangian Katz (GLK) is

$$p_{x}=\frac{1}{x!}\beta^{x}\frac{a}{c}\frac{1}{(\frac{a}{c}+x\frac{b}{c}+x)}(1-\beta)^{\frac{a}{c}+x\frac{b}{c}}\left(\frac{a}{c}+x\frac{b}{c}+1\right)_{{x}\uparrow}$$

(1)

$x=0,1,2,\ldots$, where $(x)_{{k}\uparrow}=x(x+1)\ldots(x+k-1)$ is the rising factorial with the convention that $(x)_{0}=1$, and $a>0$, $c>0$, $b\geq-c$ and $0<\beta<1$ are the parameters [Consul & Famoye, 2006]. We denote the distribution with $\mathcal{GLK}(a,b,c,\beta)$. We notice that for $-c<b<0$ some additional constraints on the parameters are needed to have all the $p_{x}\geq 0$. See the discussion at the beginning of Subsection 3.1 and Appendix B in the Supplementary. GLK distributions have probability generating function (pgf)

$$H(u)=\sum_{x=0}^{\infty}p_{x}u^{x}$$

which satisfies:

$$H(u)=(1-\beta+\beta z)^{a/c},\,\,z=u(1-\beta+\beta z)^{b/c+1},$$

(2)

or alternatively

$$H(u)=\left((1-\beta)/(1-\beta z)\right)^{a/c},\,\,z=u\left((1-\beta z)/(1-\beta)\right)^{b/c},$$

(3)

see Janardan [1998].

The GLK distribution family is very general and includes some well–known distributions and new distributions that have yet to be used in count data modelling.

• •

  The Lagrangian Katz distribution $\mathcal{LK}(a,b,\beta)$ is obtained by replacing $c$ with $\beta$ (which is called Generalized Katz in [Consul & Famoye, 2006]).

• •

  The Katz distribution $\mathcal{K}(a,\beta)$ is obtained for $b=0$ and by replacing $c$ with $\beta$, [Katz, 1965].

• •

  The Polya–Eggenberger distribution $\mathcal{PE}(a,c,\beta)$ is obtained for $b=0$, [Janardan, 1998]. Note that the Katz distribution in Consul & Famoye [2006], Tab. 2.1, is not the Katz distribution of Katz [1965], it corresponds instead to the Generalized Polya Eggenberger of the first type (GPED₁–I) of Janardan [1998] and can be obtain as the limit of the zero–truncated GLK for $a\rightarrow-c$.

• •

  The Generalized Negative Binomial distribution $\mathcal{GNB}(r,\gamma,p)$ is obtained for $c=1$, $a=r$, $b=\gamma-1$ and $\beta=p$ .

• •

  The Negative Binomial distribution $\mathcal{NB}(r,p)$ is obtained for $b=0$ $\beta=1-p$ and $r=a/c$.

• •

  The Binomial distribution $\mathcal{B}in(n,p)$ is obtained for $c=1$, $b=-1$, $a=n\in\mathbb{N}$ and $\beta=p$

• •

  The Generalized Poisson (GP) distribution $\mathcal{GP}(\theta,\lambda)$ for $c\rightarrow 0$ s.t. $b/a=\lambda$ and $a\beta/c=\theta>0$ with $0<\lambda<\theta^{-1}$. The GP limit of the GLK distribution is stated in [Consul & Famoye, 2006] without proof. In Appendix B.2 in the Supplementary Material, we provide a proof.

• •

  The Poisson distribution $\mathcal{P}(\theta)$ for $c\rightarrow 0$, $b\rightarrow 0$ s.t. $a\beta/c=\theta$.

The probability mass function of the GLK for different parameter settings is given in Fig. 1. In the top–left plot, we compare $\mathcal{K}(a,c)$, $\mathcal{LK}(a,b,\beta)$ and $\mathcal{GLK}(a,b,c,\beta)$ with the same mean. The top–right plot illustrates the sensitivity of the $\mathcal{GLK}(a,b,c,\beta)$ pmf with respect to the different parameters. All distributions have the same mean (vertical dashed line). The bottom plots illustrate the effects of the parameters on the tails (log–scale) for a $\mathcal{GLK}(a,b,c,\beta)$ with over–dispersion $VMR=50/15$ (left) and under–dispersion $VMR=13/15$ (right).

We provide in Appendix A in the Supplementary Material some useful moments of the GLK distributions, which can be used to derive the following measures of dispersion. The standard deviation to the mean ratio returns the coefficient of variation. From the results in Appendix A in the Supplementary Material it follows that the coefficient of variation is $CV=\left((1-\beta)/(a\theta\kappa)\right)^{1/2}$ where $\kappa=1-\beta-b\beta/c$ and $\theta=\beta/c$, assuming $\kappa>0$. The Fisher index is given by the variance–to–mean ratio $VMR=(1-\beta)/(\kappa^{2})$ which does not depend on the parameter $a$. For a given $\beta$, following the values of $\kappa$ ($b$ and $c$), the distribution allows for various degrees of dispersion: not dispersed ($VMR=0$), under–dispersed ($VMR<1$), equally dispersed ($VMR=1$) and over–dispersed ($VMR>1$).

We conclude this section with another important property.

The Generalized Katz INAR(1) process (GLK–INAR(1)) is defined using the binomial thinning operator, $\circ$. The binomial thinning for a non–negative discrete random variable $X$ is defined as

$$\alpha\circ X=\sum_{i=1}^{X}B_{i}(\alpha)$$

where $B_{i}(\alpha)$ are iid Bernoulli r.v.s with success probability $P(B_{i}(\alpha)=1)=\alpha$.

Figure 2 provides some trajectories of $T=100$ points each, simulated from a GLK–INAR(1) with innovation distributions given by the solid lines in the bottom plots of Fig. 1, that are $\mathcal{GLK}(3.86$, $0$, $0.60$, $0.70)$ (overdispersion) and $\mathcal{GLK}(25.00$, $0.00$, $0.70$, $0.42)$ (underdispersion). The trajectories correspond to the two parameter settings we find the empirical application to climate change discussed in Section 4, that are: (i) high persistence setting ($\alpha=0.7$, left); (ii) low persistence setting ($\alpha=0.3$, right). In all plots, the empirical mean of the observations is reported (dashed line) as a reference to illustrate the different levels of persistence in the trajectories.

Thanks to the general parametric family assumed, by setting $b=0$, $c=\beta=\theta_{1}$ and $a=\theta_{2}$, our GLK–INAR(1) nests the INARKF(1) of Kim & Lee [2017] as special case. The GLK–INAR(1) naturally nests the Poisson INAR(1) of Al-Osh & Alzaid [1987], the Negative Binomial INAR(1) of Al-Osh & Aly [1992] (NBINAR(1)), and the Generalized Poisson INAR(1) of Alzaid & Al-Osh [1993].

As for any INAR process, the GLK–INAR(1) has the following representation

$$X_{t+k}=\alpha^{k}\circ X_{t}+\sum_{j=0}^{k-1}\alpha^{j}\circ\varepsilon_{t+1-j}$$

and its conditional pgf can be written as

$$H_{X_{t+k}|X_{t}}(u)=(1-\alpha^{k}+\alpha^{k}u)^{X_{t}}\prod_{j=0}^{k-1}H(1-\alpha^{j}+\alpha^{j}u)$$

where $H(u)$ is defined in Eq. 2 or in Eq. 3. Starting from the general results on INAR processes given in Alzaid & Al-Osh [1988], one easily obtains explicit expressions for the conditional mean and variance of the GLK–INAR(1):

	$\displaystyle\mathbb{E}(X_{t+k}|X_{t})=\alpha^{k}X_{t}+\frac{1-\alpha^{k-1}}{1-\alpha}\frac{a\theta}{\kappa}$		(5)
	$\displaystyle\mathbb{V}(X_{t+k}|X_{t})=\frac{1-\alpha^{2k}}{1-\alpha^{2}}\left(\frac{a(1-\beta)\theta}{\kappa^{3}}-\frac{a\theta}{\kappa}\right)$
	$\displaystyle\hskip 70.0pt+(\alpha^{k}-\alpha^{2k})X_{t}+\frac{1-\alpha^{k}}{1-\alpha}\frac{a\theta}{\kappa}$		(6)

where $\kappa=1-\beta-b\beta/c$ and $\theta=\beta/c$.

The process $\{X_{t}\}_{t\in\mathbb{Z}}$ is a Markov Chain on $\mathbb{N}$ and the transition probability $P_{i,j}=\mathbb{P}(X_{t}=j|X_{t-1}=i)$ satisfies

	$\displaystyle P_{i,j}=\sum_{k=0}^{\min(i,j)}\mathbb{P}(\alpha\circ X_{t-1}=k|X_{t-1}=i)\mathbb{P}(\varepsilon=j-k)$
	$\displaystyle=\sum_{k=0}^{\min(i,j)}\binom{i}{k}\alpha^{k}(1-\alpha)^{i-k}p_{j-k}$

where $p_{x}$ is the pmf given in Eq. 1.

In the next Proposition, we summarize some of the asymptotic properties of a GLK–INAR(1). These properties follow from general results in Alzaid & Al-Osh [1988] and Schweer & Weiß [2014].

Since the GLK distribution satisfies the convolution property [see Janardan, 1998, Th. 8], then the GLK–INAR(1) is stable by aggregation as stated in the following

Below, we state explicit closed-form expressions for unconditional moments of the process.

From the previous proposition, under the assumption $\kappa=1-\beta-b\beta/c>0$ one obtains the unconditional variance of the process $\sigma^{2}_{X}=\mathbb{V}(X_{t})=(\sigma_{\varepsilon}^{2}+\alpha\mu_{\varepsilon})/(1-\alpha^{2})$ and the dispersion index of the process

$$VMR_{X}=\frac{\sigma^{2}_{X}}{\mu_{X}}=\frac{VMR_{\varepsilon}+\alpha}{1+\alpha}=\frac{1}{1+\alpha}\left(\alpha+\frac{1-\beta}{(1-\beta-b\beta/c)^{2}}\right)$$

where $VMR_{\varepsilon}=\sigma_{\varepsilon}^{2}/\mu_{\varepsilon}$ is the innovation index of dispersion. It follows that there is under– or over–dispersion in the marginal distribution, $VMR_{X}<1$ and $VMR_{X}>1$, if and only if there is under- or over–dispersion in the innovation, $VMR_{\varepsilon}<1$ or $VMR_{\varepsilon}<1$ respectively.

The autocorrelation function is

$$\gamma_{k}=\mathbb{C}ov(X_{t},X_{t-k})=\mathbb{E}(X_{t}X_{t-k})-\mu_{X}^{2}=\alpha^{k}\sigma^{2}_{X}$$

as in the INAR(1) process [e.g., see Al-Osh & Alzaid, 1987].

The GLK–INAR(1) process can be extended to account for various sources of model instability such as structural breaks, regimes and outliers by introducing a time–varying parameter setting [see for example Malyshkina et al., 2009]. A parsimonious approach is to assume a finite set of regimes $k=1,\ldots,K$ corresponding to different parameter configurations, i.e.

$$X_{t}=\alpha(S_{t})\circ X_{t-1}+\varepsilon_{t},\qquad t\in\mathbb{Z},$$

(7)

with $\varepsilon_{t}|S_{t}\sim\mathcal{GLK}(a(S_{t}),b(S_{t}),c(S_{t}),\beta(S_{t}))$, where the thinning coefficient and the $GLK$ parameters of the error term $\psi(S_{t})=(\alpha(S_{t}),$ $a(S_{t}),$ $b(S_{t}),c(S_{t}),$ $\beta(S_{t}))^{\prime}$ are time–varying $\psi(S_{t})=\sum_{k=1}^{K}\mathbb{I}(S_{t}=k)\psi_{k}$, where $\psi_{k}=(a(k),b(k),$ $c(k),\beta(k))^{\prime}$. The $S_{t}\in\{1,\ldots,K\}$ for $t\in\mathbb{Z}$ denotes a hidden Markov–chain process with transition probabilities $\mathbb{P}(S_{t}=j|S_{t-1}=i)=\pi_{ij}$ for $i,j\in\{1,\ldots,K\}$. From now on, this extension is denoted with MS–GLK–INAR(1).

A special case, which is relevant for a common issue in count data series, is the large proportion of zeros [e.g., Maiti et al., 2015]. The excess of zero, which leads to over–dispersion, can be handled by assuming a zero–inflated GLK–INAR(1). This can be defined by assuming that in one of the regimes, e.g. $S_{t}=1$, there is complete thinning $\alpha_{1}\to 0$, and the error distribution is a Dirac centred at zero. Some alternatives where the GLK collapses to a Dirac include the (Negative) Binomial distribution with parameters $c_{1}=1$, $b_{1}=-1$ ($b_{1}=0$), $a_{1}=1$ and $\beta_{1}\to 0$.

The transition probabilities provide information on the persistence of these events. If the probability is independent on past information, i.e. $\mathbb{P}(S_{t}=j|S_{t-1}=i)=\pi_{j}$ for all $i,j\in\{1,2\}$, then the zero–inflated regime is transitory. If the zero–inflated regime is persistent, then the duration of the regime is captured by a large $\pi_{11}$. Other regimes ($S_{t}\neq 1$) with low mean and/or large variance can also generate zeroes. This is convenient in some applications, such as in epidemiology where zeroes from $S_{t}\neq 1$ can be interpreted as under–reported cases of a particular disease [e.g., Douwes-Schultz & Schmidt, 2022].

The GLK–INAR can be extended to include more general auto–correlation structures and to the multivariate setting. The process can be modified to allow for multiple lags building on the specification strategy used in Neal & Subba Rao [2007]. In particular, a GLK integer–valued ARMA of order $p$ and $q$, i.e. GLK–INARMA$(p,q)$ can be specified using independent thinning operators.

Compared to GLK–INAR(1), in the GLK–INARMA$(p,q)$, a further restriction on the autoregressive parameters is required for stationarity, although a weaker condition can be used. For alternatives specification strategies such as combined INAR (CINAR) see for example McKenzie [2003], Weiß [2008].

For the case of random integer vectors, a multivariate GLK–INAR(1) (GLK–MINAR(1)) can be used by introducing a thinning matrix operator.

The independence between the innovation terms of the different equations allows us to estimate them separately. This assumption can be relaxed by adding common GLK errors to the equations, or under some special cases of the GLK, a joint distribution can be introduced, such as the bivariate Katz’s or Poisson distribution [Pedeli & Karlis, 2011, Diafouka et al., 2022].

With the construction in Eq. 1, the constraint $\sum_{x\geq 0}p_{x}=1$ is guaranteed by the condition $H(1)=1$. Nevertheless, some constraints on the parameters are needed to have all the $p_{x}:>0$. Three different cases are discussed below (details are given in Appendix B in the Supplementary Material.).

• •

  For parameter values $a>0,b\geq 0,c>0$ the pmf are positive. Moreover, for $a/c,b/c\in\mathbb{N}$ the extended binomial coefficient $((a+bx)/c+1)_{x\uparrow}/x!$ coincides with the standard binomial coefficient ${\frac{a+bx}{c}+x\choose x}$ [Consul & Famoye, 2006, p. 8].

• •

  For $-c<b<0$, $a/c,b/c\in\mathbb{N}$ and $(c-a)/(c+b)\leq(a+c)/|b|$, the pmf are positive for $x<x^{*}=(a+c)/|b|$, while $p_{x}=0$ for $x\geq x^{*}$.

• •

  If $-c<b<0$ but the additional constraints of the previous point are not satisfied, the terms appearing in the product $({(a+bx)}/{c}+1)_{x\uparrow}$ change sign, and there is no guarantee that the result is positive. Indeed, for all the $x$ such that $x>\max\{(a+1)/|b|,(c-a)/(c+b),2\}$ one has $(a+bx)/c+1<0$ and $(a+bx)/c+x-1>0$ and hence there is an integer $q=q_{x}$ such that $(a+bx)/c+m<0$ for $1\leq m\leq q$ and $(a+bx)/c+m>0$ for $q+1\leq m\leq k-1$. Hence

  	$\displaystyle\frac{1}{\frac{a+bx}{c}+x}\Big{(}\frac{(a+bx)}{c}+1\Big{)}_{x\uparrow}=(-1)^{q}\prod_{m=1}^{q}\Big{|}\frac{(a+bx)}{c}$
  	$\displaystyle+m\Big{|}\prod_{m=q+1}^{k-1}\Big{|}\frac{(a+bx)}{c}+m\Big{|}$

  which is negative whenever $q_{x}$ is odd. For example take $a=10$, $b=-1$ and $c=2$, for $x=20$ one has $q_{20}=5$, which shows that $p_{20}<0$ which clearly is impossible.

In a Bayesian framework, the parameter constraints can be easily included in the inference process through a suitable choice of the prior distributions. We assume:

	$\displaystyle\alpha\sim\mathcal{B}e(\kappa_{\alpha},\tau_{\alpha}),\quad a\sim\mathcal{G}a(\kappa_{a},\tau_{a}),\quad b\sim\mathcal{G}a(\kappa_{b},\tau_{b}),$
	$\displaystyle c\sim\mathcal{G}a(\kappa_{c},\tau_{c}),\quad\beta\sim\mathcal{B}e(\kappa_{\beta},\tau_{\beta})$

where $\mathcal{B}e(\kappa,\tau)$ is the beta distribution with shape parameters $\kappa$ and $\tau$ and $\mathcal{G}a(\kappa,\tau)$ the gamma distribution with shape and scale parameters $\kappa$ and $\tau$, respectively. In the empirical applications we assume a non-informative hyper-parameter setting for $\alpha$ and $\beta$, that is $\kappa_{\alpha}=\tau_{\alpha}=\kappa_{\beta}=\tau_{\beta}=1$ and an informative prior for $a,b$ and $c$ with $\kappa_{a}=\tau_{a}=1$, $\kappa_{b}=\kappa_{c}=2$ and $\tau_{b}=\tau_{c}=1/2$.

In the case of Markov–switching specification of the GLK–INAR(1), the same prior is assumed for the regime–specific parameters $h(\psi_{k})=\mathcal{B}e(\kappa_{\alpha},\tau_{\alpha})$ $\mathcal{G}a(\kappa_{a},\tau_{a})$ $\mathcal{G}a(\kappa_{b},\tau_{b})$ $\mathcal{G}a(\kappa_{c},\tau_{c})\mathcal{B}e(\kappa_{\beta},\tau_{\beta})$ for $k=1,\ldots,K$. For the transition probabilities of the allocation variable, we assume a symmetric Dirichlet prior for each row of the transition matrix, i.e. $\pi_{i\cdot}\sim\mathcal{D}(1/K,\ldots,1/K)$ with concentration parameter $1/K$, where $\pi_{i\cdot}=(\pi_{i1},\ldots,\pi_{iK})$ for $i=1,\ldots,K$.

Let $x_{1},\ldots,x_{T}$ be a sequence of observations for the GLK–INAR(1) process, then the joint posterior distribution is given by

$$f(\psi|x_{1},\ldots,x_{T})\propto f(\psi)\prod_{t=1}^{T}\prod_{i=0}^{\infty}\prod_{j=0}^{\infty}P_{ij}(\psi)^{\mathbb{I}(x_{t}-j)\mathbb{I}(x_{t-1}-i)}$$

where $\psi=(\alpha,a,b,c,\beta)$ is the parameter vector $f(\psi)$ the joint prior and

$$P_{ij}(\psi)=\sum_{k=0}^{\min(i,j)}d_{ijk}\binom{\frac{a+b(j-k)}{c}+j-k}{j-k}\alpha^{k}(1-\alpha)^{i-k}\beta^{j-k}(1-\beta)^{\frac{a+b(j-k)}{c}}$$

where $d_{ijk}=\binom{i}{k}((a/c)/((a+bx)/(c)+j-k)$.

Following the discussion above in this section, if the parameter constraint $c>0$ is not imposed, the coefficients of the Lagrangian expansion can be negative. In this case, a truncated GLK can be used, similarly to what is proposed in McCabe & Skeels [2020] for the Katz distribution, and the inference procedure can be easily extended to include this type of distribution. The truncation can be imposed by using the following recursion for the transition probability:

$$p_{i}(\psi)=p_{0}\prod_{j=0}^{i-1}\max\left\{0,\frac{U(\psi)+V(\psi)j}{a+j}\right\}$$

where $U(\psi)=a\beta/c$, $V(\psi)=U(b+c)/(a+b)$ and

$$p_{0}=\left(1+\sum_{j=1}^{\infty}\prod_{k=0}^{j-1}\max\left\{0,\frac{U(\psi)+V(\psi)j}{a+j}\right\}\right)^{-1}.$$

The probability $p_{i}$ becomes null for $i>j$ if $U(\psi)+V(\psi)j<0$ at $j$.

Since the joint posterior is not tractable, we follow a Markov Chain Monte Carlo (MCMC) framework for posterior approximation. See Robert & Casella [2013] for an introduction to MCMC methods. We overcome the difficulties in tuning the parameters of the MCMC procedure by applying the Adaptive MCMC sampler (AMCMC) proposed in Andrieu & Thoms [2008]. Following a standard procedure, the following reparametrization is considered to impose constraints on the parameters of the GLK–INAR(1). Let $\eta=(\eta_{1},\ldots,\eta_{5})$ be the 5-dimensional parameter vector obtained by the transformation $\eta=\varphi(\psi)$ with $\eta_{1}=\log(\psi_{1}/(1-\psi_{1}))$, $\eta_{2}=\log(\psi_{2})$, $\eta_{3}=\log(\psi_{3})$, $\eta_{4}=\log(\psi_{4})$, and $\eta_{5}=\log(\psi_{5}/(1-\psi_{5}))$ and let $f(\eta|x_{1},\ldots,x_{T})=f(\varphi^{-1}(\eta)|x_{1},\ldots,x_{T})J(\eta)$ be the posterior of $\eta$, with $J(\eta)=\psi_{1}\psi_{2}\psi_{3}\psi_{4}\psi_{5}(1-\psi_{1})(1-\psi_{5})$ the Jacobian of the transformation $\varphi$ given above. Given the adaptation parameters $\mu^{j}$ and $\Sigma^{(j)}$, at the $j$-th iterations, the AMCMC consists of the following three steps. First, a candidate $\eta^{\ast}$ is generated from the random walk proposal: $\eta^{\ast}=\eta^{(j-1)}+\lambda^{(j)}w^{(j)},\quad w^{(j)}\sim\mathcal{N}_{q}(\mathbf{0},\Sigma^{(j)})$. Second, the candidate is accepted with probability $\rho^{(j)}=\rho(\eta^{(j-1)},\eta^{\ast})$, where

$$\rho(\eta^{(j-1)},\eta^{\ast})=\min\left(1,\frac{f(\varphi^{-1}(\eta^{\ast})|x_{1},\ldots,x_{T})J(\eta^{\ast})}{f(\varphi^{-1}(\eta^{(j-1)})|x_{1},\ldots,x_{T})J(\eta^{(j-1)})}\right)$$

and third, the adaptive parameters are updated as follows:

	$\displaystyle\mu^{(j+1)}$	$\displaystyle=$	$\displaystyle\mu^{(j)}+\gamma^{(j)}(\mu^{(j)}-\eta^{(j)})$
	$\displaystyle\Sigma^{(j+1)}$	$\displaystyle=$	$\displaystyle\Sigma^{(j)}+\gamma^{(j)}((\mu^{(j)}-\eta^{(j)})(\mu^{(j)}-\eta^{(j)})^{\prime}-\Sigma^{(j)})$
	$\displaystyle\log\lambda^{(j+1)}$	$\displaystyle=$	$\displaystyle\log\lambda^{(j)}+\gamma^{(j)}(\rho^{(j)}-\rho^{\ast}),$

where $\rho^{*}$ is the target acceptance probability and $\gamma^{(j)}=j^{-a}$, $a>0$ is the adaptive scale [Andrieu & Thoms, 2008, , Algorithm 4]. Following the suggestions in Roberts et al. [1997] we set $\rho^{*}=0.44$.

The latent allocation variables in the Markov–switching specification of the GLK–INAR(1) are sampled the Forward–Filtering Backward sampling procedure (FFBS). The prediction and filtering probabilities are given by

		$\displaystyle\mathbb{P}(S_{t}=k|\mathcal{X}_{t-1})=\sum_{\ell=1}^{K}\pi_{\ell k}\mathbb{P}(S_{t-1}=\ell|\mathcal{X}_{t-1})$
		$\displaystyle\mathbb{P}(S_{t}=k|\mathcal{X}_{t})\propto f(x_{t}|\psi_{k},x_{t-1},S_{t}=k)\mathbb{P}(S_{t}=k|\mathcal{X}_{t-1}),$

where $\mathcal{X}_{t-1}=(x_{1},\ldots,x_{t-1})^{\prime}$, $f(x_{t}|\psi_{k},x_{t-1},S_{t}=k)=$ $\prod_{i=0}^{\infty}$ $\prod_{j=0}^{\infty}$ $P_{ij}(\psi_{k})^{\mathbb{I}(x_{t}-j)\mathbb{I}(x_{t-1}-i)}$ for $k,\ell\in\{1,\ldots,K\}$. Notice that the conditioning on the parameters $\psi$ is included only in the likelihood but not in the probabilities to simplify the notation. The filtered probabilities can be smoothed by considering all the information available, i.e.

$$\mathbb{P}\left(S_{1:T}|\mathcal{X}_{T}\right)=\mathbb{P}\left(S_{T}|\mathcal{X}_{T}\right)\prod_{t=1}^{T-1}\mathbb{P}\left(S_{t}|S_{t+1},\mathcal{X}_{t}\right)$$

where $\mathbb{P}\left(S_{t}|S_{t+1},\mathcal{X}_{t}\right)\propto\pi_{S_{t}S_{t+1}}\mathbb{P}\left(S_{t}|\mathcal{X}_{t}\right)$ and $S_{1:T}=(S_{1},\ldots,S_{T})^{\prime}$. The allocation variables are sampled directly from these smoothed probabilities.

The conditional posterior distribution of the transition probabilities of the Markov chain $S_{t}$ is conditionally conjugate and can be sampled directly from $\pi_{\cdot k}|S_{1:T}\sim\mathcal{D}(d_{1},\ldots,d_{K}),$ where $d_{k}=1/K+\sum_{t=1}^{T}\mathbb{I}_{k}(S_{t})$ for $k=1,\ldots,K$.

For the possible extensions of the GLK–INAR, such as the GLK–INRMA(p,q) and the GLK–MINAR(1), data augmentation techniques can be used to improve the efficiency of the MCMC [Neal & Subba Rao, 2007, Marques et al., 2022]. For instance, in the case of the GLK–INRMA(p,q), conditional conjugacy of the thinning parameters can be obtained by assuming each autoregressive (moving average) component is a latent variable following a binomial distribution. In the case of GLK–MINAR(1), a similar strategy can be followed, see for instance Soyer & Zhang [2022].

We illustrate the Bayesian procedure’s effectiveness in recovering the parameters’ true value and the MCMC procedure’s efficiency through some simulation experiments. We test the algorithm’s efficiency in two different settings, commonly found in the data: low persistence and high persistence (see trajectories in Fig. 2). The true values of the parameters are: $\alpha=0.3$, $a=5.3239$, $b=0.0592$, $c=0.6$, $\beta=0.5917$ in the low persistence setting, and $\alpha=0.7$, $a=5.3239$,$b=0.0592$, $c=0.6$, $\beta=0.5917$ in the high persistence setting. For each setting, we run the Gibbs sampler for 50,000 iterations on each dataset, discard the first 10,000 draws to remove dependence on initial conditions, and apply a thinning procedure with a factor of 10 to reduce the dependence between consecutive draws.

For illustrative purposes, in Figure C.1 in the Supplementary Material we show the MCMC posterior approximation for the parameter $\alpha$ (first row), the unconditional mean of the process (second row), and the marginal likelihood (last row), in one of our experiments for the high- and low-persistence settings. Each plot represents the true value (solid black line) and the Bayesian estimates. Posterior estimated are approximated by using 4,000 MCMC samples after thinning and burn-in removal (dashed red line). Figures C.2-C.3 and C.4-C.5 in the Supplementary Material exhibit 10,000 MCMC posterior draws and the MCMC approximation of the posterior distribution for all the parameters, in the high- and low-persistence settings.

In our experiments, the acceptance rate is in the range of 40%-53% for both parameter settings (see Figure C.6 in the Supplementary Material). Table C.1 in the Supplementary Material shows, for all the parameters the autocorrelation function (ACF), effective sample size (ESS), inefficiency factor (INEFF) and Geweke’s convergence diagnostic (CD) before (BT subscript) and after thinning (AT subscript). The numerical standard errors are evaluated using the nse package [Geyer, 1992, Ardia & Bluteau, 2017, Ardia et al., 2018].

The thinning procedure is effective in reducing the autocorrelation levels and in increasing the ESS. The p-values of the CD statistics indicate that the null hypothesis that two sub-samples of the MCMC draws have the same distribution is always accepted. The efficiency of the MCMC after the thinning procedure is generally improved. After thinning, on average, the inefficiency measures (5.83), the p-values of the CD statistics (0.36) and the NSE (0.02) achieved the values recommended in the literature [e.g., see Roberts et al., 1997].

It is important to underline that the persistence parameter estimation and the forecast are highly sensitive to the innovation distributional assumption. An illustration is presented in the left plot of Figure 3, where the data generating process corresponds to a GLK–INAR(1) with large overdispersion (VRM=8.6). The standard model for count data is the Poisson INAR(1) model (PINAR(1)), which cannot capture overdispersion. This misspecified model entails an underestimation of the persistence parameter (medium gray histogram). The NBINAR(1) captures the overdispersion and provides reliable persistence estimates (light gray) comparable with the one of GLK–INAR (dark gray). Nevertheless, in the case of underdispersion (VRM=0.4, right plot of Figure 3), both NBINAR(1) and PINAR(1) return an estimation bias in the persistence parameter, while the INAR–GLK gives a good approximation of the true persistence. In summary, the INAR–GLK(1) model nests standard models, such as Generalized Poisson and Negative Binomial INARs, and allows for different degrees of underdispersion and overdispersion. Hence, it can be used without preliminary testing of the dispersion features of the series.

Similarly, to exemplify the effectiveness and efficiency of the estimation procedure in different scenarios we considered: i) high and low persistence regimes with the same parameter configuration of the settings presented before and ii) a large mean regime and an zero–inflated regime where $\alpha\to 0$, $a=1$, $b\to 0$, $c=1$, $\beta\to 0$. The simulated trajectories are shown in Figure C.7 (C.8) in the Supplementary Material together with the estimated allocations of the regimes, represented by the shaded areas, with an accuracy of 97% (100%) for the two regimes (zero–inflated) scenario. Moreover, the parameters are successfully retrieved, see in Figures C, C and C in the Supplementary Material. Notice that the zero–inflated parameters are not estimated but set by default to approximate the Dirac distribution.

In conclusion, the Gibbs sampler is computationally efficient and can retrieve the true parameter values of the MS–GLK–INAR in different settings, including the single–regime and the zero–inflated specifications. The MCMC for the GLK–INAR takes 0.5 minutes for a sample size of $T=260$ observations and for 30,000 MCMC iterations. This is comparable with the Negative Binomial INAR (0.4 minutes). The method is scalable and can be applied to datasets with thousands of observations. For larger-size datasets, the theoretical moment of the process can be used to devise alternative estimation procedures, such as the method of moments. The moments of the distribution are provided in closed form in Proposition 5 in the Supplementary Material.

We used Google Trends data to measure the changes in public concern about climate change. Google Trend represents a source of big data [Choi & Varian, 2012, Scott & Varian, 2014] which have been used in many studies, for example, Anderberg et al. [2021] studied domestic violence during covid-19, Yang et al. [2021] studied influenza trends, Schiavoni et al. [2021] and Yi et al. [2021] presented applications to unemployment and Yu et al. [2019] studied oil consumption. In this study, we follow Lineman et al. [2015] and use Google search volumes as a proxy for public concern about “Climate Change” (CC) and “Global Warming” (GW). The search volume is the traffic for the specific combination of keywords relative to all queries submitted in Google Search in the world or a given region over a defined period. The indicator ranges from 0 to 100, with 100 corresponding to the largest relative search volume during the period of interest. The search volume is sampled weekly from 4th December 2016 until 21st November 2021. We analysed the dynamics at the global and country level. Countries with an excess of zeros above 95% in the search volume series have been excluded. The final dataset includes 65 countries of the about 200 countries provided by Google Trends. For illustration purposes, we report in the top plots of Fig. 4 the series of the world volume. The CC global volume exhibits overdispersion with $\widehat{VMR}=102/27.33=3.73$, skewness and kurtosis $\widehat{S}=2.09$ and $\widehat{K}=13.47$, respectively. The GW global volume has over–dispersion $\widehat{VMR}=170.42/48.56=3.51$, skewness $\widehat{S}=0.27$ and kurtosis $\widehat{K}=3.22$ (see also the histograms in the bottom plots). The country-specific indexes exhibit different levels of persistence and over–dispersion.