An Application of Fractional Differential Equations to Risk Theory

This paper defines a new class of fractional differential operators alongside a family of random variables whose density functions solve fractional differential equations equipped with these operators. These equations can be further used to construct fractional integro-differential equations for the ruin probabilities in collective renewal risk models, with inter-arrival time distributions from the aforementioned family. Gamma-time risk models and fractional Poisson risk models are two specific cases among them, whose ruin probabilities have explicit solutions, when claim sizes distributions exhibit rational Laplace transforms.


Introduction
The concept of first passage time is widely used in financial mathematics and actuarial science. It could model various things, from the time to dividends payment of a stock, to the exercise date of an American put option, or the ruin probability of an insurance company. In this paper we focus on the ruin time of an insurance business, namely the first time in which the business surplus (capital) becomes negative. Our analysis is aimed at solving equations for the probability of ruin expressed as a function of the initial capital (surplus) of the risk process.
Motivated by risk theory applications, we consider a new class of risk processes, while extending those from [28,2,7] into a fractional derivative framework. It has been proved that ruin probabilities are exponential functions when claim sizes follow an exponential distribution, for various inter-arrival time distributions [4]. This paper will derive explicit ruin probabilities in risk models with claim sizes whose distributions have rational Laplace transforms, and inter-arrival time densities solving fractional differential equations. Gammatime risk models and fractional Poisson risk models are two particular cases among them. All the results are obtained due to the introduction of a new class of fractional differential operators, which extends those from [5,34]. These operators generalize the results from [2] to a fractional derivative framework, in which their explicit results concerning ruin probabilities become particular cases. Some existed ruin probability results are retrieved (see Example 1 and 3 for details), and new results are derived. For instance, in the gamma-time risk model with Erlang(2) distributed claim sizes, the ruin probability has the form where A 1 ,B 1 ,A 2 and B 2 are constants that can be calculated on a case-by-case basis (see Example 2).
The classical collective insurance risk model describes the surplus R(t) of an insurance company over time, where u > 0 is the initial capital and c > 0 is the premium rate. The claims occur randomly. The positive random variable X i describes the size of the i-th claim, which happened after waiting T i units of time since the last claim. The process N (t) gives the number of claims that have happened up to time t. In the classical model (1), dating back to [31,32,12], all random variables are assumed independent and identically distributed. Moreover, the waiting times are usually assumed to be exponentially distributed, with the resulting counting process N (t) thus being a Poisson process. The ruin probability of this compound Poisson risk model, for an initial capital u, is defined as The net profit condition is imposed to ensure that ruin does not happen with certainty. Various generalizations of the classical risk model (1) have been considered over time. In [38], Sparre Andersen defined the renewal risk model. This model accounts for claim number processes N (t) not necessarily Poisson, but verifying the renewal property. The ruin probabilities ψ(u) in renewal models still solve integral equations derived from the renewal property, namely (3) with the universal boundary condition lim u→∞ ψ(u) = 0, as in [18]. Here f T and F X denote the probability density of the waiting time, and the distribution function of the claim size, respectively. This notation will be used throughout the paper.
There is a large actuarial literature analyzing renewal risk processes. Expressions for the Laplace transform of the ruin probability for risk models with Erlang(2, β) or mixture of 2-exponential waiting times were derived in [13,14,15] as solutions of second-order differential equations. [30] calculated the joint and marginal moments of the time of ruin, the surplus before ruin, and the deficit at ruin, whenever the inter-arrival times distributions have rational Laplace-Stieltjes transform. Subsequently, [17] computed the Laplace transform of the non-ruin probability for inter-arrival times distributions exhibiting rational Laplace transforms. [28] used a similar approach as [20] to derive a defective renewal equation for the expected discounted penalty due at ruin in a risk model with Erlang(n) inter-arrival times. Finally, [9] derived linear ordinary differential equations for ruin probabilities in Poisson jumpdiffusion processes, with phase-type jumps and obtained explicit results in a few instances. The common thread of these paper consists on deriving the ruin probabilities as solutions of (integro-)differential equations.
In an attempt to develop a general method, [35,36] introduced two algebraic structures for treating integral operators in conjunction with derivatives, integro-differential operators and integro-differential polynomials. Their method allows the description of the associated differential equations, boundary conditions and solution operators (Green's operator) in a uniform yet formal language. Their algebraic symbolic structures have immediate applications in ruin theory. For instance, as an extension of the Erlang risk model, [2] transformed the integral equation for the expected-discounted-penalty-due-at-ruin function into an integro-differential equation whenever the inter-arrival time distributions have rational Laplace transforms. Rational Laplace transforms densities are equivalent to densities that are solutions of ordinary differential equations with constant coefficients. If the claim size distributions also have rational Laplace transforms, these integro-differential equations can be further reduced to linear boundary value problems. Their symbolic computation approach permits extensions to models with premium dependent on reserves (also discussed in [16] regarding the upper and lower bounds of finite ruin probabilities), the associated boundary problems involving then linear ordinary differential equations with variable coefficients [1]. A similar duality idea has been studied in [25] and the reference therein.
We show that the probability density function of a sum of independent, heterogeneous gamma and Mittag-Leffler random variables satisfies a fractional differential equation, which we write in an operator/symbolic form. As an application, we consider a family of risk models with inter-arrival times from this family of distributions, and derive the corresponding fractional integrodifferential equations satisfied by the corresponding ruin probabilities. We consider the case of claim sizes described by sums of heterogeneous gamma random variables and show that the corresponding ruin probabilities solve fractional differential equations with constant coefficients. These equations contain both left and right fractional differential operators. We annotate here that these equations can describe other physical phenomena exhibiting anomalous diffusion, as in [23] where the "claim sizes" are height losses of the granular material contained in a silo over time [27]. For other applications, we refer to [19,22,29,40] and the references therein. We also remark that Equation (18) presented in this paper can be seen as generalized cases of the fractional boundary problems treated by [24] where critical point theory is used to analyze fractional differential equations with Dirichlet boundary conditions. The gamma-time risk model considered here is the first generalization of the case of Erlang(n)-distributed waiting times considered in [28], to that of waiting times distributed as Γ (r, λ), r being now any positive real number. This is of significance since, in practice, parameter estimation methods usually yield non-integer-valued shape parameters for the gamma distributions that best fit the available data. It thus becomes necessary to study the ruin theory related to real-valued gamma-distributed random variables. [39] dealt with a special non-integer shape gamma Γ (1/b, 1/b), b > 1 distributed claims case, and [11] provided three equivalent expressions for ruin probabilities in a Cramér-Lundberg model with gamma distributed claims. Prior to this work, as far as we know, there are no results for non-integer shape gamma-time risk model in the ruin theory literature. The fractional Poisson risk model has been previously treated in [6] and [7] for exponential claim sizes, but here, via this fractional calculus approach, we are able to derive expressions for the ruin probability for a larger class of claim sizes in fractional Poisson models.
The paper is organized as follows. In Section 2 we introduce the concept of fractional integro-differential operators. In Section 3 we present the main result and finally, in Section 4, we perform some illustrative numerical calculations and compare the behavior of the ruin probabilities as a function of the model parameters, for both, the gamma distributed waiting time, and the fractional Poisson risk models. Appendix A contains all necessary background on fractional calculus.

Fractional Integro-Differential Operators
Let L(y) denote the n-th degree polynomial y n +p 1 y n−1 +· · ·+p n−1 y +p n and consider the following associated homogeneous ordinary differential equation with constant coefficients Suppose further that Equation (4) can be expressed in the form for positive real numbers λ j and integers k j , j = 1, . . . , m. In (4) and henceforth, denotes left-composition of operators, namely The solution f (x) to (5) is the probability density function of either a sum of Erlang random variables or a mixed Erlang random variable, depending on the boundary conditions (see [2]). We would like generalize Equation (5), and characterise its solutions in the case where the exponents k j are no longer integers.

Left and Right Fractional Differential Operators
In order to generalize expression (4), it is necessary to explore the world of fractional calculus. Solving fractional differential equations has become an essential issue as fractional-order models appear to be more adequate than previously used integer-order models in various fields. A large number of available analytical methods for solving fractional order integral and differential equations is discussed in [34], including the Mellin transform method, the power series method, and the symbolic method. The symbolic method was first introduced in [5] and generalizes the Laplace transform method: it uses a specific expansion (e.g., binomial or geometric) on the differential operator and writes it as an infinite sum of fractional derivatives. However, it is always necessary to check the validity of the formal expansion since the interchange of infinite summation and integration requires justification. It is nevertheless a powerful tool for determining the possible form of the solution. Numerous examples of the application of this method to heat and mass transfer problems are discussed in [5].
In this section we define a new family of operators based on the binomial expansion. All of the related definitions and propositions of fractional calculus can be found in Appendix A. The important motivation underlying the following definition comes from realising that for positive integer n and α ∈ R, and similarly for (− d dx + α) n . We thus define the following operators as the natural generalization in terms of fractional derivatives: and the right fractional differential operator (RFDO) α The domain of definition of α a R r x and α x R r b are those of the left Riemann-Liouville fractional derivatives a D r x and right Caputo fractional derivatives C x D r b respectively, which are given in Definition 4 and Definition 6.
In the case a = 0, integration by parts yields the following characterisation of the formal adjoint of α 0 R r x . Along with the integration by parts formula in Proposition 10, this is the key calculation needed for the proof of our main result.
Proposition 1 Let α ∈ R and r > 0. The formal adjoint with respect to integration by parts of the for appropriate functions f and g (see Proposition 10).
Note that the LFDO can be used to construct differential equations for probability density functions. Consider a gamma probability density function with shape parameter r ∈ R + and rate parameter λ ∈ R + , namely When r is not an integer, instead of an ordinary differential equation, the gamma density function solves the fractional differential equation with boundary conditions λ Another distribution related to the LFDO is the Mittag-Leffler distribution, which is the waiting time distribution in the fractional Poisson process (see in Appendix C). The Mittag-Leffler probability density function with parameter µ ∈ (0, 1] and λ ∈ R + is and solves the following fractional differential equation with the boundary condition 0 D µ−1 Here, the function E µ,µ is called two-parameter Mittag-Leffler function, which is defined in Equation (26).

A generalized family of random variables
The next theorem introduces the family of random variables to which the approach presented in this paper applies to. In its full generality, we consider random variables that can be written as finite sums of independent heterogeneous gamma and Mittag-Leffler random variables. At the moment, there is no known explicit formula for the probability density function of such a random variable, but one can always express it in a convolution form. Notice that if only gamma random variables with integer shape parameters are involved in the summation, this random variable is the generalized integer gamma distribution (GIG) [10]. We now characterise the fractional boundary value problem satisfied by the density function of such random variables.

Theorem 1 Consider a random variable T defined by
in terms of gamma random variables Then the density function f m,n T (t) of T solves the following fractional differential equation with boundary conditions (when n = 0) Here and subsequently Λ m,n denotes Proof We defer the proof of Theorem 1 to Appendix B.

Remark 1
We further assume that all λ 1,i are different, i.e., λ 1,i = λ 1,k for all i = k. In other words, each variable Y i has the gamma distribution with different rate parameters. The uniqueness of the λ 1,i , rate parameter of the gamma random variable could be realized without any loss of generality. Whenever we have λ 1,i = λ 1,k , i = k, we would consider the sum of their corresponding random variables, which is still a gamma random variable.
Remark 2 One can show that the boundary conditions in Theorem 1 have various equivalent expressions. For any positive integer number k m i=1 r i + n j=1 µ j , by choosing non-negative integers k 1,i and k 2,j such that m i=1 k 1,i + n j=1 k 2,j = k, we have the boundary conditions of Equation (11) as Remark 3 Equation (11) along with its boundary conditions can be regarded as the generalization of a pair of boundary problems discussed in [36]. When the fractional differential algebra is properly defined, these fractional-order boundary problems can be factorised and further solved by obtaining their corresponding Green's operators.
The solution to Equation (11) depends on the boundary condition. When different boundary conditions are given, we may obtain density functions for other possible random variables. For instance, let us consider the differential equation with two distinct sets of boundary conditions. First, if we impose the solution is the Elang(2, λ) density function f 2,0 T (t) = λ 2 te −λt which belongs to the random variable family considered in Equation (10). However, the solution to the above equation would become f 2,0 This solution is the density function of a mixture of an exponential and an Erlang random variable, and the associated distribution does not satisfy Equation (10).

Main Results
The LFDO and RFDO give us the ability to study a very general family of distributions that may find applications in various areas, e.g, queuing theory, risk theory and control theory. Although many of the available techniques for the analysis of the associated equations are numerical or asymptotic, the fractional differential approach still offers analytic insights to the related problems. In this section, we aim at accomplishing this with particular problems arising in risk theory. A special family of renewal risk models of the form (1) will be considered, including the Erlang(n) and fractional Poisson risk models. We will show that the ruin probabilities in these models solve fractional integro-differential equations involving the LFDO and RFDO operators. Before moving on to the main result, we introduce a lemma that allows us to change the argument of our operators on a bivariate function under certain circumstances.
Lemma 1 For positive real numbers α, r and c, the following identity holds where x and y are real numbers and α x R r ∞ is defined in Equation (7).
Proof We start from the left-hand side of Equation (13). By definition we have which is the right-hand side of Equation (13).
Now we are able to generalize the result from [28,2,7] to a risk model with inter-arrival times of the form of (10). The main result of this paper is the following: where the inter-arrival times T k are assumed to be a finite sum of independent gamma random variables Y i ∼ Γ (r i , λ 1,i ) and Mittag-Leffler random variables Z j ∼ ML(µ j , λ 2,j ) as in (10). Then the ruin probability ψ(u) under model R m,n , satisfies the following fractional integro-differential equation with the universal boundary condition lim u→∞ ψ(u) = 0. Here, the constant Λ m,n is given by Equation (12) and A * m,n is the formal adjoint of A m,n (see Equation (11)) and is given by Proof For a general renewal risk model, the ruin probability solves the renewal equation (3) (see [18]). Denoting the terms in parentheses of (3) by we now apply A * m,n (c d du ) on both sides of the renewal equation and use Lemma 1 to obtain The fractional integration by parts rule in Equation (27) is applicable here, The boundary condition term evaluated at t = 0 could be computed by using the initial value theorem of Laplace transforms, Another boundary condition term evaluated at t = ∞ also equals zero due to the fact that the definition of the right Caputo fractional derivative is an integral from t to ∞. Analogously, we are able to move the first n operators n j=1 ( C t D µj ∞ + λ 2,j ) from function h to f m,n T with all boundary conditions vanishing, which leads to Now we use the integration by parts formula in Proposition 1 to take the first RFDO λ1,1 t R r1 ∞ off of h. Furthermore it can be shown that its adjoint λ1,1 0 R r1 t commutes with ( 0 D µj t + λ 2,j ) for all j = 1, . . . , n when applied on the density function f m,n T . We therefore get the right-hand side equal to: The boundary condition at t = 0 can be computed by applying the initial value theorem .
We continue to iteratively use the initial value theorem on the terms , until it eventually gives us which tends to zero when s → ∞. The boundary condition term evaluated at t = ∞ yields zero since the right Caputo derivatives vanish at infinity. Analogously, we are able to move the rest of the operators m i=1 λ1,i t R ri ∞ from the function h to f m,n T with all boundary conditions vanishing, which leads to .
Since the inter-arrival time density satisfies Equation (11), the integral term of the above equation vanishes. The boundary conditions of f m,n T ensure that the lower summand is, at t = 0, This completes the proof.

Corollary 1
The non-ruin probability φ(u) = 1 − ψ(u) for the risk model in Theorem 2 satisfies the following fractional integro-differential equation with the universal boundary condition lim u→∞ φ(u) = 1 (see Equation (12) and (15) for the definitions of the constant Λ m,n and the operator A * m,n (c d du )). Theorem 2 characterises a fractional integro-diferential equation satisfied by the ruin probability ψ for a large class of waiting time distributions. The solvability of this fractional integro-differential equation depends on the particular form of the claim size distribution function F X .
We now restrict the rest of the analysis to claim sizes X i distributed as a sum of an arbitrary number of independent gamma random variables. The next theorem shows that under such assumption, Equation (14) can be written as a boundary value problem with only fractional derivatives. It is important to note that if the claim sizes include any Mittag-Leffler components, as it is the case of T in Theorem 2, we would have E[X i ] = ∞ and ruin would happen with probability one since the net profit condition is violated.

Theorem 3
Consider the renewal risk model in Theorem 2. Assume further that the claim sizes X i are each distributed as a sum of l independent Γ (s k , α k ) distributed random variables for some s k , α k > 0, k = 1, . . . , l i.e., f X satisfies with certain boundary conditions (see Theorem 1). Let A * m,n (c d du ) and Λ m,n be as defined in Equation (15) and (12) respectively. Then the non-ruin probability φ(u) satisfies with the universal boundary condition lim u→∞ φ(u) = 1 and initial values for k = 1, . . . , l k=1 s k − 1.
Proof Taking the operator A l ( d dy ) on both sides of Equation (16) leads to Recall from Theorem 1 that the non-ruin probability function φ(u) is supported on [0, ∞), so the identity holds in this case, giving Equation (18). For the boundary conditions, we compute where f 1 stands for the density function of Γ (s 1 , α 1 ). Using Equation (8), one has , which equals to zero for k = 1, . . . , l k=1 s k − 1. This completes the proof.

The characteristic equation method
Our next goal is solving the fractional differential boundary value problem in Theorem 3 via a characteristic equation starting from the ansatz φ(u) = e −zu . The main technical difficulty in dealing with the full generality of Theorem 3 arises from the fact that the operators in Equation (18) combine two different types of differential operators: A * m,n (c d du ) is a composition of right Caputo fractional derivatives, while the operators in A l ( d du ) are LFDOs which are ultimately defined in terms of left Riemann-Liouville fractional derivatives (see Equation (17), (15) and (1)). The proposed ansatz is an eigenfunction only for the operators in A * m,n (c d du ) (see Proposition 11 and Proposition 12). When restricting to the case of s k ∈ N , k = 1, . . . , l, we simplify things greatly, since reduces to a combination of ordinary differential operators. Note that assuming s k ∈ N , k = 1, . . . , l in Equation (17) is equivalent to assuming that the claim sizes X i are each distributed as a sum of l independent Erlang random variables. Moreover, under this case, the operator A l ( d du )A * m,n (c d du ) on the left hand side of Equation (18) is a composition of right Caputo fractional derivatives. Furthermore, with the ansatz φ(u) = e −zu , Equation (18) yields the following characteristic equation for z: Note that from the definition of Λ m,n in Equation (12), z = 0 is always a root of (20). If Equation (20) has N > 0 additional distinct complex roots with positive real part, say z 1 , . . . , z N , then the non-ruin probability φ that solves Equation (18) is where the constants K p , p = 1, . . . , N are to be determined from the boundary conditions in Equation (19), which are characterized in the following result.

0, if k > L(p).
Proof We consider the p-th boundary condition where f L(p) stands for the density function of a Γ (s L(p) , α L(p) ) random variable and f L(p)+ for the density function of a sum of random varibales with distributions Γ (s k , α k ), k = L(p) + 1, . . . , L. Let Φ = φ * f L(p)+ and apply Proposition 8 to compute Note that s p,L(p)−1 < s L(p) and we have Since this holds for all 1 p N , we complete the proof.

Explicit Expressions for Ruin Probabilities in Gamma-time and Fractional Poisson Risk Models
The class of models considered in Theorem 2 is very general. In this section, we thus focus on two specific models which might be of interest to applications, and where explicit forms of ruin (non-ruin) probabilities can be derived.

Remark 4
It has been shown [4] that for any renewal risk model, the ruin probability always has an exponential form when the claim distribution is exponential. However, the fractional differential equation approach bridges a solid connection between classical risk model and a class of renewal models, which might be applied in a more sophisticated model.

Gamma-time Risk Model
A gamma-time risk model, describes the reserve process R r (t) of an insurance company by replacing the Poisson process N (t) in the classical model (1) with a renewal counting process N r (t) with Γ (r, λ 1 ) distributed waiting times. This is a natural extension of Erlang(n) risk model consiered by [28]. Being a special case of Theorem 2, the equation for ruin probability ψ r (u) in gamma-time risk model is When claim sizes in this model have rational Laplace transforms, one could use the characteristic equation method mentioned in Section 3.1 to derive explicit ruin probabilities.
Example 1 In the gamma-time risk model with Γ (r, λ 1 ) distributed interarrival times and Exp(α) distributed claim sizes, the ruin probability equals to where x 2 > λ1 c is the larger root of equation where M X and M T are moment generating functions of claim sizes and interarrival times. This means that x 2 − λ1 c is the unique positive solution γ of the Lundberg's fundamental equation. This finding coincides with the result from [4] for renewal risk models with exponential claims.
In order to compare the classical compound Poisson with a gamma-time risk model, in Figure 1a we show numerically computed ruin probabilities in the case of Example 1 with different combinations of r and λ 1 when the mean claim inter-arrival time is fixed to r/λ 1 = 1.  (24)) for the ruin probability in Example 1 with continuously varying parameters r, λ 1 . The claim sizes have fixed exponential distribution with mean α = 1 and premium rate c = 1.2. The dotted line limits the region where the net profit condition r/λ 1 < c holds (see Equation (2)).
Note the substantial impact on ψ r (u) when changing the Poisson assumption (r = 1). Ruin is more likely to happen in the gamma-time risk model with a larger shape parameter r of inter-arrival times, and vice versa. The reason is that in this case, the expected inter-arrival time r/λ 1 is fixed whereas the variance of inter-arrival time r/λ 2 1 decreases as r increases, which means that the chance of having shorter waiting periods between claims will decrease. Since ruin is usually caused by not enough capital, the model with a larger shape parameter r is more likely to survive. Figure 1a coincides with the finding from [28], which focuses on Erlang(n) risk models.
In Figure 1b we illustrate the sensitivity to the parameters r and λ 1 of the ruin probability ψ r (u) in Example 1. In order to do this, we define the statistic Namely, u 5 is the minimum capital needed to ensure a ruin probability less than 5%. Note that any combinations of r and λ 1 on or above the dashed line, marking the net profit condition, will make the ruin certain. The value of u 5 tends to infinity as the parameters approach the dashed line since the safety loading c E(T ) E(X) − 1 tends to zero. When r takes large enough values or λ 1 take small enough values (in bluer areas), the ruin probability might be less than 5% even with zero initial capital. Note that along contour lines, dλ 1 ≈ 1 c dr, so the sensitivity of the ruin probabilities to its parameters depends almost exclusively on c.
The next example goes a step further and assumes gamma distributions for both the inter-arrival times and the claim sizes. This case is simple enough that the two positive roots of the characteristic equation can be bounded.
Example 2 In the gamma-time risk model with Γ (r, λ 1 ) distributed interarrival times and Γ (2, α) distributed claim sizes, the ruin probability equals to where z 3 > λ1 c + α > z 2 > λ1 c are the two largest roots of the equation

