Fitting the Seven-Parameter Generalized Tempered Stable Distribution to Financial Data

Abstract

This paper proposes and implements a methodology to fit a seven-parameter Generalized Tempered Stable (GTS) distribution to financial data. The nonexistence of the mathematical expression of the GTS probability density function makes maximum-likelihood estimation (MLE) inadequate for providing parameter estimations. Based on the function characteristic and the fractional Fourier transform (FRFT), we provide a comprehensive approach to circumvent the problem and yield a good parameter estimation of the GTS probability. The methodology was applied to fit two heavy-tailed data (Bitcoin and Ethereum returns) and two peaked data (S&P 500 and SPY ETF returns). For each historical data, the estimation results show that six-parameter estimations are statistically significant except for the local parameter, $\mu$. The goodness of fit was assessed through Kolmogorov–Smirnov, Anderson–Darling, and Pearson’s chi-squared statistics. While the two-parameter geometric Brownian motion (GBM) hypothesis is always rejected, the GTS distribution fits significantly with a very high p-value and outperforms the Kobol, Carr–Geman–Madan–Yor, and bilateral Gamma distributions.

1. Introduction

Modeling high-frequency asset return with the normal distribution is the underlying assumption in many financial tools, such as the Black–Scholes–Merton option pricing model and the risk metric variance–covariance technique for value at risk (VAR). However, substantial empirical evidence rejects the normal distribution for various asset classes and financial markets. The symmetric and rapidly decreasing tail properties of the normal distribution cannot describe the skewed and fat-tailed properties of the asset return distribution.

The $\alpha$-stable distribution has been proposed (Nolan 2020; Sato 1999) as an alternative to the normal distribution for modeling asset return and many types of physical and economic systems. The theoretical and empirical argument is that the stable distribution generalizes the Central Limit Theorem regardless of the variance nature (finite or infinite) (Nzokem 2024; Rachev et al. 2011). There are two major drawbacks (Borak et al. 2005; Nolan 2020): firstly, the lack of closed formulas for densities and distribution functions, except for the normal distribution ($\alpha =2$), Cauchy distribution ($\alpha =1$), and Lévy distribution ($\alpha ={\displaystyle \frac{1}{2}}$) (Tsallis 1997); secondly, most of the moments of the stable distribution are infinite. An infinite variance of asset return leads to an infinite price for derivative instruments such as options.

The Generalized Tempered Stable (GTS) distribution was developed to overcome the shortcomings of the two distributions, and the tails of the GTS distribution are heavier than the normal distribution but thinner than the stable distribution (Grabchak and Samorodnitsky 2010; Kim et al. 2009). The general form of the GTS distribution can be defined by the following Lévy measure ($V(dx)$) (1):

$$\begin{array}{c}\hfill V(dx)=\left({\displaystyle \frac{{\alpha }_{+}{e}^{-{\lambda }_{+}x}}{x^{1+{\beta }_{+}}}}{\mathbf{1}}_{x>0}+{\displaystyle \frac{{\alpha }_{-}{e}^{-{\lambda }_{-}\left|x\right|}}{{\left|x\right|}^{1+{\beta }_{-}}}}{\mathbf{1}}_{x<0}\right)dx\end{array}$$

(1)

where $0\le {\beta }_{+}\le 1$, $0\le {\beta }_{-}\le 1$, ${\alpha }_{+}\ge 0$, ${\alpha }_{-}\ge 0$, ${\lambda }_{+}\ge 0$, and ${\lambda }_{-}\ge 0$. More details on the Tempered Stable distribution are provided in (Küchler and Tappe 2013; Rachev et al. 2011).

The rich class of the GTS distribution (1) has a myriad of applications ranging from financial to mathematical physics and economic systems. However, few studies (Fallahgoul and Loeper 2021; Massing 2024; Nzokem and Montshiwa 2022) have covered the methods and techniques to estimate the parameters of the GTS distribution. This study aims to contribute to the literature by providing a methodology for fitting the seven-parameter GTS distribution. As illustrations, the study used four historical prices: two heavy-tailed data (Bitcoin and Ethereum returns) and two peaked data (S&P 500 and SPY ETF returns). The GTS distribution is fitted to the underlying distribution of each data index and goodness-of-fit analysis is carried out. The main disadvantage of the GTS distribution is the lack of the closed forms of the density, cumulative, and derivative functions. We use a computational algorithm, called the enhanced fast FRFT scheme (Nzokem 2023a), to circumvent the problem.

The rest of the paper is organized as follows: Section 2 provides some theoretical framework of the GTS distribution. Section 3 presents the multivariate maximum-likelihood (ML) method and the analytic version of the two-parameter normal distribution. Section 4 presents the results of the GTS parameter estimations, along with the associated statistical tests for the heavy-tailed Bitcoin and Ethereum returns. Section 5 fits the GTS distribution to the traditional indices S&P 500 and SPY ETF returns, while Section 6 presents the results of the goodness-of-fit test. Section 7 provides the concluding remarks.

2. Generalized Tempered Stable (GTS) Distribution

The Lévy measure of the GTS distribution ($V(dx)$) is defined in (2) as a product of a tempering function $q(x)$ and a Lévy measure of the $\alpha$-stable distribution $V_{stable}(dx)$:

$$\begin{array}{cc}\hfill q(x)& =e^{-{\lambda }_{+}x}{\mathbf{1}}_{x>0}+e^{-{\lambda }_{-}\left|x\right|}{\mathbf{1}}_{x<0}\hfill \\ \hfill V_{stable}(dx)& =\left({\alpha }_{+}{\displaystyle \frac{1}{x^{1+{\beta }_{+}}}}{\mathbf{1}}_{x>0}+{\alpha }_{-}{\displaystyle \frac{1}{{\left|x\right|}^{1+{\beta }_{-}}}}{\mathbf{1}}_{x<0}\right)dx\hfill \\ \hfill V(dx)=q(x)V_{stable}(dx)& =\left({\alpha }_{+}{\displaystyle \frac{e^{-{\lambda }_{+}x}}{x^{1+{\beta }_{+}}}}{\mathbf{1}}_{x>0}+{\alpha }_{-}{\displaystyle \frac{e^{-{\lambda }_{-}\left|x\right|}}{{\left|x\right|}^{1+{\beta }_{-}}}}{\mathbf{1}}_{x<0}\right)dx\hfill \end{array}$$

(2)

where $0\le {\beta }_{+}\le 1$, $0\le {\beta }_{-}\le 1$, ${\alpha }_{+}\ge 0$, ${\alpha }_{-}\ge 0$, ${\lambda }_{+}\ge 0$ and ${\lambda }_{-}\ge 0$.

The six parameters that appear have important interpretations. ${\beta }_{+}$ and ${\beta }_{-}$ are the indexes of stability bounded below by 0 and above by 2 (Borak et al. 2005). They capture the peakedness of the distribution similarly to the $\beta$-stable distribution, but the distribution tails are tempered. If $\beta$ increases (decreases), then the peakedness decreases (increases). ${\alpha }_{+}$ and ${\alpha }_{-}$ are the scale parameters, also called the process intensity (Boyarchenko and Levendorskii 2002); they determine the arrival rate of jumps for a given size. ${\lambda }_{+}$ and ${\lambda }_{-}$ control the decay rate on the positive and negative tails. Additionally, ${\lambda }_{+}$ and ${\lambda }_{-}$ are also skewness parameters. If ${\lambda }_{+}>{\lambda }_{-}$ (${\lambda }_{+}<{\lambda }_{-}$), then the distribution is skewed to the left (right), and if ${\lambda }_{+}={\lambda }_{-}$, then it is symmetric (Fallahgoul et al. 2019; Rachev et al. 2011). $\alpha$ and $\lambda$ are related to the degree of peakedness and thickness of the distribution. If $\alpha$ increases (decreases), the peakedness and the thickness decrease (increase). Similarly, if $\lambda$ increases (decreases), then the peakedness increases (decreases) and the thickness decreases (increases) (Bianchi et al. 2019). For more details on the tempering function and Lévy measure of the tempered stable distribution, refer to (Küchler and Tappe 2013; Rachev et al. 2011).

The activity process of the GTS distribution can be studied from the integral (3) of the Lévy measure (2):

