# Fractional Skellam process of order $k$

Kuldeep Kumar Kataria  Kuldeep Kumar Kataria, Department of Mathematics,
Indian Institute of Technology Bhilai, Raipur 492015, India.
<EMAIL_ADDRESS>and  Mostafizar Khandakar  Mostafizar Khandakar,
Department of Mathematics, Indian Institute of Technology Bhilai, Raipur
492015, India<EMAIL_ADDRESS>
(Date: March 16, 2021)

###### Abstract.

We introduce and study a fractional version of the Skellam process of order
$k$ by time-changing it with an independent inverse stable subordinator. We
call it the fractional Skellam process of order $k$ (FSPoK). An integral
representation for its one-dimensional distributions and their governing
system of fractional differential equations are obtained. We derive the
probability generating function, mean, variance and covariance of the FSPoK
which are utilized to establish its long-range dependence property. Later, we
considered two time-changed versions of the FSPoK. These are obtained by time-
changing the FSPoK by an independent Lévy subordinator and its inverse. Some
distributional properties and particular cases are discussed for these time-
changed processes.

###### Key words and phrases:

Skellam distribution; Poisson process of order $k$; Lévy subordinator; LRD
property.

###### 2010 Mathematics Subject Classification:

Primary: 60G22; Secondary: 60G55

##  1\. Introduction

Let $\\{N(t)\\}_{t\geq 0}$ be a Poisson process with intensity $k\lambda$ and
$\\{X_{i}\\}_{i\geq 1}$ be a sequence of independent and identically
distributed (iid) discrete uniform random variables with support
$S=\\{1,2,\dots,k\\}$. Consider the following compound Poisson process:

$N^{k}(t)=\sum_{i=1}^{N(t)}X_{i},$

where $\\{N(t)\\}_{t\geq 0}$ is independent of $\\{X_{i}\\}_{i\geq 1}$. The
process $\\{N^{k}(t)\\}_{t\geq 0}$ is known as the Poisson process of order
$k$ (PPoK) (see Kostadinova and Minkova (2013)). For $k=1$, the PPoK reduces
to the Poisson process $\\{N(t)\\}_{t\geq 0}$ with intensity $\lambda$. In
ruin theory, the PPoK is used to model the number of claims where these arrive
in groups of size $k$.

Recently, Gupta et al. (2020) introduced and studied a Lévy process
$\\{S^{k}(t)\\}_{t\geq 0}$ by considering the difference of two PPoK, that is,

$S^{k}(t)=N_{1}^{k}(t)-N_{2}^{k}(t),$ (1.1)

where $\\{N_{1}^{k}(t)\\}_{t\geq 0}$ and $\\{N_{2}^{k}(t)\\}_{t\geq 0}$ are
independent PPoK with intensities $\lambda_{1}>0$ and $\lambda_{2}>0$,
respectively. It is called the Skellam process of order $k$ (SPoK) whose state
probabilities $p^{k}(n,t)=\mathrm{Pr}\\{S^{k}(t)=n\\}$, $n\in\mathbb{Z}$
satisfy the following system of differential equations:

$\frac{\mathrm{d}}{\mathrm{d}t}p^{k}(n,t)=-k(\lambda_{1}+\lambda_{2})p^{k}(n,t)+\lambda_{1}\sum_{j=1}^{k}p^{k}(n-j,t)+\lambda_{2}\sum_{j=1}^{k}p^{k}(n+j,t),$
(1.2)

with initial conditions $p^{k}(0,0)=1$ and $p^{k}(n,0)=0,\ n\neq 0$. For
$k=1$, the SPoK reduces to an integer-valued Lévy process, namely, the Skellam
process $\\{S(t)\\}_{t\geq 0}$ (see Barndorff-Nielsen et al. (2012)). For each
$t\geq 0$, the random variable $S(t)$ has the Skellam distribution (see
Skellam (1946)). It has several real life applications, for example, it is
used as a sensor noise model for cameras (see Hwang et al. (2007)) and in
modeling the score differences in a game of two competing teams (see Karlis
and Ntzoufras (2009)).

In the past two decades, the time-changed point processes attracted the
interest of several researchers due to their potential applications in
different fields such as finance, hydrology, econometrics, etc. The Poisson
process time-changed by a stable subordinator and by the inverse of a stable
subordinator leads to the space fractional Poisson process (see Orsingher and
Polito (2012)) and the time fractional Poisson process (see Meerschaert et al.
(2011)), respectively. For other time-changed version of the Poisson process,
we refer the reader to Beghin (2012), Orsingher and Toaldo (2015), Aletti et
al. (2018), etc., and the references therein.

A subordinator $\\{D_{f}(t)\\}_{t\geq 0}$ is a one-dimensional Lévy process
whose Laplace transform is given by (see Applebaum (2009), Section 1.3.2)

$\mathbb{E}\left(e^{-sD_{f}(t)}\right)=e^{-tf(s)},$

where the function

$f(s)=bs+\int_{0}^{\infty}\left(1-e^{-sx}\right)\,\mu(\mathrm{d}x),\ \ s>0,$
(1.3)

is called the Bernštein function. Here, $b\geq 0$ is the drift coefficient and
$\mu$ is a non-negative Lévy measure that satisfies $\mu([0,\infty))=\infty$
and $\displaystyle\int_{0}^{\infty}(x\wedge 1)\,\mu(\mathrm{d}x)<\infty$ where
$s\wedge t=\min\\{s,t\\}$. The sample paths of a subordinator
$\\{D_{f}(t)\\}_{t\geq 0}$ are non-decreasing and $D_{f}(0)=0$ almost surely
(a.s.). The first hitting time of $\\{D_{f}(t)\\}_{t\geq 0}$ is called the
inverse subordinator. It is defined as

$H_{f}(t)\coloneqq\inf\\{r\geq 0:D_{f}(r)>t\\},\ \ t\geq 0.$

A driftless subordinator, i.e., $b=0$ with $f(s)=s^{\alpha}$, $0<\alpha<1$ is
known as the stable subordinator. Its first hitting time
$\\{Y_{\alpha}(t)\\}_{t\geq 0}$ is called the inverse stable subordinator.

Recently, Gupta et al. (2020) introduced a time-changed SPoK by time-changing
it with an independent subordinator. They obtained its probability mass
function (pmf), mean, variance and covariance. Here, we introduce and study a
time-changed version of the SPoK by time-changing it with an independent
inverse stable subordinator. It is defined as

$\mathcal{S}_{\alpha}^{k}(t)=\begin{cases}S^{k}\left(Y_{\alpha}(t)\right),\ \
0<\alpha<1,\\\ S^{k}(t),\ \ \alpha=1.\end{cases}$

We call the process $\\{\mathcal{S}_{\alpha}^{k}(t)\\}_{t\geq 0}$ as the
fractional Skellam process of order $k$ (FSPoK). For $k=1$, the FSPoK reduces
to the fractional Skellam process (FSP), for more details on FSP and its
application to high frequency financial data set we refer the reader to Kerss
et al. (2014). For $\alpha=k=1$, the FSPoK reduces to the classical Skellam
process.

We obtain an integral representation for the one-dimensional distributions of
FSPoK. Also, the governing system of fractional differential equations for its
state probabilities is obtained. The probability generating function (pgf),
mean, variance and covariance of the FSPoK are derived and its long-range
dependence (LRD) property is established. Later, we explore its time-changed
versions. We consider the FSPoK time-changed by an independent Lévy
subordinator and its inverse. We compute the mean, variance and covariance for
these time-changed FSPoK. For the first time-changed version we establish the
law of iterated logarithm and the LRD property under suitable restrictions on
the Lévy subordinator. Some particular cases of these time-changed versions
are also considered by taking specific subordinators such as the gamma
subordinator, tempered stable subordinator and inverse Gaussian subordinator.
Also, we obtain the associated system of governing differential equations for
these particular cases. The results obtained in this paper generalize and
complement the results of Kerss et al. (2014) and Gupta et al. (2020).

##  2\. Fractional Skellam process of order $k$

In this section, we introduce and study a stochastic process, namely, the
fractional Skellam process of order $k$ (FSPoK) that is obtained by time-
changing the SPoK by an independent inverse stable subordinator. Let
$\\{Y_{\alpha}(t)\\}_{t\geq 0},\ 0<\alpha<1$, be an inverse stable
subordinator which is independent of SPoK $\\{S^{k}(t)\\}_{t\geq 0}$. The
FSPoK $\\{\mathcal{S}_{\alpha}^{k}(t)\\}_{t\geq 0}$ is defined as

$\mathcal{S}_{\alpha}^{k}(t)=\begin{cases}S^{k}\left(Y_{\alpha}(t)\right),\ \
0<\alpha<1,\\\ S^{k}(t),\ \ \alpha=1.\end{cases}$ (2.1)

For $k=1$, the process defined in (2.1) reduces to the fractional Skellam
process (FSP) which is introduced and studied by Kerss et al. (2014).
Moreover, for $\alpha=k=1$ the FSPoK reduces to the classical Skellam process.

In the following result we derive the system of governing differential
equations that is satisfied by the state probabilities
$p_{\alpha}^{k}(n,t)=\mathrm{Pr}\\{\mathcal{S}^{k}_{\alpha}(t)=n\\}$ of FSPoK.

######  Proposition 2.1.

The state probabilities $p_{\alpha}^{k}(n,t),\ n\in\mathbb{Z}$ of FSPoK solves
the following system of fractional differential equations:

$\partial_{t}^{\alpha}p^{k}_{\alpha}(n,t)=-k(\lambda_{1}+\lambda_{2})p^{k}_{\alpha}(n,t)+\lambda_{1}\sum_{j=1}^{k}p^{k}_{\alpha}(n-j,t)+\lambda_{2}\sum_{j=1}^{k}p^{k}_{\alpha}(n+j,t),$
(2.2)

