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 (Liesenfeld et al., 2006;Aknouche et al., 2021), economics (Freeland & McCabe, 2004;Berry & West, 2020), social sciences (Pedeli & Karlis, 2011), sports (Shahtahmassebi & Moyeed, 2016), image processing (Afrifa-Yamoah & Mueller, 2022) and oceanography (Cunha et al., 2018). Among the modelling approaches, integer-valued autoregressive processes (INAR), introduced independently by Al-Osh & Alzaid (1987) and McKenzie (1985), become very popular. The stochastic construction of the INAR relies on the binomial thinning operator and the properties of the model on the discrete selfdecomposability of the stationary distribution of the process (Steutel & van Harn, 1979). See Scotto et al. (2015) for a review. The original INAR model has been studied further in Al-Osh & Alzaid (1987) and extended in different directions. (McKenzie, 1986) introduced an INAR model with negative-Binomial and geometric marginal distributions, Jin-Guan & Yuan (1991) extended the INAR(1) model of Al-Osh & Alzaid (1987) to the higher order INAR(p). Al-Osh & Aly (1992) 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 (e.g., see Bourguignon et al., 2016). INAR models with values in the set of signed integers have been propose firstly by Kim & Park (2008) and generalised by Alzaid & Omair (2014) and Andersson & Karlis (2014). Freeland (2010) 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. Alzaid & Al-Osh (1993) propose integer-valued ARMA process with Generalized Poisson marginals and Kim & Lee (2017) 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 (e.g., see Consul & Famoye, 2006, ch. 12). It has been extended further by Janardan (1998) and Janardan (1999), which introduced the four-parameter generalized Pólya-Eggenberger (GPED) distributions of the first and second kind. Janardan (1998) showed that both families contain the Lagrangian Katz distribution as a special case. 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. 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 Weiß & Kim (2013) for a review. Despite the popularity gained in recent years by Bayesian methods, the applications to count data models are still limited (e.g., see McCabe & Martin, 2005;Neal & Subba Rao, 2007;Drovandi et al., 2016;Shang & Zhang, 2018;Garay et al., 2020b). 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 (Chen & Lee, 2016). The posterior distribution of the parameters quantifies uncertainty in the estimation (Chen & Lee, 2017), which can be included in the prediction. The inference from the Bayesian perspective may result in richer inferences in the case of small samples (Garay et al., 2020a) and extra-sample information and in robust inference in the presence of outliers (Fried et al., 2015). Finally, model selection for both nested and non-nested models can be easily carried out. 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 (Frondel et al., 2017), winegrowers (Battaglini et al., 2009), farmers (Fahad & Wang, 2018), mountain peoples (Ullah et al., 2018). Most of these studies rely on surveys conducted in a specific geographical area and sector of society, with a few exceptions. For example, Ziegler (2017) proposed a cross-country analysis of climate change beliefs and attitudes. Lineman et al. (2015) 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, . . ., where (x) k↑ = x(x + 1) . . .
(x + k -1) is the rising factorial with the convention that (x) 0 = 1, and a > 0, c > 0, b ≥ -c and 0 < β < 1 are the parameters (Consul & Famoye, 2006). We denote the distribution with GLK(a, b, c, β). We notice that for -c < b < 0 some additional constraints on the parameters are needed to have all the p x ≥ 0. See the discussion at the beginning of Subsection 3.1 and Appendix Appendix B in the Supplementary. GLK distributions have probability generating function (pgf)

which satisfies:

or alternatively

see Janardan (1998). Remark 1. Building on the Lagrangian expansion, Janardan (1998) introduced the Generalized Polya Eggenberger distribution. (Consul & Famoye, 2006) argued that since the distribution is unrelated to the Polya, it should be named Generalized Lagrangian Katz distribution. As shown in Appendix Appendix B in the Supplementary Material, it is possible to derive the Generalized Polya Eggenberger / Generalized Lagrangian Katz as a particular "generalized Lagrangian distribution". 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 LK(a, b, β) is obtained by replacing c with β (which is called Generalized Katz in (Consul & Famoye, 2006)). • The Katz distribution K(a, β) is obtained for b = 0 and by replacing c with β, (Katz, 1965). • The Polya-Eggenberger distribution PE(a, c, β) 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 1 -I) of Janardan (1998) and can be obtain as the limit of the zero-truncated GLK for a → -c.