$$\begin{array}{c}\hfill {\int }_{-\infty }^{+\infty }V(dx)=\left\{\begin{array}{c}+\infty \mathrm{if}0\le {\beta }_{+}<1\wedge 0\le {\beta }_{-}<1\hfill \\ {\alpha }_{+}{{\lambda }_{+}}^{{\beta }_{+}}\mathsf{\Gamma }(-{\beta }_{+})+{\alpha }_{-}{{\lambda }_{-}}^{{\beta }_{-}}\mathsf{\Gamma }(-{\beta }_{-})\mathrm{if}{\beta }_{+}<0\wedge {\beta }_{-}<0.\hfill \end{array}\right.\end{array}$$

(3)

As shown in (3), if ${\beta }_{+}<0$ and ${\beta }_{-}<0$, $GTS({\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$ is of a finite activity process and can be written as a compound Poisson (Barndorff-Nielsen and Shephard 2002). When $0\le {\beta }_{+}<1$ and $0\le {\beta }_{-}<1$, this Lévy density ($V(dx)$) is not integrable as it goes off to infinity too rapidly as x goes to zero (Barndorff-Nielsen and Shephard 2002), which means in practice that there will be a large number of very small jumps. As shown in (3), $GTS({\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$ is an infinite activity process with infinite jumps in any given time interval.

In addition to the infinite activities process, the variation in the process can be studied through the following integral:

$$\begin{array}{cc}\hfill {\int }_{-1}^1\left|x\right|V(dx)& ={\int }_{-1}^0\left|x\right|V(dx)+{\int }_0^1\left|x\right|V(dx)\hfill \\ \hfill & ={\alpha }_{-}{\lambda }_{-}^{{\beta }_{-}-1}\gamma (1-{\beta }_{-},{\lambda }_{-})+{\alpha }_{+}{\lambda }_{+}^{{\beta }_{+}-1}\gamma (1-{\beta }_{+},{\lambda }_{+})\hfill \end{array}$$

where $\gamma (s,x)={\int }_0^x{y}^{s-1}{e}^{-y}dy$ is the lower incomplete gamma function.

And we have:

$$\begin{array}{c}\hfill {\int }_{-1}^1\left|x\right|V(dx)<+\infty \mathrm{if}0<{\beta }_{-}\le 1\&0<{\beta }_{+}\le 1.\end{array}$$

(4)

As shown in (4), $GTS({\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$ generates a finite variance process, which is contrary to the Brownian motion process. $GTS({\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$ generates a type B Lévy process (Ken-Iti 2001), which is a purely non-Gaussian infinite activity Lévy process of finite variation whose sample paths have an infinite number of small jumps and a finite number of large jumps in any finite time interval.

The GTS distribution can be denoted by $X\sim GTS({\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$ and $X=X_{+}-X_{-}$ with $X_{+}\ge 0$, $X_{-}\ge 0$. $X_{+}\sim TS({\beta }_{+},{\alpha }_{+},{\lambda }_{+})$ and $X_{-}\sim TS({\beta }_{-},{\alpha }_{-},{\lambda }_{-})$. By adding the location parameter, the GTS distribution becomes $GTS(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$, and we have (5):

$$\begin{array}{c}\hfill Y=\mu +X=\mu +X_{+}-X_{-},Y\sim GTS(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-}).\end{array}$$

(5)

2.1. GTS Distribution and Characteristic Exponent

Theorem 1. 

Consider a variable $Y\sim GTS(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$. The characteristic exponent can be written as:

$$\begin{array}{c}\hfill \Psi (\xi )=\mu \xi i+{\alpha }_{+}\Gamma (-{\beta }_{+})\left({({\lambda }_{+}-i\xi )}^{{\beta }_{+}}-{{\lambda }_{+}}^{{\beta }_{+}}\right)+{\alpha }_{-}\Gamma (-{\beta }_{-})\left({({\lambda }_{-}+i\xi )}^{{\beta }_{-}}-{{\lambda }_{-}}^{{\beta }_{-}}\right).\end{array}$$

(6)

Proof. 

$V(dx)$ in (2) is a Lévy measure. The following relation is satisfied from (4):

$$\begin{array}{c}\hfill {\int }_{-\infty }^{+\infty }Min(1,|x|)V(dx)<+\infty .\end{array}$$

More details on the proof are provided in (Nzokem and Maposa 2024).

The Lévy–Khintchine representation (Barndorff-Nielsen and Shephard 2002) for non-negative Lévy process is applied on Y. $Y=\mu +X=\mu +X_{+}-X_{-}$ and we have:

$$\begin{array}{cc}\hfill \Psi (\xi )& =Log\left(E{e}^{iY\xi }\right)=i\mu \xi +Log\left(E{e}^{i{X}_{+}\xi }\right)+Log\left(E{e}^{-i{X}_{-}\xi }\right)\hfill \\ \hfill & =i\mu \xi +{\int }_0^{+\infty }\left(e^{iy\xi }-1\right){\displaystyle \frac{{\alpha }_{+}{e}^{-{\lambda }_{+}y}}{y^{1+{\beta }_{+}}}}dy+{\int }_0^{+\infty }\left(e^{-iy\xi }-1\right){\displaystyle \frac{{\alpha }_{-}{e}^{-{\lambda }_{-}y}}{y^{1+{\beta }_{-}}}}dy,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\int }_0^{+\infty }\left(e^{iy\xi }-1\right){\displaystyle \frac{{\alpha }_{+}{e}^{-{\lambda }_{+}y}}{y^{1+{\beta }_{+}}}}dy& ={\alpha }_{+}{\lambda }_{+}^{{\beta }_{+}}\Gamma (-{\beta }_{+})\sum _{k=1}^{+\infty }{\displaystyle \frac{\Gamma (k-{\beta }_{+})}{\Gamma (-{\beta }_{+})k!}}{({\displaystyle \frac{i\xi }{{\lambda }_{+}}})}^k\hfill \\ \hfill & ={\alpha }_{+}{\lambda }_{+}^{{\beta }_{+}}\Gamma (-{\beta }_{+})\sum _{k=1}^{+\infty }\left({\displaystyle \genfrac{}{}{0pt}{}{{\beta }_{+}}{k}}\right){(-{\displaystyle \frac{i\xi }{{\lambda }_{+}}})}^k\hfill \\ \hfill & ={\alpha }_{+}\Gamma (-{\beta }_{+})\left({({\lambda }_{+}-i\xi )}^{{\beta }_{+}}-{\lambda }_{+}^{{\beta }_{+}}\right).\hfill \end{array}$$

(8)

Similarly, we have:

$$\begin{array}{c}\hfill {\int }_0^{+\infty }\left(e^{-iy\xi }-1\right){\displaystyle \frac{{\alpha }_{-}{e}^{-{\lambda }_{-}y}}{y^{1+{\beta }_{-}}}}dy={\alpha }_{-}\Gamma (-{\beta }_{-})\left({({\lambda }_{-}+i\xi )}^{{\beta }_{-}}-{\lambda }_{-}^{{\beta }_{-}}\right).\end{array}$$

(9)

The expression in (7) becomes:

$$\begin{array}{c}\hfill \Psi (\xi )=i\mu \xi +{\alpha }_{+}\Gamma (-{\beta }_{+})\left({({\lambda }_{+}-i\xi )}^{{\beta }_{+}}-{\lambda }_{+}^{{\beta }_{+}}\right)+{\alpha }_{-}\Gamma (-{\beta }_{-})\left({({\lambda }_{-}+i\xi )}^{{\beta }_{-}}-{\lambda }_{-}^{{\beta }_{-}}\right).\end{array}$$

□

Theorem 2. 

Consider a variable $Y\sim GTS(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$.

If $({\beta }_{-},{\beta }_{+})\to (0,0)$, GTS becomes a bilateral Gamma distribution with the following characteristic exponent:

$$\begin{array}{c}\hfill \Psi (\xi )=\mu \xi i-{\alpha }_{+}log\left(1-{\displaystyle \frac{1}{{\lambda }_{+}}}i\xi \right)-{\alpha }_{-}log\left(1+{\displaystyle \frac{1}{{\lambda }_{-}}}i\xi \right).\end{array}$$

(10)

In addition to $({\beta }_{-},{\beta }_{+})\to (0,0)$, if ${\alpha }_{-}={\alpha }_{+}=\alpha$, GTS becomes a Variance-Gamma (VG) distribution with parameter $(\mu ,\delta ,\sigma ,\alpha ,\theta )$

$$\begin{array}{c}\hfill \delta ={\lambda }_{-}-{\lambda }_{+}\sigma =1\alpha ={\alpha }_{-}={\alpha }_{+}\theta ={\displaystyle \frac{1}{{\lambda }_{-}{\lambda }_{+}}}\end{array}$$

and the following characteristic exponent:

$$\begin{array}{c}\hfill \Psi (\xi )=\mu \xi i-\alpha log\left(1-{\displaystyle \frac{{\lambda }_{-}-{\lambda }_{+}}{{\lambda }_{+}{\lambda }_{-}}}i\xi +{\displaystyle \frac{1}{{\lambda }_{+}{\lambda }_{-}}}{\xi }^2\right).\end{array}$$

(11)

Proof. 

$$\begin{array}{cc}\hfill \Gamma (-{\beta }_{+})& =-{\displaystyle \frac{\Gamma (1-{\beta }_{+})}{{\beta }_{+}}}\hfill \\ \hfill \underset{{\beta }_{+}\to 0}{lim}\Gamma (-{\beta }_{+})\left({({\lambda }_{+}-i\xi )}^{{\beta }_{+}}-{\lambda }_{+}^{{\beta }_{+}}\right)& =-log\left(1-{\displaystyle \frac{1}{{\lambda }_{+}}}i\xi \right).\hfill \end{array}$$

(12)

Similarly, (12) works for ${\beta }_{-}\to 0$, and we have the characteristic exponent (10).

In addition, if ${\alpha }_{-}={\alpha }_{+}=\alpha$, from (10), the characteristic exponent becomes:

$$\begin{array}{c}\hfill \Psi (\xi )=\mu \xi i-\alpha log\left(1-{\displaystyle \frac{{\lambda }_{-}-{\lambda }_{+}}{{\lambda }_{+}{\lambda }_{-}}}i\xi +{\displaystyle \frac{1}{{\lambda }_{+}{\lambda }_{-}}}{\xi }^2\right),\end{array}$$

which is a Variance-Gamma (VG) distribution with parameter $(\mu ,{\lambda }_{-}-{\lambda }_{+},1,\alpha ,{\displaystyle \frac{1}{{\lambda }_{-}{\lambda }_{+}}})$. For more details on the VG model, refer to (Madan et al. 1998; Nzokem 2023c). □

Theorem 3 

(Cumulants ${\kappa }_k$). Consider a variable $Y\sim GTS(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$. The cumulants ${\kappa }_k$ of the GTS distribution are defined as follows:

$$\begin{array}{cc}\hfill {\kappa }_0& =0\hfill \\ \hfill {\kappa }_1& =\mu +{\alpha }_{+}{\displaystyle \frac{\Gamma (1-{\beta }_{+})}{{\lambda }_{+}^{1-{\beta }_{+}}}}-{\alpha }_{-}{\displaystyle \frac{\Gamma (1-{\beta }_{-})}{{\lambda }_{-}^{1-{\beta }_{-}}}}\hfill \\ \hfill {\kappa }_k& ={\alpha }_{+}{\displaystyle \frac{\Gamma (k-{\beta }_{+})}{{\lambda }_{+}^{k-{\beta }_{+}}}}+{(-1)}^k{\alpha }_{-}{\displaystyle \frac{\Gamma (k-{\beta }_{-})}{{\lambda }_{-}^{k-{\beta }_{-}}}}\forall k\in \mathbb{N}\backslash \{0,1\}.\hfill \end{array}$$

(13)

Proof. 

We reconsider the characteristic exponent $\Psi (\xi )$ in (7):

$$\begin{array}{cc}\hfill \Psi (\xi )& =i\mu \xi +{\int }_0^{+\infty }\left(e^{iy\xi }-1\right){\displaystyle \frac{{\alpha }_{+}{e}^{-{\lambda }_{+}y}}{y^{1+{\beta }_{+}}}}dy+{\int }_0^{+\infty }\left(e^{-iy\xi }-1\right){\displaystyle \frac{{\alpha }_{-}{e}^{{\lambda }_{-}y}}{y^{1+{\beta }_{-}}}}dy\hfill \\ \hfill & =i\mu \xi +{\alpha }_{+}\sum _{k=1}^{+\infty }{\displaystyle \frac{\Gamma (k-{\beta }_{+})}{{\lambda }_{+}^{k-{\beta }_{+}}}}{\displaystyle \frac{{(i\xi )}^k}{k!}}+{\alpha }_{-}\sum _{k=1}^{+\infty }{\displaystyle \frac{\Gamma (k-{\beta }_{-})}{{\lambda }_{-}^{k-{\beta }_{-}}}}{\displaystyle \frac{{(-i\xi )}^k}{k!}}\hfill \\ \hfill & =i\mu \xi +\sum _{k=1}^{+\infty }{\displaystyle \frac{1}{k!}}\left({\alpha }_{+}{\displaystyle \frac{\Gamma (k-{\beta }_{+})}{{\lambda }_{+}^{k-{\beta }_{+}}}}+{\alpha }_{-}{\displaystyle \frac{\Gamma (k-{\beta }_{-})}{{\lambda }_{-}^{k-{\beta }_{-}}}}{(-1)}^k\right){(i\xi )}^k\hfill \\ \hfill & =\sum _{k=0}^{+\infty }{\displaystyle \frac{{\kappa }_k}{k!}}{(i\xi )}^k.\hfill \end{array}$$

(14)

Hence, the k-th order cumulant ${\kappa }_k$ is given by comparing the coefficients of both polynomial functions in $i\xi$. For more details on the relationship between the characteristic exponent and cumulant functions, refer to (Feller 1971; Kendall 1945). □

2.2. GTS Distribution and Lévy Process

Corollary 1. 

Let $Y=\left(Y_t\right)$ be a Lévy process on ${\mathbb{R}}^{+}$ generated by $GTS(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$, and then

$$\begin{array}{c}\hfill Y_t\sim GTS(t\mu ,{\beta }_{+},{\beta }_{-},t{\alpha }_{+},t{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})\forall t\in {\mathbb{R}}^{+}.\end{array}$$

(15)

Proof. 

Let $\Psi (\xi ,t)$ be the characteristic exponent of the Lévy process $Y=\left(Y_t\right)$. By applying the infinitely divisible property, we have:

$$\begin{array}{cc}\hfill \Psi (\xi ,t)& =Log\left(E{e}^{i{Y}_t\xi }\right)=tLog\left(E{e}^{iX\xi }\right)\hfill \\ \hfill & =t\mu \xi i+t{\alpha }_{+}\Gamma (-{\beta }_{+})\left({({\lambda }_{+}-i\xi )}^{{\beta }_{+}}-{\lambda }_{+}^{{\beta }_{+}}\right)+t{\alpha }_{-}\Gamma (-{\beta }_{-})\left({({\lambda }_{-}+i\xi )}^{{\beta }_{-}}-{\lambda }_{-}^{{\beta }_{-}}\right)\hfill \end{array}$$

and we deduce that $Y_t\sim GTS(t\mu ,{\beta }_{+},{\beta }_{-},t{\alpha }_{+},t{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$. □

Theorem 4 

(Asymptotic distribution of Generalized Tempered Stable distribution process). Let $Y=Y_t$ be a Lévy process on $\mathbb{R}$ generated by $GTS(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$. Then, $Y_t$ converges in distribution to a Lévy process driving by a normal distribution with mean ${\kappa }_1$ and variance ${\kappa }_2$

$$\begin{array}{cc}\hfill Y_t\stackrel{d}{\to }N(t{\kappa }_1,t{\kappa }_2)& \mathrm{as}t\to +\infty \hfill \end{array}$$

(16)

where

$$\begin{array}{cc}\hfill {\kappa }_1& =\mu +{\alpha }_{+}{\displaystyle \frac{\Gamma (1-{\beta }_{+})}{{\lambda }_{+}^{1-{\beta }_{+}}}}-{\alpha }_{-}{\displaystyle \frac{\Gamma (1-{\beta }_{-})}{{\lambda }_{-}^{1-{\beta }_{-}}}}\hfill \\ \hfill {\kappa }_2& ={\alpha }_{+}{\displaystyle \frac{\Gamma (2-{\beta }_{+})}{{\lambda }_{+}^{2-{\beta }_{+}}}}+{\alpha }_{-}{\displaystyle \frac{\Gamma (2-{\beta }_{-})}{{\lambda }_{-}^{2-{\beta }_{-}}}}.\hfill \end{array}$$

Proof. 

The proof relies on the cumulant-generating function. As in (14), the characteristic exponent ($\Psi (\xi )$) can be written as follows:

$$\begin{array}{c}\hfill \Psi (\xi )=Log\left(E{e}^{iY\xi }\right)=\sum _{j=0}^{+\infty }{\kappa }_j{\displaystyle \frac{{(i\xi )}^j}{j!}}.\end{array}$$

(17)

Let $\varphi (\xi ,t)$ be the characteristic function of the stochastic process ${\displaystyle \frac{Y_t-t{\kappa }_1}{\sqrt{t{\kappa }_2}}}$ and we have:

$$\begin{array}{cc}\hfill \varphi (\xi ,t)& =E\left(e^{i{\displaystyle \frac{Y_t-t{\kappa }_1}{\sqrt{t{\kappa }_2}}}\xi }\right)=e^{-i{\displaystyle \frac{t{\kappa }_1}{\sqrt{t{\kappa }_2}}}\xi }E\left(e^{i{\displaystyle \frac{\xi }{\sqrt{t{\kappa }_2}}}{Y}_t}\right)\hfill \\ \hfill & =e^{-i{\displaystyle \frac{t{\kappa }_1}{\sqrt{t{\kappa }_2}}}\xi }{e}^{t\Psi ({\displaystyle \frac{\xi }{\sqrt{t{\kappa }_2}}})}=e^{-i{\displaystyle \frac{t{\kappa }_1}{\sqrt{t{\kappa }_2}}}\xi }{e}^{{\sum }_{j=0}^{+\infty }{\displaystyle \frac{t{\kappa }_j}{j!}}{(i{\displaystyle \frac{\xi }{\sqrt{t{\kappa }_2}}})}^j}\hfill \\ \hfill & =e^{-{\displaystyle \frac{{\xi }^2}{2}}+{\sum }_{j=3}^{+\infty }{\displaystyle \frac{t{\kappa }_j}{j!}}{(i{\displaystyle \frac{\xi }{\sqrt{t{\kappa }_2}}})}^j},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \underset{t\to +\infty }{lim}\sum _{j=3}^{+\infty }{\displaystyle \frac{t{\kappa }_j}{j!}}{(i{\displaystyle \frac{\xi }{\sqrt{t{\kappa }_2}}})}^j=0\underset{t\to +\infty }{lim}\varphi (\xi ,t)& =\underset{t\to +\infty }{lim}{e}^{-{\displaystyle \frac{{\xi }^2}{2}}+{\sum }_{j=3}^{+\infty }{\displaystyle \frac{t{\kappa }_j}{j!}}{(i{\displaystyle \frac{\xi }{\sqrt{t{\kappa }_2}}})}^j}\hfill \\ \hfill & =e^{-{\displaystyle \frac{1}{2}}{\xi }^2}.\hfill \end{array}$$

(19)

□

3. Multivariate Maximum-Likelihood Method

3.1. Maximum-Likelihood Method: Numerical Approach

From a probability density function $f(x,V)$ with parameter $V=(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$ and sample data $x={\left(x_j\right)}_{1\le j\le m}$, we define the likelihood function and its first and second derivatives as follows:

$$\begin{array}{cc}\hfill L_m(x,V)& =\prod _{j=1}^{m}f(x_j,V),l_m(x,V)=\sum _{j=1}^{m}log(f(x_j,V))\hfill \\ \hfill {\displaystyle \frac{d{l}_m(x,V)}{d{V}_j}}& =\sum _{i=1}^m{\displaystyle \frac{\frac{df(x_i,V)}{d{V}_j}}{f(x_i,V)}}\hfill \\ \hfill {\displaystyle \frac{d^2{l}_m(x,V)}{d{V}_{k}d{V}_j}}& =\sum _{i=1}^m\left({\displaystyle \frac{\frac{d^{2}f(x_i,V)}{d{V}_{k}d{V}_j}}{f(x_i,V)}}-{\displaystyle \frac{\frac{df(x_i,V)}{d{V}_k}}{f(x_i,V)}}{\displaystyle \frac{\frac{df(x_i,V)}{d{V}_j}}{f(x_i,V)}}\right).\hfill \end{array}$$

(20)

To perform the maximum of the likelihood function ($L_m(x,V)$), we need the gradient of the likelihood function (${\displaystyle \frac{d{l}_m(x,V)}{dV}}$), also known as the score function, and the Hessian matrix (${\displaystyle \frac{d^2{l}_m(x,V)}{dV{dV}^{\prime }}}$), which is the variance–covariance matrix generated by the likelihood function.

Given the parameters $V=(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$ and the sample data set X, we have the following quantities (21) from the previous development:

$$\begin{array}{c}\hfill I_m^{\prime }(X,V)={\left({\displaystyle \frac{d{l}_m(x,V)}{d{V}_j}}\right)}_{1\le j\le p},I_m^{″}(X,V)={\left({\displaystyle \frac{d^2{l}_m(x,V)}{d{V}_{k}d{V}_j}}\right)}_{\begin{array}{c}1\le k\le p\\ 1\le j\le p\end{array}}.\end{array}$$

(21)

We use a computational algorithm built as a composite of a standard FRFT to compute the likelihood function and its derivatives (20) in the optimization process. More details on applying the composite of FRFTs for parameter estimations are provided in (Nzokem 2021b, 2021c; Nzokem and Montshiwa 2022, 2023); for other computations (such as probability density and cumulative functions), see (Cherubini et al. 2010; Eberlein 2014; Eberlein et al. 2010; Nzokem 2023b; Nzokem and Maposa 2024).

The computational algorithm yields a local solution, V, and a negative semi-definite matrix, $I_m^{″}(x,V)$, when the following two conditions are satisfied:

$$\begin{array}{c}\hfill I_m^{\prime }(x,V)=0,U^T{\mathbf{I}}_{\mathbf{m}}^{″}(\mathbf{X},\mathbf{V})U\le 0,\forall U\in {\mathbb{R}}^p.\end{array}$$

(22)

The solutions, V, in (22) are provided by the Newton–Raphson iteration algorithm Formula (23):

$$\begin{array}{c}\hfill V^{n+1}=V^n-{\left(I_m^{″}(x,V^n)\right)}^{-1}{I}_m^{\prime }(x,V^n).\end{array}$$

(23)

More details on the maximum-likelihood and Newton–Raphson iteration procedures are provided in (Giudici et al. 2013).

3.2. Asymptotic Distribution of the Maximum-Likelihood Estimator (MLE)

Theorem 5 

(Cramer-Rao). Let $T=T(X_1,\dots ,X_m)$ be a statistic and write $E\left[T\right]=k(\theta )$. Then, under suitable (smoothness) assumptions,

$$\begin{array}{c}\hfill Var\left[T\right]\ge {\displaystyle \frac{{(\frac{dE\left[T\right]}{d\theta })}^2}{mI(\theta )}}.\end{array}$$

(24)

For the proof of Theorem 5, refer to (Casella and Berger 2024; Van den Bos 2007).

Theorem 6 

(Consistency Estimator). Let $X_1,\dots ,X_m$ be independent and identically distributed (i.i.d) random variables with density $f(x|\theta )$ satisfying some regularity conditions (Lehmann 1999). Let θ be the true parameter; then, there exists a sequence ${\widehat{\theta }}_m={\theta }_m(X_1,\dots ,X_m)$ of local maxima of the likelihood function $L_m(\theta )$ which is consistent, that is, which satisfies

$$\begin{array}{cc}\hfill {\widehat{\theta }}_m\stackrel{a.s.}{\to }\theta & \mathrm{as}m\to +\infty .\hfill \end{array}$$

(25)

More details on the proof of Theorem 6 are provided in (Casella and Berger 2024; Lehmann 1999).

Theorem 7 

(Asymptotic Efficiency and Normality). Let $X_1,\dots ,X_m$ be independent and identically distributed (i.i.d) random variables with density $f(x|\theta )$ satisfying some regularity conditions in (Lehmann 1999). There exists a solution ${\widehat{\theta }}_m={\theta }_m(X_1,\dots ,X_m)$ of the likelihood equations which is consistent, and any such solution satisfies:

$$\begin{array}{cc}\hfill {\widehat{\theta }}_m-\theta \stackrel{d}{\to }N\left(0,I_m^{-1}(\theta )\right)& \mathrm{as}m\to +\infty ,\hfill \end{array}$$

(26)

where $\theta =({\theta }_1,\dots ,{\theta }_k)$ is the actual parameter and $I_m(\theta )$ is the Fisher information matrix.

More details on the proof of Theorem 7 are provided in (Hall and Oakes 2023; Lehmann 1999; Olive 2014).

Theorem 8 

(Likelihood Ratio Test). Suppose the assumptions of Theorem 7 hold and that $({\widehat{\theta }}_{1n},\dots ,{\widehat{\theta }}_{kn})$ are consistent roots of the likelihood equations for $\theta =({\theta }_1,\dots ,{\theta }_k)$. In addition, suppose that the corresponding assumptions hold for the parameter vector $({\theta }_{r+1},\dots ,{\theta }_k)$ when $r<k$ and that ${\widehat{\widehat{\theta }}}_{r+1,n},\dots ,{\widehat{\widehat{\theta }}}_{kn}$ are consistent roots of the likelihood equations for $({\theta }_{r+1},\dots ,{\theta }_k)$ under the null hypothesis. We consider the likelihood ratio statistic

$$\begin{array}{c}\hfill {\displaystyle \frac{l_m(x,\widehat{\theta })}{l_m(x,\widehat{\widehat{\theta }})}}\end{array}$$

(27)

where $\widehat{\widehat{\theta }}=({\theta }_1,\dots ,{\theta }_r,{\widehat{\widehat{\theta }}}_{r+1,n},\dots ,{\widehat{\widehat{\theta }}}_{kn})$. Then under the null hypothesis $H_0$, if

$$\begin{array}{c}\hfill {\Delta }_n=l_m(x,\widehat{\theta })-l_m(x,\widehat{\widehat{\theta }}),\end{array}$$

(28)

the statistic $2{\Delta }_n$ has a limiting ${\chi }_r^2$ distribution.

More details on the proof of Theorem 8 are provided in (Lehmann 1999; Vuong 1989).

3.3. Asymptotic Test and Confidence Interval

The above results allow us to construct an asymptotically efficient estimator ${\widehat{\theta }}_m=({\widehat{\theta }}_{1m},\dots ,{\widehat{\theta }}_{km})$ of $\theta =({\theta }_1,\dots ,{\theta }_k)$ such that

$$\begin{array}{c}\hfill ({\widehat{\theta }}_{1m}-{\theta }_1,\dots ,{\widehat{\theta }}_{km}-{\theta }_k)\end{array}$$

(29)

has a joint multivariate limit distribution with mean $(0,\dots ,0)$ and covariance matrix $I_m^{-1}(\theta )=(J_{ij})$. In particular, we have:

$$\begin{array}{c}\hfill {\widehat{\theta }}_{jm}-{\theta }_j\stackrel{d}{\to }N(0,J_{jj})\mathrm{as}m\to +\infty .\end{array}$$

(30)

One approach to constructing an asymptotically valid confidence interval for the parameters is via the asymptotic distribution of the ML estimator (27). An approximate $(1-{\displaystyle \frac{\alpha }{2}})$ confidence interval for ${\widehat{\theta }}_{jm}$ can be written as follows:

$$\begin{array}{c}\hfill {\widehat{\theta }}_{jm}\pm z({\displaystyle \frac{\alpha }{2}})*\sqrt{J_{jj}}\mathrm{as}m\to +\infty ,\end{array}$$

(31)

where $z({\displaystyle \frac{\alpha }{2}})$ is the ${\displaystyle \frac{\alpha }{2}}$ quantile of the standard normal distribution.

3.4. Applications of the Log-Likelihood Estimator to the Normal Distribution

We suppose the sample data $x={\left(x_j\right)}_{1\le j\le m}$ are independent observations and have a normal distribution (Mensah et al. 2023) with parameter $V(\mu ,{\sigma }^2)$, that is, $y\sim \mathcal{N}(\mu ,{\sigma }^2)$; then, the density is

$$\begin{array}{c}\hfill f(y|V)={(2\pi {\sigma }^2)}^{-{\displaystyle \frac{1}{2}}}exp\left(-{\displaystyle \frac{{(y-\mu )}^2}{2{\sigma }^2}}\right).\end{array}$$

(32)

The log-likelihood function in (20) becomes

$$\begin{array}{c}\hfill l_m(x|V)=\sum _{j=1}^{m}log(f(x_j|V))=-{\displaystyle \frac{m}{2}}log(2\pi {\sigma }^2)-{\displaystyle \frac{1}{2{\sigma }^2}}\sum _{j=1}^m{(x_j-\mu )}^2.\end{array}$$

(33)

The first-order derivatives of the log-likelihood function with respect to $\mu$ and ${\sigma }^2$ in (20) become

$$\begin{array}{c}\hfill I_m^{\prime }(X,V)=\left(\begin{array}{c}{\displaystyle \frac{d{l}_m(x,V)}{d\mu }}\\ {\displaystyle \frac{d{l}_m(x,V)}{d{\sigma }^2}}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1}{{\sigma }^2}}\sum _{j=1}^m(x_j-\mu )\\ {\displaystyle \frac{1}{2{\sigma }^4}}\sum _{j=1}^m{(x_j-\mu )}^2-{\displaystyle \frac{m}{2{\sigma }^2}}.\end{array}\right)\end{array}$$

(34)

By setting $I_m^{\prime }(X,V)=0$, we have

$$\begin{array}{c}\hfill \widehat{\mu }={\displaystyle \frac{1}{m}}\sum _{j=1}^m{x}_j{\widehat{\sigma }}^2={\displaystyle \frac{1}{m}}\sum _{j=1}^m{(x_j-\widehat{\mu })}^2.\end{array}$$

(35)

The second-order derivative of the log-likelihood function with respect to $\mu$ and ${\sigma }^2$ in (20) becomes

$$\begin{array}{cc}\hfill I_m^{″}(X,V)& =\left(\begin{array}{cc}{\displaystyle \frac{d^2{l}_m(x,V)}{d{\mu }^2}}& {\displaystyle \frac{d{l}_m(x,V)}{d\mu d{\sigma }^2}}\\ {\displaystyle \frac{d{l}_m(x,V)}{d{\sigma }^{2}d\mu }}& {\displaystyle \frac{d^2{l}_m(x,V)}{{(d{\sigma }^2)}^2}}\end{array}\right)\hfill \\ \hfill & =\left(\begin{array}{cc}-{\displaystyle \frac{m}{{\sigma }^2}}& -{\displaystyle \frac{{\sum }_{j=1}^m(x_j-\mu )}{{\sigma }^4}}\\ -{\displaystyle \frac{{\sum }_{j=1}^m(x_j-\mu )}{2{\sigma }^4}}& -{\displaystyle \frac{1}{{\sigma }^6}}\sum _{j=1}^m{(x_j-\mu )}^2+{\displaystyle \frac{m}{2{\sigma }^4}}.\end{array}\right)\hfill \end{array}$$

(36)

Refer to (Casella and Berger 2024) for more details.

We have the Fisher information matrix and the inverse:

$$\begin{array}{c}\hfill I_m(V)=-E\left(I_m^{″}(X,V)\right)=\left(\begin{array}{cc}{\displaystyle \frac{m}{{\sigma }^2}}& 0\\ 0& {\displaystyle \frac{m}{2{\sigma }^4}},\end{array}\right)I_m^{-1}(V)=\left(\begin{array}{cc}{\displaystyle \frac{{\sigma }^2}{m}}& 0\\ 0& {\displaystyle \frac{2{\sigma }^4}{m}}\end{array}\right).\end{array}$$

(37)

Corollary 2. 

The limiting distribution of the MLE is given by:

$$\begin{array}{c}\hfill \left(\begin{array}{c}\widehat{\mu }\\ {\widehat{\sigma }}^2\end{array}\right)\stackrel{d}{\to }N\left(\left(\begin{array}{c}\mu \\ {\sigma }^2\end{array}\right),\left(\begin{array}{cc}{\displaystyle \frac{{\sigma }^2}{m}}& 0\\ 0& {\displaystyle \frac{2{\sigma }^4}{m}}\end{array}\right)\right),\mathrm{as}m\to +\infty .\end{array}$$

(38)

The proof of Corollary 2 comes from Theorem 7, Equation (26).

4. Fitting Tempered Stable Distribution to Cryptocurrencies: Bitcoin (BTC) and Ethereum

4.1. Data Summaries

Bitcoin was the first cryptocurrency created in 2009 by Satoshi Nakamoto. The idea behind Bitcoin was to create a peer-to-peer electronic payment system that allows online payments to be sent directly from one party to another without going through a financial institution (Nakamoto 2008). Since its inception, Bitcoin has grown in popularity and adoption and is now viewed as a viable legal tender in some countries. Bitcoin is currently used more as an investment tool, a risk-diversified tool, and less as a medium of exchange, a store of value, or a unit of account (Nzokem and Maposa 2024).

Bitcoin (BTC) and Ethereum (ETH) prices were extracted from CoinMarketCap. The period spans from 28 April 2013 to 4 July 2024 for Bitcoin and from 7 August 2015 to 4 July 2024 for Ethereum.

The daily price dynamics are provided in Figure 1. The prices have an increasing trend, even after having major significant increases and decreases over the studied period. Figure 1a,b show that Bitcoin outperforms Ethereum, which is the second-largest cryptocurrency by market capitalization after Bitcoin.

Let m be the number of observations and $S_j$ be the daily observed price on the day $t_j$ with $j=1,\dots ,m$. The daily return $(y_j)$ is computed as follows:

$$\begin{array}{c}\hfill y_j=log(S_j/S_{j-1})j=2,\dots ,m.\end{array}$$

(39)

As shown in Figure 2a,b, the daily return reaches the lowest level ($-46\%$ for Bitcoin and $-55\%$ for Ethereum) in the first quarter of 2020 amid the coronavirus pandemic and massive disruptions in the global economy. Nine values were identified as outliers and removed from the dataset to avoid a negative impact on the GTS model estimation and the empirical statistics.

4.2. Multidimensional Estimation Results for Cryptocurrencies

The results of the GTS parameter estimation are summarized in

Table 1. Maximum-likelihood GTS parameter estimation for Bitcoin.

Model	Parameter	Estimate	Std Err	z	$\mathit{Pr}(\mathit{Z}>|\mathit{z}|)$	[95% Conf.Interval]
GTS	$\mu$	−0.121571	(0.375)	−0.32	7.5 × ${10}^{-01}$	−0.856	0.613
	${\beta }_{+}$	0.315548	(0.136)	2.33	2.0 × ${10}^{-02}$	0.050	0.581
	${\beta }_{-}$	0.406563	(0.117)	3.48	4.9 × ${10}^{-04}$	0.178	0.635
	${\alpha }_{+}$	0.747714	(0.047)	15.76	6.2 × ${10}^{-56}$	0.655	0.841
	${\alpha }_{-}$	0.544565	(0.037)	14.56	4.8 × ${10}^{-48}$	0.471	0.618
	${\lambda }_{+}$	0.246530	(0.036)	6.91	4.9 × ${10}^{-12}$	0.177	0.316
	${\lambda }_{-}$	0.174772	(0.026)	6.69	2.2 × ${10}^{-11}$	0.124	0.226
	Log(ML)	−10,606
	AIC	21,227
	BIK	21,271
GBM	$\mu$	0.151997	(0.060)	2.51	1.2 × ${10}^{-02}$	0.033	0.271
	$\sigma$	3.865132	(0.330)	11.69	7.2 × ${10}^{-32}$	3.217	4.513
	Log(ML)	−11,313
	AIC	22,630
	BIK	22,638

for Bitcoin and

Table 2. Maximum-likelihood GTS parameter estimation for Ethereum.

Model	Param	Estimate	Std Err	z	$\mathit{P}\mathit{r}(\mathit{Z}>|\mathit{z}|)$	[95% Conf.Interval]
GTS	$\mu$	−0.4854	(1.008)	−0.48	6.3 × ${10}^{-01}$	−2.461	1.491
	${\beta }_{+}$	0.3904	(0.164)	2.38	1.7 × ${10}^{-02}$	0.069	0.712
	${\beta }_{-}$	0.4045	(0.210)	1.93	5.4 × ${10}^{-02}$	−0.007	0.816
	${\alpha }_{+}$	0.9582	(0.106)	9.01	1.1 × ${10}^{-19}$	0.750	1.167
	${\alpha }_{-}$	0.8005	(0.110)	7.25	4.2 × ${10}^{-13}$	0.584	1.017
	${\lambda }_{+}$	0.1667	(0.029)	5.72	1.1 × ${10}^{-08}$	0.110	0.224
	${\lambda }_{-}$	0.1708	(0.036)	4.71	2.5 × ${10}^{-06}$	0.110	0.242
	Log(ML)	−9552
	AIC	19,119
	BIK	19,162
GBM	$\mu$	0.267284	(0.091)	2.93	3.4 × ${10}^{-03}$	0.088	0.446
	$\sigma$	5.205539	(0.672)	7.74	1.0 × ${10}^{-14}$	3.887	6.524
	Log(ML)	−9960
	AIC	19,925
	BIK	19,933

for Ethereum data. The brackets are the asymptotic standard errors computed using the inverse of the Hessian matrix built in (20). The ML estimate of $\mu$ is negative for both Bitcoin and Ethereum, while others are positive, as expected in the literature. The asymptotic standard error for $\mu$ is quite large and suggests that $\mu$ is not statistically significant at 5%.

The log-likelihood, Akaike’s information Criteria (AIC), and Bayesian information criteria (BIK) statistics show that the GTS distribution with seven parameters performs better than the two-parameter normal distribution (GBM). A comprehensive and detailed examination of the statistical significance of the results will be carried out in Section 6.

Table 1 summarizes the estimation results for Bitcoin returns. The skewness parameters (${\lambda }_{+}$, ${\lambda }_{-}$) are statistically significant at 5%. The difference is positive and statistically significant, which proves that the Bitcoin return is asymmetric and skewed to the left. The process intensity parameters (${\alpha }_{+}$, ${\alpha }_{-}$) are statistically significant at 5%. Similarly, the difference is positive and statistically significant, showing that the Bitcoin is more likely to produce positive returns than negative ones. The index of stability parameters (${\beta }_{+}$, ${\beta }_{-}$) are both statistically significant at 5%. However, the difference is positive but not statistically significant.

The GTS distribution with $\beta ={\beta }_{+}={\beta }_{-}$, called the Kobol distribution, was fitted to the Bitcoin data as well, and the estimation results are presented in Appendix B.1. As shown in

Table A6. Kobol maximum-likelihood estimation for Bitcoin return data.

Model	Parameter	Estimate	Std Err	z	$\mathit{Pr}(\mathit{Z}>|\mathit{z}|)$	[95% Conf.Interval]
GTS	$\mu$	−0.292833	(0.126)	−2.32	2.1 × ${10}^{-02}$	−0.541	−0.045
	$\beta$	0.367074	(0.086)	4.27	1.9 × ${10}^{-05}$	0.199	0.535
	${\alpha }_{+}$	0.755914	(0.047)	16.02	4.7 × ${10}^{-58}$	0.663	0.848
	${\alpha }_{-}$	0.535121	(0.034)	15.68	9.6 × ${10}^{-56}$	0.468	0.602
	${\lambda }_{+}$	0.235266	(0.027)	8.87	3.6 × ${10}^{-19}$	0.183	0.287
	${\lambda }_{-}$	0.181602	(0.023)	7.94	9.8 × ${10}^{-16}$	0.137	0.226
	Log(ML)	−10,607
	AIC	21,226
	BIK	21,264

, all the parameters are statistically significant at 5%, and have the expected positive sign. However, the likelihood ratio test in Table 6 shows that the GTS distribution in Table 1 is not significantly different from the Kobol distribution as the p-value ($69.6\%$) is large. Refer to (Boyarchenko and Levendorskii 2002) for more details on the Kobol distribution.

As shown in Table 2, the parameters for Ethereum returns data are statistically significant at 5%, except $\mu$ and ${\beta }_{-}$. The difference (${\lambda }_{+}-{\lambda }_{-}$) in skewness parameters is negative and not statistically significant, showing that the Ethereum return is asymmetric and skewed to the right. Similarly, the difference (${\alpha }_{+}-{\alpha }_{-}$) in the intensity parameters is positive and not statistically significant, as shown the confidence interval. Contrary to the Bitcoin return, the Ethereum return has a larger process intensity, which provides evidence that Ethereum has a lower level of peakedness and a higher level of thickness.

We consider the following constraints $\lambda ={\lambda }_{+}={\lambda }_{-}$ and $\beta ={\beta }_{+}={\beta }_{-}$, which are the Carr–Geman–Madan–Yor (CGMY) distribution, also called the Classical Tempered Stable Distribution. The CGMY distribution was fitted as well, and the estimation results are presented in Appendix B.2. As shown in

Table A8. CGMY maximum-likelihood estimation for Ethereum return data.

Model	Parameter	Estimate	Std Err	z	$\mathit{Pr}(\mathit{Z}>|\mathit{z}|)$	[95% Conf.Interval]
GTS	$\mu$	−0.147089	(0.079)	−1.86	6.3 × ${10}^{-02}$	−0.302	−0.008
	$\beta$	0.398418	(0.127)	3.12	1.8 × ${10}^{-03}$	0.148	0.649
	$\alpha$	0.887161	(0.058)	15.22	1.2 × ${10}^{-52}$	0.773	1.001
	${\lambda }_{+}$	0.155369	(0.023)	6.56	5.2 × ${10}^{-11}$	0.109	0.202
	${\lambda }_{-}$	0.185991	(0.025)	7.29	2.9 × ${10}^{-13}$	0.136	0.236
	Log(ML)	−9554
	AIC	19,118
	BIK	19,149

, all the parameters are statistically significant at 5%, and have the expected positive sign. However, the likelihood ratio test in Table 6 shows, with a high p-value ($35.3\%$), that the GTS distribution is not significantly different from the CGMY distribution, and the null hypothesis cannot be rejected. Refer to (Carr et al. 2003; Rachev et al. 2011) for more details on the CGMY distribution.

Table 1 and Table 2 summarize the last row of

Table A1. Convergence of the GTS parameter for Bitcoin return data.

$\mathit{Iterations}$	$\mathit{\mu }$	${\mathit{\beta }}_{+}$	${\mathit{\beta }}_{-}$	${\mathit{\alpha }}_{+}$	${\mathit{\alpha }}_{-}$	${\mathit{\lambda }}_{+}$	${\mathit{\lambda }}_{-}$	$\mathit{Log}(\mathit{ML})$	$\left|\right|\frac{\mathit{dLog}(\mathit{ML})}{\mathit{dV}}\left|\right|$	$\mathit{MaxEigenValue}$
1	−0.7369246	0.4613783	0.2671787	0.8100173	0.5173470	0.2156289	0.1919378	−10609.058	282.6765666	3.6240151
2	−0.7977019	0.4654390	0.2169392	0.7846817	0.4905332	0.2164395	0.2049523	−10607.253	26.7522215	−1.6194299
3	−0.4455841	0.3884721	0.3213867	0.7758150	0.5193395	0.2340187	0.1883953	−10607.001	50.1355291	3.0916011
4	−0.7634445	0.4521878	0.2217702	0.7935129	0.4959371	0.2218253	0.2055181	−10607.210	4.8235882	−2.6390063
5	−0.4906746	0.4146531	0.3404176	0.7722729	0.5222110	0.2269202	0.1846457	−10607.059	67.6646338	9.0971871
6	−0.5515834	0.4434827	0.3335905	0.7724484	0.5190619	0.2197566	0.1853686	−10607.022	17.4476962	−0.4021102
7	−0.4914586	0.4327714	0.3503012	0.7686883	0.5235361	0.2216450	0.1826269	−10606.991	16.2838831	−0.1781480
8	−0.2900908	0.3885350	0.3956186	0.7563357	0.5370260	0.2300772	0.1754994	−10606.864	12.0116477	−2.4090216
9	−0.2752698	0.3832660	0.3969704	0.7555456	0.5377367	0.2312224	0.1753571	−10606.847	11.4457840	−2.5487401
10	−0.2609339	0.3780400	0.3982456	0.7547812	0.5384209	0.2323632	0.1752258	−10606.832	10.8628213	−2.6874876
11	−0.2085409	0.3576927	0.4025762	0.7519966	0.5408864	0.2368544	0.1748113	−10606.782	8.3600783	−3.4356438
12	−0.1970109	0.3528575	0.4034002	0.7513923	0.5414138	0.2379362	0.1747455	−10606.772	7.6818408	−3.6954428
13	−0.1761733	0.3436416	0.4046818	0.7503191	0.5423414	0.2400174	0.1746675	−10606.756	6.2516380	−4.2766527
14	−0.1668421	0.3392794	0.4051522	0.7498492	0.5427438	0.2410120	0.1746529	−10606.750	5.5002876	−4.5807361
15	−0.1581860	0.3350854	0.4055256	0.7494209	0.5431090	0.2419740	0.1746517	−10606.745	4.7262048	−4.8824015
16	−0.1501600	0.3310612	0.4058166	0.7490311	0.5434404	0.2429024	0.1746615	−10606.741	3.9306487	−5.1742197
17	−0.1209376	0.3159301	0.4066945	0.7476393	0.5446224	0.2464122	0.1747326	−10606.734	2.8592342	−6.2251311
18	−0.1216487	0.3155707	0.4064438	0.7477179	0.5445608	0.2465247	0.1747753	−10606.734	0.0014787	−6.2014232
19	−0.1215714	0.3155483	0.4064635	0.7477142	0.5445652	0.2465296	0.1747719	−10606.734	1.82 × ${10}^{-06}$	−6.2026532
20	−0.1215714	0.3155483	0.4064635	0.7477142	0.5445652	0.2465296	0.1747719	−10606.734	9.80 × ${10}^{-10}$	−6.2026530

and

Table A2. Convergence of the GTS parameter for Ethereum return data.

$\mathit{Iterations}$	$\mathit{\mu }$	${\mathit{\beta }}_{+}$	${\mathit{\beta }}_{-}$	${\mathit{\alpha }}_{+}$	${\mathit{\alpha }}_{-}$	${\mathit{\lambda }}_{+}$	${\mathit{\lambda }}_{-}$	$\mathit{Log}(\mathit{ML})$	$\left|\right|\frac{\mathit{dLog}(\mathit{ML})}{\mathit{dV}}\left|\right|$	$\mathit{MaxEigenValue}$
1	−0.1215714	0.3155483	0.4064635	0.7477142	0.5445652	0.2465296	0.1747719	−9745.171	2673.428257	206.013602
2	−0.1724835	0.3319505	0.4091022	0.7364129	0.5479934	0.2227870	0.1684568	−9700.715	2388.609394	180.884105
3	−0.2041418	0.3384742	0.4118929	0.7338794	0.5531083	0.2083203	0.1632896	−9669.986	2139.267660	157.699659
4	−0.4006157	0.3530035	0.4393474	0.7513784	0.6172425	0.1135743	0.1221930	−9586.115	1471.570475	32.410140
5	−0.6485551	0.4493817	0.4404508	0.9247887	0.7210031	0.1412949	0.1482307	−9556.026	380.605737	56.584055
6	−0.6290525	0.4371402	0.4359516	0.9780784	0.7824777	0.1582340	0.1608694	−9553.005	24.905322	−0.719221
7	−0.5545412	0.3994778	0.3918188	0.9627486	0.7936571	0.1652438	0.1724287	−9552.866	5.834338	−0.847574
8	−0.4744837	0.3913982	0.4093404	0.9582366	0.8022858	0.1665103	0.1699928	−9552.862	2.963350	−0.933466
9	−0.4825586	0.3902160	0.4051365	0.9580755	0.8007651	0.1667400	0.1706850	−9552.862	0.214871	−0.931142
10	−0.4853678	0.3904369	0.4044899	0.9582486	0.8004799	0.1667119	0.1707853	−9552.862	0.004754	−0.931872
11	−0.4853800	0.3904362	0.4044846	0.9582487	0.8004779	0.1667121	0.1707862	−9552.862	2.96 × ${10}^{-07}$	−0.931836
12	−0.4853800	0.3904362	0.4044846	0.9582487	0.8004779	0.1667121	0.1707862	−9552.862	1.18 × ${10}^{-10}$	−0.931836
13	−0.4853800	0.3904362	0.4044846	0.9582487	0.8004779	0.1667121	0.1707862	−9552.862	1.27 × ${10}^{-11}$	−0.931836

, respectively, in Appendix A.1, which describes the convergence process of the GTS parameter for Bitcoin and Ethereum data. The convergence process was obtained using the Newton–Raphson iteration algorithm (23). Each row has eleven columns made of the iteration number, the seven parameters $\mu$, ${\beta }_{+}$, ${\beta }_{-}$, ${\alpha }_{+}$, ${\alpha }_{-}$, ${\lambda }_{+}$, ${\lambda }_{-}$, and three statistical indicators, the log-likelihood ($Log(ML)$), the norm of the partial derivatives ($\left|\right|{\displaystyle \frac{dLog(ML)}{dV}}\left|\right|$), and the maximum value of the eigenvalues ($MaxEigenValue$). The statistical indicators aim at checking if the two necessary and sufficient conditions described in (22) are all met. $Log(ML)$ displays the value of the Naperian logarithm of the likelihood function $L(x,V)$, as described in (20); $\left|\right|{\displaystyle \frac{dLog(ML)}{dV}}\left|\right|$ displays the value of the norm of the first derivatives (${\displaystyle \frac{dl(x,V)}{d{V}_j}}$) described in (21); and $MaxEigenValue$ displays the maximum value of the seven eigenvalues generated by the Hessian matrix (${\displaystyle \frac{d^{2}l(x,V)}{d{V}_{k}d{V}_j}}$), as described in (21).

Similarly,

Table A7. Convergence of the Kobol parameter for Bitcoin return data.

$\mathit{Iterations}$	$\mathit{\mu }$	$\mathit{\beta }$	${\mathit{\alpha }}_{+}$	${\mathit{\alpha }}_{-}$	${\mathit{\lambda }}_{+}$	${\mathit{\lambda }}_{-}$	$\mathit{Log}(\mathit{ML})$	$\left|\right|\frac{\mathit{dLog}(\mathit{ML})}{\mathit{dV}}\left|\right|$	$\mathit{MaxEigenValue}$
1	−0.1215714	0.3155483	0.7477142	0.5445652	0.2465296	0.17477186	−10614.93879	450.0556974	25.6678081
2	−0.255172	0.36516958	0.73215119	0.53194253	0.22955186	0.17909072	−10607.01058	51.74982347	−46.893383
3	−0.2912276	0.37096854	0.75410108	0.53529439	0.23387716	0.18070591	−10606.81236	1.484798964	−53.728563
4	−0.2922408	0.36574333	0.7559582	0.53508403	0.23560819	0.18189913	−10606.81041	0.258928464	−53.391237
5	−0.2928641	0.36714311	0.75591147	0.53512239	0.23524801	0.18158588	−10606.81025	0.01286122	−53.406734
6	−0.292837	0.36708319	0.75591382	0.53512107	0.23526357	0.18159941	−10606.81025	0.00174219	−53.40643
7	−0.2928328	0.36707373	0.75591419	0.53512086	0.23526603	0.18160154	−10606.81025	1.18 × ${10}^{-05}$	−53.406384
8	−0.2928328	0.36707379	0.75591419	0.53512086	0.23526602	0.18160153	−10606.81025	1.60 × ${10}^{-06}$	−53.406384
9	−0.2928328	0.36707379	0.75591419	0.53512086	0.23526601	0.18160153	−10606.81025	2.18 × ${10}^{-07}$	−53.406384
10	−0.2928328	0.36707379	0.75591419	0.53512086	0.23526601	0.18160153	−10606.81025	1.09 × ${10}^{-08}$	−53.406384

and

Table A9. Convergence of the CGMY parameter for Ethereum return data.

$\mathit{Iterations}$	$\mathit{\mu }$	$\mathit{\beta }$	$\mathit{\alpha }$	${\mathit{\lambda }}_{+}$	${\mathit{\lambda }}_{-}$	$\mathit{Log}(\mathit{ML})$	$\left|\right|\frac{\mathit{dLog}(\mathit{ML})}{\mathit{dV}}\left|\right|$	$\mathit{MaxEigenValue}$
1	−0.48538	0.39043616	0.95824875	0.16671208	0.17078617	−9596.2658	1653.57149	−140.21456
2	−0.0545131	0.40148247	0.88205674	0.15875317	0.18060554	−9554.6834	80.0921993	−19.637869
3	−0.1479632	0.39049434	0.88271998	0.15631408	0.18704084	−9553.9065	3.84112398	−44.76325
4	−0.1465893	0.40345482	0.88868927	0.15450383	0.18507239	−9553.9036	0.4436942	−55.029934
5	−0.1472464	0.39683622	0.88667597	0.15564017	0.18628001	−9553.9026	0.14094418	−51.274906
6	−0.1470247	0.39907668	0.88736036	0.15525581	0.18587143	−9553.9025	0.05563606	−52.523819
7	−0.1471148	0.39816841	0.88708569	0.15541227	0.18603772	−9553.9025	0.02143506	−52.017698
8	−0.1470898	0.39842098	0.88716234	0.15536883	0.18599155	−9553.9025	0.00019543	−52.158334
9	−0.14709	0.39841855	0.88716161	0.15536924	0.18599199	−9553.9025	1.16 × ${10}^{-05}$	−52.156981
10	−0.14709	0.39841867	0.88716164	0.15536922	0.18599197	−9553.9025	1.78 × ${10}^{-06}$	−52.157046
11	−0.14709	0.39841869	0.88716165	0.15536922	0.18599197	−9553.9025	2.71 × ${10}^{-07}$	−52.157055
12	−0.14709	0.39841869	0.88716165	0.15536922	0.18599197	−9553.9025	4.14 × ${10}^{-08}$	−52.157057

describe the convergence process of the Kobol distribution parameter for Bitcoin returns and the CGMY distribution parameter for Ethereum returns.

GTS parameter estimations in Table 1 and Table 2 are used to evaluate the impact of each parameter on the GTS probability density function. As shown in Figure 3 and Figure 4, the effect of the GTS parameters on the probability density function has the same patterns on Bitcoin and Ethereum returns. However, the magnitudes are different. As shown in Figure 3a,b, ${\beta }_{-}$ (${\alpha }_{-}$) has a higher effect on the probability density function (pdf) than ${\beta }_{+}$ (${\alpha }_{+}$). However, ${\lambda }_{-}$ and ${\lambda }_{+}$ in both graphs seem symmetric and have the same impact.

4.3. Evaluation of the Method of Moments

The method of moments estimates the parameters of the GTS distribution by equating empirical moments and the theoretical moments of the GTS distribution. We empirically estimate the kth moments ($m_k=E(x^k)$), based on sample data $x={\left(x_j\right)}_{1\le j\le m}$ as follows:

$$\begin{array}{c}\hfill {\widehat{m}}_k={\displaystyle \frac{1}{M}}\sum _{j=1}^m{x}_j^k\mathrm{for}k=1,\dots ,7.\end{array}$$

(40)

On the other side, the cumulants (${\kappa }_k$) in Theorem 4 can be related to the moment of the GTS distribution by the following relationship (Poloskov 2021; Rota and Shen 2000; Smith 1995):

$$\begin{array}{c}\hfill m_k=E(x^k)=\sum _{j=1}^{k-1}\left({\displaystyle \genfrac{}{}{0pt}{}{k-1}{j-1}}\right){\kappa }_j{m}_{k-j}+{\kappa }_k\mathrm{for}k=1,\dots ,7.\end{array}$$

(41)

The method of moments estimator for $V=(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},{\lambda }_{-})$ is defined as the solution to the following system of equations:

$$\begin{array}{c}\hfill {\widehat{m}}_k=m_k\mathrm{for}k=1,\dots ,7.\end{array}$$

(42)

The system of Equation (42) is often not analytically solvable. For the conditions of existence and uniqueness of the solution, refer to (Küchler and Tappe 2013).

The maximum-likelihood GTS parameter estimations in Table 1 and Table 2 are used to evaluate the system of equations in (42). As shown in

Table 3. Evaluation of the method of moments.

	Bitcoin BTC			Ethereum
	Empirical(1)	Theoretical(2)	${\displaystyle \frac{(\mathbf{1})-(\mathbf{2})}{\mathbf{2}}}$	Empirical(1)	Theoretical(2)	${\displaystyle \frac{(\mathbf{1})-(\mathbf{2})}{\mathbf{2}}}$
Sample size	4083			3246
${\widehat{m}}_1$	0.152	0.152	0.0%	0.267	0.267	0.0%
${\widehat{m}}_2$	14.960	15.020	0.4%	27.161	27.388	0.8%
${\widehat{m}}_3$	−11.320	−15.640	27.6%	55.363	57.867	4.3%
${\widehat{m}}_4$	2033	2256	9.8%	5267	6307	16.5%
${\widehat{m}}_5$	−5823	−15,480	62.3%	22,368	32518	31.2%
${\widehat{m}}_6$	670,695	1,123,215	40.2%	2,114,788	4,361,562	51.5%
${\widehat{m}}_7$	−1,997,196	−19,777,988	89.9%	12,411,809	39,253,001	68.3%
Standard deviation 1	3.865	3.873	0.2%	5.206	5.226	0.4%
Skewness 2	−0.314	−0.387	18.8%	0.238	0.252	5.2%
Kurtosis 3	9.154	10.082	9.2%	7.112	8.385	15.2%
Max value	28.052			29.013
Min value	−26.620			−29.174

¹ $\sigma =\sqrt{{\kappa }_2}$; ² Skewness is estimated as ${\displaystyle \frac{{\kappa }_3}{{\kappa }_2^{3/2}}}$; ³ Kurtosis is estimated as $3+{\displaystyle \frac{{\kappa }_4}{{\kappa }_2^2}}$; ${\kappa }_1$, ${\kappa }_2$ and ${\kappa }_2$ are defined in (13).

, the solution of the maximum-likelihood method satisfies to a certain extent the equations for the first four moments: ${\widehat{m}}_1$, ${\widehat{m}}_2$, ${\widehat{m}}_3$, ${\widehat{m}}_4$ in the system (42). The seventh-moment equation has the highest relative error: 89.9% for Bitcoin (BTC) and 68.3% for Ethereum. Therefore, the maximum likelihood GTS parameter estimation is not the same as the GTS parameter estimation from the method of moments.

In addition to the method of moments estimations, the lower relative errors in Table 3 show that empirical and theoretical standard deviation ($\sigma$), skewness, and kurtosis seem to be consistent for Bitcoin and Ethereum. The empirical and theoretical statistics show that the average Ethereum daily return is greater and more volatile than the Bitcoin daily returns. Both assets are thicker than the normal distribution. However, the daily return of Bitcoin is skewed to the left, whereas the daily return of Ethereum is skewed to the right.

5. Fitting Tempered Stable Distribution to Traditional Indices: S&P 500 and SPY EFT

5.1. Data Summaries

The Standard & Poor’s 500 Composite Stock Price Index, also known as the S&P 500, is a stock index that tracks the share prices of 500 of the largest public companies with stocks listed on the New York Stock Exchange (NYSE) and the Nasdaq in the United States. It was introduced in 1957 and is often treated as a proxy for describing the overall health of the stock market or the United States (US) economy. The SPDR S&P 500 ETF (SPY), also known as the SPY ETF, is an Exchange-Traded Fund (ETF) that tracks the performance of the S&P 500. SPY ETF provides a mutual fund’s diversification, the stock’s flexibility, and lower trading fees. The data were extracted from Yahoo Finance. The historical prices span from 04 January 2010 to 22 July 2024 and were adjusted for splits and dividends.

The daily price dynamics are provided in Figure 5. Prices have an increasing trend, even after being temporally disrupted in the first quarter of 2020 by the coronavirus pandemic. The S&P 500 is priced in thousands of US dollars, whereas the SPY ETF is in hundreds of US dollars. The SPY ETF is cheaper and provides all the attributes of the S&P 500 index, as shown in Figure 5a,b.

Let the number of observations be m and the daily observed price be $S_j$ on day $t_j$ with $j=1,\dots ,m$; $t_1$ is the first observation date (4 January 2010) and $t_m$ is the last observation date (22 July 2024). The daily return, $y_j$, is computed as in (43):

$$\begin{array}{c}\hfill y_j=log(S_j/S_{j-1})j=2,\dots ,m.\end{array}$$

(43)

The SPY ETF aims to mirror the performance of the S&P 500. Figure 6a,b look similar, which is consistent with the goal of the SPY ETF. As shown in Figure 6a,b, the daily return reaches the lowest level ($-12.7\%$ for the S&P 500 and $-11.5\%$ for the SPY ETF) in the first quarter of 2020 amid the coronavirus pandemic and massive disruptions in the global economy. Nine values were identified as outliers and removed from the dataset to avoid a negative impact on GTS model estimation and the empirical statistics.

5.2. Multidimensional Estimation Results for Traditional Indices

The estimation results are provided in

Table 4. Maximum-likelihood GTS parameter estimation for S&P 500 index.

Model	Param	Estimate	Std Err	z	$\mathit{Pr}(\mathit{Z}>|\mathit{z}|)$	[95% Conf.Interval]
GTS	$\mu$	−0.249408	(0.208)	−1.20	2.3 × ${10}^{-01}$	−0.658	0.159
	${\beta }_{+}$	0.328624	(0.308)	1.07	2.9 × ${10}^{-01}$	−0.275	0.932
	${\beta }_{-}$	0.088640	(0.176)	0.50	6.1 × ${10}^{-01}$	−0.256	0.433
	${\alpha }_{+}$	0.792426	(0.350)	2.26	2.4 × ${10}^{-02}$	0.106	1.479
	${\alpha }_{-}$	0.542250	(0.107)	5.09	3.6 × ${10}^{-07}$	0.333	0.751
	${\lambda }_{+}$	1.279743	(0.348)	3.68	2.4 × ${10}^{-04}$	0.597	1.962
	${\lambda }_{-}$	0.937133	(0.144)	6.50	8.0 × ${10}^{-11}$	0.655	1.220
	Log(ML)	−4920
	AIC	9851
	BIK	9898
GBM	$\mu$	0.044875	(0.018)	2.51	1.2 × ${10}^{-02}$	0.010	0.080
	$\sigma$	1.081676	(0.027)	39.53	0.000	1.028	1.135
	Log(ML)	−5330
	AIC	10,665
	BIK	10,677

for S&P 500 return data and

Table 5. Maximum-likelihood GTS parameter estimation for SPY EFT data.

Model	Param	Estimate	Std Err	z	$\mathit{Pr}(\mathit{Z}>|\mathit{z}|)$	[95% Conf.Interval]
GTS	$\mu$	−0.260643	(0.135)	−1.94	5.3 × ${10}^{-02}$	−0.524	0.003
	${\beta }_{+}$	0.340880	(0.189)	1.80	7.1 × ${10}^{-02}$	−0.030	0.711
	${\beta }_{-}$	0.022212	(0.212)	0.10	9.2 × ${10}^{-01}$	−0.393	0.437
	${\alpha }_{+}$	0.787757	(0.225)	3.50	4.6 × ${10}^{-04}$	0.347	1.229
	${\alpha }_{-}$	0.597110	(0.141)	4.22	2.4 × ${10}^{-05}$	0.320	0.874
	${\lambda }_{+}$	1.288555	(0.226)	5.70	1.2 × ${10}^{-08}$	0.846	1.731
	${\lambda }_{-}$	1.014353	(0.177)	5.74	9.4 × ${10}^{-09}$	0.668	1.361
	Log(ML)	−4893
	AIC	9800
	BIK	9843
GBM	$\mu$	0.054344	(0.017)	3.13	1.8 × ${10}^{-03}$	0.020	0.088
	$\sigma$	1.050217	(0.026)	40.71	0.000	1.000	1.101
	Log(ML)	−54,275
	AIC	10,554
	BIK	10,566

for SPY EFT return data. As previously, the log-likelihood, AIC, and BIK statistics suggest that the GTS distribution with seven parameters performs better than the two-parameter normal distribution (GBM).

As shown in both Table 4 and Table 5, the ML estimate of $\mu$ is negative, while the others are positive, as expected in the literature. The asymptotic standard error for $\mu$, ${\beta }_{+}$ and ${\beta }_{-}$ are pretty large and result in $\mu$, ${\beta }_{+}$ and ${\beta }_{-}$ not being significantly different from zero.

However, other parameters have larger t-statistics ($\left|z\right|>2$) and are statistically significant at 5%. Except for the index of stability parameters (${\beta }_{+}$, ${\beta }_{-}$), the estimation results for the S&P 500 and SPY ETF indexes show that the difference in skewness parameters (${\lambda }_{+}$, ${\lambda }_{-}$) and intensity parameters (${\alpha }_{+}$, ${\alpha }_{-}$) are positive but are not statistically significant.

The hypothesis with ${\beta }_{+}={\beta }_{-}=0$ was considered by fitting the S&P 500 and SPY ETF indexes to the bilateral Gamma distribution. The estimation results are summarized in Appendices Appendix C.1 and Appendix C.2. As shown in

Table A10. BG maximum-likelihood estimation for S&P 500 return data.

Model	Parameter	Estimate	Std Err	z	$\mathit{Pr}(\mathit{Z}>|\mathit{z}|)$	[95% Conf.Interval]
GTS	$\mu$	−0.031467	(0.010)	−3.07	2.1 × ${10}^{-03}$	−0.052	−0.011
	${\alpha }_{+}$	1.092741	(0.058)	18.98	2.6 × ${10}^{-80}$	0.980	1.206
	${\alpha }_{-}$	0.701784	(0.042)	16.80	2.3 × ${10}^{-63}$	0.620	0.784
	${\lambda }_{+}$	1.539690	(0.064)	22.82	3.1 × ${10}^{-115}$	1.407	1.672
	${\lambda }_{-}$	1.110737	(0.050)	22.07	6.6 × ${10}^{-108}$	1.012	1.209
	Log(ML)	−4925
	AIC	9859
	BIK	9890

and

Table A12. BG maximum-likelihood estimation for SPY EFT return data.

Model	Parameter	Estimate	Std Err	z	$\mathit{Pr}(\mathit{Z}>|\mathit{z}|)$	[95% Conf.Interval]
GTS	$\mu$	0.015048	(0.012)	1.28	2.0 × ${10}^{-01}$	−0.008	0.038
	${\alpha }_{+}$	1.068239	(0.067)	16.02	8.6 × ${10}^{-58}$	0.938	1.199
	${\alpha }_{-}$	0.764449	(0.044)	17.33	3.0 × ${10}^{-67}$	0.678	0.851
	${\lambda }_{+}$	1.525718	(0.073)	20.98	1.1 × ${10}^{-97}$	1.383	1.668
	${\lambda }_{-}$	1.156439	(0.052)	22.15	1.1 × ${10}^{-108}$	1.054	1.259
	Log(ML)	−4899
	AIC	9807
	BIK	9838

, the skewness parameters (${\lambda }_{+}$, ${\lambda }_{-}$) are positive and statistically significant, and the difference (${\lambda }_{+}-{\lambda }_{-}$) is also positive and statistically significant, which proves that the S&P 500 and SPY ETF returns are skewed to the left. We have the same statistical features for the intensity parameters (${\alpha }_{+}$, ${\alpha }_{-}$), and both indexes are more likely to produce positive returns than negative returns. Refer to (Küchler and Tappe 2008; Nzokem 2021a) for more details on the bilateral Gamma distribution.

The likelihood ratio test in

Table 6. Likelihood ratio test statistic and p-value.

		GTS	GTS Variants	${\mathit{\chi }}^2$-Value	df	p-Value
	Log(ML)	−10,606.73	−10,606.81	0.1525	1	0.6962
Bitcoin	AIC	21,227.47	21,225.62
	BIK	21,271.67	21,263.51
	Log(ML)	−9552.86	−9553.90	2.0810	2	0.3533
Ethereum	AIC	19,119.72	19,117.81
	BIK	19,162.32	19,148.23
	Log(ML)	−4920.52	−4924.62	8.1828	2	0.0167
S&P 500	AIC	9851.06	9859.24
	BIK	9898.49	9890.26
	Log(ML)	−4893.21	−4898.67	10.9234	2	0.0042
SPY ETF	AIC	9800.42	9807.34
	BIK	9843.84	9838.36

shows that, even with non-statistically significant parameters, the GTS distribution fits significantly better than the bilateral Gamma distribution for both the S&P 500 and SPY ETF indexes. Contrary to the AIC statistics, the BIK statistics do not provide the same information. A comprehensive and detailed examination of the statistical significance of the results is carried out in Section 6.

Table 4 and Table 5 summarize the last row of

Table A3. Convergence of the GTS parameter for S&P 500 return data.

$\mathit{Iterations}$	$\mathit{\mu }$	${\mathit{\beta }}_{+}$	${\mathit{\beta }}_{-}$	${\mathit{\alpha }}_{+}$	${\mathit{\alpha }}_{-}$	${\mathit{\lambda }}_{+}$	${\mathit{\lambda }}_{-}$	$\mathit{Log}(\mathit{ML})$	$\left|\right|\frac{\mathit{dLog}(\mathit{ML})}{\mathit{dV}}\left|\right|$	$\mathit{MaxEigenValue}$
1	−0.2606426	0.34087979	0.02221141	0.78775729	0.59711061	1.28855513	1.01435308	−4921.0858	147.214541	−0.476265
2	−0.2747887	0.37848567	0.02517846	0.72538248	0.594628	1.22107935	1.01081205	−4920.9765	107.910271	−12.169518
3	−0.2852743	0.34562742	0.01628972	0.78353361	0.58024658	1.27423544	0.9888729	−4920.6236	23.70873	11.9588258
4	−0.2971254	0.37985815	0.05392593	0.74068472	0.55972179	1.22737986	0.96278568	−4920.5493	4.21443356	0.29705471
5	−0.3415082	0.42600675	0.0432239	0.69783497	0.56106365	1.18286494	0.966753	−4920.5722	37.0642417	−1.7903876
6	−0.2995817	0.40315129	0.12236507	0.7168274	0.522172	1.20383351	0.9117451	−4920.574	3.07232514	−0.7101089
7	−0.2944623	0.3977257	0.12218751	0.72174351	0.52260032	1.20899714	0.9121201	−4920.5701	2.63567879	−1.0187469
8	−0.2767429	0.37561063	0.11561097	0.7427615	0.52696799	1.23067384	0.91742165	−4920.5511	1.83311761	−2.1103436
9	−0.274204	0.37177939	0.11355883	0.74659763	0.52814335	1.2345524	0.91893546	−4920.5477	1.75839181	−2.177405
10	−0.2559812	0.34147926	0.09643312	0.77815581	0.53784221	1.26594249	0.93144448	−4920.5308	1.33811298	−2.6954121
11	−0.2496977	0.32954013	0.08928069	0.79125494	0.54186044	1.27868642	0.93662846	−4920.5291	0.79520373	−2.8166517
12	−0.2494237	0.32866495	0.08869445	0.79238161	0.54221561	1.27970094	0.93708759	−4920.5291	0.00166731	−2.6765739
13	−0.2494072	0.32862462	0.08864569	0.79242579	0.54224632	1.27974278	0.93712865	−4920.5291	0.00013552	−2.6768326
14	−0.2494082	0.32862428	0.08864047	0.79242619	0.54224944	1.27974312	0.93713293	−4920.5291	1.47 × ${10}^{-05}$	−2.6766945
15	−0.2494083	0.32862424	0.08863992	0.79242624	0.54224977	1.27974315	0.93713338	−4920.5291	1.57 × ${10}^{-06}$	−2.67668
16	−0.2494083	0.32862424	0.08863985	0.79242624	0.54224981	1.27974316	0.93713344	−4920.5291	1.89 × ${10}^{-09}$	−2.6766783
17	−0.2494083	0.32862424	0.08863985	0.79242624	0.54224981	1.27974316	0.93713344	−4920.5291	2.09 × ${10}^{-10}$	−2.6766783

and

Table A4. Convergence of the GTS parameter for SPY EFT return data.

$\mathit{Iterations}$	$\mathit{\mu }$	${\mathit{\beta }}_{+}$	${\mathit{\beta }}_{-}$	${\mathit{\alpha }}_{+}$	${\mathit{\alpha }}_{-}$	${\mathit{\lambda }}_{+}$	${\mathit{\lambda }}_{-}$	$\mathit{Log}(\mathit{ML})$	$\left|\right|\frac{\mathit{dLog}(\mathit{ML})}{\mathit{dV}}\left|\right|$	$\mathit{MaxEigenValue}$
1	−0.0518661	0.1161846	0.2186548	1.04269292	0.52712574	1.52244991	0.91168779	−4894.2279	14.7801725	−6.5141947
2	−0.1102477	0.18491276	0.17478472	0.94756655	0.52844271	1.4399315	0.91415148	−4893.8278	29.8166141	−1.9290981
3	−0.2094204	0.29377592	0.0891446	0.84029122	0.56054563	1.34797271	0.96500981	−4893.3554	16.9940095	4.33892902
4	−0.1985564	0.29758208	0.13230013	0.83156167	0.53656079	1.33833078	0.93230856	−4893.4206	10.9048744	1.04588745
5	−0.078883	0.25939922	0.39611543	0.84865673	0.40365522	1.35595932	0.7410936	−4895.8806	241.028178	94.6293224
6	−0.0753571	0.26704857	0.33754158	0.84120908	0.45446164	1.3452823	0.80751063	−4894.3899	25.1995505	−2.805571
7	−0.196642	0.31624372	0.20068543	0.80509106	0.50322368	1.30837612	0.88967028	−4893.888	140.257551	34.770691
8	−0.1898283	0.3045047	0.15900291	0.81380451	0.52672075	1.31341912	0.91775259	−4893.4694	6.29433991	−4.6080872
9	−0.2275214	0.32940996	0.10770535	0.79340215	0.55020025	1.29360449	0.95260474	−4893.3049	7.34361008	−8.1891832
10	−0.2726283	0.34972465	0.01601222	0.78061523	0.59736004	1.28153304	1.01757433	−4893.2192	14.0784211	−3.6408772
11	−0.2499816	0.32645286	0.01851524	0.80217301	0.60018154	1.30243672	1.01792703	−4893.2125	6.27794755	−4.6455546
12	−0.2575953	0.33832596	0.02643321	0.79008383	0.59450637	1.29085215	1.01101001	−4893.208	1.23298227	−6.8035318
13	−0.2607071	0.34052644	0.02075252	0.78817075	0.59805376	1.28895161	1.01555438	−4893.2076	0.07363298	−6.71708
14	−0.2606368	0.34088815	0.02227012	0.78774693	0.59707082	1.28854532	1.01430383	−4893.2076	0.00156771	−6.6908109
15	−0.2606432	0.34087911	0.02220633	0.78775813	0.59711397	1.28855593	1.01435731	−4893.2076	0.00010164	−6.6915902
16	−0.2606426	0.34087985	0.02221188	0.78775721	0.5971103	1.28855506	1.01435268	−4893.2076	8.45 × ${10}^{-06}$	−6.6915177
17	−0.2606426	0.34087979	0.02221141	0.78775729	0.59711061	1.28855513	1.01435308	−4893.2076	7.21 × ${10}^{-07}$	−6.6915239

, respectively, in Appendix A.1, which describes the convergence process of the GTS parameter for S&P 500 index, and SPY ETF return data. The convergence process was obtained using the Newton–Raphson iteration algorithm (23). Each row has eleven columns made of the iteration number; the seven parameters $\mu$, ${\beta }_{+}$, ${\beta }_{-}$, ${\alpha }_{+}$, ${\alpha }_{-}$, ${\lambda }_{+}$, ${\lambda }_{-}$; and three statistical indicators, the log-likelihood ($Log(ML)$), the norm of the partial derivatives ($\left|\right|{\displaystyle \frac{dLog(ML)}{dV}}\left|\right|$), and the maximum value of the eigenvalues ($MaxEigenValue$). The statistical indicators aim at checking if the two necessary and sufficient conditions described in (22) are all met. $Log(ML)$ displays the value of the Naperian logarithm of the likelihood function $L(x,V)$, as described in (20); $\left|\right|{\displaystyle \frac{dLog(ML)}{dV}}\left|\right|$ displays the value of the norm of the first derivatives (${\displaystyle \frac{dl(x,V)}{d{V}_j}}$) described in Equation (21); and $MaxEigenValue$ displays the maximum value of the seven eigenvalues generated by the Hessian matrix (${\displaystyle \frac{d^{2}l(x,V)}{d{V}_{k}d{V}_j}}$), as described in (21).

Similarly,

Table A11. Convergence of the BG parameter for S&P 500 return data.

$\mathit{Iterations}$	$\mathit{\mu }$	${\mathit{\alpha }}_{+}$	${\mathit{\alpha }}_{-}$	${\mathit{\lambda }}_{+}$	${\mathit{\lambda }}_{-}$	$\mathit{Log}(\mathit{ML})$	$\left|\right|\frac{\mathit{dLog}(\mathit{ML})}{\mathit{dV}}\left|\right|$	$\mathit{MaxEigenValue}$
1	0	0.79242624	0.54224981	1.27974316	0.93713344	−4951.1439	1138.53458	−265.251
2	−0.0038447	0.93153413	0.64461254	1.41132138	1.05504868	−4931.7583	549.025405	−171.22804
3	−0.0103214	1.03062846	0.70198868	1.49426667	1.10555794	−4926.8412	286.215345	−126.18156
4	−0.0186317	1.07922912	0.71421996	1.53377391	1.11392285	−4925.393	135.694287	−113.12071
5	−0.0279475	1.09450205	0.70493092	1.54418172	1.10795103	−4924.7065	38.0551545	−116.58686
6	−0.0313951	1.09325663	0.70162581	1.54038766	1.10996346	−4924.621	1.54417452	−120.06271
7	−0.0314671	1.09274119	0.70178365	1.53969027	1.11073682	−4924.6205	0.02788236	−120.35435
8	−0.0314664	1.09276971	0.70182788	1.53971127	1.11079928	−4924.6205	0.00198685	−120.34482
9	−0.0314663	1.09277213	0.70183158	1.5397131	1.11080431	−4924.6205	0.00016039	−120.34394
10	−0.0314662	1.09277232	0.70183188	1.53971325	1.11080472	−4924.6205	1.29 × ${10}^{-05}$	−120.34387
11	−0.0314662	1.09277234	0.7018319	1.53971326	1.11080475	−4924.6205	1.04 × ${10}^{-06}$	−120.34387
12	−0.0314662	1.09277234	0.7018319	1.53971326	1.11080476	−4924.6205	8.43 × ${10}^{-08}$	−120.34387
13	−0.0314662	1.09277234	0.7018319	1.53971326	1.11080476	−4924.6205	6.80 × ${10}^{-09}$	−120.34387
14	−0.0314662	1.09277234	0.7018319	1.53971326	1.11080476	−4924.6205	5.63 × ${10}^{-10}$	−120.34387
15	−0.0314662	1.09277234	0.7018319	1.53971326	1.11080476	−4924.6205	5.73 × ${10}^{-11}$	−120.34387

and

Table A13. Convergence of the BG parameter for SPY EFT return data.

$\mathit{Iterations}$	$\mathit{\mu }$	${\mathit{\alpha }}_{+}$	${\mathit{\alpha }}_{-}$	${\mathit{\lambda }}_{+}$	${\mathit{\lambda }}_{-}$	$\mathit{Log}(\mathit{ML})$	$\left|\right|\frac{\mathit{dLog}(\mathit{ML})}{\mathit{dV}}\left|\right|$	$\mathit{MaxEigenValue}$
1	0	0.78775729	0.59711061	1.28855513	1.01435308	−4918.7331	406.35365	−252.28104
2	0.02867773	0.97562263	0.67572827	1.46822249	0.97596762	−4908.5992	226.190986	−116.11753
3	0.02727407	1.05127517	0.78618306	1.51846501	1.17883595	−4899.0798	45.2281041	−96.788275
4	0.00834089	1.07692251	0.75226045	1.5348003	1.14577232	−4898.9955	131.637516	−107.77617
5	0.01126962	1.07242358	0.7568497	1.53011913	1.1494212	−4898.751	48.0005418	−103.03258
6	0.01386478	1.06933921	0.76167668	1.52688303	1.15363987	−4898.6763	11.3246873	−100.1483
7	0.01492745	1.06823047	0.76397541	1.52573245	1.15588409	−4898.6693	0.98802026	−99.171136
8	0.01504464	1.06821119	0.76439389	1.52569528	1.15636567	−4898.6693	0.02529683	−99.040575
9	0.01504742	1.06823539	0.76444163	1.5257152	1.15642976	−4898.6693	0.00300624	−99.030178
10	0.01504762	1.06823881	0.76444764	1.52571803	1.15643796	−4898.6693	0.00038489	−99.028926
11	0.01504764	1.06823925	0.76444841	1.52571839	1.15643901	−4898.6693	4.92 × ${10}^{-05}$	−99.028765
12	0.01504765	1.06823931	0.76444851	1.52571844	1.15643915	−4898.6693	6.30 × ${10}^{-06}$	−99.028745
13	0.01504765	1.06823932	0.76444852	1.52571845	1.15643917	−4898.6693	1.32 × ${10}^{-08}$	−99.028742
14	0.01504765	1.06823932	0.76444852	1.52571845	1.15643917	−4898.6693	1.69 × ${10}^{-09}$	−99.028742
15	0.01504765	1.06823932	0.76444852	1.52571845	1.15643917	−4898.6693	2.23 × ${10}^{-10}$	−99.028742

describe the convergence process of the bilateral Gamma distribution parameter for S&P 500 index and SPY ETF return data.

The GTS parameter estimations in Table 4 and Table 5 were used to evaluate the impact of the parameters on the GTS probability density function. As shown in Figure 7 and Figure 8, the effect of the GTS parameters on the probability density function generated by the S&P 500 and SPY ETF have the same patterns. However, the magnitudes are different. As shown in Figure 7a,b on the S&P 500 return data, ${\beta }_{+}$ (${\alpha }_{+}$) has a higher effect on the probability density function than ${\beta }_{-}$ (${\alpha }_{-}$). However, ${\lambda }_{-}$ and ${\lambda }_{+}$ in Figure 7c are symmetric and have the same impact.

5.3. Evaluation of the Methods of Moments

The maximum-likelihood GTS parameter estimations in Table 4 and Table 5 are used to evaluate the system of equations in (42). As shown in

Table 7. Evaluation of the methods of moments.

	S&P 500 Index			SPY ETF
	Empirical(1)	Theoretical(2)	${\displaystyle \frac{(\mathbf{1})-(\mathbf{2})}{\mathbf{2}}}$	Empirical(1)	Theoretical(2)	${\displaystyle \frac{(\mathbf{1})-(\mathbf{2})}{\mathbf{2}}}$
Sample size	3656			3655
${\widehat{m}}_1$	0.045	0.045	−0.5%	0.054	0.054	0.0%
${\widehat{m}}_2$	1.069	1.083	−1.3%	1.053	1.044	0.8%
${\widehat{m}}_3$	−0.447	−0.341	31.2%	−0.214	−0.351	−39.0%
${\widehat{m}}_4$	8.371	9.764	−14.3%	8.197	7.691	6.6%
${\widehat{m}}_5$	−16.386	−11.128	47.3%	−3.969	−12.717	−68.8%
${\widehat{m}}_6$	193.563	247.811	−21.9%	157.645	162.048	−2.7%
${\widehat{m}}_7$	−840.097	−547.882	53.3%	−85.003	−602.447	−85.9%
Standard deviation 1	1.082	1.033	4.7%	1.050	1.021	2.9%
Skewness 2	−0.432	−0.535	−19.2%	−0.358	−0.490	−26.9%
Kurtosis 3	8.413	7.435	13.1%	7.495	7.177	4.4%
Max value	6.797			6.501
Min value	−7.901			−6.734

¹ $\sigma =\sqrt{{\kappa }_2}$; ² Skewness is estimated as ${\displaystyle \frac{{\kappa }_3}{{\kappa }_2^{3/2}}}$; ³ Kurtosis is estimated as $3+{\displaystyle \frac{{\kappa }_4}{{\kappa }_2^2}}$; ${\kappa }_1$, ${\kappa }_2$ and ${\kappa }_2$ are defined in (13).

, the solution of the maximum-likelihood method satisfies to a certain extent the equations (42) for the following first four moments: ${\widehat{m}}_1$, ${\widehat{m}}_2$, ${\widehat{m}}_4$, ${\widehat{m}}_5$. As for Bitcoin and Ethereum, the seventh-moment equation has the highest relative error: 53.3% for S&P 500 index and −85.9% for SPY ETF. Therefore, the maximum-likelihood GTS parameter estimation is not the GTS parameter estimation from the method of moments.

In addition to the moment estimations in Table 7, the empirical and theoretical standard deviation ($\sigma$), skewness, and kurtosis are consistent with lower relative errors for both S&P 500 and SPY ETF. The empirical and theoretical statistics show that both assets are skewed to the left and also thicker than the normal distribution.

6. Goodness-of-Fit Analysis

6.1. Kolmogorov–Smirnov (KS) Analysis

Given the sample of daily return $\{y_1,y_2\dots y_m\}$ of size m and the empirical cumulative distribution function, $F_m(x)$, for each index, the Kolmogorov–Smirnov (KS) test is performed under the null hypothesis, $H_0$, that the sample $\{y_1,y_2\dots y_m\}$ comes from the GTS distribution, $F(x)$. The cumulative distribution function of the theoretical distribution, $F(x)$, needs to be computed. The density function, $f(x)$, does not have a closed form, the same for the cumulative function, $F(x)$, in (45). However, we know the closed form of the Fourier of the density function, $\mathcal{F}\left[f\right]$, and the relationship in (46) provides the Fourier of the cumulative distribution function, $\mathcal{F}\left[F\right]$. The GTS distribution function, $F(x)$, was computed from the inverse of the Fourier of the cumulative distribution, $\mathcal{F}\left[F\right]$, in (47):

$$\begin{array}{cc}\hfill Y& \sim GTS(\mu ,{\beta }_{+},{\beta }_{-},{\alpha }_{+},{\alpha }_{-},{\lambda }_{+},f{\lambda }_{-})\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill F(x)& ={\int }_{-\infty }^{x}f(t)\mathrm{d}tf\mathrm{is}\mathrm{the}\mathrm{density}\mathrm{function}\mathrm{of}Y\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill \mathcal{F}\left[F\right](x)& ={\displaystyle \frac{\mathcal{F}\left[f\right](x)}{ix}}+\pi \mathcal{F}\left[f\right](0)\delta (x)\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill F(x)& ={\displaystyle \frac{1}{2\pi }}{\int }_{-\infty }^{+\infty }{\displaystyle \frac{\mathcal{F}\left[f\right](y)}{iy}}{e}^{ixy}\mathrm{d}y+{\displaystyle \frac{1}{2}}\hfill \end{array}$$

(47)

See Appendix A in (Nzokem 2021a) for (46) proof.

The two-sided KS goodness-of-fit statistic ($D_m$) is defined as follows:

$$\begin{array}{c}\hfill D_m=\underset{x}{sup}|F(x)-F_m(x)|,\end{array}$$

(48)

where m is the sample size, $F_m(x)$ denotes the empirical cumulative distribution of $\{y_1,y_2\dots y_m\}$.

The distribution of Kolmogorov’s goodness-of-fit measure $D_m$ has been studied extensively in the literature. It was shown (Massey 1951) that the $D_m$ distribution is independent of the theoretical distribution, $F(x)$, under the null hypothesis, $H_0$. The discrete, mixed, and discontinuous distributions case has also been studied (Dimitrova et al. 2020). Under the null hypothesis, $H_0$, that the sample $\{y_1,y_2\dots y_m\}$ of size m comes from the hypothesized continuous distribution, it was shown (An 1933) that the asymptotic statistic $\sqrt{n}{D}_n$ converges to the Kolmogorov distribution.

The limiting form for the distribution function of Kolmogorov’s goodness-of-fit measure $D_m$ is

$$\begin{array}{c}\hfill \underset{m\to +\infty }{lim}Pr(\sqrt{m}{D}_m\le x)=1-2\sum _{k=1}^{+\infty }{(-1)}^{k-1}{e}^{-2k^2{x}^2}={\displaystyle \frac{\sqrt{2\pi }}{x}}\sum _{k=1}^{+\infty }{e}^{-{\displaystyle \frac{{(2k-1)}^2{\pi }^2}{8x^2}}}.\end{array}$$

(49)

The first representation was given in (An 1933), and the second came from a standard relation for theta functions (Marsaglia et al. 2003).

As shown in Figure 9, the asymptotic statistic, $\sqrt{n}{D}_n$, is a positively skewed distribution with a mean and a standard deviation (Marsaglia et al. 2003) as follows.

$$\begin{array}{c}\hfill \mu =\sqrt{{\displaystyle \frac{\pi }{2}}}log(2)\sim 0.8687,\sigma =\sqrt{{\displaystyle \frac{{\pi }^2}{12}}-{\mu }^2}\sim 0.2603.\end{array}$$

(50)

At a $5\%$ risk level, the risk threshold is $d=1.3581$ and represents the area in the shaded area under the probability density function.

The p-value of the test statistic, $D_m$, is computed based on (49) as follows:

$$\begin{array}{c}\hfill p\_value=Pr(D_m>{\widehat{D}}_m|H_0)=1-Pr(\sqrt{m}{D}_m\le \sqrt{m}{\widehat{D}}_m).\end{array}$$

(51)

A p-value is defined as the probability that values are even more extreme or more in the tail than our test statistic. A small p-value leads to a rejection of the null hypothesis, $H_0$, because the test statistic, $D_m$, is already extreme. We reject the hypothesis if the p-value is less than our level of significance, which we take to be equal to 0.05.

${\widehat{D}}_m$ is a realization value of the KS estimator $D_m$ computed from the sample $\{y_1,y_2\dots y_m\}$. ${\widehat{D}}_m$ is estimated (Krysicki et al. 1999) as follows:

$$\begin{array}{c}\hfill {\widehat{D}}_m=Max(\underset{0\le j\le P}{sup}|F(x_j)-F_m(x_j)|,\underset{1\le j\le P}{sup}|F(x_j)-F_m(x_{j-1})|).\end{array}$$

(52)

The following computations were performed for Bitcoin (BTC) data, and the quantity ${\widehat{D}}_m$ was obtained:

$$\begin{array}{cc}\hfill & \underset{0\le j\le P}{sup}|F(x_j)-F_m(x_j)|=0.01300\hfill \\ \hfill & \underset{1\le j\le P}{sup}|F(x_j)-F_m(x_{j-1})|=0.00538\hfill \\ \hfill & {\widehat{D}}_m=0.01300\hfill \\ \hfill & p\_value=prob(\sqrt{m}{D}_m>0.6903|H_0)=49.48\%.\hfill \end{array}$$

(53)

For each index, KS statistics (${\widehat{D}}_m$) and associated p-values were computed and summarized in

Table 8. Kolmogorov–Smirnov statistics and p-values.

	GTS			GBN			GTS Variants			Sample Size
Index	${\widehat{\mathit{D}}}_{\mathit{m}}$	$\sqrt{\mathit{m}}{\widehat{\mathit{D}}}_{\mathit{m}}$	p-Value	${\widehat{\mathit{D}}}_{\mathit{m}}$	$\sqrt{\mathit{m}}{\widehat{\mathit{D}}}_{\mathit{m}}$	p-Value	${\widehat{\mathit{D}}}_{\mathit{m}}$	$\sqrt{\mathit{m}}{\widehat{\mathit{D}}}_{\mathit{m}}$	p-Value	m
Bicoin BTC	0.013	0.830	0.494	0.106	6.803	0.000	0.014 1	0.863	0.445	4083
Ethereum	0.012	0.721	0.674	0.092	5.249	0.000	0.013 2	0.749	0.627	3246
S&P 500	0.012	0.750	0.627	0.091	5.550	0.000	0.014 3	0.897	0.395	3656
SPY ETF	0.014	0.869	0.436	0.089	5.438	0.000	0.016 3	1.010	0.258	3655

¹ Kobol distribution ($\beta ={\beta }_{-}={\beta }_{+}$); ² Carr–Geman–Madan–Yor (CGMY) distributions ($\beta ={\beta }_{-}={\beta }_{+}$; $\alpha ={\alpha }_{-}={\alpha }_{+}$); ³ bilateral Gamma distribution (${\beta }_{-}={\beta }_{+}=0$).

, along with the index sample size, m.

The asymptotic statistics, $\sqrt{n}{D}_n$, produced from the two-parameter geometric Brownian motion (GBM) hypothesis, have high values and show that the GBM hypothesis is always rejected. On the other hand, the high p-values generated by the asymptotic statistics suggest insufficient evidence to reject the assumption that the data were randomly sampled from a GTS. The same observations work for the GTS variants: the Kobol, CGMY, and bilateral Gamma distributions. In addition, as shown the p-value indicator in Table 8, the GTS distribution outperforms the bilateral Gamma distribution for the S&P 500 and SPY ETF indexes. However, the Kobol and CGMY distributions, respectively, for Bitcoin and Ethereum have almost the same performance as the GTS distribution.

6.2. Anderson–Darling Test Analysis

The Anderson–Darling test (Anderson 2008) is a goodness-of-fit test that allows the control of the hypothesis that the distribution of a random variable observed in a sample follows a certain theoretical distribution. The Anderson–Darling statistic belongs to the class of quadratic EDF statistics (Stephens 1974) based on the empirical distribution function. The quadratic EDF statistics measure the distance between the hypothesized distribution ($F(x)$) and empirical distribution. It is defined as

$$\begin{array}{c}\hfill m{\int }_{-\infty }^{+\infty }{\left(F_m(x)-F(x)\right)}^{2}w(x)dFx,\end{array}$$

(54)

where m is the number of elements in the sample, $w(x)$ is a weighting function, and $F_m(x)$ is the empirical distribution function defined on the sample of size m.

When the weighting function is $w(x)=1$, the statistic (54) is the Cramér–-Von Mises statistic, while the Anderson–Darling statistic is obtained by choosing the weighting function $w(x)=F(x)\left(1-F(x)\right)$. Compared with the Cramér–Von Mises statistic, the Anderson–Darling statistic places more weight on the tails of the distribution.

The Anderson–Darling statistic is

$$\begin{array}{c}\hfill A_m^2=m{\int }_{-\infty }^{+\infty }{\displaystyle \frac{\left(F_m(x)-F(x)\right)}{F(x)\left(1-F(x)\right)}}dF(x).\end{array}$$

(55)

It can be shown that the asymptotic distribution of the Anderson–Darling statistic, $A_m^2$, is independent of the theoretical distribution under the null hypothesis. The asymptotic distribution (Lewis 1961; Marsaglia and Marsaglia 2004) is defined as follows:

$$\begin{array}{cc}\hfill G(x)& =\underset{m\to \infty }{lim}Pr\left[A_m^2<x\right]=\sum _{j=0}^{+\infty }{a}_j{(x{b}_j)}^{-{\displaystyle \frac{1}{2}}}exp(-{\displaystyle \frac{b_j}{x}}){\int }_0^{+\infty }{f}_j(y)exp(-y^2)dy\hfill \\ \hfill f_j(y)& =exp\left({\displaystyle \frac{1}{8}}{\displaystyle \frac{x{b}_j}{y^{2}x+b_j}}\right),a_j={\displaystyle \frac{{(-1)}^j{(2)}^{\frac{1}{2}}(4j+1)\Gamma (j+\frac{1}{2})}{j!}}\hfill \\ \hfill b_j& ={\displaystyle \frac{1}{2}}{(4j+1)}^2{\pi }^2.\hfill \end{array}$$

(56)

As shown in Figure 10, the asymptotic distribution of the Anderson–Darling statistic ($A_m^2$) is a positively skewed distribution with a mean and a standard deviation (Anderson 2011) as follows

$$\begin{array}{c}\hfill \mu =1,\sigma =\sqrt{{\displaystyle \frac{2}{3}}({\pi }^2-9)}\sim 0.761.\end{array}$$

(57)

At a $5\%$ risk level, the risk threshold is $d=2.4941$ and represents the area in the shaded area under the probability density function.

The p-value of the test statistic, $A_m^2$, is defined as follows:

$$\begin{array}{c}\hfill p\text{-}value=prob(A_m^2>{\widehat{A}}_m^2|H_0)=1-G({\widehat{A}}_m^2).\end{array}$$

(58)

In order to compute the Anderson–Darling statistic, $A_m^2$, in (55), the sample of daily return $\{y_1,y_2\dots y_m\}$ of size m is arranged in ascending order: $y_{(1)}<y_{(2)}<\cdots <y_{(m)}$. The Anderson–Darling statistic (Lewis 1961) then becomes

$$\begin{array}{c}\hfill A_m^2=-m-{\displaystyle \frac{1}{m}}\sum _{j=1}^m\left[(2j-1)log(F(y_{(j)}))+(2(n-j)+1)log(F(y_{(j)}))\right].\end{array}$$

(59)

For each index, the Anderson–Darling statistic (59) is computed, along with the p-value statistic.

Table 9. Anderson–Darling statistics and p-values.

	GTS		GBN		GTS Variants		Sample Size
Index	${\widehat{\mathit{A}}}_{\mathit{m}}^{\mathbf{2}}$	p-Value	${\widehat{\mathit{A}}}_{\mathit{m}}^{\mathbf{2}}$	p-Value	${\widehat{\mathit{A}}}_{\mathit{m}}^{\mathbf{2}}$	p-Value	m
Bicoin BTC	0.1098	0.9999	99.706	0.0000	0.1105 1	0.9999	4083
Ethereum	0.1018	0.9999	59.157	0.0001	0.2123 2	0.9866	3246
S&P 500	0.3007	0.9376	54.304	0.0001	0.5010 3	0.7458	3656
SPY ETF	0.3017	0.9368	51.516	0.0001	0.6684 3	0.5857	3655

¹ Kobol distribution ($\beta ={\beta }_{-}={\beta }_{+}$); ² Carr–Geman–Madan–Yor (CGMY) distributions ($\beta ={\beta }_{-}={\beta }_{+}$; $\alpha ={\alpha }_{-}={\alpha }_{+}$); ³ bilateral Gamma distribution (${\beta }_{-}={\beta }_{+}=0$).

shows the KS statistics ($A_m^2$) and p-values for the GTS, GBM, and GTS variant distributions. While the two-parameter GBM hypothesis is always rejected, the GTS hypothesis is accepted and yields a very high p-value.

In addition, as shown by the p-value indicator in Table 9, the GTS distribution outperforms the bilateral Gamma distribution for the S&P 500 and SPY ETF indexes. However, the Kobol and CGMY distributions for Bitcoin and Ethereum, respectively, have almost the same performance as the GTS distribution.

6.3. Pearson’s Chi-Squared Test Analysis

Pearson’s chi-squared test (Schoutens 2003) counts the number of sample points falling into certain intervals and compares them with the expected number under the null hypothesis. Under the null hypothesis, $H_0$, a random sample $\{y_1,y_2\dots y_m\}$ comes from the GTS distribution, which has seven parameters estimated in Section 5. Suppose that m observations in the sample from a population are classified into K mutually exclusive classes with respective observed numbers of observations $N_j$ (for $j=1,2,\dots ,K$), and a null hypothesis gives the probability ${\mathsf{\Pi }}_j=F(x_j)-F(x_{j-1})$ (47) that an observation falls into the jth class.

The following Pearson statistic calculates the value of the chi-squared goodness-of-fit test:

$$\begin{array}{c}\hfill {\chi }^2(K-1-p)=\sum _{j=1}^K{\displaystyle \frac{{\left(N_j-m{\mathsf{\Pi }}_j\right)}^2}{m{\mathsf{\Pi }}_j}}.\end{array}$$

(60)

Under the null hypothesis assumption, as m goes to $+\infty$, the limiting distribution (Schoutens 2003) of the Pearson statistic (60) follows the ${\chi }^2(K-1-p)$ distribution with $K-1-p$ degrees of freedom, and p is the number of estimated parameters.

Table 10. Pearson statistics and p-values.

	GTS		GBN		GTS Variants			Class Number
Index	${\widehat{\mathit{\chi }}}^{\mathbf{2}}(\mathit{K}-\mathbf{8})$	p-Value	${\widehat{\mathit{\chi }}}^{\mathbf{2}}(\mathit{K}-\mathbf{3})$	p-Value	${\widehat{\mathit{\chi }}}^{\mathbf{2}}(\mathit{K}-\mathit{p}-\mathbf{1})$	$\mathit{p}$	p-Value	K
Bicoin BTC	12.234	0.508	1375	0.000	12.549 1	6	0.562	21
Ethereum	6.910	0.863	805	0.000	8.618 2	5	0.854	20
S&P 500	9.886	0.703	574	0.000	12.844 3	5	0.614	21
SPY ETF	13.955	0.377	605	0.000	18.228 3	5	0.251	21

¹ Kobol distribution ($\beta ={\beta }_{-}={\beta }_{+}$); ² Carr–Geman–Madan–Yor (CGMY) distributions ($\beta ={\beta }_{-}={\beta }_{+}$; $\alpha ={\alpha }_{-}={\alpha }_{+}$); ³ bilateral Gamma distribution (${\beta }_{-}={\beta }_{+}=0$).

shows the Pearson chi-squared statistics (${\widehat{\chi }}^2(K-1-p)$), p-values, and class number for the GTS, GBM, and GTS variant distributions. While the two-parameter GBM hypothesis is always rejected, the GTS hypothesis is accepted and yields a high p-value.

In addition, as shown by the p-value indicator in Table 10, the GTS distribution outperforms the bilateral Gamma distribution for the S&P 500 and SPY ETF indexes. However, the Kobol and CGMY distributions for Bitcoin and Ethereum, respectively, have almost the same performance as the GTS distribution. For more details on the estimation of the Pearson statistic inputs under the GTS distribution, refer to

Table A5. Observed versus expected statistics under GTS distribution.

	Bitcoin			Ethereum			sp500			SPY EFT
k	x(k)	n*${\mathbf{\Pi }}_{\mathbf{k}}$	n(k)	x(k)	n*${\mathbf{\Pi }}_{\mathbf{k}}$	n(k)	x(k)	n*${\mathbf{\Pi }}_{\mathbf{k}}$	n(k)	x(k)	n*${\mathbf{\Pi }}_{\mathbf{k}}$	n(k)
1	−18.988	7.512	8	−20.861	7.531	6	−4.341	10.264	12	−4.405	8.327	11
2	−17.080	4.144	7	−18.321	5.583	6	−3.935	5.456	5	−4.007	4.562	3
3	−15.172	6.603	7	−15.781	10.018	15	−3.529	8.442	7	−3.608	7.116	7
4	−13.264	10.678	9	−13.241	18.331	20	−3.123	13.138	15	−3.210	11.144	12
5	−11.356	17.586	13	−10.700	34.424	29	−2.717	20.588	20	−2.811	17.538	18
6	−9.448	29.661	32	−8.160	66.980	68	−2.311	32.543	30	−2.413	27.775	27
7	−7.540	51.657	47	−5.620	137.268	134	−1.905	52.023	48	−2.015	44.350	41
8	−5.632	94.188	107	−3.080	305.591	305	−1.499	84.479	89	−1.616	71.617	73
9	−3.724	184.486	168	−0.540	769.951	769	−1.093	140.406	147	−1.218	117.552	122
10	−1.816	411.503	419	2.000	965.210	966	−0.687	242.564	244	−0.819	198.101	186
11	0.092	1195.725	1186	4.540	458.955	466	−0.281	455.971	456	−0.421	351.660	348
12	2.000	1150.470	1159	7.080	219.873	222	0.126	896.809	896	−0.023	725.476	730
13	3.908	473.177	469	9.620	111.760	101	0.532	749.106	733	0.376	867.735	867
14	5.816	217.783	227	12.160	59.253	60	0.938	430.841	426	0.774	541.022	522
15	7.724	107.387	102	14.700	32.379	32	1.344	234.692	260	1.173	300.464	325
16	9.632	55.272	51	17.241	18.099	21	1.750	126.820	137	1.571	163.491	189
17	11.540	29.294	39	19.781	10.296	12	2.156	68.688	57	1.969	88.862	82
18	13.448	15.861	14	22.321	5.939	5	2.562	37.387	33	2.368	48.491	36
19	15.356	8.729	9	24.861	3.465	5	2.968	20.460	15	2.766	26.600	23
20	17.264	4.866	4		5.091	4	3.374	11.256	12	3.165	14.669	13
21		6.419	6					14.067	14		18.450	20

in Appendix A.2.

7. Conclusions

This study provides a methodology for fitting the rich class of the seven-parameter GTS distribution to financial data. Four historical prices were considered in the methodology application: two heavy-tailed data (Bitcoin and Ethereum returns) and two peaked data (S&P 500 and SPY ETF returns). The study used each historical data to fit the seven-parameter GTS distribution to the underlying data return distribution. The advanced fast FRFT scheme, based on the classic fast FRFT algorithm and the 11-point composite Newton–Cotes rule, was used to perform the maximum-likelihood estimation of seven parameters of the GTS distribution. The maximum likelihood estimate results show that, for each index, the location parameter, $\mu$, is negative, while the others are positive, as expected in the literature. The statistical significance of the parameters was analyzed. The non-statistical significance of the index of stability parameters (${\beta }_{+}$, ${\beta }_{-}$) has led to the fitting of the Kobol, CGMY, and bilateral Gamma distributions. The goodness of fit was assessed through Kolmogorov–Smirnov, Anderson–Darling, and Pearson’s chi-squared statistics. While the two-parameter GBM hypothesis is always rejected, the goodness-of-fit analysis shows that the GTS distribution fits significantly the four historical data with a very high p-value.

As a main limitation of the study, the applied methodology is computation-intensive, and the researchers need good skills in computer programming. In future work, the estimated parameter of the GTS distribution will be used in the Ornstein–Uhlenbeck-type process to simulate the daily cumulative returns of financial assets.