# Tail Index Estimation for Discrete Heavy-Tailed Distributions with Application to Statistical Inference for Regular Markov Chains

Patrice Bertail, Stephan Clémençon, Carlos Fernández

► To cite this version:

Patrice Bertail, Stephan Clémençon, Carlos Fernández. Tail Index Estimation for Discrete Heavy-Tailed Distributions with Application to Statistical Inference for Regular Markov Chains. Test, 2025, 34, pp.691-713. 10.1007/s11749-025-00975-9. hal-04299756v2

HAL Id: hal-04299756

https://hal.science/hal-04299756v2

Submitted on 22 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Tail index estimation for discrete heavy-tailed distributions with application to statistical inference for regular markov chains

Patrice Bertail<fnref n="1" /> · Stephan Clémençon<fnref n="2" /> · Carlos Fernández<fnref n="1" />,<fnref n="2" />

Received: 18 September 2024 / Accepted: 6 May 2025 / Published online: 2 June 2025
© The Author(s) 2025

## Abstract

It is the purpose of this paper to investigate the issue of estimating the regularity index $\beta > 0$ of a discrete heavy-tailed r.v. $S$, i.e. a r.v. $S$ valued in $\mathbb{N}^*$ such that $\mathbb{P}(S > n) = L(n) \cdot n^{-\beta}$ for all $n \ge 1$, where $L: \mathbb{R}_+^* \to \mathbb{R}_+$ is a slowly varying function. Such discrete probability laws, referred to as generalized Zipf’s laws sometimes, are commonly used to model rank-size distributions after a preliminary range segmentation in a wide variety of areas such as e.g. quantitative linguistics, social sciences or information theory. As a first go, we consider the situation where inference is based on independent copies $S_1, \dots, S_n$ of the generic variable $S$. The estimator $\hat{\beta}$ we propose can be derived by means of a suitable reformulation of the regularly varying condition, replacing $S$’s survivor function by its empirical counterpart. Under mild assumptions, a non-asymptotic bound for the deviation between $\hat{\beta}$ and $\beta$ is established, as well as limit results (consistency and asymptotic normality). Beyond the i.i.d. case, the inference method proposed is extended to the estimation of the regularity index of a regenerative $\beta$-null-recurrent Markov chain. Since the parameter $\beta$ can be then viewed as the tail index of the (regularly varying) distribution of the return time of the chain $X$ to any (pseudo-) regenerative set, in this case, the estimator is constructed from the successive regeneration times. Because the durations between consecutive regeneration times are asymptotically independent, we can prove that the consistency of the estimator promoted is preserved. In addition to the theoretical analysis carried

P. Bertail, S. Clémençon and C. Fernández: These authors contributed equally to this work

<figure><img src="image_2.png" /></figure>

Carlos Fernández
fernandez@telecom-paris.fr

Patrice Bertail
patrice.bertail@parisnanterre.fr

Stephan Clémençon
stephan.clemençon@telecom-paris.fr

<footnote n="1">¹ MODAL'X, UMR CNRS 9023, Université Paris Nanterre, 200 Avenue de la République, 92000 Nanterre, Ile-de-France, France</footnote>

<footnote n="2">² LTCI, Telecom Paris, Institut Polytechnique de Paris, 19 place Marguerite Perey, 91123 Palaiseau, Ile-de-France, France</footnote>
out, simulation results provide empirical evidence of the relevance of the inference technique proposed.

**Keywords** Generalized discrete Pareto distribution · Nonparametric estimation · Null-recurrent Markov chain · Regularity index · Zipf's law

**Mathematics Subject Classification** 60K35

# 1 Introduction

This article is devoted to the study of the problem of estimating the regularity index $\beta > 0$ of a generalized discrete Pareto distribution, namely the probability distribution of a random variable $S$ defined on a probability space $(\Omega, \mathcal{F}, \mathbb{P})$, taking its values in $\mathbb{N}^*$ and such that:

$$ \mathbb{P}(S > n) = n^{-\beta} L(n) \quad \text{for all } n = 1, 2, \dots, \qquad (1) $$

where $L: \mathbb{R}_+ \to \mathbb{R}$ is a slowly varying function, i.e. such that $L(\lambda z)/L(z) \to +1$ as $z \to +\infty$ for any $\lambda > 0$, see Bingham et al. (1987). Such discrete power law probability distributions, also referred to as generalized Zipf's laws sometimes, are often used to model the distribution of discrete data exhibiting a specific rank-frequency relationship, namely when the logarithm of the frequency and that of the rank order are nearly proportional. Such a phenomenon has been empirically observed in many ranking systems: in quantitative linguistics (i.e. when analysing word frequency law in natural language, see e.g. Manning and Schütze 1999) in the first place, as well as in a very wide variety of situations, too numerous to be exhaustively listed here. One may refer to Sidra et al. (2018), Lazzardi et al. (2021) or Zanette (2006) among many others. In this paper, we first consider the issue of estimating the parameter $\beta$ involved in (1) (supposedly unknown, like the function $L$) in the classic (asymptotic) i.i.d. statistical setting, i.e. based on an increasing number $n \ge 1$ of independent copies $S_1, \dots, S_n$ of the generic r.v. $S$. Statistical inference for discrete heavy-tailed distributions has not received much attention in the literature and still poses methodological problems that this article seeks to resolve. Most of the very few dedicated methods documented either deal with very specific cases as in e.g. Goldstein et al. (2004), Matsui et al. (2013) and Clauset et al. (2009) or else consist in empirically applying techniques originally designed for continuous heavy-tailed distributions to discrete data after first adding independent uniform noise (Voitalov et al. 2019; Kim and Kokoszka 2020, 2023).

The vast majority of the regular variation index estimators proposed in the literature, Hill or Pickand estimators in particular (Hill 1975; Pickands 1975), are based on order statistics, which causes obvious difficulties in the discrete case because of the possible occurrence of many ties. Indeed, in that case, many spacings are equal to 0 which causes these estimators to behave erroneously. This is stressed in Matsui et al. (2013) (see their figures 1 and 2). More importantly, Matsui et al. (2013) constructs an example of a family of discrete value heavy-tailed distributions that do not satisfy any second-order condition so that the asymptotic normality of the Hill estimator

<figure><img src="image_3.png" /></figure>
cannot be proved directly. In contrast, the estimator under study here is based on the analysis of the probability of exponentially separated tail events. It simply rests on the fact that, as can be immediately deduced from (1), we have $\ln(p_k) - \ln(p_{k+1}) = \beta + \ln(L(e^k)/L(e^{k+1}))$, where $\ln(x)$ denotes the natural logarithm of any real number $x > 0$ and $p_k = \mathbb{P}(S > e^k)$ for all $l \in \mathbb{N}$, and that $L(e^{k+1})/L(e^k)$ is expected to be very close to 1 for $k \in \mathbb{N}$ chosen sufficiently large. A natural (plug-in) inference technique can be then devised by replacing the tail probabilities $p_k$ with their empirical versions $\hat{p}_l^{(n)} = (1/n) \sum_{i=1}^n \mathbb{I}\{S_i > e^k\}$ for $k \in \mathbb{N}$, where $\mathbb{I}\{A\}$ means the indicator function of any event $A$. This yields the estimator

$$ \hat{\beta}_n(k) = \ln(\hat{p}_k^{(n)}) - \ln(\hat{p}_{k+1}^{(n)}), \quad (2) $$

provided that $\hat{p}_{k+1}^{(n)} > 0$ (as shall be seen, this occurs with large probability if $n$ is sufficiently large). By convention, we set $\hat{\beta}_n(k) = 0$ when $\hat{p}_{k+1}^{(n)} = 0$. We point out that it has exactly the same form as that proposed and analysed in Carpentier and Kim (2015) in a different context, that of approximately Pareto distributions<fnref n="1" /> namely. In the discrete generalized Pareto framework, we prove that for an appropriate choice of the hyperparameter $k = k_n$ (typically chosen of order $\ln(n)$), the estimator (2) is strongly consistent and asymptotically normal as $n \to +\infty$. We also show that non-asymptotic upper confidence bounds for the absolute deviations between $\hat{\beta}_n(k)$ and $\beta$ are also established here.

Although estimation of the parameter $\beta$ in (1) in the discrete i.i.d. setting is an important issue in itself, the present paper also finds its motivation in the problem of recovering statistically the regularity index of a regenerative regular $\beta$-null-recurrent Markov chain $X = (X_n)_{n \in \mathbb{N}}$, based on the observation of a finite sample path $X_1, \dots, X_n$ with $n \ge 1$. As explained in Chen (1999, 2000), for regular Markov chains, the regularity index $\beta$ controls the (sublinear) rate at which the number of visits to any given Harris set increases with observation time $n$, no matter the initial distribution. These Markov chains are strongly non-stationary and are not mixing in any sense. Typical examples of these chains are random walks and their variations, for instance Bessel random walks (see Alexander 2011). Such processes appear also naturally in many applications: for population dynamics or stochastic growth models in demography, see for instance Adam (2016) and the references therein; for branching processes, Galton-Watson processes with immigration (or population dependent) also exhibit a fat tail behaviour of the time returns of the chain, see (Pakes 1971; Zubkov 1972; Klebaner 1993).

