Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
63,245 characters · 17 sections · 93 citation commands
First--order integer--valued autoregressive processes with Generalized Katz innovations
\doublespacing
In the recent years there has been a large interest in discrete--time integer--valued models, also due to increased availability of count data in very diverse fields including finance LieMolPoh2006,AknoucheAlmohaimeedDimitrakopoulos2021, economics FreCab2004,Berri2020, social sciences PedKar2011, sports ShaMoy2016, image processing Afrifa2022 and oceanography CunVasBou2018. Among the modelling approaches, integer--valued autoregressive processes (INAR), introduced independently by AlOAlz1987 and McK1985, become \textcolor{black}{very} popular. The stochastic construction of the INAR relies on the binomial thinning operator and the properties of the model on the discrete self--decomposability of the stationary distribution of the process Ste79. See ScoWeiGou2015 for a review.
The original INAR model has been studied further in AlOAlz1987 and extended in different directions. McK1986 introduced an INAR model with negative--Binomial and geometric marginal distributions, DuLi1991 extended the INAR(1) model of AlOAlz1987 to the higher order INAR$(p)$. AlOAly1992 introduced a negative--binomial INAR with a new iterated thinning operator. Other extensions of the INAR process have been made to include a seasonal structure in the model Bouetal16. INAR models with values in the set of signed integers have been propose firstly by KimPar2008 and generalised by AlzOma2014 and AndKar2014. Fre2010 proposed a true integer--valued autoregressive model (TINAR(1)). More flexible INAR models have been introduced by assuming more flexible distributions for the innovations terms. AlzAlO1993 propose integer--valued ARMA process with Generalized Poisson marginals and KimLee2017 introduced INAR with Katz innovations.
This paper introduces a general class of INARs with Generalized Lagrangian Katz innovations. The Lagrangian Katz family is a flexible distribution and naturally arises as first crossing probabilities, which is a common problem in actuarial mathematics, e.g. claim number distribution in cascading processes or ruin probability in discrete--time risk models ConFam2006. It has been extended further by jan98 and jan99, which introduced the four--parameter generalized P\'{o}lya--Eggenberger (GPED) distributions of the first and second kind. jan98 showed that both families contain the Lagrangian Katz distribution as a special case. \textcolor{black}{We consider the four-parameters GPED of the first kind, also known as Generalized Lagrangian Katz (GLK). The resulting process family provides a flexible modelling framework for count data, allowing for under and over--dispersion, asymmetry, and excess of kurtosis and includes standard INAR models such as Generalized Poisson and Negative Binomial as special cases.} \textcolor{black}{Further extensions are provided, such as the Markov--Switching and the zero--inflated GLK--INARs, to account for different sources of model instability and excess of zeros. }
Various approaches to inference have traditionally been presented for count data models, such as the conditional likelihood approach, generalized method of moments and Yule--Walker approach. See WeiKim2013 for a review. Despite the popularity gained in recent years by Bayesian methods, the applications to count data models are still limited CabMar05, NeSu07, Droetal16, ShaZha18, Garetal20. Thus, we provide a Bayesian inference procedure for our model and illustrate the procedure's efficiency on a synthetic dataset. The Bayesian approach to inference entirely considers parameter uncertainty in the prior knowledge about a random process. It allows for imposing parameter restrictions by specifying the prior distribution Chen2016. The posterior distribution of the parameters quantifies uncertainty in the estimation Chen17, which can be included in the prediction. The inference from the Bayesian perspective may result in richer inferences in the case of small samples Garayetal20 and extra--sample information and in robust inference in the presence of outliers fried2015retrospective. Finally, model selection for both nested and non--nested models can be easily carried out.
\textcolor{black}{We illustrate the model's flexibility with an application to an original Google Trend dataset of 130 time--series measuring the public concern about climate change in different countries. The contrasting features of the series, such as excess of zeros, outliers, and regimes, are common in count data and provide a challenging and diversified ground for illustrating the robustness and flexibility of the GLK--INAR model}. Assessing public awareness and knowledge of a specific topic and understanding the dynamics of social consciousness allows for designing more effective public policies. For this reason, researchers measured and studied the level of awareness about the effects of climate change in different sectors of society such as households Frondel2017, winegrowers bat2009, farmers Fahad2018, mountain peoples Ullah2018. Most of these studies rely on surveys conducted in a specific geographical area and sector of society, with a few exceptions. For example, Zieg2017 proposed a cross--country analysis of climate change beliefs and attitudes. lineman2015talking provided a broader and global perspective by exploiting the potentiality of big data provided by Google Trend. This extended climate change perception literature along two lines. First, we consider a multi--country dataset, including country--specific measures to capture worldwide heterogeneity in public awareness. Moreover, we offer a model--based approach and an inference procedure to analyze these measures.
The paper is organized as follows. Section 2 introduces the GLK family and INAR process with some extensions such as the Markov--Switching GLK--INAR. Section 3 proposes a Bayesian inference procedure and provides some simulation results. Section 4 provides some illustrations on a multi--country Google Trend dataset related to climate change. Section 5 concludes.
The probability mass function (pmf), $P(X=x)=p_x$, of the Generalized Lagrangian Katz (GLK) is
$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 ConFam2006. 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 (ref) and Appendix (ref) in the Supplementary. GLK distributions have probability generating function (pgf)
which satisfies:
or alternatively
see jan98.
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 probability mass function of the GLK for different parameter settings is given in Fig. (ref). 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 (ref) 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 (ref) 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
where $B_{i}(\alpha)$ are iid Bernoulli r.v.s with success probability $P(B_{i}(\alpha)=1)=\alpha$.
Figure (ref) 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. (ref), 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 (ref), 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 KimLee2017 as special case. The GLK--INAR(1) naturally nests the Poisson INAR(1) of AlOAlz1987, the Negative Binomial INAR(1) of AlOAly1992 (NBINAR(1)), and the Generalized Poisson INAR(1) of AlzAlO1993.
As for any INAR process, the GLK--INAR(1) has the following representation
and its conditional pgf can be written as
where $H(u)$ is defined in Eq. (ref) or in Eq. (ref). Starting from the general results on INAR processes given in AlOAlz1988, one easily obtains explicit expressions for the conditional mean and variance of the GLK--INAR(1):
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
where $p_x$ is the pmf given in Eq. (ref).
In the next Proposition, we summarize some of the asymptotic properties of a GLK--INAR(1). These properties follow from general results in AlOAlz1988 and Schweer2014.
Since the GLK distribution satisfies the convolution property jan98, 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
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
as in the INAR(1) process AlOAlz1987. {\color{black}
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 malyshkina2009markov. A parsimonious approach is to assume a finite set of regimes $k=1,\ldots, K$ corresponding to different parameter configurations, i.e.
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))'$ 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))'$. 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 maiti2015time. 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\to0$.
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 douwes2022zero.
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 NeSu07. 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 mckenzie2003ch,weiss2008combined.
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 PedKar2011,diafouka2022bivariate. }
With the construction in Eq. (ref), 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 (ref) in the Supplementary Material.).
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:
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$.
{\color{black} 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
where $\psi=(\alpha,a,b,c,\beta)$ is the parameter vector $f(\psi)$ the joint prior and
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 mccabe2020distributions 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:
where $U(\psi)=a\beta/c$, $V(\psi)=U(b+c)/(a+b)$ and
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 RobCas2013 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 Andr08. 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
and third, the adaptive parameters are updated as follows:
where $\rho^{*}$ is the target acceptance probability and $\gamma^{(j)}=j^{-a}$, $a> 0$ is the adaptive scale Andr08. Following the suggestions in Robetal1997 we set $\rho^{*} = 0.44$.
{\color{black} 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
where $\mathcal{X}_{t-1}=(x_1,\ldots,x_{t-1})'$, $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.
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)'$. 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$. }
{\color{black} 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 NeSu07,c2022bayesian. 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 soyer2022bayesian. }
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. (ref)). 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 (ref) 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 (ref)-(ref) and (ref)-(ref) 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 (ref) in the Supplementary Material). Table (ref) 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 Geyer92, ArdBlu17, Ardiaetal18.
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 Robetal1997.
{\color{black} 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 (ref), 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 (ref)), 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 (ref) ((ref)) 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 (ref), (ref) and (ref) 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 (ref) 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 choi2012,scott2014 which have been used in many studies, for example, AndHel21 studied domestic violence during covid-19, Yaetal21 studied influenza trends, Schetal2020 and Dinetal21 presented applications to unemployment and Yuetal19 studied oil consumption. In this study, we follow lineman2015talking 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. (ref) 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.
The posterior distribution of the autoregressive coefficient is given in Fig. (ref). The coefficient estimate and posterior credible interval (in parenthesis) are $\widehat{\alpha}=0.56$ $(0.50,0.62)$ and $\widehat{\alpha}=0.62$ $(0.56,0.67)$ for the GW and the CC dataset, respectively (see also the approximation to the posterior distribution of the parameters in Figures (ref) and (ref) in the Supplementary Material). This result indicates that the public concern about climate risk is persistent over time worldwide at an aggregate level. The estimated parameter of the innovation process and their 0.95% credible intervals (in parenthesis) are $\widehat{a}=3.53$ $(1.56,6.08)$, $\widehat{b}=0.04$ $(0.01, 0.11)$, $\widehat{c}=0.21$ $(0.05,0.47)$ and $\widehat{\beta}=0.48$ $(0.20,0.65)$ for the GW dataset and $\widehat{a}=3.26$ $(1.44,5.72)$, $\widehat{b}=0.12$ $(0.021,0.310)$, $\widehat{c}=0.26$ $(0.032,0.726)$ and $\widehat{\beta}=0.35$ $(0.067,0.623)$ for the CC one.
The results indicate a deviation from the Negative Binomial model. Thus we apply the DIC criterion $DIC=-4\mathbb{E}(\log f(X|\psi)|y)+2\log f(X|\widehat{\psi})$ to compare GLK--INAR(1) and NB--INAR(1). The $DIC$ is computed following Spiegel2002:
where $f(X|\psi)$ is the likelihood of the model, $\psi^{(j)}$ $j=1,\ldots,N$ the MCMC draws after thinning and burn-in sample removal, and $\widehat{\psi}$ is the parameter estimate. The DICs for the GLK (NB) INARs fitted on the aggregate CC and GW series are $1.6743\cdot 10^3$ $(1.6862\cdot 10^3)$ and $1.8735\cdot 10^3$ $(1.8834\cdot 10^3)$, respectively.
{\color{black} Given the high kurtosis levels and the multi--modality in the empirical distribution of both series (see Figure (ref)), the Markov Switching INAR is used to deal with outliers and parameter instability. We use DIC and RMSE to select the number of regimes and the model (see Table (ref) in the Supplementary Material). We find that GLK--INAR with two or three regimes present the best fit in--sample and out--sample for both the CC and GW datasets. The results with three--regimes are presented in Figure (ref). The three regimes identify different persistence levels: high (right plot), medium (middle) and low (left). Some of the regimes also have different unconditional mean levels (see Figure (ref) left plots). In terms of one--step--ahead forecasting, in both datasets, the model can reproduce the upward trend at the end of the sample and effectively cover the true values within their 90% credible intervals (see Figure (ref) right plots).}
{\color{black}
We run the analysis at a disaggregate level. The results are given in Figures (ref)-(ref) and Tables (ref)-(ref) in the Supplementary Material. Figure (ref) provides evidence of an inverse relationship between estimated persistence $\widehat{\alpha}$ and dispersion $\widehat{VMR}$ cross countries (reference lines in the left plot). There is evidence of this inverse relationship in both the CC (blue dots) and GW (red dots) datasets. The plot on the right indicates an inverse (direct) relationship between the estimated unconditional mean $\widehat{\mu_\varepsilon}/(1-\widehat{\alpha})$ and the dispersion index $\widehat{VMR}$ for the GW (CC). In the same picture, we indicate the parameter estimates for the world volume of searches (stars). }
The terms “Climate Change" and “Global Warming" are used interchangeably. Nevertheless, they describe different phenomena and can be used to determine the public's level of understanding about these two parallel concepts lineman2015talking. We investigate the relationships in the search volumes through the lens of our GLK--INAR(1) model. The left plot in Fig. (ref) shows the unconditional mean of the search volumes for the two concepts in all countries (dots). In public attention, the two concepts are connected in the long run. We find a positive association for both countries with large (percentage of zeros $<$ 21%) and low search volumes (percentage of zeros $>$ 21%). There is an asymmetric effect in the overdispersion (right plot), and in all countries, the GW search volume has a larger VMR than the CC volume. This can be explained by the larger variability induced by the changes in the use of the GW term in official communications.
Comparing the coefficients across the rows of Tables (ref)-(ref) in the Supplementary Material, we find evidence of two types of series, one with high persistence and the other with low persistence. Moreover, for each country, the level of persistence is similar across the two datasets (compare columns of Tables (ref)-(ref) in the Supplementary Material).
Tables (ref)-(ref) in the Supplementary Material report the marginal likelihood of the GLK--INAR(1) and Lagrangian Katz INAR(1) in columns GLK and LK, respectively. We find evidence of a better fitting of the GLK--INAR(1) for some countries and variables, e.g. CC searches in India and CC and GW searches in South Africa. To get further insights into the results, we study the relationship between the dynamic and dispersion properties of the series and the actual level of climate risk of the countries. We consider the Global Climate Risk Index (CRI), which ranks countries and regions following the impacts of extreme weather events (such as storms, hurricanes, floods, heatwaves, etc.). The lower the index value, the larger the climate risk is. Following the values of the CRI for 2021, based on the events recorded from 2000 to 2019, our dataset includes some of the countries most exposed to climate risk, such as Japan, Philippines, Germany, South Africa, India, Sri Lanka and Canada eckstein2021.
The left plot in Fig. (ref) shows the unconditional mean against the CRI. There is evidence of a positive relationship between the public interest in climate-related topics and the actual level of climatic risk. The lower the CRI level, the larger the Google search volumes are (see dashed lines). For example, India has a high risk (CRI equal to 7) and a very high long-run level of public attention.
The right plot reports the coefficient of variation against the CRI for all countries in the “Climate Change" (blue) and “Global Warming" (red) datasets. The dashed lines represent linear regressions estimated on the data. There is evidence of a negative relationship between the dispersion of public concern and climatic risk; in countries with more significant risk levels, the Google search volumes are less over--dispersed.
{\color{black} To deal with the excess of zeros, which are very frequent in more than 62% of the series, we apply the MS--GLK--INAR(1) with two states, where the first state represents a prolonged absence of searches on Google and the second a persistent search activity. The MS--GLK--INAR(1) performs better than MS--NBINAR(1) in 119 out of the 130 CC and GW time series following the DIC (see Table (ref) and (ref) in the Supplementary Material). Accounting for the excess of zeros allows for improving the estimation of the persistence and provides an estimate of the probability $\hat{\pi}_{11}$ to stay in an inactive search regime. The findings on the persistence parameter discussed in this section for the GLK--INAR(1) are confirmed by the MS--GLK--INAR. Furthermore, there is evidence of a positive relationship between the probability $\hat{\pi}_{11}$ and the CRI, consistent with the results on the Google search persistence.}
\textcolor{black}{A novel integer--valued autoregressive process is proposed with Generalized Lagrangian Katz innovations (GLK--INAR). Theoretical properties of the model, such as stationarity, moments, and semi--self--decomposability, are provided. To deal with parameter instability and excess of zeroes, we also propose a Markov--Switching GLK--INAR. A Bayesian approach to inference and an efficient Gibbs sampling procedure have been proposed, which naturally account for uncertainty when forecasting. The modelling framework is applied to a Google Trend dataset measuring the public concern about climate change in 65 countries. The greater flexibility of the GLK--INAR allows for a superior fitting compared to the standard INAR models and a better comprehension of the heterogeneity in public perception. More specifically, new evidence is provided about the long-run level of public attention, its persistence and dispersion in countries with low and high levels of climate risk. The Markov-switching GLK-INAR identified regimes with the absence of searches and changes in the dynamic features of the series.}
The authors acknowledge support from: the MUR– PRIN project ‘Discrete random structures for Bayesian learning and prediction’ under g.a. n. 2022CLTYP4 and the Next Generation EU -- ‘GRINS– Growing Resilient, INclusive and Sustainable’ project (PE0000018), National Recovery and Resilience Plan (NRRP)-- PE9 -- Mission 4, C2, Intervention 1.3. The views and opinions expressed are only those of the authors and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
This supplement consists of four Appendices. Appendix A provides some properties of the GLK distributions. Appendix B provides proof of the paper's results. Appendix C contains the simulation results, while Appendix D includes more details on the empirical application.
\setcounter{figure}{0} \setcounter{equation}{0} \setcounter{table}{0}