with the initial conditions $p^{k}_{\alpha}(0,0)=1$ and
$p^{k}_{\alpha}(n,0)=0,\ n\neq 0$. Here, $\partial_{t}^{\alpha}$ denotes the
Caputo fractional derivative defined as (see Kilbas et al. (2006))

$\partial_{t}^{\alpha}f(t):=\left\\{\begin{array}[]{ll}\dfrac{1}{\Gamma{(1-\alpha)}}\displaystyle\int^{t}_{0}(t-s)^{-\alpha}f^{\prime}(s)\,\mathrm{d}s,\
\ 0<\alpha<1,\\\ \\\ f^{\prime}(t),\ \ \alpha=1.\end{array}\right.$

###### Proof.

From (2.1), we have

$p^{k}_{\alpha}(n,t)=\int_{0}^{\infty}p^{k}(n,u)h_{\alpha}(u,t)\,\mathrm{d}u,$
(2.3)

where $h_{\alpha}(u,t)$ is the pdf of $\\{Y_{\alpha}(t)\\}_{t\geq 0}$. For
$\gamma\geq 0$, the Riemann-Liouville (R-L) derivative is defined as (see
Kilbas et al. (2006))

$D_{t}^{\gamma}f(t):=\left\\{\begin{array}[]{ll}\dfrac{1}{\Gamma{(m-\gamma)}}\displaystyle\frac{\mathrm{d}^{m}}{\mathrm{d}t^{m}}\int^{t}_{0}\frac{f(s)}{(t-s)^{\gamma+1-m}}\,\mathrm{d}s,\
\ m-1<\gamma<m,\\\ \\\
\displaystyle\frac{\mathrm{d}^{m}}{\mathrm{d}t^{m}}f(t),\ \
\gamma=m,\end{array}\right.$ (2.4)

where $m$ is a positive integer. Using the following result (see Meerschaert
and Straka (2013))

$D_{t}^{\alpha}h_{\alpha}(u,t)=-\frac{\partial}{\partial u}h_{\alpha}(u,t),$

we get

$\displaystyle D_{t}^{\alpha}p^{k}_{\alpha}(n,t)$
$\displaystyle=-\int_{0}^{\infty}p^{k}(n,u)\frac{\partial}{\partial
u}h_{\alpha}(u,t)\mathrm{d}u$
$\displaystyle=p^{k}(n,0)h_{\alpha}(0+,t)+\int_{0}^{\infty}h_{\alpha}(u,t)\frac{\mathrm{d}}{\mathrm{d}u}p^{k}(n,u)\,\mathrm{d}u$
$\displaystyle=p^{k}(n,0)\frac{t^{-\alpha}}{\Gamma(1-\alpha)}+\int_{0}^{\infty}h_{\alpha}(u,t)\frac{\mathrm{d}}{\mathrm{d}u}p^{k}(n,u)\,\mathrm{d}u,$
(2.5)

where in the last step we used $h_{\alpha}(0+,t)=t^{-\alpha}/\Gamma(1-\alpha)$
(see Meerschaert and Straka (2013)). The following relation holds:

$\partial_{t}^{\alpha}p^{k}_{\alpha}(n,t)=D_{t}^{\alpha}p^{k}_{\alpha}(n,t)-p^{k}_{\alpha}(n,0)\frac{t^{-\alpha}}{\Gamma(1-\alpha)}.$
(2.6)

From (2.3), we have $p^{k}_{\alpha}(n,0)=p^{k}(n,0)$ as
$h_{\alpha}(u,0)=\delta_{0}(u)$. Substituting (1.2) in (2) and then using it
in (2.6), we get

$\displaystyle\partial_{t}^{\alpha}p^{k}_{\alpha}(n,t)$
$\displaystyle=\int_{0}^{\infty}h_{\alpha}(u,t)\frac{\mathrm{d}}{\mathrm{d}u}p^{k}(n,u)\,\mathrm{d}u$
$\displaystyle=\int_{0}^{\infty}\left(-k(\lambda_{1}+\lambda_{2})p^{k}(n,u)+\lambda_{1}\sum_{j=1}^{k}p^{k}(n-j,u)+\lambda_{2}\sum_{j=1}^{k}p^{k}(n+j,u)\right)h_{\alpha}(u,t)\,\mathrm{d}u$
$\displaystyle=-k(\lambda_{1}+\lambda_{2})p^{k}_{\alpha}(n,t)+\lambda_{1}\sum_{j=1}^{k}p^{k}_{\alpha}(n-j,t)+\lambda_{2}\sum_{j=1}^{k}p^{k}_{\alpha}(n+j,t).$

This completes the proof. ∎

On substituting $k=1$ in (2.2), we get the system of governing differential
equations for the state probabilities
$p_{\alpha}(n,t)=\mathrm{Pr}\\{\mathcal{S}_{\alpha}(t)=n\\},\ n\in\mathbb{Z}$
of FSP (see Kerss et al. (2014), Eq. (3.5)) as follows:

$\partial_{t}^{\alpha}p_{\alpha}(n,t)=-(\lambda_{1}+\lambda_{2})p_{\alpha}(n,t)+\lambda_{1}p_{\alpha}(n-1,t)+\lambda_{2}p_{\alpha}(n+1,t),$

with $p_{\alpha}(0,0)=1$ and $p_{\alpha}(n,0)=0,\ n\neq 0$.

To state the next result we need two special functions, namely, the Wright
function $M_{\alpha}(\cdot)$, $0<\alpha<1$ and the modified Bessel function of
first kind $I_{n}(\cdot)$. These are defined as follows (see Mainardi (2010),
Abramowitz and Stegun (1972)):

$M_{\alpha}(z)=\sum_{m=0}^{\infty}\frac{(-z)^{m}}{m!\Gamma(1-m\alpha-\alpha)}$

and

$I_{n}(z)=\sum_{m=0}^{\infty}\frac{(z/2)^{2m+n}}{m!(m+n)!}.$

######  Theorem 2.1.

For $n\in\mathbb{Z}$, the state probability $p_{\alpha}^{k}(n,t)$ of FSPoK is
given by

$p_{\alpha}^{k}(n,t)=\frac{1}{t^{\alpha}}\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{n/2}\int_{0}^{\infty}e^{-ku(\lambda_{1}+\lambda_{2})}I_{|n|}\left(2uk\sqrt{\lambda_{1}\lambda_{2}}\right)M_{\alpha}\left(\frac{u}{t^{\alpha}}\right)\mathrm{d}u.$

###### Proof.

For $n\in\mathbb{Z}$, the state probability $p^{k}(n,t)$ of SPoK is given by
(see Gupta et al. (2020), Eq. (38)):

$p^{k}(n,t)=e^{-kt(\lambda_{1}+\lambda_{2})}\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{n/2}I_{|n|}\left(2tk\sqrt{\lambda_{1}\lambda_{2}}\right).$
(2.7)

The following result holds (see Meerschaert et al. (2015), Section 3):

$h_{\alpha}(u,t)=\frac{1}{t^{\alpha}}M_{\alpha}\left(\frac{u}{t^{\alpha}}\right).$
(2.8)

The result follows on substituting (2.7) and (2.8) in (2.3). ∎

The three-parameter Mittag-Leffler function is defined as (see Kilbas et al.
(2006), p. 45)

$E_{\alpha,\beta}^{\gamma}(x)\coloneqq\frac{1}{\Gamma(\gamma)}\sum_{k=0}^{\infty}\frac{\Gamma(\gamma+k)x^{k}}{k!\Gamma(k\alpha+\beta)},\
\ x\in\mathbb{R},$ (2.9)

where $\alpha>0$, $\beta>0$ and $\gamma>0$. For $\gamma=1$, it reduces to two-
parameter Mittag-Leffler function. For $\gamma=\beta=1$, it further reduces to
the Mittag-Leffler function.

In the next result we obtain the pgf of FSPoK.

######  Proposition 2.2.

The pgf
$G^{k}_{\alpha}(\theta,t)=\mathbb{E}\left(\theta^{\mathcal{S}^{k}_{\alpha}(t)}\right)$,
$0<\theta<1$, of FSPoK is given by

$G^{k}_{\alpha}(\theta,t)=E_{\alpha,1}\left(-\left(k(\lambda_{1}+\lambda_{2})-\lambda_{1}\sum_{j=1}^{k}\theta^{j}-\lambda_{2}\sum_{j=1}^{k}\theta^{-j}\right)t^{\alpha}\right).$

###### Proof.

The pgf of SPoK is given by (see Gupta et al. (2020), Eq. (39))

$G^{k}(\theta,t)=\exp\left(-t\left(k(\lambda_{1}+\lambda_{2})-\lambda_{1}\sum_{j=1}^{k}\theta^{j}-\lambda_{2}\sum_{j=1}^{k}\theta^{-j}\right)\right).$

Using the above result, we get

$\displaystyle G^{k}_{\alpha}(\theta,t)$
$\displaystyle=\int_{0}^{\infty}G^{k}(\theta,u)h_{\alpha}(u,t)\,\mathrm{d}u$
$\displaystyle=\int_{0}^{\infty}\exp\left(-u\left(k(\lambda_{1}+\lambda_{2})-\lambda_{1}\sum_{j=1}^{k}\theta^{j}-\lambda_{2}\sum_{j=1}^{k}\theta^{-j}\right)\right)h_{\alpha}(u,t)\,\mathrm{d}u$
$\displaystyle=E_{\alpha,1}\left(-\left(k(\lambda_{1}+\lambda_{2})-\lambda_{1}\sum_{j=1}^{k}\theta^{j}-\lambda_{2}\sum_{j=1}^{k}\theta^{-j}\right)t^{\alpha}\right).$

This completes the proof. ∎

Using the fact that the Mittag-Leffler function is an eigenfunction of the
Caputo fractional derivative, it follows that

$\partial_{t}^{\alpha}G^{k}_{\alpha}(\theta,t)=-\left(k(\lambda_{1}+\lambda_{2})-\lambda_{1}\sum_{j=1}^{k}\theta^{j}-\lambda_{2}\sum_{j=1}^{k}\theta^{-j}\right)G^{k}_{\alpha}(\theta,t),\
\ G^{k}_{\alpha}(\theta,0)=1.$

######  Proposition 2.3.

The one-dimensional distributions of FSPoK
$\\{\mathcal{S}^{k}_{\alpha}(t)\\}_{t\geq 0}$ are not infinitely divisible.

###### Proof.

From (1.1) and (2.1), we get

$\displaystyle\mathcal{S}^{k}_{\alpha}(t)$
$\displaystyle=N_{1}^{k}\left(Y_{\alpha}(t)\right)-N_{2}^{k}\left(Y_{\alpha}(t)\right)$
$\displaystyle\overset{d}{=}N_{1}^{k}\left(t^{\alpha}Y_{\alpha}(1)\right)-N_{2}^{k}\left(t^{\alpha}Y_{\alpha}(1)\right),$

where $\overset{d}{=}$ means equal in distribution. Here, we have used the
self-similarity property of $\\{Y_{\alpha}(t)\\}_{t\geq 0}$. The following
result holds for PPoK (see Sengar et al. (2020), Eq. (9)):

$\frac{N^{k}(t)}{t}\to\frac{k(k+1)}{2}\lambda,\ \ \mathrm{in\ probability}\
\text{as}\ t\to\infty.$ (2.10)

Thus, $N^{k}(t)/t\to k(k+1)\lambda/2$, in distribution as $t\to\infty$.
Therefore,

$\displaystyle\lim_{t\to\infty}\frac{\mathcal{S}^{k}_{\alpha}(t)}{t^{\alpha}}$
$\displaystyle=\lim_{t\to\infty}\frac{N_{1}^{k}\left(t^{\alpha}Y_{\alpha}(1)\right)-N_{2}^{k}\left(t^{\alpha}Y_{\alpha}(1)\right)}{t^{\alpha}},\
\ \mathrm{in\ distribution}$
$\displaystyle=Y_{\alpha}(1)\left(\lim_{t\to\infty}\frac{N_{1}^{k}\left(t^{\alpha}Y_{\alpha}(1)\right)}{t^{\alpha}Y_{\alpha}(1)}-\lim_{t\to\infty}\frac{N_{2}^{k}\left(t^{\alpha}Y_{\alpha}(1)\right)}{t^{\alpha}Y_{\alpha}(1)}\right)$
$\displaystyle=\frac{k(k+1)}{2}(\lambda_{1}-\lambda_{2})Y_{\alpha}(1),\ \
\mathrm{in\ distribution}.$

It is known that $Y_{\alpha}(1)$ is not infinitely divisible (see Vellaisamy
and Kumar (2018)). Let $\mathcal{S}^{k}_{\alpha}(t)$ be infinitely divisible.
Then, it follows that $\mathcal{S}^{k}_{\alpha}(t)/t^{\alpha}$ is infinitely
divisible (see Steutel and van Harn (2004), Proposition 2.1). As the limit of
a sequence of infinitely divisible random variables is infinitely divisible
(see Steutel and van Harn (2004), Proposition 2.2), it follows that
$Y_{\alpha}(1)$ is infinitely divisible. This leads to a contradiction. ∎

Let us assume that

$r_{1}=k(k+1)(\lambda_{1}-\lambda_{2})/2\ \ \text{and}\ \
r_{2}=k(k+1)(2k+1)(\lambda_{1}+\lambda_{2})/6.$ (2.11)

The mean, variance and covariance function of $\\{S^{k}(t)\\}_{t\geq 0}$ are
given by (see Gupta et al. (2020))

$\mathbb{E}\left(S^{k}(t)\right)=r_{1}t,\ \
\operatorname{Var}\left(S^{k}(t)\right)=r_{2}t\ \ \text{and}\ \
\operatorname{Cov}\left(S^{k}(s),S^{k}(t)\right)=r_{2}s,\ \ 0<s\leq t.$ (2.12)

The mean, variance and covariance of FSPoK are

$\displaystyle\mathbb{E}\left(\mathcal{S}^{k}_{\alpha}(t)\right)$
$\displaystyle=r_{1}\mathbb{E}\left(Y_{\alpha}(t)\right),$ (2.13)
$\displaystyle\operatorname{Var}\left(\mathcal{S}^{k}_{\alpha}(t)\right)$
$\displaystyle=r_{2}\mathbb{E}\left(Y_{\alpha}(t)\right)+r_{1}^{2}\operatorname{Var}\left(Y_{\alpha}(t)\right),$
(2.14)
$\displaystyle\operatorname{Cov}\left(\mathcal{S}^{k}_{\alpha}(s),\mathcal{S}^{k}_{\alpha}(t)\right)$
$\displaystyle=r_{2}\mathbb{E}\left(Y_{\alpha}(s)\right)+r_{1}^{2}\operatorname{Cov}\left(Y_{\alpha}(s),Y_{\alpha}(t)\right),$
(2.15)

which are obtained by using Theorem 2.1 of Leonenko et al. (2014).

######  Remark 2.1.

On substituting $k=1$ in (2.13)-(2.15), we get the mean, variance and
covariance of FSP (see Kerss et al. (2014), Remark 3.2).

Next, we show that the FSPoK possesses the LRD property. The following
definition will be used (see D’Ovidio and Nane (2014), Maheshwari and
Vellaisamy (2016)):

######  Definition 2.1.

Let $s>0$ be fixed and $\\{X(t)\\}_{t\geq 0}$ be a stochastic process such
that

$\lim_{t\to\infty}\frac{\operatorname{Cov}(X(s),X(t))}{\sqrt{\operatorname{Var}(X(s))}\sqrt{\operatorname{Var}(X(t))}t^{-\gamma}}=c(s).$

If $\gamma\in(0,1)$ the process $\\{X(t)\\}_{t\geq 0}$ has the LRD property,
and if $\gamma\in(1,2)$ it has the SRD property.

######  Theorem 2.2.

The FSPoK has the LRD property.

###### Proof.

For a fixed $s>0$, we have

$\lim_{t\to\infty}\frac{\operatorname{Cov}\left(\mathcal{S}^{k}_{\alpha}(s),\mathcal{S}^{k}_{\alpha}(t)\right)}{\sqrt{\operatorname{Var}\left(\mathcal{S}^{k}_{\alpha}(t)\right)}\sqrt{\operatorname{Var}\left(\mathcal{S}^{k}_{\alpha}(t)\right)}t^{-\alpha}}=c(s),$
(2.16)

where

$c(s)=\left(\frac{1}{\Gamma(2\alpha)}-\frac{1}{\alpha(\Gamma(\alpha))^{2}}\right)^{-1}\left(\frac{\alpha
r_{2}}{\Gamma(1+\alpha)r_{1}^{2}}+\frac{\alpha
s^{\alpha}}{\Gamma(1+2\alpha)}\right).$

The result in (2.16) follows by using a result on p. $10$ of Leonenko et al.
(2014), and the mean and variance of SPoK which are given in (2.12). Thus, the
FSPoK exhibits the LRD property as $0<\alpha<1$. ∎

######  Remark 2.2.

For a fixed $h>0$, the increment of FSPoK is defined as

$Z^{k}_{h}(t)=\mathcal{S}^{k}_{\alpha}(t+h)-\mathcal{S}^{k}_{\alpha}(t),\ \
t\geq 0.$

It can be shown that the increment process $\\{Z^{k}_{h}(t)\\}_{t\geq 0}$
exhibits the SRD property. The proof follows similar lines to that of Theorem
1 of Maheshwari and Vellaisamy (2016), and thus it is omitted.

##  3\. The FSPoK time-changed by a Lévy subordinator

In this section, we consider a time-changed version of the FSPoK. We call it
the time-changed fractional Skellam process of order $k$ (TCFSPoK) and denote
it by $\\{\mathcal{Z}^{f}_{\alpha}(t)\\}_{t\geq 0}$, $0<\alpha\leq 1$. It is
defined as the FSPoK time-changed by an independent Lévy subordinator
$\\{D_{f}(t)\\}_{t\geq 0}$ with $\mathbb{E}\left(D^{r}_{f}(t)\right)<\infty$
for all $r>0$. Thus,

$\mathcal{Z}^{f}_{\alpha}(t)\coloneqq\mathcal{S}^{k}_{\alpha}(D_{f}(t)),\ \
t\geq 0,$ (3.1)

where $\\{\mathcal{S}^{k}_{\alpha}(t)\\}_{t\geq 0}$ is independent of
$\\{D_{f}(t)\\}_{t\geq 0}$.

For $\alpha=1$, the TCFSPoK reduces to a time-changed version of the SPoK (see
Gupta et al. (2020)), that is,

$\mathcal{Z}^{f}(t)\coloneqq\mathcal{S}^{k}(D_{f}(t)),\ \ t\geq 0.$

Its pmf $p^{f}(n,t)=\mathrm{Pr}\\{\mathcal{Z}^{f}(t)=n\\}$ is given by (see
Gupta et al. (2020), Eq. (47))

$p^{f}(n,t)=\sum_{x=\max(0,-n)}^{\infty}\frac{(k\lambda_{1})^{n+x}(k\lambda_{2})^{x}}{(n+x)!x!}\mathbb{E}\left(e^{-k(\lambda_{1}+\lambda_{2})D_{f}(t)}D_{f}^{2n+x}(t)\right),\
\ n\in\mathbb{Z}.$

The mean and covariance of inverse stable subordinator are given by (see
Leonenko et al. (2014), Eq. (8) and Eq. (10))

$\mathbb{E}\left(Y_{\alpha}(t)\right)=\frac{t^{\alpha}}{\Gamma(\alpha+1)}$
(3.2)

and

$\operatorname{Cov}\left(Y_{\alpha}(s),Y_{\alpha}(t)\right)=\frac{1}{\Gamma^{2}(\alpha+1)}\left(\alpha
s^{2\alpha}B(\alpha,\alpha+1)+F(\alpha;s,t)\right),\ \ 0<s\leq t,$ (3.3)

where $F(\alpha;s,t)=\alpha t^{2\alpha}B(\alpha,\alpha+1;s/t)-(ts)^{\alpha}$.
Here, $B(\alpha,\alpha+1)$ and $B(\alpha,\alpha+1;s/t)$ denote the beta
function and the incomplete beta function, respectively.

Let $l_{1}=r_{1}/\Gamma(\alpha+1)$, $l_{2}=r_{2}/\Gamma(\alpha+1)$ and
$d=l_{1}^{2}\alpha B(\alpha,\alpha+1)$ where $r_{1}$ and $r_{2}$ are given by
(2.11).

The mean of TCFSPoK is obtained as follows:

$\mathbb{E}\left(\mathcal{Z}^{f}_{\alpha}(t)\right)=\mathbb{E}\left(\mathbb{E}\left(\mathcal{S}^{k}_{\alpha}(D_{f}(t))|D_{f}(t)\right)\right)=l_{1}\mathbb{E}\left(D_{f}^{\alpha}(t)\right).$
(3.4)

On substituting (3.2) and (3.3) in (2.15), we get

$\mathbb{E}\left(\mathcal{S}^{k}_{\alpha}(s)\mathcal{S}^{k}_{\alpha}(t)\right)=l_{2}s^{\alpha}+ds^{2\alpha}+\alpha
l_{1}^{2}\left(t^{2\alpha}B(\alpha,\alpha+1;s/t)\right).$

Thus,

$\displaystyle\mathbb{E}\left(\mathcal{Z}^{f}_{\alpha}(s)\mathcal{Z}^{f}_{\alpha}(t)\right)$
$\displaystyle=\mathbb{E}\left(\mathbb{E}\left(\mathcal{S}^{k}_{\alpha}(D_{f}(s))\mathcal{S}^{k}_{\alpha}(D_{f}(t))|D_{f}(s),D_{f}(t)\right)\right)$
$\displaystyle=l_{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)+d\mathbb{E}\left(D_{f}^{2\alpha}(s)\right)+\alpha
l_{1}^{2}\mathbb{E}\left(D_{f}^{2\alpha}(t)B\left(\alpha,\alpha+1;D_{f}(s)/D_{f}(t)\right)\right).$

Hence, the covariance of TCFSPoK is given by

$\displaystyle\operatorname{Cov}\left(\mathcal{Z}^{f}_{\alpha}(s),\mathcal{Z}^{f}_{\alpha}(t)\right)$
$\displaystyle=l_{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)+d\mathbb{E}\left(D_{f}^{2\alpha}(s)\right)-l_{1}^{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)\mathbb{E}\left(D_{f}^{\alpha}(t)\right)$
$\displaystyle\ \ +\alpha
l_{1}^{2}\mathbb{E}\left(D_{f}^{2\alpha}(t)B\left(\alpha,\alpha+1;D_{f}(s)/D_{f}(t)\right)\right).$
(3.5)

Also, its variance is obtained by substituting $s=t$ in (3) which is given by

$\operatorname{Var}\left(\mathcal{Z}^{f}_{\alpha}(t)\right)=\mathbb{E}\left(D_{f}^{\alpha}(t)\right)\left(l_{2}-l_{1}^{2}\mathbb{E}\left(D_{f}^{\alpha}(t)\right)\right)+2d\mathbb{E}\left(D_{f}^{2\alpha}(t)\right).$
(3.6)

######  Remark 3.1.

On substituting $k=1$ and taking $\lambda_{2}\to 0$ in (3.4)-(3.6), we get the
mean, covariance and variance of a time-changed version of the fractional
Poisson process (see Maheshwari and Vellaisamy (2019), Theorem 3.2).

Next we show that the TCFSPoK exhibits the LRD property under certain
restrictions on the Lévy subordinator. For positive functions $f$ and $g$, the
notation $f\sim g$ stands for $f(t)/g(t)\to 1$ as $t\to\infty$.

######  Theorem 3.1.

Let $\\{D_{f}(t)\\}_{t\geq 0}$ be a Lévy subordinator with $f(\cdot)$ being
its associated Bernštein function defined in (1.3) such that
$\mathbb{E}\left(D^{r}_{f}(t)\right)<\infty$ for all $r>0$. The TCFSPoK
$\\{\mathcal{Z}^{f}_{\alpha}(t)\\}_{t\geq 0}$, $0<\alpha<1$ exhibits the LRD
property if

$\mathbb{E}\left(D_{f}^{i\alpha}(t)\right)\sim k_{i}t^{i\rho},\ \ i=1,2,$
(3.7)

for some $0<\rho<1$, and positive constants $k_{1}$ and $k_{2}$ such that
$k_{2}\geq k_{1}^{2}$.

###### Proof.

Let $0<s<t$. In (3), we use the following asymptotic result for large $t$ (see
Maheshwari and Vellaisamy (2019), Theorem 3.3)

$\alpha\mathbb{E}\left(D_{f}^{2\alpha}(t)B\left(\alpha,\alpha+1;D_{f}(s)/D_{f}(t)\right)\right)\sim\mathbb{E}\left(D_{f}^{\alpha}(s)\right){\mathbb{E}}\left(D_{f}^{\alpha}(t-s)\right),$

to obtain

$\displaystyle\operatorname{Cov}\left(\mathcal{Z}^{f}_{\alpha}(s),\mathcal{Z}^{f}_{\alpha}(t)\right)$
$\displaystyle\sim
l_{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)+d\mathbb{E}\left(D_{f}^{2\alpha}(s)\right)-l_{1}^{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)\mathbb{E}\left(D_{f}^{\alpha}(t)\right)$
$\displaystyle\ \ \
+l_{1}^{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right){\mathbb{E}}\left(D_{f}^{\alpha}(t-s)\right)$
$\displaystyle\sim
l_{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)+d\mathbb{E}\left(D_{f}^{2\alpha}(s)\right)-l_{1}^{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)k_{1}(t^{\rho}-(t-s)^{\rho})$
$\displaystyle\sim
l_{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)+d\mathbb{E}\left(D_{f}^{2\alpha}(s)\right)-l_{1}^{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)k_{1}s\rho
t^{\rho-1},$

where we have used (3.7) in the penultimate step. Again by using (3.7) in
(3.6), we get

$\displaystyle\operatorname{Var}\left(\mathcal{Z}^{f}_{\alpha}(t)\right)$
$\displaystyle\sim
l_{2}k_{1}t^{\rho}-l_{1}^{2}k^{2}_{1}t^{2\rho}+2dk_{2}t^{2\rho}$
$\displaystyle\sim\left(2dk_{2}-k^{2}_{1}l_{1}^{2}\right)t^{2\rho}$
$\displaystyle=\frac{r_{1}^{2}}{\alpha}\left(\frac{k_{2}}{\Gamma(2\alpha)}-\frac{k_{1}^{2}}{\alpha\Gamma^{2}(\alpha)}\right)t^{2\rho}.$

Thus, for large $t$, we have

$\displaystyle\operatorname{Corr}\left(\mathcal{Z}^{f}_{\alpha}(s),\mathcal{Z}^{f}_{\alpha}(t)\right)$
$\displaystyle\sim\frac{l_{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)+d\mathbb{E}\left(D_{f}^{2\alpha}(s)\right)-l_{1}^{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)k_{1}s\rho
t^{\rho-1}}{\sqrt{\operatorname{Var}\left(\mathcal{Z}^{f}_{\alpha}(s)\right)}\sqrt{\left(2dk_{2}-k^{2}_{1}l_{1}^{2}\right)t^{2\rho}}}$
$\displaystyle\sim c_{1}(s)t^{-\rho},$

where

$c_{1}(s)=\frac{l_{2}\mathbb{E}\left(D_{f}^{\alpha}(s)\right)+d\mathbb{E}\left(D_{f}^{2\alpha}(s)\right)}{\sqrt{\operatorname{Var}\left(\mathcal{Z}^{f}_{\alpha}(s)\right)\left(2dk_{2}-k^{2}_{1}l_{1}^{2}\right)}}.$

Hence, the TCFSPoK exhibits the LRD property as $0<\rho<1$. ∎

######  Remark 3.2.

In a similar way it can be shown that $\\{\mathcal{Z}^{f}(t)\\}_{t\geq 0}$
exhibits the LRD property.

The following result will be used to prove the law of iterated logarithm (LIL)
for TCFSPoK (see Bertoin (1996), Theorem 14, p. 92).

######  Lemma 3.1.

Let $\\{D_{f}(t)\\}_{t\geq 0}$ be a Lévy subordinator whose associated
Bernštein function $f$ is regularly varying at $0+$ with index $0<\gamma<1$,
i.e., $\lim_{x\rightarrow 0+}f(\lambda x)/f(x)=\lambda^{\gamma}$, $\lambda>0$.
Also, let

$g(t)=\frac{\log\log t}{\phi(t^{-1}\log\log t)},\ \ t>e,$ (3.8)

where $\phi$ is the inverse of $f$. Then,

$\liminf_{t\to\infty}\frac{D_{f}(t)}{g(t)}=\gamma(1-\gamma)^{(1-\gamma)/\gamma},\
\ \text{a.s.}$ (3.9)

######  Theorem 3.2.

Let the Bernštein function $f(\cdot)$ associated with Lévy subordinator
$\\{D_{f}(t)\\}_{t\geq 0}$ be regularly varying at 0+ with index $0<\gamma<1$.
Then,

$\liminf_{t\rightarrow\infty}\frac{\mathcal{Z}^{f}_{\alpha}(t)}{(g(t))^{\alpha}}\stackrel{{\scriptstyle
d}}{{=}}\frac{k(k+1)}{2}(\lambda_{1}-\lambda_{2})Y_{\alpha}(1)\gamma^{\alpha}\left(1-\gamma\right)^{\alpha(1-\gamma)/\gamma},$
(3.10)

where $g(t)$ is given in (3.8).

###### Proof.

Note that $D_{f}(t)\to\infty$ as $t\to\infty$, a.s. (see Bertoin (1996), p.
73). From (1.1), (2.1) and (3.1), we get

$\displaystyle\mathcal{Z}^{f}_{\alpha}(t)$
$\displaystyle=N_{1}^{k}(Y_{\alpha}(D_{f}(t)))-N_{2}^{k}(Y_{\alpha}(D_{f}(t)))$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}N_{1}^{k}(D^{\alpha}_{f}(t)Y_{\alpha}(1))-N_{2}^{k}(D^{\alpha}_{f}(t)Y_{\alpha}(1)),$

where we have used the self-similarity property of $\\{Y_{\alpha}(t)\\}_{t\geq
0}$. Hence,

$\displaystyle\liminf_{t\rightarrow\infty}\frac{\mathcal{Z}^{f}_{\alpha}(t)}{(g(t))^{\alpha}}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\liminf_{t\rightarrow\infty}\frac{N_{1}^{k}(D^{\alpha}_{f}(t)Y_{\alpha}(1))-N_{2}^{k}(D^{\alpha}_{f}(t)Y_{\alpha}(1))}{(g(t))^{\alpha}}$
$\displaystyle=\liminf_{t\rightarrow\infty}\left(\frac{N_{1}^{k}(D^{\alpha}_{f}(t)Y_{\alpha}(1))-N_{2}^{k}(D^{\alpha}_{f}(t)Y_{\alpha}(1))}{D^{\alpha}_{f}(t)Y_{\alpha}(1)}\right)\frac{D^{\alpha}_{f}(t)Y_{\alpha}(1)}{(g(t))^{\alpha}}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{k(k+1)}{2}(\lambda_{1}-\lambda_{2})Y_{\alpha}(1)\left(\liminf_{t\rightarrow\infty}\frac{D_{f}(t)}{g(t)}\right)^{\alpha},\
\ (\text{using}\ (\ref{limit}))$ $\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{k(k+1)}{2}(\lambda_{1}-\lambda_{2})Y_{\alpha}(1)\gamma^{\alpha}\left(1-\gamma\right)^{\alpha(1-\gamma)/\gamma},$

where the last step follows from (3.9). ∎

######  Remark 3.3.

On substituting $k=1$ and taking $\lambda_{2}\to 0$ in (3.10), we get the LIL
for a time-changed version of the fractional Poisson process (see Maheshwari
and Vellaisamy (2019), Theorem 3.5).

###  3.1. Some special cases of the TCFSPoK

In this subsection, we time-change the FSPoK and the SPoK by three specific
Lévy subordinators, namely, the gamma subordinator, the tempered stable
subordinator and the inverse Gaussian subordinator. We obtain the governing
systems of differential equations for their one-dimensional distributions.

####  3.1.1. FSPoK time-changed by gamma subordinator

Let $\\{Z(t)\\}_{t\geq 0}$ be a gamma subordinator with the following
probability density function (pdf):

$g(x,t)=\frac{a^{bt}}{\Gamma(bt)}x^{bt-1}e^{-ax},\ \ x>0,$

where $a>0$ and $b>0$. The Bernštein function $f_{1}(s)$ associated with
$\\{Z(t)\\}_{t\geq 0}$ is given by $f_{1}(s)=b\log(1+s/a),\ s>0$ (see
Applebaum (2009), p. 55).

The FSPoK time-changed by an independent gamma subordinator is defined as

$\mathcal{Z}^{f_{1}}_{\alpha}(t)\coloneqq\mathcal{S}^{k}_{\alpha}(Z(t)),\ \
t\geq 0.$ (3.11)

The following result will be used to obtain the governing system of
differential equations for its pmf
$p_{\alpha}^{f_{1}}(n,t)=\mathrm{Pr}\\{\mathcal{Z}^{f_{1}}_{\alpha}(t)=n\\}$,
$n\in\mathbb{Z}$.

######  Lemma 3.2 (Vellaisamy and Maheshwari (2018)).

For any $\gamma\geq 1$, the pdf $g(x,t)$ of gamma subordinator solves

$\displaystyle D_{t}^{\gamma}g(x,t)$
$\displaystyle=bD_{t}^{\gamma-1}\left(\log(ax)-\psi(bt)\right)g(x,t),\ \ \
x>0,$ $\displaystyle g(x,0)$ $\displaystyle=0.$

Here, $\psi(x)\coloneqq\Gamma^{\prime}(x)/\Gamma(x)$ is the digamma function
and $D_{t}^{\gamma}$ is the R-L fractional derivative defined in (2.4).

######  Theorem 3.3.

Let $\gamma\geq 1$ and $\psi(x)$ be the digamma function. Then, the pmf
$p_{\alpha}^{f_{1}}(n,t)$ solves the following equation:

$D_{t}^{\gamma}p_{\alpha}^{f_{1}}(n,t)=bD_{t}^{\gamma-1}\left(\log(a)-\psi(bt)\right)p_{\alpha}^{f_{1}}(n,t)+b\int_{0}^{\infty}p^{k}_{\alpha}(n,x)\log(x)D_{t}^{\gamma-1}g(x,t)\,\mathrm{d}x.$

###### Proof.

From (3.11), we have

$p_{\alpha}^{f_{1}}(n,t)=\int_{0}^{\infty}p^{k}_{\alpha}(n,x)g(x,t)\,\mathrm{d}x.$
(3.12)

Taking the R-L fractional derivative in (3.12) and using Lemma 3.2, we get

$\displaystyle D_{t}^{\gamma}p_{\alpha}^{f_{1}}(n,t)$
$\displaystyle=\int_{0}^{\infty}p^{k}_{\alpha}(n,x)D_{t}^{\gamma}g(x,t)\,\mathrm{d}x$
$\displaystyle=b\int_{0}^{\infty}p^{k}_{\alpha}(n,x)D_{t}^{\gamma-1}\left(\log(ax)-\psi(bt)\right)g(x,t)\,\mathrm{d}x$
$\displaystyle=bD_{t}^{\gamma-1}\log(a)\int_{0}^{\infty}p^{k}_{\alpha}(n,x)g(x,t)\,\mathrm{d}x+b\int_{0}^{\infty}p^{k}_{\alpha}(n,x)\log(x)D_{t}^{\gamma-1}g(x,t)\,\mathrm{d}x$
$\displaystyle\ \
-bD_{t}^{\gamma-1}\psi(bt)\int_{0}^{\infty}p^{k}_{\alpha}(n,x)g(x,t)\,\mathrm{d}x$
$\displaystyle=bD_{t}^{\gamma-1}\left(\log(a)-\psi(bt)\right)p_{\alpha}^{f_{1}}(n,t)+b\int_{0}^{\infty}p^{k}_{\alpha}(n,x)\log(x)D_{t}^{\gamma-1}g(x,t)\,\mathrm{d}x.$

This completes the proof. ∎

Next we discuss two particular cases of the time-changed SPoK (see Gupta et
al. (2020)).

####  3.1.2. SPoK time-changed by tempered stable subordinator

Let $\\{\mathscr{D}_{\eta,\nu}(t)\\}_{t\geq 0}$ denote the tempered stable
subordinator (TSS) with stability index $0<\nu<1$ and the tempering parameter
$\eta>0$. The Bernštein function $f_{2}(s)$ associated with TSS is given by

$f_{2}(s)=(\eta+s)^{\nu}-\eta^{\nu},\ s>0.$ (3.13)

The SPoK time-changed by an independent TSS is defined as

$\mathcal{Z}^{f_{2}}(t)\coloneqq\mathcal{S}^{k}(\mathscr{D}_{\eta,\nu}(t)),\ \
t\geq 0.$ (3.14)

######  Proposition 3.1.

The pmf $p^{f_{2}}(n,t)=\mathrm{Pr}\\{\mathcal{Z}^{f_{2}}(t)=n\\}$,
$n\in\mathbb{Z}$, is the solution of the following differential equation:

$\left(\eta^{\nu}-\frac{\mathrm{d}}{\mathrm{d}t}\right)^{1/\nu}p^{f_{2}}(n,t)=\eta
p^{f_{2}}(n,t)+k(\lambda_{1}+\lambda_{2})p^{f_{2}}(n,t)-\lambda_{1}\sum_{j=1}^{k}p^{f_{2}}(n-j,t)-\lambda_{2}\sum_{j=1}^{k}p^{f_{2}}(n+j,t).$

###### Proof.

From (3.14), we have

$p^{f_{2}}(n,t)=\int_{0}^{\infty}p^{k}(n,x)h_{\eta,\nu}(x,t)\,\mathrm{d}x.$

Here, $h_{\eta,\nu}(x,t)$ denotes the pdf of TSS. It is known that (See Beghin
(2015), Eq. (15))

$\frac{\partial}{\partial x}h_{\eta,\nu}(x,t)=-\eta
h_{\eta,\nu}(x,t)+\left(\eta^{\nu}-\frac{\partial}{\partial
t}\right)^{1/\nu}h_{\eta,\nu}(x,t),$

with initial conditions $h_{\eta,\nu}(x,0)=\delta_{0}(x)$ and
$h_{\eta,\nu}(0,t)=0$. Here, $\delta_{0}(x)$ is the Dirac delta function.
Using the following results: $\lim_{x\to
0}h_{\eta,\nu}(x,t)=\lim_{x\to\infty}h_{\eta,\nu}(x,t)=0$, we get

$\displaystyle\left(\eta^{\nu}-\frac{\mathrm{d}}{\mathrm{d}t}\right)^{1/\nu}p^{f_{2}}(n,t)$
$\displaystyle=\int_{0}^{\infty}p^{k}(n,x)\left(\frac{\partial}{\partial
x}h_{\eta,\nu}(x,t)+\eta h_{\eta,\nu}(x,t)\right)\,\mathrm{d}x$
$\displaystyle=\eta
p^{f_{2}}(n,t)-\int_{0}^{\infty}h_{\eta,\nu}(x,t)\frac{\mathrm{d}}{\mathrm{d}x}p^{k}(n,x)\,\mathrm{d}x$
$\displaystyle=\eta
p^{f_{2}}(n,t)-\int_{0}^{\infty}\bigg{(}-k(\lambda_{1}+\lambda_{2})p^{k}(n,x)+\lambda_{1}\sum_{j=1}^{k}p^{k}(n-j,x)$
$\displaystyle\ \
+\lambda_{2}\sum_{j=1}^{k}p^{k}(n+j,x)\bigg{)}h_{\eta,\nu}(x,t)\,\mathrm{d}x,\
\ {(\text{using (\ref{skellam gov})})}$ $\displaystyle=\eta
p^{f_{2}}(n,t)+k(\lambda_{1}+\lambda_{2})p^{f_{2}}(n,t)$ $\displaystyle\ \
-\lambda_{1}\sum_{j=1}^{k}p^{f_{2}}(n-j,t)-\lambda_{2}\sum_{j=1}^{k}p^{f_{2}}(n+j,t).$

This completes the proof. ∎

If $\nu^{-1}=m\geq 2$ is an integer then the pmf $p^{f_{2}}(n,t)$ solves

$\sum_{k=1}^{m}(-1)^{k}\binom{m}{k}\eta^{(1-k/m)}\frac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}p^{f_{2}}(n,t)=k(\lambda_{1}+\lambda_{2})p^{f_{2}}(n,t)-\lambda_{1}\sum_{j=1}^{k}p^{f_{2}}(n-j,t)-\lambda_{2}\sum_{j=1}^{k}p^{f_{2}}(n+j,t).$

