Alpha Estimation via Sample Splitting: A Two-Sample Framework for Stable-like Distributions

TL;DR

This paper introduces a novel semiparametric estimator for the stability index alpha in stable distributions, using a two-sample approach via random splitting to improve robustness and computational efficiency.

Contribution

It proposes a new alpha estimator leveraging sample splitting and empirical quantiles, avoiding complex likelihood calculations and enhancing robustness and efficiency.

Findings

Estimator is consistent and asymptotically normal.

Performs well in small samples and heavy-tailed scenarios.

Offers significant computational advantages over maximum likelihood methods.

Abstract

Stable distributions provide a flexible framework for modeling heavy-tailed and skewed data, with the stability index $\alpha$ quantifying tail heaviness. We propose a new semiparametric estimator for $\alpha$ that leverages the two-sum closure property of stable distributions within a location-scale framework. The method transforms a single sample into two pseudo-independent samples via repeated random splitting and estimates $\alpha$ using weighted least squares applied to empirical quantiles. This approach avoids intractable likelihood calculations, offers computational advantages over maximum likelihood estimation, and remains robust to skewness. We establish consistency and asymptotic properties of the estimator and assess its finite-sample performance via simulation. Results indicate competitive accuracy, particularly in small samples and heavy-tailed settings, with substantial computational savings.

Tables

Table 1. Table 1: Bias and RMSE of the split–sample estimates when ( α , β ) = ( 1 , 0 ) 𝛼 𝛽 1 0 (\alpha,\beta)=(1,0)

		${\hat{\alpha }}_1$			${\hat{\alpha }}_2$			${\hat{\alpha }}_3$
		Bias	RMSE		Bias	RMSE		Bias	RMSE
$n+2m=150$	$B=1$	0.102	0.515		0.102	0.515		0.102	0.515
	$B=10$	-0.036	0.201		0.090	0.209		0.005	0.223
	$B=100$	-0.064	0.148		0.091	0.153		-0.025	0.139
	$B=250$	-0.065	0.145		0.092	0.149		-0.026	0.134
	$B=500$	-0.066	0.144		0.092	0.148		-0.027	0.132
$n+2m=300$	$B=1$	0.105	0.398		0.105	.398		0.105	0.398
	$B=10$	-0.010	0.131		0.097	0.174		0.020	0.156
	$B=100$	-0.025	0.090		0.095	0.131		0.003	0.092
	$B=250$	-0.026	0.087		0.094	0.127		0.002	0.087
	$B=500$	-0.026	0.086		0.094	0.126		0.001	0.085
$n+2m=600$	$B=1$	0.041	0.253		0.041	0.253		0.041	0.253
	$B=10$	-0.009	0.086		0.049	0.110		0.006	0.110
	$B=100$	-0.013	0.061		0.051	0.082		0.001	0.063
	$B=250$	-0.012	0.058		0.052	0.080		0.001	0.059
	$B=500$	-0.012	0.057		0.052	0.079		0.001	0.057

Table 2. Table 2: Bias and RMSE of the split–sample estimates when ( α , β ) = ( 1.95 , 0 ) 𝛼 𝛽 1.95 0 (\alpha,\beta)=(1.95,0)

		${\hat{\alpha }}_1$			${\hat{\alpha }}_2$			${\hat{\alpha }}_3$
		Bias	RMSE		Bias	RMSE		Bias	RMSE
$n+2m=150$	$B=1$	-0.285	0.520		-0.285	0.520		-0.285	0.520
	$B=10$	-0.127	0.234		-0.285	0.321		-0.148	0.251
	$B=100$	-0.092	0.155		-0.281	0.289		-0.084	0.154
	$B=250$	-0.088	0.147		-0.281	0.286		-0.076	0.140
	$B=500$	-0.087	0.143		-0.280	0.285		-0.074	0.133
$n+2m=300$	$B=1$	-0.174	0.338		-0.174	0.338		-0.174	0.338
	$B=10$	-0.076	0.162		-0.173	0.200		-0.075	0.162
	$B=100$	-0.046	0.100		-0.173	0.180		-0.027	0.091
	$B=250$	-0.044	0.094		-0.174	0.179		-0.021	0.081
	$B=500$	-0.043	0.092		-0.174	0.179		-0.019	0.077
$n+2m=600$	$B=1$	-0.129	0.254		-0.129	0.254		-0.129	0.254
	$B=10$	-0.044	0.113		-0.122	0.143		-0.050	0.121
	$B=100$	-0.027	0.072		-0.124	0.129		-0.015	0.069
	$B=250$	-0.024	0.066		-0.124	0.128		-0.011	0.061
	$B=500$	-0.024	0.065		-0.124	0.128		-0.011	0.059

Table 3. Table 3: Average percentage (%) of times α b ^ ^ subscript 𝛼 𝑏 {\hat{\alpha_{b}}} is on the boundary of the parameter space for B = 250 𝐵 250 B=250 .

	$n+2m=150$			$n+2m=300$			$n+2m=600$
	${\hat{\alpha }}_b=0$	${\hat{\alpha }}_b=2$		${\hat{\alpha }}_b=0$	${\hat{\alpha }}_b=2$		${\hat{\alpha }}_b=0$	${\hat{\alpha }}_b=2$