In the regenerative case (i.e. when the chain $X$ possesses an accessible atom, a Harris set on which the transition probability is constant), the distribution of the regenerative time, the return time to the atom, is a discrete generalized Pareto (1) and the parameter $\beta$ is its tail index. Due to the non-standard behaviour of traditional estimators in this context, statistical inference for null-recurrent Markov chains is very poorly documented in the literature (Karlsen and Tjostheim 2001; Karlsen et al. 2010; Myklebust

<footnote n="1">¹ The distribution of a real-valued r.v. $X$ is said to be approximately Pareto with tail index $\beta > 0$ if its survivor function is of the form: $\forall x > 0, \mathbb{P}(X > x) = L(x)x^{-\beta}$, where $L$ is asymptotically constant at infinity, i.e. there exists $C \in (0, \infty)$ s.t. $L(x) \to C$ as $x \to +\infty$.</footnote>
et al. 2012; Gao et al. 2013) and, to the best of our knowledge, estimation of the key quantity $\beta$ has not received much attention. It is also the goal of this article to extend the use of the estimator (2) to the case where the $S_i$'s are the successive durations between the consecutive regeneration times up to time $n$. The main difficulty naturally arises from the fact that the number $1 + N_n \ge 0$ of regeneration times (and thus the number of durations) is now random, and it is not independent of the variables $S_1, \dots, S_{N_n}$ (in particular, $S_1 + \dots + S_{N_n} \le n$ by construction). We show that the estimator remains strongly consistent in this scenario and we provide non-asymptotic upper confidence bounds for the absolute deviations between $\hat{\beta}_{N_n}(k)$ and $\beta$. For illustration purposes, numerical experiments have been carried out, providing empirical evidence of the relevance of the estimation method promoted. Extension to the general case of (pseudo-regenerative) null-recurrent chains is also discussed; the difficulties inherent in applying the methodology originally proposed in Bertail and Clémençon (2006b) in the positive recurrent case to mimic regenerative Nummelin extensions (Nummelin 1984) being explained at length.

The paper is organized as follows. A thorough analysis of the behaviour of the estimator (2) in the i.i.d. case, illustrated by numerical experiments, is first carried out in section 2. Some of the results thus established (consistency and non-asymptotic bounds) are next extended in section 3 to the regenerative regular Markovian setup, when the estimator is computed based on a single finite-length trajectory of the atomic chain. Experimental results are also displayed, and the main barrier to the extension of the methodology promoted to general (i.e. pseudo-regenerative) regular null-recurrent chains is also discussed therein. Technical proofs are deferred to the Supplementary material.

## 2 Tail Index Estimation—The Discrete Heavy-Tailed i.i.d. Case

Throughout this section, $S_1, \dots, S_n$ are independent copies of a generic discrete generalized Pareto r.v. $S$, i.e. a r.v. $S$ with survivor function of type (1), where the parameter $\beta > 0$ and the slowly varying function $L$ are supposedly unknown. As a first go, we start to investigate the (asymptotic) behaviour of the estimator (2) in this basic general framework and next develop the analysis in particular situations, i.e. when the function $L$ has a specific form.

### 2.1 Main Results—Confidence Bounds and Limit Theorems

As explained in the Introduction section, the estimator (2) can be viewed as an empirical counterpart of the quantity

$$ \beta(k) = \ln(p_k) - \ln(p_{k+1}) = \beta + \ln \left( \frac{L(e^k)}{L(e^{k+1})} \right), \quad (3) $$

see (1), which tends to $\beta$ as $k \to \infty$ by virtue of the slow variation property of $L$. It was initially proposed in Carpentier and Kim (2015), where only the case when $L$ is

<figure><img src="image_4.png" /></figure>
asymptotically constant is analysed. Finite sample inequalities and strong consistency are established therein, as well as uniform rates of convergence under the stronger assumption that $F$ is second-order Pareto (i.e. $L(x) = C_0+O(x^{-\alpha\beta})$ for some positive constants $\alpha$ and $C_0$). In what follows, we extend their finite sample inequalities and strong consistency results to the entire class of distributions with regularly varying survival functions and prove the asymptotic normality of the estimator.

The proposition below provides an upper confidence bound for the absolute deviations between (2) and $\beta$.

**Proposition 2.1** Let $\delta \in (0, 1/2)$ and set $u_n(\delta) = \ln(2/\delta)/n$ for all $n \ge 1$. If $k \ge 1$ is such that $p_{k+1} \ge 16u_n(\delta)$, then, with probability at least $1-2\delta$, we have:

$$|\hat{\beta}_n(k) - \beta| \le 6 \sqrt{\frac{u_n(\delta)}{p_{k+1}}} + \left|\ln\left(\frac{L(e^k)}{L(e^{k+1})}\right)\right|. \quad (4)$$

The bound (4) reveals some sort of “bias-variance” trade-off controlled by the hyperparameter $k > 0$. The second term on the right hand can be viewed as the bias of the inference method, as the estimator (2) serves as an empirical counterpart to the approximation (3). This bias term diminishes as $k$ increases towards infinity. Conversely, the first term, whose presence is due to the random nature of the estimator, gets bigger as $k$ increases (although it cannot get too big due to the condition $p_{k+1} \ge 16u_n(\delta)$). It is important to note that second-order slow variation conditions (Goldie and Smith 1987) are required to bound the (vanishing) bias term in (4), as shall be explained in subsection 2.2. The following result reveals that for an appropriate choice of $k = k_n$, the estimator (2) is strongly consistent.

**Theorem 2.1** (Strong consistency). Suppose that, as $n \to +\infty$, we have $k_n \to +\infty$ and $(\ln n) \exp(k_n\beta)/n = o(L(\exp(k_n)))$. Then, we have:

$$\hat{\beta}_n(k_n) \to \beta \text{ almost-surely, as } n \to +\infty.$$

In particular, as stated below, strong consistency is guaranteed when $k_n$ is of logarithmic order.

**Corollary 2.1** Let $0 < A < 1/\beta$. Then, we have:

$$\hat{\beta}_n(A \ln n) \to \beta \text{ almost-surely, as } n \to +\infty.$$

Now, the following results establish the asymptotic normality of the deviation between (2) and $\beta(k_n)$, when appropriately normalized.

**Theorem 2.2** (Asymptotic normality). Suppose that $k_n$ satisfies the conditions of Theorem 2.1, then

(i) Then, as $n \to +\infty$, we have the convergence in distribution:

$$\sqrt{np_{k_n}} (\hat{\beta}_n(k_n) - \beta(k_n)) \Rightarrow \mathcal{N}(0, e^\beta - 1).$$
(ii) In addition, asymptotic normality holds true for the “standardized” deviation:

$$ \frac{\sqrt{n \hat{p}_{k_n}^{(n)} (\hat{\beta}_n(k_n) - \beta(k_n))}}{\sqrt{e^{\hat{\beta}_n(k_n)} - 1}} \Rightarrow \mathcal{N}(0, 1), \text{ as } n \to +\infty. $$

The asymptotic normality results above can be extended to the deviation between (2) and $\beta$, provided that the bias term $\beta(k_n) - \beta$ vanishes at an appropriate rate, as stated below.

**Corollary 2.2** Suppose that the conditions of Theorem 2.2 are fulfilled. In addition, assume that $k_n$ is such that

$$ \sqrt{np_{k_n}} \left( 1 - \frac{L(e^{k_n})}{L(e^{k_n+1})} \right) \to 0, \text{ as } n \to +\infty. \quad (5) $$

(i) Then, we have the convergence in distribution

$$ \sqrt{np_{k_n}} (\hat{\beta}_n(k_n) - \beta) \Rightarrow \mathcal{N}(0, e^{\beta} - 1) \text{ as } n \to +\infty. $$

(ii) In addition, the “studentized” version is asymptotically normal:

$$ \frac{\sqrt{n \hat{p}_{k_n}^{(n)} (\hat{\beta}_n(k_n) - \beta)}}{\sqrt{e^{\hat{\beta}_n(k_n)} - 1}} \Rightarrow \mathcal{N}(0, 1) \text{ as } n \to +\infty. $$

Of course, the condition (5) on $k_n$ can be difficult to check in practice. This is a common issue in tail estimation and in the statistical analysis of extreme values more generally. The choice of the hyperparameter $k$ somehow governs the (asymptotic) bias-variance trade-off: The estimator (2) is expected to have a large variance when $k$ is large and to have a large bias if $k$ is too small. As depicted in Fig. 1a, to choose $k$, one may use the same approach as that originally proposed for the Hill estimator (see e.g. Resnick 2007), which consists of plotting the values of (2) for a range of values of $k$ and choosing $k$ in a region where the estimator shows some stability. Section 3.2 of Carpentier and Kim (2015) suggests an adaptive algorithm for selecting $k$ for $\hat{\beta}_n$, assuming the stricter condition that $F$ belongs to the class of second-order Pareto distributions.

