The calculation of the distribution function of a strictly stable law at large
X The calculation of a strictly stable law at $x\to\infty$ V. V. Saenko V. V.
Saenko Ulyanovsk State University, S.P. Kapitsa Research Institute of
Technology, 42, Leo Tolstoy St., Ulyanovsk, 432017, email address:
<EMAIL_ADDRESS>The paper considers the problem of calculating the
distribution function of a strictly stable law at $x\to\infty$. To solve this
problem, an expansion of the distribution function in a power series was
obtained, and an estimate of the remainder term was also obtained. It was
shown that in the case $\alpha<1$ this series was convergent for any $x$, in
the case $\alpha=1$ the series was convergent at $N\to\infty$ in the domain
$|x|>1$, and in the case $\alpha>1$ the series was asymptotic at $x\to\infty$.
The case $\alpha=1$ was considered separately and it was demonstrated that in
that case the series converges to the generalized Cauchy distribution. An
estimate for the threshold coordinate $x_{\varepsilon}^{N}$ was obtained which
determined the area of applicability of the obtained expansion. It was shown
that in the domain $|x|\geqslant x_{\varepsilon}^{N}$ this power series could
be used to calculate the distribution function, which completely solved the
problem of calculating the distribution function at large $x$.

##  1 Introduction

The main method for calculating the probability density and the distribution
function of stable laws is the use of integral representations of these
quantities. The reason for this situation is the impossibility of obtaining
expressions for these quantities in elementary functions in the general case.
The exception comprises only five cases: the Levy distribution
($\alpha=1/2,\theta=1$), the symmetric Levy distribution
($\alpha=1/2,\theta=-1$), the Cauchy distribution ($\alpha=1,\theta=0$),
Gaussian distribution ($\alpha=2,\theta=0$) and generalized Cauchy
distribution ($\alpha=1,-1\leqslant\theta\leqslant 1$).

When performing the inverse Fourier transform of the characteristic function,
it is possible to obtain two types of integral representations. The first type
includes representations expressing the probability density and the
distribution function in terms of an improper integral of the oscillating
function. The works [1, 2] are devoted to obtaining and studying such
representations. The second type includes integral representations expressing
the probability density and distribution function in terms of a definite
integral of a monotone function. The works[3, 4, 5] are devoted to obtaining
and studying integral representations in the parameterization “B”, the works
[6, 7], are devoted to obtaining and studying integral representations in the
parameterization “M” and the paper [8] to devoted to the parameterization “C”.
Here, to determine various parametrizations of the characteristic function of
the stable law, the notation was used which was introduced in the book [4].
Further in the text, we will continue adhering to these notations.

Integral representations of the second type are most widely used due to the
convenience of their use. The method of the inverse Fourier transform, which
leads to this type of integral representations, is called the stationary phase
method. The convenience of using such representations lies in the fact that
the integrand is a monotonic function and in a wide range of coordinates and
parameters there are no difficulties in calculating the definite integral of
such a function. The integral representations for the parameterization “M” of
the characteristic function gained in popularity. These integral
representations served as the basis for the development of several software
products [9, 10, 11, 12, 13].

Both the first and the second type of the integral representation of a stable
law have their disadvantages. The main difficulty in using the integral
representation of the first type is the oscillating integrand. In some cases,
numerical methods are cannot calculate the integral of such a function. In
particular, the work [1] indicates the following problems for the integral
representation in the parameterization “M”: 1) in the case $\alpha<0.75$ the
integration domain becomes very large, which leads to difficulties in
numerical integration; the integration domain becomes very large which leads
to difficulties in numerical integration; 2) if $\beta\neq 0$ and
$0<|\alpha-1|<0.001$ there are calculation problems in calculating the term
with $(\tan(\pi\alpha/2)(t-t^{\alpha})$; 3) when $x$ is very large, the
integrand oscillates very quickly. In the paper [2] the authors propose to
modernize the standard quadrature numerical integration algorithm to adapt the
calculation of integrals of an oscillating function. This gave an opportunity
to reduce the lower limit of the parameter $\alpha$ from 0.75 to 0.5. To
calculate the probability density at large $x$ it is proposed to use the
expansion of the probability density in a power series. However, the paper
points out that the proposed scheme is not applicable for symmetric
distributions in the case of $\alpha<0.5$ and for asymmetric distributions in
the cases of $\alpha<0.5$ and $0.9<\alpha<1.1$.

The second type of integral representations also has some features that lead
to difficulties in numerical integration. The cause of the calculation
difficulties is the behavior of the integrand at very small and very large
values of the coordinate $x$. In the case of the integral representation for
the probability density in these two cases, the integrand has the form of a
very narrow peak. As a result, numerical integration algorithms cannot
determine this peak and give an incorrect integration result. This behavior of
the integrand is pointed out in the articles [6, 14, 10, 2]. To eliminate this
problem various numerical algorithms are proposed to use in the papers [6, 10,
11] However, all these algorithms increase the accuracy of calculations, but
do not eliminate the problem completely.

To solve this problem, it is expedient to use the approaches not possessing
any specific features in these areas that can lead to calculation
difficulties. The most suitable idea is to use expansions of the probability
density and distribution function in power series at $x\to 0$ and
$x\to\infty$. The articles [15, 16] show that the use of expansions of stable
laws in power series at $x\to 0$ and $x\to\infty$ makes it possible to solve
the problem of calculating stable laws completely at very small and very large
$x$. However, in these articles the problem of calculating the distribution
function of a strictly stable law in the case of $x\to\infty$ was left out of
consideration. Therefore, the main purpose of this paper is to fill this gap.

This paper considers the problem of calculating the distribution function in
the case of $x\to\infty$ with the characteristic function

$\hat{g}(t,\alpha,\theta,\lambda)=\exp\left\\{-\lambda|t|^{\alpha}\exp\\{-i\tfrac{\pi}{2}\alpha\theta\sign
t\\}\right\\},\quad t\in\mathbf{R},$ (1)

where $\alpha\in(0,2]$, $|\theta|\leqslant\min(1,2/\alpha-1)$, $\lambda>0$.
According to the terminology introduced in the book [4], this characteristic
function corresponds to the parameterization “C”. In the paper [8] the inverse
Fourier transform of this characteristic function was performed using the
stationary phase method, and integral representations for the probability
density and distribution function were obtained (see Appendix A). As one can
see the formula (32) expresses the distribution function in terms of a
definite integral, and belongs to the second type of integral representations.
In the general case, the integrand in (33) is a monotonic function varying
from 0 to 1 as the integration variable changes from the lower limit point
$\varphi=-\pi\theta/2$ to the upper limit point $\varphi=\pi/2$. In the case
$0<\alpha<1$ it is a decreasing function, and in the case of
$1<\alpha\leqslant 2$ it is an increasing function. However, at very small and
very large values of $x$ the change in the function from 0 to 1 occurs so fast
that numerical integration algorithms cannot recognize it. As a result, this
leads to an incorrect calculation of the integral and points to the fact that
in this range of coordinates it is no longer possible to use the integral
representation (32) to calculate the distribution function. In this paper, to
calculate the distribution function in the specified range of coordinates, it
is proposed to use the expansion of the distribution function in a power
series at $x\to\infty$. To do this, such an expansion of the distribution
function will be obtained and the conditions for the applicability of this
expansion will be determined.

##  2 Representation of the distribution function as a power series

We will obtain the expansion of the distribution function in a series at
$x\to\infty$ for a strictly stable law with a characteristic function (1).
Without loss of generality, we will assume that the scale parameter is
$\lambda=1$. It is generally accepted to call strictly stable laws with the
scale parameter $\lambda=1$ standard strictly stable laws and shorthand
notations are used for them. The characteristic function is usually denoted as
$\hat{g}(t,\alpha,\theta,1)\equiv\hat{g}(t,\alpha,\theta)$, the probability
density - $g(x,\alpha,\theta,1)\equiv g(x,\alpha,\theta)$, the distribution
function - $G(x,\alpha,\theta,1)\equiv G(x,\alpha,\theta)$, a strictly stable
random quanitity $Y(\alpha,\theta,1)\equiv Y(\alpha,\theta)$. Further in the
text we will use this notation. It should be noted that to transform a
standard strictly stable law into a strictly stable law with an arbitrary
$\lambda$ one can use remark 5 and remark 7 from the paper [8], (see also[4,
17]).

We also need the inversion property, which for a standard strictly stable law
with the characteristic function (1) has the form

For any admissible parameters $(\alpha,\theta)$

$Y(\alpha,-\theta)\stackrel{{\scriptstyle d}}{{=}}-Y(\alpha,\theta).$

The proof of this property was given in the paper [8] (see also [4, 5]). In
the terms of the distribution function $G(x,\alpha,\theta)$ this property
takes the form

$G(-x,\alpha,\theta)=1-G(x,\alpha,-\theta).$ (2)

The convenience of this property lies in the fact that, when studying the
distribution function, it gives us an opportunity to confine ourselves to
considering only the case $x\geqslant 0$. Expressions for the case $x<0$ are
obtained using this formula.

To solve the stated problem, we need to expand the probability density into a
series at $x\to\infty$. A similar expansion was obtained in the article [16],
where the following theorem was proved.

######  Theorem 2.1.

In the case $x\to\pm\infty$ for any admissible set of parameters
$(\alpha,\theta)$ except for the values $\theta=\pm 1$ for the probability
density $g(x,\alpha,\theta)$ the representation in the form of a power series
is valid

$g(x,\alpha,\theta)=g_{N}^{\infty}(|x|,\alpha,\theta^{*})+R_{N}^{\infty}(|x|,\alpha,\theta^{*}),$

where $\theta^{*}=\theta\sign(x)$ and

