Fractional Poisson Processes and Related Planar Random Motions

We present three different fractional versions of the Poisson process and some related results concerning the distribution of order statistics and the compound Poisson process. The main version is constructed by considering the difference-differential equation governing the distribution of the standard Poisson process, N ( t ) , t > 0, and by replacing the time-derivative with the fractional Dzerbayshan-Caputo derivative of order ν ∈ ( 0,1 ] . For this process, denoted by N ν ( t ) , t > 0, we obtain an interesting probabilistic representation in terms of a composition of the standard Poisson process with a random time, of the form N ν ( t ) = N ( T 2 ν ( t )) , t > 0. The time argument T 2 ν ( t ) , t > 0, is itself a random process whose distribution is related to the fractional diffusion equation. We also construct a planar random motion described by a particle moving at ﬁnite velocity and changing direction at times spaced by the fractional Poisson process N ν . For this model we obtain the distributions of the random vector representing the position at time t , under the condition of a ﬁxed number of events and in the unconditional case. For some speciﬁc values of ν ∈ ( 0,1 ] we show that the random position has a Brownian behavior (for ν = 1 / 2) or a cylindrical-wave structure (for ν = 1).


Introduction
Attempts to construct fractional versions of the Poisson process have been undertaken by Repin and Saichev (2000), Jumarie (2001) and Laskin (2003).
In analogy with the well-known fractional Brownian motion, the idea of passing from the usual Poisson process to a possible fractional version has inspired the papers by Wang et al. (2006Wang et al. ( )-(2007. The process proposed by them is constructed as a stochastic integral with respect to the Poisson measure.
We follow here a completely different approach and we construct various types of fractional Poisson processes.
We note that, for ν = 1, (1.1) coincides with the equation governing the homogeneous Poisson process, so that our results generalize the well-known distributions holding in the standard case.
The solution to (1.1)-(1.2) is expressed as the distribution of a process, which we will denote by ν (t), t > 0. The presence of the fractional derivative in (1.1) implies a series of consequences, the most important being that ν has dependent increments.
We will obtain various forms of the solution to (1.1): the first one is expressed as for k ≥ 0, t > 0.
The distribution (1.4) coincides with formula (25) of Laskin (2003), despite the latter being obtained as a solution to equation (1.1) with the fractional derivative defined in the Riemann-Liouville sense (instead of (1.3)). The additional impulse term appearing in equation (25) of Laskin (2003) justifies this coincidence.
In Jumarie (2001) the fractional derivative is also meant in the Riemann-Liouville sense, but it is subject to a fractional initial condition, thus obtaining again formula (1.4). For a general reference on fractional calculus see Samko et al. (1993).
An alternative form of the distribution is given in terms of the generalized Mittag-Leffler function, which is defined (see Podlubny (1999), p.17) as , ν, µ > 0, x ∈ . (1.5) For small values of k we have meaningful expression: for example, for k = 0, we get that Pr ν (t) = 0 = E ν,1 (−λt ν ), (1.6) which represents the distribution of the waiting time of the first fractional Poisson event.
The relationship (1.7) means that the fractional Poisson process ν has the following representation ν (t) = N ( 2ν (t)) (1.9) and thus it can be considered as a homogeneous Poisson process stopped at a random time 2ν (t).
On the basis of the previous results we construct, in Section 3, a planar random motion described by a particle moving at finite velocity c and changing direction at times spaced by the fractional Poisson process ν . At each change of direction (at fractional Poisson instants) the new direction is chosen with uniform law in [0, 2π] . This has been our first motivation of the research of this paper.
In the standard case, when the process governing the changes of directions is the standard Poisson process, this model has been studied in Kolesnik and Orsingher (2005). The conditional distribution of the random vector (X (t), Y (t)) , t > 0 (representing the position of the moving particle at time t) is given by for t > 0, x, y ∈ C c t = x, y : x 2 + y 2 ≤ c 2 t 2 .
We obtain here the analogue to (1.10) for the fractional case where ν ∈ (0, 1]. In particular, for the case ν = 1/2, the fractional process governing the changes of directions is 1/2 (t) = N (|B(t)|), since, in this case, equation (1.8) reduces to the heat-equation. We prove that the conditional distribution of the random vector (X (|B(t)|), Y (|B(t)|)) , t > 0 can be written as where B α, β denotes a Beta function of parameters α, β and B(t), t > 0, is a standard Brownian motion, independent of (X (t), Y (t)). This means that the planar motion with a Brownian time can be regarded as a planar Brownian motion, whose volatility is itself random and possesses a Beta distribution depending on the number of changes of directions.
The unconditional version of the distribution of the planar motion takes interesting forms, especially in some particular cases. In general, for any ν ∈ (0, 1], we have that where u 2ν = u 2ν (x, y, t; λ 2 ) is the fundamental solution to the fractional planar wave equation (see (3.16) below).
In the special case of ν = 1, the previous expression reduces to the following very interesting form: (1.14) The distribution (1.14) is the marginal of the unconditional distribution of the planar random motion (X (t), Y (t)) , t > 0 (performed at velocity 1), given in formula (1.3) of Orsingher and De Gregorio (2007).
The kernel function in (1.13) emerges in the study of cylindrical waves, that is the analysis of the following problem  (1.15) Section 4 is devoted to the investigation of alternative forms of fractional Poisson processes. The first one, denoted by N ν (t), is defined by generalizing to the case ν ∈ (0, 1] the Poisson distribution as follows so that the probability generating function reads As we will see below, for ν = 1, formulae (1.16) and (1.17) reduce to the corresponding results holding for the standard Poisson case. Also this version of fractional Poisson process does not possess the property of independence of increments.
We prove that the probability generating function (1.17) is a solution to the fractional equation Based on the model above, we construct a fractional counterpart of the Poisson compound process, defined (for any fixed t) as the sum of a random number N ν (t) of i.i.d. r.v.'s and we derive its characteristic function.
Moreover, for the first two versions of fractional Poisson processes, we obtain some results concerning the maximum, minimum and order statistics of flows of i.i.d. r.v.'s paced by ν and N ν .
In hydrology and seismology the sequence of catastrophic events can be successively represented by a fractional Poisson process (which is devoid of the lack of memory property of the usual Poisson process). The evaluation of the distribution of the maximum of i.i.d. r.v.'s, whose number is represented by a fractional Poisson process, is therefore of practical use in determining the probability of extremal events.

A first form of the fractional Poisson process
We construct the first type of fractional Poisson process (which we will denote by ν (t), t > 0) by replacing, in the differential equations governing the distribution of the classical Poisson process, the standard derivatives by the fractional derivatives, defined in (1.3), i.e. in the Dzerbayshan-Caputo sense.
Therefore we are interested in the solution of The use of definition (1.3) permits us to avoid fractional derivatives in the initial conditions and to obtain different and more explicit expressions of the solution, from a probabilistic point of view.
The equation governing the probability generating function G in the standard case, i.e.
is then replaced by the following equation subject to the initial condition G(u, 0) = 1.
The Laplace transform of the solution to (2.4) becomes The inverse Laplace transform of (2.5) yields the following probability generating function The uniqueness of the solutions to (2.4) can be proved by taking two solutions (G 1 and G 2 ) and by solving the Cauchy problem for h(u, t) = G 1 (u, t) − G 2 (u, t), with initial condition h(u, 0) = 0. The existence of the solution is a consequence of our analysis.
Much information can be extracted from (2.6): some corollaries are immediate, for example Var ν (t) = λt ν Γ(ν + 1) We note that the usual equality between the mean value and the variance does not hold for this model, while, for ν = 1, it can be checked that 1 (t) = Var 1 (t). The most interesting result is the distribution function related to (2.6), which can be presented as follows For small values of k it is possible to write the distribution (2.10) in terms of Mittag-Leffler functions as follows Pr We derive now an explicit representation of the solution to (2.1)-(2.2), which suggests an interesting probabilistic interpretation of the related process. The next result represents also an important link between the fractional Poisson process and the fractional diffusion equations.
It is well-known that v 2ν ( y, t) = 1 where be the folded solution to (2.12), then the probability generating function of the fractional Poisson process ν (t), t > 0 can be written as The corresponding distribution reads, for k ≥ 0, where 2ν (t), t > 0 represents a random time with transition density given in (2.15).
Proof The Laplace transform (2.5) can be written as We note that where f ν (·; y) is a stable law S ν (µ, β, σ) of order ν, with parameters µ = 0, β = 1 and σ = y cos πν Therefore, by using together (2.18) and (2.19), we obtain the inverse Laplace transform as (2.20) It has been proved in Orsingher and Beghin (2004) that the solution v 2ν to (2.12) can be expressed as We note that, in particular for ν = 1/2, the random time 1 (t), t > 0 becomes a reflecting Brownian motion, since, in this case, equation (2.12) reduces to the heat-equation and the solution v 1 ( y, t) is the density of a Brownian motion B(t), t > 0 with infinitesimal variance 2. After folding up the solution, we find the following density so that 1/2 coincides with a standard Poisson process at a Brownian time.

Remark 2.2
We can check that, by replacing (2.14) into (2.16) we retrieve (2.6). We apply the following representations in terms of integrals on the Hankel contour H a of the Wright and Mittag-Leffler functions, respectively: For a plot of the Hankel contour see Fig.1. Therefore we achieve

Remark 2.3
We derive from (2.17) an alternative expression of the distribution of the fractional Poisson process (in view of (2.14)), as follows: Pr where the dependence on k is limited to the binomial coefficient. It is easy to check in this form that the distribution (2.25) sums up to one with respect to k, as follows: Pr where, in the last step, we applied the contour-integral representation of the reciprocal of the Gamma function, i.e. 1 The expression (2.25) suggests another useful representation of the solution in terms of derivatives of Mittag-Leffler functions, which generalizes the particular cases obtained in (2.11): Pr where, in the last step we have applied (5.1) in Appendix.
For k = 0, 1, 2, 3 we obtain again formulae (2.11): for example, for k = 2 and, by considering result (5.2) in Appendix, we have that The case k = 3 can be derived by performing similar steps and applying (5.1) and (5.3).

Remark 2.4
We can evaluate the Laplace transform of (2.17) as follows A further check stems from the fact that which coincides with (2.5).

Remark 2.5
We note that the fractional Poisson process ν loses the lack of memory property of the homogeneous Poisson process: indeed it is easy to check that it possesses dependent increments, since, from (2.6), we get We consider now a sequence of a random number of i.i.d. non-negative random variables X i , i ≥ 1 with distribution function F (w) = Pr X i < w , i ≥ 1, under the assumption that this random number is represented by ν (t). We will obtain the distribution of the k-th order statistic X ν (t) (k) , k ≥ 1, conditioned on the event ν (t) ≥ k.
, t > 0 be a fractional Poisson process of parameter λ, then the conditional Proof We recall that in the standard case the k-th order statistics possess the following distribution, therefore in this case we have that We can evaluate the numerator as By inserting (2.31) into (2.30) we get (2.28), in view of (2.10).

Remark 2.6
We note that for z → ∞ the previous expression converges to 1, as happens for the standard case (2.29).
As a particular case we can obtain the distribution of the minimum of a random number of i.i.d. random variables X i , i ≥ 1, by taking k = 1 in (2.28): Again, for z → ∞, (2.32) converges to 1.
As far as the distribution of the maximum is concerned we can evaluate it as follows: We can verify that for ν = 1, i.e. in the standard case, result (2.33) yields the well known formula: However if we consider that and this distribution has a jump of height e −λt in z = 0. If N (t) = 0 we assume that max 1≤ j≤N (t) X j = −∞.
In the fractional case we have analogously that The limit for z → ∞ is 1, while for z → 0 there is a jump equal to Pr ν (t) = 0 .

Planar random motions changing direction at fractional Poisson times
Let us consider a planar random motion described by a particle moving at finite velocity c and changing direction at Poisson time instants (for details on this process see Kolesnik and Orsingher (2005)). At each Poisson event, the new direction is chosen with uniform law in [0, 2π] . Let us denote the current position of the moving particle by (X (t), Y (t)) , t > 0. We recall the formulae (11) and (18) in Kolesnik and Orsingher (2005) that provide the conditional distribution of the planar random motion for t > 0, x, y ∈ C c t = x, y : x 2 + y 2 ≤ c t and the corresponding characteristic function where J k 2 is the Bessel function of order k/2.
We will examine here the same planar motion (taken at a Brownian time) when the process governing the changes of direction is replaced by a fractional Poisson process. In particular, we will refer to the process ν (t), t > 0 analyzed in Section 2.
We start by considering the case where ν = 1 2 , since it permits us to exploit the representation (2.23). Indeed in this case, in view of Theorem 2.1, the fractional process 1 2 coincides in distribution with a Poisson process at the random time 1 . As we have already seen in Remark 2.1, 1 reduces to the reflecting Brownian motion, i.e. 1 (t) = |B(t)|, t > 0, where again B(t), t > 0 is a Brownian motion with infinitesimal variance 2.
We present the distribution of (X (|B(t)|), Y (|B(t)|)), t > 0, conditional on a fixed number of events At time |B(t)| the set of possible positions of the process is represented by the whole circle with radius c|B(t)|, if k ≥ 1, and coincides with the circumference of the same circle, if k = 0. In both cases, since the radius is the random variable c|B(t)| distributed on the positive real line, the support covers the whole plane.
In the next theorem we will show that the distribution of the random vector under the condition of a number k of events possesses a surprisingly explicit form: it coincides with a bivariate Gaussian with variance (itself random) equal to 2c 2 tW k , where W k is a Beta random variable with parameters 1 2 , k+1 2 .

Theorem 3.1
Under the condition N1 2 (t) = k the conditional distribution of the vector Proof For k ≥ 1, in view of (3.1) we have that We now take the Fourier transform of (3.4) and apply (3.2), so that we get 2 e iαx+iβ y d x d y By inverting the characteristic function (3.5) we get formula (3.3), for k ≥ 1.
For k = 0, by performing similar steps we can check that result (3.3) is valid also in this case.
We can analogously consider the planar random motion changing direction at fractional Poisson times of order ν = 1 4 , by applying the results obtained in Orsingher and Beghin (2009). In this paper it is proved that the solution to the fractional diffusion equation of order 1 can be written as Formula (3.7) coincides with the density of the n-times iterated Brownian motion where the j 's are standard independent Brownian motions (i.e. with infinitesimal variance 1).
In order to apply Theorem 2.1, we must suitably adapt the previous results, taking into account that, , (3.9) and this differs from (3.6) for the coefficient representing the infinitesimal variance. Thus, instead of (3.7), we must use the solution to (3.9), which reads Therefore we can interpret the random time as the iterated reflected Brownian motion 2 (t) = |I 1 (t)| = |B 1 (|B 2 (t)|)|, where B 1 and B 2 are independent Brownian motions with infinitesimal variance 2. We study the random motion under the condition that the changes of directions are governed by a fractional Poisson process 1 4 (t) = N (|B 1 (|B 2 (t)|)|) = k, k ≥ 0. 1 4 (t) = k, the conditional distribution of the vector
Proof As in (3.4) we get, for k ≥ 1, By taking the Fourier transform of (3.12) we get so that we derive (3.11).
We can conclude that the random vector, under the condition of a number k of Poisson events, possesses a bivariate Gaussian distribution with variance (itself random) equal to 2c 2 |B(t)|W k , where W k is a Beta random variable with parameters 1 2 , k+1 2 and B(t) is a Brownian motion with infinitesimal variance 2.  = d x d y for k ≥ 0 and thus the process possesses a conditional distribution coinciding with the bivariate Gaussian law with variance equal to 2c 2 |B 1 (...|B n−1 (t)|...)|W k . It is easy to check that in the case ν = 1 4 , i.e. for n = 2, formula (3.14) reduces to (3.11).
We consider now the general case 0 < ν ≤ 1: the random process representing the time argument is 2ν (t), t > 0 and the condition is therefore equal to ν (t) = N ( 2ν (t)) = k.

Theorem 3.3 Under the condition ν (t) = k, the conditional characteristic function of the vector
Proof In this case we must evaluate the following integral (which is analogous to (3.5) and (3.13)) from which we derive (3.15).
We note that, for ν = 1 2 , formula (3.15) reduces to (3.5). Theorems 3.1, 3.2 and 3.3, as well as Remark 3.1, show that the assumption of random time for the planar motion at finite velocity (with direction changing at fractional Poisson times) implies a drastic change in its structure. This is transparent in (3.5) where the resulting motion is a planar Brownian motion, whose infinitesimal variance is 2c 2 W k and W k ∼ Bet a 1 2 , k+1 2 . Therefore the larger the number k of changes of direction, the more concentrated becomes the distribution of the random position. Clearly W k = 1 k+2 and this should convince the reader that the planar motion considered has sample paths that coil up around the starting point.
We now derive the law of the planar random motion in the general case 0 < ν ≤ 1. We will see that it can be expressed explicitly in terms of the planar process governed by the following twodimensional fractional differential equation: with the additional condition Also in this case the time argument of the solution u 2ν = u 2ν (x, y, t; λ 2 ) will be itself random and will be represented in terms of a r.v. W k ∼ Bet a 1 2 , k+1 2 .

Theorem 3.4 Under the condition
where u 2ν (x, y, t; c 2 w) is the solution to the fractional two-dimensional heat-wave equation (3.16), where λ 2 = c 2 w.
Proof The inverse Fourier transform of E 2ν,1 (−w t 2ν c 2 (α 2 + β 2 )) can be calculated as follows so that we get the first form of (3.18). In (3.19) we used the following integral representation of the Bessel function: We observe that (3.19) coincides with the solution u 2ν (x, y, t; λ 2 ) of the fractional equation (3.16), for λ 2 = c 2 w, and then the second form of (3.18) follows.
We can show that the two expressions in (3.21) coincide by evaluating the Fourier transform of the second line: The characteristic function (3.22) can be inverted as follows which coincides with the first line of (3.21).
We consider now the unconditional distribution of X ( 2ν (t)), Y ( 2ν (t)) , t > 0, that, for any ν ∈ (0, 1], is presented in the next result: , t > 0 has the following representation for ν ∈ (0, 1]. Proof In view of formulae (2.10) and (3.18) we have that The double sum appearing in (3.24) can be further treated as follows: By inserting (3.25) into (3.24) we have that and result (3.23) easily follows.

Remark 3.4
For ν = 1 formula (3.23) can be rewritten in a more suitable form as where u 2 (x, y, t; c 2 w) is given in (3.21).
The expression in square brackets can be rewritten as follows so that, by the change of variable m + 1 = k, (3.29) becomes It can be easily checked that (3.30) integrates to one, by considering that u 2 (x, y, t; c 2 w) is a transition density and then 2 u 2 (x, y, t; c 2 w)d x d y = 1. By inserting formula (3.21) into (3.30) we get We can recognize in (3.31) the following expression which represents, for |u| < t, the distribution at time t of the projection on the horizontal (or, equivalently, vertical) axis of the planar motion changing directions at Poisson times, performed at speed equal to 1 (see formula (1.3) in Orsingher and De Gregorio (2007), for c = 1). This motion can also be regarded as a telegraph process changing velocity, randomly, at Poisson times (see, for details, Stadje and Zacks (2004)).
Therefore we finally get The distribution (3.32) is concentrated in a circle of radius c t and therefore has a different nature than all the other distributions (3.23), which are positive-valued on the whole plane. In effect (3.32) represents a cylindrical wave expanding on a circle with random radius, whose distribution is related to a telegraph process with random velocity.

Alternative forms of fractional Poisson processes
We present now two different versions of fractional Poisson processes obtained by following an approach alternative to that used in Section 2.
The first process we construct here, denoted by N ν (t), t > 0, is obtained by replacing the factorial function in the Poisson distribution with the Gamma function and writing the probabilities as follows: where λ > 0, 0 < ν ≤ 1.
For ν = 1 the distribution (4.1) coincides with that of the homogeneous Poisson process. The process with one-dimensional distribution (4.1) can be viewed as a fractional version of the Poisson process because the Mittag-Leffler function emerges as the counterpart to the exponential function in the analysis of equations with fractional order derivatives. Nevertheless some basic properties of the Poisson process are lost when considering (4.1): in particular the independence of increments on non-overlapping intervals.
It is easily ascertained that where T ν 1 = inf z : N ν (z) = 1 and thus (4.1) loses the lack of memory property of the Poisson distribution.
We note that the distribution (4.1) can be expressed in terms of the k-th factorial moments of the first-type fractional Poisson process ν , given in (2.9), as follows: The probability generating function of (4.1) reads which can be compared with (2.6) of the first model.
We note that the function G ν (u ν , t) is a solution to the fractional equation with the initial condition G ν (u ν , 0) = 1.
Indeed we can check that The mean value of the distribution (4.1) is while the variance reads We note that, for ν = 1, the expressions (4.6) and (4.7) coincide and thus as happens in the standard Poisson case.
We will consider now a sequence of a random number of non-negative i.i.d. random variables with distribution function F (w) = Pr X i < w , i ≥ 1, under the assumption that this random number is represented by N ν (t). The distributions of the maximum and the minimum of this sequence can be written as follows Pr max X 1 , ..., X N ν (t) < w = E ν,1 (λt F (w)) E ν,1 (λt) (4.8) Pr min X 1 , ..., We can compare result (4.8) with the analogous distribution obtained for the fractional Poisson process of the first type, ν , in (2.34): we note that they display similar structures and coincide for ν = 1.
Moreover we define a fractional counterpart of the compound Poisson process as follows and we obtain explicitly the characteristic function that satisfies the following relationship The distribution of the waiting time of the first event T ν 1 is equal to For the waiting time of the k-th event T ν k we can write that We construct now the third version of the fractional Poisson process (which we will denote by ν (t), t > 0), by assuming that its interarrival times U j coincide in distribution with the waiting time of the first event of the first model, i.e. T ν 1 = inf t : ν (t) = 1 . Therefore we deduce the distribution of the r.v.'s representing the interarrival times, U j , j ≥ 1, from (2.11), as follows and then (in view of (5.2) in Appendix) for any j ≥ 1, λ > 0, t > 0 and 0 < ν ≤ 1.
By taking the Laplace transform of (4.14) we get where g(t; λ, m) is a Gamma density of parameters λ, m.
We can obtain the distribution of the new process ν (t), t > 0, from that of the interarrival times U j , as follows: Pr ν (t) = m (4.17) From the analysis presented above we see that the r.v.'s representing the increments between two successive occurrences, i.e. ν (T ν m ) − ν (T ν m−1 ), m ≥ 1, are independent. However, since, for arbitrary time instants t j−1 , t j , the r.v.'s ν (t j ) − ν (t j−1 ) are not independent, the process is not endowed with independent increments.
This coincides with the result of Remark 2.4 and thus the one-dimensional distributions of the first and third models coincide.
Another representation of the distribution of the fractional Poisson process is given in the next theorem.  Proof In order to write down explicitly the distribution (4.17) we recall that the factor e −µ ν t appearing in (4.16) is the Laplace transform of the stable law f ν (·; t) (see (2.19)), therefore (4.16) can be rewritten as  where, in the last step, we have used result (2.23).
If we now take the first derivative of (4.26) we obtain

Appendix
We note that the following identity for the Mittag-Leffler function hold: for 0 < ν < 1 and x = 0. In order to prove (5.1) we use the contour integral representation of the inverse of the Gamma function, so that we can write H a e t t −1+ν k+ν d t = 1 x 1 2πi H a e t t −1+ν 1 1 − t ν 2πi H a e t t −1+ν x − t ν d t = −E ν,1 (x) .