We point out that, contrary to (variants of) the Hill inference method, the integer hyperparameter $k$ does not refer here to the number of order statistics involved in the estimator. Instead, it determines a threshold that defines the tail probability we want to estimate. This implies that the range of $k$ is the set of integers between 1 and $\max_{1 \le i \le n} \ln S_i - 1$, while, for classical tail index estimators, the hyperparameter varies between 1 and $n$.

**Remark 1** (QUANTILES VS TAIL PROBABILITIES) In extreme value theory, it is customary to estimate the tail index by looking at the largest values, i.e. a fraction of the order statistics. The general idea consists of fixing the tail probabilities and then using
<figure><img src="image_5.png" /><figcaption>Fig. 1 (a) Behaviour of $\hat{\beta}_n(k)$ for different values of $k$, to estimate the parameter $\beta = 0.8$ based on a dataset of $10^6$ independent realizations of a Zeta distribution (A discrete r.v. $W$ follows a Zeta distribution with parameter $\beta$ if $P(W = k) = 1/(k^{\beta-1}\zeta(\beta-1))$ where $\zeta$ is the Riemann zeta function. Its cdf satisfies $P(W \ge k) \sim k^{\beta}/(\beta\zeta(\beta-1))$. This distribution is also known as Zipf's distribution due to its relationship with Zipf's law.). (b) Behaviour of the Hill estimator for the same dataset. Notice that the range of the x-axis in both graphs is different; in the case of $\hat{\beta}_n(k)$, the range of $k$ is the set of integers between 1 and $\max_{1 \le i \le n} \ln S_i - 1$ while in the case of the Hill method, the x-axis is the number of order statistics involved in the estimator; hence, it goes from 1 to $n$</figcaption></figure>

the order statistics (i.e. empirical quantiles) to estimate the quantiles (De Haan and Ferreira 2007, Chapter 3). The estimator proposed here follows the converse approach: we fix the thresholds/quantiles, next use the data to estimate the tail probabilities and finally use these estimations in order to mimic the limit behaviour of (3). In general, estimators based on order statistics are difficult to analyse in a non-asymptotic way, whereas our estimate is based on probabilities of well-chosen tail events.

**Remark 2 (AVERAGED VERSIONS)** From a practical point of view, rather than picking a single value for $k$, another natural approach would consist in averaging the estimators (2) over a range of values for the hyperparameter. Let $k$ and $m$ be such that $k > m$ and define

$$ \beta(k, m) = \frac{1}{2m+1} \sum_{j=-m}^{m} \beta(k+j), \quad \hat{\beta}_n(k, m) = \frac{1}{2m+1} \sum_{j=-m}^{m} \hat{\beta}_n(k+j). $$

One may easily check that

$$ \beta(k, m) = \beta + \frac{1}{2m+1} \left| \ln \left( \frac{L(e^{k-m})}{L(e^{k+m+1})} \right) \right|. \quad (6) $$

The non-asymptotic bound in Proposition 2.1 can be extended to the averaged version, as revealed by the analysis carried out in Section B of the Supplementary Material, as well as the strong consistency and asymptotic normality results. However, the asymptotic variance of the averaged version is shown to increase with $m$.
**Remark 3** (SECOND- ORDER CONDITIONS) According to (De Haan and Ferreira 2007, Remark 3.2.6), under the following second-order condition:

$$ \forall t > 0, \alpha > 0, \quad \lim_{x \to +\infty} \left( \frac{\bar{F}(tx)}{\bar{F}(x)} - t^{-\beta} \right) x^{\alpha} = 0 \qquad (7) $$

the Hill estimator, properly standardized, converges in distribution to a centred normal random variable. We point out that if the distribution of S satisfies (7), the sequence $k_n$ satisfies the hypothesis of Theorem 2.2, and there is a positive constant C such that $n \le e^{Ck_n}$ for n large enough, then the conditions of Corollary 2.2 are satisfied.

As explained in Example 1, it is possible to build a discrete heavy tail distribution that do not satisfy any second order condition (as highlighted in Matsui et al. (2013)) but for which asymptotic normality of (2) remarkably holds true.

**Example 1** Let $l \in \mathbb{N}$, $\beta > 0$ and U be a uniform distribution in $[0, 1]$. Consider the random variable $X = 10^{-l} \lfloor 10^l U^{-1/\beta} \rfloor$. This transformation of $U^{-1/\beta}$ turns all but the first $l$ digits behind the comma into zeros. Equation (1.6) to Matsui et al. (2013) shows that the survival function of this random variable is given by

$$ \bar{F}(x) = \frac{(10^l)^\beta}{(\lfloor 10^l x \rfloor + 1)^\beta}, \qquad (8) $$

and that $\bar{F}(x) \sim x^{-\beta}$ as $x \to \infty$. Moreover, the authors showed that this distribution does not satisfy the second order condition (7). Define $L(x) = (10^l x)^\beta (\lfloor 10^l x \rfloor + 1)^{-\beta}$, then $L(x) \to 1$ as $x \to \infty$, hence it is slowly varying, and equation (8) can be written as $\bar{F}(x) = x^{-\beta} L(x)$. In section A.5 of the Supplementary Material, we show that there exists a positive constant K such that, for x big enough

$$ \left| \frac{L(x)}{L(ex)} - 1 \right| \le K x^{-1}, \qquad (9) $$

therefore, if $k_n \to +\infty$, for n large enough, we have

$$ \sqrt{np_{k_n}} \left| 1 - \frac{L(e^{k_n})}{L(e^{k_n+1})} \right| \le K \sqrt{\frac{nL(e^{k_n})}{e^{k_n}\beta}} (e^{-k_n}) \le 2K \sqrt{n e^{-k_n}(\beta+2)}. \qquad (10) $$

This shows that if $k_n$ satisfies the conditions of Theorem 2.1 and $ne^{-k_n(\beta+2)}$ goes to 0, then the conditions of Corollary 2.2 are satisfied and $\sqrt{np_{k_n}}(\hat{\beta}_n(k_n) - \beta)$ converges in distribution to a centred normal random variable. Moreover, if we take $k_n = A \ln n$ (as in Corollary 2.1), then the asymptotic normality is guaranteed as long as $1/(\beta+2) < A < 1/\beta$.

In the next subsection, we discuss further how the behaviour of the slowly varying function L impacts the “bias-variance” contributions revealed by the bound (4).

<figure><img src="image_6.png" /></figure>
## 2.2 Refined 'Bias vs Variance' Analysis-Examples

We now consider several specific cases of distributions of type (1) (i.e. several instances of the slowly varying functions $L$) to explicit the asymptotic order of magnitude of the terms $1/\sqrt{np_{k+1}}$ and $|\ln(L(e^k)/L(e^{k+1}))|$ involved in the bound (4), when $k_n$ is picked as in Corollary 2.1: $k_n = A \ln n$ with $0 < A < 1/\beta$.

* **The logarithmic case:** Suppose that $L(n) = C \ln n$, where $C > 0$. In this situation, we have $|\ln(L(e^{k_n})/L(e^{k_n+1}))| \sim 1/(A \ln n)$ as $n \to +\infty$, whereas $1/\sqrt{np_{k+1}} = O(1/\sqrt{n^{1-A\beta}\ln n})$.

* **The inversely logarithmic case:** Consider now the situation where $L(n) = C/\ln n$ with $C > 0$. Then, we still have we have $|\ln(L(e^{k_n})/L(e^{k_n+1}))| \sim 1/(A \ln n)$, while $1/\sqrt{np_{k+1}} = O(\sqrt{(\ln n)/n^{1-A\beta}})$ as $n \to +\infty$.

We point out that, in the two examples above, the conditions of Corollary 2.2 are not met, the bias being too big to get asymptotic normality (centred at $\beta$).

* **The asymptotically constant case:** Suppose that $L(n) = C_0(1 + \epsilon(n))$ where $C_0 > 0$ and $\epsilon(n) \to 0$ as $n \to +\infty$. In this case, $|\ln(L(e^{k_n})/L(e^{k_n+1}))| = O(\epsilon(n^A))$ and $1/\sqrt{np_{k+1}} = O(1/\sqrt{n^{1-A\beta}})$. Hence, if $|\epsilon(n)| = O(n^{-\lambda})$ for some $\lambda > 0$, then the conditions of Corollary 2.2 are satisfied if we take $k_n = A \ln n$ with $\frac{1}{\beta+2\lambda} < A < 1/\beta$.

* **Second-order Pareto distributions:** A distribution function belongs to the class of second-order Pareto distributions, if there exist positive constants $C_0, C_1, \alpha, \beta$ such that $|\bar{F}(x) - C_0x^{-\beta}| \le C_1x^{-\beta(1+\alpha)}$. Clearly, if $F$ is second-order Pareto, then $\bar{F}$ is regularly varying with index $\beta$ and $L(x) = C_0 + O(x^{-\alpha\beta})$. By our previous example, we now have that, if we take $k_n = A \ln n$, then $|\ln(L(e^{k_n})/L(e^{k_n+1}))| = O(n^{-A\alpha\beta})$ and $1/\sqrt{np_{k+1}} = O(1/\sqrt{n^{1-A\beta}})$. Moreover, choosing $A$ such that $\frac{1}{\beta(1+2\alpha)} < A < 1/\beta$ guarantees that the conditions of Corollary 2.2 are satisfied.

