Parameter estimation for stable distributions with application to commodity futures log-returns

Abstract

This paper explores the theory behind the rich and robust family of $\mathit{\alpha }$-stable distributions to estimate parameters from financial asset log-returns data. We discuss four-parameter estimation methods including the quantiles, logarithmic moments method, maximum likelihood (ML), and the empirical characteristics function (ECF) method. The contribution of the paper is two-fold: first, we discuss the above parametric approaches and investigate their performance through error analysis. Moreover, we argue that the ECF performs better than the ML over a wide range of shape parameter values, $\mathit{\alpha }$ including values closest to 0 and 2 and that the ECF has a better convergence rate than the ML. Secondly, we compare the t location-scale distribution to the general stable distribution and show that the former fails to capture skewness which might exist in the data. This is observed through applying the ECF to commodity futures log-returns data to obtain the skewness parameter.

Keywords:

• stable distribution
• parameter estimation
• density estimation

AMS Subject Classifications:

• 62G05
• 62G07
• 62G32

Public Interest Statement

This paper is entitled parameter estimation for stable distribution with applications to commodity future log-returns. The paper is useful to individuals interested in investing their wealth in financial markets. It provides essential information on how historical asset prices can inform future market movements via parameter estimation. This is crucial to portfolio managers, speculators, and hedgers. It’s imperative that the most accurate estimation method is established. Market data distribution deviates from the normal distribution, it exhibits skews, high or low peaks, and fat or skinny tails. The current paper is geared towards establishing the best estimation method among known methods in economic and financial analysis for skewed data.

1. Introduction

The motivation for this paper derives from the fact that parameter estimation from historical data is an important analysis to financial market participants. It provides useful information for portfolio managers, speculators, and hedgers. It is therefore, imperative that the most accurate estimation method is established. It is a known fact that in general, market data deviates from the Gaussian distribution, its distribution is either skewed, high or low peaked, and/or with fat or skinny tails. The current paper is geared towards establishing a better parameter estimation method among the commonly known ECF, ML, quantile, and logarithm moments methods used in economic and financial analysis for skewed data assumed to flow stable distributions.

The application of stable distributions in finance is traced way back in the late 50s when Mandelbrot (1959,1962,1963) developed a hypothesis that revolutionalized the way economists viewed and interpreted prices in speculative markets such as grains and securities markets. The hypothesis suggested that prices were not Gaussian as it had been previously believed by market participants based on Bachelier (1900). Mandelbrot’s hypothesis was therefore, an extension of the widely embraced breakthrough of Bachelier (1900).

In the following years Zolotarev (1964) developed integral representations of stable laws and the results have been used to develop parameter estimation techniques for the stable laws. Fama (1963) reviewed the validity of Mandelbrot’s hypothesis and came up with statistical tools suitable for dealing with speculative prices. Dumouchel (1971) employs this class of distributions in statistical inference for long-tailed data. Graphical representation of their densities and the estimation of their parameters via interpolation appear in Holt and Crow (1973) and in Koutrouvelis (1980) using regression. Parameter estimation methods based on quantile methods are presented in Fama and Roll (1971) for symmetric stable distributions but this approach faces a problem of discontinuity of the traditional location parameter in the asymmetrical cases when the exponent parameter passes unity. A remedy and generalization of the quantile approach is later introduced by McCulloch (1986).

A different parameter estimation technique based on fractional lower order moments (FLOM) appears in Ma and Nikias (1995) where the authors develop new methods for estimating parameters in impulsive signal environments. However, their methods only cover symmetric stable distributions. There was a need to extend the method to asymmetric systems. This came through by Kuruoğlu (2001) where a generalized FLOM method is introduced. Generally, FLOM methods pose a challenge of having to estimate the Sinc function and this in turn affects the accuracy of the results. As a consequence a better estimation approach referred to as logarithmic moments method (LM) is proposed by Kuruoğlu (2001) to avoid having to compute the Sinc.

The third estimation method utilizes the maximum likelihood (ML). It is known that the ML approach is widely favored in economic and financial applications due to its generality and asymptotic efficiency (see for instance, Yu, 2004). However, there are cases where the ML method can be unreliable especially when the likelihood function is not tractable, or its not bounded over the parameter space or does not have a closed form representation. For instance, in this current paper the densities considered do not have closed form expressions. However, since there is a one-to-one correspondence between the density function and its Fourier transform it could be worth exploiting the latter since it always exists and its bounded. This leads us to next estimation method.

The fourth estimation approach is the empirical characteristic function (ECF) method discussed in Yu (2004). Although the likelihood function can be unbounded, its Fourier transform is always bounded and, while the likelihood function might not be tractable or could not be of a closed form, the Fourier transform could have a closed form expression. The Fourier transform of the density function is the characteristic function (CF), hence the name empirical characteristic function (ECF) method. In this paper we aim to show that this approach performs better than all the previously mentioned methods. A useful software package that can be used to estimate stable distributions is provided in Nolan (1997). A more theoretical approach to statistical estimation of the parameters of stable laws is extensively discussed in Zolotarev (1980). Readers interested in how to simulate stable process can refer to two excellent literatures of Weron and Weron (1995) and Zolotarev (1986).

This paper explores the theory behind the rich and robust family of $\mathit{\alpha }$-stable distributions to estimate parameters from financial asset log-returns data. We discuss four-parameter estimation methods including the quantiles, logarithmic moments method, ML, and the empirical characteristics function (ECF) method. The contribution of the paper is two-fold: first, we discuss the above parametric approaches and investigate their performance through error analysis. Moreover, we argue that the ECF performs better than the ML over a wide range of shape parameter values, $\mathit{\alpha }$ including values closest to 0 and 2 and that the ECF has a better convergence rate than the ML. Secondly, we compare the t location-scale distribution to the general stable distribution and show that the former fails to capture skewness which might exist in the data. This is observed through applying the ECF to commodity futures log-returns data to obtain the stable parameters.

The rest of the paper is organized as follows: in Section 2 we define a stable process and its construction from independent and identically distributed random variables based on a generalized central limit theorem and discuss its characterization. In Section 3 we study the density and distribution properties of stable processes through their characteristic functions. Section 4 explains how the four-parameter estimation methods discussed in this paper work and provides an analysis on their accuracy. In Section 5 we study and analyze some commodity data and show that the data deviates from the normal distribution hypothesis. We use the ECF to obtain the four stable parameters from the data and in addition, fit the data to various distributions to determine the closest shape of the data which turns out to be the t location-scale distribution for all our data. This distribution is suited for data that is highly peaked and heavily tailed with outliers. However, we propose stable distribution fitting to check for any existing tails. Section 6 concludes.

2. Stable processes

Stable also known as alpha-stable (or equivalently $\mathit{\alpha }$-stable) processes belong to a general class of Lévy distributions. They are limiting distributions with a definitive exponent parameter $\mathit{\alpha }$ that determines their shape.

2.1. Definition and construction

Definition 2.1

Let $X_1,X_2,\dots ,X_n$ be independent and identically distributed random variables and suppose a random variable S defined by(1) $$\begin{array}{c}\hfill S\Rightarrow \frac{1}{a_n}\left(\sum _{i=1}^n{X}_i-b_n\right),\end{array}$$(1)

where “$\Rightarrow$” represents weak convergence in distribution, $a_n$ is a positive constant and $b_n$ is real. Then S is a stable process and the constants $a_n$ and $b_n$ need not to be finite.

Definition 2.1 allows modeling of a number of natural phenomenon beyond normality using stable distributions. The fact that $a_n$ and $b_n$ do not necessarily have to be finite provides the generalized central limit theorem.

Definition 2.2

(Generalized Central Limit Theorem Rachev (2003))    Suppose $X_1,X_2,\dots$ denotes a sequence of independent and identically distributed random variables and let sequences $a_n\in \mathbb{R}$ and $b_n\in {\mathbb{R}}^{+}$. Then we can define a sequence(2) $$\begin{array}{c}\hfill Z_n:=\frac{1}{b_n}\left(\sum _{i=1}^n{X}_i-a_n\right)\end{array}$$(2)

of sums $Z_n$ such that their distribution functions weakly converge to some limiting distribution:(3) $$\begin{array}{c}\hfill \mathtt{P}(Z_n<x)\Rightarrow H(x),n⟶\infty ,\end{array}$$(3)

where H(x) is some limiting distribution.

