An Application of Fractional Differential Equations to Risk Theory

TL;DR

This paper introduces fractional differential operators and their application to risk models, deriving explicit ruin probability solutions for specific cases with rational Laplace transforms.

Contribution

It develops a new class of fractional differential operators and applies them to risk theory, providing explicit solutions for ruin probabilities in fractional risk models.

Findings

Explicit ruin probability solutions for gamma-time risk models.

Explicit ruin probability solutions for fractional Poisson risk models.

Applicable when claim size distributions have rational Laplace transforms.

Abstract

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.

Equations

$$A_1{e}^{-B_{1}u}+A_2{e}^{-B_{2}u},u>0,$$

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

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

$$c\mathbb{E}\left[T_i\right]>\mathbb{E}\left[X_i\right]$$

$$\psi (u)$$

$$\mathcal{L}\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)[f](x)=f^{(n)}(x)+p_1{f}^{(n-1)}(x)+\cdots +p_{n-1}{f}^{\prime }(x)+p_{n}f(x)=0.$$

$$\underset{j=0}{\overset{m}{⨀}}{\left(\frac{\mathrm{d}}{\mathrm{d}x}+{\lambda }_j\right)}^{k_j}[f](x)=0$$

$$\underset{j=1}{\overset{m}{⨀}}{\mathcal{L}}_j[f]:=({\mathcal{L}}_m\circ \cdots \circ {\mathcal{L}}_1)[f].$$

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

$$\text{prescript}\alpha a{\mathrm{R}}_x^r\left[f\right](x)\text{vbox}..=e^{-\alpha x}\text{prescript}a{\mathrm{D}}_x^r\left[e^{\alpha x}f(x)\right]$$

$$\text{prescript}\alpha x{\mathrm{R}}_b^r\left[g\right](x)\text{vbox}..=e^{\alpha x}\text{prescript}\mathrm{C}x{\mathrm{D}}_b^r\left[e^{-\alpha x}g(x)\right].$$

$${\int }_0^{\infty }\text{prescript}\alpha 0{\mathrm{R}}_x^r[f](x)g(x)\mathrm{d}x={\int }_0^{\infty }f(x)\text{prescript}\alpha x{\mathrm{R}}_{\infty }^r[g](x)\mathrm{d}x,$$

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

$$\text{prescript}\lambda 0{\mathrm{R}}_x^r[f_r](x)=e^{-\lambda x}\text{prescript}0{\mathrm{D}}_x^r\left[e^{\lambda x}{f}_r(x)\right]=0,x>0,$$

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

$$\left(\text{prescript}00{\mathrm{R}}_x^{\mu }+\lambda \right)[f_{\mu }](x)=(\text{prescript}0{\mathrm{D}}_x^{\mu }+\lambda )[f_{\mu }](x)=0,x>0,$$

$$T=\sum _{i=1}^m{Y}_i+\sum _{j=1}^n{Z}_j,$$

$${\mathcal{A}}_{m,n}\left(\frac{\mathrm{d}}{\mathrm{d}t}\right)\left[f_T^{m,n}\right](t)\text{vbox}..=\underset{j=1}{\overset{n}{⨀}}\left(\text{prescript}0{\mathrm{D}}_t^{{\mu }_j}+{\lambda }_{2,j}\right)\underset{i=1}{\overset{m}{⨀}}\text{prescript}{\lambda }_{1,i}0{\mathrm{R}}_t^{r_i}\left[f_T^{m,n}\right](t)=0,$$

$$\text{prescript}0{\mathrm{D}}_t^{{\mu }_1-1}\underset{j=2}{\overset{n}{⨀}}\left(\text{prescript}0{\mathrm{D}}_t^{{\mu }_j}+{\lambda }_{2,j}\right)\underset{i=1}{\overset{m}{⨀}}\text{prescript}{\lambda }_{1,i}0{\mathrm{R}}_t^{r_i}{[f_T^{m,n}](t)|}_{t=0}={\Lambda }_{m,n},$$