Moreover, if we take $k_n = \frac{1}{\beta(2\alpha+1)}\ln(n)$, as in Theorem 2 of Carpentier and Kim (2015), then $\sqrt{np_{k_n}} \sim \sqrt{C}n^{-\alpha/(1+2\alpha)}$ and $\lim_n \sqrt{np_{k_n}} \ln \frac{L(e^{k_n})}{L(e^{k_n+1})} = K$ for some positive constant $K$, which combined with our Theorem 2.2 gives that $\sqrt{np_{k_n}}(\beta_n(k_n) - \beta) \Rightarrow \mathcal{N}(K, e^\beta - 1)$, which gives the same rate of convergence obtained in Carpentier and Kim (2015).

* **Slow variation with a remainder (SR2):** Consider the case where the slowly varying function satisfies the condition SR2 introduced in Bingham et al. (1987): there exist two real-valued functions $k$ and $g$ defined on $\mathbb{R}_+$ such that, for all $\lambda > 0$,

$$ \frac{L(\lambda x)}{L(x)} - 1 \sim \kappa(\lambda)g(x), \text{ as } x \to +\infty, \quad (11) $$

where $\kappa(\lambda) = c \int_1^\lambda \theta^{\rho-1}d\theta$, $c > 0$ and $g$ is regularly varying with index $\rho \le 0$, i.e. $g(x) = x^\rho U(x)$ where $U$ is a slowly varying function. Under the additional assumption that $g$ has positive decrease, Corollary 3.12.3 in Bingham et al. (1987) gives the following representation:
<figure><img src="image_7.png" /><figcaption>Fig. 2 Monte Carlo average and 95% confidence interval for $\hat{\beta}_{10^4}(k)$, as $k$ is varying. The true value of $\beta$ in both cases is 0.5</figcaption></figure>

$$L(x) = C (1 - c|\rho|^{-1} g(x) + o(g(x))), \quad \text{as } x \to +\infty, \quad (12)$$

where $C$ is a finite constant. The result below provides precise control of the bias of the estimation method in this case.

**Lemma 2.1** Suppose that conditions (11) and (12) are fulfilled. Then, as $n \to +\infty$, we have:

$$\ln \left( \frac{L(n^A)}{L(en^A)} \right) = -c|\rho|^{-1} n^{-A|\rho|} \left( U(n^A) - e^{-|\rho|} U(en^A) \right) + o \left( n^{-A|\rho|} U(n^A) \right).$$

In this situation, the bias of the method is thus of order $O(n^{-A|\rho|})$, while $1/\sqrt{np_{k+1}}$ is of order $O(n^{-(1-A\beta)/2})$. Hence, if $1/(\beta + 2|\rho|) \le A < 1/\beta$, the conditions of Corollary 2.2 are satisfied with $k_n = A \ln n$.

To illustrate this trade-off, we present the following Monte Carlo experiment: we generate $10^4$ samples of a heavy-tailed distribution and calculate $\hat{\beta}_{10^4}(k)$ for all admissible values of $k$, we repeat this experiment 100 times, and then we calculate the mean and the 95% confidence interval of $\hat{\beta}_{10^4}(k)$ for each value of $k$. The results of these simulations, for the cases where $L(n)$ is asymptotically constant and $L(n)$ is logarithmic, are presented in Figures 2a and 2b. As expected, the behaviour of the estimator is way better in the former case than in the latter.

### 3 Regular Null-Recurrent Chains—Regularity Index Estimation

We start by setting out the notations used throughout this section, now standard in the Markov chain literature and listing first the properties supposedly satisfied by the class of Markov chains under study. One may refer to Meyn et al. (2009) for an excellent account of the Markov chain theory. The concept of $\beta$-regularity for describing how fast a Harris chain returns to Harris sets is then recalled, together with related asymptotic properties, invoked in the subsequent statistical analysis, for clarity's sake. Then, the main results of this paper, related to the inference of the
parameter $\beta$ and the extended use of the estimator (2) in the (regenerative) Markovian case, are established and discussed. Here, $X = (X_n)_{n \in \mathbb{N}}$ denotes a time-homogeneous Markov chain, with state space $E$, equipped with a countably generated $\sigma$-field $\mathcal{E}$, and transition probability $\Pi(x, dy)$. For any probability distribution $\nu$ on $E$, we denote by $\mathbb{P}_\nu$ the probability distribution on the underlying space such that $X_0 \sim \nu(dx)$ and by $\mathbb{E}_\nu[\cdot]$ the corresponding expectation. For notational convenience, we shall write $\mathbb{P}_x$ and $\mathbb{E}_x[\cdot]$ when $\nu$ is the Dirac mass at $x \in E$. In the following, we also denote by $t \in \mathbb{R} \mapsto \lfloor t \rfloor$ the floor function and by $\Gamma(z) = \int_{t \ge 0} t^{z-1} e^{-t} dt$ the Gamma function.

## 3.1 Background and Preliminaries

Throughout the section, we suppose that the chain $X$ is $\psi$-irreducible, meaning that there exists some $\sigma$-finite measure $\psi$ on $(E, \mathcal{E})$ such that any measurable set $B \subset E$, weighted by $\psi$, can be reached by the chain with positive probability in a finite number of steps, i.e. $\sum_{n \ge 1} \Pi_n(x, B) > 0$, no matter the starting point $x \in E$, denoting by $\Pi_n(x, dy)$ the $n$-th iterate of the transition probability $\Pi(x, dy)$. Recall that an irreducibility measure is said to be maximal if it dominates any other irreducibility measure. We also assume that $X$ is aperiodic (rather than replacing $\Pi$ by an iterate) and Harris recurrent, i.e. that, with probability one, it visits an infinite number of times any measurable subset $B \subset E$, weighted by maximal irreducibility measures, whatever the initial state: $\forall x \in E$, $\mathbb{P}_x(\sum_{n=1}^\infty \mathbb{I}\{X_n \in B\} = \infty) = 1$. When Harris recurrent, a transition kernel $\Pi(x, dy)$ has a nonzero invariant (positive) measure $\mu(dx)$ (i.e. such that $\int_{x \in E} \mu(dx)\Pi(x, dy) = \mu(dy)$), that is unique up to a multiplicative factor (notice incidentally that $\mu(dx)$ is a maximal irreducibility measure). Measurable sets weighted by $\mu$ are said to be Harris. For Harris recurrent chains, recall that the following strong ratio limit theorem holds. We have indeed, as $n \to \infty$,

$$ \frac{\sum_{i=1}^{n} \mathbb{I}\{X_i \in B\}}{\sum_{i=1}^{n} \mathbb{I}\{X_i \in C\}} \to \frac{\mu(B)}{\mu(C)} \quad \mathbb{P}_\nu\text{-almost-surely,} \quad (13) $$

for any initial distribution $\nu$ and any measurable sets $B$ and $C$ s.t. $\mu(C) > 0$. When the measure $\mu(dx)$ is finite, the chain is said to be positive recurrent and, by convention, rather than considering $\mu(dx)/\mu(E)$, by $\mu(dx)$ we mean the stationary probability measure in this case.

**Regular chains.** For a wide class of Harris Markov chains, the regularity index describes how fast the occupation time related to a Harris set $B$ (i.e. the number of visits to $B$)

$$ \Sigma_n(B) = \sum_{i=1}^{n} \mathbb{I}\{X_i \in B\} $$

increases with time $n$. When $X$ is positive recurrent, it follows from the strong law of large numbers that occupation times of Harris sets grow in a linear fashion with the observation time: as $n \to \infty$, $\Sigma_n(B) \sim \mu(B)n \mathbb{P}_\nu$-almost-surely. In the general Harris case, some technical assumptions are required in order to be able to specify
the growing rate. In order to formulate them rigorously, further concepts are required. Recall that a special set (also referred to as a D-set sometimes (Chen 1999)) for the chain X is any Harris set D such that $\mu(D) < \infty$ and $\sup_{x \in E} \mathbb{E}_x[\sum_{i=1}^{\tau_B} \mathbb{I}\{X_i \in D\}] < \infty$, for any Harris set $B \subset E$, denoting by $\tau_B = \inf\{i \ge 1 : X_i \in B\}$ the hitting time to B. We recall that special sets not only exist but there are many of them: actually, any Harris set contains a special set at least, see Proposition 5.13 in Nummelin (1984). For any special set D and initial distribution $\nu$, consider the so-termed truncated Green function:

$$G_{v,D}(t) = \frac{1}{\mu(D)} \sum_{n=1}^{[t]} \nu \Pi_n(D)$$

where $\nu \Pi_n(B) = \int_{x \in E} \nu(dx) \Pi_n(x, B) = \mathbb{P}_\nu(X_n \in B)$ for any $B \in \mathcal{E}$. Harris recurrence entails that $G_{v,D}(t) \to +\infty$ as $t \to \infty$.