The traditional central limit theorem assumes finite mean $a:=\mathtt{E}[X_i]$ and finite variance ${\mathit{\sigma }}^2:=\mathtt{V}ar[X_i]$ and defines the sequence of sums(4) $$\begin{array}{c}\hfill Z_n:=\frac{1}{\mathit{\sigma }\sqrt{n}}\left(\sum _{i=1}^n{X}_i-na\right),\end{array}$$(4)

such that the distribution functions of $Z_n$ weakly converge to $h^{sG}(x)$:(5) $$\begin{array}{c}\hfill \mathtt{P}(x_1<Z_n<x_2)\Rightarrow \underset{x_1}{\overset{x_2}{\int }}{h}^{sG}(x)dx,n⟶\infty \end{array}$$(5)

where $h^{sG}(x)$ denotes the standard Gaussian distribution(6) $$\begin{array}{c}\hfill h^{sG}(x)=\frac{1}{\sqrt{2\mathit{\pi }}}exp(-x^2/2).\end{array}$$(6)

Suppose the independent and identically distributed random variables $X_i$ equal to a positive constant c almost surely and the sequences $a_n$ and $b_n$ in (2) are defined by $a_n=(n-1)c$ and $b_n=1$, then $Z_n$ is also equal to c for all $n>0$ almost surely. In this case the random variables $X_i$ are mutually independent and as a result, the limiting distribution for the sums $Z_n$ belong to the stable family of distributions by definition. This is one reason why they are regarded as stable.

2.2. Parametrization

Definition 2.3

A stable distribution is a four-parameter family denoted by $S(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\mu })$:

(1)	$\mathit{\alpha }\in (0,2]$ is the characteristic exponent responsible for the shape of the distribution.
(2)	$\mathit{\beta }\in [-1,1]$ is responsible for skewness of the distribution.
(3)	$\mathit{\nu }>0$ is the scale parameter (it narrows or extends the distribution around $\mathit{\mu }$).
(4)	$\mathit{\mu }\in \mathbb{R}$ is the location parameter (it shifts the distribution to the left or the right).

Suppose a random variable s follows a stable distribution $S(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\mu })$ then the random variable $z=(s-\mathit{\mu })/\mathit{\nu }$ has the same-shaped distribution as s but with the location parameter $\mathit{\mu }=0$ and the scale parameter $\mathit{\nu }=1$. This is another reason why they are referred to as stable, the shape is maintained after any rescaling.

Densities of $\mathit{\alpha }$-stable distributions do not have closed-form representations except for the case of a Gaussian ($\mathit{\alpha }=2$), Cauchy ($\mathit{\alpha }=1,\mathit{\beta }=0$) and Inverse Gaussian or Pearson ($\mathit{\alpha }=0.5,\mathit{\beta }=\pm 1$) distributions.

(1)	Gaussian distribution $N(\mathit{\mu },{\mathit{\sigma }}^2)$: $S\left(2,0,\frac{\mathit{\sigma }}{\sqrt{2}},\mathit{\mu }\right)$. $$\begin{array}{c}\hfill h^G(x)=\frac{1}{\mathit{\sigma }\sqrt{2\mathit{\pi }}}exp\left(-\frac{{(x-\mathit{\mu })}^2}{2{\mathit{\sigma }}^2}\right);-\infty <x<\infty .\end{array}$$
(2)	Cauchy distribution: $S(1,0,\mathit{\nu },\mathit{\mu })$. $$\begin{array}{c}\hfill h^C(x)=\frac{1}{\mathit{\pi }}\frac{\mathit{\nu }}{{\mathit{\nu }}^2+{(1-x)}^2};-\infty <x<\infty .\end{array}$$
(3)	Levy distribution (Inverse-Gaussian or Pearson): $S(1/2,1,\mathit{\nu },\mathit{\mu })$. $$\begin{array}{c}\hfill h^L(x)=\sqrt{\frac{\mathit{\nu }}{2\mathit{\pi }}}{(x-\mathit{\mu })}^{-3/2}exp\left(-\frac{\mathit{\nu }}{2(x-\mathit{\mu })}\right);\mathit{\mu }<x<\infty .\end{array}$$

The densities are generally computed using characteristic functions through transformations such as the Fourier¹. One can also refer to the work of Zolotarev (1964,1980,1986) for straight-forward and easy-to-compute integral representations of stable distribution and density functions. The distribution functions for the different $\mathit{\alpha }$ values have been tabulated in Dumouchel (1971), Fama and Roll (1968) and Holt and Crow (1973).

3. Density and distribution properties

3.1. Special case

Let $(X_t,t\ge 0)$ denote a Lévy process. The characterization of $X_t$ is deduced from the Lévy-Khintchine formula.

Definition 3.1

(Lévy-Khintchine & Applebaum, 2004)    Let $X={(X_t)}_{t\ge 0}$ be a Lévy process. There exist $b\in \mathbb{R}$, $\mathit{\sigma }\ge 0$ such that the characteristic function of X is given by(7) $$\begin{array}{c}\hfill \mathrm{\Phi }(t):=\mathtt{E}[e^{itX}]=exp\left(itb-\frac{1}{2}{\mathit{\sigma }}^2{t}^2+\underset{\mathbb{R}-\{0\}}{\int }(e^{itx}-1-itx{1}_{|x|<1})m(dx)\right),\end{array}$$(7)

where $1_{\{·\}}$ is an indicator function and m is a $\mathit{\sigma }$-finite measure satisfying the constraint(8) $$\begin{array}{c}\hfill \underset{\mathbb{R}-\{0\}}{\int }{min(1,|x|}^2)m(dx)<\infty ;\text{alternatively}\underset{\mathbb{R}-\{0\}}{\int }\frac{{|x|}^2}{1+{|x|}^2}m(dx)<\infty .\end{array}$$(8)

Definition 3.2

(The Lévy-Itô Decomposition Applebaum (2004))    If $X_t$ is a Lévy process, there exist $b\in \mathbb{R}$, a Brownian motion $B_{\mathit{\sigma }}(t)$ with variance $\mathit{\sigma }\in {\mathbb{R}}^{+}$ and an independent Poisson random measure N on ${\mathbb{R}}^{+}\times (\mathbb{R}-\{0\})$ such that, for each $t\ge 0$,(9) $$\begin{array}{c}\hfill X_t=bt+B_{\mathit{\sigma }}(t)+\underset{|x|<1}{\int }x\tilde{N}(t,dx)+\underset{|x|\ge 1}{\int }xN(t,dx),\end{array}$$(9)

where(10) $$\begin{array}{c}\hfill b=\mathtt{E}\left[X_1-\underset{|x|\ge 1}{\int }xN(1,dx)\right].\end{array}$$(10)

The compensated compound Poisson random measure is defined by $\tilde{N}=N-t\mathit{\lambda }$ to preserve the martingale property. The Lévy measure $\mathit{\lambda }$ satisfies (8).

A stable distribution can be constructed by setting $\mathit{\sigma }$ to zero in (7) or the second term on the right of (9) to zero and the Lévy measure in (8) to(11) $$\begin{array}{c}\hfill m(dx)=\frac{C}{{|x|}^{1+\mathit{\alpha }}}\text{d}x;C>0,\end{array}$$(11)

This gives a pure jump Lévy process which is a simple example of a stable family of distributions. We discuss a general case in the following.

3.2. General case