$g_{N}^{\infty}(x,\alpha,\theta)=\frac{1}{\pi}\sum_{n=0}^{N-1}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n+1)\sin\left(\tfrac{\pi}{2}\alpha n(1+\theta)\right)x^{-\alpha n-1},\quad
x>0,$ (3) $|R_{N}^{\infty}(x,\alpha,\theta)|\leqslant\frac{x^{-\alpha
N-1}}{\pi N!}\left(\Gamma(\alpha
N+1)+x^{-\alpha}\Gamma(\alpha(N+1)+1)\right),\quad x>0.$

Using this theorem, one can obtain an expansion of the distribution function
at $x\to\infty$. As a result, the following theorem turns out to be true.

######  Theorem 2.2.

For any admissible values of parameters $(\alpha,\theta)$ except for the
values $\theta=\pm 1$ at $x\to\pm\infty$ for the distribution function
$G(x,\alpha,\theta)$ the representation in the form of a power series is valid

$G(x,\alpha,\theta)=\tfrac{1}{2}(1+\sign(x))-\sign(x)\left(G_{N}^{\infty}(|x|,\alpha,\theta^{*})+\mathcal{R}_{N}^{\infty}(|x|,\alpha,\theta^{*})\right),$
(4)

where $\theta^{*}=\theta\sign(x)$,

$\displaystyle G_{N}^{\infty}(x,\alpha,\theta)$
$\displaystyle=\frac{1}{\pi}\sum_{n=1}^{N-1}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin\left(\tfrac{\pi}{2}\alpha n(1+\theta)\right)x^{-\alpha n},\quad x>0,$
(5) $\displaystyle|\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)|$
$\displaystyle\leqslant\frac{x^{-\alpha N}}{\pi N!}\left(\Gamma(\alpha
N)+x^{-\alpha}\Gamma(\alpha(N+1))\right),\quad x>0.$ (6)

######  Proof 2.3.

To prove it we will use theorem 2.1. Without loss of generality, we will
consider the case $x>0$. The case $x<0$ can be obtained using the inversion
property (2). By definition, at $x>0$ the distribution function has the form

$G^{(+)}(x,\alpha,\theta)=1-\int_{x}^{\infty}g(\xi,\alpha,\theta)d\xi,\quad
x>0,$ (7)

where $g(x,\alpha,\theta)$ is the probability density of a strictly stable law
and the superscript “$(+)$” shows that this expression determines the
distribution function on the positive part of the semiaxis.

Thus, the expansion of the distribution function in a series is determined by
the expansion of the probability density in a series. It is known that the
expansion of any function in a Taylor series consists of the $N$-th partial
sum and the remainder term. Consequently, the expansion of the probability
density $g(x,\alpha,\theta)$ can be written in the form

$g(x,\alpha,\theta)=g_{N}^{\infty}(x,\alpha,\theta)+R_{N}^{\infty}(x,\alpha,\theta),\quad
x>0,$ (8)

where $g_{N}^{\infty}(x,\alpha,\theta)$ is the $N$-th partial sum and
$R_{N}^{\infty}(x,\alpha,\theta)$ is the remainder term of the series. In the
case $x\to\infty$ the first summand is determined by the expression (3) and
for the remainder term, we use the expression obtained in the article [16]

$R_{N}^{\infty}(x,\alpha,\theta)=\frac{1}{\pi x}\Re
ie^{-i\frac{\pi}{2}\theta}\int_{0}^{\infty}\exp\left\\{-\tau
e^{-i\frac{\pi}{2}\theta}\right\\}R_{N}\left(-\left(\tfrac{i\tau}{x}\right)^{\alpha}\right)d\tau,\quad
x>0,$ (9)

where $R_{N}(y)=\frac{y^{N}}{N!}e^{y\zeta}$, $(0<\zeta<1)$ is the remainder
term in the Lagrange form.

It should be noted that in the case of $x\to\infty$ the expression (8), as
well as theorem 2.1 and the expression (9) are valid when the condition
$\tau/x\to 0$ is satisfied, where $\tau$ is the integration variable in the
inverse Fourier transform formula. In particular, the integration variable in
the formula (9). See the paper [16] for detail. Hence, here, and further in
the text, we will assume everywhere that $\tau/x\to 0$.

To obtain the expansion of the distribution function in a power series at
$x\to\infty$ we will substitute the expression (8) in the expression (7). As a
result, we get

$G^{(+)}(x,\alpha,\theta)=1-\int_{x}^{\infty}g_{N}^{\infty}(\xi,\alpha,\theta)d\xi-\int_{x}^{\infty}R_{N}^{\infty}(x,\alpha,\theta)d\xi,\quad
x>0,$ (10)

where $g_{N}^{\infty}(x,\alpha,\theta)$ has the form (3), and
$R(x,\alpha,\theta)$ is determined by the expression (9).

Interchanging the order of integration and summation in the second summand, we
obtain

$G_{N}^{\infty}(x,\alpha,\theta)=\int_{x}^{\infty}g_{N}^{\infty}(\xi,\alpha,\theta)d\xi=\frac{1}{\pi}\sum_{n=0}^{N-1}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n+1)\sin\left(\tfrac{\pi}{2}\alpha
n(1+\theta)\right)\int_{x}^{\infty}\xi^{-\alpha n-1}d\xi=\\\
=\frac{1}{\pi}\sum_{n=0}^{N-1}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin\left(\tfrac{\pi}{2}\alpha n(1+\theta)\right)x^{-\alpha n},\quad x>0.$

We should pay attention that at $n=0$ the correspondent summand in the sum is
equal to zero. Therefore, the summation can be started with $n=1$. As a
result, we come to the expression (5).

Now we will obtain an expression for the remainder term. Using the expression
(9) and changing the order of integration in some places, for the third
summand in (10) we obtain

$\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)=\int_{x}^{\infty}R_{N}^{\infty}(\xi,\alpha,\theta)d\xi=\frac{1}{\pi}\Re
ie^{-\frac{\pi}{2}\theta}\int_{x}^{\infty}\frac{d\xi}{\xi}\int_{0}^{\infty}\exp\left\\{-\tau
e^{-\frac{\pi}{2}\theta}\right\\}R_{N}\left(-\left(\frac{i\tau}{\xi}\right)^{\alpha}\right)d\tau=\\\
=\frac{1}{\pi}\Re
ie^{-\frac{\pi}{2}\theta}\int_{x}^{\infty}\frac{d\xi}{\xi}\int_{0}^{\infty}\exp\left\\{-\tau
e^{-i\frac{\pi}{2}\theta}\right\\}\frac{1}{N!}\left(-\left(\frac{i\tau}{\xi}\right)^{\alpha}\right)^{N}\exp\left\\{-\left(\frac{i\tau}{\xi}\right)^{\alpha}\zeta\right\\}d\xi=\\\
=\frac{1}{\pi N!}\Re
ie^{-\frac{\pi}{2}\theta}\int_{0}^{\infty}\exp\left\\{-\tau
e^{-i\frac{\pi}{2}\theta}\right\\}(-(i\tau)^{\alpha})^{N}d\tau\int_{x}^{\infty}\xi^{-\alpha
N-1}\exp\left\\{-\left(\frac{i\tau}{\xi}\right)^{\alpha}\zeta\right\\}d\xi$
(11)

Unfortunately, we cannot calculate the integral since the exact value of the
variable $\zeta$ is unknown. It is only known that this variable takes values
from the interval $0<\zeta<1$. Nevertheless, it is possible to estimate the
value of this integral. We will consider
$|\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)|$. Taking into account that the
case $\tau/x\to 0$ is being considered, we can expand the multiplier
$\exp\left\\{-\left(\tfrac{i\tau}{\xi}\right)^{\alpha}\zeta\right\\}$ in a
Taylor series and leave only the summands of the first order of smallness. We
have

$\exp\left\\{-\left(\frac{i\tau}{x}\right)^{\alpha}\zeta\right\\}=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!}\left(\left(\frac{i\tau}{x}\right)^{\alpha}\zeta\right)^{k}\approx
1-\zeta\left(\frac{i\tau}{x}\right)^{\alpha}$

To calculate the obtained integral we will also need the following formula
given in [18] (see §1.5, the formula (31))

$\int_{0}^{\infty}t^{\gamma-1}e^{-ct\cos\beta-
ict\sin\beta}dt=\Gamma(\gamma)c^{-\gamma}e^{-i\gamma\beta},\
-\frac{\pi}{2}<\beta<\frac{\pi}{2},\ \Re\gamma>0\ \text{or}\
\beta=\pm\frac{\pi}{2},\ 0<\Re\gamma<1.$

If we use the Euler formula $\cos\beta+i\sin\beta=e^{i\beta}$, then this
integral can be represented in the form

$\int_{0}^{\infty}t^{\gamma-1}e^{-ct\exp\\{i\beta\\}}dt=\Gamma(\gamma)c^{-\gamma}e^{-i\gamma\beta},\quad-\frac{\pi}{2}<\beta<\frac{\pi}{2},\
\Re\gamma>0\ \text{or}\ \beta=\pm\frac{\pi}{2},\ 0<\Re\gamma<1.$ (12)

Taking into account that $x>0$ for (11) the following estimates turn out to be
valid