$$\text{prescript}0{\mathrm{D}}_t^{{\mu }_1-k}\underset{j=2}{\overset{n}{⨀}}\left(\text{prescript}0{\mathrm{D}}_t^{{\mu }_j}+{\lambda }_{2,j}\right)\underset{i=1}{\overset{m}{⨀}}\text{prescript}{\lambda }_{1,i}0{\mathrm{R}}_t^{r_i}{[f_T^{m,n}](t)|}_{t=0}=0,$$

$${\Lambda }_{m,n}\text{vbox}..=\prod _{i=1}^m{\lambda }_{1,i}^{r_i}\prod _{j=1}^n{\lambda }_{2,j}.$$

$$\left(\bigodot_{j=1}^{n}\left(\prescript{}{0}{\mathrm{D}}_{t}^{\mu_{j}-k_{2,j}}+\lambda_{2,j}\prescript{}{0}{\mathrm{I}}_{t}^{k_{2,j}}\right)\bigodot_{i=1}^{m}\prescript{\lambda_{1,i}}{0}{\mathrm{R}}_{t}^{r_{i}-k_{1,i}}\right)\left.[f_{T}^{m,n}](t)\right|_{t=0}=\begin{cases*}\Lambda_{m,n},\,k=1\\ \quad\\ 0,\quad\,\,\,\,k>1.\end{cases*}$$

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

$$\{\begin{array}{l}\left(\frac{\mathrm{d}}{\mathrm{d}x}+\lambda \right)f_T^{2,0}(t){|}_{t=0}={\lambda }^2,\\ \\ \lambda f_T^{2,0}(t){|}_{t=0}=0,\end{array}$$

$$\{\begin{array}{l}\left(\frac{\mathrm{d}}{\mathrm{d}x}+\lambda \right)f_T^{2,0}(t){|}_{t=0}=\frac{1}{2}{\lambda }^2,\\ \\ \lambda f_T^{2,0}(t){|}_{t=0}=\frac{1}{2}{\lambda }^2.\end{array}$$

$$\text{prescript}\alpha x{\mathrm{R}}_{\infty }^r[f(x+cy)](x,y)=c^{-r}\text{prescript}\alpha cy{\mathrm{R}}_{\infty }^r[f(x+cy)](x,y),$$

$$\text{prescript}\alpha x{\mathrm{R}}_{\infty }^r[f(x+cy)](x,y)=e^{\alpha x}\frac{1}{\Gamma (n-r)}{\int }_x^{\infty }(t-x{)}^{n-r-1}\frac{{\mathrm{d}}^n}{\mathrm{d}{t}^n}\left[e^{-\alpha t}f(t+cy)\right]\mathrm{d}t.$$

$$\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,$$

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

$${\mathcal{A}}_{m,n}^{\ast }\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)$$

$${\mathcal{A}}_{m,n}^{\ast }\left(c\frac{\mathrm{d}}{\mathrm{d}u}\right)\text{vbox}..=\underset{j=1}{\overset{n}{⨀}}\left(c^{{\mu }_j}\text{prescript}\mathrm{C}u{\mathrm{D}}_{\infty }^{{\mu }_j}+{\lambda }_{2,j}\right)\underset{i=1}{\overset{m}{⨀}}\left(c^{r_i}\text{prescript}{\lambda }_{1,i}/cu{\mathrm{R}}_{\infty }^{r_i}\right).$$

Full text

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 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 li2004ruin ; albrecher2010algebraic ; biard2014fractional 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 AsmAlb2010 . 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. 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 babenko1986heat ; podlubny1998fractional . These operators generalize the results from albrecher2010algebraic 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

[TABLE]

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,

[TABLE]

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 lundberg1903approximerad ; lundberg1926forsakringsteknisk ; Cramer , 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

[TABLE]

The net profit condition

[TABLE]

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 andersen1957collective , 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 $\psi(u)$ in renewal models still solve integral equations derived from the renewal property, namely

[TABLE]

with the universal boundary condition $\lim_{u\rightarrow\infty}\psi(u)=0$, as in feller2008introduction . 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, $\beta$) or mixture of 2-exponential waiting times were derived in dickson1998class ; dickson1998ruin ; dickson2001time as solutions of second-order differential equations. lin2000moments 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, dufresne2002general computed the Laplace transform of the non-ruin probability for inter-arrival times distributions exhibiting rational Laplace transforms. li2004ruin used a similar approach as gerber1998time 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, chen2007ode 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 paper consists on deriving the ruin probabilities as solutions of (integro-)differential equations.

