1 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 dividend payments 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 Li and Garrido [23], Albrecher et al. [2] and Biard and Saussereau [6] 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; see e.g. Asmussen and Albrecher [4, Theorem 2.1]. The present paper derives explicit ruin probabilities in risk models with claim sizes whose distributions have rational Laplace transforms, and with inter-arrival time densities solving fractional differential equations. Gamma-time 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 Podlubny [28, Sect. 6.3.2]. These operators generalise the results from Albrecher et al. [2] to a fractional derivative framework, in which their explicit results concerning ruin probabilities become particular cases. Some existing ruin probability results are retrieved (see Examples 4.2 and 4.5 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

$$ A_{1} e^{-B_{1} u}+A_{2} e^{-B_{2} u}, \quad u>0, $$

where \(A_{1},B_{1},A_{2}\) and \(B_{2}\) are constants that can be calculated on a case-by-case basis (see Example 4.4).

The classical collective insurance risk model describes the surplus\(R(t) \) of an insurance company over time as

$$ R(t)=u+ct-\sum _{i=1}^{N(t)}X_{i}, \quad t>0, $$
(1.1)

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 quantity \(N(t) \) gives the number of claims that have happened up to time \(t \). In the classical model (1.1), dating back to Lundberg [25, 26], Cramér [11], 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\) thus being a Poisson process. The ruin probability of this compound Poisson risk model, for an initial capital \(u\), is defined as

$$ \psi (u)=\mathbb{P}\left [\inf \{R(t):t>0\}< 0 \mid\! R(0)=u\right ]. $$

The net profit condition

$$ c\,\mathbb{E}\left [T_{i}\right ]>\mathbb{E}\left [X_{i}\right ] $$
(1.2)

is imposed to ensure that ruin does not happen with certainty. Various generalisations of the classical risk model (1.1) have been considered over time. Sparre Andersen [32] introduced the renewal risk model. This model accounts for claim number processes \(N\) not necessarily Poisson, but verifying the renewal property. The ruin probabilities \(\psi (u)\) in renewal models still solve integral equations derived from the renewal property, namely

$$\begin{aligned} \psi (u) &=\int _{0}^{\infty }f_{T}(t)\bigg(\int _{0}^{u+ct}\psi (u+ct-y) \,\mathrm{d}F_{X}(y)+\int _{u+ct}^{\infty }\,\mathrm{d}F_{X}(y)\bigg) \,\mathrm{d}t \end{aligned}$$
(1.3)

with the universal boundary condition \(\lim _{u\rightarrow \infty } \psi (u)=0\), as in Feller [17, Sect. 6.7]. 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 is used throughout the paper.

There is a large actuarial literature analysing renewal risk processes. Expressions for the Laplace transform of the ruin probability for risk models with \(\text{Erlang}(2, \beta )\) or mixture of 2-exponential waiting times were derived in Dickson [12] and Dickson and Hipp [13, 14] as solutions of second-order differential equations. Lin and Willmot [24] calculated the joint and marginal moments of the time of ruin, the surplus before ruin and the deficit at ruin, whenever the inter-arrival time distributions have rational Laplace–Stieltjes transforms. Subsequently, Dufresne [16] computed the Laplace transform of the non-ruin probability for inter-arrival time distributions exhibiting rational Laplace transforms. Li and Garrido [23] used a similar approach as Gerber and Shiu [18] to derive a defective renewal equation for the expected discounted penalty due at ruin in a risk model with \(\text{Erlang}(n)\) inter-arrival times. Finally, Chen et al. [8] derived linear ordinary differential equations for ruin probabilities in Poisson jump-diffusion processes with phase-type jumps and obtained explicit results in a few instances. The common thread of these papers consists of deriving the ruin probabilities as solutions of (integro-)differential equations.

In an attempt to develop a general method, Rosenkranz and Regensburger [29, 30] 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 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, Albrecher et al. [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 transform 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 premia dependent on reserves (also discussed in Djehiche [15] regarding the upper and lower bounds of finite-time ruin probabilities), the associated boundary problems then involving linear ordinary differential equations with variable coefficients; see Albrecher et al. [1]. A similar duality idea has been studied in Kolokoltsov and Lee [21] and the references 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 integro-differential 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 also remark that Eq. (3.6) presented in this paper can be seen as a generalised case of the fractional boundary problems treated by Jin and Liu [20], where critical point theory is used to analyse the fractional differential equations with Dirichlet boundary conditions.

The gamma-time risk model considered here is the first generalisation of the case of \(\text{Erlang}( n )\)-distributed waiting times considered in Li and Garrido [23] to that of waiting times distributed as \(\Gamma (r, \lambda )\), \(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. Thorin [33] dealt with a special case of claims with a non-integer shape gamma \(\Gamma (1/b,1/b)\), \(b>1\), distribution, and Constantinescu et al. [10] provided three equivalent expressions for ruin probabilities in a Cramér–Lundberg model with gamma-distributed claims. The fractional Poisson risk model has been previously treated in Beghin and Macci [5] and Biard and Saussereau [6] 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 organised as follows. In Sect. 2, we introduce the concept of fractional integro-differential operators. In Sect. 3, we present the main result, and finally, in Sect. 4, we perform some illustrative numerical calculations and compare the behaviour of the ruin probabilities as a function of the model parameters, for both the gamma-distributed waiting times and the fractional Poisson risk models. Appendix A contains all necessary background on fractional calculus.

2 Fractional integro-differential operators

Let \(\mathcal{L}(y)\) denote the \(n\)th degree polynomial \(y^{n}+p_{1}y ^{n-1}+\cdots +p_{n-1}y+p_{n}\) and consider the associated homogeneous ordinary differential equation with constant coefficients, given by

$$ \mathcal{L}\bigg(\frac{\mathrm{d}}{\mathrm{d}x}\bigg)f(x)=f^{(n)}(x)+p _{1}f^{(n-1)}(x)+\cdots +p_{n-1}f'(x)+p_{n}f(x)=0. $$
(2.1)

Suppose further that (2.1) can be expressed in the form

$$ \bigodot _{j=0}^{m} \bigg(\frac{\mathrm{d}}{\mathrm{d}x}+\lambda _{j} \bigg)^{k_{j}}f(x)=0 $$
(2.2)

for positive real numbers \(\lambda _{j} \) and integers \(k_{j} \), \(j=1,\dots , m \). In (2.1) and henceforth, ⨀ denotes left-composition of operators, namely

$$ \bigodot _{j=1}^{m} \mathcal{L}_{j} f := (\mathcal{L}_{m} \circ \cdots \circ \mathcal{L}_{1})f. $$

The solution \(f(x)\) to (2.2) 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 Albrecher et al. [2]). We should like to generalise (2.2) and characterise its solutions in the case where the exponents \(k_{j} \) are no longer integers.

2.1 Left and right fractional differential operators

In order to generalise (2.1), 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 Podlubny [28, Chaps. 5 and 6], including the Mellin transform method, the power series method and the symbolic method.

The symbolic method generalises 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, and there are numerous examples of the application of this method to heat and mass transfer problems.

In this section, we define a new family of operators based on the binomial expansion. All relevant definitions and results 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 \(\alpha \in \mathbb{R} \),

$$ \bigg(\frac{\mathrm{d}}{\mathrm{d}x}+\alpha \bigg)^{n} [f](x) = e^{- \alpha x} \frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}} \left (e^{\alpha x} f(x)\right ), $$