$|\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)|=\frac{1}{\pi N!}\left|\Re
ie^{-\frac{\pi}{2}\theta}\int_{0}^{\infty}\exp\left\\{-\tau
e^{-i\frac{\pi}{2}\theta}\right\\}(-(i\tau)^{\alpha})^{N}d\tau\int_{x}^{\infty}\xi^{-\alpha
N-1}\exp\left\\{-\left(\frac{i\tau}{\xi}\right)^{\alpha}\zeta\right\\}d\xi\right|\\\
\leqslant\frac{1}{\pi N!}\left|\int_{0}^{\infty}\exp\left\\{-\tau
e^{-i\frac{\pi}{2}\theta}\right\\}(-(i\tau)^{\alpha})^{N}d\tau\int_{x}^{\infty}\xi^{-\alpha
N-1}\left(1-(i\tau)^{\alpha}\xi^{-\alpha}\zeta\right)d\xi\right|=\\\
=\frac{1}{\pi N!}\left|\int_{0}^{\infty}\exp\left\\{-\tau
e^{-i\frac{\pi}{2}\theta}\right\\}(-(i\tau)^{\alpha})^{N}\left(\frac{x^{-\alpha
N}}{\alpha N}-(i\tau)^{\alpha}\frac{\zeta
x^{-\alpha(N+1)}}{\alpha(N+1)}\right)d\tau\right|\leqslant\\\
\leqslant\frac{1}{\pi N!}\frac{x^{-\alpha N}}{\alpha
N}\left|\int_{0}^{\infty}\exp\left\\{-\tau
e^{-i\frac{\pi}{2}\theta}\right\\}\tau^{\alpha N}d\tau\right|+\frac{1}{\pi
N!}\frac{\zeta
x^{-\alpha(N+1)}}{\alpha(N+1)}\left|\int_{0}^{\infty}\exp\left\\{-\tau
e^{-i\frac{\pi}{2}\theta}\right\\}\tau^{\alpha(N+1)}d\tau\right|=\\\
=\frac{x^{-\alpha N}}{\pi N!}\left(\frac{\Gamma(\alpha N+1)}{\alpha N}+\zeta
x^{-\alpha}\frac{\Gamma(\alpha(N+1)+1)}{\alpha(N+1)}\right)\leqslant\frac{x^{-\alpha
N}}{\pi N!}\left(\Gamma(\alpha N)+x^{-\alpha}\Gamma(\alpha(N+1))\right).$

Here, when passing in the last equality, the formula (12) was used to
calculate the integrals, and when passing in the last inequality, it was
assumed that $\zeta=1$.

To substantiate the validity of using the formula (12) when calculating the
integrals in this expression, we examine the range of the argument
$-\tfrac{\pi}{2}\theta$. The range of admissible values for the parameter
$\theta$ is determined by the inequality
$|\theta|\leqslant\min(1,2/\alpha-1)$. From this it follows that if
$\alpha\leqslant 1$, then $-1\leqslant\theta\leqslant 1$, and if
$1<\alpha\leqslant 2$, then $-(2/\alpha-1)\leqslant\theta\leqslant
2/\alpha-1$. Thus, for any $0<\alpha\leqslant 2$ we obtain
$-\tfrac{\pi}{2}\leqslant-\tfrac{\pi}{2}\theta\leqslant\tfrac{\pi}{2}$. The
extreme values of this interval $\pm\tfrac{\pi}{2}$ are attained at the values
$\alpha\leqslant 1$ and $\theta=\mp 1$. Now we will compare the integral in
(12) with the integrals in the expression during the passage in the last
equality. We see that the integral (12) coincides with these integrals except
for the case $-\tfrac{\pi}{2}\theta=\pm\tfrac{\pi}{2}$. These two points are
out of the range of admissible values for the argument $\beta$ in the formula
(12). Therefore, they should be excluded from consideration.

Now getting back to (10), we obtain

$G^{(+)}(x,\alpha,\theta)=1-G_{N}^{\infty}(x,\alpha,\theta)-\mathcal{R}_{N}^{\infty}(x,\alpha,\theta),\quad
x>0,$ (13)

where for $\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)$ the estimate is valid

$\left|\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)\right|\leqslant\frac{x^{-\alpha
N}}{\pi N!}\left(\Gamma(\alpha N)+x^{-\alpha}\Gamma(\alpha(N+1))\right),\quad
x>0.$

Since the case $x>0$ has been considered so far, these expressions are valid
at $x>0$. To obtain the expansion of the distribution function at $x<0$ we
will use the inversion property. Using in the formula (2) the expression (13)
we obtain

$G^{(-)}(-x,\alpha,\theta)=G_{N}^{\infty}(x,\alpha,-\theta)+\mathcal{R}_{N}^{\infty}(x,\alpha,-\theta),\quad
x>0.$

If we now introduce the notation $\theta^{*}=\theta\sign(x)$ and take the
coordinate $x$ in absolute value, then we can combine the formulas for
$G^{(+)}(x,\alpha,\theta)$ and $G^{(-)}(x,\alpha,\theta)$ into one formula. As
a result, we obtain the expression (4). Thus, the theorem is proved.

$\Box$

The proved theorem determines the expansion of the distribution function of a
strictly stable law with characteristic function (1) into a power series at
$x\to\infty$. Now we examine the issue of the convergence of the obtained
expansion. Since this expansion was obtained by integrating the expansion for
the probability density, taking into account the results of Corollary 1,
proved in the paper [16], one can state that this series converges in the case
of $\alpha<1$ for all $x$, in the case $\alpha=1$, only for $|x|>1$, and in
the case $\alpha>1$ the series is asymptotic one at $x\to\infty$. A more
precise formulation is given by the following corollary.

######  Corollary 2.4.

In the case $\alpha<1$ the series (5) converges for any $x$ at $N\to\infty$.
In this case for the distribution function $G(x,\alpha,\theta)$ for any
admissible $\theta$ the representation is valid in the form of an infinite
series

$G(x,\alpha,\theta)=\tfrac{1}{2}(1+\sign(x))-\frac{\sign(x)}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin\left(\tfrac{\pi}{2}\alpha n(1+\theta^{*})\right)|x|^{-\alpha n}.$ (14)

In the case $\alpha=1$ the series (5) converges for the values $|x|>1$ at
$N\to\infty$. In this case the representation is valid in the form of an
infinite series for the distribution function $G(x,1,\theta)$ for any
admissible $\theta$

$G(x,1,\theta)=\tfrac{1}{2}(1+\sign(x))-\frac{\sign(x)}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n}\sin\left(\tfrac{\pi}{2}n(1+\theta^{*})\right)|x|^{-n},\quad|x|>1.$
(15)

In the case $\alpha>1$ the series (5) diverges for any $x$ at $N\to\infty$. In
this case the asymptotic expansion is valid for the distribution function
$G(x,\alpha,\theta)$ for any admissible $\theta$

$G(x,\alpha,\theta)\sim\tfrac{1}{2}(1+\sign(x))-\frac{\sign(x)}{\pi}\sum_{n=1}^{N-1}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin\left(\tfrac{\pi}{2}\alpha n(1+\theta^{*})\right)|x|^{-\alpha n},\quad
x\to\pm\infty.$ (16)

Here, everywhere $\theta^{*}=\theta\sign(x)$.

######  Proof 2.5.

Without loss of generality, we first consider the case $x>0$. The expansion
for the case $x<0$ will be obtained using the inversion property (2). It was
previously obtained that at positive $x$ the representation (13) is valid.
From this expression and also from (6) it follows that

$|G^{(+)}(x,\alpha,\theta)-1+G_{N}^{\infty}(x,\alpha,\theta)|\leqslant\frac{x^{-\alpha
N}}{\pi N!}\left(\Gamma(\alpha N)+x^{-\alpha}\Gamma(\alpha(N+1))\right),\quad
x>0.$ (17)

We examine the convergence of the series (5). Since this series is sign-
alternating, the inequalities are valid

$G_{N}^{\infty}(x,\alpha,\theta)\leqslant|G_{N}^{\infty}(x,\alpha,\theta)|\leqslant\frac{1}{\pi}\sum_{n=1}^{N-1}\left|\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin\left(\tfrac{\pi}{2}\alpha n(1+\theta)\right)x^{-\alpha n}\right|\\\
\leqslant\frac{1}{\pi}\sum_{n=1}^{N-1}\frac{\Gamma(\alpha
n)}{\Gamma(n+1)}x^{-\alpha n},\quad x>0.$

We will make use of the Cauchy criterion in the limiting form and the Stirling
formula

$\Gamma(z)\sim e^{-z}z^{z-\frac{1}{2}}\sqrt{2\pi},\quad z\to\infty,|\arg
z|<\pi.$ (18)

As a result, we obtain

$\lim_{n\to\infty}\left(\frac{1}{\pi}\frac{\Gamma(\alpha
n)}{\Gamma(n+1)}x^{-\alpha
n}\right)^{1/n}=\lim_{n\to\infty}\left(\frac{e^{-\alpha n}(\alpha n)^{\alpha
n-1/2}\sqrt{2\pi}x^{-\alpha n}}{\pi
e^{-n-1}(n+1)^{n+1-1/2}\sqrt{2\pi}}\right)^{1/n}=\lim_{n\to\infty}\frac{e^{-\alpha}(\alpha
n)^{\alpha-\frac{1}{2n}}x^{-\alpha}}{\pi^{\frac{1}{n}}e^{-1-\frac{1}{n}}(n+1)^{1+\frac{1}{2n}}}\\\
=e^{1-\alpha}\alpha^{\alpha}x^{-\alpha}\lim_{n\to\infty}n^{\alpha}(n+1)^{-1}=e^{1-\alpha}\alpha^{\alpha}x^{-\alpha}\lim_{n\to\infty}n^{\alpha-1}=\begin{cases}0,&\mbox{if}\
\alpha<1,\\\ x^{-1},&\mbox{if}\ \alpha=1,\\\ \infty,&\mbox{if}\
\alpha>1.\end{cases}$ (19)

From this it is clear that in the case $\alpha<1$ the series (5) is convergent
for any $x$, in the case $\alpha=1$ this series converges at $x>1$ and
diverges at $x\leqslant 1$. In the case $\alpha>1$ this series diverges for
any $x$.