In the following, we restrict our attention to a specific class of Harris chains for which the rate at which $G_{v,D}(t)$ grows to infinity as $t \to \infty$ can be characterized. Notice that, in such cases, the rate would be independent from the pair $(\nu, D)$. Indeed, by virtue of Theorem 7.3 in Nummelin (1984), we have $G_{v_1, D_1}(t)/G_{v_2, D_2}(t) \to 1$ as $t$ goes to infinity, for any distributions $v_1$ and $v_2$ and any special sets $D_1$ and $D_2$. One may thus give the following definition, see Chen (1999, 2000).

**Definition 1** ($\beta$-REGULAR MARKOV CHAIN). Let $\beta \in [0, 1]$. A Harris chain X is said to be $\beta$-regular if there exists a special set D and a distribution $\nu$ such that the function $G_{v,D}$ is $\beta$-regularly varying: $\forall t > 0$,

$$\lim_{\lambda \to \infty} \frac{G_{v,D}(\lambda t)}{G_{v,D}(\lambda)} = t^{\beta}. \qquad (14)$$

We point out that property (14) can be rephrased as follows: there exists a slowly varying function $L_{v,D}(t)$ such that $G_{v,D}(t) = L_{v,D}(t)t^{\beta}$. Notice incidentally that “$\beta$-regularity” is called “$\beta$-null recurrence” in Karlsen and Tjostheim (2001) when $\beta < 1$, while $\beta = 1$ corresponds to the positive recurrent case. The parameter $\beta$ thus rules the “frequency” at which a (supposedly regular) Harris chain X recurs (Chen 1999) and it is the purpose of this section to investigate the issue of estimating it with asymptotic guarantees, based on the observation of a single path $X_1, \dots, X_n$ of size $n \to +\infty$. In particular, we shall focus in subsection 3.3 on the case of regenerative chains, for which an extension of the estimator (2) can be used with statistical guarantees.

**Regenerative regular chains.** Recall that a Markov chain is regenerative when it possesses an accessible atom, i.e. a measurable set A such that $\psi(A) > 0$ and $\Pi(x, \cdot) = \Pi(y, \cdot)$ for all $(x, y) \in A^2$. By $\tau_A = \tau_A(1) = \inf\{n \ge 1, X_n \in A\}$ is meant the hitting time to A and we denote by $\tau_A(j) = \inf\{n > \tau_A(j-1), X_n \in A\}$, for $j \ge 2$, the successive return times to A, by $\mathbb{P}_A$ the probability measure on the underlying space such that $X_0 \in A$ and by $\mathbb{E}_A[\cdot]$ the $\mathbb{P}_A$-expectation. In the regenerative case, it results from the strong Markov property that the blocks of observations in between consecutive visits to the atom

$$\mathcal{B}_1 = (X_{\tau_A(1)+1}, \dots, X_{\tau_A(2)}), \dots, \mathcal{B}_j = (X_{\tau_A(j)+1}, \dots, X_{\tau_A(j+1)}), \dots \quad (15)$$

<figure><img src="image_8.png" /></figure>
form a collection of i.i.d. random variables, taking their values in the torus $T = \cup_{n=1}^{\infty} E^n$, and the sequence $\{\tau_A(j)\}_{j \ge 1}$, corresponding to successive times at which the chain forgets its past is a (possibly delayed) renewal process.

The class of regenerative chains includes a wide variety of Markov processes, including all chains with countable state spaces (where any recurrent state serves as an accessible atom). This class also incorporates many Markov models frequently employed in operations research and queueing theory (see, for example, Asmussen 2010). Further examples of (regular) regenerative chains are presented in Examples 2 and 3 below.

**Example 2** (BESSEL RANDOM WALKS) A Bessel random walk with drift $\delta \in [-1, +\infty)$ is a Markov chain with $N$ as state space, jumps in $\{-1, +1\}$, reflecting at 0 and with transition probabilities of the form:

$$ \Pi(0, 1) = 1 \quad \text{and} \quad 1 - \Pi(k, k-1) = \Pi(k, k+1) = \frac{1+h(k)-\delta/(2k)}{2} \quad \forall k \ge 1, $$

where $h(k) \in (-1 + \delta/(2k), 1 + \delta/(2k))$ and $h(k) = o(1/k)$ as $k \to +\infty$. It is recurrent when $\delta > -1$, positive recurrent when $\delta > 1$ and transient when $\delta = -1$. For $\delta = 1$, it is either recurrent or else transient, depending on the function $h(x)$. In the null-recurrent case, the chain is $\beta$-regular with $\beta = (1 + \delta)/2$, see Theorem 2.1 in Alexander (2011). Of course, when $\delta = 0$ and $h \equiv 0$, this chain corresponds to a simple reflected random walk with $p = 1/2$.

**Example 3** (NULL RECURRENT, NOT NECESSARILY REGULAR, CHAINS) By means of the model below, originally presented in Myklebust et al. (2012), one can generate $\beta$-null-recurrent chains for any $\beta > 0$, as well as null-recurrent chains that are not regular. Let $\{\eta_n\}_{n \in N}$ be a sequence of i.i.d. real-valued random variables. Consider the chain defined by:

$$ X_n = (X_{n-1} - 1)\mathbb{I}\{X_{n-1} > 1\} + \eta_n \mathbb{I}\{X_{n-1} \in [0, 1]\} \quad n \ge 1. $$

This chain is regenerative, with the interval $[0, 1]$ as atom. In addition, we have $\mathbb{P}_X(\tau_{[0,1]} > n) = \mathbb{P}(\lfloor \eta_1 \rfloor > n)$. Hence, $X$ is null-recurrent iff $\mathbb{E}[\lfloor \eta_1 \rfloor] = \infty$ and $\beta$-regular with $\beta \in (0, 1]$ iff the r.v. $\lfloor \eta_1 \rfloor$ has generalized discrete Pareto distribution with tail index $\beta$.

In the regenerative setting, all stochastic stability properties may be expressed in terms of speed of return to the atom. For instance, when $X$ is Harris recurrent, see Theorem 10.0.1 in Meyn et al. (2009), the invariant measure is equal to the occupation measure between two consecutive visits to the atom (up to a multiplicative factor): $\forall B \in \mathcal{E}, \mu(B) \sim \mathbb{E}_A[\sum_{i=1}^{\tau_A} \mathbb{I}\{X_i \in B\}]$. For instance, the chain is positive recurrent if and only if the expected return time to the atom is finite<fnref n="2" />, i.e. $\mathbb{E}_A[\tau_A] < \infty$, see Theorem 10.2.2 in Meyn et al. (2009). More generally, the $\beta$-regularity property can be characterized by the heaviness of the tail of the probability distribution of the regeneration times in the atomic case, as the following result shows.

<footnote n="2">² Its (unique) invariant probability distribution $\mu$ is then given by $\mu(B) = (1/\mathbb{E}_A[\tau_A])\mathbb{E}_A[\sum_{i=1}^{\tau_A} \mathbb{I}\{X_i \in B\}]$, for all $B \in \mathcal{E}$.</footnote>
**Proposition 3.1** (Karlsen and Tjostheim 2001, Theorem 3.1) Suppose that X is regenerative Harris recurrent. Let A be an atom for X and $\beta \in [0, 1]$. The following assertions are equivalent.

(i) The chain X is $\beta$-regular.

(ii) There exists a slowly varying function $L_A: \mathbb{R}_+ \to \mathbb{R}_+$ such that: $\forall m \ge 1$,

$$ \mathbb{P}_A(\tau_A \ge m) = L_A(m) \cdot m^{-\beta}. \quad (16) $$

As a direct consequence of Proposition 3.1, we have that if X is regenerative and $\beta$-regular, then $\beta = \sup\{\theta \in [0, 1] : \mathbb{E}_A[\tau_A^\theta] < +\infty\}$.

Based on the decomposition (15) of the whole sample path, limit theorems for regenerative Markov chains can be derived from the application of their i.i.d. counterparts to the sequence of blocks $(B_k)_{k \ge 1}$, see e.g. Meyn et al. (2009). This approach is usually referred to as the regenerative method and is extensively used to establish the asymptotic results stated in subsection 3.2 in the atomic case. Notice however that the regenerative blocks $B_1, \dots, B_{N_n}$, where $\Sigma_n(A) = N_n - 1$ denotes the (random) number of regenerations before time $n$, forming the truncated trajectory up to time $n$ are not independent (the sum of their length being less than $n$ in particular), which causes technical difficulties when establishing higher-order or non-asymptotic results, see e.g. Bertail and Clémençon (2006b) or Bertail and Clémençon (2004) and the references therein.

The inference technique for the regularity index $\beta$ of the chain X developed in subsection 3.3 is based on characterization (ii): the parameter $\beta$ is the tail index of a discrete generalized Pareto r.v., the regeneration time namely, i.e. the conditional survivor function of $\tau_A$ given $X_0 \in A$. Incidentally, notice that the parameter $\beta$ does not depend on the atom A considered (in contrast to the estimator analysed in subsection 3.3).

Based on a (random) number $N_n$ of (dependent) realizations of the regenerative time, namely

