Skellam Type Processes of Order k and Beyond

In this article, we introduce the Skellam process of order k and its running average. We also discuss the time-changed Skellam process of order k. In particular, we discuss the space-fractional Skellam process and tempered space-fractional Skellam process via time changes in Skellam process by independent stable subordinator and tempered stable subordinator, respectively. We derive the marginal probabilities, Lévy measures, governing difference-differential equations of the introduced processes. Our results generalize the Skellam process and running average of Poisson process in several directions.


Introduction
The Skellam distribution is obtained by taking the difference between two independent Poisson distributed random variables, which was introduced for the case of different intensities λ 1 , λ 2 by (see [1]) and for equal means in [2]. For large values of λ 1 + λ 2 , the distribution can be approximated by the normal distribution and if λ 2 is very close to 0, then the distribution tends to a Poisson distribution with intensity λ 1 . Similarly, if λ 1 tends to 0, the distribution tends to a Poisson distribution with non-positive integer values. The Skellam random variable is infinitely divisible, since it is the difference of two infinitely divisible random variables (see Proposition 2.1 in [3]). Therefore, one can define a continuous time Lévy process for Skellam distribution, which is called Skellam process.
The Skellam process is an integer valued Lévy process and it can also be obtained by taking the difference of two independent Poisson processes. Its marginal probability mass function (pmf) involves the modified Bessel function of the first kind. The Skellam process has various applications in different areas, such as to model the intensity difference of pixels in cameras (see [4]) and for modeling the difference of the number of goals of two competing teams in a football game [5]. The model based on the difference of two point processes is proposed in (see [6][7][8][9]).
Recently, the time-fractional Skellam process has been studied in [10], which is obtained by time-changing the Skellam process with an inverse stable subordinator. Further, they provided the application of time-fractional Skellam process in modeling of arrivals of jumps in high frequency trading data. It is shown that the inter arrival times between the positive and negative jumps follow Mittag-Leffler distribution rather then the exponential distribution. Similar observations are observed in the case of Danish fire insurance data (see [11]). Buchak and Sakhno, in [12], have also proposed the governing equations for time-fractional Skellam processes. Recently, [13] introduced time-changed Poisson process of order k, which is obtained by time changing the Poisson process of order k (see [14]) by general subordinators.
In this paper, we introduce Skellam process of order k and its running average. We also discuss the time-changed Skellam process of order k. In particular, we discuss space-fractional Skellam process

Preliminaries
In this section, we collect relevant definitions and some results on Skellam process, subordinators, space-fractional Poisson process, and tempered space-fractional Poisson process. These results will be used to define the space-fractional Skellam processes and tempered space-fractional Skellam processes.

Skellam Process
In this section, we revisit the Skellam process and also provide a characterization of it. Let S(t) be a Skellam process, such that where N 1 (t) and N 2 (t) are two independent homogeneous Poisson processes with intensity λ 1 > 0 and λ 2 > 0, respectively. The Skellam process is defined in [8] and the distribution has been introduced and studied in [1], see also [2]. This process is only symmetric when λ 1 = λ 2 . The pmf s k (t) = P(S(t) = k) of S(t) is given by (see e.g., [1,10]) where I k is modified Bessel function of first kind (see [19], p. 375), The pmf s k (t) satisfies the following differential difference equation (see [10]) with initial conditions s 0 (0) = 1 and s k (0) = 0, k = 0. For a real-valued Lévy process Z(t) the characteristic function admits the form where the function ψ Z is called characteristic exponent and it admits the following Lévy-Khintchine representation (see [20]) Here, a ∈ R, b ≥ 0 and π Z is a Lévy measure. If π Z (dx) = ν Z (x)dx for some function ν Z , then ν Z is called the Lévy density of the process Z. The Skellam process is a Lévy process, its Lévy density ν S is a linear combination of two Dirac delta functions, ν S (y) = λ 1 δ 1 (y) + λ 2 δ −1 (y) and the corresponding characteristic exponent is given by The moment generating function (mgf) of Skellam process is With the help of mgf, one can easily find the moments of Skellam process. In the next result, we give a characterization of Skellam process, which is not available in literature as per our knowledge. For a function h, we write h(δ) = o(δ) if lim δ→0 h(δ)/δ = 0. Theorem 1. Suppose that an arrival process has the independent and stationary increments and it also satisfies the following incremental condition, then the process is Skellam.
m > n, m = n + 1; Proof. Consider the interval [0,t], which is discretized with n sub-intervals of size δ each, such that nδ = t. For k ≥ 0, we have by taking n → ∞. The result follows now by using the definition of modified Bessel function of first kind I k . Similarly, we prove when k < 0.