In an attempt to develop a general method, rosenkranz2008integro ; rosenkranz2008solving 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, albrecher2010algebraic 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 djehiche1993large regarding the upper and lower bounds of finite ruin probabilities), the associated boundary problems involving then linear ordinary differential equations with variable coefficients albrecher2013exact . A similar duality idea has been studied in kolokoltsov2013stochastic 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 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 annotate here that these equations can describe other physical phenomena exhibiting anomalous diffusion, as in jiao2012existence where the “claim sizes” are height losses of the granular material contained in a silo over time leszczynski2011modeling . For other applications, we refer to fix2004least ; jiao2011existence ; li2013existence ; torres2014mountain 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 jin2016eigenvalue 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 li2004ruin , 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. MR0381218 dealt with a special non-integer shape gamma $\Gamma(1/b,1/b)$, $b>1$ distributed claims case, and cons2017 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 beghin2013large and biard2014fractional 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.

2 Fractional Integro-Differential Operators

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

[TABLE]

Suppose further that Equation (4) can be expressed in the form

[TABLE]

for positive real numbers $\lambda_{j}$ and integers $k_{j}$, $j=1,\dots,m$. In (4) and henceforth, $\bigodot$ denotes left-composition of operators, namely

[TABLE]

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 albrecher2010algebraic ). We would like generalize Equation (5), 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 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 podlubny1998fractional , including the Mellin transform method, the power series method, and the symbolic method.

The symbolic method was first introduced in babenko1986heat 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 babenko1986heat .

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 $\alpha\in\mathbb{R}$,

[TABLE]

and similarly for $(-\frac{\mathrm{d}}{\mathrm{d}x}+\alpha)^{n}$. We thus define the following operators as the natural generalization 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) $\prescript{\alpha}{a}{\mathrm{R}}_{x}^{r}$ is defined by

[TABLE]

and the right fractional differential operator (RFDO) $\prescript{\alpha}{x}{\mathrm{R}}_{b}^{r}$ by

[TABLE]

The domain of definition of $\prescript{\alpha}{a}{\mathrm{R}}_{x}^{r}$ and $\prescript{\alpha}{x}{\mathrm{R}}_{b}^{r}$ are those of the left Riemann-Liouville fractional derivatives $\prescript{}{a}{\mathrm{D}}_{x}^{r}$ and right Caputo fractional derivatives $\prescript{\mathrm{C}}{x}{\mathrm{D}}_{b}^{r}$ 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 $\prescript{\alpha}{0}{\mathrm{R}}_{x}^{r}$. 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 $\alpha\in\mathbb{R}$ and $r>0$. The formal adjoint with respect to integration by parts of the LFDO $\prescript{\alpha}{0}{\mathrm{R}}_{x}^{r}$ is the RFDO $\prescript{\alpha}{x}{\mathrm{R}}_{\infty}^{r}$, namely,

[TABLE]

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\in\mathbb{R}^{+}$ and rate parameter $\lambda\in\mathbb{R}^{+}$, namely

[TABLE]

When $r$ is not an integer, instead of an ordinary differential equation, the gamma density function solves the fractional differential equation