and similarly for \((-\frac{\mathrm{d}}{\mathrm{d}x}+\alpha )^{n} \). We thus define the following operators as the natural generalisation in terms of fractional derivatives.

Definition 1

Let \(r>0\), \(\alpha \in \mathbb{R}\), \(a\in [-\infty ,\infty )\) and \(b\in (-\infty ,\infty ]\). The left fractional differential operator (LFDO) R x r a α is defined by

R x r a α [f](x):= e α x a D x r ( e α x f(x)),

and the right fractional differential operator (RFDO) R b r x α by

R b r x α [g](x):= e α x x C D b r ( e α x g(x)).
(2.3)

The domains of definition of R x r a α and R b r x α are those of the left Riemann–Liouville fractional derivatives D x r a and right Caputo fractional derivatives D b r x C , respectively, which are given in Appendix A, points 2) and 4).

In the case \(a=0 \), integration by parts yields the following characterisation of the formal adjoint of R x r 0 α . Along with the integration by parts formula in (A.3), this is the key calculation needed for the proof of our main result.

Proposition 2

Let\(\alpha \in \mathbb{R} \)and\(r>0 \). The formal adjoint with respect to integration by parts of the LFDO R x r 0 α is the RFDO R r x α , i.e.,

0 0 α R x r [f](x)g(x)dx= 0 f ( x ) x α R r [g](x)dx,

for appropriate functions\(f\)and\(g\); see (A.3).

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\in \mathbb{R}_{+}\) and rate parameter \(\lambda \in \mathbb{R}_{+}\), i.e.,

$$ f_{r}(x) =\frac{\lambda ^{r}}{\Gamma (r)}x^{r-1}e^{-\lambda x}, \quad x>0. $$

When \(r\) is not an integer, instead of an ordinary differential equation, the gamma density function solves the fractional differential equation (see (A.1))

R x r 0 λ [ f r ](x)= e λ x 0 D x r ( e λ x f r ( x ) ) =0,x>0,
(2.4)

together with the boundary conditions R x r 1 0 λ [ f r ](0)= λ r and R x r k 0 λ [ f r ](0)=0 for \(k=2,\dots ,\lceil {r}\rceil \). Another distribution related to the LFDO is the Mittag-Leffler distribution, which is the waiting time distribution in the fractional Poisson process (see Appendix C). The Mittag-Leffler probability density function with parameters \(\mu \in (0,1]\) and \(\lambda \in \mathbb{R}_{+}\) is

$$ f_{\mu }(x) =\lambda x^{\mu -1}E_{\mu ,\mu }(-\lambda x^{\mu }), \quad t>0, $$

and solves the fractional differential equation

( 0 0 R x μ +λ)[ f μ ](x)= ( 0 D x μ +λ)[ f μ ](x)=0,x>0,
(2.5)

with the boundary condition D x μ 1 0 [f](0)=λ. Here, the function \(E_{\mu ,\mu }\) is called two-parameter Mittag-Leffler function; it is defined in (C.1).

2.2 A generalised family of random variables

The next theorem introduces the family of random variables to which the approach presented in this paper applies. 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 generalised integer gamma distribution (GIG) [9]. We now characterise the fractional boundary value problem satisfied by the density function of such random variables.

Theorem 3

Consider a random variable\(T\)defined by

$$ T=\sum _{i=1}^{m}Y_{i}+\sum _{j=1}^{n}Z_{j}, $$
(2.6)

in terms of gamma random variables\(Y_{i}\sim \Gamma (r_{i},\lambda _{1,i})\)and Mittag-Leffler random variables\(Z_{j}\sim \mathrm{ML}( \mu _{j},\lambda _{2,j})\), all independent of each other, where\(r_{i}\), \(\lambda _{1,i}\), \(\lambda _{2,j}\in \mathbb{R}_{+}\)and\(\mu _{j}\in (0,1]\). Then the density function\(f_{T}^{m,n}(t)\)of\(T\)solves the fractional differential equation

A m , n ( d d t ) [ f T m , n ](t):= j = 1 n ( 0 D t μ j + λ 2 , j ) i = 1 m 0 λ 1 , i R t r i [ f T m , n ](t)=0,
(2.7)

with boundary conditions (when\(n\neq 0\))

D t μ 1 1 0 j = 2 n ( 0 D t μ j + λ 2 , j ) i = 1 m 0 λ 1 , i R t r i [ f T m , n ] ( t ) | t = 0 = Λ m , n , D t μ 1 k 0 j = 2 n ( 0 D t μ j + λ 2 , j ) i = 1 m 0 λ 1 , i R t r i [ f T m , n ] ( t ) | t = 0 = 0 ,

for\(k=2,\dots ,\lceil \sum _{j=1}^{n} \mu _{j} +\sum _{i=1}^{m} r_{i} \rceil \). Here and subsequently, \(\Lambda _{m,n}\)denotes

$$ \Lambda _{m,n} :=\prod _{i=1}^{m} \lambda _{1,i}^{r_{i}}\prod _{j=1}^{n} \lambda _{2,j}. $$
(2.8)