Poisson Process of Order k (PPoK)
In this section, we recall the definition and some important properties of Poisson process of order k (PPoK). Kostadinova and Minkova (see [14]) introduced and studied the PPoK. Let x 1 , x 2 , · · · , x k be non-negative integers and Additionally, let {N k (t)} t≥0 , represent the PPoK with rate parameter λt, then probability mass function (pmf) is given by The pmf of N k (t) satisfies the following differential-difference equations (see [14]) with initial condition p N k 0 (0) = 1 and p N k n (0) = 0 and n ∧ k = min{k, n}. The characteristic function of where i = √ −1. The process PPoK is Lévy, so it is infinite divisible i.e. φ N k (t) (u) = (φ N k (1) (u)) t . The Lévy density for PPoK is easy to derive and it is given by where δ j is the Dirac delta function concentrated at j. The transition probability of the PPoK {N k (t)} t≥0 is also given by Kostadinova and Minkova [14], The probability generating function (pgf) G N k (s, t) is given by (see [14]) The mean, variance and covariance function of the PPoK are given by

Subordinators
Let D f (t) be real valued Lévy process with non-decreasing sample paths and its Laplace transform has the form is the integral representation of Bernstein functions (see [21]). The Bernstein functions are C ∞ , non-negative and such that (−1) m d m dx m f (x) ≤ 0 for m ≥ 1 in [21]. Here, π denote the non-negative Lévy measure on the positive half line, such that ∞ 0 (x ∧ 1)π(dx) < ∞, π([0, ∞)) = ∞, and b is the drift coefficient. The right continuous inverse E f (t) = inf{u ≥ 0 : D f (u) > t} is the inverse and first exist time of D f (t), which is non-Markovian with non-stationary and non-independent increments. Next, we analyze some special cases of Lévy subordinators with drift coefficient b = 0, which is, It is worth noting that, among the subordinators given in (14), all of the integer order moments of α-stable subordinators are infinite and others subordinators have all finite moments.

The Space-Fractional Poisson Process
In this section, we discuss main properties of space-fractional Poisson process (SFPP). We also provide the Lévy density for SFPP, which is not discussed in the literature. The SFPP N α (t) was introduced by (see [22]), as follows where D α (t) is an α-stable subordinator, which is independent of homogeneous Poisson process N(t).
The probability generating function (pgf) of this process is of the form The pmf of SFPP is where h ψ i (z) is the Fox Wright function (see formula (1.11.14) in [23]). It was shown in [22] that the pmf of the SFPP satisfies the following fractional differential-difference equations with initial conditions where δ k,0 is the Kronecker delta function, given by The fractional difference operator is defined in [24], where B is the backward shift operator. The characteristic function of SFPP is Proposition 1. The Lévy density ν N α (x) of SFPP is given by Proof. We use Lévy-Khintchine formula (see [20]), which is the characteristic exponent of SFPP from Equation (23).

Tempered Space-Fractional Poisson Process
The tempered space-fractional Poisson process (TSFPP) can be obtained by subordinating homogeneous Poisson process N(t) with the independent tempered stable subordinator D α,µ (t) (see [25]) This process have finite integer order moments due to the tempered α-stable subordinator. The pmf of TSFPP is given by (see [25]) The governing difference-differential equation is given by The characteristic function of TSFPP, While using a standard conditioning argument, the mean and variance of TSFPP are given by Proof. Using (28), the characteristic exponent of TSFPP is given by . We find the Lévy density with the help of Lévy-Khintchine formula (see [20]), hence proved.

Definition 1.
A stochastic process X(t) is over-dispersed, equi-dispersed or under-dispersed [18], if the Fisher index of dispersion, given by (see e.g., [17]) is more than 1, equal to 1, or smaller than 1, respectively, for all t > 0.

Running Average of PPoK
In this section, we first introduced the running average of PPoK and their main properties. These results will be used further to discuss the running average of SPoK.
Definition 2 (Running average of PPoK). We define the running average process N k A (t), t ≥ 0 by taking time-scaled integral of the path of the PPoK (see [26]), We can write the differential equation with initial condition N k Which shows that it has continuous sample paths of bounded total variation. We explored the compound Poisson representation and distribution properties of running average of PPoK. The characteristic of N k A (t) is obtained using the Lemma 1 of [26]. We recall Lemma 1 from [26] for ease of reference.