We examine the behavior of the remainder term (6) in the case $N\to\infty$.
Using the Stirling’s formula (18) and taking into account that $N+1\approx N$
at $N\to\infty$ we obtain

$\lim_{N\to\infty}\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)\leqslant\lim_{N\to\infty}\frac{x^{-\alpha
N}}{\pi N!}\left(\Gamma(\alpha N)+x^{-\alpha}\Gamma(\alpha(N+1))\right)\\\
=\frac{1}{\pi}\lim_{N\to\infty}\frac{x^{-\alpha N}\Gamma(\alpha
N)+x^{-\alpha(N+1)}\Gamma(\alpha(N+1))}{\Gamma(N+1)}\approx\frac{2}{\pi}\lim_{N\to\infty}x^{-\alpha
N}\frac{\Gamma(\alpha N)}{\Gamma(N)}\\\
=\frac{2}{\pi}\lim_{N\to\infty}x^{-\alpha N}\frac{e^{-\alpha N}(\alpha
N)^{\alpha
N-1/2}\sqrt{2\pi}}{e^{-N}N^{N-1/2}\sqrt{2\pi}}=\frac{2}{\pi}\lim_{N\to\infty}x^{-\alpha
N}e^{N(1-\alpha)}\alpha^{\alpha N-1/2}N^{N(\alpha-1)}\\\
=\frac{2}{\pi\sqrt{\alpha}}\lim_{N\to\infty}\exp\left\\{N(1-\alpha)(1-\ln
N)+\alpha N(\ln\alpha-\ln x)\right\\}=\begin{cases}0,&\mbox{if }\alpha<1\\\
\infty,&\mbox{if }\alpha=1,x\leqslant 1,\\\ 0,&\mbox{if }\alpha=1,x>1,\\\
\infty,&\mbox{if }\alpha>1.\end{cases}$ (20)

We will consider the case $\alpha<1$. Generalizing the results obtained, we
see that in this case the series (5) is convergent, and the limit of the
remainder term $\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)$ is equal to zero.
From this it follows that the right part (17) is an element of an
infinitesimal sequence. In turn, this means that for any fixed $x$ the
sequence $1-G_{N}^{\infty}(x,\alpha,\theta)$ converges to the distribution
function $G^{(+)}(x,\alpha,\theta)$ at $N\to\infty$. Consequently, in the case
$\alpha<1$, for the distribution function $G^{(+)}(x,\alpha,\theta)$ the
representation is valid in the form of an infinite series

$G^{(+)}(x,\alpha,\theta)=1-\frac{1}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin(\tfrac{\pi}{2}\alpha n(1+\theta))x^{-\alpha n},\quad x>0$

Using the inversion property (2) for negative $x$ we obtain

$G^{(-)}(x,\alpha,\theta)=\frac{1}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin(\tfrac{\pi}{2}\alpha n(1-\theta))(-x)^{-\alpha n},\quad x<0.$

If now we introduce the parameter $\theta^{*}=\theta\sign(x)$ and take the
coordinate $x$ in absolute value, then we can combine the last two formulas.
As a result, we get the formula (14). This proves the first item of the
corollary.

Now we consider the case $\alpha=1$ and $x>1$. As it follows from the
expression (19), in this case the series (5) is convergent. It also follows
from the expression (20) that
$\lim_{N\to\infty}\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)=0$. Therefore, the
right part of the expression (17) is an element of an infinitesimal sequence.
In turn, this means that for any fixed $x>1$ the sequence
$1-G_{N}^{\infty}(x,1,\theta)$ converges to the distribution function
$G^{(+)}(x,1,\theta)$. Therefore, the representation in the form of an
infinite series is valid for the distribution function $G^{(+)}(x,1,\theta)$

$G^{(+)}(x,1,\theta)=1-\frac{1}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n}\sin\left(\tfrac{\pi}{2}n(1+\theta)\right)x^{-n},\quad
x>1.$ (21)

Using the formula (2) for negative $x$ we obtain the formula

$G^{(-)}(x,1,\theta)=\frac{1}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n}\sin\left(\tfrac{\pi}{2}n(1-\theta)\right)(-x)^{-n},\quad
x<-1.$

Using now the notation $\theta^{*}=\theta\sign(x)$ and taking the coordinate
$x$ in absolute value we can combine the last two formulas in one expression.
As a result, we obtain the expression (15). This proves the second item of the
corollary.

Now we consider the case $\alpha>1$. From the expression (19) it follows that
in this case the series (5) is divergent at $N\to\infty$. However, from the
expression (6) it is clear that at some fixed $N$ the estimate is valid

$\mathcal{R}_{N}^{\infty}(x,\alpha,\theta)=O(x^{-\alpha N}),\quad x\to\infty.$

Thus, for each fixed $N$ at $x>0$ from the expression (13) we obtain

$G^{(+)}(x,\alpha,\theta)=1-\frac{1}{\pi}\sum_{n=1}^{N-1}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin(\tfrac{\pi}{2}\alpha n(1+\theta))x^{-\alpha n}+O\left(x^{-\alpha
N}\right),\quad x\to\infty.$

Using the formula (2) we obtain the representation for negative $x$

$G^{(-)}(x,\alpha,\theta)=\frac{1}{\pi}\sum_{n=1}^{N-1}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin(\tfrac{\pi}{2}\alpha n(1-\theta))(-x)^{-\alpha n}+O\left((-x)^{-\alpha
N}\right),\quad x\to-\infty.$

Introducing now the notation $\theta^{*}=\theta\sign(x)$ and taking the
coordinate $x$ in absolute value we can combine the last two expressions into
one. As a result, we obtain

$G(x,\alpha,\theta)=\tfrac{1}{2}(1+\sign(x))-\frac{\sign(x)}{\pi}\sum_{n=1}^{N-1}\frac{(-1)^{n+1}}{n!}\Gamma(\alpha
n)\sin(\tfrac{\pi}{2}\alpha n(1+\theta^{*}))|x|^{-\alpha
n}+O\left(|x|^{-\alpha N}\right),\quad x\to\pm\infty.$

Thus, we have obtained the definition of an asymptotic series. Therefore, this
expression can be written as (16). This proves the third item of the
corollary. It should be noted that we do not consider the case $\alpha=1$ and
$|x|<1$ since in this case the series (5) diverges. Thus, the corollary is
completely proved.

$\Box$

It should be noted that the case $\alpha=1$ is related to one of those few
cases when both the probability density and the distribution function are
expressed in terms of elementary functions. In this case, the distribution
function is expressed by the formula (34). The derivation of this formula, as
well as the proof of corollary A.2 (see Appendix A) can be found in the
paper[8]. In the article [16] it was shown that the expansion of the
probability density in a series in the case $\alpha=1$ at $x\to\infty$
converges to the density
$g(x,1,\theta)=\frac{\cos(\pi\theta/2)}{\pi(x^{2}-2x\sin(\pi\theta/2)+1)}$ at
$N\to\infty$. Similarly, for the distribution function one can show that the
expansion (15) converges to the distribution function (34) at $|x|>1$. We
formulate this result in the form of a remark.

######  Remark 2.6.

In the case $\alpha=1$ for any $-1<\theta<1$ in the domain $|x|>1$ the series
(15) converges to the distribution function (34).

######  Proof 2.7.

To prove this, we will consider the distribution function (34) and show that
the expansion of this distribution function in a Taylor series at $x\to\infty$
has the form (15). Using the reduction formulas
$\cos\left(\tfrac{\pi}{2}\theta\right)=\sin\left(\tfrac{\pi}{2}+\tfrac{\pi}{2}\theta\right)$,
$\sin\left(\tfrac{\pi}{2}\theta\right)=-\cos\left(\tfrac{\pi}{2}+\tfrac{\pi}{2}\theta\right)$,
we will write down the distribution function (34) in the form

$G(x,1,\theta)=\frac{1}{2}+\frac{1}{\pi}\arctan\left(\frac{x+\cos\left(\tfrac{\pi}{2}(1+\theta)\right)}{\sin\left(\tfrac{\pi}{2}(1+\theta)\right)}\right)$
(22)

Further, since we need to obtain the expansion of the distribution function in
a Taylor series at $x\to\infty$, then in this expression we will substitute
the variable $x=1/y$ in this expression and start considering the case
$x\geqslant 0$. The relationship obtained in this way we will denote as
$G^{(+)}(y,1,\theta)$. As a result, we obtain

$G^{(+)}(y,1,\theta)=\frac{1}{2}+\frac{1}{\pi}\arctan\left(\frac{\tfrac{1}{y}+\cos\left(\tfrac{\pi}{2}(1+\theta)\right)}{\sin\left(\tfrac{\pi}{2}(1+\theta)\right)}\right),\quad
y\geqslant 0.$ (23)

Hence it is clear that the behavior of the function (23) at $y\to 0$
corresponds to the behavior of the function (22) at $x\to\infty$. Therefore,
expanding the function (23) in a Taylor series in the vicinity of the point
$y=0$ and getting back to the variable $x$, we obtain the expansion of the
distribution function (22) into a power series at $x\to\infty$.

We will take into account that the function $\arctan(x)$ is infinitely
differentiable. Consequently, the expansion of the function (23) in a Taylor
series in the vicinity of the point $y=0$ has the form

$G^{(+)}(y,1,\theta)=G^{(+)}(0,1,\theta)+\sum_{n=1}^{\infty}\frac{1}{n!}\left.\frac{d^{n}G^{(+)}(y,1,\theta)}{dy^{n}}\right|_{y=0}y^{n}.$
(24)

At the beginning, we will calculate the first derivative of the function (23).
We obtain

$\frac{dG^{(+)}(y,1,\theta)}{dy}=-\frac{\sin\left(\tfrac{\pi}{2}(1+\theta)\right)}{\pi\left(y^{2}+2y\cos\left(\tfrac{\pi}{2}(1+\theta)\right)+1\right)}.$