Fractional Poisson Risk Model
The fractional (compound) Poisson risk model is a special case of the classic compound Poisson risk model (1) where the counting process is chosen as fractional Poisson process N µ (t). Namely, the inter-arrival times are Mittag-Leffler distributed T ∼ ML(µ, λ 2 ), with λ 2 > 0, 0 < µ 1. Since, when µ = 1, the fractional Poisson process reduces into a Poisson process, we need the net profit condition to compute the ruin probability.
The following examples are derived under the assumption 0 < µ < 1 (in the fractional Poisson risk model). Note that in this case E[T i ] = ∞, so the net profit condition (2) holds whenever E[X i ] < ∞. It follows from Theorem 2 that the ruin probability ψ µ of a fractional Poisson risk model satisfies the following fractional integrodifferential equation with the universal boundary condition lim u→∞ ψ µ (u) = 0. Explicit expressions for ruin probabilities in the fractional Poisson risk model with exponential claims has been derived in [7]. The same result can be obtained via the fractional differential equation approach introduced in this paper.

Example 3
In the fractional Poisson risk model with T ∼ ML(µ, λ 2 ) and exponentially distributed claim sizes with parameter α, the ruin probability equals where x 2 is the unique positive solution of c µ x − αc µ + λ 2 x 1−µ = 0. Figure 2a shows the ruin probability ψ µ (u) for various combinations of the parameters λ 2 , µ, with fixed exponential claim size distribution.  (24)) for the ruin probability in Example 3 with continuously varying parameters µ, λ 2 . The claim sizes have fixed exponential distribution with mean α = 1 and premium rate c = 1.2.
Note the substantial impact on ψ µ (u) when the classical Poisson assumption (µ = 1) is changed. Increasing either λ 2 or µ increases the chances for ruin to happen. The reason is that, for large enough t, the expected number of jumps, before time t, in the fractional Poisson process (see Equation (28)), is an increasing function of both λ 2 and µ. Figure 2b shows the values of natural logarithm of u 5 as a function of µ and λ 2 . Note that the contour lines in this plot are not parallel to each other. As the values of µ decrease, the parameter λ 2 plays a less significant role in the ruin probability function.
Notice that the operator C u D µ ∞ tends to the identity operator, when µ → 0+. Thus, we obtain the following result.