Lemma 1. If X t is a Lévy process and Y t its Riemann integral is defined by
Proposition 3. The characteristic function of N k A (t) is given by Proof. Using the Equation (10), we have Using (32) and (31), we have Proposition 4. The running average process has a compound Poisson representation, such that where X i = 1, 2, . . . are independent, identically distributed (iid) copies of X random variables, independent of N(t) and N(t) is a Poisson process with intensity kλ. Subsequently, Further, the random variable X has following pdf where V i follows discrete uniform distribution over (0, k) and U i follows continuous uniform distribution over (0, i), i = 1, 2, . . . , k. (45), the characteristic function of X is given by which is equal to the characteristic function of PPoK that is given in (33). Hence, by the uniqueness of characteristic function, the result follows.
Using the definition the first two moments for random variable X given in Proposition (4) Corollary 1. Putting k = 1, the running average of PPoK N k A (t) reduces to the running average of standard Poisson process N A (t) (see Appendix in [26]). (2k + 1). If k = 1 the process is under-dispersed and for k > 1 it is over-dispersed.
Next we discuss the long-range dependence (LRD) property of running average of PPoK. We recall the definition of LRD for a non-stationary process.
Definition 3 (Long range dependence (LRD)). Let X(t) be a stochastic process that has a correlation function for s ≥ t for fixed s, that satisfies, We say that, if d ∈ (0, 1), then X(t) has the LRD property and if d ∈ (1, 2) it has short-range dependence (SRD) property [27]. Proof. Let 0 ≤ s < t < ∞, then the correlation function for running average of PPoK N k

Skellam Process of Order k (SPoK)
In this section, we introduce and study the Skellam process of order k (SPoK).
Proof. For m ≥ 0, using the pmf of PPoK that is given in (8), it follows using the multinomial theorem and modified Bessel function given in (2). Similarly, it follows for m < 0.

Proposition 7. The Lévy density for SPoK is
Proof. The proof follows by using the independence of two PPoK used in the definition of SPoK.

Remark 3.
Using (12), the pgf of SPoK is given by Further, the characteristic function of SPoK is given by

SPoK as a Pure Birth and Death Process
In this section, we provide the transition probabilities of SPoK at time t + δ, given that we started at time t. Over such a short interval of length δ → 0, it is nearly impossible to observe more than k event; in fact, the probability to see more than k event is o(δ).

Proposition 8. The transition probabilities of SPoK are given by
m > n, m = n + i, i = 1, 2, . . . , k; Basically, at most k events can occur in a very small interval of time δ. Additionally, even though the probability for more than k event is non-zero, it is negligible.
Proof. Note that S k (t) = N k 1 (t) − N k 2 (t). We call N k 1 (t) as the first process and N k 2 (t) as the second process. For i = 1, 2, · · · , k, we have P(the first process has i+j arrivals and the second process has j arrivals) + P(the first process has i arrivals and the second process has 0 arrivals) + o(δ) Similarly, for i = 1, 2, · · · , k, we have P(the first process has j arrivals and the second process has i+j arrivals) + P(the first process has 0 arrivals and the second process has i arrivals) + o(δ) Further, P(S k (t + δ) = n|S k (t) = n) = k ∑ j=1 P(the first process has j arrivals and the second process has j arrivals) + P(the first process has 0 arrivals and the second process has 0 arrivals) + o(δ)

Remark 4.
The pmf R m (t) of SPoK satisfies the following difference differential equation with initial condition R 0 (0) = 1 and R m (0) = 0 for m = 0. Let B be the backward shift operator defined in (22) and F be the forward shift operator defined by F j X(t) = X(t + j), such that (1 − F) α = ∑ ∞ j=0 ( α j )F j . Multiplying by s m and summing for all m in (42), we obtain the following differential equation for the pgf The mean, variance and covariance of SPoK can be easily calculated by using the pgf, Remark 5. For the SPoK, when ] > 1 and hence S k (t) exhibits over-dispersion. For λ 1 < λ 2 , the process is under-dispersed.
Next, we show the LRD property for SPoK. Proposition 9. The SPoK has LRD property defined in Definition 3.

Proof. The correlation function of SPoK satisfies
Hence, SPoK exhibits the LRD property.

Running Average of SPoK
In this section, we introduce and study the new stochastic Lévy process, which is the running average of SPoK.