$$ S_j = \tau_A(j + 1) - \tau_A(j) \text{ for } j = 1, \dots, N_n, $$

one may naturally compute the estimator (2). As will be shown in Theorem 3.2, in spite of the dependence structure between the $S_j$'s, the consistency property is preserved in the Markovian framework.

## 3.2 Limit Theorems of the occupation times for Regular Markov Chains

We now recall the limit results related to the behaviour of the random occupation times $\Sigma_n(.)$ for regular Markov chains and discuss their limitations regarding their possible use to infer the regularity index $\beta$ with (asymptotic) guarantees. The latter essentially reveals that the empirical occupation measures $\Sigma_n(B)$ of Harris sets B grow at the sublinear rate $n^\beta$ (up to a slowly varying factor). As shall be seen, however, due to the great dispersion of their (asymptotic) distribution, the empirical occupation measures can hardly be used directly to estimate the key parameter $\beta$.

The result stated below claims that the logarithm of the occupation time of any Harris set provides a strongly consistent estimator of the regularity index $\beta$ when
appropriately normalized. It corresponds to the comment in Remark 3.7 of Karlsen and Tjostheim (2001).

**Proposition 3.2** Suppose that the chain X is β-regular with $\beta \in (0, 1)$ and B is a Harris set. Let v be its initial probability distribution. Then, we have:

$$ \tilde{\beta}_n(B) \to \beta \mathbb{P}_v - \text{a.s., as } n \to +\infty, \quad (17) $$

where

$$ \tilde{\beta}_n(B) = \frac{\ln \Sigma_n(B)}{\ln n}. \quad (18) $$

It was pointed out in Remark 3.7 of Karlsen and Tjostheim (2001) that this estimator is of limited practical use due to its slow rate of convergence, although no specific rate was given therein. Equation (4.2) in Chen (1999) may suggest that $|\tilde{\beta}_n(B) - \beta|$ is almost-surely $O(\ln K_{v,D}(n)/\ln n)$ as $n \to +\infty$. However, this has not been proven. In the following result, we obtain the limit distribution of the estimator and show that under the additional assumption that $L_{v,D}$ has a finite nonzero limit, the estimator has a logarithmic rate of convergence.

**Theorem 3.1** Suppose that the chain X is β-regular with $\beta \in (0, 1)$, B is a Harris set, and v is its initial probability distribution. Then, as $n \to +\infty$,

$$ \ln(n) \left( \tilde{\beta}_n(B) - \beta - \frac{\ln L_{v,D}}{\ln n} \right) \Rightarrow \ln \left( \frac{\mu(B)}{(Z_\beta)^\beta} \right) \text{ in } \mathbb{P}_v\text{-distribution.} $$

where $Z_\beta$ is a stable random variable with Laplace transform

$$ \psi_\beta(t) = \exp(-t^\beta / \Gamma(\beta + 1)), \quad t \ge 0. $$

In addition, if $\lim_{n \to +\infty} L_{v,D}(n)$ exists and is not 0, then, there exists a constant $\kappa > 0$ such that

$$ \ln(n) (\tilde{\beta}_n(B) - \beta) \Rightarrow \ln (\kappa/(Z_\beta)^\beta) \text{ in } \mathbb{P}_v\text{-distribution.} $$

The almost sure convergence of $\tilde{\beta}_n(B)$ towards $\beta$ suggests that, for $n$ large enough, $\ln \Sigma_n(B) \approx \beta \ln(n)$ and that the log-log plot of $\Sigma_n(B)$ and $n$ should look like a linear function with slope $\beta$, which could be possibly used to infer the value of $\beta$. Unfortunately, the dispersion of such a plot (and that of the process $\sigma_n(B)$, asymptotically described by Theorem D.1 of the Supplementary Material) is way too large in practice. To illustrate this, we simulated a simple symmetric random walk ($\beta = 0.5$) with $n = 10^5$ points, and we computed $\ln \Sigma_{\lfloor nt \rfloor}(B)$ for $0.1 \le t \le 1$ (choosing $B = \{0\}$, which is an atom for this regenerative regular chain). The outcomes of this simulation are depicted in Fig. 3.

This simulation illustrates in particular the slow convergence of $\tilde{\beta}_n(B)$ described in Corollary 3.1. Hence, in the regenerative case, it is more suitable to exploit the tail behaviour of the regenerative times (cf Proposition 3.1) to estimate the regularity index $\beta$, as shall be investigated in the next subsection.
<figure><img src="image_9.png" /><figcaption>Fig. 3 Log-log plot of ln Σ$_{[nt]}$({0}) in dark blue, the orange line representing the plot of the linear function x ↦→ 0.5 × x, while the blue area represents the 95% confidence interval for ln Σ$_{[nt]}$({0}) calculated from 100 independent trajectories. In this experiment n = 10⁵ and 0.1 ≤ t ≤ 1 (color figure online)</figcaption></figure>

3.3 Regularity Index of a Regular Chain—Statistical Inference

Assume that X is a regenerative regular chain with atom A and unknown regularity
index β ∈ (0, 1), and suppose that a sample path X₁, ..., n of length n ≥ 1 is
observed. Because the chain is Harris recurrent, the number of observed regeneration
times Σₙ(A) almost-surely tends to +∞ as n → +∞. Hence, with probability 1, we
have Nₙ ≥ 1 for n large enough and one can define

$$
\hat{p}_l^{(N_n)} = \frac{1}{N_n} \sum_{i=1}^{N_n} \mathbb{I}\{S_i > e^l\} \quad \text{for any } l \in \mathbb{R} \qquad (19)
$$

and form the statistic

$$
\hat{\beta}_{N_n}(k) = \ln \left( \hat{p}_k^{(N_n)} / \hat{p}_{k+1}^{(N_n)} \right), \qquad (20)
$$

provided that $\hat{p}_{k+1}^{(N_n)} > 0$. We point out that, due to the dependence between $N_n$ and the $S_i$'s, (19) is a biased (though strongly consistent) estimator of $p_l = \mathbb{P}_A(\tau_A > e^l)$ for any $l \in \mathbb{R}$. The estimator (20) is of the same form as (2) except that the number $N_n$ of observations is random, and it is not independent of the sequence $\{S_i\}$ of observations (in particular $S_1 + \dots + S_{N_n} \le n$).

Obtaining properties of this new estimator requires understanding the behaviour of $N_n$, to this, the key is the following result, proved in Karlsen and Tjostheim (2001), which gives the asymptotic distribution of $N_n$.
**Proposition 3.3** (Karlsen and Tjostheim 2001, Theorem 3.2) If X is an atomic β-regular Markov chain with initial probability distribution v, then

$$N_n \overline{F}(n) \Rightarrow \frac{M_\beta(1)}{\Gamma(1-\beta)} \text{ in } \mathbb{P}_v\text{-distribution.}$$

where $M_β(1)$ is a Mittag-Leffler distribution with parameter $β$.

From Proposition 3.3 we can deduce the following corollary which provides a deterministic control of $N_n$ with high probability.

**Corollary 3.1** For any $\delta \in (0, 1)$, let $H_\delta > R_\delta > 0$ be such that

$$\mathbb{P}\left(\frac{M_{\beta}(1)}{\Gamma(1-\beta)} \in [R_{\delta}, H_{\delta}]\right) \geq 1 - \delta/2.$$

Then, there exists $n_\delta$ such that

$$\forall n \geq n_{\delta} \quad \mathbb{P}_{v}\left(N_{n} \in\left[\frac{R_{\delta}}{\overline{F}(n)}, \frac{H_{\delta}}{\overline{F}(n)}\right]\right)>1-\delta.$$

Equipped with Corollary 3.1, we can provide the following result, which is the analogue of Proposition 2.1 in the Markovian case.

**Proposition 3.4** Let X be an atomic Markov chain, β-regular with $β \in (0, 1)$ and with initial probability distribution v. Let $\delta \in (0, 1/2)$, take $H_\delta$, $R_\delta$ and $n_\delta$ as in Corollary 3.1 and set $v_n(\delta) = \ln(36/\delta)\overline{F}(n)/H_\delta$. Then, as soon as $n \geq n_\delta$ and $p_{k+1} \geq (40H_\delta R_\delta^{-1})^2 v_n(\delta)$, we have with probability larger than $1 - 2\delta$:

$$|\hat{\beta}_{N_n}(k) - \beta| \leq \frac{60H_\delta}{R_\delta} \sqrt{\frac{v_n(\delta)}{p_{k+1}}} + \left|\ln\left(\frac{L(e^k)}{L(e^{k+1})}\right)\right|.$$

**Theorem 3.2** Suppose that the atomic chain X is β-regular with $β \in (0, 1)$. Let v be its initial probability distribution. If $k_n \to +\infty$ s.t. $(\ln n) \exp(k_n\beta)/n = o(L_A(\exp k_n))$ as $n \to +\infty$, then the estimator (20) is strongly consistent:

$$\hat{\beta}_{N_n}(k_{N_n}) \to \beta \quad \text{a.s as } n \to +\infty. \qquad (21)$$