[TABLE]

with boundary conditions $\prescript{\lambda}{0}{\mathrm{R}}_{x}^{r-1}[f_{r}](0)=\lambda^{r}$ and $\prescript{\lambda}{0}{\mathrm{R}}_{x}^{r-k}[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 in Appendix C). The Mittag-Leffler probability density function with parameter $\mu\in(0,1]$ and $\lambda\in\mathbb{R}^{+}$ is

[TABLE]

and solves the following fractional differential equation

[TABLE]

with the boundary condition $\prescript{}{0}{\mathrm{D}}_{x}^{\mu-1}[f](0)=\lambda$. Here, the function $E_{\mu,\mu}$ is called two-parameter Mittag-Leffler function, which is defined in Equation (26).

2.2 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) coelho1998generalized . We now characterise the fractional boundary value problem satisfied by the density function of such random variables.

Theorem 2.1

Consider a random variable $T$ defined by

[TABLE]

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. Here $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 following fractional differential equation

[TABLE]

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

[TABLE]

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

[TABLE]

Proof

We defer the proof of Theorem 2.1 to Appendix B.

Remark 1

We further assume that all $\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 the gamma distribution with different rate parameters. The uniqueness of the $\lambda_{1,i}$, rate parameter of the gamma random variable could be realized without any loss of generality. Whenever we have $\lambda_{1,i}=\lambda_{1,k}$, $i\neq 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 2.1 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 non-negative 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 Equation (11) as

[TABLE]

Remark 3

Equation (11) along with its boundary conditions can be regarded as the generalization of a pair of boundary problems discussed in rosenkranz2008solving . 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

[TABLE]

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

[TABLE]

the solution is the Elang($2,\lambda$) density function $f^{2,0}_{T}(t)=\lambda^{2}te^{-\lambda t}$ which belongs to the random variable family considered in Equation (10). 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 are changed to

[TABLE]

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).

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) 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 $\alpha$, $r$ and $c$, the following identity holds

[TABLE]

where $x$ and $y$ are real numbers and $\prescript{\alpha}{x}{\mathrm{R}}_{\infty}^{r}$ is defined in Equation (7).

Proof

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

[TABLE]

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

[TABLE]

which is the right-hand side of Equation (13).

Now we are able to generalize the result from li2004ruin ; albrecher2010algebraic ; biard2014fractional to a risk model with inter-arrival times of the form of (10). The main result of this paper is the following:

Theorem 3.1

Consider a renewal risk model

[TABLE]

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 (10). Then the ruin probability $\psi(u)$ under model $R_{m,n}$, satisfies the following fractional integro-differential equation

[TABLE]

with the universal boundary condition $\lim_{u\rightarrow\infty}\psi(u)=0$. Here, the constant $\Lambda_{m,n}$ is given by Equation (12) and $\mathcal{A}_{m,n}^{*}$ is the formal adjoint of $\mathcal{A}_{m,n}$ (see Equation (11)) and is given by

[TABLE]

Proof

For a general renewal risk model, the ruin probability solves the renewal equation (3) (see feller2008introduction ). Denoting the terms in parentheses of (3) by

[TABLE]

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

[TABLE]

The fractional integration by parts rule in Equation (10) is applicable here,

[TABLE]

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

[TABLE]

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 $\bigodot_{j=1}^{n}(\prescript{\mathrm{C}}{t}{\mathrm{D}}_{\infty}^{\mu_{j}}+\lambda_{2,j})$ from function $h$ to $f_{T}^{m,n}$ with all boundary conditions vanishing, which leads to

[TABLE]