Proof

We defer the proof of Theorem 2.3 to Appendix B. □

Remark 4

We further assume that all the \(\lambda _{1,i}\) are different, i.e., \(\lambda _{1,i}\neq \lambda _{1,k}\) for all \(i\neq k\). In other words, each variable \(Y_{i}\) has a gamma distribution with a different rate parameter. The uniqueness of the \(\lambda _{1,i}\) can be realised without any loss of generality. Whenever we have \(\lambda _{1,i}=\lambda _{1,k}\) for some \(i\neq k\), we can consider the sum of their corresponding random variables, which is still a gamma random variable with the same rate parameter.

Remark 5

One can show that the boundary conditions in Theorem 2.3 have various equivalent expressions. For any positive integer number \(k\leqslant \lceil \sum _{i=1}^{m} r_{i}+\sum _{j=1}^{n} \mu _{j}\rceil \), by choosing nonnegative integers \(k_{1,i}\) and \(k_{2,j}\) such that \(\sum _{i=1}^{m} k_{1,i} +\sum _{j=1}^{n} k_{2,j} =k\), we have the boundary conditions of (2.7) as

( j = 1 n ( 0 D t μ j k 2 , j + λ 2 , j 0 I t k 2 , j ) i = 1 m 0 λ 1 , i R t r i k 1 , i )[ f T m , n ](t) | t = 0 = { Λ m , n , k = 1 , 0 , k > 1 .

Remark 6

Equation (2.7) along with its boundary conditions can be regarded as the generalisation of a pair of boundary problems discussed in Rosenkranz and Regensburger [30]. When the fractional differential algebra is properly defined, these fractional-order boundary problems can be factorised and further solved by obtaining their corresponding Green operators.

The solution to (2.7) 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

$$ \left (\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right )^{2} f^{2,0} _{T}(t)=0 $$

with two distinct sets of boundary conditions. First, if we impose

$$ \textstyle\begin{cases} \displaystyle \left (\frac{\mathrm{d}}{\mathrm{d}x}+\lambda \right ) f^{2,0}_{T}(t) \bigg|_{t=0}=\lambda ^{2}, \\ \displaystyle \lambda f^{2,0}_{T}(t)\big|_{t=0}=0, \end{cases} $$

the solution is the \(\text{Erlang}(2,\lambda )\) density function \(f^{2,0} _{T}(t)=\lambda ^{2}te^{-\lambda t}\), which belongs to the random variable family considered in (2.6). However, the solution to the above equation would become \(f^{2,0}_{T}(t)=\frac{1}{2}\lambda e ^{-\lambda t}+\frac{1}{2}\lambda ^{2}te^{-\lambda t}\) if the boundary conditions were changed to

$$ \textstyle\begin{cases} \displaystyle \bigg(\frac{\mathrm{d}}{\mathrm{d}x}+\lambda \bigg)f^{2,0}_{T}(t) \bigg|_{t=0}= \displaystyle \frac{1}{2} \lambda ^{2}, \\ \displaystyle \lambda f^{2,0}_{T}(t)\big|_{t=0}= \displaystyle \frac{1}{2} \lambda ^{2}. \end{cases} $$

This solution is the density function of a mixture of an exponential and an Erlang random variable, and the associated distribution does not satisfy (2.6).

3 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.1) is considered, including the \(\text{Erlang}(n)\) and fractional Poisson risk models. We shall 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\(\alpha \), \(r\)and\(c\), we have the identity

R r x α [f(x+cy)](x,y)= c r y α c R r [f(x+cy)](x,y),
(3.1)

where\(x\)and\(y\)are real numbers and R r x α is defined in (2.3).

Proof

We start from the left-hand side of (3.1). By definition, we have

R r x α [f(x+cy)](x,y)= e α x 1 Γ ( n r ) x ( t x ) n r 1 d n d t n ( e α t f ( t + c y ) ) dt.

Letting \(s=\frac{1}{c}(t-x)+y\) leads to

$$\begin{aligned} \frac{1}{\Gamma (n-r)}\int _{y}^{\infty }e^{\alpha cy}(s-y)^{n-r-1}c ^{-r}\frac{\mathrm{d}^{n}}{\mathrm{d}y^{n}}\left (e^{-\alpha cs}f(cs+x)\right ) \,\mathrm{d}s, \end{aligned}$$

which is the right-hand side of (3.1). □

Now we are able to generalise the results from [23, 2, 6] to a risk model with inter-arrival times of the form (2.6). The main result of this paper is the following.

Theorem 2

Consider a renewal risk model

$$ R_{m,n}(t)=u+ct-\sum _{i=1}^{N_{m,n}(t)}X_{i}, \quad t>0, $$

where the inter-arrival times\(T_{k}\)are assumed to be a finite sum of independent gamma random variables\(Y_{i}\sim \Gamma (r_{i},\lambda _{1,i})\)and Mittag-Leffler random variables\(Z_{j}\sim \mathrm{ML}( \mu _{j},\lambda _{2,j})\)as in (2.6). Then the ruin probability\(\psi (u)\)under the model\(R_{m,n} \)satisfies the fractional integro-differential equation

$$\begin{aligned} &\mathcal{A}_{m,n}^{*}\left (c\frac{\mathrm{d}}{\mathrm{d}u}\right )[ \psi ](u) =\Lambda _{m,n}\left (\int _{0}^{u}\psi (u-y)\,\mathrm{d}F _{X}(y)+\int _{u}^{\infty }\,\mathrm{d}F_{X}(y)\right ) \end{aligned}$$
(3.2)

with the universal boundary condition\(\lim _{u\rightarrow \infty } \psi (u)=0\). Here, the constant\(\Lambda _{m,n} \)is given by (2.8), and\(\mathcal{A}_{m,n}^{*} \)is the formal adjoint of\(\mathcal{A}_{m,n} \) (see (2.7)) and is given by

A m , n ( c d d u ) := j = 1 n ( c μ j u C D μ j + λ 2 , j ) i = 1 m ( c r i u λ 1 , i / c R r i ).
(3.3)

Proof