####  3.1.3. SPoK time-changed by inverse Gaussian subordinator

Let $\\{Y(t)\\}_{t\geq 0}$ be an inverse Gaussian (IG) subordinator with the
following pdf (see Applebaum (2009), Eq. (1.27))

$q(x,t)=(2\pi)^{-1/2}\delta tx^{-3/2}e^{\delta\gamma
t-\frac{1}{2}(\delta^{2}t^{2}x^{-1}+\gamma^{2}x)},\ \ x>0,\ \delta>0,\
\gamma>0$

and the associated Bernštein function

$f_{3}(s)=\delta\left(\sqrt{2s+\gamma^{2}}-\gamma\right),\ \ s>0.$ (3.15)

The SPoK time-changed by an independent IG subordinator is defined as

$\mathcal{Z}^{f_{3}}(t)\coloneqq\mathcal{S}^{k}(Y(t)),\ \ t\geq 0.$ (3.16)

######  Proposition 3.2.

The pmf $p^{f_{3}}(n,t)=\mathrm{Pr}\\{\mathcal{Z}^{f_{3}}(t)=n\\}$,
$n\in\mathbb{Z}$, solves the following differential equation:

$\displaystyle\left(\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}-2\delta\gamma\frac{\mathrm{d}}{\mathrm{d}t}\right)p^{f_{3}}(n,t)$
$\displaystyle=2\delta^{2}\Bigg{(}k(\lambda_{1}+\lambda_{2})p^{f_{3}}(n,t)$
$\displaystyle\hskip
85.35826pt-\lambda_{1}\sum_{j=1}^{k}p^{f_{3}}(n-j,t)-\lambda_{2}\sum_{j=1}^{k}p^{f_{3}}(n+j,t)\Bigg{)}.$