• • The Generalized Poisson (GP) distribution GP(θ, λ) for c → 0 s.t. b/a = λ and aβ/c = θ > 0 with 0 < λ < θ -1 . The GP limit of the GLK distribution is stated in (Consul & Famoye, 2006) without proof. In Appendix Appendix B.2 in the Supplementary Material, we provide a proof. The probability mass function of the GLK for different parameter settings is given in We provide in Appendix 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 Appendix A in the Supplementary Material it follows that the coefficient of variation is CV = ((1 -β)/(aθκ)) 1/2 where κ = 1 -β -bβ/c and θ = β/c, assuming κ > 0. The Fisher index is given by the variance-to-mean ratio

) which does not depend on the parameter a. For a given β, following the values of κ (b and c), the distribution allows for various degrees of dispersion: not

We conclude this section with another important property. Proof. From the pgf of a GLK given in Eq. 3

which is the pgf of the sum of n independent GLKs with distribution GLK(a/n, b, c, β)

where a/n > 0 according to the definition of GLK. The Generalized Katz INAR(1) process (GLK-INAR(1)) is defined using the binomial thinning operator, •. The binomial thinning for a non-negative discrete random variable

where B i (α) are iid Bernoulli r.v.s with success probability P (B i (α) = 1) = α. Definition 1 (GLK-INAR process). For α ∈ (0, 1), the GLK-INAR(1) process is defined by

where ε t are iid random variables with Generalized Lagrangian Katz distribution GLK(a, b, c, β),  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. 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):

where κ = 1 -β -bβ/c and θ = β/c. Remark 2. Setting b = 0, c = β = θ 1 and a = θ 2 the results in Kim & Lee (2017) Th. 2.2 are obtained. The process {X t } t∈Z is a Markov Chain on N and the transition probability P i,j = P(X t = j|X t-1 = i) satisfies

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). Proposition 2. Assume that {X t } t∈Z is a GLK-INAR(1). (i) The process {X t } t∈Z is an irreducible, aperiodic and positive recurrent Markov chain. Hence there is a unique stationary distribution for the process {X t } t∈Z . (ii) The marginal distribution of the stationary process {X t } t∈Z is infinitely divisible. Proof. Alzaid (1987). An alternative derivation of (ii) is as follows. Since at stationarity the process satisfies

where ε ∼ GLK(a, b, c, β), and the innovation terms are infinite divisible by Proposition 1, the stationary distribution satisfies the definition of discrete semi-self-decomposability given in Bouzar (2008). Theorem 2 in Bouzar (2008) yields that it is also infinitely divisible. 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 Proposition 3. Let {X jt } t∈Z with j = 1, 2, . . . , J be a sequence of independent GLK-INAR(1) which satisfy:

The process Y t = X 1t + . . . + X Jt is GLK-INAR(1) which satisfies:

Below, we state explicit closed-form expressions for unconditional moments of the process. ε and σ 2 ε the mean, second order non-central moment and variance given in Prop. 5 for a GLK(a, b, c, β).
For a GLK-INAR(1) process, the following unconditional moments can be derived:

(iv) Higher-order non-central moments can be derived using the formula:

ε where s(m, k) and S(m, k) denote the Stirling's numbers of the I and II kind, respectively. Proof. First-and second-order moments are known from Al-Osh & Alzaid (1987) for general INAR. Specifying the GLK innovations gives (i)-(iii). High-order moments can be computed similarly. See e.g. Weiß (2013). For the sake of completeness, details are given in Appendix Appendix B.3 in the Supplementary Material. From the previous proposition, under the assumption κ = 1 -β -bβ/c > 0 one obtains the unconditional variance of the process σ 2

the dispersion index of the process

ε /µ ε is the innovation index of dispersion. It follows that there is under-or over-dispersion in the marginal distribution, V M R X < 1 and V M R X > 1, if and only if there is under-or over-dispersion in the innovation,

