**1. Introduction**

The problem of calculating the density of a strictly stable law with the characteristic function

$$\vec{g}(t, u, \theta, \lambda) = \exp\{-\lambda |t|^{\alpha} \exp\{-i\frac{\pi}{2} \alpha \theta \operatorname{sign} t\}\}, \quad t \in \mathbb{R}, \tag{1}$$

where 0 < *α* 2, |*θ*| min(1, 2/*α* − <sup>1</sup>), *λ* > 0 is considered in the paper. One of the reasons why it became necessary to calculate the density of a strictly stable law with this characteristic function is the need to calculate the density of a fractionally stable law which is defined by the expression

$$g(\mathbf{x}, \mathbf{a}, \nu, \theta) = \int\_0^\infty g(\mathbf{x} y^{\nu/a}, \mathbf{a}, \theta) g(y, \nu, \mathbf{1}) y^{\nu/a} dy. \tag{2}$$

Here, *g*(*<sup>x</sup>*, *α*, *θ*) and *g*(*y*, *ν*, 1) are the densities of strictly stable and one-sided strictly stable laws with the characteristic function in Equation (1) and parameter *λ* = 1. For the first time, the density in Equation (2) was obtained in the article [1]. The density in Equation (2) go<sup>t</sup> its name in the work [2], since the random variable *<sup>Z</sup>*(*<sup>α</sup>*, *ν*, *θ*) distributed by this law is defined by the ratio *<sup>Z</sup>*(*<sup>α</sup>*, *ν*, *θ*) = *<sup>Y</sup>*(*<sup>α</sup>*, *<sup>θ</sup>*)[*V*(*<sup>ν</sup>*, <sup>1</sup>)]−*<sup>ν</sup>*/*<sup>α</sup>*. Here the random variables *<sup>Y</sup>*(*<sup>α</sup>*, *θ*) and *<sup>V</sup>*(*<sup>ν</sup>*, 1) are distributed by the laws *g*(*<sup>x</sup>*, *α*, *θ*) and *g*(*y*, *ν*, <sup>1</sup>), respectively.

The density Equation (2) appears as a limit distribution with the following random walk scheme. Let the particle be at the origin *x* = 0 at the initial time *t* = 0 and it stays at this point during random time *T*1. Then, it instantly moves with an equal probability to the right or left at random distance *X*1 and it stays at rest again random time *T*2. Then, the whole process is repeated in the same way. Values *Xi*, *i* = 1, 2, ... are independent identically distributed random variables belonging to the domain of normal attraction of a strictly stable law *g*(*<sup>t</sup>*, *α*, *θ*) with the characteristic function in Equation (1). Values *Ti*, *i* = 1, 2, ... are independent both between themselves and of the sequence {*Xi*} by identically distributed random variables belonging to the domain of normal attraction of a strictly stable law *g*(*<sup>t</sup>*, *ν*, 1) with the characteristic function in Equation (1) and 0 < *ν* 1. We will form the sum of these random variables *<sup>S</sup>*(*τ*) = ∑ *<sup>N</sup>*(*τ*) *i*=1 *Xi*, *τ* > 0, where *<sup>N</sup>*(*τ*) is the counting process: *<sup>N</sup>*(*τ*) = max {*n* : ∑*n i*=1 *Ti <sup>τ</sup>*}. The physical interpretation of the sum *<sup>S</sup>*(*τ*) is the coordinate of particle *x* at time *τ*. In the works [1,2] it has been shown that the asymptotic (at *τ* → ∞) distribution of the sum *<sup>S</sup>*(*τ*) is described by the distribution Equation (2).

The described random walk scheme is called Continuous Time Random Walk (CTRW). For the first time it was considered in the work [3]. Later, it was described in the works [4–6]. For more detailed familiarity with this model one can look through the overviews [7,8]. In the work [1] it has been shown that the asymptotic (at *t* → ∞) distribution of the CTRW process is described with the distribution (2). In the work [9], it has been shown that the CTRW process in large time asymptotics is described with the fractional-differential equation of diffusion. The solution of this equation is expressed through fractional-stable distributions in Equation (2). This is one of the factors determining the interest in studying the class of the fractional-stable laws.

Another factor is the appearance of these distributions in various processes occurring in a plasma [10,11] or in biology processes [12–14]. In particular, the fractional-stable distributions were used to describe a distribution of the gene expression in cells of tissues of various organisms in the following papers [12–14]. It is known that the distribution of the gene expression is described by laws with the power decrease in the density [15–17]. Since the density in Equation (2) decreases as *x*<sup>−</sup>*α*−<sup>1</sup> at *x* → <sup>∞</sup>, therefore this class of distributions was used to describe the gene expression distribution. In the articles [12,13] the fractional-stable distributions were used to describe the gene expression obtained with the microarray technology. In the paper [14], these distributions were used to describe the results obtained with the Next Generation Sequence technology. In the papers [12–14] the Monte Carlo method was used to calculate the density *q*(*<sup>x</sup>*, *α*, *ν*, *<sup>θ</sup>*). To estimate the parameters (*<sup>α</sup>*, *ν*, *θ*) of the fractional-stable law a method described in [18] was used which is also based on the Monte Carlo method. However, to construct more effective estimators of the parameters (*<sup>α</sup>*, *ν*, *<sup>θ</sup>*), for example, the maximum likelihood estimation, one should be able to calculate the density *q*(*<sup>x</sup>*, *α*, *ν*, *<sup>θ</sup>*). As a result we come again to the necessity of calculating the integral of Equation (2).

As we can see from Equation (2), the density of a fractional-stable law is defined using the Mellin convolution of two strictly stable densities with the characteristic function in Equation (1). Hence, to calculate the density in Equation (2) it is necessary to be able to calculate densities *g*(*<sup>x</sup>*, *α*, *θ*) for any admissible set of parameters (*<sup>α</sup>*, *<sup>θ</sup>*). It should be pointed out that the problem of calculating densities of stable laws at present is well studied. The solution to this problem is based on the inverse Fourier transform of the characteristic function of a stable law. There are several methods for performing the inverse Fourier transform: a direct calculation of the inverse Fourier transform [19–26], the use of the fast Fourier transform algorithm [27,28], the use of the inversion formula followed by the numerical calculation of the integral [29,30], and the use of the inversion method by V. Zolotarev [31–34].

Direct implementation of the inverse Fourier transform of the characteristic function of the stable law leads to the appearance of special functions. As a rule, such a transformation can be implemented if the shift parameter of the stable law *γ* = 0. Therefore, practically all cases when it is possible to express the density of a stable law through special functions are referred to strictly stable laws. In addition, the density of a strictly stable law can be obtained only for rational values of characteristic exponent *α* and parameter of skewness. For instance, in the works [19,20], representations were obtained for the densities of stable laws through the Fox H-function. Representations for the densities of stable laws through an incomplete hypergeometric function were obtained in the article [21]. Later, the results of the article [21]

were generalized in [22] in which representations were obtained through Meijer's G-function. In the work [23] the expressions for density were obtained through the Fox *H*-function and hypergeometric function. In the work [24] expressions for the density of a one-sided stable law were obtained at 0 < *α* < 1 through the hypergeometric function. In the subsequent work [25] using the law of duality and the Mellin transform, the authors generalized the result to the case of two-sided distributions (− ∞ < *x* < ∞) and to the range of values 0 < *α* 2, where *α* is a rational number. It has been mentioned earlier that it is possible to directly implement the inverse Fourier transform of the characteristic function if the shift parameter *γ* = 0. An exception may be represented by the work [26]. In this work, the author was able to invert the characteristic function of the stable law for arbitrary values of the shift parameter and scale. As a result, for *α* ∈ (1, 2] it was possible to express the density of stable laws through a generalization of the Srivastava-Daoust of Kampé de Fériet two-variable hypergeometric function.

In the work [27], to invert the characteristic function of the stable law, the fast Fourier transform (FFT) algorithm is used. Using the FFT algorithm allows one to quickly reverse the characteristic function of the stable law and obtain numerical values of the density. However, the FFT algorithm allows one to calculate density values only on a grid of equally spaced coordinate values. This is not always convenient, since one should use interpolation methods to calculate density values at intermediate points. In the paper [28], standard quadrature numerical integration algorithms are redefined to invert the characteristic function. In the proposed approach, the FFT algorithm is used to calculate the value of the integrand at the nodes of the grid. This approach makes it possible to reduce the approximation error in the central part of the distribution. To calculate the density in the tails of the distribution, the Bergström expansion of the density of a stable law in a series is used [35] (see also § 2.4 in [32]). However, the accuracy of the proposed method depends on the values *α* and *β* and turns out to be effective only with values *α* ∈ (1, 2].

In the papers [29,30], the inversion formula is used to calculate the density of a stable law

$$\log(\mathbf{x}, a, \boldsymbol{\beta}) = \frac{1}{\pi} \Re \int\_0^\infty e^{it\mathbf{x}} \boldsymbol{\hat{g}}(t, a, -\boldsymbol{\beta}) dt = \frac{1}{\pi} \int\_0^\infty \cos(h(t, \mathbf{x}, a, \boldsymbol{\beta})) e^{-t^\mathbf{x}} dt,\tag{3}$$

where *g*<sup>ˆ</sup>(*<sup>t</sup>*, *α*, *β*) is the characteristic function of a stable law with the scale parameter *λ* = 1 and shift parameter *γ* = 0. In this case, the density *g*(*<sup>x</sup>*, *α*, *β*) is expressed through the improper integral of real variables. To calculate it, one can use standard algorithms of numerical integration. This approach was used in the work [29] where the characteristic function was chosen as

$$\mathcal{G}(t, a, \beta, \lambda, \gamma) = \begin{cases} \exp\left\{it\lambda\gamma - \lambda|t|^a + it\lambda\left(|t|^{a-1} - 1\right)\beta\tan(\pi a/2)\right\}, & a \neq 1, \\\exp\left\{it\lambda\gamma - \lambda|t|^a - it\lambda\beta(2/\pi)\ln|t|\right\}, & a = 1. \end{cases} \tag{4}$$

Here the parameters vary within 0 < *α* 2, −1 *β* 1, − ∞ < *γ* < <sup>∞</sup>, *λ* > 0. As it was pointed out in the paper, this approach does not have difficulty with the values *α* > 1.1. Difficulties with calculation arise at *α* < 0.75, *α* ≈ 1, and *β* -= 0. In addition to it, at greater values of *x* the integrand begins to oscillate fast which leads to difficulties in numerical integration. In the paper [30], it is proposed to use an optimized generalized Gaussian scheme of numerical integration to calculate the integral of Equation (3) with the characteristic function in Equation (4). In this work the constants *<sup>B</sup>*<sup>∞</sup>80 and *<sup>B</sup>*<sup>∞</sup>40 were introduced (more detailed information about the definition of these constants see [30]). If *β* = 0 the proposed integration scheme is effective at 0.5 *α* 2 and |*x*| *<sup>B</sup>*<sup>∞</sup>40. If *β* -= 0 the scheme is effective for values *α* ∈ [0.5, 0.9] ∪ [1.1, 2.0] and |*x* − *ζ*| *<sup>B</sup>*<sup>∞</sup>80. With the values of |*x*| > *<sup>B</sup>*<sup>∞</sup>40 and |*x* − *ζ*| > *<sup>B</sup>*<sup>∞</sup>80 an asymptotic expansion of the density is used in a series. With the values *α* ∈ (0.9, 1.1), *β* -= 0 and *α* ∈ (0, 0.5), *β* ∈ [−1, 1] the scheme is not applicable.

The use of Equation (3) leads to the appearance of fast oscillating functions under the sign of the integral. To ge<sup>t</sup> around this problem, in the paper [31], Zolotarev V.M. developed a method of inverting the characteristic function of a stable law. Using this method, in the paper [31] (see also § 2.2 in [32], § 4.4 in [36]) an integral representation was obtained for the density of a stable law with the characteristic function

$$\mathcal{G}(t, a, \beta, \lambda, \gamma) = \begin{cases} \exp\left\{it\lambda\gamma - \lambda|t|^a \exp\left\{-i\frac{\pi}{2}\beta \mathcal{K}(a)\operatorname{sign}t\right\}\right\}, & a \neq 1, \\\exp\left\{it\lambda\gamma - \lambda|t|^a \left(\frac{\pi}{2} + i\beta \log|t|\operatorname{sign}t\right)\right\}, & a = 1, \end{cases} \tag{5}$$

where 0 < *α* 2, −1 *β* 1, *λ* > 0, − ∞ < *γ* < <sup>∞</sup>, *<sup>K</sup>*(*α*) = *α* − 1 + sign (1 − *<sup>α</sup>*). The obtained integral representation of the density of the stable law is expressed through a definite integral. It is not possible to calculate this integral analytically. However, using the methods of numerical integration, it is possible to calculate and obtain the probability density and the distribution function of a stable law. Using the specified method of inverting the characteristic function in the work [33] integral representations for probability density and distribution functions of a stable law with the characteristic function in Equation (4) were obtained. In the paper [37], a slight modification of the characteristic function in Equation (4) is considered and it is noted that for calculation purposes it is more convenient to use this particular modification. Subsequently, this integral representation formed the basis of various software packages for calculating the probability density and distribution function of stable laws [38–42]. In the paper [43], it is indicated that difficulties in calculating the integral in the integral representation obtained in [33] arise with (1) small values of *α* and *x* → 0, (2) *x* → ∞ and (3) *α* close either to 1 or 2. In this paper, the authors proposed a method of solving the last two problems for symmetric stable laws and note that using the proposed approach, it is possible to calculate the densities of stable laws for values *α* close to either 1 or 2 as well as at *x* → 0 and *x* → ∞.

Having slightly modified Zolotarev's method [31,32] of inverting the characteristic function in the paper [34] expansions were obtained for the density of stable laws in power series. Investigating trans-stable distributions, the authors obtained expansions in the power series of densities of stable laws for the cases 0 < *α* < 1 and 1 < *α* < 2. In each of the ranges 0 < *α* < 1 and 1 < *α* < 2 expansions are represented in the form of "internal" (*x* → 0) and "external" (*x* → ∞) expansions. To describe the behavior of the density of a stable law in the whole range of values 0 < *x* < ∞ these two expansions are put together.

Thus, all the results related to obtaining expressions for the probability density of stable laws were obtained for laws with characteristic functions in Equations (4) and (5). However, to calculate the density in Equation (2), it is necessary to have an expression for the probability density *g*(*<sup>x</sup>*, *α*, *θ*) with a characteristic function Equation (1). It should be emphasized that an integral representation for the density of a stable law with the characteristic function in Equation (1) is presented in the paper [44]. However, the expression cited is valid only for *x* > 0 and *α* -= 1. In this paper, we will obtain an integral representation for the density and distribution function of a stable law with a characteristic function Equation (1) for arbitrary *x* and any admissible values of parameters *α* and *θ*.