Definition 5.
The following stochastic process defined by taking the time-scaled integral of the path of the SPoK, is called the running average of SPoK.
Next, we provide the compound Poisson representation of running average of SPoK.
Proof. By using the Lemma 3.1 to Equation (40) after scaling by 1/t.

Remark 6.
It is easily observable that Equation (43) has removable singularity at u = 0. To remove that singularity, we can define φ S k A (t) (0) = 1.

Proposition 11. Let Y(t) be a compound Poisson process
where N(t) is a Poisson process with rate parameter k(λ 1 + λ 2 ) > 0 and {J n } n≥1 are iid random variables with mixed double uniform distribution function p j , which are independent of N(t). Subsequently, The random variable J 1 being a mixed double uniformly distributed has density where V j follows discrete uniform distribution over (0, k) with pmf p V j (x) = P(V j = x) = 1 k , j = 1, 2, . . . k, and U i be doubly uniform distributed random variables with density Further, 0 < w < 1 is a weight parameter and 1(·) is the indicator function. Here, we obtained the characteristic of J 1 using the Fourier transform of (45), putting the characteristic function φ J 1 (u) in the above expression yields the characteristic function of S k A (t), which completes the proof.

Time-Changed Skellam Process of Order k
We consider time-changed SPoK, which can be obtained by subordinating SPoK S k (t) with the independent Lévy subordinator Note that the stable subordinator does not satisfy the condition E[D f (t)] c < ∞. The mgf of time-changed SPoK Z f (t) is given by .

Theorem 2.
The pmf H f (t) = P(Z f (t) = m) of time-changed SPoK is given by Proof. Let h f (x, t) be the probability density function of Lévy subordinator. Using conditional argument The mean and covariance of time changed SPoK are given by,

Space Fractional Skellam Process and Tempered Space Fractional Skellam Process
In this section, we introduce time-changed Skellam processes where time-change are stable subordinator and tempered stable subordinator. These processes give the space-fractional version of the Skellam process similar to the time-fractional version of the Skellam process introduced in [10].
Next, we derive the mgf of SFSP. We use the expression for marginal (pmf) of SFPP that is given in (17) to obtain the marginal pmf of SFSP.
In the next result, we obtain the state probabilities of the SFSP. Theorem 3. The pmf H k (t) = P(S α 1 ,α 2 (t) = k) of SFSP is given by for k ∈ Z.
In the next theorem, we discuss the governing differential-difference equation of the marginal pmf of SFSP. Theorem 4. The marginal distribution H k (t) = P(S α 1 ,α 2 (t) = k) of SFSP satisfies the following differential difference equations with initial conditions H 0 (0) = 1 and H k (0) = 0 for k = 0.
Proof. The proof follows by using pgf.

Tempered Space-Fractional Skellam Process (TSFSP)
In this section, we present the tempered space-fractional Skellam process (TSFSP). We discuss the corresponding fractional difference-differential equations, marginal pmfs, and moments of this process. Definition 7 (TSFSP). The TSFSP is obtained by taking the difference of two independent tempered space fractional Poisson processes. Let D α 1 ,µ 1 (t), D α 2 ,µ 2 (t) be two independent TSS (see [28]) and N 1 (t), N 2 (t) be two independent Poisson processes that are independent of TSS. Subsequently, the stochastic process S µ1,µ2 is called the TSFSP.

Remark 10.
We use this expression to calculate the marginal distribution of TSFSP. The mgf is obtained using the conditioning argument. Let f α,µ (x, t) be the density function of D α,µ (t). Subsequently, Using (54), the mgf of TSFSP is

Simulation of SFSP and TSFSP
We present the algorithm to simulate the sample trajectories for SFSP and TSFSP. We use Python 3.7 and its libraries Numpy and Matplotlib for the simulation purpose. It is worth mentioning that Python is an open source and freely available software.
Use the first two steps of previous algorithm for generating the increments of α-stable subordinator D α 1 (t).
For α 1 = 0.6 and α 2 = 0.9 and fixed t, it is more probable that the value of the rv D α 1 (t) is higher than the rv D α 2 (t). Thus, for same intensity parameter λ for Poisson process the process N(D α 1 (t)) will have generally more arrivals than the process N(D α 2 (t)) until time t. This is evident from the trajectories of the SFSP in Figure 1, because the trajectories are biased towards positive side. The TSFPP is a finite mean process, however SFPP is an infinite mean process and hence SFSP paths are expected to have large jumps, since there could be a large number of arrivals in any interval.