###### Proof.

From (3.16), we have

$p^{f_{3}}(n,t)=\int_{0}^{\infty}p^{k}(n,x)q(x,t)\,\mathrm{d}x.$ (3.17)

On taking derivatives, we get

$\frac{\mathrm{d}}{\mathrm{d}t}p^{f_{3}}(n,t)=\int_{0}^{\infty}p^{k}(n,x)\frac{\partial}{\partial
t}q(x,t)\,\mathrm{d}x$

and

$\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}p^{f_{3}}(n,t)=\int_{0}^{\infty}p^{k}(n,x)\frac{\partial^{2}}{\partial
t^{2}}q(x,t)\,\mathrm{d}x.$

The following results hold for $q(x,t)$ (see Vellaisamy and Kumar (2018), Eq.
(3.3))

$\frac{\partial^{2}}{\partial
t^{2}}q(x,t)-2\delta\gamma\frac{\partial}{\partial
t}q(x,t)=2\delta^{2}\frac{\partial}{\partial x}q(x,t)$ (3.18)

and $\lim_{x\to\infty}q(x,t)=\lim_{x\to 0}q(x,t)=0$. From (3.17), we get

$\displaystyle\left(\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}-2\delta\gamma\frac{\mathrm{d}}{\mathrm{d}t}\right)p^{f_{3}}(n,t)$
$\displaystyle=\int_{0}^{\infty}p^{k}(n,x)\left(\frac{\partial^{2}}{\partial
t^{2}}-2\delta\gamma\frac{\partial}{\partial t}\right)q(x,t)\,\mathrm{d}x$
$\displaystyle=2\delta^{2}\int_{0}^{\infty}p^{k}(n,x)\frac{\partial}{\partial
x}q(x,t)\,\mathrm{d}x,\ \ (\text{using (\ref{qxt212})})$
$\displaystyle=-2\delta^{2}\int_{0}^{\infty}q(x,t)\frac{\mathrm{d}}{\mathrm{d}x}p^{k}(n,x)\,\mathrm{d}x$
$\displaystyle=-2\delta^{2}\int_{0}^{\infty}\bigg{(}-k(\lambda_{1}+\lambda_{2})p^{k}(n,x)+\lambda_{1}\sum_{j=1}^{k}p^{k}(n-j,x)$
$\displaystyle\ \
+\lambda_{2}\sum_{j=1}^{k}p^{k}(n+j,x)\bigg{)}q(x,t)\,\mathrm{d}x,\ \
(\text{using (\ref{skellam gov})})$
$\displaystyle=2\delta^{2}\Bigg{(}k(\lambda_{1}+\lambda_{2})p^{f_{3}}(n,t)-\lambda_{1}\sum_{j=1}^{k}p^{f_{3}}(n-j,t)$
$\displaystyle\hskip
113.81102pt-\lambda_{2}\sum_{j=1}^{k}p^{f_{3}}(n+j,t)\Bigg{)}.$