Now we use the integration by parts formula in Proposition 1 to take the first RFDO $\prescript{\lambda_{1,1}}{t}{\mathrm{R}}_{\infty}^{r_{1}}$ off of $h$. Furthermore it can be shown that its adjoint $\prescript{\lambda_{1,1}}{0}{\mathrm{R}}_{t}^{r_{1}}$ commutes with $(\prescript{}{0}{\mathrm{D}}_{t}^{\mu_{j}}+\lambda_{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:

[TABLE]

The boundary condition at $t=0$ can be computed by applying the initial value theorem

[TABLE]

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

[TABLE]

until it eventually gives us

[TABLE]

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 $\bigodot_{i=1}^{m}\prescript{\lambda_{1,i}}{t}{\mathrm{R}}_{\infty}^{r_{i}}$ from the function $h$ to $f_{T}^{m,n}$ with all boundary conditions vanishing, which leads to

[TABLE]

Since the inter-arrival time density satisfies Equation (11), 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$,

[TABLE]

This completes the proof.

Corollary 1

The non-ruin probability $\phi(u)=1-\psi(u)$ for the risk model in Theorem 3.1 satisfies the following fractional integro-differential equation

[TABLE]

with the universal boundary condition $\lim_{u\rightarrow\infty}\phi(u)=1$ (see Equation (12) and (15) 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.1 characterises a fractional integro-diferential 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 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 3.1, we would have $\mathbb{E}[X_{i}]=\infty$ and ruin would happen with probability one since the net profit condition is violated.

Theorem 3.2

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

[TABLE]

with certain boundary conditions (see Theorem 2.1). Let $\mathcal{A}_{m,n}^{*}(c\frac{\mathrm{d}}{\mathrm{d}u})$ and $\Lambda_{m,n}$ be as defined in Equation (15) and (12) respectively. Then the non-ruin probability $\phi(u)$ satisfies

[TABLE]

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

[TABLE]

for $k^{\prime}=1,\,\dots,\,\left\lceil\sum_{k=1}^{l}s_{k}\right\rceil-1.$

Proof

Taking the operator $\mathcal{A}_{l}(\frac{\mathrm{d}}{\mathrm{d}y})$ on both sides of Equation (16) leads to

[TABLE]

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

[TABLE]

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

[TABLE]

where $f_{1}$ stands for the density function of $\Gamma(s_{1},\alpha_{1})$. Using Equation (8), one has

[TABLE]

which equals to zero for $k^{\prime}=1,\,\dots,\,\lceil\sum_{k=1}^{l}s_{k}\rceil-1$. This completes the proof.

3.1 The characteristic equation method

Our next goal is solving the fractional differential boundary value problem in Theorem 3.2 via a characteristic equation starting from the ansatz $\phi(u)=e^{-zu}$. The main technical difficulty in dealing with the full generality of Theorem 3.2 arises from the fact that the operators in Equation (18) 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}_{l}(\frac{\mathrm{d}}{\mathrm{d}u})$ 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 $\mathcal{A}_{m,n}^{*}(c\frac{\mathrm{d}}{\mathrm{d}u})$ (see Proposition 11 and Proposition 12). When restricting to the case of $s_{k}\in\mathbb{N}$ , $k=1,\dots,l$, we simplify things greatly, since

[TABLE]

reduces to a combination of ordinary differential operators.

Note that assuming $s_{k}\in\mathbb{N}$ , $k=1,\dots,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 $\mathcal{A}_{l}(\frac{\mathrm{d}}{\mathrm{d}u})\mathcal{A}_{m,n}^{*}(c\frac{\mathrm{d}}{\mathrm{d}u})$ on the left hand side of Equation (18) is a composition of right Caputo fractional derivatives. Furthermore, with the ansatz $\phi(u)=e^{-zu}$, Equation (18) yields the following characteristic equation for $z$:

[TABLE]

Note that from the definition of $\Lambda_{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},\dots,z_{N}$, then the non-ruin probability $\phi$ that solves Equation (18) is

[TABLE]

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

Proposition 2

Suppose $s_{k}\in\mathbb{N}$, $k=1,\dots,l$, in Theorem 3.2. The number of initial-value boundary conditions of $\phi(u)$ is $N=\sum_{k=1}^{l}s_{k}$ and they are given explicitly by:

[TABLE]

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

[TABLE]

and define

[TABLE]

Proof

We consider the $p$-th boundary condition

[TABLE]

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 varibales with distributions $\Gamma(s_{k},\alpha_{k}),$ $k=L(p)+1,\dots,L$. Let $\Phi=\phi*f_{L(p)+}$ and apply Proposition 8 to compute

[TABLE]

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

[TABLE]

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

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

[TABLE]

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

[TABLE]

4 Explicit Expressions for Ruin Probabilities in Gamma-time and Fractional Poisson Risk Models

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

4.1 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 $\Gamma(r,\lambda_{1})$ distributed waiting times. This is a natural extension of Erlang$(n)$ risk model consiered by li2004ruin .

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

[TABLE]

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 $\Gamma(r,\lambda_{1})$ distributed inter-arrival times and Exp$(\alpha)$ distributed claim sizes, the ruin probability equals to

[TABLE]

where $x_{2}>\frac{\lambda_{1}}{c}$ is the larger root of equation

[TABLE]

Remark 5

Let $s=x_{2}-\frac{\lambda_{1}}{c}$ in Equation (23), one has

[TABLE]

where $M_{X}$ and $M_{T}$ are moment generating functions of claim sizes and inter-arrival times. This means that $x_{2}-\frac{\lambda_{1}}{c}$ is the unique positive solution $\gamma$ of the Lundberg’s fundamental equation. This finding coincides with the result from AsmAlb2010 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 $\lambda_{1}$ when the mean claim inter-arrival time is fixed to $r/\lambda_{1}=1$.

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 of inter-arrival time $r/\lambda_{1}^{2}$ 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 li2004ruin , which focuses on Erlang($n$) risk models.

In Figure 1b we illustrate the sensitivity to the parameters $r$ and $\lambda_{1}$ of the ruin probability $\psi_{r}(u)$ in Example 1. In order to do this, we define the statistic

[TABLE]

Namely, $u_{5}$ is the minimum capital needed to ensure a ruin probability 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}$ 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, $\mathrm{d}\lambda_{1}\approx\frac{1}{c}\,\mathrm{d}r$, 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 $\Gamma(r,\lambda_{1})$ distributed inter-arrival times and $\Gamma(2,\alpha)$ distributed claim sizes, the ruin probability equals to

[TABLE]

where $z_{3}>\frac{\lambda_{1}}{c}+\alpha>z_{2}>\frac{\lambda_{1}}{c}$ are the two largest roots of the equation

[TABLE]

4.2 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_{\mu}(t)$. 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, when $\mu=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<\mu<1$ (in the fractional Poisson risk model). Note that in this case $\mathbb{E}[T_{i}]=\infty$, so the net profit condition (2) holds whenever $\mathbb{E}[X_{i}]<\infty$. It follows from Theorem 3.1 that the ruin probability $\psi_{\mu}$ of a fractional Poisson risk model satisfies the following fractional integro-differential equation

[TABLE]

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 has been derived in biard2014fractional . 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\sim\mathrm{ML}(\mu,\lambda_{2})$ and exponentially distributed claim sizes with parameter $\alpha$, the ruin probability equals

[TABLE]

where $x_{2}$ is the unique positive solution of $c^{\mu}x-\alpha c^{\mu}+\lambda_{2}x^{1-\mu}=0.$

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

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 Equation (28)), is an increasing function of both $\lambda_{2}$ and $\mu$. Figure 2b shows the values of 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 $\prescript{\mathrm{C}}{u}{\mathrm{D}}_{\infty}^{\mu}$ tends to the identity operator, when $\mu\rightarrow 0+$. Thus, we obtain the following result.