$\alpha =1.0$	2.67	12.67		0.15	5.81		0.00	1.10
$\alpha =1.5$	2.04	27.93		0.15	22.43		0.00	13.90
$\alpha =1.95$	2.65	45.53		0.27	47.30		0.00	45.97

Table 4. Table 4: RMSE for different estimators for α = 1 𝛼 1 \alpha=1 and β ∈ { 0 , 0.75 } 𝛽 0 0.75 \beta\in\{0,0.75\}

	$n=150$			$n=300$			$n=600$
	$\beta =0$	$\beta =0.75$		$\beta =0$	$\beta =0.75$		$\beta =0$	$\beta =0.75$
MLE	0.089	0.090		0.059	0.064		0.042	0.045
MQE	0.124	0.176		0.083	0.121		0.060	0.088
SSE3	0.127	0.157		0.090	0.110		0.065	0.074
SSE9	0.144	0.127		0.087	0.081		0.058	0.053
SSE19	N/A	N/A		0.121	0.112		0.062	0.057

Table 5. Table 5: Ratio of SSE 9 RMSE over MQE RMSE for n = 600 𝑛 600 n=600 and different choices of ( α , β ) 𝛼 𝛽 (\alpha,\beta)

	Monte Carlo Relative Efficiency
$\alpha$	$0.5$	$0.75$	$1.0$	$1.25$	$1.5$	$1.75$	$2.0$
$\beta =0$	$0.467$	$0.944$	$0.976$	$1.005$	$0.878$	$0.684$	$0.627$
$\beta =0.75$	$1.141$	$0.654$	$0.601$	$0.656$	$0.687$	$0.663$	$0.627$

Equations