This completes the proof. ∎

##  4\. The FSPoK time-changed by inverse subordinator

The first hitting time of subordinator $\\{D_{f}(t)\\}_{t\geq 0}$ is called
the inverse subordinator. It is defined as

$H_{f}(t)\coloneqq\inf\\{r\geq 0:D_{f}(r)>t\\},\ \ t\geq 0.$

It is known that $\mathbb{E}\left(H^{r}_{f}(t)\right)<\infty$ for all $r>0$
(see Aletti et al. (2018), Section 2.1).

Let

$\bar{\mathcal{Z}}^{f}_{\alpha}(t)\coloneqq\mathcal{S}^{k}_{\alpha}(H_{f}(t)),\
\ t\geq 0,$ (4.1)

be the FSPoK time-changed by an independent inverse subordinator.

For $\alpha=1$, the process defined in (4.1) reduces to a time-changed version
of the SPoK, that is,

$\bar{\mathcal{Z}}^{f}(t)\coloneqq\mathcal{S}^{k}(H_{f}(t)),\ \ t\geq 0.$

The pmf $\bar{p}^{f}(n,t)=\mathrm{Pr}\\{\bar{\mathcal{Z}}^{f}(t)=n\\}$ of
$\\{\bar{\mathcal{Z}}^{f}(t)\\}_{t\geq 0}$ is given by