Corollary 2

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 an integral equation

[TABLE]

with the universal boundary condition $\lim_{u\rightarrow\infty}\psi_{0}(u)=0$.

Substituting $u=0$ into Equation (25) gives $\psi_{0}(0)=\frac{\lambda_{2}}{\lambda_{2}+1}$, which only depends on the value of $\lambda_{2}$. Taking Laplace transform both sides with respect to $u$ leads to

[TABLE]

which can be explicitly inverted back in some cases.

Appendix 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 podlubny1998fractional . 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 mittag1903nouvelle as a generalization of the exponential function.

Definition 2

The two-parameter Mittag-Leffler function is defined as

[TABLE]

Proposition 3

The Laplace transform of $z^{\alpha k+\beta-1}E^{(k)}_{\alpha,\beta}(\pm az^{\alpha})$ is (see Equation (1.80) in podlubny1998fractional )

[TABLE]

A.2 Fractional Integrals and Derivatives

As per hilfer2008threefold , we define and denote:

Definition 3

The left Riemann-Liouville fractional integral of order $r>0$ with lower limit $a\in\mathbb{R}$ is defined on locally integrable functions $f$ as

[TABLE]

and the right Riemann-Liouville fractional integral of order $r>0$ with upper limit $b\in\mathbb{R}$ is defined as