In particular, strong consistency holds for $\hat{\beta}_{N_n}$ ($A \ln N_n$) with $A < 1$.

**Remark 4** (ON INVESTIGATING CONVERGENCE RATES) Due to the impossibility of tightly controlling the sequence $N_n$ by a deterministic quantity in probability and the non-linearity of the estimator, we have not been able to extend to the Markovian case the asymptotic normality results of Theorem 2.2 and Corollary 2.2 via Anscombe's theorem (Gut 2013, Theorem 7.3.2). Heuristically, if in Theorem 2.2 we take $k_n = \ln N_n$ and replace $N_n$ by its approximate expectation $n^\beta L_{v,A}(n)$ (Karlsen and Tjostheim 2001, Lemma 3.3), we would get a convergence rate of order $n^{-\beta(1-\beta)/2}L_1(n)$, where
$L_1(n)$ is the slowly varying function given by $\sqrt{L_{v,A}(n^\beta L_{v,A}(n))/L_{v,A}(n)^{1-\beta}}$. This suggests a convergence rate of order $n^{-\beta(1-\beta)/2}$ when $L_{v,A} \sim C > 0$. However, we have not been able to prove this claim.

**Remark 5 (TRAJECTORIES OF RANDOM LENGTH)** Suppose that the trajectory is observed until the $N$-th regeneration, i.e. $n = \tau_A(N)$, with $N \ge 2$. In this case, we will obtain a sequence of $N$ i.i.d blocks whose sizes follow the heavy-tailed distribution described in (1), and therefore, the results of Section 2.1 can be applied directly to this sequence. Notice that in this case, the total number of points observed in the chain (i.e. the amount of time we need to wait to collect the $N$ blocks) is a random variable that, while finite with probability one, has infinite expectation.

**Remark 6 (THE ATOMIC POSITIVE RECURRENT CASE)** When the chain is positive recurrent (or equivalently 1-regular), the estimator (20) can be naturally used to estimate the tail index $\beta' \ge 1$ of the regeneration time, when the latter has a regularly varying distribution. Dedicated theoretical results can be found in section C of the Supplementary Material.

## 3.4 Simulation Experiments

In order to analyse empirically the finite sample behaviour of our estimator in the Markovian case, we dedicate this subsection to a simulation example. As a generative model for our experiments, we use the Bessel random walks, defined in Example 2, as it will allow us to study the accuracy of our estimator, not just for different sample sizes, but also for different values of $\beta$.

For different values of $\beta$ and $n$, we have performed the following twofold experiment:

1. First, we have generated one path of a $\beta$-regular Bessel random walk of length $n$. We have then computed the estimator $\tilde{\beta}_{N_n}(k)$ based on this trajectory and plotted the results in order to determine the value $k_{\beta,n}$ where it first stabilizes. Using the same trajectory,

2. We have simulated 2000 independent paths of a $\beta$-regular Bessel random walk of length $n$. Then, we have computed the estimator $\tilde{\beta}_{N_n}(k_{\beta,n})$ based on each of these trajectories, as well as $\tilde{\beta}_n(\{0\})$ and the Hill estimator.

The results of this experiment are presented in Table 1. They show that the accuracy of our estimator greatly increases as the length $n$ gets larger (notice that the bias and variance are both divided by two in order of magnitude when $n$ increases from $10^3$ to $10^6$). The histograms shown in Figure 4 hint a possible asymptotic normality of the estimator $\tilde{\beta}_{N_n}$. The simulation study also suggests that our estimator outperforms the Hill estimator for this type of discrete datasets (as it has less bias and similar variance). Regarding the comparison with $\tilde{\beta}_n$, the experiments show that for moderately large sample sizes ($n \ge 10^5$) our estimator $\tilde{\beta}_{N_n}$ outperforms $\tilde{\beta}_n$ as it has less bias and comparable variance, while for smaller sample sizes $\tilde{\beta}_n$ works better. This behaviour can be attributed to the nature of our estimator, which focuses on estimating the probabilities of tail events $\{S > e^k\}$ to approximate the quantity $\beta(k)$ and use the fact

<figure><img src="image_10.png" /></figure>
<table><caption>**Table 1** Bias-Variance results of the simulation example</caption><thead><tr><th rowspan="2">β</th><th rowspan="2">n</th><th rowspan="2">k<sub>β,n</sub></th><th colspan="3">|Bias|<sup>2</sup></th><th colspan="3">Variance</th></tr><tr><th>β̂<sub>Nn</sub></th><th>β̂<sub>n</sub>(0)</th><th>Hill</th><th>β̂<sub>Nn</sub></th><th>β̂<sub>n</sub>(0)</th><th>Hill</th></tr></thead><tbody><tr><td rowspan="5">0.5</td><td>10<sup>3</sup></td><td>2</td><td>0.04738</td><td><strong>0.00690</strong></td><td>0.33076</td><td>0.12041</td><td><strong>0.02003</strong></td><td>0.74502</td></tr><tr><td>10<sup>4</sup></td><td>2</td><td>0.01649</td><td><strong>0.00429</strong></td><td>0.05514</td><td>0.05976</td><td><strong>0.01281</strong></td><td>0.21975</td></tr><tr><td>10<sup>5</sup></td><td>3</td><td><strong>0.00042</strong></td><td>0.00326</td><td>0.02041</td><td>0.03321</td><td><strong>0.00947</strong></td><td>0.09476</td></tr><tr><td>10<sup>6</sup></td><td>4</td><td><strong>0.00037</strong></td><td>0.00202</td><td>0.00629</td><td>0.02071</td><td><strong>0.00619</strong></td><td>0.03244</td></tr><tr><td>10<sup>7</sup></td><td>4</td><td><strong>0.00003</strong></td><td>0.00148</td><td>0.00216</td><td>0.00691</td><td><strong>0.00482</strong></td><td>0.03207</td></tr><tr><td rowspan="5">0.7</td><td>10<sup>3</sup></td><td>2</td><td>0.04194</td><td><strong>0.01796</strong></td><td>0.15864</td><td>0.12603</td><td><strong>0.01607</strong></td><td>0.55112</td></tr><tr><td>10<sup>4</sup></td><td>3</td><td><strong>0.00103</strong></td><td>0.01113</td><td>0.02915</td><td>0.05633</td><td><strong>0.01018</strong></td><td>0.10492</td></tr><tr><td>10<sup>5</sup></td><td>4</td><td><strong>0.00089</strong></td><td>0.00661</td><td>0.00634</td><td>0.02536</td><td><strong>0.00546</strong></td><td>0.02779</td></tr><tr><td>10<sup>6</sup></td><td>4</td><td><strong>0.00020</strong></td><td>0.00526</td><td>0.00283</td><td>0.00763</td><td><strong>0.00459</strong></td><td>0.09319</td></tr><tr><td>10<sup>7</sup></td><td>5</td><td><strong>0.00000</strong></td><td>0.00368</td><td>0.00040</td><td><strong>0.00292</strong></td><td>0.00319</td><td>0.00388</td></tr><tr><td rowspan="5">0.9</td><td>10<sup>3</sup></td><td>2</td><td>0.07861</td><td><strong>0.04123</strong></td><td>0.16550</td><td>0.15024</td><td><strong>0.00856</strong></td><td>0.50321</td></tr><tr><td>10<sup>4</sup></td><td>3</td><td><strong>0.00115</strong></td><td>0.02677</td><td>0.02007</td><td>0.04543</td><td><strong>0.00486</strong></td><td>0.09275</td></tr><tr><td>10<sup>5</sup></td><td>4</td><td><strong>0.00091</strong></td><td>0.01857</td><td>0.00329</td><td>0.01743</td><td><strong>0.00277</strong></td><td>0.01223</td></tr><tr><td>10<sup>6</sup></td><td>5</td><td><strong>0.00005</strong></td><td>0.01375</td><td>0.00105</td><td>0.00793</td><td><strong>0.00199</strong></td><td>0.00672</td></tr><tr><td>10<sup>7</sup></td><td>5</td><td><strong>0.00002</strong></td><td>0.01032</td><td>0.00068</td><td><strong>0.00095</strong></td><td>0.00137</td><td>0.00540</td></tr></tbody></table>

that, as $k$ tends to infinity, $\beta(k)$ converges to $\beta$ (see (3)). When $n$ is small, there are fewer sample points available in the tail, which forces us to select smaller values of $k$ to obtain accurate estimations of the tail probability, which worsens the estimation of $\beta$.

### 3.5 Perspectives—Extension to the Pseudo-regenerative Case

Harris chains are not necessarily regenerative, of course. However, the construction proposed in Nummelin (1978, 1984), referred to as the *Nummelin splitting technique*, permits to build a regenerative extension of any Harris chain (see section E of the Supplementary Material for a detailed description). Essentially, given a Harris recurrent Markov chain $X$, this technique allows the construction an atomic chain, called the “split chain” that has a very similar communication structure as the original chain. Moreover, if $X$ is $\beta$-regular, then the “split chain” is also $\beta$-regular with the same value of $\beta$ (Chen 1999, pp. 19).