The autocorrelation function is

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, . . . , K corresponding to different parameter configurations, i.e. with ε t |S t ∼ GLK(a(S t ), b(S t ), c(S t ), β(S t )), where the thinning coefficient and the GLK parameters of the error term ψ(S t ) = (α(S t ), a(S t ), b(S t ), c(S t ), β(S t )) ′ are time-varying

The S t ∈ {1, . . . , K} for t ∈ Z denotes a hidden Markov-chain process with transition probabilities P(S t = j|S t-1 = i) = π ij for i, j ∈ {1, . . . , K}. From now on, this extension is denoted with

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 α 1 → 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β 1 → 0. The transition probabilities provide information on the persistence of these events. If the probability is independent on past information, i.e. P(S t = j|S t-1 = i) = π j for all i, j ∈ {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 π 11 . Other regimes (S t ̸ = 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 ̸ = 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. Definition 2 (GLK-INARMA(p, q) process). Let α ℓ ∈ (0, 1) for ℓ = 1, . . . , p and ζ r ∈ (0, 1) for r = 1, . . . , q, the GLK-INARMA(p,q) process is defined as

where ε t are iid random variables with Generalized Lagrangian Katz distribution GLK(a, b, c, β)

and p ℓ=1 α ℓ < 1 and q r=1 ζ r < 1. 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. Definition 3 (GLK-MINAR(1) process). Let X t = (X 1t , . . . , X Jt ) ′ be a random integer

, where α ij ∈ (0, 1), the GLK-MINAR(1) process can be defined by

is the thinning matrix operator, such that for each i = 1, . . . , J,

and α ij • X jt-1 refers to the binomial thinning operator. 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 x≥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

Appendix B in the Supplementary Material.). • For parameter values a > 0, b ≥ 0, c > 0 the pmf are positive. Moreover, for Consul & Famoye, 2006, p. 8). • 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↑ 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 ≤ m ≤ q and

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. Remark 4. It should be noted that alternative definitions for -c < b < 0 can be considered. For example one can set to 0 the p x < 0, i.e. when x > max{(a + 1)/|b|. In this case re-scaling the p x is necessary to get x * x=0 p x = 1. The resulting pmf is not a generalized Lagrangian distribution (due to the truncation and rescaling), and the normalizing constant is not in closed form. See, for example, McCabe & Skeels (2020) for a discussion on the parameter values for the Katz distributions. 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 Be(κ, τ ) is the beta distribution with shape parameters κ and τ and Ga(κ, τ ) the gamma distribution with shape and scale parameters κ and τ , respectively. In the empirical applications we assume a non-informative hyper-parameter setting for α and β, that is κ α = τ α = κ β = τ β = 1 and an informative prior for a, b and c with κ a = τ a = 1,

In the case of Markov-switching specification of the GLK-INAR(1), the same prior is assumed for the regime-specific parameters h

For the transition probabilities of the allocation variable, we assume a symmetric Dirichlet prior for each row of the transition matrix, i.e. Let x 1 , . . . , x T be a sequence of observations for the GLK-INAR(1) process, then the joint posterior distribution is given by

where ψ = (α, a, b, c, β) is the parameter vector f (ψ) the joint prior and

where d ijk = 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:

. The probability p i becomes null for i > j if U (ψ) + V (ψ)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 η = (η 1 , . . . , η 5 ) be the 5-dimensional parameter vector obtained by the transformation

), and

) the Jacobian of the transformation φ given above. Given the adaptation parameters µ j and Σ (j) , at the j-th iterations, the AMCMC consists of the following three steps. First, a candidate η * is generated from the random walk proposal: η * = η (j-1) + λ (j) w (j) , w (j) ∼ N q (0, Σ (j) ). Second, the candidate is accepted with probability ρ (j) = ρ(η (j-1) , η * ), where

and third, the adaptive parameters are updated as follows:

where ρ * is the target acceptance probability and γ (j) = j -a , a > 0 is the adaptive scale (Andrieu & Thoms, 2008, , Algorithm 4). Following the suggestions in Roberts et al. (1997) we set ρ * = 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