$\bar{p}^{f}(n,t)=\sum_{x=\max(0,-n)}^{\infty}\frac{(k\lambda_{1})^{n+x}(k\lambda_{2})^{x}}{(n+x)!x!}\mathbb{E}\left(e^{-k(\lambda_{1}+\lambda_{2})H_{f}(t)}H_{f}^{2n+x}(t)\right),\
\ n\in\mathbb{Z}.$

The proof follows along the similar lines to that of Theorem 2 of Gupta et al.
(2020).

The mean, variance and covariance of
$\\{\bar{\mathcal{Z}}^{f}_{\alpha}(t)\\}_{t\geq 0}$ are given as follows:

Let $l_{1}=r_{1}/\Gamma(\alpha+1)$, $l_{2}=r_{2}/\Gamma(\alpha+1)$ and
$d=l_{1}^{2}\alpha B(\alpha,\alpha+1)$ where $B(\alpha,\alpha+1)$ is the beta
function, and $r_{1}$ and $r_{2}$ are given by (2.11). Then,

$(i)\ \
\mathbb{E}\left(\bar{\mathcal{Z}}^{f}_{\alpha}(t)\right)=l_{1}\mathbb{E}\left(H_{f}^{\alpha}(t)\right)$,
$(ii)\
\operatorname{Var}\left(\bar{\mathcal{Z}}^{f}_{\alpha}(t)\right)=\mathbb{E}\left(H_{f}^{\alpha}(t)\right)\left(l_{2}-l_{1}^{2}\mathbb{E}\left(H_{f}^{\alpha}(t)\right)\right)+2d\mathbb{E}\left(H_{f}^{2\alpha}(t)\right)$,
$(iii)\
\operatorname{Cov}\left(\bar{\mathcal{Z}}^{f}_{\alpha}(s),\bar{\mathcal{Z}}^{f}_{\alpha}(t)\right)=l_{2}\mathbb{E}\left(H_{f}^{\alpha}(s)\right)+d\mathbb{E}\left(H_{f}^{2\alpha}(s)\right)-l_{1}^{2}\mathbb{E}\left(H_{f}^{\alpha}(s)\right)\mathbb{E}\left(H_{f}^{\alpha}(t)\right)$
$\hskip 130.88284pt+\alpha
l_{1}^{2}\mathbb{E}\left(H_{f}^{2\alpha}(t)B\left(\alpha,\alpha+1;H_{f}(s)/H_{f}(t)\right)\right)$,