We represent this expression in the form

$\frac{dG^{(+)}(y,1,\theta)}{dy}=-\frac{\sin\left(\tfrac{\pi}{2}(1+\theta)\right)}{\pi}f(g(y)),$

where

$f\equiv f(g)=1/g,\quad g\equiv
g(y)=y^{2}+2y\cos\left(\tfrac{\pi}{2}(1+\theta)\right)+1.$ (25)

Thus, for the derivative of the series of $n$ of the function
$G^{(+)}(y,1,\theta)$ we get

$\frac{d^{n}G^{(+)}(y,1,\theta)}{dy^{n}}=\frac{d^{n-1}}{dy^{n-1}}\frac{dG^{(+)}(y,1,\theta)}{dy}=-\frac{\sin\left(\tfrac{\pi}{2}(1+\theta)\right)}{\pi}\frac{d^{n-1}f(g(y))}{dy^{n-1}}\\\
=-\frac{\sin\left(\tfrac{\pi}{2}(1+\theta)\right)}{\pi}\sum_{k=0}^{\left[\tfrac{n-1}{2}\right]}\frac{(-1)^{n-1-k}(n-1)!(n-1-k)!}{k!(n-1-2k)!}\frac{\left(2y+2\cos\left(\tfrac{\pi}{2}(1+\theta)\right)\right)^{n-1-2k}}{\left(y^{2}+2y\cos\left(\tfrac{\pi}{2}(1+\theta)\right)+1\right)^{n-k}},$

where the formula was used

$\frac{d^{n}f(g(y))}{dy^{n}}=\sum_{k=0}^{\left[\tfrac{n}{2}\right]}\frac{(-1)^{n-k}n!(n-k)!}{k!(n-2k)!}\frac{\left(2y+2\cos\left(\tfrac{\pi}{2}(1+\theta)\right)\right)^{n-2k}}{\left(y^{2}+2y\cos\left(\tfrac{\pi}{2}(1+\theta)\right)+1\right)^{n-k+1}},$

which was obtained in the article [16].

At the point $y=0$ this derivative has the value

$\left.\frac{d^{n}G^{(+)}(y,1,\theta)}{dy^{n}}\right|_{y=0}=-\frac{\sin\left(\tfrac{\pi}{2}(1+\theta)\right)}{\pi}(-1)^{n-1}(n-1)!\\\
\times\sum_{k=0}^{\left[\tfrac{n-1}{2}\right]}\frac{(-1)^{k}(n-k-1)!}{k!(n-2k-1)!}\left(2\cos\left(\tfrac{\pi}{2}(1+\theta)\right)\right)^{n-2k-1}=\frac{(n-1)!}{\pi}(-1)^{n-1}\sin\left(\tfrac{\pi}{2}n(1+\theta)\right).$
(26)

Here it was taken into account that $(-1)^{-k}=(-1)^{k}$, and also the formula
was used for $\sin(n\varphi)$ (see, for example, [19])

$\sin(n\varphi)=\sin\varphi\sum_{k=0}^{\left[\frac{n-1}{2}\right]}(-1)^{k}\frac{(n-k-1)!}{k!(n-2k-1)!}(2\cos\varphi)^{n-2k-1}.$

Now we substitute the expression (26) in the formula (24) and take into
account that $G^{(+)}(0,1,\theta)=\frac{1}{2}+\arctan(\infty)=1$. As a result,
the expression (24) becomes

$G^{(+)}(y,1,\theta)=1-\frac{1}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}\sin\left(\tfrac{\pi}{2}(1+\theta)\right)y^{n}.$

Thus, we have obtained the expansion of the function (23) in a Taylor series
in the vicinity of the point $y=0$.

To obtain the expansion of the distribution function (22) at $x\to\infty$, one
must return to the variable $x$ in the last expression. By substituting the
variable $y=1/x$, we obtain

$G^{(+)}(x,1,\theta)=1-\frac{1}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}\sin\left(\tfrac{\pi}{2}(1+\theta)\right)x^{-n}.$

This expression is the expansion of the distribution function (22) at
$x\to\infty$. As one can see, it completely coincides with the expansion
obtained earlier (21). Corollary 2.4 shows that this series converges at
$n\to\infty$ in the domain $x>1$. Hence

$G^{(+)}(x,1,\theta)=1-\frac{1}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}\sin\left(\tfrac{\pi}{2}(1+\theta)\right)x^{-n},\quad
x>1.$

To obtain the expansion of the distribution function for negative $x$ we use
the inversion property and, in particular, the formula (2). As a result, we
get

$G^{(-)}(x,1,\theta)=\frac{1}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n}\sin\left(\tfrac{\pi}{2}n(1-\theta)\right)(-x)^{-n},\quad
x<-1.$

If we now introduce the parameter $\theta^{*}=\theta\sign(x)$ and take the
variable $x$ in absolute value then it is possible to combine the formulas for
$G^{(+)}(x,1,\theta)$ and $G^{(-)}(x,1,\theta)$ into one formula. As a result,
we obtain

$G(x,1,\theta)=\tfrac{1}{2}(1+\sign(x))-\frac{\sign(x)}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n}\sin\left(\tfrac{\pi}{2}n(1+\theta^{*})\right)|x|^{-n},\quad|x|>1.$

This formula coincides completely with the expression (15). Thus, the
expansion of the distribution function (22) coincides exactly with the
expansion (15)which was obtained earlier. Therefore, in the domain $|x|>1$ the
series (15) converges to the distribution function (34). The remark has been
proved.

$\Box$

Theorem 2.2 gives an opportunity to calculate the distribution function for
large values of the coordinate $x$ using a power series (4). As mentioned in
the Introduction, for large values of the coordinate $x$ the use of the
integral representation (32) no longer allows one to calculate the
distribution function correctly, and here it is necessary to apply other
calculation methods. The most obvious and probably the only way to calculate
the distribution function at large values of $x$ is to use the expansions (4).

However, before we use this expansion, it is necessary to obtain a criterion
that makes it possible to determine the coordinate $x$ in which, for a
specified value of $\alpha$ and the number of summands $N$ in the sum (4), the
specified accuracy of calculating the distribution function will be achieved.
Such a criterion can be obtained by evaluating the remainder term (6). From
the formula (4) and (6) it follows that

$\left|G(x,\alpha,\theta)-\tfrac{1}{2}(1+\sign(x))+\sign(x)G_{N}^{\infty}(|x|,\alpha,\theta^{*})\right|\leqslant\frac{|x|^{-\alpha
N}}{\alpha N!}\left(\Gamma(\alpha N)+|x|^{-\alpha}\Gamma(\alpha(N+1))\right).$

If now for the specified $\alpha$ and $N$ we give the value of the absolute
error $\varepsilon$, the value of which should not exceed the true absolute
error of the calculation using the expansion (4), i.e.

$\left|G(x,\alpha,\theta)-\tfrac{1}{2}(1+\sign(x))+\sign(x)G_{N}^{\infty}(|x|,\alpha,\theta^{*})\right|\leqslant\varepsilon,$

then this gives an opportunity to introduce the threshold coordinate
$x_{\varepsilon}^{N}$. The value of the threshold coordinate can be found from
the solution to the equation

$\frac{\left|x_{\varepsilon}^{N}\right|^{-\alpha N}}{\alpha
N!}\left(\Gamma(\alpha
N)+\left|x_{\varepsilon}^{N}\right|^{-\alpha}\Gamma(\alpha(N+1))\right)=\varepsilon.$
(27)

Unfortunately, it is impossible to solve this equation analytically and obtain
an explicit expression for the threshold coordinate $x_{\varepsilon}^{N}$.
Nevertheless, numerical methods without much difficulty help to find a
solution to this equation for specified $\alpha$, $N$ and $\varepsilon$. As a
result, we obtain the condition

$\left|G(x,\alpha,\theta)-\tfrac{1}{2}(1+\sign(x))+\sign(x)G_{N}^{\infty}(|x|,\alpha,\theta^{*})\right|\leqslant\varepsilon,\quad|x|\geqslant
x_{\varepsilon}^{N}.$ (28)

It means that at the values of $|x|\geqslant x_{\varepsilon}^{N}$ the absolute
error in calculating the distribution function using a power series (4) will
not exceed the specified accuracy level $\varepsilon$. Here $\alpha,N$ and
$\varepsilon$ are specified and $x_{\varepsilon}^{N}$ is found from the
solution to the equation (27).

Taking account of all the foregoing, we obtain formulas for calculating the
distribution function

$G_{N}(x,\alpha,\theta)=\tfrac{1}{2}(1+\sign(x))-\sign(x)G_{N}^{\infty}(|x|,\alpha,\theta^{*}),\quad|x|\geqslant
x_{\varepsilon}^{N}.$ (29)

Here $G_{N}^{\infty}(x,\alpha,\theta)$ is determined by the expression (5),
and the threshold coordinate $x_{\varepsilon}^{N}$ is found from the solution
to the equation (27). The accuracy level $\varepsilon$ and the number of
summands $N$ are specified beforehand. In this case, it can be guaranteed that
in the domain $|x|\geqslant x_{\varepsilon}^{N}$ at a specified $N$ the
absolute error in calculating the distribution function using this formula
will not exceed $\varepsilon$.

Figures 1, 2, and 3 show the calculation results of the distribution function
using the formula (29) and the results of calculating the quantity of the
absolute error for the values $\alpha=0.7,1,1.3$. Figures 1a, 2a, and 3a
demonstrate the results of calculating the distribution function. In these
figures, the solid curve corresponds to the exact value of the distribution
function. In the case $\alpha\neq 1$ the integral representation (32) was used
for calculations, and in the case $\alpha=1$ the formula (34) was used. The
dashed-dotted curves in these figures correspond to the results of the
calculation using the expansion (29) for the values $N=3,10,30,60,90$. The
circles in these figures show the position of the threshold coordinate
$x_{\varepsilon}^{N}$ for the chosen values of $N$ and the specified accuracy
level $\varepsilon=10^{-5}$. The value of the threshold coordinate for each
$N$ and selected $\varepsilon$ was found by solving the equation (27).