[TABLE]

Definition 4

The left Riemann-Liouville fractional derivative of order $r>0$ with lower limit $a$ is defined as the integer order derivatives of fractional integrals as follows,

[TABLE]

where $n=\lfloor r\rfloor+1$, and $\lfloor r\rfloor$ denotes the floor function. Similarly, the right Riemann-Liouville fractional derivative of order $r>0$ with upper limit $b$ is defined as

[TABLE]

These two operators are well defined on the Lebesgue space $L^{\lceil r\rceil}([a,b])$ (see Definition 2.1 in samko1993fractional ). Here, $\lceil r\rceil$ denotes the ceiling function.

Proposition 4

The Riemann-Liouville fractional derivatives are the left inverse operators of the corresponding fractional integrals (see Section 3.2 in valerio2013fractional )

[TABLE]

Proposition 5

The left Riemann-Liouville fractional integrals $\prescript{}{a}{\mathrm{I}}^{r}_{x}$ and left fractional derivative $\prescript{}{a}{\mathrm{D}}^{r}_{x}$ of the power function $(x-a)^{p}$ are (see Equation (2.117) in podlubny1998fractional )

[TABLE]

Proposition 6

The left fractional derivative $\prescript{}{0}{\mathrm{D}}^{r}_{x}$ of the two-parameter Mittag-Leffler functions satisfies (see Equation (1.82) in podlubny1998fractional )

[TABLE]

Proposition 7

The left Riemann-Liouville fractional derivatives $\prescript{}{0}{\mathrm{D}}^{r}_{x}$ of an integral depending on a parameter $t\in\mathbb{R}$ is given by (see Equation (2.212) in podlubny1998fractional )

[TABLE]

Proposition 8

The left Riemann-Liouville fractional derivatives $\prescript{}{0}{\mathrm{D}}^{r}_{x}$ of the (positive density) convolution integral equals to (see Equation (2.213) in podlubny1998fractional )

[TABLE]

Definition 5

The Weyl-Liouville fractional derivatives samko1993fractional ; butzer2000introduction are special cases of the Riemann-Liouville derivatives, whenever $a$ is replaced by $-\infty$ or $b$ is replaced by $\infty$ in Definition 4. The right Weyl-Liouville fractional derivative is defined for functions $f\in L^{\lceil r\rceil}([a,b])$ as

[TABLE]

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\in L^{\lceil r\rceil}([a,b])$ as

[TABLE]

Proposition 9

The Caputo fractional derivatives are the left inverse operators of their corresponding fractional integrals (see Section 3.2 in valerio2013fractional )

[TABLE]

Proposition 10

The Caputo and left Riemann-Liouville fractional derivatives are related by the following integration by parts formula (see Section 2.1 in almeida2011necessary )