In the following, ${(S_t)}_{t\ge 0}$ will represent a stable process. Its characteristic function $\mathrm{\Phi }$ is obtained using the definition of domain of attraction of stable random variables and the Lévy-Khinchine representation formula in Definition 3.1 (see Applebaum, 2004):(12) $$\begin{array}{c}\hfill \mathrm{\Phi }(\mathit{\theta })=\left\{\begin{array}{c}exp\left(-{\mathit{\nu }}^{\mathit{\alpha }}{|\mathit{\theta }|}^{\mathit{\alpha }}\left(1-i\mathit{\beta }\text{sign}(\mathit{\theta })tan\left(\frac{\mathit{\pi }\mathit{\alpha }}{2}\right)\right)+i\mathit{\mu }\mathit{\theta }\right);\text{for}\mathit{\alpha }\ne 1.\hfill \\ exp\left(-\mathit{\nu }|\mathit{\theta }|\left(1+i\mathit{\beta }\text{sign}(\mathit{\theta })\frac{2}{\mathit{\pi }}log|\mathit{\theta }|\right)+i\mathit{\mu }\mathit{\theta }\right);\text{for}\mathit{\alpha }=1.\hfill \end{array}\right.\end{array}$$(12)

Alternative forms of parametrization are discussed in McCulloch (1986) for easier numerical analysis. More discussion on this to follow in Section 3.4.

The density of $S_t$ is computed from (12) using the Fourier transform:(13) $$\begin{array}{c}\hfill h_{S_t}(s)=\frac{1}{2\mathit{\pi }}\underset{-\infty }{\overset{\infty }{\int }}{e}^{-ist}\mathrm{\Phi }(t)dt.\end{array}$$(13)

Figure 1 shows density graphs for different exponent parameter values. The density is defined over the whole real line and for application purposes in finance log-returns data is usually used instead of raw asset prices to fit this family of distributions.

Figure 1. $\mathit{\alpha }$-stable densities for $\mathit{\alpha }\in (0,2]$.

The drawback in approximating (13) is that elementary techniques such as expressing the integral in terms of simple functions or using infinite polynomial expressions of the density function are not sufficient for meaningful numerical analysis. Some authors propose a standard parameterized integral expression of the density given by (see Ament & O’Neal, 2016)(14) $$\begin{array}{c}\hfill h_{S_t}(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\mu })=\frac{1}{\mathit{\sigma }\mathit{\pi }}\underset{0}{\overset{\infty }{\int }}{e}^{-t^{\mathit{\alpha }}}·cos\left(t·\left(\frac{s-\mathit{\mu }}{\mathit{\sigma }}\right)-\mathit{\beta }{t}^{\mathit{\alpha }}tan\left(\frac{\mathit{\pi }\mathit{\alpha }}{2}\right)\right)dt.\end{array}$$(14)

However, this representation consists of an oscillating integrand which in turn leads to another alternative approach presented in Zolotarev (1986) where the density of $S_t$ is given by(15) $$\begin{array}{cc}\hfill h_{S_t}(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\mu })& =\left\{\begin{array}{cc}\frac{\mathit{\alpha }{\left|\frac{s-\mathit{\mu }}{\mathit{\sigma }}\right|}^{\frac{1}{\mathit{\alpha }-1}}}{2\mathit{\sigma }|\mathit{\alpha }-1|}\underset{-\mathit{\theta }}{\overset{1}{\int }}{U}_{\mathit{\alpha }}(\mathit{\phi },\mathit{\theta })exp\left(-{\left|\frac{s-\mathit{\mu }}{\mathit{\sigma }}\right|}^{\frac{\mathit{\alpha }}{\mathit{\alpha }-1}}{U}_{\mathit{\alpha }}(\mathit{\phi },\mathit{\theta })\right)d\mathit{\phi };\hfill & \text{if}s\ne \mathit{\mu }\hfill \\ \frac{1}{\mathit{\pi }\mathit{\sigma }}·\mathrm{\Gamma }\left(1+\frac{1}{\mathit{\alpha }}\right)·cos\left(\frac{1}{\mathit{\alpha }}arctan\left(\mathit{\beta }·tan\left(\frac{\mathit{\pi }\mathit{\alpha }}{2}\right)\right)\right);\hfill & \text{if}s=\mathit{\mu }\hfill \end{array}\right.\hfill \end{array}$$(15) (16) $$\begin{array}{cc}\hfill U_{\mathit{\alpha }}(\mathit{\phi },\mathit{\vartheta })& ={\left(\frac{sin\left(\frac{\mathit{\pi }}{2}\mathit{\alpha }(\mathit{\phi }+\mathit{\vartheta })\right)}{cos\left(\frac{\mathit{\pi }\mathit{\phi }}{2}\right)}\right)}^{\frac{\mathit{\alpha }}{1-\mathit{\alpha }}}·\left(\frac{cos\left(\frac{\mathit{\pi }}{2}(\mathit{\alpha }-1)\mathit{\phi }+\mathit{\alpha }\mathit{\vartheta }\right)}{cos\left(\frac{\mathit{\pi }\mathit{\phi }}{2}\right)}\right),\hfill \end{array}$$(16)

where $\mathit{\theta }=arctan\left(\mathit{\beta }tan\frac{\mathit{\pi }\mathit{\alpha }}{2}\right)\frac{2}{\mathit{\alpha }\mathit{\pi }}\text{sign}(s-\mathit{\mu })$.

3.3. Some properties of stable distribution functions

Firstly, recall that for any two admissible sets of parameters of stable distributions we can find two unique numbers $a>0$ and b such that(17) $$\begin{array}{c}\hfill S(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\mu })\stackrel{d}{=}aS(\mathit{\alpha },\mathit{\beta },{\mathit{\nu }}^{\prime },{\mathit{\mu }}^{\prime }),\end{array}$$(17)

where(18) $$\begin{array}{c}\hfill a=\frac{\mathit{\nu }}{{\mathit{\nu }}^{\prime }},b=\left\{\begin{array}{cc}\mathit{\mu }-{\mathit{\mu }}^{\prime }\frac{\mathit{\nu }}{{\mathit{\nu }}^{\prime }},\hfill & \mathit{\alpha }\ne 1\hfill \\ \mathit{\mu }-{\mathit{\mu }}^{\prime }\frac{\mathit{\nu }}{{\mathit{\nu }}^{\prime }}+\mathit{\nu }\mathit{\beta }\frac{2}{\mathit{\pi }}log\frac{\mathit{\nu }}{{\mathit{\nu }}^{\prime }},\hfill & \mathit{\alpha }=1.\hfill \end{array}\right.\end{array}$$(18)

The intuition is that a general stable distribution can be expressed in terms of a standard stable distribution. That is, we can write $S(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\mu })\stackrel{d}{=}aS(\mathit{\alpha },\mathit{\beta },1,0)+b$ where(19) $$\begin{array}{c}\hfill a=\mathit{\nu },b=\left\{\begin{array}{cc}\mathit{\mu };\hfill & \mathit{\alpha }\ne 1\hfill \\ \mathit{\mu }+\mathit{\nu }\mathit{\beta }\frac{2}{\mathit{\pi }}log\mathit{\nu };\hfill & \mathit{\alpha }=1.\hfill \end{array}\right.\end{array}$$(19)

Secondly, suppose h, H and $\mathrm{\Phi }$ denote the respective probability, cumulative density and characteristic functions of a stable random variable, S, where$$\begin{array}{c}\hfill h(s,\mathit{\alpha },\mathit{\beta })=\frac{1}{2\mathit{\pi }}\underset{-\infty }{\overset{\infty }{\int }}(cosst-isinst)\mathrm{\Phi }(t,\mathit{\alpha },\mathit{\beta })dt,\end{array}$$

then it is readily seen that the following properties hold:

(1)	$h(-s,\mathit{\alpha },\mathit{\beta })=h(s,\mathit{\alpha },-\mathit{\beta })$.
(2)	$H(-s,\mathit{\alpha },\mathit{\beta })=1-H(s,\mathit{\alpha },-\mathit{\beta })$.
(3)	$\mathrm{\Phi }(-s,\mathit{\alpha },\mathit{\beta })=\mathrm{\Phi }(s,\mathit{\alpha },-\mathit{\beta })$.

The above three relations can be verified by trigonometric properties.

3.4. Simulating $\mathit{\alpha }$-stable random variables

The two excellent references for simulating stable processes are Zolotarev (1986) and Chambers, Mallows, and Stuck (1976).

Definition 3.3

Suppose $S_t$ is a stable process with parameters $(\mathit{\alpha },{\mathit{\beta }}_2,{\mathit{\nu }}_2,\mathit{\mu })$, the characteristic function is given by(20) $$\begin{array}{c}\hfill ln\mathrm{\Phi }(t)=\left\{\begin{array}{cc}i\mathit{\mu }t-{\mathit{\nu }}_2^{\mathit{\alpha }}{|t|}^{\mathit{\alpha }}exp(-i{\mathit{\beta }}_{2}sign(t))\frac{\mathit{\pi }}{2}K(\mathit{\alpha })),\hfill & \mathit{\alpha }\ne 1;\hfill \\ i\mathit{\mu }t-{\mathit{\nu }}_2|t|\left(\frac{\mathit{\pi }}{2}+i{\mathit{\beta }}_{2}sign(t))ln|t|\right);\hfill & \mathit{\alpha }=1;\hfill \end{array}\right.\end{array}$$(20)