Figures 1b, 2b, and 3b give the absolute error of calculating the distribution
function using the expansion (29). In these figures, the solid curve is the
exact value of the absolute error
$|G(x,\alpha,\theta)-G_{N}(x,\alpha,\theta)|$. Here to calculate
$G(x,\alpha,\theta)$ the integral representation (32) was used in the case
$\alpha\neq 1$, and the formula (34) was applied in the case $\alpha=1$. To
calculate $G_{N}(x,\alpha,\theta)$ the formula (29) was used. The dashed-
dotted curves are the estimate of the remainder term (6). The dotted line
shows the position of the specified accuracy level $\varepsilon$, and the
circles are the position of the threshold coordinate $x_{\varepsilon}^{N}$ for
each value of $N$.

Figure 1: a) The distribution function $G(x,\alpha,\theta)$ for the parameter
values shown in the figure. The solid curve is the integral representation
(32), the dash-dotted curves are the representation in the form of a power
series (29) for different values of the number of summands $N$ in the sum. The
circles show the position of the threshold coordinate $x_{\varepsilon}^{N}$
for corresponding values of $N$ and the specified accuracy level
$\varepsilon$. (b) Graph of the absolute error of calculating the distribution
function $G(x,\alpha,\theta)$ with the use of a power series (29). The solid
curves are the exact value of the absolute error
$|G(x,\alpha,\theta)-G_{N}(x,\alpha,\theta)|$, the dash-dotted curves are the
residual term estimate (6), the dotted line is the specified accuracy level
$\varepsilon$, the circles show the position of the threshold coordinate
$x_{\varepsilon}^{N}$

Figure 2: a) The distribution function $G(x,\alpha,\theta)$ for the parameter
values shown in the figure. The solid curve is the formula (34), the dash-
dotted curves are the power series representation (29) for different values of
the number of summands $N$ in the sum. The circles show the position of the
threshold coordinate $x_{\varepsilon}^{N}$ for the corresponding values of $N$
and the specified level of accuracy $\varepsilon$. (b) Graph of the absolute
error of calculating the distribution function $G(x,\alpha,\theta)$ with the
use of a power series (29). The solid curves are the exact value of the
absolute error $|G(x,\alpha,\theta)-G_{N}(x,\alpha,\theta)|$, the dash-dotted
curves are the residual term estimate (6), the dotted line is the specified
accuracy level $\varepsilon$, the circles show the position of the threshold
coordinate $x_{\varepsilon}^{N}$

Figure 3: a) The distribution function $G(x,\alpha,\theta)$ for the parameter
values shown in the figure. The solid curve is the integral representation
(32), the dash-dotted curves are the representation in a power series (29) for
different values of the number of summands $N$ in the sum. The circles show
the position of the threshold coordinate $x_{\varepsilon}^{N}$ for the
corresponding values $N$ and the specified accuracy level $\varepsilon$. (b)
Graph of the absolute error of calculating the distribution function
$G(x,\alpha,\theta)$ using a power series (29). The solid curves are the exact
value of the absolute error $|G(x,\alpha,\theta)-G_{N}(x,\alpha,\theta)|$, the
dash-dotted curves are the residual term estimate (6), the dotted line is the
specified accuracy level $\varepsilon$, the circles show the position of the
threshold coordinate $x_{\varepsilon}^{N}$

Figures 1b, 2b, and 3b clearly show that in the domain $x>x_{\varepsilon}^{N}$
the exact value of the absolute error (solid curves) turns out to be less than
the estimate of the remainder term (6) (dashed-dotted curves). It is also
clearly seen that in the domain $x>x_{\varepsilon}^{N}$, both the exact value
of the absolute error and the residual term estimate are less than the value
of the selected accuracy level $\varepsilon$. This means that in the domain
$|x|>x_{\varepsilon}^{N}$ the formula (29) can be used to calculate the
distribution function. At the same time, it can be guaranteed that the
absolute error in calculating the distribution function using this formula
will not exceed the selected accuracy level $\varepsilon$, and in reality it
will be much less than this value.

If we now analyze the behavior of the threshold coordinate
$x_{\varepsilon}^{N}$ from the number of terms $N$, then we can see that for
each of the considered cases $\alpha<1$, $\alpha=1$ and $\alpha>1$ this
behavior differs. It can be seen from fig. 1 that in the case $\alpha=0.7$ as
$N$ increases the value $x_{\varepsilon}^{N}$ decreases. This behavior of the
threshold coordinate is the result of corollary 2.4. Indeed, in the first item
of this corollary it is proved that in the case $\alpha<1$ at $N\to\infty$ the
series (14) converges for any $x$. This means that as the number of terms in
the formula (29) increases, the accuracy of calculating the distribution
function at some fixed point $x$ will increase. In turn, this leads to the
fact that the range of coordinates $x$ under which the condition (28) is
satisfied will expand. Therefore, as $N$ increases the value of the threshold
coordinate $x_{\varepsilon}^{N}$ will decrease.

To calculate the limit value of the coordinate $x_{\varepsilon}^{N}$ at
$N\to\infty$ we will consider the equation (27) and assume that $N\to\infty$.
Taking account that $N+1\approx N$ at $N\to\infty$ this equation takes the
form $2|x_{\varepsilon}^{N}|^{-\alpha N}\Gamma(\alpha
N)=\varepsilon\pi\Gamma(N)$. From here we find

$|x_{\varepsilon}^{N}|=\left(\frac{2\Gamma(\alpha
N)}{\pi\varepsilon\Gamma(N)}\right)^{\frac{1}{\alpha N}}.$

We now find the limit of this expression at $N\to\infty$. Using the Stirling
formula (18), we obtain

$\lim_{N\to\infty}|x_{\varepsilon}^{N}|=\lim_{N\to\infty}\left(\frac{2\Gamma(\alpha
N)}{\pi\varepsilon\Gamma(N)}\right)^{\frac{1}{\alpha
N}}=\lim_{N\to\infty}\left(\frac{2}{\pi\varepsilon}\frac{e^{-\alpha N}(\alpha
N)^{\alpha
N-\frac{1}{2}}\sqrt{2\pi}}{e^{-N}(N)^{N-\frac{1}{2}}\sqrt{2\pi}}\right)^{\frac{1}{\alpha
N}}\\\
=e^{\frac{1}{\alpha}-1}\alpha^{-1}\lim_{N\to\infty}\left(\frac{2}{\pi\varepsilon\alpha^{2}}\right)^{\frac{1}{\alpha
N}}N^{1-\frac{1}{\alpha}}=\begin{cases}0,&\alpha<1\\\ 1,&\alpha=1\\\
\infty,&\alpha>1.\end{cases}$ (30)

Thus, in the case $\alpha<1$ we get
$\lim_{N\to\infty}|x_{\varepsilon}^{N}|=0$. This result is a consequence of
the convergence of the series (5) in the case $\alpha<1$.

A similar behavior of the threshold coordinate is also observed in the case
$\alpha=1$. Figure 2 shows that as the number of $N$ summands in the formula
(29) increases the value of the threshold coordinate decreases. However,
unlike the previous case, it follows from the expression (30) that in this
case $\lim_{N\to\infty}|x_{\varepsilon}^{N}|=1$. Such behavior of the
threshold coordinate is the result of corollary 2.4, where in the second item
it was proved that in the case $\alpha=1$ the series (15) converges at
$N\to\infty$ in the domain $|x|>1$.

In the case $\alpha>1$ the behavior of the threshold coordinate changes.
Figure 3a shows that as the number of summands $N$ in the formula (29)
increases, the value of the threshold coordinate first decreases:
$x_{\varepsilon}^{3}>x_{\varepsilon}^{10}>x_{\varepsilon}^{30}$. However, a
further increase in $N$ leads to an increase in the threshold coordinate:
$x_{\varepsilon}^{30}<x_{\varepsilon}^{60}<x_{\varepsilon}^{90}$. Such
behavior of the threshold is in full accordance with corollary 2.4, where in
the third item it was proved that in the case $\alpha>1$ the series (5)
diverged at $N\to\infty$. The cause of the divergence of this series lies in
the presence of the multiplier $\Gamma(\alpha n)/\Gamma(n+1)$. Hence, it is
clear that at $\alpha>1$ this multiplier is more than 1 and as $n$ increases
the value of this multiplier only increases. Such behavior of this multiplier
is the reason for the divergence of the series (5). In this series there is
also a multiplier $x^{-\alpha n}$. As one can see, as the value of $x$
increases, the value of this multiplier decreases. The competition between
these two factors leads to the observed behavior of the threshold coordinate
$x_{\varepsilon}^{N}$.

Indeed, the threshold coordinate is found as a result of solving the equation
$|G(x,\alpha,\theta)-G_{N}(x,\alpha,\theta)|=\varepsilon$. Therefore, first
the increase in the number of summands $N$ in the sum (29) leads to a decrease
in the coordinate $x_{\varepsilon}^{N}$. This is testified by the fact that
$x_{\varepsilon}^{3}>x_{\varepsilon}^{10}>x_{\varepsilon}^{30}$. However,
further increase in $N$ leads to the fact that the factor $\Gamma(\alpha
n)/\Gamma(n+1)$ starts growing rapidly, and to compensate for this growth, it
is necessary to increase $x$ in the multiplier $x^{-\alpha n}$. This is what
leads to a shift of the threshold coordinate towards larger values of $x$ as
$N$ increases. Fig. 3 demonstrates such behavior from which it is clear that
$x_{\varepsilon}^{30}<x_{\varepsilon}^{60}<x_{\varepsilon}^{90}$. In case, if
$n\to\infty$, then the factor $\Gamma(\alpha n)/\Gamma(n)\to\infty$.
Therefore, to compensate for this growth, it is necessary that $x\to\infty$ in
the multiplier $x^{-\alpha n}$. Thus, the obtained conclusion is in full
accordance with the expression (30), which shows that
$\lim_{N\to\infty}|x_{\varepsilon}^{N}|=\infty$, if $\alpha>1$.