[TABLE]

Proposition 11

The eigenfunction of left fractional derivative $\prescript{}{0}{\mathrm{D}}^{r}_{x}$ (or $\prescript{\mathrm{C}}{0}{\mathrm{D}}^{r}_{x}$) is $x^{1-\alpha}E_{\alpha,\alpha}(\lambda x^{\alpha})$ with eigenvalue $\lambda\in\mathbb{R}$ (see Section 2.2.3 in hilfer2008threefold ).

Proposition 12

The eigenfunction of right fractional derivative $\prescript{}{x}{\mathrm{D}}^{r}_{\infty}$ (or $\prescript{\mathrm{C}}{x}{\mathrm{D}}^{r}_{\infty}$) is $e^{-\lambda x}$ with eigenvalue $\lambda^{r}$, where $\lambda\in\mathbb{R}^{+}$ (see Section 4 in valerio2013fractional ).

Proposition 13

The Laplace transform of the left Riemann-Liouville fractional derivative of order $r>0$ is (see Equation (2.245) in podlubny1998fractional )

[TABLE]

Appendix B Proof of Theorem 2.1

Proof

We will use induction on both variables to validate (11) together with the extra statement: for any function $g$ supported on $[0,\infty)$, $\mathcal{A}_{m,n}\left(\frac{\mathrm{d}}{\mathrm{d}t}\right)[f_{T}^{m,n}*g](t)=\Lambda_{m,n}g(t)$. Base step: when $m=1,n=0$ or $m=0,n=1$, from Equation (8) and (9) we have $\mathcal{A}_{1,0}(\frac{\mathrm{d}}{\mathrm{d}t})[f_{T}^{1,0}](t)=0$ and $\mathcal{A}_{0,1}(\frac{\mathrm{d}}{\mathrm{d}t})[f_{T}^{0,1}](t)=0.$ Furthermore, a simple calculation yields

[TABLE]

Inductive step: for a non-negative $m$ and $n$, we assume that the statements

[TABLE]

hold. We then compute,

[TABLE]

thereby showing $m+1$ and $n+1$ cases are true. To validate the boundary conditions, we compute

[TABLE]

which equals to $\Lambda_{m,n}$ when $k=1$, and [math] for $k>1$. This completes the proof.

Appendix C Review of Fractional Poisson Process

The fractional Poisson process, denoted by $N_{\mu}(t)$, $t>0$, $\mu\in(0,1]$, is a fractional non-Markovian generalisation of Poisson process $N(t)$, $t>0$. The distribution of fractional Poisson process $P_{\mu}(n,t)=\mathbb{P}[N_{\mu}(t)=n]$ is defined by solving a fractional generalisation of the Kolmogorov-Feller equation laskin2003fractional

[TABLE]

where $\lambda$ is the intensity parameter and $\delta_{n,0}$ is the Kronecker symbol.

Moreover, laskin2003fractional showed the inter-arrival times of a fractional Poisson process have probability density function $f_{\mu}(t)=\lambda t^{\mu-1}E_{\mu,\mu}(-\lambda t^{\mu}),$ $t>0.$ The Laplace transform of the inter-arrival time density $f_{\mu}(t)$ is $\mathcal{L}\left\{f_{\mu}(t);s\right\}=\hat{f}_{\mu}(s)=\frac{\lambda}{s^{\mu}+\lambda}.$ The mean and variance of $N_{\mu}(t)$ are

[TABLE]

respectively ${\mathbb{V}\mathrm{ar}}\,N_{\mu}(t)=\frac{2(\lambda t^{\mu})^{2}}{\Gamma(2\mu+1)}-\frac{(\lambda t^{\mu})^{2}}{\left[\Gamma(\mu+1)\right]^{2}}+\frac{\lambda t^{\mu}}{\Gamma(\mu+1)},$ as in laskin2003fractional .