Notice that the conditioning on the parameters ψ 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 P (S t |S t+1 , X t ) ∝ π StS t+1 P (S t |X t ) and S 1:T = (S 1 , . . . , 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 π

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: α = 0.3, a = 5.3239, b = 0.0592, c = 0.6, β = 0.5917 in the low persistence setting, and α = 0.7, a = 5.3239,b = 0.0592, c = 0.6, β = 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. 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 Supplemen-tary 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 fore-

The posterior distribution of the autoregressive coefficient is given in

The results indicate a deviation from the Negative Binomial model. Thus we apply the DIC criterion DIC = -4E(log f (X|ψ)|y) + 2 log f (X| ψ) to compare GLK-INAR(1)

and NB-INAR(1). The DIC is computed following (Spiegelhalter et al., 2002):

where f (X|ψ) is the likelihood of the model, ψ Given the high kurtosis levels and the multi-modality in the empirical distribution of both series (see Figure 4), 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 D.3 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 6. 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 7 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 7 right plots). We run the analysis at a disaggregate level. The results are given in Figures 89and Tables D.1-D.2 in the Supplementary Material. Figure 8 provides evidence of an inverse relationship between estimated persistence α and dispersion V M R 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 µ ε /(1-α) and the dispersion index V M R 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 Lineman et al. (2015). We investigate the relationships in the search volumes through the lens of our GLK-INAR(1) model. The left plot in Fig. 9 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 D.1-D.2 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 D.1-D.2 in the Supplementary Material). Tables D.1-D.2 in the Supplementary Material report the marginal likelihood of the

(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 (Eckstein et al., 2021, see). The left plot in Fig. 10 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. 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  The skewness and the kurtosis of the distribution are

respectively. For a given value of θ, there is negative skewness if β < (1 + ξ)/(2κ + ξ)

with ξ = 3θ(b + c)a -1/2 and positive otherwise. where we set ℓ = j to get the result and the following equivalent Lagrangian expansion used in Janardan ( 1998)

In particular, if f (1) = g(1) = 1, the function u → f (z(u)) defines the pgf of the "generalized Lagrangian distribution"

, the expressions in ( 2) and ( 1) follows after some algebra as detailed in the following. The expansion coefficients become

Hence

while, for k ≥ 2, the k-th coefficient of the Lagrangian expansion in Eq. B.2 is

where ξ k = a c + k b c + 1 and (x) k↑ = x(x + 1) . . . (x + k -1) is the rising factorial. We now discuss conditions for which p k ≥ 0 for all k ≥ 0. (iv) Let us denote with (x) m = x(x -1) . . . (x -m +1) the falling factorial and with µ (k) = E((X) k ) the m-order falling factorial moment of a random variable X. The following two results will be used. The relationships between non-central moments and falling factorial moments are

where s(m, k) and S(m, k) are the Stirling numbers of the I and II kind, respectively (e.g., see Consul & Famoye, 2006, p. 18). Let X and Y be two random variables then

which can be proved by induction. Let α • X a binomial thinning with X a discrete random variable, then

where |κ| = κ 1 + . . . + κ X . Using the results given above and stationarity (i.e. E(

which implies the m-order falling factorial moment of a INAR(1) is

and the m-order moment is and Geweke's convergence diagnostic (CD) of the posterior MCMC samples for the two settings: low persistence and high persistence. We ran the proposed MCMC algorithm for 50,000 iterations and evaluate the statistics before (subscript BT) and after (subscript AT) removing the first 10,000 burn-in samples, and applying a thinning procedure with a factor of 10. In parenthesis the p-values of the Geweke's convergence diagnostic. High persistence

Figure 6: Posterior approximation of the persistence parameter for a three-state Markov Switching model: high (right plot), medium (middle) and low (left) persistence, for the Climate Change (top) and Global Warming (bottom) datasets. .1: Autocorrelation function (ACF), effective sample size (ESS), inefficiency factor (INEFF), numerical standard errors (NSE)

Preprint submitted to ElsevierDecember 18, 2024