It should be pointed out, if in the considered case ($\alpha>1$) we fix some
arbitrary $N$, then in view of the presence of the factor $x^{-\alpha N}$ with
an increase in the value of $x$ one can achieve any preset calculation
accuracy $\varepsilon$. Consequently, at $\alpha>1$ the formula (29) is
asymptotic at $|x|\to\infty$, which is the result of corollary 2.4.

##  3 Calculation of the distribution function at $x\to\infty$

We return to the question of calculating the distribution function of a
strictly stable law in the case of large values of the coordinate $x$. The
main approach to calculate the distribution function is to use the integral
representation (32). In theory, this integral representation is valid for all
values of the parameters $\alpha,\theta$ (except for the value $\alpha=1$) and
all $x$. However, in practice, it is not always possible to calculate the
integral numerically included in this integral representation. Problems arise
at small and large values of the coordinate $x$. The cause of the difficulties
that arise is the behavior of the integrand in the formula (32).

Figure 4:  The relationship between the integrand of the integral
representation for the distribution function (32) and the value of the
integration variable $\varphi$. The figure shows graphs of the integrand for
the values of the parameters $\alpha=1.1,\theta=0$ and the specified values of
the coordinate $x$

Fig. 4 presents a graph of the integrand of the integral representation (32)
for the parameters $\alpha=1.1,\theta=0$ and the specified values of
coordinates $x$ depending on the integration variable $\varphi$. The variable
$\varphi$ changes within the range from $-\pi\theta/2$ to $\pi/2$. It is clear
from the figure that at very small and large values of $x$ the integrand in
(32) increases very sharply from 0 to the value of 1. In the case $\alpha<1$
the behavior of the integrand will be inverse. In this case, the function is
decreasing, and therefore it will sharply decrease from 1 to 0. With a further
decrease or increase in the value of $x$ the steep increase (in the case
$\alpha>1$) or decrease (in the case $\alpha<1$) in sections will increase. As
a result, at some $x$ numerical integration algorithms cannot recognize the
monotonic nature of the function and begin to produce an incorrect result.

Figure 5:  Distribution function $G(x,\alpha,\theta)$ a strictly stable law.
The figure on the left is the case $\alpha<1$, the figure on the right is the
case $\alpha>1$. The values of the indicator $\alpha$ are given in the figures
for all graphs $\theta=0$. The solid curves are integral representation (32),
dash-dotted curves are the representation in a power series (29)

The most suitable method for calculating the distribution function in the case
of small and large values of $x$ is to use asymptotic expansions. The problem
of calculating the distribution function in the case $x\to 0$ was considered
in the article [15]. In the case $x\to\infty$ it is expedient to use theorem
2.2 and, in particular, the formula (29). Figure 5 shows the results of
calculating the distribution function $G(x,\alpha,\theta)$ using the integral
representation (32) (solid curves) and the formula (29) (dash-dotted curves)
at large values of $x$. The left figure shows the case $\alpha<1$, the right
one shows the case $\alpha>1$. To calculate the integral in the formula (32)
the Gauss-Kronrod algorithm was used. One can see from the figures that at
large values of $x$ the numerical integration algorithm used is incapable of
calculating the integral in (32) and starts giving an incorrect result. It can
also be seen from the figure that the value of the critical coordinate
$x_{\mbox{\scriptsize cr}}$, at which the numerical integration algorithm
starts calculating the integral incorrectly, depends on the value of $\alpha$.
For the value $\alpha=0.5$ the value is $x_{\mbox{\scriptsize cr}}\approx
3.2\cdot 10^{13}$, for the value $\alpha=0.7$ the value is
$x_{\mbox{\scriptsize cr}}\approx 3\cdot 10^{9}$, for the value $\alpha=0.9$
the value is $x_{\mbox{\scriptsize cr}}\approx 10^{7}$, for the value
$\alpha=1.1$ the value is $x_{\mbox{\scriptsize cr}}\approx 4\cdot 10^{5}$,
for the value $\alpha=1.4$ the value is $x_{\mbox{\scriptsize cr}}\approx
3.5\cdot 10^{4}$, for the value $\alpha=1.7$ the value is
$x_{\mbox{\scriptsize cr}}\approx 6\cdot 10^{3}$. It is clear that as the
value $\alpha$ decreases, the value $x_{\mbox{\scriptsize cr}}$ increases.
Thus, at $|x|>x_{\mbox{\scriptsize cr}}$ other methods of calculating the
distribution function should be used.

The use of the formula (29) to calculate the distribution function at
$|x|>x_{\mbox{\scriptsize cr}}$ solves the problem completely. It is clear
that in the domain $x<x_{\mbox{\scriptsize cr}}$ the calculation results using
the integral representation (32) and the formula (29) coincide completely. It
should be noted that such a coincidence will be observed in the domain
$x_{\varepsilon}^{N}\leqslant x\leqslant x_{\mbox{\scriptsize cr}}$. In the
domain $x>x_{\mbox{\scriptsize cr}}$ the numerical integration algorithm no
longer makes it possible to obtain the correct value of the distribution
function using the integral representation (32), whereas the use of the
formula (29) does not lead to any calculation difficulties. It should be noted
that the number of summands $N=30$ was used to calculate the distribution
function with the help of the formula (29). The threshold coordinate
$x_{\varepsilon}^{N}$ for the accuracy level $\varepsilon=10^{-5}$ for the
presented graphs has the following values: at
$\alpha=0.5,x_{\varepsilon}^{30}=0.088$, at
$\alpha=0.7,x_{\varepsilon}^{30}=0.402$, at
$\alpha=0.9,x_{\varepsilon}^{30}=1.000$, at
$\alpha=1.1,x_{\varepsilon}^{30}=1.860$, at
$\alpha=1.4,x_{\varepsilon}^{30}=3.552$ and at
$\alpha=1.7,x_{\varepsilon}^{30}=5.612$. As one can see, the values of the
threshold coordinate for each of the graphs presented in the figure are
significantly less than the range of values of $x$, which is given in the
figures. Consequently, it can be asserted that in the domain $|x|\geqslant
x_{\varepsilon}^{N}$ to calculate the distribution function, one can use the
formula (29). In this case, the absolute error of calculating the distribution
function will not exceed the specified accuracy level $\varepsilon$ and as $x$
increases the absolute error of calculation will only decrease. Thus, the use
of theorem 2.2 and, in particular, the formula (29) solves the problem of
calculating the distribution function completely at $x\to\infty$.

It should be noted that the presented results are related to standard strict-
stable laws, i.e. to laws with scale parameter $\lambda=1$. To transform the
distribution function of a standard strictly stable law into a distribution
function of a strictly stable law with an arbitrary $\lambda$ one can use
remark 7 from the article [8], (see also[4, 17]).

##  4 Conclusion

The article considers the problem of calculating the distribution function of
a strictly stable law with the characteristic function (1) at large values of
the coordinate $x$. The need to solve this problem is dictated by the
inability of numerical integration algorithms to calculate the integral
correctly in the integral representation (32) at large $x$. The cause of such
difficulties lies in the behavior of the integrand. In this regard, the use of
the integral representation to calculate the distribution function at large
values of the coordinate is no longer possible, and it is necessary to apply
other approaches to solve this problem.

To solve it, it was proposed to use the expansion of the distribution function
in a power series at $x\to\infty$. In the article, such an expansion was
obtained, as well as an estimate for the remainder term. The results are
formulated in theorem 2.2. The convergence of this series has been studied and
it has been shown that in the case of $\alpha<1$ the series is convergent, in
the case $\alpha>1$ it is asymptotic, and in the case $\alpha=1$ the series is
convergent at $|x|>1$. These results are formulated as a corollary 2.4. It
should be noted that the results formulated in this corollary for the cases
$\alpha<1$ and $\alpha>1$ are not new and generalize the known results related
to the convergence of the expansion of the distribution function in a series
(see, for example[4, 5]). Nevertheless, the study of the convergence of the
series for the expansion of the distribution function in the case $\alpha=1$
was carried out for the first time. The study of this case showed that the
obtained series is convergent at $|x|>1$. In addition, we managed to show that
in this case at $N\to\infty$ this series converges to the distribution
function of the generalized Cauchy distribution (34).

The estimate of the remainder term obtained in theorem 2.2, turned out to be
very useful in the problem of calculating the distribution function. Using
this estimate, we managed to obtain the equation (27) for the threshold
coordinate $x_{\varepsilon}^{N}$. The threshold coordinate makes it possible
to determine the range of coordinates $x$ in which the absolute calculation
error does not exceed the required level of accuracy $\varepsilon$ at
specified values of $\alpha$ and the number of summands $N$ in the expansion
(5). As a result, the formula (29) is valid for calculating the distribution
function. The calculations performed showed that when using this formula, the
absolute error in calculating the distribution function in the domain
$|x|>x_{\varepsilon}^{N}$ does not exceed the required level of accuracy
$\varepsilon$, and in reality is much less than this value. As $|x|$
increases, the absolute calculation error will only decrease. This makes it
possible to use the formula (29) to calculate the distribution function even
for those values of $x$, for which the use of the integral representation (32)
turns out to be impossible. Indeed, the calculations have shown that for the
integral representation (32) there is a critical value of the coordinate
$x_{\mbox{\scriptsize cr}}$ at which numerical integration algorithms can no
longer calculate the integral correctly (see Fig. 5). At the same time, using
the formula (29) does not lead to any calculation difficulties. Thus, using
this formula to calculate the distribution function of a strictly stable law
in the coordinate range $|x|>x_{\mbox{\scriptsize cr}}$ solves the problem of
calculating the distribution function at large $x$.