where(21) $$\begin{array}{cc}\hfill K(\mathit{\alpha })& =\mathit{\alpha }-1+\text{sign}(1-\mathit{\alpha })=\left\{\begin{array}{cc}\mathit{\alpha };\hfill & \mathit{\alpha }\ne 1\hfill \\ \mathit{\alpha }-2;\hfill & \mathit{\alpha }=1\hfill \end{array}\right.\hfill \end{array}$$(21) (22) $$\begin{array}{cc}\hfill ({\mathit{\beta }}_2,{\mathit{\nu }}_2)& =\left\{\begin{array}{cc}\frac{2}{\mathit{\pi }K(\mathit{\alpha })}{tan}^{-1}\left(\mathit{\beta }tan\frac{\mathit{\pi }\mathit{\alpha }}{2}\right),\mathit{\nu }{\left(1+{\mathit{\beta }}^2{tan}^2\frac{\mathit{\pi }\mathit{\alpha }}{2}\right)}^{\frac{1}{2\mathit{\alpha }}};\hfill & \mathit{\alpha }\ne 1\hfill \\ \mathit{\beta },\frac{2}{\mathit{\pi }}\mathit{\nu };\hfill & \mathit{\alpha }=1\hfill \end{array}\right.\hfill \end{array}$$(22)

Lemma 3.4

Let $\mathit{\gamma }\in \left[-\frac{\mathit{\pi }}{2},\frac{\mathit{\pi }}{2}\right]$ be a uniformly distributed random variable and let W be an independent exponential random variable with mean 1. Then(23) $$\begin{array}{c}\hfill S=\left\{\begin{array}{cc}\frac{sin\mathit{\alpha }\left(\mathit{\gamma }+\frac{\mathit{\pi }}{2}{\mathit{\beta }}_2\frac{K(\mathit{\alpha })}{\mathit{\alpha }}\right)}{{(cos\mathit{\gamma })}^{\frac{1}{\mathit{\alpha }}}}{\left(\frac{cos\left(\mathit{\gamma }-\mathit{\alpha }\left(\mathit{\gamma }+\frac{\mathit{\pi }}{2}{\mathit{\beta }}_2\frac{K(\mathit{\alpha })}{\mathit{\alpha }}\right)\right)}{W}\right)}^{\frac{1-\mathit{\alpha }}{\mathit{\alpha }}};\hfill & \mathit{\alpha }\ne 1\hfill \\ \left(\frac{\mathit{\pi }}{2}+{\mathit{\beta }}_2\mathit{\gamma }\right)tan\mathit{\gamma }-{\mathit{\beta }}_{2}log\left(\frac{Wcos\mathit{\gamma }}{\frac{\mathit{\pi }}{2}+{\mathit{\beta }}_2\mathit{\gamma }}\right);\hfill & \mathit{\alpha }=1\hfill \end{array}\right.\end{array}$$(23)

is a standard $\mathit{\alpha }$-stable process with parameters $(\mathit{\alpha },{\mathit{\beta }}_2,1,0)$.

Proof

See Zolotarev (1986). $\square$

A stable random variable can be easily generated using Lemma 3.4. Programming languages such as R or MATLAB can be utilized to generate a uniformly distributed random variable U on the interval $\left(-\frac{\mathit{\pi }}{2},\frac{\mathit{\pi }}{2}\right)$ and an independent exponential random variable E with mean 1². Then the stable random variable would be generated by computing(24) $$\begin{array}{c}\hfill S=\left\{\begin{array}{cc}{A}_{\mathit{\alpha },\mathit{\beta }}\frac{sin(\mathit{\alpha }(U+B_{\mathit{\alpha },\mathit{\beta }}))}{{(cosU)}^{\frac{1}{\mathit{\alpha }}}}{\left(\frac{cos(U-\mathit{\alpha }(U+B_{\mathit{\alpha },\mathit{\beta }}))}{E}\right)}^{\frac{1-\mathit{\alpha }}{\mathit{\alpha }}};\hfill & \mathit{\alpha }\ne 1\hfill \\ \frac{2}{\mathit{\pi }}\left(\left(\frac{\mathit{\pi }}{2}+\mathit{\beta }U\right)tanU-\mathit{\beta }log\left(\frac{\frac{\mathit{\pi }}{2}EcosU}{\frac{\mathit{\pi }}{2}+\mathit{\beta }U}\right)\right);\hfill & \mathit{\alpha }=1\hfill \end{array}\right.\end{array}$$(24)

where $A_{\mathit{\alpha },\mathit{\beta }}={\left(1+{\mathit{\beta }}^2{tan}^2\frac{\mathit{\pi }\mathit{\alpha }}{2}\right)}^{\frac{1}{2\mathit{\alpha }}}$ and $B_{\mathit{\alpha },\mathit{\beta }}=\frac{{tan}^{-1}\left(\mathit{\beta }tan\frac{\mathit{\pi }\mathit{\alpha }}{2}\right)}{\mathit{\alpha }}$.

3.5. Moments of stable processes

Statistical moments $\mathtt{E}[|·{|}^k]$ of stable distributions are finite only when $k\le \mathit{\alpha }$. Moreover, for $\mathit{\alpha }<2$ the variance is infinite, for $\mathit{\alpha }\in (0,1]$ the mean does not exist and the mean is zero when $\mathit{\alpha }\in (1,2)$. This is not always the case for symmetric stable distributions where $\mathit{\beta }=0$.

3.5.1. Fractional lower order moments

The FLOM is an alternative for computing moments of $\mathit{\alpha }$-stable random variables especially in situations where the mean and/or variance are infinite. FLOM representation formulas are discussed in Ma and Nikias (1995) for symmetric stable random data and its generalization to asymmetric stable random data in Kuruoğlu (2001). In the latter, if $S_i\sim S(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\gamma })$ and $\mathit{\alpha }\ne 1$, then$$\begin{array}{cc}\hfill \mathtt{E}[S^{<p>}]& =\frac{\mathrm{\Gamma }\left(1-\frac{p}{\mathit{\alpha }}\right)}{\mathrm{\Gamma }(1-p)}{\left|\frac{\mathit{\gamma }}{cos\mathit{\theta }}\right|}^{\frac{p}{\mathit{\gamma }}}\frac{sin\left(\frac{p\mathit{\theta }}{\mathit{\alpha }}\right)}{sin\left(\frac{p\mathit{\pi }}{2}\right)},\text{for}p\in (-2,-1)\cup (-1,\mathit{\alpha }).\hfill \\ \hfill \mathtt{E}[|S{|}^p]& =\frac{\mathrm{\Gamma }\left(1-\frac{p}{\mathit{\alpha }}\right)}{\mathrm{\Gamma }(1-p)}{\left|\frac{\mathit{\gamma }}{cos\mathit{\theta }}\right|}^{\frac{p}{\mathit{\gamma }}}\frac{cos\left(\frac{p\mathit{\theta }}{\mathit{\alpha }}\right)}{cos\left(\frac{p\mathit{\pi }}{2}\right)},\text{for}p\in (-1,\mathit{\alpha }).\hfill \end{array}$$

where $\mathit{\theta }=arctan\left(\mathit{\beta }tan\frac{\mathit{\alpha }\mathit{\pi }}{2}\right)$ and $\mathrm{\Gamma }$ denotes the Gamma function. From the above representations, moments with negative values of p are attainable. This results into the logarithmic moments approach that provides an easier way of estimating stable distribtuion parameters compared to the FLOM.

3.5.2. Logarithmic moments

This approach is as a result of the challenges encountered when using the FLOM method which requires computing Gamma functions, the inversion of the sinc function and it only works for some p. The current method suggests computing derivatives with respect to the moment order p resulting in moments of the logarithms of the stable process. We illustrate in the following.

Lemma 3.5

Let S denote a symmetric stable random variable and let $p\in \mathbb{R}$. Then(25) $$\begin{array}{c}\hfill M_n:={\mathtt{E}[(log|S|)}^n]=\underset{p\to 0}{lim}\frac{{\text{d}}^n}{\text{d}{p}^n}{\mathtt{E}[|S|}^p],n=1,2,\dots .\end{array}$$(25)

The moments follow readily for $n=1,2,\dots$. i.e.(26) $$\begin{array}{cc}\hfill M_1& =\mathtt{E}[log|S|]={\mathit{\phi }}_0\left(1-\frac{1}{\mathit{\alpha }}\right)+\frac{1}{\mathit{\alpha }}log\left|\frac{\mathit{\nu }}{cos\mathit{\theta }}\right|.\hfill \end{array}$$(26) (27) $$\begin{array}{cc}\hfill M_2& =\mathtt{E}[(log|S|-{\mathtt{E}[log|S|])}^2]={\mathit{\phi }}_1\left(\frac{1}{2}+\frac{1}{{\mathit{\alpha }}^2}\right)-\frac{{\mathit{\theta }}^2}{{\mathit{\alpha }}^2}.\hfill \end{array}$$(27) (28) $$\begin{array}{cc}\hfill M_3& =\mathtt{E}[(log|S|-{\mathtt{E}[log|S|])}^3]={\mathit{\phi }}_3\left(1-\frac{1}{{\mathit{\alpha }}^3}\right).\hfill \end{array}$$(28)

where $\mathit{\theta }=arctan(\mathit{\beta }tan\mathit{\alpha }\mathit{\pi }/2)$ and terms ${\mathit{\phi }}_k$ are given by ${\mathit{\phi }}_0=-0.57721566$, ${\mathit{\phi }}_1={\mathit{\pi }}^2/6$, $\mathit{\phi }=1.2020569$ derived from the polygamma function(29) $$\begin{array}{c}\hfill {\mathit{\phi }}_{k-1}=\frac{{\text{d}}^k}{\text{d}{x}^k}{log\mathrm{\Gamma }(x)|}_{x=1}.\end{array}$$(29)

Proof

3.6 The proof is provided in Kuruoğlu (2001). $\square$

4. Parameter estimation of stable processes

The four common methods for estimating parameters of stable processes include: quantiles method (see Fama & Roll, 1971; McCulloch, 1986,1996 ), the logarithmic moments method (see Kuruoğlu, 2001), the empirical characteristics method (see Yang, 2012), and the ML method (see Nolan, 2001). We investigate their accuracy in the following.

4.1. The quantiles method

The quantile method was pioneered by Fama and Roll (1971) but was much more appreciated through McCulloch (1986) after its extension to include asymmetric distributions and for cases where $\mathit{\alpha }\in [0.6,2]$ unlike the former approach that restricts it to $\mathit{\alpha }\ge 1$.

Suppose $\widehat{s}$ is a given data sample then the estimates for $\mathit{\alpha }$ and $\mathit{\beta }$ are given by $\widehat{\mathit{\alpha }}={\mathrm{\Theta }}_1({\widehat{\mathit{\vartheta }}}_{\mathit{\alpha }},{\widehat{\mathit{\vartheta }}}_{\mathit{\beta }})$ and $\widehat{\mathit{\beta }}={\mathrm{\Theta }}_2({\widehat{\mathit{\vartheta }}}_{\mathit{\alpha }},{\widehat{\mathit{\vartheta }}}_{\mathit{\beta }})$ where(30) $$\begin{array}{c}\hfill {\widehat{\mathit{\vartheta }}}_{\mathit{\alpha }}=\frac{{\widehat{s}}_{0.95}-{\widehat{s}}_{0.05}}{{\widehat{s}}_{0.75}-{\widehat{s}}_{0.25}},{\widehat{\mathit{\vartheta }}}_{\mathit{\beta }}=\frac{{\widehat{s}}_{0.95}+{\widehat{s}}_{0.05}-2{\widehat{s}}_{0.05}}{{\widehat{s}}_{0.95}-{\widehat{s}}_{0.05}}.\end{array}$$(30)

The notation ${\widehat{s}}_q$ represents the qth quantile of sample $\widehat{s}$ and, $\widehat{\mathit{\alpha }}$ and $\widehat{\mathit{\beta }}$ are obtained by functions ${\mathrm{\Theta }}_1({\widehat{\mathit{\vartheta }}}_{\mathit{\alpha }},{\widehat{\mathit{\vartheta }}}_{\mathit{\beta }})$ and ${\mathrm{\Theta }}_2({\widehat{\mathit{\vartheta }}}_{\mathit{\alpha }},{\widehat{\mathit{\vartheta }}}_{\mathit{\beta }})$ given in Tables III and IV in McCulloch (1986) through linear interpolation. Consequently, the scale parameter is given by(31) $$\begin{array}{c}\hfill \widehat{\mathit{\nu }}=\frac{{\widehat{s}}_{0.75}-{\widehat{s}}_{0.25}}{{\mathrm{\Theta }}_3(\widehat{\mathit{\alpha }},\widehat{\mathit{\beta }})},\end{array}$$(31)

where ${\mathrm{\Theta }}_3(\widehat{\mathit{\alpha }},\widehat{\mathit{\beta }})$ is given by Table V in McCulloch (1986). The consistent estimator $\mathit{\nu }$ is then obtained through interpolation.

Finally the location parameter $\mathit{\mu }$ is estimated through a new parameter defined by(32) $$\begin{array}{c}\hfill \mathit{\zeta }=\left\{\begin{array}{cc}\mathit{\mu }+\mathit{\beta }\mathit{\gamma }tan\frac{\mathit{\pi }\mathit{\alpha }}{2};\hfill & \mathit{\alpha }\ne 1\hfill \\ \mathit{\mu };\hfill & \mathit{\alpha }=1.\hfill \end{array}\right.\end{array}$$(32)

Moreover, $\mathit{\zeta }$ is estimated by(33) $$\begin{array}{c}\hfill \widehat{\mathit{\zeta }}={\widehat{s}}_{0.5}+\widehat{\mathit{\nu }}{\mathrm{\Theta }}_5(\widehat{\mathit{\alpha }},\widehat{\mathit{\beta }}),\end{array}$$(33)

where ${\mathrm{\Theta }}_5(\widehat{\mathit{\alpha }},\widehat{\mathit{\beta }})$ is obtained from Table VII (McCulloch, 1986) by linear interpolation. The location parameter is estimated consistently by(34) $$\begin{array}{c}\hfill \widehat{\mathit{\mu }}=\widehat{\mathit{\zeta }}+\widehat{\mathit{\beta }}\widehat{\mathit{\nu }}tan\frac{\mathit{\pi }\widehat{\mathit{\alpha }}}{2}.\end{array}$$(34)

4.2. Empirical characteristic function method

Suppose a set of observable data $\{s_1,s_2,\dots ,s_N\}$ follows a stable distribution. Then we can approximate the characteristic function of this data by applying a basic Monte Carlo approach based on the law of large numbers i.e.(35) $$\begin{array}{c}\hfill \mathrm{\Phi }(u)=\mathtt{E}[e^{iu{s}_j}]\approx \widehat{\mathrm{\Phi }}(u)=\frac{1}{N}\sum _{j=1}^N{e}^{iu{s}_j}.\end{array}$$(35)

We can express the characteristic function (12) in terms of the cosine and sine function from basic trigonometric principles, i.e.(36) $$\begin{array}{c}\hfill \mathrm{\Phi }(u)=e^{-{|\mathit{\nu }u|}^{\mathit{\alpha }}}(cos\mathit{\eta }+isin\mathit{\eta }),\end{array}$$(36)

where$$\begin{array}{cc}\hfill \mathit{\eta }& =\mathit{\nu }u-{|\mathit{\nu }u|}^{\mathit{\alpha }}\mathit{\beta }\text{sign}(u)\mathit{\omega }(u,\mathit{\alpha })\hfill \\ \hfill \mathit{\omega }(u,\mathit{\alpha })& =\left\{\begin{array}{cc}tan\frac{\mathit{\pi }\mathit{\alpha }}{2},\hfill & \mathit{\alpha }\ne 1\hfill \\ \frac{2log|u|}{\mathit{\pi }},\hfill & \mathit{\alpha }=1\hfill \end{array}\right.\hfill \end{array}$$

As a result, we observe that(37) $$\begin{array}{c}\hfill |\mathrm{\Phi }(u)|=e^{-{|\mathit{\nu }u|}^{\mathit{\alpha }}}.\end{array}$$(37)

The estimated characteristic function relates to the model parameters by(38) $$\begin{array}{c}\hfill log|\widehat{\mathrm{\Phi }}(u_k)|={\mathit{\nu }}^{\mathit{\alpha }}{|u_k|}^{\mathit{\alpha }};\text{for}k=1,2,u_k>0,\mathit{\alpha }\ne 1.\end{array}$$(38)

Solving this system leads to the estimation representation formulas for the stability and variance parameters:$$\begin{array}{cc}\hfill \widehat{\mathit{\alpha }}& =\frac{log\frac{log|\widehat{\mathrm{\Phi }}(u_1)|}{log|\widehat{\mathrm{\Phi }}(u_2)|}}{log\left|\frac{u_1}{u_2}\right|}.\hfill \\ \hfill log\widehat{\mathit{\nu }}& =\frac{log|u_1|log(-log|\widehat{\mathrm{\Phi }}(u_2)|)-log|u_2|log(-log|\widehat{\mathrm{\Phi }}(u_1)|)}{log\left|\frac{u_1}{u_2}\right|}.\hfill \end{array}$$

The real and imaginary parts of the characteristic function (36) provide estimates for $\widehat{\mathit{\beta }}$ and $\widehat{\mathit{\mu }}$:(39) $$\begin{array}{c}\hfill arctan\frac{Im(\mathrm{\Phi }(u))}{Re(\mathrm{\Phi }(u))}=\mathit{\mu }u-{|\mathit{\nu }u|}^{\mathit{\alpha }}\mathit{\beta }\text{sign}(u)\mathit{\omega }(u,\mathit{\alpha }).\end{array}$$(39)

Suppose $\mathrm{{\rm Y}}(u):=arctan(Im\mathrm{\Phi }(u)/Re\mathrm{\Phi }(u))$ and choose another set of positive numbers $u_k,k=3,4$ together with $\widehat{\mathit{\alpha }}$ and $\widehat{\mathit{\nu }}$ then the estimates of the location and skewness parameters are given respectively by(40) $$\begin{array}{cc}\hfill \widehat{\mathit{\mu }}& =\frac{u_4^{\widehat{\mathit{\alpha }}}\mathrm{{\rm Y}}(u_3)-u_3^{\widehat{\mathit{\alpha }}}\mathrm{{\rm Y}}(u_4)}{u_3{u}_4^{\widehat{\mathit{\alpha }}}-u_3{u}_4^{\widehat{\mathit{\alpha }}}}.\hfill \end{array}$$(40) (41) $$\begin{array}{cc}\hfill \widehat{\mathit{\beta }}& =\frac{u_4\mathrm{{\rm Y}}(u_3)-u_3\mathrm{{\rm Y}}(u_4)}{{\widehat{\mathit{\nu }}}^{\widehat{\mathit{\alpha }}}tan\frac{\mathit{\pi }\widehat{\mathit{\alpha }}}{2}(u_4{u}_3^{\widehat{\mathit{\alpha }}}-u_3{u}_4^{\widehat{\mathit{\alpha }}})}.\hfill \end{array}$$(41)

Notice, it can be deduced from Equation (36) that$$\begin{array}{c}\hfill {log(-log(|\mathrm{\Phi }(u)|}^2))=log(2{\mathit{\nu }}^{\mathit{\alpha }})+\mathit{\alpha }log(u).\end{array}$$

This provides an alternative way to envision the regression estimation method:$$\begin{array}{c}\hfill y_k=m+\mathit{\alpha }{x}_k+{\mathit{\epsilon }}_k;k=1,2,\dots ,M;\end{array}$$

where $y_k=log(-log|\widehat{\mathrm{\Phi }}(u_k){|}^2),m=log(2{\mathit{\nu }}^{\mathit{\alpha }}),x_k=log(u_k)$ and ${\mathit{\epsilon }}_k$ is an error term. The stability parameter $\mathit{\alpha }$ and the scale parameter $\mathit{\nu }$ can be estimated by selecting $u_k=\frac{\mathit{\pi }k}{25},k=1,2,\dots ,M$; of real data (see Koutrouvelis, 1980, Table I). The estimates $\widehat{\mathit{\alpha }}$ and $\widehat{\mathit{\nu }}$ are then used to estimate $\mathit{\beta }$ and $\mathit{\mu }$ using the following relation$$\begin{array}{c}\hfill z_l={\mathit{\eta }}_l+{\mathit{\varsigma }}_l,l=1,2,\dots ,Q.\end{array}$$

where $z_l={\mathrm{{\rm Y}}}_n(u_l)+\mathit{\pi }{k}_n(u_l)$, ${\mathit{\eta }}_l={\widehat{\mathit{\nu }}}_{l}u-{|{\widehat{\mathit{\nu }}}_{l}u|}^{\widehat{\mathit{\alpha }}}\mathit{\beta }\text{sign}(u)\mathit{\omega }(u,\widehat{\mathit{\alpha }})$ and ${\mathit{\varsigma }}_l$ is some random error. The proposed real data set for Q (see Koutrouvelis, 1980, Table II) is $u_l=\frac{\mathit{\pi }l}{50},l=1,2,\dots ,Q$.

4.3. Logarithmic moments method

This approach follows the theory discussed in Section 3.5.2. The key innovation with this method is that there is no need of computing Gamma functions and the sinc function as in the FLOM. Secondly, techniques of parameter estimation for symmetric stable random variables (i.e. $\mathit{\beta }=0$) can be applied to skewed stable random variables (i.e. $\mathit{\beta }\ne 0$) and, techniques of parameter estimation for centered stable random variables (i.e. $\mathit{\mu }=0$) to non-centered ones (i.e. $\mathit{\mu }\ne 0$) through centro-symmetrization. However, this comes at a cost of losing almost half of the sample data. Therefore to obtain better estimates one has to use large sample data sets.

4.3.1. Centro-symmetrization of stable random data sets

Let $S_k$ be a sequence of n independent stable random variables distributed according to$$\begin{array}{c}\hfill S_k\sim S(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\mu }).\end{array}$$

Then the distribution of a weighted sum of the above sequence with weights $a_k$ can be estimated using their characteristic function:(42) $$\begin{array}{c}\hfill Z=\sum _{k=1}^n{a}_k{S}_k\sim S\left(\mathit{\alpha },\frac{\sum _{k=1}^n{a}_k^{<\mathit{\alpha }>}}{\sum _{k=1}^n{|a_k|}^{\mathit{\alpha }}}\mathit{\beta },\sum _{k=1}^n{|a_k|}^{\mathit{\alpha }}\mathit{\nu },\sum _{k=1}^n{a}_k\mathit{\mu }\right),\end{array}$$(42)

where the pth power of a number x is defined by$$\begin{array}{c}\hfill x^{<p>}=\text{sign}(x){|x|}^p.\end{array}$$

As a result, it is easy to obtain sequences of independent stable random variables with zero $\mathit{\mu }$, zero $\mathit{\beta }$ as well as both zero $\mathit{\mu }$ and zero $\mathit{\beta }$ for $\mathit{\alpha }\ne 1$. This yields the centred, deskewed, and symmetrized sequences:(43) $$\begin{array}{cc}\hfill S_k^C& =S_{3k}+S_{3k-1}-2S_{3k-2}\sim S\left(\mathit{\alpha },\left[\frac{2-2^{\mathit{\alpha }}}{2+2^{\mathit{\alpha }}}\right]\mathit{\beta },[2+2^{\mathit{\alpha }}]\mathit{\nu },0\right),\hfill \end{array}$$(43) (44) $$\begin{array}{cc}\hfill S_k^D& =S_{3k}+S_{3k-1}-2^{1/\mathit{\alpha }}{S}_{3k-2}\sim S(\mathit{\alpha },0,4\mathit{\nu },[2-2^{1/\mathit{\alpha }}]\mathit{\mu }),\hfill \end{array}$$(44) (45) $$\begin{array}{cc}\hfill S_k^S& =S_{2k}-S_{2k-1}\sim S(\mathit{\alpha },0,2\mathit{\nu },0).\hfill \end{array}$$(45)

4.3.2. Parameter estimation

Suppose $S_k$ is a data set assumed to be drawned from $S(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\mu })$. Then the exponent parameter $\mathit{\alpha }$ is estimated by setting $\mathit{\theta }=0$ in (27), and the log moment $M_2$ is estimated from the obverted data (45). That is,(46) $$\begin{array}{c}\hfill \widehat{\mathit{\alpha }}={\left(\frac{M_2}{{\mathit{\phi }}_1}-\frac{1}{2}\right)}^{-1/2}.\end{array}$$(46)

The estimated $\widehat{\mathit{\alpha }}$ is used to estimate $\mathit{\theta }$ using (26) where $M_1$ is estimated from the obverted data (44). That is,(47) $$\begin{array}{c}\hfill |\widehat{\mathit{\theta }}|={\left(\left(\frac{{\mathit{\phi }}_1}{2}-M_2\right){\widehat{\mathit{\alpha }}}^2+{\mathit{\phi }}_1\right)}^{1/2}.\end{array}$$(47)

From the definition of $\mathit{\theta }$, $|{\mathit{\beta }}_0|$ can be estimated by(48) $$\begin{array}{c}\hfill {\widehat{\mathit{\beta }}}_0=\frac{tan\widehat{\mathit{\theta }}}{tan\frac{\widehat{\mathit{\alpha }}\mathit{\pi }}{2}}.\end{array}$$(48)

Centering (see (43)) requires $|{\widehat{\mathit{\beta }}}_0|$ to be multiplied by $(2+2^{\mathit{\alpha }})/(2-2^{\mathit{\alpha }})$ to obtain $|\widehat{\mathit{\beta }}|$ of the original data where the sign of $\mathit{\beta }$ is determined by$$\begin{array}{c}\hfill K=\text{sign}(|S_{max}-S_{\text{md}}|-|S_{min}-S_{\text{md}}|),\text{such that}\widehat{\mathit{\beta }}=K|\widehat{\mathit{\beta }}|.\end{array}$$

where $S_{max}$, $S_{\text{md}}$, $S_{min}$ is the maximum, median and minimum of the original data.

Next we estimate the scale parameter ${\widehat{\mathit{\nu }}}_0$ using (26) where $M_1$ is estimated from the obverted data (43). That is(49) $$\begin{array}{c}\hfill {\widehat{\mathit{\nu }}}_0=|cos\widehat{\mathit{\theta }}|exp((M_1-{\mathit{\phi }}_0)\widehat{\mathit{\alpha }}+{\mathit{\phi }}_0).\end{array}$$(49)

Again centering (see (43)) gives the parameter estimate $\widehat{\mathit{\nu }}$ of the original data by $\widehat{\mathit{\nu }}={\widehat{\mathit{\nu }}}_0{(2-2^{1/\mathit{\alpha }})}^{-1}$.

Finally, the location parameter $\mathit{\mu }$ is estimated by(50) $$\begin{array}{c}\hfill \widehat{\mathit{\mu }}={\widehat{\mathit{\mu }}}_0{(2-2^{1/\mathit{\alpha }})}^{-1}.\end{array}$$(50)

where ${\mathit{\mu }}_0$ is the median or mean of the obverted data ().

4.4. Maximum likelihood method

The ML method is the most favored parameter estimation method in economic and financial applications. The method relies on the density function which in the case of stable distributions poses a closed form representation problem. In this case we propose a numerical estimation of the density function. For a vector $s=(s_1,s_2,\dots ,s_n)$ of independent identically distributed random variables assumed to follow a stable distribtion, the ML estimate of the parameter vector $\mathrm{\Theta }=(\mathit{\alpha },\mathit{\beta },\mathit{\nu },\mathit{\mu })$ is obtained by maximizing the log-likelihood function given by(51) $$\begin{array}{c}\hfill L_{\mathrm{\Theta }}(s)=\sum _{i=1}^{n}ln\tilde{h}(s_i;\mathrm{\Theta }),\end{array}$$(51)

where $\tilde{h}(s;\mathrm{\Theta })$ denotes a numerically estimated stable probability density function. It is shown for instance in Mittnik, Rachev, Doganoglu, and Chenyao (1999) that the best algorithms to compute the ML is by using Fast Fourier Transforms (FFT) or by direct integration method as in Nolan (2001). The ML algorithms require carefully chosen initial input parameters which in our case can be obtained for example, through the quantiles method described above. The FFT is faster for large data sets and the direct integral approach is suitable for smaller data sets since it can be evaluated at any arbitrary point.

In the following section, we analyze commodities and apply the empirical characteristic functions method to estimate the stable distribution parameters.

It is important to mention the restrictions on the parameters under which the different estimation methods operate.

4.5. Error analysis

In this section we simulate datasets from the stable family of distributions based on the theory in Chambers et al. (1976) and Weron and Weron (1995). Then use the above four methods to retrieve the stable parameters from the simulated data. Our focus is on the $\mathit{\alpha }$ and $\mathit{\beta }$ but the arguments extend to the other two parameters.

First, it is important to mention that all the four methods perform poorly close to the boundaries i.e. $\mathit{\alpha }\to 0$, $\mathit{\alpha }\to 2$ and $\mathit{\beta }\to \pm 1$. Moreover, literature shows that the methods operate efficiently under the parameter restrictions in Table 1.

Table 1. Estimation methods and their parameter restrictions

Estimation method	Parameter restrictions
Quantile	$\mathit{\alpha }\ge 0.1$
Logarithm moments	$\mathit{\beta }=\mathit{\mu }=0$
Maximum likelihood	$\mathit{\alpha }\ge 0.4$
Empirical characteristic function	$\mathit{\alpha }\ge 0.1$

In addition, the MLE seems the most preferred and used estimation method. However, we observe in our analysis that this method fails for particular parameter ranges and it is not robust. For instance in estimating $0.1<\mathit{\alpha }<1.0$ with respect to $\mathit{\beta }$, the MLE fails to converge and returns huge unrealistic errors. This is why we do not include it in Figure 2(a). Similarly, for $\mathit{\beta }=0.4$ estimation with respect to $\mathit{\alpha }$, the logarithm moments method returns either negative or very large $\mathit{\beta }$ values which is expected according to the constraints in Table 1. We omit its graph in Figure 2(b). Meanwhile, we notice that in both cases, the quantile and ECF methods work well with the latter providing relatively the best estimates.

Figure 2. Method comparison for $\mathit{\alpha }=0.4$ and $\mathit{\beta }=0.4$ estimation.

The graphs in Figure 3 show the error associated with estimating $1.0<\mathit{\alpha }<2.0$ for different $\mathit{\beta }$ values. Note that all the four methods work well and we still notice the ECF being relatively the most accurate and robust method. Recall that for $\mathit{\alpha }\to 1$ and $\mathit{\alpha }\to 2$ the estimation methods perform poorly. An example is Figure 3(a) (for $\mathit{\alpha }=1.4$) which was the closest for which the ML would converge but for higher $\mathit{\alpha }>1.4$ values but far less than 2.0 (see for instance, Figure 3(a) for $\mathit{\alpha }=1.7$) the methods performed relatively better except for the logarithm moments methods.

Figure 3. Method comparison for $\mathit{\alpha }=1.4$ and $\mathit{\alpha }=1.7$ estimation.

Figure 4. $\mathit{\beta }$ estimation for differing data set sizes and $\mathit{\alpha }$ values.

The graphs in Figure 4 illustrate convergence of the quantile, ECF and the MLE in estimating $\mathit{\alpha }=1.4$ and $\mathit{\alpha }=1.7$. We simulated 50,000 points and divided it into 100 sets starting with a 500-sized set and increasing it by 500 to 50,000. The logarithm moments method performed extremely poorly and incomparable to the above three methods. It is not included in Figure 4(a) and (b). The ECF is seen to perform better than the quantile and ML methods with a relatively better convergence rate. Similary Figure 5 shows the convergence rates for the quantile, ECF and ML estimation methods. The ECF still provides a better precision in both cases i.e. Figure 5(a) and (b).

Figure 5. $\mathit{\alpha }$ estimation for differing data set sizes for $\mathit{\beta }=0.4$ values.

In summary the empirical characteristic function method outperforms all the three other methods discussed in this paper in the following way:

(1)	It is robust and can consistently estimate a wide range of $\mathit{\alpha }$ and $\mathit{\beta }$ parameters.
(2)	It provides a better precision compared to the quantile, logarithm moments and MLE methods for a wide range of $\mathit{\alpha }$ and $\mathit{\beta }$ parameters.
(3)	It has a better convergence rate.

Therefore the quantile, logarithm moments or the ML methods can be used to provide initial parameters for the ECF method. Similarly, the latter can be used to provide initial parameters for better future estimators.

The following section is devoted to extracting stable parameters from log-returns commodity futures data using the ECF method.

5. Commodity data

The data sets used here are obtained from Quandl Financial and Economic Data website. The sets differ in sizes and include settled prices of Corn Futures Continuous Contract C#1 from 1959-07-01 to 2017-02-10; Crude Oil Futures Continuous Contract C#1 from 1983-03-30 to 2017-02-10; Gasoline Futures Continuous Contract C#1 from 2005-10-03 to 2017-02-10; Gold Futures Continuous Contract C#1 from 1974-12-31 to 2017-02-10; Natural Gas Futures Continuous Contract C#1 from 1990-04-03 to 2017-02-10; Platinum Futures Continuous Contract C#1 from 1969-01-02 to 2017-02-10; Silver Futures Continuous Contract C#1 from 1963-06-13 to 2017-02-10; Soybeans Futures Continuous Contract C#1 from 1959-07-01 to 2017-02-10; Wheat Futures Continuous Contract C#1 from 1959-07-01 to 2017-02-10. To avoid multi-distributional effects, we work with log-returns of the data sets.

5.1. The t-location-scale distribution

The t-location-scale distribution is most suited for modeling data distributions with heavier tails, more prone to outliers than the Gaussian distribution. The distribution uses the following parameters

Table

Parameter	Description	Support
$\mathit{\mu }$	Location parameter	$-\infty <\mathit{\mu }<\infty$
$\mathit{\nu }$	Scale parameter	$\mathit{\nu }>0$
${\mathit{\alpha }}^{\ast }$	Shape parameter	${\mathit{\alpha }}^{\ast }>0$

The probability density function (pdf) of the t-location-scale distribution is given by$$\begin{array}{c}\hfill h(x)=\frac{\mathrm{\Gamma }\left(\frac{{\mathit{\alpha }}^{\ast }+1}{2}\right)}{\mathit{\nu }\sqrt{{\mathit{\alpha }}^{\ast }\mathit{\pi }}\mathrm{\Gamma }\left(\frac{{\mathit{\alpha }}^{\ast }}{2}\right)}{\left(\frac{{\mathit{\alpha }}^{\ast }+{\left(\frac{x-\mathit{\mu }}{\mathit{\nu }}\right)}^2}{{\mathit{\alpha }}^{\ast }}\right)}^{-\left(\frac{{\mathit{\alpha }}^{\ast }+1}{2}\right)},\end{array}$$

where $\mathrm{\Gamma }(·)$ denotes the gamma function. The mean of the t-location-scale distribution is given by $\mathit{\mu }$ and it is defined for ${\mathit{\alpha }}^{\ast }>1$ and undefined otherwise. The variance is given by$$\begin{array}{c}\hfill \mathbb{V}ar={\mathit{\nu }}^2\frac{{\mathit{\alpha }}^{\ast }}{{\mathit{\alpha }}^{\ast }-2}.\end{array}$$

The t-location-scale distribution approaches the Gaussian distribution as ${\mathit{\alpha }}^{\ast }$ approaches infinity and smaller values of ${\mathit{\alpha }}^{\ast }$ yield heavier tails. This distribution does not take skewness into consideration and its three parameters are usually estimated using the ML estimation method.

Using algorithms by Sheppard (2012) on our log-returns commodity futures data we obtained fittings in Figures 6–8.

Table 2. t- location-scale distribution parameters extracted from the log-returns data

		t-location scale parameters
		$\mathit{\mu }$	$\mathit{\nu }$	${\mathit{\alpha }}^{\ast }$
	Gasoline	0.000905848	0.0193048	4.45922
3*Energy	Natural gas	0.000117668	0.0181944	2.50848
	Crude oil	0.000313708	0.00912801	1.75246
	Corn	$5.52294e-05$	0.00439924	3.03782
3*Grains	Soy beans	0.000400308	0.00535733	2.43218
	Wheat	$1.63112e-05$	0.0113397	3.41863
	Gold	0.000905944	0.0108491	2.70553
3*Metals	Platinum	0.00044219	0.0110761	3.13631
	Silver	$-1.72459e-05$	0.000682631	0.512196

According to the ${\mathit{\alpha }}^{\ast }$ values, the log-returns data exhibit some tails. To determine the nature of the details one would require to run some QQ plots but this can also be observed directly from the Figures 6–8.

Figure 6. Energy: The data exhibits high peaks and skinny tails.

Figure 7. Grains: The data exhibits high peaks and skinny tails.

Figure 8. Metals: The data exhibits high peaks and skinny tails.

It is important to mention that QQ plots do not straight away provide conclusive evidence about the nature of the tails. More tests would still need to be made. For instance under the t-location scale it is not obvious to observe any skewness in the data. We however, view this effect when we fit the data to stable distribution (see Table 2) as discussed in the following section.

5.2. Stable distribution fitting

On the other hand, by assuming stable distribution for our log-returns commodity futures data, we employed the ECF method and obtained the stable parameters in Table 3.

Table 3. Stable distribution parameters extracted from the log-returns data

		Stable distribution parameters
		$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\nu }$	$\mathit{\mu }$
	Gasoline	1.7504	$-0.3806$	0.0152	$-0.0005$
3*Energy	Natural gas	1.5329	0.0371	0.015	0.0005
	Crude oil	1.2322	$-0.1526$	0.0075	$-0.0022$
	Corn	1.651	0.2117	0.0036	0.0004
3*Grains	Soy beans	1.4665	$-0.0968$	0.0043	0.0001
	Wheat	1.638	0.0929	0.0091	0.0003
	Gold	1.5007	$-0.1324$	0.0088	0.0001
3*Metals	Platinum	1.5943	$-0.1339$	0.0089	$-0.0001$
	Silver	0.4461	0.0176	0.0011	$-0.0001$

Log-returns of commodity futures are not only high peaked but they also have left and right skinny tails with extreme outliers as observed from the QQ-plots for energy commodities (i.e. Crude oil, Natural gas and Gasoline) in Figure 9, the grains commodities in Figure 10 and the precious metals in Figure 11.

Figure 9. Energy: In all, the left and right tails are skinny.

Figure 10. Grains: In all, the left and right tails are skinny.

Figure 11. Metals: In all, the left and right tails are skinny.

Table 3 shows stable distribution parameters extracted from the log-returns data using the empirical characteristic function parameter estimation method. We notice that the data exhibit a bit of skewness which is not reflected in the t-location-scale distribution fitting.

6. Conclusion

First we showed that the ECF provides the best precision in estimating a wide range of $\mathit{\alpha }$ and $\mathit{\beta }$ parameters, it is robust and provides better convergence compared to the quantile, ML, and the logarithm moments. Secondly, we have illustrated that in general, the distribution of the commodity futures log-returns data is closest to a t-location-scale distribution due to its high peaks, skinny tails and extreme outliers. Moreover, by using the ECF estimation method we realize some minor skewness effects not captured in the t-location-scale fitting. We recommend the ECF as a suitable approach for estimating parameters of any skewed financial market data and could be used to obtain initial input parameters for future and better estimation techniques.

Cover image

Source: Original image from Parameter estimation for stable distributions with application to commodity futures log-returns by M. Kateregga, S. Mataramvura and D. Taylor.

Additional information

Funding

This work was supported by funds from the National Research Foundation of South Africa (NRF), the African Institute for Mathematical Sciences (AIMS) and the African Collaboration for Quantitative Finance and Risk Research (ACQuFRR) which is the research section of the African Institute of Financial Markets and Risk Management (AIFMRM), which delivers postgraduate education and training in financial markets, risk management and quantitative finance at the University of Cape Town in South Africa.

Notes on contributors

M. Kateregga

Mr Michael Kateregga is a finishing PhD student at the University of Cape Town in South Africa. His research is in the field of mathematical finance and his PhD thesis is entitled Stable Distributions with Applications in Finance. The current paper is a chapter in his thesis which is due for submission in August, 2017. Mr Kateregga is also a researcher at the African Collaboration for Quantitative Finance and Risk Research (ACQuFRR) which is the research section of the African Institute of Financial Markets and Risk Management (AIFMRM), which delivers postgraduate education and training in financial markets, risk management and quantitative finance. Mr Kateregga also works with the African Institute for Mathematical Sciences (AIMS) in South Africa as a Research Assistant.

Notes

1 Note that characteristic functions always exist.

2 These are easily obtained from in-built functions in MATLAB