$${\psi }_X(t)=E(e^{itX})=\{\begin{array}{ll}\text{exp(−γα∣t∣α[1+iβtan (2πα)(sign t)(∣γt∣1−α−1)]+iδt)}& \text{ if α=1}\\ \text{exp(−γ∣t∣[1+iβ(π2)(sign t)(log γ∣t∣)]+iδt)}& \text{ if α=1}\end{array}$$

$$Y\sim \text{S}(\alpha ,\beta ,\sigma \gamma ,\mu +\sigma \delta ).$$

$$Y\sim \text{S}(\alpha ,\beta ,2^{1/\alpha }\gamma ,{\delta }^{\prime })$$

$${\delta }^{\prime }=\{\begin{array}{ll}2\delta +(\text{tan 2πα})\beta \gamma [2^{1/\alpha }-2]& \text{if α=1}\\ 2\delta +\frac{2}{\pi }\beta \gamma [2^{1/\alpha }\text{log (21/αγ)}-2\text{log γ}]& \text{if α=1}.\end{array}$$

$$\alpha =\frac{\log 2}{\log \sigma }.$$

$$\tilde{F}(X_{(j)})=\frac{j-1}{n-1},j=1,2,\dots ,n,$$

$$G^{-1}(t)=\mu +\sigma F^{-1}(t),0<t<1.$$

$${\mathrm{W}}_1(\theta )$$

$${\mathrm{W}}_2(\theta )$$

$$\mathrm{V}(\theta )=[{{\mathrm{W}}_1(\theta )}^{\top },{{\mathrm{W}}_2(\theta )}^{\top }{]}^{\top }$$

$$\Omega =\left(\begin{array}{cc}\lambda \Sigma & 0\\ 0& (1-\lambda )\Sigma \end{array}\right)$$

$$\tilde{Q}(\theta )={\text{V}(\theta )}^{\top }{\tilde{\Omega }}^{-1}\text{V}(\theta )$$

$$\hat{\alpha }=\{\begin{array}{ll}0& \text{if σ^∈(0,1)}\\ 2& \text{if σ^∈[1,2)}\\ \frac{\log 2}{\log \hat{\sigma }}& \text{if σ^∈[2,∞)}.\end{array}$$

$${\theta }_{\mathcal{K}}=t\left(X_{k_1},\dots ,X_{k_n},X_{k_{n+1}}+X_{k_{n+2}},\dots ,X_{k_{n+2m-1}}+X_{k_{n+2m}}\right)$$

$$\bar{\theta }=\frac{1}{K}\sum _{k=1}^K{\theta }_{{\mathcal{K}}_k}$$

$${\hat{\alpha }}_1=\{\begin{array}{ll}0& \text{if σˉ∈(0,1)}\\ 2& \text{if σˉ∈[1,2)}\\ \log 2/\log \bar{\sigma }& \text{if σˉ∈[2,∞)}.\end{array}$$

$${\left(\frac{mn}{m+n}\right)}^{1/2}\left(\hat{\alpha }-\alpha \right)\stackrel{d}{\to }N\left(0,\frac{{\alpha }^4}{(log2{)}^2}{\Gamma }_{22}\right)$$

Full text

Some Key Words: Stable Distributions, Tail Index, Characteristic Exponent, Location-Scale Model, Split-Sample Estimation, Data Permutation.

1 Introduction

The stable probability law was introduced by Paul Lévy in 1924 in his work on sums of independent and identically distributed (iid) random variables. Originating from his attempt to generalize the central limit theorem, the class of stable distributions is defined as follows: A non-degenerate random variable $X$ is said to have a stable distribution if and only if for all $n\geq 2$, there exist constants $a_{n}\in\mathbb{R}$ and $b_{n}>0$ such that $X_{1}+X_{2}+\ldots+X_{n}\stackrel{{\scriptstyle d}}{{=}}a_{n}+b_{n}X$. Here $X_{1},X_{2},\ldots,X_{n}$ are independent copies of $X$ and the symbol $\stackrel{{\scriptstyle d}}{{=}}$ is used to denote equality in distribution. The random variable $X$ is called strictly stable if and only if $a_{n}=0$ for all values of $n$. See Feller (2008) and Nolan (2013), for a comprehensive overview of existing results on stable distributions.

Stable distributions are typically described in terms of their characteristic functions. The random variable $X$ is said to have a stable distribution $\text{S}(\alpha,\beta,\gamma,\delta)$ if the characteristic function of $X$, for all real $t$, is given by

[TABLE]

where $\alpha\in(0,2]$ is commonly referred to as the stability index, $\beta\in[-1,1]$ is a skewness parameter, $\gamma>0$ is a scale parameter and $\delta\in\mathbb{R}$ is a location parameter. The index parameter $\alpha$, also called the stability index or characteristic exponent, measures the heaviness of the tails of the distribution. As $\alpha$ decreases, the tail heaviness of the distribution increases. For $\alpha<=1$, the mean of the distribution does not exist, while for $1<\alpha<2$,  the variance of the distribution does not exist either. In general, a stable random variable with index parameter $\alpha\in(0,2)$ possesses absolute moments of order $p$ where $p<\alpha$; that is, $E(|X|^{p})<\infty$ for $p<\alpha$.

A sum of iid stable random variables with common characteristic exponent $\alpha$ is again stable, retaining the characteristic exponent of the original distribution. This property is termed stability. From a practical point of view, stable distributions are an attractive option for modeling data that exhibit heavy tails and skewness as these features can easily capture by stable distributions. With three exceptions discussed below, stable distributions do not have closed-form density functions. However, all non–degenerate stable distributions are continuous with infinitely differentiable density functions. The three special cases in which there exists a closed-form expression for the density function are the normal, Cauchy and Lévy distributions. When setting $\alpha=2$, (1) corresponds to the characteristic function of a Normal$(\delta,2\gamma^{2})$ distribution. Similarly, upon setting $(\alpha,\beta)=(1,0)$, the distribution is Cauchy$(\gamma,\delta)$ while upon setting $(\alpha,\beta)=(1/2,1)$, the resulting distribution is Lévy with $\gamma$ and $\delta+\gamma$ being the scale and location parameters respectively.

The lack of a closed-form density function together with the non-existence of moments has made parameter estimation a historically challenging task. While some applications require the estimation of all four parameters, in many instances the parameter of greatest interest is the stability index $\alpha$ which determines the tail heaviness of the distribution. This paper focuses only on the estimation of the parameter $\alpha$. Many methods have been proposed to estimate $\alpha$, all of which fall in to three categories: maximum likelihood, quantile methods and characteristic function methods. Quantile methods include both methods based on sample quantiles and methods based on extreme order statistics.

DuMouchel (1973, 1975) did extensive work on using a maximum likelihood type method to estimate the parameters of a stable distribution. His method relied on grouping the data into bins and numerically maximizing an approximate log–likelihood function. Mittnik et al. (1999) used the fast Fourier transform to estimate the parameters, while Nolan (2001) developed routines for numerical computation of the integrals involved in the Fisher information matrix which can then also be used for maximum likelihood estimation. Lombardi (2007) proposed a MCMC method to estimate the stable parameters. Buckle (1995) developed Bayesian methodology for inference in stable distributions, while Peters et al. (2012) proposed likelihood-free Bayesian inference for stable models. Still, direct maximization of the likelihood function presents many challenges in practice and it is therefore not a popular approach when estimating the tail index.

Fama & Roll (1968, 1971) did early work on estimating $\alpha$ using order statistics, but their method only applied to symmetric stable distributions for $\alpha\in[1,2]$. McCulloch (1986) developed the quantile method now popular in application which works for both symmetric and skew stable distributions and $\alpha\in[0.5,2]$. The Hill estimator, see Hill (1975), is a popular measure of tail heaviness for distributions with Pareto-like tails, and can also be adapted to estimate $\alpha$. However, the Hill estimator tends to have large bias in small to moderately sized samples.

Authors that have considered estimators making direct use the characteristic function include Press (1972), Paulson et al. (1975), Koutrouvelis (1980) and Brockwell & Brown (1981). These methods have good performance properties, but some still shy away from models that rely on inference in the complex domain.

In this paper, the problem of estimating $\alpha$ is viewed through the lens of representing stable random variables as members of a location–scale family of distributions. The general framework presented relies on repeatedly partitioning the observed data into two independent samples, say $X_{1},\ldots,X_{n}$ and $Y_{1},\ldots,Y_{m}$. The $X$ sample is created by selecting values from the observed data without replacement. Thereafter, the remaining observations are randomly paired and the $Y$ sample is created by adding the data values in each pair. This “split-sample” approach allows the estimation of $\alpha$ to be treated as a two-sample estimation problem. In particular, a quantile–based method of estimation is proposed in this paper. This quantile method is a variation of a method due to Potgieter & Lombard (2012) who consider estimating location and scale parameters for two independent samples belonging to the same location-scale family.

In Section 2, the location–scale representation of stable distributions is discussed and a connection between the scale parameter in the two-sample setting and the stability index $\alpha$ of a single sample is established. In Section 3, the split-sample estimator is formally defined. Additionally, the two-sample quantile method of Potgieter and Lombard is reviewed and extended to the present framework for estimating $\alpha$. Section 4 presents results from a simulation study carried out to investigate the number of sample partitions needed to give good RMSE performance of the proposed estimator. Section 5 deals with some practical considerations, such as choosing the number of quantiles to use in the Potgieter & Lombard method. Section 6 compares the proposed estimator to other existing methods and recommendations for implementing the method are discussed in Section 7.

2 Stable Distribution Location–Scale Representation

Suppose that $X$ and $X^{\prime}$ are iid $\text{S}(\alpha,\beta,\gamma,\delta)$ random variables and define $Y=X+X^{\prime}$. As the sum of two stable random variables with common parameters $\alpha$ and $\beta$ is again stable, it follows that there are constants $\mu\in\mathbb{R}$ and $\sigma>0$ such that $Y\stackrel{{\scriptstyle d}}{{=}}\mu+\sigma X$ where $\stackrel{{\scriptstyle d}}{{=}}$ indicates equality in distribution. Thus, the distribution of $Y$ has two equivalent representations,

[TABLE]

and

[TABLE]

where $\delta^{\prime}$ is given by

[TABLE]

Here, relation (3) is derived as the characteristic function of the sum of two iid stable variables, $X+X^{\prime}$, while (2) is the characteristic function of the random variable $\mu+\sigma X$ and follows from the properties of a location-scale transformation. As (2) and (3) are equivalent, the parameters in the two cases must also be equal. Specifically, by equating the scale parameters in the two formulations of $Y$, it follows that $\sigma\gamma=2^{1/\alpha}\gamma$. Solving for $\alpha$ gives

[TABLE]

Therefore, the problem of estimating $\alpha$ is equivalent to that of estimating $\sigma$, the scale parameter relating $X$ and $Y$. This relation forms the basis of the estimation procedure proposed in this paper.

Now, let $X$ and $Y$ be two independent random variables such that $Y\stackrel{{\scriptstyle d}}{{=}}\mu+\sigma X$ for appropriate constants $\mu$ and $\sigma$. For the time being, assume that independent samples $X_{1},\ldots,X_{n}$ and $Y_{1},\ldots,Y_{m}$ are observed. Potgieter & Lombard (2012) proposed a nonparametric method called asymptotic likelihood (AL) for estimating $\mu$ and $\sigma$ from the two independent samples. Their method only assumes that the random variables $X$ and $Y$ have continuous and strictly increasing distribution functions with differentiable density functions. Despite the general lack of closed-form expressions for the density and distribution functions, stable distributions do satisfy these assumptions. The AL method can therefore be applied to the stable setting to estimate $\sigma$, and subsequently $\alpha$, from the independent $X$- and $Y$-samples. The question of obtaining these samples is further addressed in Section 3, while the remainder of this section gives a brief overview of the implementation of the AL method in the present setting.

Let $F$ and $G$ (correspondingly $f$ and $g$) denote the respective distribution (density) functions of random variables $X$ and $Y$. For iid random samples $X_{1},X_{2},\ldots,X_{n}$ and $Y_{1},Y_{2},\ldots,Y_{m}$, denote the respective order statistics by $X_{(1)}<X_{(2)}<\ldots<X_{(n)}$ and $Y_{(1)}<Y_{(2)}<\ldots<Y_{(m)}$. Let $\widetilde{F}$ and $\widetilde{G}$ denote, respectively, empirical distribution functions of $X$ and $Y$ made continuous using linear interpolation. That is, let $\widetilde{F}(x)=0$ for $x<X_{(1)}$, $\widetilde{F}(x)=1$ for $x>X_{(n)}$,

[TABLE]

and for $x\in(X_{(j)},X_{(j+1)})$, define $\widetilde{F}(x)$ to be the interpolated value between the pairs $(X_{(j)},(j-1)/(n-1))$ and $(X_{(j+1)},j/(n-1))$, $j=1,\ldots,n-1$. A similar definition holds for $\widetilde{G}(x)$. The continuous empirical quantile functions $\widetilde{F}^{-1}$ and $\widetilde{G}^{-1}$ are uniquely defined by the relation $\widetilde{F}[\widetilde{F}^{-1}(t)]=\widetilde{G}[\widetilde{G}^{-1}(t)]=t$.

Now, when the relation $Y\stackrel{{\scriptstyle d}}{{=}}\mu+\sigma X$ holds, the quantile functions $F^{-1}$ and $G^{-1}$ satisfy

[TABLE]

Define $\varepsilon_{n}(t)=\widetilde{F}^{-1}(t)-F^{-1}(t)$ and $\delta_{m}(t)=\widetilde{G}^{-1}(t)-\mu-\sigma F^{-1}(t)$. Potgieter & Lombard (2012) show that for fixed $0<t_{1}<\ldots<t_{k}<1$, the independent random vectors $\left[n^{1/2}f(F^{-1}(t_{j}))\varepsilon_{n}(t_{j})\right]^{\top}$, $j=1,\ldots,k$ and $\left[m^{1/2}g[G^{-1}(t_{j})]\delta_{m}(t_{j})\right]^{\top}$, $j=1,\ldots,k$ converge in distribution to multivariate normal distributions with common covariance matrix $\Sigma_{ij}=\text{min}(t_{i},t_{j})-t_{i}t_{j}$ as $\min(n,m)\rightarrow\infty$.

Now, let $\phi_{j}=F^{-1}(t_{j})$, $j=1,\ldots,k$ and define parameter vector $\theta=[\phi_{1},\phi_{2},\ldots,\phi_{k},\mu,\sigma]^{\top}$. The $k+2$ parameters in $\theta$ can be estimated using the established asymptotic normality. Define vectors

[TABLE]

and

[TABLE]

where $\widetilde{f}$ and $\widetilde{g}$ are kernel density estimators of $f$ and $g$ respectively. It then follows that as $\min(n,m)\rightarrow\infty$, $(nm/(n+m))^{1/2}\text{V}(\theta)$ converges in distribution to a $2k$ multivariate normal distribution with zero mean and covariance matrix $\Omega$ given by

[TABLE]

where $\lambda=\lim_{m,n\rightarrow\infty}n/(n+m)$ and the component of the asymptotic log-likelihood of $\text{V}(\theta)$ involving the $k+2$ parameters is $Q(\theta)={\text{V}(\theta)}^{\top}\Omega^{-1}\text{V}(\theta)/2$. Now define,

[TABLE]

where $\widetilde{\Omega}$ is $\Omega$ with $\lambda$ replaced by $\widetilde{\lambda}=n/(n+m)$. The estimator $\hat{\theta}$ that minimizes $\widetilde{Q}(\theta)$, cannot expressed in closed form but can be easily found using standard numerical optimization routines. The component estimators $\hat{\mu}$ and $\hat{\sigma}$ are called the AL estimators of $\mu$ and $\sigma$.

3 Split–sample estimator

The method for estimating $\sigma$ outlined in Section 2 assumes the availability of two independent samples satisfying the specified location-scale relationship, while only the equivalent of an $X$-sample is observed in practice. Suppose this sample consists of $n+2m$ observations from a stable distribution with unknown index parameter $\alpha$. It is possible to create two independent samples from the single observed sample using the following method: First select $n$ observations randomly and treat them as the $X$-sample. Next, form randomly $m$ pairs from the remaining $2m$ observations and sum the observations in each pair. Treat these $m$ sums as the $Y$-sample. By the properties of stable distributions, the relation $Y\stackrel{{\scriptstyle d}}{{=}}\mu+\sigma X$ holds for these constructed samples. Additionally, the $X$ and $Y$ samples constructed in this way are independent and therefore the AL method proposed by Potgieter & Lombard can be used to estimate $\sigma$ and then, subsequently, $\alpha$ using (4). It should be noted that although the method guarantees $\hat{\sigma}>0$, the estimator resulting from applying (4) to $\hat{\sigma}$ may not be in the interval $[0,2]$. It is therefore reasonable to define

[TABLE]

The perceived discontinuity in the definition of $\hat{\alpha}$ results from a discontinuity in $1/\log{\sigma}$ when $\sigma=1$. Since this method involves splitting the sample, we referred to the estimator $\hat{\alpha}$ as the split–sample estimator (SSE) of $\alpha$.

The estimator proposed above, of course, uses only one random permutation of the data. Ideally, all possible sample permutations would be constructed, and each permutation would be used to construct an estimate of $\sigma$ and/or $\alpha$. Finally, these estimators would then be combined to create some ensemble estimator of $\alpha$. Specifically, let $t(a_{1},\ldots,a_{n},a_{n+1},a_{n+m})$ be a function that is permutation-invariant in the first $n$ arguments $a_{1},\ldots,a_{n}$ and also permutation-invariant in the last $m$ arguments $a_{n+1},\ldots,a_{n+m}$. For $\mathcal{K}=(k_{1},\ldots,k_{n+2m})$ a random permutation of the integers $1,\ldots,n+2m$, define

[TABLE]

to be the statistic $t$ calculated for the permutation $\mathcal{K}$. Letting $K$ denote the total number of possible permutations of the data, define

[TABLE]

where $\theta_{\mathcal{K}_{k}}$ denotes the value of statistic $t$ for the kth permutation of the data. Note that the effect here is that of creating a $U$-statistic with a symmetric kernel. Whereas $\theta_{\mathcal{K}}$ is not symmetric in its arguments, $\bar{\theta}$ is such. In the present setting, the statistic $t$ evaluated for one permutation represents an AL estimate of $\sigma$ based on a single permutation of the data, whereas the ideal is to calculate $\bar{\theta}$ by evaluating $t$ for all possible data permutations. However this is not realistic, as there are $\binom{n+2m}{n}\binom{2m}{2}\binom{2m-2}{2}\ldots\binom{4}{2}=\binom{n+2m}{n}\prod\limits_{i=0}^{m-1}\binom{2m-2i}{2}$ unique possible $X$- and $Y$-samples that can be created using the proposed data-splitting method. Even when $n$ and $m$ are only moderately large, this constitutes too large a number of sample permutations to practically evaluate all of them. For example, considering the scenario $m=n$ with sample size $n+2m=30$, there are approximate $7.08\times 10^{22}$ such sample permutations and when $n+2m=150$, there are approximately $1.67\times 10^{183}$ such sample permutations. The number of possible sample permutations grows at a super-exponential rate.

Of course, the inability to evaluate all possible data permutations should not steer one towards the other extreme where only a single data permutation is used to estimate $\sigma$, as a big loss in efficiency could result. A compromise is proposed, in that data permutation process creating $X$ and $Y$ samples is repeated $B$ where $B$ is some ”large” integer and these $B$ estimates of $\sigma$ can then be combined in an appropriate manner to estimate $\alpha$.

Let $\hat{\sigma}_{1},\ldots,\hat{\sigma}_{B}$ denote the estimates of $\sigma$ resulting from randomly splitting the sample $B$ times. The question of how to combine these to create an estimate $\hat{\alpha}$ is now considered. Proposed here are three ways of combining the $B$ estimates to find an estimate of $\alpha$:

1. (i)

Define $\bar{\sigma}=\left(\sum_{j=1}^{B}\hat{\sigma}_{j}\right)/B$. That is, $\bar{\sigma}$ is the average of the $B$ values $\hat{\sigma}_{1},\ldots,\hat{\sigma}_{B}$. Using this, define estimator $\hat{\alpha}_{1}$,

[TABLE] 2. (ii)

Let $\hat{\alpha}_{j}$ denote the value obtained after applying transformation (7) to the $\hat{\sigma}_{j}$ for $\,j=1,\ldots,B$. Here, $\hat{\alpha}_{j}$ denotes the estimate of $\alpha$ for the $j^{\text{th}}$ random split of the sample. Define estimator $\hat{\alpha}_{2}=\bar{\alpha}=\left(\sum_{j=1}^{B}\hat{\alpha}_{j}\right)/B$. 3. (iii)

Estimate $\alpha$ using $\hat{\alpha}_{3}=\text{median}\left(\hat{\alpha}_{1},\ldots,\hat{\alpha}_{B}\right)$. That is, instead of taking the average as in (ii), the median of the $B$ estimates of $\alpha$ is evaluated.

Here, estimators $\hat{\alpha}_{1}$ and $\hat{\alpha}_{2}$ fall within the outlined framework of approximating a $U$-statistic with symmetric kernel. On the other hand, $\hat{\alpha}_{3}$ is outside this framework, but is included as a robust alternative. The performance of the three estimators, as well as the number splits $B$ to be used, are investigated in the simulation study presented in the next section.

4 Simulation Study

The accuracy of the estimators defined in the previous section will depend on $B$, the number of random splits used. A too small $B$ results in an estimate with large variability, while a very large $B$ detracts from the practical viability of the approach due to computational cost. Therefore, a good choice of $B$ is essential. A simulation study has been performed to assess how the choice of $B$ affects the defined estimators. In the simulation, samples of size $n+2m\in\left\{150,300,600\right\}$ were drawn from standardized stable distributions ($\gamma=1$ and $\delta=0$) for values $\alpha\in\{0.5,1,1.5,1.95\}$ and $\beta\in\left\{0,0.75\right\}$. The samples were split such that $m=n$ when evaluating the AL estimator. In addition, the AL method was implemented using $k=9$ equally spaced $t$-values, $t=[1/10,\ldots,9/10]$. The three estimators $\hat{\alpha}_{j}$, $j=1,2,3$ were evaluated for $B\in\left\{1,10,50,100,500\right\}$. A total of $N=1000$ random samples were drawn for each configuration of $(\alpha,\beta,n)$. The mean value of the estimate and the Monte Carlo RMSE was evaluated for each configuration.

Tables 1 and 2 About Here

Table 1 presents results for the case $(\alpha,\beta)=(1,0)$ and Table 2 presents results for $(\alpha,\beta)=(1.95,0)$. Generally, as the number of sample splits $B$ increases, the both the bias and the RMSE tend to decrease. Generally there is a steady decrease in RMSE when going from $B=1$ to $B=100$, but there is only a small decrease in RMSE when going from $B=100$ to $B=250$, while going from $B=250$ to $B=500$, the decrease seems to be negligible. Simulation results for parameter configurations not presented here all show a similar pattern to those seen in Tables 1 and 2. It is clear that RMSE continues to decrease as $B$ increases, but that this reduction diminishes for $B\geq 250$. To illustrate, in the context of Table 2, when further increasing $B$ to $2500$, the RMSE of $\hat{\alpha}_{1}$ decreases from $0.065$ at $B=500$ to $0.063$ at $B=2500$ and the RMSE of $\hat{\alpha}_{3}$ similarly decreases from $0.059$ to $0.057$. For practical purposes, a recommendation is made to use $B=250$. This value ensures fast computation, but already shows good performance. A practitioner who wanted to see any further improvement in RMSE would have to choose $B$ a whole order of magnitude larger.

When comparing the estimators $\hat{\alpha}_{j}$, $j=1,2,3$, it should be noted that $\hat{\alpha}_{2}$ consistently performs much worse (in terms of bias and RMSE) than both $\hat{\alpha}_{1}$ and $\hat{\alpha}_{3}$. This can be explained, at least in part, by the truncation that occurs when applying transformation (7) from the $\sigma$-scale to the $\alpha$-scale. When a large proportion of the $\hat{\sigma}_{b}$ are outside the interval $[\sqrt{2},\infty)$, that same proportion of $\hat{\alpha}_{b}$ are on the boundaries ([math] and $2$) of the parameter space. As the estimator $\hat{\alpha}_{2}$ is calculated by averaging on the $\alpha$-scale, a large proportion of boundary values can increase the bias of the estimate. On the other hand, $\hat{\alpha}_{1}$ is calculated by averaging on the $\sigma$-scale and the truncation only comes into play when $\bar{\sigma}$ is outside $[\sqrt{2},\infty)$. Similarly, since $\hat{\alpha}_{3}$ is the median of the $B$ split-sample estimates, the truncation only comes into play if more than $50\%$ of the values are truncated to a specific boundary.

To illustrate the occurrence of boundary values, Table 3 reports the results of a simulation study in which $N=2000$ samples were generated from a standard stable distribution with $\alpha\in\{1,1.5,1.95\}$ and $\beta=0$ with $n+2m\in\left\{150,300,600\right\}$. For each simulated set of data, values $\hat{\alpha_{b}}$ were calculated from $B=250$ random splits. The table reports the average percentage of $\hat{\alpha_{b}}$ that were truncated to either [math] or $2$.

Table 3 About Here

The content of Table 3 is unsurprising. As the sample size increases, the occurrence of truncation to the boundaries decreases. The one exception is when $\alpha=1.95$. This value is very close to the boundary and a very large sample size would have to be observed before there will be substantial decrease in boundary truncation.

In terms of a “best estimator”, there does not appear to be a clear choice between $\hat{\alpha}_{1}$ and $\hat{\alpha}_{3}$. In Tables 1 and 2, the RMSE of $\hat{\alpha}_{3}$ is generally smaller than that of $\hat{\alpha}_{1}$, but the RMSE values are very close to one another. Figure 1 shows the RMSE of the three estimators for $\alpha\in\left(0.5,2\right)$ with $\beta=0$, $n+2m=300$ and $k=9$ based on $N=2000$ samples for each value of $\alpha$. A simple smoother was applied to the RMSE values to enhance readability of plot. Similar plots (not shown here) were produced for settings with $\beta=0.75$ and also $n+2m=600$; the same general trends were visible in these.

Figure 1 About Here

Inspection of Figure 1 shows that $\hat{\alpha}_{1}$ and $\hat{\alpha}_{3}$ perform better than $\hat{\alpha}_{2}$ over a large part of the parameter space. The estimator $\hat{\alpha}_{2}$ performs better in the approximate range $1.4\leq\alpha\leq 1.8$, but performs very poorly outside this range. As the true $\alpha$ becomes smaller and the underlying data distribution has heavier tails, $\hat{\alpha}_{1}$ has the best performance among the three estimators. When $\alpha>1.8$, $\hat{\alpha}_{3}$ performs better than the other two estimators. Generally, the RMSE of $\hat{\alpha}_{1}$ and $\hat{\alpha}_{3}$ are very similar with $\hat{\alpha}_{1}$ having smaller RMSE over a large range of $\alpha$.

5 Choice of $k$ and $(t_{1},\ldots,t_{k})$

The companion questions of how large to choose $k$ and how to choose the values $t_{1},\ldots,t_{k}$ have not yet been addressed. Intuition might suggest that one should choose $k$ as large as possible. However, the estimator can perform very poorly when $k$ is chosen too large. To illustrate, samples of size $n+2m=300$ were generated from a Cauchy distribution ($\alpha=1$) and estimators $\hat{\alpha}_{2}^{(k)}$ were calculated for each sample using equi-spaced $t_{j}=j/(k+1)$, $j=1,\ldots,k$ for $k=9,19$ and $29$. Based on $N=1000$ simulated samples, $RMSE(\hat{\alpha}_{1}^{(k=9)})=0.0919$, $RMSE(\hat{\alpha}_{2}^{(k=19)})=0.1253$ and $RMSE(\hat{\alpha}_{2}^{(k=29)})=0.1353$. The accuracy of the estimator decreases as $k$ increases. The AL method of Potgieter & Lombard (2012) is similar to a GLS method described by Hsieh (1995). A result from the latter shows $k=o(n^{-1/6})$ to be the optimal sample-size dependent rate for $k$. Practically, this means $k$ should not be chosen too large. Based on extensive simulation work, it is recommended that choices $k=9$ or $k=19$ be used for most applications. These values performed well across a wide range of sample sizes. If the underlying distribution has very heavy tails, say $\alpha<1$, the choice $k=3$ tends to be more robust to outliers.

Once a choice of $k$ has been made, it is still unclear what the best approach is to choosing the values $t_{1},\ldots,t_{k}$. This will be investigated in terms of the asymptotic distribution of the estimation of $\alpha$ based on a single random split. Potgieter & Lombard (2012) derive an expression for the asymptotic covariance matrix of $\left(mn/(m+n)\right)^{1/2}\left[\hat{\mu}-\mu,\hat{\sigma}-\sigma\right]^{\top}$ where $\hat{\mu}$ and $\hat{\sigma}$ denote the AL estimators based on two independent samples. Denote this covariance matrix $\sigma^{2}\Gamma$. The elements of the matrix $\Gamma$ are somewhat tedious expressions involving the points $t_{1},\ldots,t_{k}$, as well as the density and quantile functions $f(x)$ and $F^{-1}(t)$ of the underlying distribution. These expressions are omitted for brevity. For $\hat{\alpha}=\log 2/\log\hat{\sigma}$, a standard application of the delta method gives

[TABLE]

where $\Gamma_{22}$ indicates the element in the $2^{nd}$ row, $2^{nd}$ column of $\Gamma$.

The asymptotic variance in (10) does provide a view of the difficulty inherent in choosing “optimal” $t$-values. Noting that $\Gamma_{22}$ is an implicit function $\alpha$ and $\beta$ for the underlying stable distribution, one could choose the $t$-values such that $\alpha^{4}\Gamma_{22}$ is a minimized. Of course, in practice this is not possible as $\alpha$ and $\beta$ are unknown. Simulation studies not reported here were done to see how such optimal values would compare against the simple choice of equi-spaced values $t_{j}=j/(k+1)$, $j=1,\ldots,k$. While improvement in RMSE was observed in some instances, the simple choice of equi-spaced was usually very competitive with (and in a few instances outperfomed) the asymptotic optimal values. Therefore, the recommendation is made here to use equally spaced values. Further study is recommended to develop an adaptive approach for choosing the $t$-values. For example, after finding an initial estimate of $\alpha$, this initial estimate can be used to update the $t$-values and then re-estimate $\alpha$.

6 Comparison of estimators

In this section, the split–sample AL approach developed is compared with two other existing methods for estimating the stability index. The first of these, maximum likelihood (ML), is computationally expensive. However, it is included here as a performance benchmark as maximum likelihood estimators are asymptotically unbiased and minimum variance estimators. Next, McCulloch’s quantile estimator (MQE) is also considered. The MQE estimator is widely used in practice and, similar to the split–sample estimator, is based on sample order statistics.The split–sample estimator $\hat{\alpha}_{1}$ based on $B=250$ random partitions and using equi-spaced $t$-values, $t_{j}=j/(k+1)$ for $j=1,\ldots,k$ and choices $k=3,9$ and $19$ (hereafter referred to as SSEk). Samples were drawn from a $S(\alpha,\beta,1,0)$ distribution with $\alpha\in\{0.5,0.75,\ldots,1.75,2\}$ and $\beta\in\{0,0.75\}$ for sample sizes $n\in\{150,300,600\}$. The RMSE was estimated from $N=2000$ samples drawn for each parameter configuration. The results for the case $\alpha=1$ are reported below in Table 4. In this table, the boldface entries correspond to the estimates with the smallest and second smallest RMSE.

Table 4 About Here

Several interesting observations can be made upon inspection of Table 4. Consider first the symmetric case where $\beta=0$. As one would expect, the ML estimator has the smallest RMSE among all the estimators considered. When the sample size is small ($n+2m<=300$), the MQE has the second smallest RMSE. It should be noted, however, that it only performs marginally better than the split–sample estimator. Specifically, SSE3 is competitive at sample size $150$ as is SSE9 is at sample size $300$. For sample $600$, SSE9 outperforms MQE. A few additional simulations were done at even larger sample sizes and this trend was also observed there. In larger samples, SSE9 always has smaller RMSE than MQE. It should also be noted that no numerical values for RMSE were reported for estimator SSE${}_{1}9$ and sample size $150$, as convergence problems were frequently encountered when calculating the estimator. This is likely an artifact of using $t_{1}=0.05$ and $t_{19}=0.95$ (corresponding to the 5th and 95th sample percentiles) to do estimation in a small sample drawn from a heavy-tailed distribution. These sample percentiles have large variance and will often be values to what some might label extreme observations. In the asymmetric case ($\beta=0.75$), the estimator SSE9 always performs better than the MQE.

These simulation results are fairly representative of what was observed for other values of $\alpha$. Generally, SSE19 becomes the preferred estimator over SSE9 as either the value of $\alpha$ increases or as the sample size increases. Additionally, the relative efficiency of the methods was also estimated in the simulation study. Table 5 reports the ratios $\mathrm{RE}=\mathrm{RMSE}(\mathrm{SSE}_{9})/\mathrm{RMSE}(\mathrm{MQE})$ for a sample size of $600$. Here, a values of $\mathrm{RE}<1$ indicates superiority of SSE9 relative to MQE.

Table 5 About Here

In Table 5, it is evident that SSE9 generally performs much better than MQE. The estimated relative efficiency is often well below $0.7$. There are two notable exceptions to this general statement. In the symmetric case, when $\alpha$ ranges from $0.75$ to $1.25$, the MQE is very competitive, even outperforming SSE9 slightly at $\alpha=1.25$. In the asymmetric case, MQE performs better than SSE9 when $\alpha=0.5$, but nowhere else.

7 Recommendations

The split-sample approach shows promise as a method for estimating the stability index $\alpha$. In a side-by-side comparison with the McCulloch quantile estimator, the split-sample approach frequently outperforms the McCulloch estimator. The split-sample approach could conceivably be further improved by choosing the design points $t_{1},\ldots,t_{k}$ in some adaptive way as suggested at the end of Section 5. Additionally, the problem of estimating the standard error of the estimator has not been considered here. The bootstrap is one option, but does suffer from computational cost in that it becomes a nested problem involving a first-level bootstrap sampling procedure and a second-level data splitting procedure. Both of these questions are being considered by the authors in ongoing research.

As a final remark, the question of computational cost does arise when considering the practical implementation of the split-sample approach. Specifically, the process of permuting the data does become more time-consuming as the sample size gets large. However, the estimator proposed is highly parallelizable. In large sample situations, this will more than compensate for any computing time required to permute the data $B$ times.

Tables and Figures