The Calculation of the Density and Distribution Functions of Strictly Stable Laws

Integral representations for the probability density and distribution function of a strictly stable law with the characteristic function in the Zolotarev’s “C” parametrization were obtained in the paper. The obtained integral representations express the probability density and distribution function of standard strictly stable laws through a definite integral. Using the methods of numerical integration, the obtained integral representations allow us to calculate the probability density and distribution function of a strictly stable law for a wide range of admissible values of parameters ( α , θ ) . A number of cases were given when numerical algorithms had difficulty in calculating the density. Formulas were given to calculate the density and distribution function with an arbitrary value of the scale parameter λ .

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 X i , i = 1, 2, . . . are independent identically distributed random variables belonging to the domain 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 g(x, α, β) = 1 π ∞ 0 e itxĝ (t, α, −β)dt = 1 π ∞ 0 cos(h(t, x, α, β))e −t α dt, whereĝ(t, α, β) is the characteristic function of a stable law with the scale parameter λ = 1 and shift parameter γ = 0. In this case, the density g(x, α, β) 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 aŝ g(t, α, β, λ, γ) = exp itλγ − λ|t| α + itλ |t| α−1 − 1 β tan(πα/2) , α = 1, exp {itλγ − λ|t| α − itλβ(2/π) ln |t|} , α = 1.
Here the parameters vary within 0 < α 2, −1 β 1, −∞ < γ < ∞, λ > 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 B ∞ 80 and B ∞ 40 were introduced (more detailed information about the definition of these constants see [30]). If β = 0 the proposed integration scheme is effective at 0. 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 get 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 g(t, α, β, λ, γ) = exp itλγ − λ|t| α exp −i π 2 βK(α) sign t , α = 1, exp itλγ − λ|t| α π 2 + iβ log |t| sign t , where 0 < α 2, −1 β 1, λ > 0, −∞ < γ < ∞, K(α) = α − 1 + sign (1 − α). 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][39][40][41][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(x, α, θ) 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 θ.

Auxiliary Results
Thus, the objective is to obtain an integral representation of the density of a strictly stable law with a characteristic function Equation (1). Without losing generality we will further assume everywhere that λ = 1. A strictly stable law with a parameter λ = 1 is commonly called the standard strictly stable law. An abbreviated notation of the characteristic function is accepted for standard strictly stable lawŝ g(t, α, θ, 1) ≡ĝ(t, α, θ), for density g(x, α, θ, 1) ≡ g(x, α, θ), for the distribution function G(x, α, θ, 1) ≡ G(x, α, θ), and random variable Y(α, θ, 1) ≡ Y(α, θ). Everywhere below, for standard strictly stable laws, we will use this notation. To obtain an integral representation, we use the method of inverting the characteristic function of a stable law for the first time proposed by V. Zolotarev in the work [31] and described in detail in his monograph [32] (see also § 4.4 in [36]). To prove the main theorem, we will need some auxiliary results.
Property 1 (Property of inversion). For any admissible set of values of parameters (α, θ) Proof. The proof of this property is simple enough. Applying the definition of the characteristic function in Equation (1) and by making the substitution of a variable t → −τ, we obtain It follows directly from here Equation (6).
In terms of the characteristic functionĝ(t, α, θ), probability density functions g(x, α, θ) and distribution functions G(x, α, θ) of a strictly stable law the property of inversion is written in the form The utility of this property consists in the fact that owing to this property it is sufficient to consider the issue of the density representation g(x, α, θ) or the distribution function G(x, α, θ) only for x 0 or for θ 0. For negative values of the argument x or the parameter θ expressions can be obtained according to the expressions given earlier.
The following property will be useful further.

Property 2.
For any two admissible sets of parameters (α, θ, λ) and (α, θ, λ ), there is such a unambiguously defined real a > 0, that Y(α, θ, λ) d = aY(α, θ, λ ). For the characteristic function in Equation (1), the value a is connected with parameters in the following way a = (λ/λ ) 1/α . This property is a full analog of property 2.1 in [32] (see also § 3.7 in [36]) formulated for strictly stable random variables with the characteristic function in Equation (1). This property is proved in the same way to the one which is performed in [32]. In the particular case that is of interest to us λ = 1, we obtain We now formulate a lemma which makes it possible to perform the inverse Fourier transform of the characteristic function and obtain the density of a strictly stable law. Lemma 1. The probability density function g(x, α, θ) for any admissible values of parameters (α, θ) and any x can be obtained with the help of the inversion formulas Proof. Performing the inverse of the characteristic functionĝ(t, α, θ) we obtain Let us consider the integral I 1 . By substituting the integration variable in this integral −t → τ, we obtain Now having calculated the sum I 1 + I 2 , we will obtain As a result, we obtained the first formula in Equation (9). In order to obtain the second formula in (9) it is necessary to subtract in the penultimate manipulation the imaginary component in Equation (10).
Later, we need analytic continuation of the characteristic functionĝ(t, α, θ) in the complex plane z. We will carry out this analytic continuation in the complex plane z with a semiaxis z = t > 0 with a cut along the negative part of the real axis arg z = −π. The resulting analytic continuation of the functionĝ(t, α, θ) with a half-line t > 0 will be designated as g + (z, α, θ). Using the characteristic function in Equation (1), we obtain The idea of analytic continuation of the integrand in the formula of inversion in Equation (9) in the complex plane z and subsequent calculation of the resulting integral Γ exp{izx}g + (z, α, θ)dz underlies the method of inverting a characteristic function developed by Zolotarev V.M. in the work [31]. This integral is calculated due to such a change in the integration contour Γ at which its real part does not change (for more details see [32]). To substantiate the change in the integration contour, we need the following lemma.
Proof. Let us consider the integral I(C R ) = C R e izx g + (z, α, −θ)dz. As an integration contour C R we will consider contour lines that represent an arc of a circle of radius R which has ϕ 1 arg z ϕ 2 or The task is to determine under what conditions imposed on the contour C R and the parameters α and θ limits lim R→0 I(C R ) = 0 and lim R→∞ I(C R ) = 0.
For any contour C R the inequality is true. Assuming in this expression z = re iϕ and taking into account that C R is an arc of a circle of radius R , we obtain where U(r, ϕ, α, θ) = exp {−rx sin ϕ − r α cos(α(ϕ + πθ/2)) + ln r} .
Taking into account that 0 < α < 1 and −1 θ 1, we deduce that −π/2 < αθπ/2 < π/2 and, therefore, cos(αθπ/2) > 0. Assuming that R → ∞ in this expression and in view of the fact that R α grows faster than ln R, we deduce lim R→∞ U(R, 0, α, θ) = 0. Without loss of generality the case ϕ = π can be excluded from consideration. Hence, ϕ 1 = 0, and ϕ 2 we represent in the form ϕ 2 = π − ε, where 0 < ε π/6 is an arbitrary fixed number. As a result, we deduce and integration contour in Equation (13) takes the form Now assuming that R → ∞ in Equation (14) and using Equation (16) we deduce From here it follows that in the case α < 1 and x > 0 the integral I(C R ) = 0, at R → ∞ where contour integration has the form (17). The first item of the lemma is proved.
(3) Now let us consider the case 0 < α 2 and R → 0. It is known that for a specified range of values of the parameter α the parameter θ can take the values |θ| min(1, 2/α − 1). Assuming that R → 0 in Equation (14) we obtain From here we can see that the behavior of this integral at R → 0 will be determined by the behavior U(R, ϕ, α, θ) at R → 0. Using Equation (15), we obtain For any values of ϕ and any x. Choosing ϕ 1 = −π and ϕ 2 = π a contour C R will take the form Now using Equations (26) and (27) in Equation (25) we will finally obtain lim R→0 |I(C R )| lim R→0 π −π U(R, ϕ, α, θ)dϕ = 0, for any admissible values of parameters α and θ and for any x. It proves the third item of the lemma.
(4) We will consider the case x = 0 and R → ∞. In this case, the expression in Equation (15) will take the form Thus, the integral in Equation (14) will tend to zero at Since no additional limitations for parameters α and θ were introduced here, then this result is true for any admissible values of these parameters. It should be pointed out that there are strict inequalities here. In fact, if ϕ = ± π 2α − πθ 2 , then cos(α(ϕ + πθ/2)) = 0 and in this case U(R, ϕ, α, θ) → ∞ at R → ∞. Choosing in Equation (14) the values ϕ 1 = −π/(2α) − πθ/2 + ε and ϕ 2 = π/(2α) − πθ/2 − ε as the limits of integration, contour integration C R at x = 0 takes the form where ε is an arbitrary small positive number. Now using Equations (28) and (30) in Equation (14), we obtain at x = 0 and any admissible values of parameters α and θ. It proves the fourth item of the lemma. (5) Let us consider the case α = 1. In this case the parameter θ can vary within the limits −1 θ 1. It is necessary to determine under which conditions the integral in Equation (14) tends to zero at R → ∞. As in previous cases, assuming that R → ∞ in Equation (14), we obtain From this it follows that these conditions are a consequence of the behavior U(R, ϕ, 1, θ) at R → ∞. Applying Equation (15), we obtain In view of the fact that R grows faster than ln R at R → ∞ we obtain that This inequality allows us to determine the conditions imposed on ϕ. It should be noted that one should exclude the case x = 1, θ = 1 from consideration. In fact, substituting these values in Equation (32), we obtain lim R→∞ U(R, ϕ, 1, 1) = ∞ for any ϕ. Thus, from Equation (31) we obtain Making some transformations the inequality in Equation (33) can be written in the form From this inequality we can see that it is necessary to consider three cases: x > 1, 0 x < 1, and x = 1.
Now using these two trigonometric identities in Equation (40), we find Thus, at x = 0 the condition in Equation (39) will take the form − π 2 − πθ 2 < ϕ π 2 , if α = 1 and θ < 0. The condition obtained is the same as the condition in Equation (29), if to limit the latter above with the value π/2 at θ < 0. Now we need to consider the case θ = θ 0 . In Figure 2 we can see that when increasing the parameter θ the point ϕ 0 (θ, x) will approach the value −π/2. With such a change of the parameter θ the function graph f (ϕ, θ, x) in Figure 1 will shift to the left and with the value θ = θ 0 will take the form that is given in Figure 1. As we can see from Figure 2, in this case the function ϕ 0 (θ, x) has a discontinuity in the point θ = θ 0 . It is connected with the fact that the principal branch of the function arctan x is investigated. In the vicinity of this point we have From this we can see that the point θ 0 is the point of discontinuity of the first kind. It is possible to eliminate the discontinuity of the function ϕ 0 (θ, x) in the point θ 0 by defining the value of this function in the given point. We will select Thus, the condition in Equation (33) is met if It should be noted that there are precisely strict inequalities here. In fact, if to take ϕ = ±π/2, then in the case considered these values will be the solution of the Equation (38) which will lead to divergence of the integral in Equation (31). Now we will consider the case θ > θ 0 . From Figure 1 one can see that by increasing the parameter θ from θ 0 to 1 the half-period of the function f (ϕ, θ, x) satisfying the condition f (ϕ, θ, x) > 0 will keep moving to the left. At the same time, the left point which is the solution of an equation f (ϕ, θ, x) = 0 will become smaller than −π/2. Since we take interest in the interval from −π/2 to π/2, then the left bound of the interval will be the value −π/2. Keeping in mind that the half-period of the function f (ϕ, θ, x) is equal to π, the right bound of the interval f (ϕ, θ, x) > 0 will be the second solution of the equation f (ϕ, θ, x) = 0. As a result, the case considered the condition in Equation (33) will be met with values ϕ satisfying the inequality The graph of the function ϕ 0 (θ, x) is given in Figure 2.
If x = 0, we have θ 0 = 0 and condition θ > θ 0 turns into a condition θ > 0. As it was pointed out earlier, in this case it is convenient to represent ϕ 0 (θ, x) in the form of Equation (40). Now using the trigonometric identities in Equations (41) and (42) in Equation (40) we obtain Thus, at x = 0 the condition in Equation (46) takes the form − π 2 ϕ < π 2 − πθ 2 , if x = 0 and θ > 0. The condition obtained is the same as the condition in Equation (29), if to limit the latter below with a value −π/2 at θ > 0.
Consider now the case x = 1. In this case the inequality in Equation (35) has the form From this it follows that 1 − sin(πθ/2) > 0, for any θ < θ 0 , where θ 0 = (2/π) arcsin 1 = 1. As a result, selecting the principal branch of the function arctan y, the inequality in Equation (48) takes the form This inequality gives a condition under which the inequality will be satisfied in Equation (35) in case x = 1. The graph of the function ϕ 0 (θ, 1) is given in Figure 2.

Remark 1.
The statements of the Lemma 2 are preserved if to substitute contours C R with their parts.
Using the lemma which was proved one can substantiate the validity of transition from an integral ∞ 0 exp{itx}ĝ(t, α, θ)dt along the real variable t to an integral Γ exp{izx}g + (z, α, θ)dz in the complex variable z along some contour Γ for inversion formulas in Equation (9). We will state this result in the form of a lemma. Lemma 3. Let us consider the family of contours {Γ} in the complex plane z with a cut along a half-line arg z = −π satisfying the following conditions: 1. Every contour starts in the point z = 0. 2. None of the contours Γ intersects the lines of the cut. 3. Moving from the point z = 0 along the contour Γ we let it tend to infinity but in such a way that starting from some place all points z ∈ Γ have values of arguments within the limits: where θ 0 = (2/π) arcsin x and ϕ 0 (θ, x) have the form in Equation (12) and ε > 0 is any arbitrary small number. Then for any contour of the specified type and any pair of admissible parameters (α, θ) and any x 0 (with the exception of the point x = 1, α = 1, θ = 1) where t is real.
Proof. From constraints imposed on contours Γ one can see that it is possible to divide the whole family of contours {Γ} into two kinds. The contours which start in the point z = 0 and tend to infinity without intersecting the positive part of a real semiaxis are referred to the contours of the first kind. The contours intersecting the positive part of a real semiaxis are referred to the contours of the second kind. We will consider, at first, contours of the first kind. We will introduce the following notation: z r is the intersection point of the contour Γ with a circle |z| = r, C r is an arc of a circle of radius r (not crossing the cut) which is formed when moving from the point z = r to the point z r and Γ r,R is a part of a contour Γ which is formed when moving from the point z r to the point z R . We form a closed contour G r,R = [r, R] ∪ C R ∪Γ r,R ∪ C r (see Figure 3). The line means that we go along the contour in the opposite direction. Since the function h(z) = e izx g + (z, α, −θ) is analytic in the region restricted by the contour G r,R , then by using the Cauchy theorem We will assume in this expression that r → 0 and R → ∞. Using Lemma 2 and Remark 1 we find that C R h(z)dz + C r h(z)dz → 0, at R → ∞ and r → 0. Therefore, the equality in Equation (59)  Now consider the contours of the second kind. These contours are characterized by the feature that they intersect the real axis. Therefore, to prove the lemma, we consider two closed auxiliary contours: the contour G r,x = [r, x] ∪Γ r,x ∪ C r and contour G x,R = [x, R] ∪ C R ∪ Γ x,R (see Figure 4). Here x is the intersection point of a contour Γ with a real axis. Since the function h(z) is analytic within the regions restricted with the contours Γ r,x and Γ x,R , then by using the Cauchy theorem Now we will assume in these expressions that r → 0 and R → ∞. Using Lemma 2 and Remark 1, we find that C r h(z)dz → 0 at r → 0 and C R h(z)dz → 0 at R → ∞. Now summing up Equations (60) and (61) in view of this result we obtain hence, the equality in Equation (59) is true. Figure 4. Auxiliary contours G r,x (heavy curve) and G x,R (dashed heavy curve).
From here it follows Equation (64). The lemma is completely proved.
In Figures 5 and 6, the graphs of the functions ϕ 0 (θ, x) and ϕ(θ, x) are given for the cases x 1 and 0 x < 1 which clearly illustrate the validity of the lemma that has just been proved.
ϕ0(θ, x) Figure 5. The graph of the function ϕ(θ, x) and ϕ 0 (θ, x) depending on the parameter θ in the case of x 1 (The graphs are plotted for x = 1.5 and x = 1).

Main Results
Now we can formulate the main theorem which gives an integral representation for the probability density of a stable law g(x, α, θ) with the characteristic function in Equation (1). Theorem 1. The distribution density g(x, α, θ) of a strictly stable law with a characteristic function as in Equation (1) can be represented in the form 1. If α = 1 and x = 0 for any values |θ| min(1, 2/α − 1) where θ * = θ sign x and U(ϕ, α, θ) = sin α ϕ + π 2 θ cos ϕ 2. If x = 0, then for any 0 < α 2 and |θ| min(1, 2/α − 1) g(0, α, θ) = 1 π cos (πθ/2) Γ (1/α + 1) 3. If α = 1, then for any |θ| 1 and any values x Proof. To obtain the expression for the probability density we use the inversion formulas in Equation (9). In principle, it makes no difference which formula to use. The result will differ then only by the sign of the parameter θ. We will use the first formula in Equation (9), we have Without loss of generality, we assume that x > 0. The density g(x, α, θ) for x < 0 can be obtained with the inversion property in Equation (7). Next, let us make the substitution of integration variable t → z, where z is complex-valued. Such a substitution means that we analytically extend the integral to the complex plane. With this analytical continuation is carried out starting with the positive part of a real semiaxis. As a result, we have thatĝ(t, α, −θ) → g + (z, α, −θ), where g + (z, α, θ) is defined by Equation (11). The improper integral becomes the integral along the contour Γ. We will define the contour Γ in such a way that e izx g + (z, α, θ) = 0, and the contour itself Γ should start in the point z = 0 tend to infinity. Since in the inversion formula of Equation (73) the variable is t > 0 then from this it follows that arg z must lie within the limits −π/2 arg z π/2. As result, the contour Γ will take the form However, the specific type of the contour has to be defined. In view of the foregoing, the expression in Equation (73) will be written in the form As a result, the problem consists in determining the contour form Γ, in proving the validity of transition from Equation (73) to Equation (75) and in calculating this integral.
1. The contours of the first group are made up of the contours with values 0 < α < 1 and θ = −1.
In Figure 7 the contour Γ 1 corresponds to this case. From the definition of the contour in Equation (80) one can see that in this case the admissible region of an angle ϕ takes the form π 2 ϕ π 2 . This means that the contour goes along the positive part of the imaginary axis: Γ ≡ Γ 1 = I + .
2. The contours of the second group include the contours with values 0 < α < 1 and −1 < θ 0. In Figure 7 the contours Γ 2 , Γ 2 , Γ 2 , Γ 2 are referred to this case. The contours of this group start in the point z = 0 at ϕ = − π 2 θ and tend to infinity at ϕ → π 2 . As one can see, in this case − π 2 θ 0, and contours of this group do not cross the real semiaxis. 3. The third group is made up of the contours with values of parameters 0 < α < 1 and 0 < θ < 1.
In Figure 7 this group consists of the contours Γ 3 , Γ 3 , Γ 3 , Γ 3 . The contours of this group start in the point z = 0 coming out at an angle ϕ = − π 2 θ, and tend to infinity at ϕ → π 2 . As we can see in this case − π 2 θ < 0. Therefore, the contours of this group approach the point z = 0 at values of ϕ < 0 which, in its turn, means that these contours intersect the positive part of the real axis. 4. The fourth group is composed of the contours 0 < α < 1, θ = 1. In Figure 7 the contour Γ 4 . corresponds to this case. From the expression in Equation (78) we can see that in this case at ϕ → − π 2 in this expression there is an indeterminate form 0/0. Evaluating this indeterminate form according to L'Hôpital's rule we get lim ϕ→− π 2 r(ϕ) = (α/x) 1/(1−α) . Thus, the contours of this group start in the point z = −i(α/x) 1/(1−α) at ϕ = − π 2 and tend to infinity at ϕ → π/2. As we can see the contours of this group also cross the positive part of the real axis. 5. The fifth group includes the contours with parameters 1 < α 2 and −(2/α − 1) θ 0. In Figure 8 the contours Γ 5 , Γ 5 , Γ 5 , Γ 5 are referred to this case. The contours of this group start in the point z = 0 at ϕ = π/2 and tend to infinity at ϕ → − π 2 θ. Since in this case θ 0, then the condition ϕ > 0 is met for all points of the contour. Therefore, the contours of this group do not cross the positive part of the real axis. 6. The sixth group consists of the contours with parameters 1 < α 2, 0 < θ < 2/α − 1. In Figure 8 the contours Γ 6 , Γ 6 correspond to this case. The contours of this groups start in the point z = 0 at ϕ = π/2 and tend to infinity at ϕ → − π 2 θ. Since in this case θ > 0, then − π 2 θ < 0 and, consequently, the contours of this group cross the positive part of the real semiaxis. 7. The seventh group comprises the contours with parameters 1 < α 2 and θ = 2/α − 1. In Figure 8 the contour Γ 7 corresponds to this case. One should pay attention that in this case at ϕ = π/2 the function r(ϕ) defined by Equation (78) has an indeterminate form 0/0. Evaluating this indeterminate form according to L'Hôpital's rule we obtain lim ϕ→π/2 r(ϕ) = (α/x) 1/(1−α) . Thus, the contours of this group start in the point z = i(α/x) 1/(1−α) at ϕ = π/2 and tend to infinity at ϕ → − π 2 θ. The contours of this group also cross the positive part of the real axis. Figure 7. The type of a contour Γ at α ∈ (0, 1) and various values of parameter θ. The contours are given for the case α = 0.6 and Γ 1 - This indicates that the contours of the first, second and third groups (0 < α < 1, −1 θ < 1) and contours of the fifth and sixth groups (1 < α 2, −(2/α − 1) θ < 2/α − 1) satisfy the conditions of the Lemma 3. In fact, the contours of these groups start in the point z = 0, do not cross the line of the cut that goes through a half-line arg z = −π and tend to infinity. The contours of the first, second, third groups tend to infinity at ϕ = π 2 . This means that the contours of these groups satisfy the condition in Equation (51). The contours of the fifth and sixth groups tend to infinity at ϕ = − π 2 θ. Thus, these contours satisfy the condition in Equation (52).
As a result, for the case 1 < α 2, θ = 2/α − 1 we get Now, applying Lemma 3 and taking account of the equalities in Equations (82) and (83) we find that where the contour Γ is defined by the expression in Equation (80). As a result, an improper integral along the positive part of the real axis in the expression of Equation (73) can be replaced with an integral along the contour Γ. Thus, in the case considered (α = 1) we showed the validity of transition from Equation (73) to Equation (75). Returning to Equation (75), taking into consideration Equation (76) and representing complex z in the form z = re iϕ , we obtain exp{−rx sin ϕ − r α cos(α(ϕ + π 2 θ))} cos rx cos ϕ − r α sin(α(ϕ + π 2 θ)) d[r cos ϕ].
We will place here the expression in Equation (78). One should pay attention that in this case the Equation (77) is valid, consequently, we obtain where A(ϕ) = −r(ϕ)x sin ϕ − r(ϕ) α cos(α(ϕ + π 2 θ)) and r(ϕ) is defined by Equation (78). We transform the function A(ϕ) where U(ϕ, α, θ) has the form of Equation (70). Now we consider the differential, we have d[r(ϕ) Using now this result in the expression for d[r(ϕ) cos ϕ], we have where U(ϕ, α, θ) has the form of Equation (70). Now using Equations (86) and (87) in Equation (85) and also taking into consideration that the motion along the contour Γ, having the form of Equation (80) now described by the parameter change ϕ, we obtain From the expression in Equation (84) one can see that the integration limits ϕ min and ϕ max should be selected in such a way that the motion along the contour Γ could correspond to a change from r = 0 to r = ∞. From the expression in Equation (81) it is clear that in the case 0 < α < 1 a change in the angle ϕ from − π 2 θ to π 2 corresponds to the motion along the contour Γ from r = 0 to r = ∞. Therefore, in this case ϕ min = − π 2 θ, ϕ max = π 2 . In the case 1 < α 2 a change in the angle ϕ from π 2 to − π 2 θ corresponds to the motion from r = 0 to r = ∞. Therefore, in this case ϕ min = π 2 , and ϕ max = − π 2 θ. Combining these two cases the expression in Equation (88) takes the form It should be pointed out that this formula was obtained on the assumption x > 0. The case x < 0 is easy to obtain using the property of inversion for the density of probabilities in Equation (7). For this it is enough to replace the parameter θ in the formula in Equation (89) with −θ. It is possible to combine these two cases if in the formula of Equation (89) to substitute the parameter θ for the parameter θ * = θ sign x and the value x is taken in absolute value. As a result, we obtain the expression in Equation (69) valid for any x = 0 and any admissible α = 1.
Now we consider the case 0 < α 2 and x = 0. In this case the inversion formula of Equation (75) takes the form where the contour Γ is determined by the expression in Equation (74), but it is necessary to determine the specific type of a contour in the case under consideration.
We represent the complex number z in the form z = re iϕ . As a result, the integrand in Equation (90) takes the form Using the condition g + (0, α, θ) = 0, we obtain This equation has two solutions: r = 0, for any ϕ, and ϕ = −πθ/2, r 0. It is clear that if r = 0, then g + (0, α, θ) = 0 for any value of ϕ, for definiteness we will select ϕ = −πθ/2 if r = 0. As a result, the contour of integration of Equation (74) Here it was taken into account that [dz] = cos(πθ/2)dr on the contour Γ. The limits of integration were selected in such a way that when moving along the contour r could change from 0 to ∞. Since the contour Γ in the case under consideration is a half-line coming out of the point z = 0 at an angle arg z = −πθ/2, then the motion along the contour Γ from 0 to ∞ corresponds to a change r from 0 to ∞. Making in this integral a substitution of a variable r α = y and using the definition of the gamma-function Γ(n) = ∞ 0 t n−1 e −t dt, n > 0 we obtain the expression Equation (71). Thus, the second item of the theorem is proved. Now we consider the case α = 1. We will make an assumption that x 0. In this case, the inversion formula in Equation (75) takes the form where the contour of integration defined by Equation (74) will be written in the form Here it should be noted that analytic continuation of the functionĝ(t, 1, −θ) from the positive part of the real axis t to the complex plane z at α = 1 is an analytic function and it has the form g + (z, 1, −θ) = exp −z exp i π 2 θ . As in the previous case, we begin by defining the form of the integration contour Γ. Consider the integrand in Equation (93) and we represent the complex number z in the form z = re iϕ . As a result, we obtain exp{izx − z exp{i π 2 θ}} = exp ixre iϕ − re iϕ exp{i π 2 θ} = exp{−r x sin ϕ + cos((ϕ + π 2 θ)) } exp ir x cos ϕ − sin((ϕ + π 2 θ)) . (95) From this expression we get that the condition e izx g + (z, 1, −θ) = 0 leads to an equation This equation has two solutions. The first solution is r = 0. The second solution we obtain from the equation x cos ϕ − sin((ϕ + π 2 θ)) = 0, r 0. Solving it with respect to ϕ we get As a result, the contour of integration in Equation (94) takes the form Thus, the contours Γ with different values θ are half-lines coming out of the origin at an angle ϕ(θ, x) and tending to infinity.
Next, it is necessary to substantiate the transition from the integral in Equation (73) to the integral along the contour in Equation (93). For this we will use Lemmas 3 and 4. From the definition in Equation (98), it is clear that for all admissible values θ and x 0 the contour Γ satisfies item 1 and 2 of the Lemma 3. There is only one thing left, to find out if the contour in Equation (98) satisfies item 3 of this lemma.
Consider the case x > 1 at first. According to item 3 of the Lemma 3 for the equality in Equation (59) to be valid the contour of Equation (98) must satisfy the condition in Equation (54). Now using the Lemma 4. According to this lemma in the case x 1 the functions ϕ(θ, x) and ϕ 0 (θ, x) are connected between each other with a ratio (63) from which it directly follows that ϕ 0 (θ, x) < ϕ(θ, x) for all θ ∈ [−1, 1]. Therefore, at x > 1 the condition in Equation (54) is met and we can move from the integral in Equation (73) to the integral in Equation (93). In the case x = 1 the contour Γ must meet the condition in Equation (58). Now using the Lemma 4 we obtain that in this case for all θ ∈ [−1, 1] the inequality ϕ 0 (θ, x) < ϕ(θ, x) is true. Therefore, in this case the condition in Equation (58) is satisfied. Now we consider the case 0 x < 1. Similar to previous case, applying the Lemma 4 namely, the formula in Equation (64) we get that in this case the contour in Equation (98) satisfies the conditions in Equations (55)-(57) Lemma 3. Thus, the contour in Equation (98) completely satisfies the conditions of the Lemma 3 and therefore, the equality 1 π ∞ 0 e itxĝ (t, 1, −θ)dt = 1 π Γ e izx g + (z, 1, −θ)dz is true. This makes it possible to replace the improper integral in the expression of Equation (73) over the real variable with the integral over the contour Γ. Thus, the possibility for the transition from the formula in Equation (73) to the formula in Equation (75) is substantiated. We need to note that in the case under consideration α = 1 the expression in Equation (75) takes the form of Equation (93). It should also be pointed out that that due to Equation (34) the case α = 1, θ = 1, x = 1 is excluded in the Lemma 3. That is why, here, this case should also be excluded from consideration. Now taking into account Equation (95), the expression of Equation (93) will be written in the following form Now using here the definition of Equation (98) we obtain that [dz] = cos(ϕ(θ, x))dr. The motion along the contour Γ should take place in such a way that it would start in the point z = 0 in the process of moving it would tend to infinity. Therefore, the motion from r = 0 to r = ∞ corresponds to such motion. Taking into consideration that on the contour Γ the condition of Equation (96) is met, we obtain g(x, 1, θ) = 1 π ∞ 0 exp{−r x sin ϕ(θ, x) + cos(ϕ(θ, x) + π 2 θ) } cos(ϕ(θ, x))dr As we can see, the integral obtained is easy to calculate. As a result, we obtain g(x, 1, θ) = 1 π cos(ϕ(θ, x)) x sin ϕ(θ, x) + cos(ϕ(θ, x) + πθ /2) .
Using now the definition of the function ϕ(θ, x) (97) after simple transformations we get Recall that consideration was carried out for the case x 0. The case x < 0 can be obtained using the inversion property for the density of probability of Equation (7). For this it is enough to replace θ with −θ. It is possible to combine these two cases if to introduce a parameter θ * = θ sign x. However, we want to note that if to perform this replacement in the expression of Equation (99), then the expression itself will not change g(x, 1, θ) → g(x, 1, θ * ) ≡ g(x, 1, θ). Consequently, the formula in Equation (99) is true for any x. Thus, the theorem is completely proved.
We will make some remarks on the proved theorem.

Remark 2.
The proof of the case α = 1 was carried out under the assumption x 0. Therefore, the formula in Equation (99) was obtained for the case x 0. The generalization of this formula for the case x < 0 was carried out using the inversion property g(−x, 1, θ) = g(x, 1, −θ). It should be noted here that in the process of proving this case the point θ = 1, x = 1 was excluded from consideration. Therefore, in view of the inversion property, the point θ = −1, x = −1 should also be excluded. In these two points the density g(x, 1, θ) has a peculiarity. In fact, substituting the value θ = ±1, we obtain However, an indeterminate form 0/0 can be evaluated in the following remark. (71) and (72) can be obtained directly from the inversion formulas in Equation (9) without resorting to analytic continuation of the characteristic functionĝ(t, α, θ) to the complex plane with the subsequent transition from the improper integral over the real variable (73) to the integral along the contour in Equation (75). We first consider the case α = 1. Using the inversion formula (the first formula in Equation (9)) we obtain

Remark 3. Equations
We determine under which conditions this integral will converge. For the integral I the inequality is valid For the integrand we have Since in the case considered (α = 1) the parameter θ varies within the limits −1 θ 1, we obtain cos(πθ/2) > 0, if − 1 < θ < 1. From here we get lim t→∞ exp t ix − exp i π 2 θ = 0, if −1 < θ < 1. Thus, the integral in Equation (101) will converge, and, therefore, and the integral I will also converge at −1 < θ < 1.
Returning to Equation (100) and by calculating the integral directly, we get Thus the behavior of the formula in Equation (106) coincides with the behavior of the integral I in Equations (104) and (105) in the cases θ = ±1. Therefore, the formula in Equation (106) is true for any −1 θ 1. The expression in Equation (107) means that the density g(x, 1, ±1) is a degenerate distribution in the point x = ±1. In other words, Thus, the obtained expression in Equation (106) completely coincides with the one previously obtained in the Theorem 1 the density in Equation (72), and the conclusion presented in this remark is an alternative way of deducing this density.

Remark 4.
By a similar method, one can obtain the density value at x = 0. Using the first formula in Equation (9) and making a substitution of the integration variable t α = τ we get To calculate the integral obtained it is necessary to use the formula (see [45], Section 1.5. the Equation (35)) It is clear that the integral in Equation (109) completely satisfies the conditions of this integral. Consequently, using it in Equation (109) we get g(0, α, θ) = 1 π Γ(1/α + 1) cos(πθ/2). The formula obtained coincides completely with the formula in Equation (71).
We will make another useful remark.

Remark 5.
The Theorem 1 formulates an integral representation for the density of a standard strictly stable law. However, it is useful to have a formula that allows one to convert the density of a standard strictly stable law to the density of a strictly stable law with arbitrary λ. The Property 2 and, in particular, the formula in Equation (8) allows one to obtain such a formula.

Now substituting here Equation (89) we find
If we combine the cases α < 1 and α > 1 we obtain This formula defines the distribution function of stable law for the case x > 0 and α = 1. The case x < 0 is reduced to the case x > 0 using the property of an inversion, namely, the formula in Equation (7) for G(x, α, θ). As a result, we get G (−) (−x, α, θ) = 1 − G (+) (x, α, −θ), x > 0. This formula gives the distribution function for negative x. Combining the formulas for G (+) (x, α, θ) and G (−) (−x, α, θ), we obtain the formula in Equation (112) which is true for any x = 0 and α = 1.
Thus, the second item of the corollary has been proved.

Remark 7.
The proved corollary gives an integral representation for the distribution function of a standard strictly stable law. In order to get the distribution function of a strictly stable law for an arbitrary λ it is necessary to use the Remark 5. By definition G(x, α, θ, λ) = x −∞ g(y, α, θ, λ)dy. Using now the relation in Equation (111) and changing the variable xλ −1/α = τ we arrive at the relation It should be noted the density in Equation (72) and distribution function in Equation (114) is not new. These formulas were deduced by V.M. Zolotarev (see § 2.3 in [32]). As we can see, this distribution is the stable law for α = 1 and any −1 θ 1 and it is expressed in terms of elementary functions.

The Calculation of the Density and Distribution Function of a Stable Law
Integral representations for the probability density and distribution function of a stable law with the characteristic function in Equation (1) were obtained in the Theorem 1 and Corollary 1. These integral representations in Equations (69) and (112) express the probability density and the distribution function in terms of a definite integral. That is why, using the methods of numerical integration, it is possible to calculate the values of these integrals without much difficulty.
In this paper, to calculate definite integrals in Equations (69) and (112) we used the adaptive Gaussian-Kronrod numerical integration algorithm for 31 points. To implement the program for calculating the functions g(x, α, θ) and G(x, α, θ) we used the implementation of this algorithm in the gsl library (GNU Scientific Library) of version 1.8 [46]. The calculation results for the functions g(x, α, θ) and G(x, α, θ) by the Equations (69) and (112) are given in Figures 9-14.
Let us analyze the results presented in more detail. We first consider the case α < 1. The results related to this case are given in . From these figures it can be seen that when θ = 1 the probability density is concentrated on the positive semiaxis. Thus, g(x, α, 1) = 0, G(x, α, 1) = 0, if x < 0, α < 1. Similarly, for the case θ = −1 and α < 1, we obtain that the negative semiaxis is the area of concentration of the probability density. Consequently, g(x, α, −1) = 0, G(x, α, −1) = 1, if x > 0, α < 1. This result is in complete agreement with remarks 3 and 4 on page 79 of theorem 2.2.3 from the book by [32].
In the introduction, it was noted that in the work [31] (see also [32]) integral representations were obtained for the probability density and distribution function of a stable law with the characteristic function in Equation (5). In order to avoid any confusion the parameters of a stable law with the characteristic function in Equation (5) will be designated as (α , β, λ , γ). For the parameters of a strictly stable law with the characteristic function in Equation (1) we keep the notation (α, θ, λ). The parameters (α , β, λ , γ) are related to the parameters (α, θ, λ) by the relations (see [32,36]): α = α, In the case of α = 1 from the relation for the parameters θ and β we get θ = β, if α < 1, and θ = β(1 − 2/α), if α > 1. Thus, at α < 1 the parameters θ and β coincide. Therefore, the corresponding properties for the probability density and the distribution function of stable laws with characteristic functions in Equations (1) and (5) coincide. The situation slightly changes if α > 1. It can be seen from the above relation that, firstly, the admitted region of the parameter θ narrows in comparison with the admitted region of the parameter β, secondly, the parameter θ changes its sign to the opposite with respect to the parameter β. Since the parameter θ has the meaning of an asymmetry parameter, then a change of the sign of this parameter when passing from α < 1 to α > 1 will affect the form of probability density. This is clearly seen in Figures 11 and 12 for densities with extreme values of the parameter θ. Comparing densities for the values α = 0.9, θ = 1 given in Figure 11 and density for the values α = 1.2, θ = 0.66 given in Figure 12 one can see that these densities are turned into different directions. This fact is a consequence of the fact that the parameter θ changed the sign compared to the sign of the parameter β. Thus, for the same sign of the parameter θ of density g(x, α, θ) at α < 1 and α > 1 will be turned into different directions. The reason for this behavior is related to the selected parameter system (α, θ, λ) of the characteristic function in Equation (1) It should be pointed out that the Theorem 1 and Corollary 1 formulate expressions for probability density and distribution functions for strict stable laws with a scale parameter λ = 1. To obtain the density and distribution function with an arbitrary value of the scale parameter λ it is necessary to use Remarks 5 and 7.

Conclusions
In this paper, integral representations for the probability density have been obtained (Theorem 1) and distribution function (Corollary 1) of a standard (λ = 1) strictly stable law with the characteristic function in Equation (1). In the general case α = 1 and x = 0 the probability density and distribution function are expressed in terms of a definite integral. In the case α = 1 for any x and in the case x = 0 for any admissible α and θ the probability density and distribution function are expressed in terms of elementary functions. Applying the method of numerical integration, the values of the density and distribution function of strictly stable laws with the characteristic function in Equation (1) were calculated. The calculations show that the numerical methods do not have any difficulties in calculating the density and distribution function for the selected parameter values.
However, this does not mean that one can calculate the density and function of distribution for all admissible parameters by using obtained integral representations. Most likely, numerical integration algorithms will have difficulty in calculating the integral for small values α, at α ≈ 1 and for bigger values of x. The results of the works in which integral representations for densities of stable laws with characteristic functions in Equations (5) and (4) were investigated testify to this. An integral representation for a stable law with the characteristic function in Equation (5) was obtained in the work [31] (see also § 2.2 in [32], § 4.4 in [36]). In the work [33], it was pointed out that when values of α close to 1 problems arise with the numerical calculation of the integral in this integral representation. An integral representation for the density of a stable law with the characteristic function in Equation (4) was obtained in [33]. In this work, it was emphasized that when calculating the density, calculation difficulties arise at values 0 < |α − 1| < 0.02 and at values α close to zero. In the works by [38,43] the same problems are mentioned when calculating the integral in the representation obtained in the work [33]. Based on this, it should be expected that, with the above parameter values, calculation difficulties will also arise with the density and distribution functions obtained in the Theorem 1 and Corollary 1. In particular, directly from the expressions in Equations (69) and (113), it can be seen that at α close to 1, but not equal to 1, problems may arise with the numerical calculation of the integral. This is indicated by the exponent α/(α − 1). It can be seen that when α → 1 this value increases unlimitedly. Most likely, in this case, one will have to look for other ways of calculating the density and distribution function of a strictly stable law.
In conclusion, we would like to point out that the integral representation of the density g(x, α, θ) formulated in the Theorem 1 was used to calculate the density in Equation (2). To calculate the improper integral in Equation (2) we used the adaptive quadrature Gaussian-Kornord numerical integration algorithm on 15 points. We used the implementation of this algorithm in the library gsl (GNU Scientific Library) version 1.8 [46]. The calculations performed in some cases show the presence of problems of numerical integration. In particular, at x close to zero, the calculated density behaves like a periodic function. In addition, in some cases, the integration algorithm generates an integration error. All this indicates the need for additional study of the integrand function in Equation (2) and adapting this expression for numerical integration algorithms. It should be noted that the most likely causes of these difficulties may be the ones described above when calculating the density g(x, α, θ). Therefore, first of all, it is necessary to find a solution to the problems described above. To calculate the density at x close to zero and for bigger values x the most promising approach is to use an expansion of the strictly stable density in the power series. The method described in the article [43] can be used to calculate the density at α → 1. However, the possibility of using this approach requires additional research.
Funding: This work was supported by the Ministry of Higher Education and Science of the Russian Federation (project No. 0830-2020-0008).