A data-driven algorithm for obtaining samples of the split chain, given a finite sample of $X$ and the transition kernel density $\pi$, has been proposed in Bertail and Clémençon (2006a). We applied this algorithm to a simple random walk on $\mathbb{R}$, defined as $X_0 = 0$, $X_{n+1} = X_n + Z_n$ for $n \ge 1$, where $\{Z_n\}$ is a sequence of independent standard normal random variables.<fnref n="3" /> After constructing the split chain, we applied our

<footnote n="3">³ The transition kernel density of this Markov chain is $\pi(x, y) = f(y - x)$, where $f$ is the standard normal density function.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 4 Histograms and KDE for the estimators $\hat{\beta}_{N_n}$, $\tilde{\beta}_{N_n} (\{0\})$ and Hill for different values of $n$ and $\beta$</figcaption></figure>

estimator to it. The results, presented in Fig. 5a, show that when using the estimator $\hat{\beta}_{N_n}(k)$ on the split chain, the estimates of $\beta$ closely approximate the true value ($\beta = 1/2$) when $k$ is chosen near $\ln N_n$.

When the kernel density is unknown, a procedure to approximate the split chain, based on an estimation $\hat{\pi}_n$ of the kernel, was presented as Algorithm 3 in Bertail and Clémençon (2006a). As indicated in Theorem 3.1 of Bertail and Clémençon (2006b), the accuracy of this construction depends on the rate at which $\pi$ is estimated by $\hat{\pi}_n$. To our knowledge, in the null-recurrent case, the sole consistent estimator of the transition density documented in the literature is the Nadaraya–Watson estimator. The proof of the consistency can be found in section 5 of Karlsen and Tjostheim (2001). However, no results regarding its rate of convergence have been established so far. Moreover, the practical choice of the bandwidth parameter involved in this estimator is a difficult and largely unresolved problem, as discussed on pp. 412 in Karlsen and Tjostheim (2001). Hence, an ambitious line of further research consists in understanding how to implement practically the approximate regenerative block construction presented in Bertail and Clémençon (2006a) in order to extend the estimation methodology studied in the

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 5** Application of the pseudo-regeneration technique to construct the split chain and estimate $\beta$ using $\hat{\beta}_{N_n}(k)$ in the random walk $X_{n+1} = X_n + Z_n$ where $Z_n$ is an i.i.d. sequence of standard normal random variables. The total number of observed points in the chain is $10^6$ and the number of pseudo-blocks ($N_n$) is 418</figcaption></figure>

previous subsection to the regular pseudo-regenerative case with theoretical guarantees.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s11749-025-00975-9.

**Author Contributions** This paper is part of Carlos Fernández PhD thesis.

**Funding** Open access funding provided by Télécom Paris. This research has been conducted as part of the project Labex MME-DII (ANR11-LBX-0023-01). It was partially funded by a grant from the Hi! Paris foundation.

**Availability of data and materials** Not applicable.

**Code Availability** The code used to generate the results in the paper is available at https://github.com/carlosds731/tail_index_estimation.
**Declarations**

**Competing Interests** The authors have no competing interests to declare that are relevant to the content of this article.

**Ethics approval** Not applicable.

**Consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**References**

Adam E (July 2016) Persistence et vitesse d'extinction pour des modèles de populations stochastiques multitypes en temps discret. Université Paris Saclay (COmUE), Theses

Alexander K (2011) Excursions and local limit theorems for bessel-like random walks. Electronic Journal of Probability 16. https://doi.org/10.1214/ejp.v16-848

Asmussen S (2010) Applied Probability and Queues, 2nd edn. Stochastic Modelling and Applied Probability, Springer, New York

Bertail P, Clémentçon S (2004) Edgeworth expansions of suitably normalized sample mean statistics for atomic markov chains. Probability theory and related fields 130(3):388-414

Bertail, P., Clémentçon, S.: Regeneration-based statistics for Harris recurrent Markov chains. In: Bertail, P., Soulier, P., Doukhan, P. (eds.) Dependence in Probability and Statistics, pp. 3–54. Springer, New York (2006a). https://doi.org/10.1007/0-387-36062-X_1

Bertail P, Clémentçon S (2006b) Regenerative block bootstrap for markov chains. Bernoulli 12(4):689–712

Bingham NH, Goldie CM, Teugels JL (1987) Regular Variation. Encyclopedia of mathematics and its applications 27. Cambridge University Press, Cambridge

Chen X (1999) How often does a harris recurrent markov chain recur? The Annals of Probability 27. https://doi.org/10.1214/aop/1022677449

Chen X (2000) On the limit laws of the second order for additive functionals of harris recurrent markov chains. Probability Theory and Related Fields 116. https://doi.org/10.1007/pl00008724

Carpentier A, Kim AKH (2015) Adaptive and minimax optimal estimation of the tail coefficient. Statistica Sinica 25(3):1133–1144

Clauset A, Shalizi CR, Newman MEJ (2009) Power-law distributions in empirical data. SIAM Review 51. https://doi.org/10.1137/070710111

De Haan L, Ferreira A (2007) Extreme Value Theory: an Introduction, 1st edn. Springer, New York. https://doi.org/10.1007/0-387-34471-3

Goldstein ML, Morris SA, Yen GG (2004) Problems with fitting to the power-law distribution. The European Physical Journal B / Condensed Matter and Complex Systems 41. https://doi.org/10.1140/epjb/e2004-00316-5

Goldie CM, Smith RL (1987) Slow variation with remainder: theory and applications. Quart. J. Math. Oxford 38(1):45–71

Gao J, Tjøstheim D, Yin J (2013) Estimation in threshold autoregressive models with a stationary and a unit root regime. Journal of Econometrics 172(1):1–13

Gut A (2013) Probability: A Graduate Course, 2nd edn. Springer texts in statistics, Springer, New York
Hill BM (1975) A simple general approach to inference about the tail of a distribution. The Annals of Statistics 3. https://doi.org/10.2307/2958370

Kim M, Kokoszka P (2020) Consistency of the hill estimator for time series observed with measurement errors. Journal of Time series Analysis 41(3):421–435

Kim M, Kokoszka P (2023) Asymptotic and finite sample properties of hill-type estimators in the presence of errors in observations. Journal of Nonparametric Statistics 35(1):1–18

Klebaner FC (1993) Population-dependent branching processes with a threshold. Stochastic Processes and their Applications 46(1):115–127. https://doi.org/10.1016/0304-4149(93)90087-K

Karlsen HA, Myklebust T, Tjøstheim D (2010) Nonparametric regression estimation in a null recurrent time series. Journal of Statistical Planning and Inference 140. https://doi.org/10.1016/j.jspi.2010.04.029

Karlsen H, Tjøstheim D (2001) Nonparametric estimation in null recurrent time series. The Annals of Statistics 29. https://doi.org/10.1214/aos/1009210546

Lazzardi S, Valle F, Mazzolini A, Scialdone A, Caselle M, Osella M (2021) Emergent statistical laws in single-cell transcriptomic data. bioRxiv https://doi.org/10.1101/2021.06.16.448706

Myklebust T, Karlsen H, Tjøstheim D (2012) Null recurrent unit root process. Econometric Theory 28. https://doi.org/10.1017/S0266466611000119

Matsui M, Mikosch T, Tafakori L (2013) Estimation of the tail index for lattice-valued sequences. Extremes 16. https://doi.org/10.1007/s10687-012-0167-9

Manning CD, Schütze H (1999) Foundations of Statistical Natural Language Processing. The MIT Press, Cambridge, Massachusetts

Meyn S, Tweedie R, Glynn P (2009) Markov Chains and Stochastic Stability. Cambridge Mathematical Library, 2nd edn. Cambridge University Press, Cambridge

Nummelin E (1978) A splitting technique for harris recurrent markov chains. Probability Theory and Related Fields 43. https://doi.org/10.1007/bf00534764

Nummelin E (1984) General Irreducible Markov Chains and Non-Negative Operators. Cambridge Tracts in Mathematics 83. Cambridge University Press, Cambridge

Pakes AG (1971) On the critical galton-watson process with immigration. Journal of the Australian Mathematical Society 1971-nov vol. 12 iss. 4 12. https://doi.org/10.1017/s1446788700010375

Pickands J (1975) Statistical Inference Using Extreme Order Statistics. The Annals of Statistics 3(1):119–131. https://doi.org/10.1214/aos/1176343003

Resnick S (2007) Heavy-tail Phenomena: Probabilistic and Statistical Modeling. Springer, New York

Sidra A, Shougeng H, Nadeem AB (2018) Zipf's law and city size distribution: A survey of the literature and future research agenda. Physica A: Statistical Mechanics and Its Applications 492:75–92

Voitalov I, Hoorn P, Hofstad R, Krioukov D (2019) Scale-free networks well done. Physical Review Research 1(3):033034

Zanette DH (2006) Zipf's law and the creation of musical context. Musicae Scientiae 10(1):3–18. https://doi.org/10.1177/102986490601000101

Zubkov AM (1972) Life-periods of a branching process with immigration. Theory of Probability & Its Applications 17(1):174–183. https://doi.org/10.1137/1117018

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.