As noted in the Introduction, the integral representation of the distribution
function (32) has two domains in which numerical methods have difficulties in
calculating the integral. These are the ranges of coordinates at $x\to 0$ and
at $x\to\infty$. This article shows that in the domain of large values of
coordinates, the formula (29) can be used for calculation. The problem of
calculating the distribution function at small values of $x$ was solved
earlier in the paper [15]. In this paper, the expansion of the distribution
function in a power series at $x\to 0$ was obtained and the area of
applicability of this expansion was determined. Thus, if to use the formula
(29) to calculate the distribution function at large values of $x$, at small
values of $x$ to use the results of the article [15], and to use in the
intermediate domain the integral representation (32), then we get an
opportunity to calculate correctly the distribution function of a strictly
stable law with the characteristic function (1) on the entire real line.

In conclusion, it should be pointed out that there are similar difficulties in
calculating the probability density of a strictly stable law. Calculation
difficulties also arise at small and large values of the coordinate $x$. The
problem of calculating the probability density is solved in the papers [15,
16]. These articles show that if in the domain of small coordinates $x$ to use
the probability density expansion from the article [15], in the domain of
large coordinates to use the expansion from the article [16], and in the
intermediate domain to use the integral representation for the density
probability obtained in the article [8], then we get an opportunity to
calculate the probability density of a strictly stable law correctly on the
entire real line. Thus, the problem of calculating the probability density and
distribution function of a strictly stable law with the characteristic
function (1) on the entire real line turns out to be solved.

The author thanks M. Yu. Dudikov for translation of the article into English.

##  Appendix A Integral representation of the distribution function

To perform the inverse Fourier transform and obtain the probability density
distribution, the following lemma is useful, which determines the inversion
formula

######  Lemma A.1.

The probability density distribution $g(x,\alpha,\theta)$ for any admissible
set of parameters $(\alpha,\theta)$ and any $x$ can be obtained using the
inverse transform formulas

$g(x,\alpha,\theta)=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-itx}\hat{g}(t,\alpha,\theta)dt=\left\\{\begin{array}[]{c}\displaystyle\frac{1}{\pi}\Re\int_{0}^{\infty}e^{itx}\hat{g}(t,\alpha,-\theta)dt,\\\
\displaystyle\frac{1}{\pi}\Re\int_{0}^{\infty}e^{-itx}\hat{g}(t,\alpha,\theta)dt.\end{array}\right.$
(31)

The proof of this lemma can be found in the article [8].

There is also an integral representation for the distribution function. For a
strictly stable law with a characteristic function (1) it was obtained in the
article [8] and formulated as a corollary

######  Corollary A.2.

The distribution function of the stable law $G(x,\alpha,\theta)$ with
characteristic function (1) can be represented in the form

  1. 1.

If $\alpha\neq 1$, then for any $|\theta|\leqslant\min(1,2/\alpha-1)$ and
$x\neq 0$

$G(x,\alpha,\theta)=\tfrac{1}{2}(1-\sign(x))+\sign(x)G^{(+)}(|x|,\alpha,\theta^{*}),$
(32)

where $\theta^{*}=\theta\sign(x)$,

$G^{(+)}(x,\alpha,\theta)=1-\frac{(1+\theta)}{4}(1+\sign(1-\alpha))\\\
+\frac{\sign(1-\alpha)}{\pi}\int_{-\pi\theta/2}^{\pi/2}\exp\left\\{-x^{\alpha/(\alpha-1)}U(\varphi,\alpha,\theta)\right\\}d\varphi,\quad
x>0,$ (33)

and $U(\varphi,\alpha,\theta)$ is determined by the expression

$U(\varphi,\alpha,\theta)=\left(\frac{\sin\left(\alpha\left(\varphi+\frac{\pi}{2}\theta\right)\right)}{\cos\varphi}\right)^{\alpha/(1-\alpha)}\frac{\cos\left(\varphi(1-\alpha)-\frac{\pi}{2}\alpha\theta\right)}{\cos\varphi}.$

  2. 2.

If $\alpha=1$, then for any $-1\leqslant\theta\leqslant 1$ and any $x$

$G(x,1,\theta)=\frac{1}{2}+\frac{1}{\pi}\arctan\left(\frac{x-\sin(\pi\theta/2)}{\cos(\pi\theta/2)}\right).$
(34)

  3. 3.

If $x=0$, then for any admissible $\alpha$ and $\theta$

$G(0,\alpha,\theta)=(1-\theta)/2.$

## References

  * [1] J. P. Nolan, An algorithm for evaluating stable densities in Zolotarev’s (M) parameterization, Mathematical and Computer Modelling 29 (10-12) (1999) 229–233.  doi:10.1016/S0895-7177(99)00105-3.
URL https://linkinghub.elsevier.com/retrieve/pii/S0895717799001053

  * [2] S. Ament, M. O’Neil, Accurate and efficient numerical calculation of stable densities via optimized quadrature and asymptotics, Statistics and Computing 28 (1) (2018) 171–185.  arXiv:1607.04247, doi:10.1007/s11222-017-9725-y.
URL https://link.springer.com/article/10.1007/s11222-017-9725-y

  * [3] V. M. Zolotarev, On the representation of stable laws by integrals, Selected Translations in Mathematical Statistics and Probability 4 (1964) 84–88.
  * [4] V. M. Zolotarev, One-dimensional stable Distributions, Amer. Mat. Soc., Providence, RI, 1986.
  * [5] V. V. Uchaikin, V. M. Zolotarev, Chance and stability Stable Distributions and their Applications, VSP, Utrecht, 1999.
  * [6] J. P. Nolan, Numerical calculation of stable densities and distribution functions, Communications in Statistics. Stochastic Models 13 (4) (1997) 759–774.  doi:10.1080/15326349708807450.
URL http://www.tandfonline.com/doi/abs/10.1080/15326349708807450

  * [7] J. P. Nolan, Computational aspects of stable distributions, WIREs Computational Statistics 14 (4) (2022) e1569.  doi:10.1002/wics.1569.
URL https://onlinelibrary.wiley.com/doi/10.1002/wics.1569

  * [8] V. V. Saenko, The Calculation of the Density and Distribution Functions of Strictly Stable Laws, Mathematics 8 (5) (2020) 775.  doi:10.3390/math8050775.
URL https://www.mdpi.com/2227-7390/8/5/775

  * [9] Y. Liang, W. Chen, A survey on computing Lévy stable distributions and a new MATLAB toolbox, Signal Processing 93 (1) (2013) 242–251.  doi:10.1016/j.sigpro.2012.07.035.
URL https://linkinghub.elsevier.com/retrieve/pii/S0165168412002654

  * [10] J. Royuela-del Val, F. Simmross-Wattenberg, C. Alberola-López, Libstable: Fast, parallel, and high-precision computation of $\alpha$-stable distributions in R, C/C++, and MATLAB, Journal of Statistical Software 78 (June) (2017).  doi:10.18637/jss.v078.i01.
URL https://www.jstatsoft.org/article/view/v078i01

  * [11] G. Julián-Moreno, J. E. López de Vergara, I. González, L. de Pedro, J. Royuela-del Val, F. Simmross-Wattenberg, Fast parallel $\alpha$-stable distribution function evaluation and parameter estimation using OpenCL in GPGPUs, Statistics and Computing 27 (5) (2017) 1365–1382.  doi:10.1007/s11222-016-9691-9.
URL http://link.springer.com/10.1007/s11222-016-9691-9

  * [12] R. H. Rimmer, J. P. Nolan, Stable Distributions in Mathematica, Mathematica journal 9 (4) (2005) 776–789.
  * [13] M. Veillette, MATLAB code: Alpha-Stable distributions (2008).
URL http://math.bu.edu/people/mveillet/research.html

  * [14] T. K. Pogány, S. Nadarajah, Remarks on the Stable S $\alpha$ ($\beta$,$\gamma$,$\mu$) Distribution, Methodology and Computing in Applied Probability 17 (2) (2015) 515–524.  doi:10.1007/s11009-014-9404-9.
URL http://link.springer.com/10.1007/s11009-014-9404-9

  * [15] V. V. Saenko, The Calculation of the Probability Density and Distribution Function of a Strictly Stable Law in the Vicinity of Zero, Mathematics 10 (20) (2022) 3861.  arXiv:2210.06920, doi:10.3390/math10203861.
URL https://www.mdpi.com/2227-7390/10/20/3861

  * [16] V. V. Saenko, The calculation of the probability density of a strictly stable law at large X (2023) 1–29arXiv:2303.03016.
URL http://arxiv.org/abs/2303.03016

  * [17] V. M. Zolotarev, V. V. Uchaikin, V. V. Saenko, Superdiffusion and stable laws, Journal of Experimental and Theoretical Physics 88 (4) (1999) 780–787.  doi:10.1134/1.558856.
URL http://www.springerlink.com/index/10.1134/1.558856

  * [18] H. Bateman, Higher Transcendental Functions., Vol. 1, McGraw-Hill Book Company, Inc, New York, 1953.
  * [19] M. R. Spiegel, S. Lipshutz, J. Liu, Schaum’s Outline Mathematical Handbook of Formulas and Tables., 5th Edition, McGraw-Hill Education, 2018.