A Basic facts from fractional calculus
The fractional calculus is the theory of integrals and derivatives of arbitrary order, which unifies and generalizes the notions of integer-order differentiation and n-fold integration [34]. The definitions of several special functions, fractional integrals and fractional derivatives used in this paper are listed in this section.

A.1 Mittag-Leffler Function
The Mittag-Leffler function was firstly introduced by [33] as a generalization of the exponential function.

A.2 Fractional Integrals and Derivatives
As per [21], we define and denote:

Definition 3
The left Riemann-Liouville fractional integral of order r > 0 with lower limit a ∈ R is defined on locally integrable functions f as a I r x f (x) = 1 Γ (r) x a (x − y) r−1 f (y) dy, x > a, and the right Riemann-Liouville fractional integral of order r > 0 with upper limit b ∈ R is defined as Definition 5 The Weyl-Liouville fractional derivatives [37,8] are special cases of the Riemann-Liouville derivatives, whenever a is replaced by −∞ or b is replaced by ∞ in Definition 4. The right Weyl-Liouville fractional derivative is defined for functions f ∈ L r ([a, b]) as (y − x) n−r−1 f (y) dy, n = r + 1.

Definition 6
The Caputo fractional derivatives are defined as fractional integrals on integer-order derivatives. The right Caputo fractional derivative is defined on functions f ∈ L r ([a, b]) as Proposition 11 The eigenfunction of left fractional derivative 0 D r x (or C 0 D r x ) is x 1−α E α,α (λx α ) with eigenvalue λ ∈ R (see Section 2.2.3 in [21]).

Proposition 12
The eigenfunction of right fractional derivative x D r ∞ (or C x D r ∞ ) is e −λx with eigenvalue λ r , where λ ∈ R + (see Section 4 in [41]). where λ is the intensity parameter and δ n,0 is the Kronecker symbol.