where $0<s\leq t$ and $B(\alpha,\alpha+1;H_{f}(s)/H_{f}(t))$ is the incomplete
beta function.

The proof of $(i)$-$(iii)$ follows similar lines to the corresponding results
of $\\{\mathcal{Z}^{f}_{\alpha}(t)\\}_{t\geq 0}$ given in the previous
section.

Next we discuss two particular cases of the time-changed process
$\\{\bar{\mathcal{Z}}^{f}(t)\\}_{t\geq 0}$.

###  4.1. SPoK time-changed by the inverse TSS

Let $\\{\mathscr{L}_{\eta,\nu}(t)\\}_{t\geq 0}$ denote the inverse TSS which
is defined as the first hitting time of TSS
$\\{\mathscr{D}_{\eta,\nu}(t)\\}_{t\geq 0}$, $0<\nu<1$, $\eta>0$. That is,

$\mathscr{L}_{\eta,\nu}(t)=\inf\\{s\geq 0:\mathscr{D}_{\eta,\nu}(s)>t\\},\ \
t\geq 0.$

The SPoK time-changed by an independent inverse TSS is defined as

$\bar{\mathcal{Z}}^{f_{2}}(t)\coloneqq\mathcal{S}^{k}(\mathscr{L}_{\eta,\nu}(t)),\
\ t\geq 0,$ (4.2)

where the associated Bernštein function $f_{2}$ is given in (3.13).

######  Proposition 4.1.

For $n\in\mathbb{Z}$, the pmf
$\bar{p}^{f_{2}}(n,t)=\mathrm{Pr}\\{\bar{\mathcal{Z}}^{f_{2}}(t)=n\\}$ solves
the following differential equation:

$\displaystyle\left(\eta+\frac{\mathrm{d}}{\mathrm{d}t}\right)^{\nu}\bar{p}^{f_{2}}(n,t)$
$\displaystyle=p^{k}(n,x)l_{\eta,\nu}(x,t)\big{|}_{x=0}-k(\lambda_{1}+\lambda_{2})\bar{p}^{f_{2}}(n,t)+\lambda_{1}\sum_{j=1}^{k}\bar{p}^{f_{2}}(n-j,t)$
$\displaystyle\ \
+\lambda_{2}\sum_{j=1}^{k}\bar{p}^{f_{2}}(n+j,t)+\eta^{\nu}\bar{p}^{f_{2}}(n,t)-t^{-\nu}E^{1-\nu}_{1,1-\nu}(-\eta
t)p^{k}(n,0),$

where $E^{1-\nu}_{1,1-\nu}(\cdot)$ is the three-parameter Mittag-Leffler
function defined in (2.9).

###### Proof.

From (4.2), we get

$\bar{p}^{f_{2}}(n,t)=\int_{0}^{\infty}p^{k}(n,x)l_{\eta,\nu}(x,t)\,\mathrm{d}x,\
\ t\geq 0,$ (4.3)

where $p^{k}(n,x)$ and $l_{\eta,\nu}(x,t)$ are the pmf of SPoK and the pdf of
$\mathscr{L}_{\eta,\nu}(t)$, respectively.

The following result will be used (See Kumar et al. (2019), Eq. (25)):

$\frac{\partial}{\partial
x}l_{\eta,\nu}(x,t)=-\left(\eta+\frac{\partial}{\partial
t}\right)^{\nu}l_{\eta,\nu}(x,t)+\eta^{\nu}l_{\eta,\nu}(x,t)-t^{-\nu}E^{1-\nu}_{1,1-\nu}(-\eta
t)\delta_{0}(x),$

where $\delta_{0}(x)=l_{\eta,\nu}(x,0)$. Also, we have
$\lim_{x\to\infty}l_{\eta,\nu}(x,t)=0$ (see Alrawashdeh et al. (2017), Lemma
4.6). Using the above results in (4.3), we get

$\displaystyle\left(\eta+\frac{\mathrm{d}}{\mathrm{d}t}\right)^{\nu}\bar{p}^{f_{2}}(n,t)$
$\displaystyle=\int_{0}^{\infty}p^{k}(n,x)\left(-\frac{\partial}{\partial
x}l_{\eta,\nu}(x,t)+\eta^{\nu}l_{\eta,\nu}(x,t)-t^{-\nu}E^{1-\nu}_{1,1-\nu}(-\eta
t)\delta_{0}(x)\right)\,\mathrm{d}x$
$\displaystyle=p^{k}(n,x)l_{\eta,\nu}(x,t)\big{|}_{x=0}+\int_{0}^{\infty}l_{\eta,\nu}(x,t)\frac{\mathrm{d}}{\mathrm{d}x}p^{k}(n,x)\,\mathrm{d}x+\eta^{\nu}\bar{p}^{f_{2}}(n,t)$
$\displaystyle\ \ -t^{-\nu}E^{1-\nu}_{1,1-\nu}(-\eta
t)\int_{0}^{\infty}p^{k}(n,x)\delta_{0}(x)\,\mathrm{d}x$
$\displaystyle=p^{k}(n,x)l_{\eta,\nu}(x,t)\big{|}_{x=0}+\int_{0}^{\infty}\bigg{(}-k(\lambda_{1}+\lambda_{2})p^{k}(n,x)+\lambda_{1}\sum_{j=1}^{k}p^{k}(n-j,x)$
$\displaystyle\ \
+\lambda_{2}\sum_{j=1}^{k}p^{k}(n+j,x)\bigg{)}l_{\eta,\nu}(x,t)\,\mathrm{d}x+\eta^{\nu}\bar{p}^{f_{2}}(n,t)-t^{-\nu}E^{1-\nu}_{1,1-\nu}(-\eta
t)p^{k}(n,0)$
$\displaystyle=p^{k}(n,x)l_{\eta,\nu}(x,t)\big{|}_{x=0}-k(\lambda_{1}+\lambda_{2})\bar{p}^{f_{2}}(n,t)+\lambda_{1}\sum_{j=1}^{k}\bar{p}^{f_{2}}(n-j,t)$
$\displaystyle\ \
+\lambda_{2}\sum_{j=1}^{k}\bar{p}^{f_{2}}(n+j,t)+\eta^{\nu}\bar{p}^{f_{2}}(n,t)-t^{-\nu}E^{1-\nu}_{1,1-\nu}(-\eta
t)p^{k}(n,0).$

This gives the required result. ∎

###  4.2. SPoK time-changed by the first hitting time of IG subordinator

Consider an IG subordinator whose associated Bernštein function $f_{3}$ is
given in (3.15). The first hitting time $\\{H(t)\\}_{t\geq 0}$ of the IG
subordinator $\\{Y(t)\\}_{t\geq 0}$ is defined as

$H(t)\coloneqq\inf\\{s\geq 0:Y(s)>t\\},\ \ t\geq 0.$

We define the following time-changed process:

$\bar{\mathcal{Z}}^{f_{3}}(t)\coloneqq\mathcal{S}^{k}(H(t)),\ \ t\geq 0,$
(4.4)

where the SPoK $\\{\mathcal{S}^{k}(t)\\}_{t\geq 0}$ is independent of
$\\{H(t)\\}_{t\geq 0}$.

######  Proposition 4.2.

For $n\in\mathbb{Z}$, the pmf
$\bar{p}^{f_{3}}(n,t)=\mathrm{Pr}\\{\bar{\mathcal{Z}}^{f_{3}}(t)=n\\}$ solves
the following differential equation:

$\displaystyle\delta\left(\gamma^{2}+2\frac{\mathrm{d}}{\mathrm{d}t}\right)^{1/2}\bar{p}^{f_{3}}(n,t)$
$\displaystyle=\left(\delta\gamma-k(\lambda_{1}+\lambda_{2})\right)\bar{p}^{f_{3}}(n,t)+\lambda_{1}\sum_{j=1}^{k}\bar{p}^{f_{3}}(n-j,t)$
$\displaystyle\ \
+\lambda_{2}\sum_{j=1}^{k}\bar{p}^{f_{3}}(n+j,t)-\delta\gamma\mathrm{Erf}\left(\frac{\gamma\sqrt{t}}{\sqrt{2}}\right)p^{k}(n,0),$

where $\mathrm{Erf}(\cdot)$ is the error function.

###### Proof.

From (4.4), we have

$\bar{p}^{f_{3}}(n,t)=\int_{0}^{\infty}p^{k}(n,x)h(x,t)\,\mathrm{d}x,\ \ t\geq
0,$ (4.5)

where $h(x,t)$ is the pdf of $H(t)$.

The following result will be used (see Wylomańska et al. (2016), Eq. (2.22))