For a general renewal risk model, the ruin probability solves the renewal equation (1.3) (see Feller [17, Sect. 6.7]). Denoting the terms in parentheses of (1.3) by

$$ h(u+ct)=\int _{0}^{u+ct}\psi (u+ct-y)\,\mathrm{d}F_{X}(y)+\int _{u+ct} ^{\infty }\,\mathrm{d}F_{X}(y), $$

we now apply \(\mathcal{A}_{m,n}^{*}(c\frac{\mathrm{d}}{\mathrm{d}u})\) on both sides of the renewal equation and use Lemma 3.1 to obtain

$$\begin{aligned} &\mathcal{A}_{m,n}^{*}\left (c\frac{\mathrm{d}}{\mathrm{d}u}\right )[ \psi ](u) =\int _{0}^{\infty }f_{T}^{m,n}(t) \mathcal{A}_{m,n}^{*} \left (\frac{\mathrm{d}}{\mathrm{d}t}\right )[h(u+ct)](u,t)\, \mathrm{d}t. \end{aligned}$$

The fractional integration by parts rule in (A.3) is applicable here, and yields

0 f T m , n ( t ) A m , n ( d d t ) [ h ( u + c t ) ] ( u , t ) d t = 0 ( 0 D t μ 1 + λ 2 , 1 ) [ f T m , n ] ( t ) A m , n 1 ( d d t ) [ h ( u + c t ) ] ( u , t ) d t + k = 0 μ 1 ( ( 1 ) μ 1 + 1 + k 0 D t μ 1 + k μ 1 1 [ f T m , n ] ( t ) = : : + ( A m , n 1 ( d d t ) [ h ( u + c t ) ] ( u , t ) | 0 ) .

The boundary condition term evaluated at \(t=0\) could be computed by using the initial value theorem of Laplace transforms,

I t 1 μ 1 0 [ f T m , n ](0)= lim s ( s μ 1 j = 1 n λ 2 , j s μ j + λ 2 , j i = 1 m ( λ 1 , i s + λ 1 , i ) r i )=0.

Another boundary condition term evaluated at \(t=\infty \) also equals zero due to the fact that the definition of the right Caputo fractional derivative is an integral from \(t\) to \(\infty \). Analogously, we are able to move the first \(n\) operators j = 1 n ( t C D μ j + λ 2 , j ) from function \(h\) to \(f_{T}^{m,n}\) with all boundary conditions vanishing, which leads to

A m , n ( c d d u ) [ ψ ] ( u ) = 0 j = 1 n ( 0 D t μ j + λ 2 , j ) [ f T m , n ] ( t ) i = 1 m t λ 1 , i R r i [ h ( u + c t ) ] ( u , t ) d t .

Now we use the integration by parts formula in Proposition 2.2 to take the first RFDO R r 1 t λ 1 , 1 off \(h\). Furthermore, it can be shown that its adjoint R t r 1 0 λ 1 , 1 commutes with ( 0 D t μ j + λ 2 , j ) for all \(j=1,\dots , n\) when applied on the density function \(f_{T}^{m,n}\). We therefore get the right-hand side equal to

0 j = 1 n ( 0 D t μ j + λ 2 , j ) 0 λ 1 , 1 R t r 1 [ f T m , n ] ( t ) i = 2 m t λ 1 , i R r i [ h ( u + c t ) ] ( u , t ) d t + k = 0 r 1 ( ( 1 ) r 1 + 1 + k i = 2 m t λ 1 , i R r i [ h ( u + c t ) ] ( u , t ) + ( j = 1 n ( 0 D t μ j + λ 2 , j ) 0 λ 1 , i R t r 1 + k r 1 1 [ f T m , n ] ( t ) ) | 0 .

The boundary condition at \(t=0\) can be computed by applying the initial value theorem to give

j = 1 n ( 0 D t μ j + λ 2 , j ) 0 λ 1 , 1 R t r 1 + k r 1 1 [ f T m , n ] ( 0 ) = j = 1 n λ 2 , j lim s ( λ 1 , 1 r 1 s ( s + λ 1 , 1 ) r 1 + 1 k i = 2 m ( λ 1 , i s + λ 1 , i ) r i = : j = 1 n λ 2 , j lim s ( s = 0 k 1 ( s + λ 1 , 1 ) ( 0 D t r 1 + k r 1 2 ( e λ 1 , 1 f T m , 0 ( t ) ) ) | t = 0 ) .

We continue to iteratively use the initial value theorem on the terms

s ( s + λ 1 , 1 ) ( 0 D t r 1 + k r 1 2 ( e λ 1 , 1 t f T m , 0 (t))) | t = 0 ,

until it eventually gives us

s ( s + λ 1 , 1 ) r 1 1 ( 0 I t r 1 + 1 r 1 ( e λ 1 , 1 t f T m , 0 (t))) | t = 0 =s ( s + λ 1 , 1 ) r 1 2 i = 1 m ( λ 1 , i s ) r i ,

which tends to zero when \(s\rightarrow \infty \). The boundary condition term evaluated at \(t=\infty \) yields zero since the right Caputo derivatives vanish at infinity. Analogously, we are able to move the rest of the operators i = 1 m t λ 1 , i R r i from the function \(h\) to \(f_{T}^{m,n}\) with all boundary conditions vanishing, which leads to

$$\begin{aligned} \mathcal{A}_{m,n}^{*}\left (c\frac{\mathrm{d}}{\mathrm{d}u}\right )[ \psi ](u) &=\int _{0}^{\infty }\mathcal{A}_{m,n}\left (\frac{ \mathrm{d}}{\mathrm{d}t}\right ) [f_{T}^{m,n}](t) [h(u+ct)](u,t)\, \mathrm{d}t \\ & \phantom{=:}+\bigg(h(u+ct) \mathcal{A}_{m-1,n}\Big(\frac{\mathrm{d}}{ \mathrm{d}t}\Big) [f_{T}^{m,n}](t)\bigg|_{t=0}\bigg).\vadjust{\eject} \end{aligned}$$

Since the inter-arrival time density satisfies (2.7), the integral term of the above equation vanishes. The boundary conditions of \(f_{T}^{m,n}\) ensure that the lower summand is, at \(t=0 \),

h(u) j = 1 n ( 0 D t μ j + λ 2 , j ) 0 λ 1 , n R t r n 1 i = 1 m 1 0 λ 1 , i R t r i [ f T m , n ](0)= Λ m , n h(u).

This completes the proof. □

Corollary 3

The non-ruin probability\(\phi (u) = 1-\psi (u)\)for the risk model in Theorem3.2satisfies the fractional integro-differential equation

$$\begin{aligned} &\mathcal{A}_{m,n}^{*}\left (c\frac{\mathrm{d}}{\mathrm{d}u}\right )[ \phi ](u) =\Lambda _{m,n}\left (\int _{0}^{u}\phi (u-y)\,\mathrm{d}F _{X}(y)\right ) \end{aligned}$$
(3.4)

with the universal boundary condition\(\lim _{u\rightarrow \infty } \phi (u)=1\) (see (2.8) and (3.3) for the definitions of the constant\(\Lambda _{m,n} \)and the operator\(\mathcal{A}_{m,n}^{*}(c\frac{\mathrm{d}}{\mathrm{d}u}) \)).

Theorem 3.2 characterises a fractional integro-differential equation satisfied by the ruin probability \(\psi \) 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 this assumption, (3.2) 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 is the case of \(T\) in Theorem 3.2, we have \(\mathbb{E}[X_{i}] = \infty \) and ruin happens with probability one since the net profit condition is violated.

Theorem 4

Consider the renewal risk model in Theorem3.2. Assume further that the claim sizes\(X_{i}\)are each distributed as a sum of\(\ell \)independent\(\Gamma (s_{k},\alpha _{k})\)-distributed random variables for some\(s_{k},\alpha _{k} >0 \), \(k=1,\dots ,\ell \), i.e., \(f_{X}\)satisfies

A ( d d u ) [ f X ] (u):= k = 1 0 α k R u s k [ f X ] (u)=0,
(3.5)

with boundary conditions (see Theorem2.3)

( 0 α 1 R u s 1 1 k = 2 0 α k R u s k ) [ f X ] ( u ) | u = 0 = k = 1 α k s k , ( 0 α 1 R u s 1 q k = 2 0 α k R u s k ) [ f X ] ( u ) | u = 0 = 0 ,

for\(q=2,\dots ,\lceil \sum _{k=1}^{\ell }s_{k} \rceil \). Let\(\mathcal{A}_{m,n}^{*}(c\frac{\mathrm{d}}{\mathrm{d}u}) \)and\(\Lambda _{m,n} \)be as in (3.3) and (2.8), respectively. Then the non-ruin probability\(\phi (u)\)satisfies

$$\begin{aligned} &\mathcal{A}_{\ell }\left (\frac{\mathrm{d}}{\mathrm{d}u}\right ) \mathcal{A}_{m,n}^{*}\left (c\frac{\mathrm{d}}{\mathrm{d}u}\right )[ \phi ](u) =\Lambda _{m,n} \prod _{k=1}^{\ell }\alpha _{k}^{s_{k}} \phi (u) \end{aligned}$$
(3.6)

with the universal boundary condition\(\lim _{u\rightarrow \infty } \phi (u)=1\)and initial values

( 0 α 1 R u s 1 k k = 2 0 α k R u s k j = 1 n ( c μ j u C D μ j + λ 2 , j ) i = 1 m ( c r i u λ 1 , i / c R r i ))[ϕ](0)=0,
(3.7)

for\(k'=1,\,\dots ,\,\lceil \sum _{k=1}^{\ell }s_{k} \rceil -1\).

Proof

Taking the operator \(\mathcal{A}_{\ell }( \frac{\mathrm{d}}{\mathrm{d}y})\) on both sides of (3.4) leads to

$$\begin{aligned} &\mathcal{A}_{\ell }\left (\frac{\mathrm{d}}{\mathrm{d}u}\right ) \mathcal{A}_{m,n}^{*}\left (c\frac{\mathrm{d}}{\mathrm{d}u}\right )[ \phi ](u) = \Lambda _{m,n} \mathcal{A}_{\ell }\left (\frac{\mathrm{d}}{ \mathrm{d}u}\right )\left (\int _{0}^{u}\phi (u-y)f_{X}(y)\,\mathrm{d}y\right ). \end{aligned}$$

Recall from Theorem 2.3 that the non-ruin probability function \(\phi (u)\) is supported on \([0,\infty )\); so the identity

A ( d d u ) ( 0 u ϕ ( u y ) f X ( y ) d y ) = k = 1 0 α k R u s k [ϕ f X ](u)= k = 1 α s k ϕ(u)

holds in this case, giving (3.6). For the boundary conditions, we compute

( 0 α 1 R u s 1 k k = 2 0 α k R u s k j = 1 n ( c μ j u C D μ j + λ 2 , j ) i = 1 m ( c r i u λ 1 , i / c R r i ) ) [ ϕ ] ( 0 ) = Λ m , n k = 2 α k s k 0 α 1 R u s 1 k ( ϕ ( u ) f 1 ( u ) ) | u = 0 ,

where \(f_{1}\) stands for the density function of \(\Gamma (s_{1},\alpha _{1})\). Using (A.2), one has

Λ m , n k = 2 α k s k e α 1 u 0 D u s 1 k ( 0 u e α 1 ( u y ) ϕ ( u y ) e α 1 y f 1 ( y ) d y ) | u = 0 = Λ m , n k = 2 α k s k ( e α 1 u ( e α 1 u ϕ ( u ) α 1 s 1 Γ ( k ) u k 1 ) | u = 0 + ϕ ( 0 ) α 1 s 1 Γ ( k + 1 ) y k | y = 0 ) ,

which equals zero for \(k'=1,\,\dots ,\,\lceil \sum _{k=1}^{\ell }s_{k} \rceil -1\). This completes the proof. □

3.1 The characteristic equation method

Our next goal is to solve the fractional differential boundary value problem in Theorem 3.4 via a characteristic equation starting from the ansatz \(\phi (u) = e^{-z u}\). The main technical difficulty in dealing with the full generality of Theorem 3.4 arises from the fact that the operators in (3.6) combine two different types of differential operators: \(\mathcal{A}_{m,n}^{*}(c\frac{ \mathrm{d}}{\mathrm{d}u}) \) is a composition of right Caputo fractional derivatives, while the operators in \(\mathcal{A}_{\ell }(\frac{ \mathrm{d}}{\mathrm{d}u}) \) are LFDOs which are ultimately defined in terms of left Riemann–Liouville fractional derivatives (see (3.5), (3.3) and (2.1)). The proposed ansatz is an eigenfunction only for the operators in \(\mathcal{A}_{m,n} ^{*}(c\frac{\mathrm{d}}{\mathrm{d}u}) \) (see (A.4)). When restricting to the case of \(s_{k} \in \mathbb{N}\), \(k=1,\dots , \ell \), we simplify things greatly since

A ( d d u ) = k = 1 0 α k R u s k = k = 1 ( d d u + α k ) s k

reduces to a combination of ordinary differential operators.

Note that assuming \(s_{k} \in \mathbb{N}\), \(k=1,\dots ,\ell \), in (3.5) is equivalent to assuming that the claim sizes \(X_{i} \) are each distributed as a sum of \(\ell \) independent Erlang random variables. Moreover, in this case, the operator \(\mathcal{A} _{\ell }(\frac{\mathrm{d}}{\mathrm{d}u}) \mathcal{A}_{m,n}^{*}(c\frac{ \mathrm{d}}{\mathrm{d}u}) \) on the left-hand side of (3.6) is a composition of right Caputo fractional derivatives. Furthermore, with the ansatz \(\phi (u) = e^{-z u} \), (3.6) yields for \(z \) the characteristic equation

$$ \prod _{k=1}^{\ell }(-z+\alpha _{k})^{s_{k}}\prod _{j=1}^{n}(c^{\mu _{j}}z ^{\mu _{j}}+\lambda _{2,j})\prod _{i=1}^{m}(cz+\lambda _{1,i})^{r_{i}}= \Lambda _{m,n} \prod _{k=1}^{\ell }\alpha _{k}^{s_{k}}. $$
(3.8)

Note that from the definition of \(\Lambda _{m,n} \) in (2.8), \(z=0 \) is always a root of (3.8). If (3.8) has \(N>0\) additional distinct complex roots with positive real parts, say \(z_{1},\dots ,z_{N}\), then the non-ruin probability \(\phi \) that solves (3.6) is

$$ \phi (u) = 1+ \sum _{p=1}^{N} K_{p} e^{-z_{p} u}, $$
(3.9)

where the constants \(K_{p}\), \(p=1,\dots ,N \), are to be determined from the boundary conditions in (3.7), which are characterised in the following result.

Proposition 5

Suppose\(s_{k} \in \mathbb{N}\), \(k=1,\dots ,\ell \), in Theorem3.4. The number of initial-value boundary conditions of\(\phi (u)\)is\(N=\sum _{k=1}^{\ell }s_{k}\), and they are given explicitly by

k = 1 0 α k R u s p , k j = 1 n ( c μ j u C D μ j + λ 2 , j ) i = 1 m ( c r i u λ 1 , i / c R r i )[ϕ](0)=0,p=1,,N,
(3.10)

where the values of\(s_{p,k}\)are to be computed as follows: let

$$ L(p)=\inf \bigg\{ \ell \in \mathbb{N}: \sum _{k=1}^{\ell }s_{k} \leqslant p\bigg\} , \quad p=1,\dots ,N, $$

and define

$$ s_{p,k}= \textstyle\begin{cases} s_{k}, & \textit{if } k< L(p), \\ p-\sum _{i=1}^{L(p)-1}s_{i}-1, &\textit{if } k=L(p), \\ 0, & \textit{if } k>L(p). \end{cases} $$

Proof

We consider the \(p\)th boundary condition

k = 1 L ( p ) 0 α k R u s p , k A m , n ( c d d u ) [ ϕ ] ( 0 ) = Λ m , n k = 1 L ( p ) 1 α k s k 0 α L ( p ) R u s p , L ( p ) [ ϕ f L ( p ) f L ( p ) + ] ( 0 ) ,

where \(f_{L(p)}\) stands for the density function of a \(\Gamma (s_{L(p)}, \alpha _{L(p)})\) random variable and \(f_{L(p)+}\) for the density function of a sum of random variables with distributions \(\Gamma (s_{k},\alpha _{k})\), \(k = L(p)+1,\dots , L\). Let \(\Phi =\phi *f_{L(p)+}\) and apply (A.2) to compute

R u s p , L ( p ) 0 α L ( p ) [ Φ f L ( p ) ] ( u ) = Φ ( u ) 0 D y s p , L ( p ) 1 ( e α L ( p ) y f L ( p ) ( y ) ) | y = 0 = : + e α L ( p ) u ( e α L ( p ) ( u ) Φ ( u ) 0 D u s p , L ( p ) e α L ( p ) u f L ( p ) ( u ) ) .

Note that \(s_{p,L(p)-1}< s_{L(p)}\) and we have

R u s p , L ( p ) 0 α L ( p ) [Φ f L ( p ) ](0)= 0 u Φ ( u y ) 0 α L ( p ) R y s p , L ( p ) f L ( p ) (y)dy | u = 0 =0.

Since this holds for all \(1\leqslant p\leqslant N\), we complete the proof. □

Substituting the expression in (3.9) for \(\phi (u) \) into the boundary conditions (3.10) yields explicit linear equations for the unknown constants \(K_{p} \), \(p=1,\dots ,N \). First, denote

$$ \Delta _{p} :=\prod _{j=1}^{n}(c^{\mu _{j}}z_{p}^{\mu _{j}}+ \lambda _{2,j})\prod _{i=1}^{m}(cz_{p}+\lambda _{1,i})^{r_{i}}, \quad p =1,\dots ,N. $$

Then the constants \(K_{p} \), \(p=1,\dots ,N \) in (3.9) satisfy

$$ \textstyle\begin{cases} {\Lambda _{m,n}+\sum _{p=1}^{N}\Delta _{p} K_{p}=0}, \\ {\alpha _{1}\Lambda _{m,n}+\Delta \sum _{p=1}^{N}(-z_{p}+\alpha _{1}) K_{p}=0}, \\ {\quad \vdots } \\ {\alpha _{1}^{s_{1}}\Lambda _{m,n}+\sum _{p=1}^{N}\Delta _{p}(-z_{p}+\alpha _{1})^{s_{1}} K_{p}=0}, \\ {\alpha _{1}^{s_{1}}\alpha _{2}\Lambda _{m,n}+\sum _{p=1}^{N}\Delta _{p}(-z_{p}+\alpha _{1})^{s_{1}}(-z_{p}+\alpha _{2}) K_{p}=0}, \\ {\quad \vdots } \\ {\alpha _{1}^{s_{1}}\alpha _{2}^{s_{2}}\Lambda _{m,n}+\sum _{p=1}^{N}\Delta _{p}(-z_{p}+\alpha _{1})^{s_{1}}(-z_{p}+\alpha _{2})^{s_{2}} K_{p}=0}, \\ {\quad \vdots } \\ {\prod _{k=1}^{\ell -1} \alpha _{k}^{s_{k}}\alpha _{\ell }^{s_{\ell }-1}\Lambda _{m,n}+\sum _{p=1}^{N}\Delta _{p}\prod _{k=1}^{\ell -1}(-z_{p}+\alpha _{k})^{s_{k}}(-z_{p}+\alpha _{\ell })^{s_{\ell }-1} K_{p}=0.} \end{cases} $$

4 Explicit expressions for ruin probabilities in gamma-time and fractional Poisson risk models

The class of models considered in Theorem 3.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 1

It has been shown in Asmussen and Albrecher [4, Theorem 2.1] 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 the classical risk model and a class of renewal models which might be applied in a more sophisticated model.

4.1 Gamma-time risk model

A gamma-time risk model describes the reserve process \(R_{r}\) of an insurance company by replacing the Poisson process \(N\) in the classical model (1.1) with a renewal counting process \(N_{r}\) with \(\Gamma (r,\lambda _{1})\)-distributed waiting times. This is a natural extension of the \(\text{Erlang}(n)\) risk model considered by Li and Garrido [23].

Being a special case of Theorem 3.2, the equation for the ruin probability \(\psi _{r}(u)\) in the gamma-time risk model is

c r e λ 1 c u u C D r ( e λ 1 c u ψ r (u))= λ 1 r ( 0 u ψ r ( u y ) d F X ( y ) + u d F X ( y ) ) .

When the claim sizes in this model have rational Laplace transforms, one could use the characteristic equation method mentioned in Sect. 3.1 to derive explicit ruin probabilities.

Example 2

In the gamma-time risk model with \(\Gamma (r,\lambda _{1})\)-distributed inter-arrival times and \(\text{Exp}(\alpha )\)-distributed claim sizes, the ruin probability equals

$$ \psi _{r}(u)=\left (\frac{\lambda _{1}}{cx_{2}}\right )^{r}e^{-(x_{2}-\frac{ \lambda _{1}}{c}) u}, \quad u>0, $$
(4.1)

where \(x_{2}>\frac{\lambda _{1}}{c}\) is the larger root of the equation

$$ c^{r} x^{r}\bigg(x-\Big(\frac{\lambda _{1}}{c}+\alpha \Big)\bigg)+ \alpha \lambda _{1}^{r}=0. $$

Remark 3

Letting \(s=x_{2}-\frac{\lambda _{1}}{c}\) in (4.1), one has

$$\begin{aligned} \big(M_{X}(s)M_{T}(-cs)\big)^{-1}-1 &=\bigg(1-\frac{s}{\alpha }\bigg) \bigg(1+\frac{cs}{\lambda _{1}}\bigg)^{r}-1 \\ &=\frac{c^{r}}{\lambda _{1}^{r}\alpha }\bigg(\Big(\alpha +\frac{ \lambda _{1}}{c}-x_{2}\Big)x_{2}^{r}-\frac{\lambda _{1}^{r}}{c^{r}} \bigg) \\ &=\frac{-1}{\lambda _{1}^{r}\alpha }\big(c^{r} x_{2}^{r+1}- (c^{r-1} \lambda _{1}+\alpha c^{r} )x_{2}^{r}+\alpha \lambda _{1}^{r}\big)=0, \end{aligned}$$

where \(M_{X}\) and \(M_{T}\) are the moment-generating functions of the claim sizes and inter-arrival times. This means that \(x_{2}-\frac{ \lambda _{1}}{c}\) is the unique positive solution \(\gamma \) of Lundberg’s fundamental equation. This finding coincides with the result from [4, Theorem 2.1] for renewal risk models with exponential claims.

In order to compare the classical compound Poisson with a gamma-time risk model, we show in Fig. 1(a) numerically computed ruin probabilities in the case of Example 4.2 with different combinations of \(r \) and \(\lambda _{1} \) when the mean claim inter-arrival time is fixed to \(r/\lambda _{1}=1\).

Fig. 1
figure 1

(a) Ruin probabilities in the case of Example 4.2 for \(\lambda _{1}=r \in \{ 0.5, 1, 1.5, 2, 2.5\}\). Claim sizes are exponentially distributed with mean \(\alpha = 1\) and \(c = 1.2\) in order to ensure the net profit condition. (b) Natural logarithm of \(u_{5} \) (see (4.2)) for the ruin probability in Example 4.2 with continuously varying parameters \(r, \lambda _{1} \). The claim sizes have a fixed exponential distribution with mean \(\alpha = 1\) and premium rate \(c = 1.2\). The dotted line limits the region where the net profit condition \(r/\lambda _{1} < c\) holds (see (1.2))

Note the substantial impact on \(\psi _{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/\lambda _{1}\) is fixed whereas the variance \(r/\lambda _{1}^{2}\) of the inter-arrival time 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 1(a) coincides with the findings from [23], which focuses on \(\text{Erlang}(n)\) risk models.

In Fig. 1(b), we illustrate the sensitivity to the parameters \(r\) and \(\lambda _{1} \) of the ruin probability \(\psi _{r}(u) \) in Example 4.2. In order to do this, we define the statistic

$$ u_{5} :=\inf \left \{u \geq 0: \psi _{r}(u) < 0.05\right \}. $$
(4.2)

Namely, \(u_{5} \) is the minimum capital needed to ensure a ruin probability of less than 5%. Note that any combinations of \(r\) and \(\lambda _{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 \(\frac{c\,\mathbb{E}[T]}{\mathbb{E}[X]}-1\) tends to zero. When \(r\) takes large enough values or \(\lambda _{1}\) takes small enough values (in darker blue areas), the ruin probability might be less than 5% even with zero initial capital. Note that along contour lines, \(\mathrm{d}\lambda _{1} \approx \frac{1}{c} \, \mathrm{d} r \) so that the sensitivity of the ruin probability 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 4

In the gamma-time risk model with \(\Gamma (r,\lambda _{1})\)-distributed inter-arrival times and \(\Gamma (2,\alpha )\)-distributed claim sizes, the ruin probability equals

$$ \psi _{r}(u)=\frac{\frac{\lambda _{1}}{c}-z_{3}}{z_{2}-z_{3}}\left (\frac{ \lambda _{1}}{cz_{2}}\right )^{r}e^{(\frac{\lambda _{1}}{c}-z_{2}) u}+\frac{\frac{ \lambda _{1}}{c}-z_{2}}{z_{3}-z_{2}}\left (\frac{\lambda _{1}}{cz_{3}}\right ) ^{r}e^{(\frac{\lambda _{1}}{c}-z_{3}) u}, \quad u>0, $$

where \(z_{3}>\frac{\lambda _{1}}{c}+\alpha >z_{2}> \frac{\lambda _{1}}{c}\) are the two largest roots of the equation

$$ c^{r} z^{r} \bigg(z-\Big(\frac{\lambda _{1}}{c}+\alpha \Big)\bigg)^{2}- \alpha ^{2} \lambda _{1}^{r}=0. $$

4.2 Fractional Poisson risk model

The fractional (compound) Poisson risk model is a special case of the classic compound Poisson risk model (1.1) where the counting process is chosen as fractional Poisson process \(N_{\mu }\). Namely, the inter-arrival times are Mittag-Leffler-distributed \(T \sim \text{ML}(\mu , \lambda _{2})\) with \(\lambda _{2} >0 \), \(0 < \mu \leqslant 1 \). Since for \(\mu =1\), the fractional Poisson process reduces to a Poisson process, we need the net profit condition to compute the ruin probability. The following examples are derived under the assumption \(0<\mu <1\) (in the fractional Poisson risk model). Note that in this case \(\mathbb{E} [T_{i}] = \infty \), so the net profit condition (1.2) holds whenever \(\mathbb{E} [X_{i}] < \infty \). It follows from Theorem 3.2 that the ruin probability \(\psi _{\mu }\) of a fractional Poisson risk model satisfies the fractional integro-differential equation

c μ u C D μ ψ μ (u)+ λ 2 ψ μ (u)= λ 2 ( 0 u ψ μ ( u y ) d F X ( y ) + u d F X ( y ) ) ,

with the universal boundary condition \(\lim _{u\rightarrow \infty } \psi _{\mu }(u)=0\). Explicit expressions for ruin probabilities in the fractional Poisson risk model with exponential claims have been derived in [6]. The same result can be obtained via the fractional differential equation approach introduced in this paper.

Example 5

In the fractional Poisson risk model with \(T \sim \mathrm{ML}(\mu , \lambda _{2}) \) and exponentially distributed claim sizes with parameter \(\alpha \), the ruin probability equals

$$ \psi _{\mu }(u)=\bigg(1-\frac{x_{2}}{\alpha }\bigg)e^{-x_{2} u}, \quad u>0, $$

where \(x_{2}\) is the unique positive solution of \(c^{\mu }x-\alpha c ^{\mu }+\lambda _{2} x^{1-\mu }=0\).

Figure 2(a) shows the ruin probability \(\psi _{\mu }(u) \) for various combinations of the parameters \(\lambda _{2}, \mu \), with fixed exponential claim size distribution.

Fig. 2
figure 2

(a) Ruin probabilities in the case of Example 4.5 for different combinations of \(\lambda _{2}, \mu \). Claim sizes are taken exponentially distributed with mean \(\alpha = 1\) and \(c = 1.2\). (b) Natural logarithm of \(u_{5} \) (see Eq. (4.2)) for the ruin probability in Example 4.5 with continuously varying parameters \(\mu , \lambda _{2} \). The claim sizes have fixed exponential distribution with mean \(\alpha = 1\) and premium rate \(c = 1.2\)

Note the substantial impact on \(\psi _{\mu }(u)\) when the classical Poisson assumption (\(\mu = 1\)) is changed. Increasing either \(\lambda _{2}\) or \(\mu \) 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 (C.2)) is an increasing function of both \(\lambda _{2}\) and \(\mu \). Figure 2(b) shows the values of the natural logarithm of \(u_{5}\) as a function of \(\mu \) and \(\lambda _{2}\). Note that the contour lines in this plot are not parallel to each other. As the values of \(\mu \) decrease, the parameter \(\lambda _{2}\) plays a less significant role in the ruin probability function.

Notice that the operator D μ u C tends to the identity operator when \(\mu \rightarrow 0+\). Thus we obtain the following result.

Corollary 6

In the fractional Poisson risk model, the ruin probability\(\psi _{ \mu }(u)\)converges to a function\(\psi _{0}(u)\)as\(\mu \rightarrow 0\). Moreover, the function\(\psi _{0}(u)\)satisfies the integral equation

$$ \left (1+\lambda _{2}\right )\psi _{0}(u)=\lambda _{2}\int _{0}^{u}\psi _{0}(u-y)\,\mathrm{d}F_{X}(y)+\lambda _{2}\int _{u}^{\infty }\mathrm{d}F _{X}(y), $$
(4.3)

with the universal boundary condition\(\lim _{u\rightarrow \infty }\psi _{0}(u)=0\).

Substituting \(u=0\) into (4.3) gives \(\psi _{0}(0)=\frac{\lambda _{2}}{\lambda _{2}+1}\) which only depends on the value of \(\lambda _{2}\). Taking Laplace transforms of both sides with respect to \(u\) leads to

$$ \hat{\psi }_{0}(s)=\frac{1-\hat{f}(s)}{(\lambda _{2}+1)s-\lambda _{2} s \hat{f}(s)}, $$

which can be explicitly inverted back in some cases.