$\frac{\partial}{\partial
x}h(x,t)=-\delta\left(\gamma^{2}+2\frac{\partial}{\partial
t}\right)^{1/2}h(x,t)+\delta\gamma
h(x,t)-\frac{\delta\sqrt{2}e^{-\gamma^{2}t/2}}{\sqrt{\pi t}}\delta_{0}(x),$

with the initial condition $h(x,0)=\delta_{0}(x)$. Using the above result in
(4.5), we get

$\displaystyle\delta\left(\gamma^{2}+2\frac{\mathrm{d}}{\mathrm{d}t}\right)^{1/2}\bar{p}^{f_{3}}(n,t)$
$\displaystyle=\int_{0}^{\infty}p^{k}(n,x)\left(-\frac{\partial}{\partial
x}h(x,t)+\delta\gamma h(x,t)-\frac{\delta\sqrt{2}e^{-\gamma^{2}t/2}}{\sqrt{\pi
t}}\delta_{0}(x)\right)\mathrm{d}x$
$\displaystyle=p^{k}(n,0)h(0,t)+\int_{0}^{\infty}h(x,t)\frac{\mathrm{d}}{\mathrm{d}x}p^{k}(n,x)\mathrm{d}x+\delta\gamma\bar{p}^{f_{3}}(n,t)$
$\displaystyle\ \ -\frac{\delta\sqrt{2}e^{-\gamma^{2}t/2}}{\sqrt{\pi
t}}p^{k}(n,0)$
$\displaystyle=p^{k}(n,0)h(0,t)+\int_{0}^{\infty}\bigg{(}-k(\lambda_{1}+\lambda_{2})p^{k}(n,x)+\lambda_{1}\sum_{j=1}^{k}p^{k}(n-j,x)$
$\displaystyle\ \
+\lambda_{2}\sum_{j=1}^{k}p^{k}(n+j,x)\bigg{)}h(x,t)\mathrm{d}x+\delta\gamma\bar{p}^{f_{3}}(n,t)-\frac{\delta\sqrt{2}e^{-\gamma^{2}t/2}}{\sqrt{\pi
t}}p^{k}(n,0)$
$\displaystyle=p^{k}(n,0)h(0,t)-k(\lambda_{1}+\lambda_{2})\bar{p}^{f_{3}}(n,t)+\lambda_{1}\sum_{j=1}^{k}\bar{p}^{f_{3}}(n-j,t)$
$\displaystyle\ \
+\lambda_{2}\sum_{j=1}^{k}\bar{p}^{f_{3}}(n+j,t)+\delta\gamma\bar{p}^{f_{3}}(n,t)-\frac{\delta\sqrt{2}e^{-\gamma^{2}t/2}}{\sqrt{\pi
t}}p^{k}(n,0)$
$\displaystyle=\left(\delta\gamma-k(\lambda_{1}+\lambda_{2})\right)\bar{p}^{f_{3}}(n,t)+\lambda_{1}\sum_{j=1}^{k}\bar{p}^{f_{3}}(n-j,t)$
$\displaystyle\ \
+\lambda_{2}\sum_{j=1}^{k}\bar{p}^{f_{3}}(n+j,t)-\delta\gamma\mathrm{Erf}\left(\frac{\gamma\sqrt{t}}{\sqrt{2}}\right)p^{k}(n,0).$

The last step follows by using the following result (see Vellaisamy and Kumar
(2018), Proposition 2.2):

$\lim\limits_{x\to 0}h(x,t)=h(0,t)=\delta
e^{-\gamma^{2}t/2}\left(\sqrt{\frac{2}{\pi t}}-\gamma
e^{\gamma^{2}t/2}\mathrm{Erf}\left(\frac{\gamma\sqrt{t}}{\sqrt{2}}\right)\right).$

This completes the proof. ∎

##  5\. Concluding remarks

In this paper, we introduced and studied a fractional version of the SPoK,
namely, the FSPoK. It is obtained by time-changing the SPoK with an
independent inverse stable subordinator. Its one-dimensional distribution,
pgf, mean, variance and covariance are obtained. It is shown that the FSPoK
exhibits the LRD property. Also, we considered two time-changed versions of
the FSPoK where time-change is done using a Lévy subordinator and its inverse.
Some distributional properties and particular cases are discussed for these
time-changed processes. Kerss et al. (2014) introduced and studied a
fractional version of the Skellam process, namely, the FSP. They showed its
applications to high frequency financial data set. As the FSPoK is a
generalized version of the FSP, we expect the FSPoK to have potential
applications in finance and related fields.

## References

  * [1] Abramowitz, M. and Stegun, I. A. (1972). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. National Bureau of Standards Applied Mathematics Series, vol. 55, 10th printing.
  * [2] Aletti, G., Leonenko, N. and Merzbach, E. (2018). Fractional Poisson fields and martingales. J. Stat. Phys. 170(4), 700-730.
  * [3] Alrawashdeh, M. S., Kelly, J. F., Meerschaert, M. M. and Scheffler, H. -P. (2017). Applications of inverse tempered stable subordinators. Comput. Math. Appl. 73(6), 892-905.
  * [4] Applebaum, D. (2009). Lévy Processes and Stochastic Calculus. Cambridge University Press, Cambridge.
  * [5] Barndorff-Nielsen, O. E., Pollard, D. G. and Shephard, N. (2012). Integer-valued Lévy processes and low latency financial econometrics. Quant. Finance 12(4), 587-605.
  * [6] Beghin, L. (2012). Random-time processes governed by differential equations of fractional distributed order. Chaos Solitons Fractals 45(11), 1314-1327.
  * [7] Beghin, L. (2015). On fractional tempered stable processes and their governing differential equations. J. Comput. Phys. 293, 29-39.
  * [8] Bertoin, J. (1996). Lévy Processes. Cambridge University Press, Cambridge.
  * [9] D’Ovidio, M. and Nane, E. (2014). Time dependent random fields on spherical non-homogeneous surfaces. Stochastic Process. Appl. 124(6), 2098-2131.
  * [10] Gupta, N., Kumar, A. and Leonenko, N. (2020). Skellam type processes of order $K$ and beyond. Entropy. 22(11).
  * [11] Hwang, Y., Kim, J. and Kweon, I. (2007). Sensor noise modeling using the Skellam distribution: Application to the color edge detection. In Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, 1-8.
  * [12] Karlis, D. and Ntzoufras, I. ( 2009). Bayesian modelling of football outcomes: Using the Skellam’s distribution for the goal difference. IMA J. Manag. Math. 20, 133-145.
  * [13] Kerss, A., Leonenko, N. N. and Sikorskii, A. (2014). Fractional Skellam processes with applications to finance. Fract. Calc. Appl. Anal. 17(2), 532-551.
  * [14] Kostadinova, K. Y. and Minkova, L. D. (2013). On the Poisson process of order $k$. Pliska Stud. Math. Bulgar. 22, 117-128.
  * [15] Kilbas, A. A., Srivastava, H. M. and Trujillo, J. J. (2006). Theory and Applications of Fractional Differential Equations. Elsevier, Amsterdam.
  * [16] Kumar, A., Gajda, J., Wylomańska, A. and Poloczański, R. (2019). Fractional Brownian motion delayed by tempered and inverse tempered stable subordinators. Methodol. Comput. Appl. Probab. 21(1), 185-202.
  * [17] Leonenko, N. N., Meerschaert, M. M., Schilling, R. L. and Sikorskii, A. (2014). Correlation structure of time-changed Lévy processes. Commun. Appl. Ind. Math. 6(1), e-483, 22.
  * [18] Mainardi, F. (2010). Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models. Imperial College Press, London.
  * [19] Maheshwari, A. and Vellaisamy, P. (2016). On the long-range dependence of fractional Poisson and negative binomial processes. J. Appl. Probab. 53(4), 989-1000.
  * [20] Maheshwari, A. and Vellaisamy, P. (2019). Fractional Poisson process time-changed by Lévy subordinator and its inverse. J. Theoret. Probab. 32(3), 1278-1305.
  * [21] Meerschaert, M. M., Schilling, R. L. and Sikorskii, A. (2015). Stochastic solutions for fractional wave equations. Nonlinear Dynam. 80(4), 1685-1695.
  * [22] Meerschaert, M. M., Nane, E. and Vellaisamy, P. (2011). The fractional Poisson process and the inverse stable subordinator. Electron. J. Probab. 16(59), 1600-1620.
  * [23] Meerschaert, M. M. and Straka, P. (2013). Inverse stable subordinators. Math. Model. Nat. Phenom. 8(2), 1-16.
  * [24] Orsingher, E. and Polito, F. (2012). The space-fractional Poisson process. Statist. Probab. Lett. 82(4), 852-858.
  * [25] Orsingher, E. and Toaldo, B. (2015). Counting processes with Bernštein intertimes and random jumps. J. Appl. Probab. 52(4), 1028-1044.
  * [26] Sengar, A. S., Maheshwari, A. and Upadhye, N. S. (2020). Time-changed Poisson processes of order k. Stoch. Anal. Appl. 38(1), 124-148.
  * [27] Skellam, J. G. (1946). The frequency distribution of the difference between two Poisson variates belonging to different populations. J. Roy. Statist. Soc. (N.S.) 109, 296.
  * [28] Steutel, F. W. and van Harn, K. (2004). Infinite Divisibility of Probability Distributions on the Real Line. Marcel Dekker, New York.
  * [29] Vellaisamy, P. and Kumar, A. (2018). First-exit times of an inverse Gaussian process. Stochastics 90(1), 29-48.
  * [30] Vellaisamy, P. and Maheshwari, A. (2018). Fractional negative binomial and Pólya processes. Probab. Math. Statist. 38(1), 77-101.
  * [31] Wylomańska, A., Kumar, A., Poloczański, R. and Vellaisamy, P. (2016). Inverse Gaussian and its inverse process as the subordinators of fractional Brownian motion. Phys. Rev. E 94, 042128.