Fractional non-homogeneous Poisson and Pólya-Aeppli processes of order k and beyond

Abstract

We introduce two classes of point processes: a fractional non-homogeneous Poisson process of order k and a fractional non-homogeneous Pólya-Aeppli process of order k. We characterize these processes by deriving their non-local governing equations. We further study the covariance structure of the processes and investigate the long-range dependence property.

Keywords:

• non-homogeneous fractional Poisson process of order k non-homogeneous fractional Pólya-Aeppli process of order k long range dependence
• Caputo fractional derivative
• α-stable Lévy subordinators
• fractional integro-differential difference equations

AMS classification:

• 60G55
• 26A33
• 60G05
• 60G51

1. Introduction

Fractional Poisson processes (FPP) enjoy the property of non-stationarity and long range dependence, which makes them an attractive modeling tool. These processes are widely used in statistics, finance, meteorology, physics and network science, see for instance (Baleanu et al. 2012) p. 332, and (Kumar, Leonenko, and Pichler 2020).

Fractional Poisson processes were introduced as renewal processes in Mainardi, Gorenflo, and Scalas (2004). The authors generalized the characterization of the Poisson process as the counting process for epochs defined as sum of independent non-negative exponential random variables, and, instead of the exponential, the authors used a Mittag-Leffler distribution. The theory of FPP was further developed by Beghin and Orsingher (2009, 2010) and by Meerschaert, Nane, and Vellaisamy (2011).

In particular, Meerschaert, Nane, and Vellaisamy (2011) defined FPP by means of a time-change for the Poisson process N(t), where the time variable t is replaced by the inverse α-stable subordinator $Y_{\alpha }(t).$ Remarkably, they could prove the equality in distribution between $N(Y_{\alpha }(t))$ and the counting process defined by (Mainardi, Gorenflo, and Scalas 2004).

Leonenko, Scalas, and Trinh (2017) used the same time-change technique to introduce a non-homogeneous fractional Poisson process (NFPP) by replacing the time variable in the FPP with an appropriate function of time.

In a recent paper, Gupta, Kumar, and Leonenko (2020) and Gupta and Kumar (2021) generalize the results available on fractional Poisson processes using the z-transform technique.

Kostadinova and Minkova (2019) introduced a Poisson process of order $k$ with insurance modeling in mind. This process models the claim arrival in groups of size k, where the number of arrivals in a group is uniformly distributed over k points.

The Pólya-Aeppli process of order k was studied in Chukova and Minkova (2015) and later by Kostadinova and Lazarova (2019). In this process, the uniform distribution on the integers $1,\dots ,k$ is replaced by the truncated geometric distribution of parameter ρ.

To deal with dependent inter-arrival times, a generalization of Poisson processes of order k was proposed by Sengar, Maheshwari, and Upadhye (2020). These authors extended the Poisson process of order k by means of time change with a general Lévy subordinator as well as an inverse Lévy subordinator.

Here, we combine the compound Poisson processes of order k and fractional Poisson processes, namely we study a fractional non-homogeneous Poisson process of order k and a fractional non-homogeneous Pólya-Aeppli process of order k (see the definitions below). First, we generalize the results of Kostadinova and Minkova (2019) by considering a non-homogeneous Poisson process of order k. Then, we generalize the results of Sengar, Maheshwari, and Upadhye (2020) by introducing the time non-homogeneity in the fractional Poisson process of order k. Finally, we study a non-homogeneous fractional Pólya-Aeppli process of order k.

This paper is organized as follows. Section 2 collects some known results from the theory of subordinators and provides the definition of the compound distributions of order k. In Section 3, we consider a non-homogeneous fractional Poisson process of order k. We obtain the governing equations and calculate the moments and the covariance function of the process. Section 4 is devoted to a non-homogeneous fractional Pólya-Aeppli process of order k. We derive the non-local governing equations for the marginal distributions of these processes, using non-local operators known as Caputo derivatives. The moments and the covariance structure of the processes are derived, as well.

2. Preliminaries

This section presents known results in the theory of subordinators and provides the definition of the compound distributions of order k.

2.1. Compound distributions of order k

Consider a random variable that can be represented as a random sum $N=X_1+X_2+\dots +X_Y,$ where ${\left\{X_i\right\}}_{i=1}^{\infty }$ is a sequence of independent identically distributed random variables (i.i.d. r.v’s), independent of a non-negative integer-valued random variable Y. The probability distribution of N is called compound distribution and the distribution of X₁ is called compounding distribution.

A well-known and widely used example is the compound Poisson distribution, where Y has a Poisson distribution. If $X_i\in \{1,2,\dots ,k\},$ then the random variable N has a compound discrete distribution of order k.

Compound discrete distributions of order k were studied by Philippou (1983) and Philippou, Georghiou, and Philippou (1983).

As mentioned previously, in this paper we deal with two types of compounding distributions: the discrete uniform distribution and the truncated geometric distribution. They respectively induce the Poisson distribution of order k and the Pólya-Aeppli distribution of order k as will be shown in the following. We say that the random variable X is uniformly distributed over k points if its probability mass function (pmf) is of the form (1) $$\mathbb{P}[X=m]=\frac{1}{k},\hspace{1em}m=1,\dots ,k.$$(1)

Its probability generating function (pgf) is $$G_X(u)=\mathbb{E}\left[u^X\right]=\frac{1}{k}(u+u^2+\dots +u^k)=\frac{u}{k}·\frac{1-u^k}{1-u},\hspace{1em}u\in (0,1).$$

The random variable X has a truncated geometric distribution with parameter ρ and with success probability $1-\rho$ if (2) $$\mathbb{P}[X=m]=\frac{1-\rho }{1-{\rho }^k}{\rho }^{m-1},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}m=1,2,\dots ,k,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\rho \in [0,1).$$(2)

Consequently, the pgf of X is given by (3) $$G_X(u)=\mathbb{E}\left[u^X\right]=\frac{(1-\rho )u}{1-{\rho }^k}\frac{1-{\rho }^k{u}^k}{1-\rho u},\hspace{1em}u\in (0,1).$$(3)

Note, that for $k\to \infty ,$ the truncated geometric distribution asymptotically coincides with the geometric distribution with parameter $1-\rho .$

We can now define the Poisson distribution of order k.

Definition 1

(Poisson distribution of order k). The random variable N has Poisson distribution of order k with parameter Λ if $N=X_1+X_2+\dots +X_Y,$ where:

(1) ${\left\{X_i\right\}}_{i\ge 1}$ are the i.i.d. r.v’s with the uniform distribution; (2) Y has Poisson distribution with parameter $\Lambda >0;$ (3) Y and ${\left\{X_i\right\}}_{i\ge 1}$ are independent.

Note that $$\mathbb{P}[N=m]=e^{-\Lambda k}\sum _{(n_1,\dots n_k)\in \Omega (k,m)}\frac{{\Lambda }^{n_1+\dots +n_k}}{n_1!·\dots ·n_k!}=e^{-\Lambda k}\sum _{\Omega (k,m)}\frac{{\Lambda }^{z_k}}{{\Pi }_k!},$$ where $z_k=n_1+n_2+\dots +n_k,\hspace{1em}{\Pi }_k!=n_1!·n_2!·\dots ·n_k!,$ and (4) $$\Omega (k,m)=\{(n_1,\dots n_k):n_1+2n_2+\dots k{n}_k=m\}.$$(4)

The pgf of the Poisson distribution of order k is (5) $$G_N(u)=\mathbb{E}\left[u^N\right]=\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\Lambda \left(k-\sum _{j=1}^k{u}^j\right)\right\}.$$(5)

Note that $N\stackrel{d}{=}{\sum }_{j=1}^{k}j{Y}_j,$ where $Y_j,j=1,\dots ,k$ are independent copies of Poisson random variable Y with parameter Λ, and “$\stackrel{d}{=}$“stands for equality in distributions. We now introduce the Pólya-Aeppli distribution of order k as a compound Poisson distribution (see (Minkova 2010)).

Definition 2

(Pólya-Aeppli distribution of order k). The random variable N has Pólya-Aeppli distribution of order k with parameter $1-\rho$ if $N=X_1+X_2+\dots +X_Y,$ where: (i) ${\left\{X_i\right\}}_{i\ge 1}$ are the i.i.d. r.v’s with the truncated geometric distribution with parameter $1-\rho ,$ given by (2); (ii) Y has Poisson distribution with parameter $\Lambda ;$ (iii) Y and ${\left\{X_i\right\}}_{i\ge 1}$ are independent.

Note that the probability generating function of N is $G_N(u)={\mathrm{e}}^{-\Lambda (1-G_{X_1}(u))},$ where $G_{X_1}$ is given by (3).

The probability mass function of Pólya-Aeppli distribution of order k is defined by (see Minkova 2010, Theorem 3.1): (6) $$\mathbb{P}[N=m]=q_m(\Lambda ),m=0,1,2,\dots ,$$(6) where $$\begin{array}{c}{q}_0(\Lambda )=e^{-\Lambda }\\ q_m(\Lambda )=e^{-\Lambda }\sum _{j=1}^m\left(\begin{array}{c}m-1\\ j-1\end{array}\right)\frac{Q^J}{j!}{\rho }^{m-j},m=1,2,\dots ,k\\ q_m(\Lambda )=e^{-\Lambda }[\sum _{j=1}^m\left(\begin{array}{c}m-1\\ j-1\end{array}\right)\frac{Q^J}{j!}{\rho }^{m-j}-\sum _{n=1}^l{(-1)}^{n-1}\frac{{(Q{\rho }^k)}^n}{n!}\times \\ \times \sum _{j=0}^{m-n(k+1)}\left(\begin{array}{c}m-n(k+1)+n-1\\ j+n-1\end{array}\right)\frac{Q^J}{j!}{\rho }^{m-j-n(k+1)}],\end{array}$$ and $$Q=\frac{\Lambda (1-\rho )}{1-{\rho }^k},m=l(k+1)+r,r=0,1,..,k,l=1,2,\dots .$$

2.2. Inverse $\alpha$-stable subordinator

Let $L_{\alpha }=\{L_{\alpha }(t);t\ge 0\}$ be a $\alpha$-stable Lévy subordinator, that is Lévy process with Laplace transform: $$\mathbb{E}\left[{\mathrm{e}}^{-s{L}_{\alpha }(t)}\right]=e^{-t{s}^{\alpha }},\hspace{1em}0<\alpha <1,\hspace{1em}s\ge 0.$$

Then the inverse α-stable subordinator $\{Y_{\alpha }(t);t\ge 0\}$ (see e.g., Meerschaert and Sikorskii 2019, 103) is defined as the first passage time of $L_{\alpha }:$ (7) $$Y_{\alpha }(t)=\inf \{u>0:L_{\alpha }(u)>t\},\hspace{1em}t\ge 0.$$(7)

We will use the following properties of the inverse α-stable subordinator:

1. The density of $Y_{\alpha }(t)$ is of the form (see (Meerschaert and Sikorskii 2019) p.113): (8) $$h_{\alpha }(t,x)=\frac{d}{dx}\mathbb{P}[Y_{\alpha }(t)\le x]=\frac{t}{\alpha }{x}^{-1-\frac{1}{\alpha }}{g}_{\alpha }\left(t{x}^{-\frac{1}{\alpha }}\right),\hspace{1em}x>0,\hspace{1em}t>0,$$(8)

   where $$g_{\alpha }(x)=\frac{1}{\pi }\sum _{k=1}^{\infty }{(-1)}^{k+1}\frac{\Gamma (\alpha k+1)}{k!}\frac{1}{x^{\alpha k+1}}\text{sin}\hspace{0.17em}(\pi k\alpha )$$ is the density of $L_{\alpha }(1)$ (see e.g. (Kataria and Vellaisamy 2018)).

2. The Laplace transform (9) $${\tilde{h}}_{\alpha }(s,x)=\underset{0}{\overset{\infty }{\int }}{\mathrm{e}}^{-st}{h}_{\alpha }(t,x)dt=s^{\alpha -1}{e}^{-x{s}^{\alpha }},\hspace{1em}s\ge 0.$$(9)

3. The moments of the inverse α-stable subordinator are as follows: (10) $$\mathbb{E}[Y_{\alpha }^{\nu }(t)]=\frac{\Gamma (\nu +1)}{\Gamma (\alpha \nu +1)}{t}^{\alpha \nu },\nu >0,\hspace{1em}\text{Var}[Y_{\alpha }(t)]=t^{2\alpha }\left[\frac{2}{\Gamma (2\alpha +1)}-\frac{1}{{(\Gamma (\alpha +1))}^2}\right].$$(10)

   (see e.g. (Kataria and Vellaisamy 2018, 1640).

4. The covariance function (see (Leonenko et al. 2014; Leonenko, Scalas, and Trinh 2017)) is (11) $$\text{Cov}[Y_{\alpha }(t),Y_{\alpha }(s)]=\frac{1}{\Gamma (1+\alpha )\Gamma (\alpha )}\underset{0}{\overset{\min (t,s)}{\int }}((t-\tau {)}^{\alpha }+{(s-\tau )}^{\alpha }){\tau }^{\alpha -1}d\tau -\frac{{(st)}^{\alpha }}{{\Gamma }^2(1+\alpha )}.$$(11)

3. Poisson processes of order k

The Poisson process of order k was introduced in Kostadinova and Minkova (2019), see also (Sengar, Maheshwari, and Upadhye 2020).

Definition 3.

The Poisson process of order k (PPk) $N=\{N(t);t\ge 0\}$ is defined as a compound Poisson process with the compounding discrete uniform distribution: (12) $$N(t)=X_1+\dots +X_{N_1(t)},$$(12) where (1) X_i are independent copies of a discrete uniform random variable distributed over k points given by (1); (2) $N_1=\{N_1(t);t\ge 0\}$ is the Poisson process with parameter $k\lambda ;$ (3) N₁ and ${\left\{X_i\right\}}_{i\ge 1}$ are independent.

The following Kolmogorov forward equations are valid for $p_m(t)=\mathbb{P}[N(t)=m]:$ (13) $$\frac{d}{dt}{p}_0(t)=-k\lambda p_0(t)$$(13) (14) $$\frac{d}{dt}{p}_m(t)=-k\lambda p_m(t)+\lambda \sum _{j=1}^{m\wedge k}{p}_{m-j}(t),\hspace{1em}m=1,2,\dots$$(14) with the initial condition $p_0(0)=1,:p_m(0)=0,m\ge 1,$ and $m\wedge k=\min (m,k).$ The pgf is of the form: $$G_{N(t)}(u)=\mathbb{E}\left[u^{N(t)}\right]=\hspace{0.17em}\text{exp}\hspace{0.17em}\{\lambda t(u+\dots +u^k-k)\},$$ and the first two moments are given by (15) $$\mathbb{E}[N(t)]=\frac{k(k+1)}{2}\lambda t,\hspace{1em}\text{Cov}[N(t),N(s)]=\frac{k(k+1)(2k+1)}{6}\lambda \min (s,t).$$(15)

3.1. Fractional Poisson process of order k

In this sub-section we shall derive governing equations for a fractional Poisson process of order k and we shall investigate its long-range dependence properties. It is worth noting that Sengar, Maheshwari, and Upadhye (2020) studied the Poisson process of order k time-changed by a general Lévy subordinator and its inverse. However, among their examples, they did not explicitly consider the governing equations for the inverse α-stable subordinator (this particular process is studied in Gupta and Kumar (2021)). That is why we specify some formulae of Sengar, Maheshwari, and Upadhye (2020) that will be used in the following sub-sections. In particular, below, we use Equation (10)(10) $$\mathbb{E}[Y_{\alpha }^{\nu }(t)]=\frac{\Gamma (\nu +1)}{\Gamma (\alpha \nu +1)}{t}^{\alpha \nu },\nu >0,\hspace{1em}\text{Var}[Y_{\alpha }(t)]=t^{2\alpha }\left[\frac{2}{\Gamma (2\alpha +1)}-\frac{1}{{(\Gamma (\alpha +1))}^2}\right].$$(10) to derive the marginal distributions of the fractional Poisson process of order k.

Definition 4.

(Fractional Poisson process of order k). The process $N_{\alpha }(t)$ is called fractional Poisson process of order k (FPPk) if (16) $$N_{\alpha }(t)=N(Y_{\alpha }(t)),\hspace{1em}0<\alpha <1,$$(16) where (1) $Y_{\alpha }(t)$ is the inverse α-stable subordinator, given by (7); (2) N is the Poisson process of order k, given by (12); (3) $Y_{\alpha }(t)$ and N are independent.

The marginal distributions of the FPPk process is given by $$\begin{array}{c}{p}_m^{\alpha }(t)=\mathbb{P}[N_{\alpha }(t)=m]=\sum _{\Omega (k,m)}\frac{{\lambda }^{z_k}}{{\Pi }_k!}\sum _{n=0}^{\infty }\frac{{(-k\lambda )}^n}{n!}\mathbb{E}\left[{(Y_{\alpha }(t))}^{z_k+n}\right]=\\ =\sum _{\Omega (k,m)}\frac{{\lambda }^{z_k}}{{\Pi }_k!}\sum _{n=0}^{\infty }\frac{{(-k\lambda )}^n}{n!}\frac{\Gamma (z_k+n+1)}{\Gamma (\alpha (z_k+n)+1)}{t}^{\alpha (z_k+n)},m=0,1,..\end{array}$$ where $z_k=n_1+n_2+\dots +n_k,\hspace{1em}{\Pi }_k!=n_1!n_2!\dots n_k!,$ and $\Omega (k,m)$ is defined in (4).

Also $$\begin{array}{c}\mathbb{E}[N_{\alpha }(t)]=k\lambda \frac{(k+1)}{2}\mathbb{E}[Y_{\alpha }(t)],\\ \text{Var}[N_{\alpha }(t)]=k\lambda \frac{(k+1)(2k+1)}{6}\mathbb{E}[Y_{\alpha }(t)]+{\left(k\lambda \frac{(k+1)}{2}\right)}^2\text{Var}(Y_{\alpha }(t)),\\ \text{Cov}[N_{\alpha }(t),N_{\alpha }(s)]=\frac{k(k+1)(2k+1)\lambda {(\min (t,s))}^{\alpha }}{6\Gamma (1+\alpha )}+{\left(\frac{k\lambda (k+1)}{2}\right)}^2\text{Cov}(Y_{\alpha }(s),Y_{\alpha }(t)),\end{array}$$ where the variance and covariance of the process $Y_{\alpha }(t)$ are given by (10) and (11).

3.1.1. Correlation structure and long-range dependence

There exist many definitions of the long-range dependence property. Here, we shall use the definition given in Biard and Saussereau (2014).

Definition 5.

The process $\{X(t);t\ge 0\}$ has a long-range dependence property (LRD) if for fixed s and some c(s) and $\alpha \in (0,1):\underset{t\to \infty }{\lim }[\text{Corr}(X(s),X(t))/t^{-\alpha }]=c(s),$ where $\text{Corr}$ is the correlation function of the process $X.$

We now investigate the asymptotic behavior of the correlation function of the FPPk process defined by (16).

Theorem 3.1.

The process $N_{\alpha }(t)$ has the LRD property.

Proof.

Using the result of Leonenko et al. (2014) we have that for a fixed s > 0 $$\text{Corr}[N_{\alpha }(t),N_{\alpha }(s)]\sim t^{-\alpha }C(\alpha ,s)\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}t\to \infty ,$$ where $C(\alpha ,s)={(\frac{1}{\Gamma (2\alpha )}-\frac{1}{\alpha {(\Gamma (\alpha ))}^2})}^{-1}\left[\frac{\alpha \text{Var}[N(1)]}{\Gamma (1+\alpha ){(\mathbb{E}[N(1)])}^2}+\frac{\alpha s^{\alpha }}{\Gamma (1+2\alpha )}\right],$ and $\mathbb{E}[N(1)]$ and $\text{Var}[N(1)]$ are given by (15).

3.1.2. Governing equations

In the sequel we will employ the fractional Caputo (or Caputo-Djrbashian) derivative which is defined as follows (see e.g., (Meerschaert and Sikorskii 2019, 30) (17) $$D_t^{\alpha }f(t)=\{\begin{array}{cc}\hfill & \frac{1}{\Gamma (\alpha )}\underset{0}{\overset{t}{\int }}\frac{df(u)}{du}\frac{du}{{(t-u)}^{\alpha }},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}0<\alpha <1,\hfill \\ \hfill & \frac{df(u)}{du},\alpha =1.\hfill \end{array}$$(17)

Theorem 3.2.

The governing fractional difference-differential equations for

$p_m^{\alpha }(t),t\ge 0$ are given by (18) $$D_t^{\alpha }{p}_0^{\alpha }(t)=-k\lambda p_0^{\alpha }(t)$$(18) (19) $$D_t^{\alpha }{p}_m^{\alpha }(t)=-k\lambda p_m^{\alpha }(t)+\lambda \sum _{j=1}^{m\wedge k}{p}_{m-j}^{\alpha }(t),\hspace{1em}m=1,2,\dots$$(19) with the initial condition $$p_m^{\alpha }(0)={\delta }_{m,0}=\{\begin{array}{c}1,\hspace{1em}m=0\hfill \\ 0,\hspace{1em}m\ge 1.\hfill \end{array}$$

Note, that by setting $\alpha =1,$ we get the governing equations of the Poisson process of order k given in Equation (13)(13) $$\frac{d}{dt}{p}_0(t)=-k\lambda p_0(t)$$(13) .

Proof.

Note that (20) $$D_t^{\alpha }{h}_{\alpha }(t,u)=-\frac{\partial }{\partial u}{h}_{\alpha }(t,u)$$(20) and remember that (21) $$p_n^{\alpha }(t)=\underset{0}{\overset{\infty }{\int }}{p}_n(u)h_{\alpha }(t,u)du\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}n=0,1,2\dots$$(21)

We first consider the case $n\ge 1.$ By taking the fractional Caputo derivative of both sides (21) and using property (20), we get $$\begin{array}{c}{D}_t^{\alpha }{p}_m^{\alpha }(t)=-\underset{0}{\overset{\infty }{\int }}{p}_m(u)\frac{\partial }{\partial u}{h}_{\alpha }(t,u)du=\\ =\underset{0}{\overset{\infty }{\int }}[-k\lambda p_m(u)+\lambda \sum _{j=1}^{m\wedge k}{p}_{m-j}(u)]h_{\alpha }(t,u)du-p_m(u)h_{\alpha }(t,u){|}_0^{\infty }=\\ =-k\lambda p_m^{\alpha }(t)+\lambda \sum _{j=1}^{m\wedge k}{p}_{m-j}^{\alpha }(t).\end{array}$$

For n = 0 we have $$D_t^{\alpha }{p}_0^{\alpha }(t)=-\underset{0}{\overset{\infty }{\int }}{p}_0(u)\frac{\partial }{\partial u}{h}_{\alpha }(t,u)du=\underset{0}{\overset{\infty }{\int }}[-k\lambda p_0(u)]h_{\alpha }(t,u)du=-k\lambda p_0^{\alpha }(t).$$

Remark 1.

Note that Sengar, Maheshwari, and Upadhye (2020) derived governing equations in which the Caputo derivative is replaced by a more general non-local operator. We present the proof of Theorem 3.2 for the sake of completeness.

3.2. Non-homogeneous fractional Poisson process of order k

We now generalize the fractional Poisson process of order k by introducing a deterministic, time dependent intensity or rate function $\lambda (t):[0,\infty )\to [0,\infty ),$ such that for every fixed $t>0,$ the cumulative rate function follows the following equation $$\Lambda (t)=\underset{0}{\overset{t}{\int }}\lambda (u)du<\infty$$

Denote $\Lambda (s,t)={\int }_s^t\lambda (u)du=\Lambda (t)-\Lambda (s),0\le s<t.$ Let $N_1^1(t);t\ge 0$ be a homogeneous Poisson process (HPP) of unit intensity, and $N_1^1(\Lambda (t)),t\ge 0,$ be a non-homogeneous Poisson process (NPP) with rate function $\lambda (t),$ then $$N^n(t)=X_1+\dots +X_{N_1^1(k\Lambda (t))},\hspace{1em}t\ge 0,$$ is non-homogeneous Poisson process of order k (NPPk), with rate function $\lambda (t),t\ge 0,$ where ${\left\{X_i\right\}}_{i\ge 1}$ are the i.i.d.r.v’s with the uniform distribution on $\{1,2,\dots ,k\},$ independent of $N_1^1(\Lambda (t)).$ The mgf of the process Nⁿ is of the form: $$G_{N^n(t)}(u)=\mathbb{E}\left[u^{N^n(t)}\right]=\hspace{0.17em}\text{exp}\hspace{0.17em}\{\Lambda (t)(u+\dots +u^k-k)\}.$$

The process Nⁿ has the following distributions of its increments: (22) $$\begin{array}{c}{p}_m^n(t,u)=\mathbb{P}[N^n(t+u)-N^n(u)=m]=\\ =e^{-k\Lambda (u,t+u)}\sum _{\Omega (k,m)}\frac{{[\Lambda (u,u+t)]}^{n_1+\dots +n_k}}{n_1!\dots n_k!},\hspace{1em}m=0,1,\dots \end{array}$$(22)

Incidentally, this model includes Weibull’s rate function: $\Lambda (t):=\Lambda (0,t)={(\frac{t}{b})}^c,\hspace{1em}\lambda (t)=\frac{c}{b}{(\frac{t}{b})}^{c-1},\hspace{1em}c\ge 0,b>0;$ Makeham’s rate function: $\Lambda (t)=\frac{c}{b}{e}^{bt}-\frac{c}{b}+\mu t,\hspace{1em}\lambda (t)=c{e}^{bt}+\mu ,\hspace{1em}c>0,b>0,\mu \ge 0,$ and many others.

We define a non-homogeneous fractional Poisson process of order k (FNPPk) as (23) $$N_{\alpha }^{\ast }(t)=N^n(Y_{\alpha }(t)),\hspace{1em}t\ge 0,\hspace{1em}0<\alpha <1,$$(23) where $Y_{\alpha }(t)$ is the inverse α-stable subordinator (7), independent of the NPPk process $N^n.$

3.2.1. Marginal distributions

Define the increment process: $I_{\alpha }(t,v)=N(\Lambda (Y_{\alpha }(t)+v))-N(\Lambda (v)).$ Its marginal distributions can be written as follows: (24) $$p_m^{\ast }(t,v)=\mathbb{P}[I_{\alpha }(t,v)=m]=\underset{0}{\overset{\infty }{\int }}{p}_m^n(u,v)h_{\alpha }(t,u)du,$$(24) where $h_{\alpha }(t,u)$ is the density of the inverse α-stable subordinator (8) and $p_x^n(u,v)$ is given by (22). Consequently the marginal distributions of $N_{\alpha }^{\ast }(t)$ are given by $$\mathbb{P}[N_{\alpha }^{\ast }(t)=m]=p_m^{\ast }(t,0)=\underset{0}{\overset{\infty }{\int }}{p}_m^n(u,0)h_{\alpha }(t,u)du.$$

For the NFPP $N_1^1(\Lambda (Y_{\alpha }(t));t\ge 0,$ of order $k=1,$ Leonenko et al. (Leonenko, Scalas, and Trinh 2017) derived the governing equations for the marginal distributions $\mathbb{P}[I_{\alpha }^1(t,v)=m]$ of the corresponding increment process $I_{\alpha }^1(t,v)=N_1^1(\Lambda (Y_{\alpha }(t)+v))-N_1(\Lambda (v))$ of NFPP (of order k = 1), where $N_1^1$ is the homogeneous Poisson process of intensity 1. We shall derive the governing equations for the marginal distributions $p_x^{\ast }(t,v)$ of FNPPk.

Theorem 3.3.

The marginal distributions $p_x^{\ast }(t,v)$ satisfy the following fractional differential-difference integral equations (25) $$\begin{array}{c}{D}_t^{\alpha }{p}_0^{\ast }(u,v)=-k\underset{0}{\overset{\infty }{\int }}\lambda (u+v)p_0^n(u,v)h_{\alpha }(t,u)du\hspace{1em}0\le v<u\\ D_t^{\alpha }{p}_m^{\ast }(u,v)=\underset{0}{\overset{\infty }{\int }}[-k\lambda (u+v)p_m^n(u,v)+\lambda (u+v)\sum _{j=1}^{m\wedge k}{p}_{m-j}^n(u,v)]h_{\alpha }(t,u)du,\hspace{1em}m=1,2,\dots \mathrm{ }\mathrm{ }\end{array}$$(25) with the initial condition: $p_m^{\ast }(0,v)={\delta }_{m,0},$ where $p_m^n(u,v)$ is given by (22).

Proof.

Note that the mgf of $p_m^n(u,v)$ is of the form $${\widehat{p}}_s^n(u,v)=\mathbb{E}\left[s^{N^n(v+u)-N^n(v)}\right]=\hspace{0.17em}\text{exp}\hspace{0.17em}\{\Lambda (v,u+v)(s+\dots +s^k-k)\},$$ while the Laplace transform with respect to t of $h_{\alpha }(t,u)$ is given by (9). Taking both the mgf and the Laplace transform in (24) as above, we have (26) $${\overline{p}}_s^{\ast }(r,v)=\underset{0}{\overset{\infty }{\int }}{\widehat{p}}_s(u,v){\tilde{h}}_{\alpha }(r,u)du=r^{\alpha -1}\underset{0}{\overset{\infty }{\int }}\hspace{0.17em}\text{exp}\hspace{0.17em}\{\Lambda (v,u+v)(s+\dots +s^k-k)\}{e}^{-u{r}^{\alpha }}du.$$(26)

Note that for $U(u)=\hspace{0.17em}\text{exp}\hspace{0.17em}\{\Lambda (v,u+v)(s+\dots +s^k-k)\},$ we have (27) $$\frac{d}{du}U(u)=(s+s^2+\dots +s^k-k)\lambda (u+v)\hspace{0.17em}\text{exp}\hspace{0.17em}\{\Lambda (v,u+v)(s+s^2+\dots +s^k-k)\}.$$(27)

Thus, integrating (26) by parts with U as above, and $V=-e^{-u{r}^{\alpha }}/r^{\alpha },$ we get (28) $$\begin{array}{c}{\overline{p}}_s^{\ast }(r,v)=r^{\alpha -1}\{\left[-\frac{1}{r^{\alpha }}({\mathrm{e}}^{-u{r}^{\alpha }}{\mathrm{e}}^{\Lambda (v,u+v)(s+\dots s^k-k)}{|}_0^{\infty }\right]+\\ +\frac{1}{r^{\alpha }}(s+s^2+\dots +s^k-k)\underset{0}{\overset{\infty }{\int }}k\lambda (v,u+v)\hspace{0.17em}\text{exp}\hspace{0.17em}\{\Lambda (v,u+v)(s+s^2+\dots +s^k-k)\}{\mathrm{e}}^{-u{r}^{\alpha }}du\}=\\ =\frac{1}{r^{\alpha }}\left[r^{\alpha -1}+(s+s^2+\dots +s^k-k)\underset{0}{\overset{\infty }{\int }}\lambda (u+v)\hspace{0.17em}\text{exp}\hspace{0.17em}\{\Lambda (v,u+v)(s+s^2+\dots +s^k-k)\}{r}^{\alpha -1}{\mathrm{e}}^{-u{r}^{\alpha }}du\right].\end{array}$$(28)

We shall use the following property of the Caputo derivative: $$L_r\left\{D_t^{\alpha }f\right\}(r)=r^{\alpha }L\left\{f\right\}(r)-r^{\alpha -1}f(0^{+}),$$ where $L\{(f)\}(r)$ stands for the Laplace transform of function $f.$ Note that $p_y^{\ast }(0^{+},v)=1,$ since $Y_{\alpha }(0)=0$ a.s. Hence, by (28) $$\begin{array}{c}{r}^{\alpha }{\overline{p}}_s^{\ast }(r,v)-r^{\alpha -1}{\overline{p}}_s^{\ast }(0,v)=L_r\{D_t^{\alpha }{\overline{p}}_s^{\ast }(r,v)\}(r)=\\ =(s+s^2+\dots +s^k-k)\underset{0}{\overset{\infty }{\int }}\lambda (u+v)\hspace{0.17em}\text{exp}\hspace{0.17em}\{\Lambda (v,u+v)(s+s^2+\dots +s^k-k)\}{r}^{\alpha -1}{\mathrm{e}}^{-u{r}^{\alpha }}du.\end{array}$$

Inverting the Laplace transform yields $$\begin{array}{c}{D}_t^{\alpha }{\widehat{p}}_s^{\ast }(t,v)=(s+s^2+\dots +s^k-k)\underset{0}{\overset{\infty }{\int }}\lambda (u+v)\hspace{0.17em}\text{exp}\hspace{0.17em}\{\Lambda (v,u+v)(s+s^2+\dots +s^k-k)\}{h}_{\alpha }(t,u)du=\\ =\underset{0}{\overset{\infty }{\int }}\lambda (u+v)(s+s^2+\dots +s^k-k){\widehat{p}}_s(u,v)h_{\alpha }(t,u)du,\end{array}$$ where the mgf $${\widehat{p}}_s(u,v)=\sum _m{s}^m{p}_m(u,v).$$

Finally, by inverting the mgf $(s+s^2+\dots +s^k-k){\widehat{p}}_s(u,v),$ we obtain: $$D_t^{\alpha }{p}_m^{\ast }(u,v)=\underset{0}{\overset{\infty }{\int }}\lambda (u+v)\left[-k{p}_m(u,v)+\sum _{j=1}^{m\wedge k}{p}_{m-j}(u,v)\right]h_{\alpha }(t,u)du,$$ since the mgf of $$-k{p}_m(u,v)+\sum _{j=1}^{m\wedge k}{p}_{m-j}(u,v)$$ is equal to $$\sum _m{s}^m[-k{p}_m(u,v)+\sum _{j=1}^{m\wedge k}{p}_{m-j}(u,v)]=(s+s^2+\dots +s^k-k){\widehat{p}}_s(u,v).$$

3.2.2. Covariance structure

One can show that for NPPk $\mathbb{E}[N^n(t)]=\frac{k(k+1)}{2}\Lambda (t),$ and its covariance function is $$\text{Cov}[N^n(t),N^n(s)]=\frac{k(k+1)(2k+1)}{6}\Lambda (\min (s,t)).$$

Then the mean and covariance function of FNPPk are given by $$\mathbb{E}[N_{\alpha }^{\ast }(t)]=\frac{k(k+1)}{2}\mathbb{E}[\Lambda (Y_{\alpha }(t)]$$ $$\text{Cov}[N_{\alpha }^{\ast }(t),N_{\alpha }^{\ast }(s)]=\frac{k(k+1)(2k+1)}{6}\mathbb{E}[\Lambda (Y_{\alpha }(\min (s,t))+{\left(\frac{k(k+1)}{2}\right)}^2\text{Cov}[\Lambda (Y_{\alpha }(t)),\Lambda (N_{\alpha }(s))].$$

4. Pólya-Aeppli process of order k

The Pólya-Aeppli process of order k was defined and studied in the context of ruin problems in Chukova and Minkova (2015) and later by Kostadinova and Lazarova (2019). Related pure fractional birth processes were studied in Orsingher and Polito (2010).

Definition 6.

The process $N_{\mathit{PAk}}(t)$ is said to be the Pólya-Aeppli process of order k (PAk) if $$N_{\mathit{PAk}}(t)=X_1+\dots +X_{N_1(t)},$$ where (i) the random variables X_i are i.i.d with the truncated geometric distribution of parameter $\rho \in [0,1),$ given by (2); (ii) $N=\{N(t);t\ge 0\}$ is a homogeneous Poisson process (HPP) with intensity $\lambda >0,$ independent of ${\left\{X_i\right\}}_{i=1}^{\infty }.$

The following Kolmogorov forward equations are valid for the marginal distributions $p_m(t)=\mathbb{P}[N_{\mathit{PAk}}(t)=m]:$ (29) $$\begin{array}{c}\frac{d}{dt}{p}_0(t)=-\lambda p_0(t)\\ \frac{d}{dt}{p}_m(t)=-\lambda p_m(t)+\lambda \frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^{m\wedge k}{\rho }^{j-1}{p}_{m-j}(t),\end{array}$$(29) where $p_m(0)={\delta }_{m,0.}$

The marginal distributions of the PAk process are given by (30) $$p_m(t):=\mathbb{P}[N_{\mathit{PAk}}(t)=m]=q_m(\lambda t),m=0,1,2,..,$$(30) where q_m are given by (6).

More explicit expressions for $p_m(t)$ can be found in Minkova (2010). The expectation ad variance are as follows: (31) $$\begin{array}{c}\mathbb{E}[N_{\mathit{PAk}}(t)]=\lambda t\frac{1+\rho +\dots +{\rho }^{k-1}-k{\rho }^k}{1-{\rho }^k},\\ \text{Var}[N_{\mathit{PAk}}(t)]=\frac{\lambda t}{1-{\rho }^k}\left[1+3\rho +5{\rho }^2+\dots +(2k-1){\rho }^{k-1}-k^2{\rho }^k\right].\end{array}$$(31)

Note that PAk process is a compound Poisson process with the pgf $$G_{N_{\mathit{PAk}}(t)}(u)=\mathbb{E}\left[u^{N_{\mathit{PAk}}(t)}\right]=\mathbb{P}[N_{\mathit{PAk}}(t)=m]=e^{-\lambda t(1-G_X(u))},$$ where $G_X(u)=\mathbb{E}\left[u^X\right]$ is given by (3).

4.1. Non-homogeneous Pólya-Aeppli process of order k

We now consider a non-homogeneous version by introducing a deterministic time dependent intensity function $\lambda (t)$ as above, and $\Lambda (s,s+t)=\Lambda (s+t)-\Lambda (s),\Lambda (t)={\int }_0^t\lambda (u)du.$

Definition 7.

(Non-homogeneous Pólya-Aeppli process of order k). We define a non-homogeneous Pólya-Aeppli process of order k with cumulative rate function $\Lambda (t)$ and parameter ρ as (32) $$N_{\mathit{PAk}}^n(t)=X_1+\dots +X_{N_1^n(t)},$$(32) where (i) $\{N_1^n(t);t\ge 0\}$ is a non-homogeneous Poisson process (NPP) with cumulative rate function $\Lambda (t);$ (ii) X_i are i.i.d. r.v’s following the truncated geometric distribution with parameter $\rho ,$ given by (3); (iii) $\{N_1^n(t);t\ge 0\}$ is independent from $X_i,i=1,2,\dots$

Note, that the random variable $N_{\mathit{PAk}}^n(t+s)-N_{\mathit{PAk}}^n(s),:s,:t\ge 0$ has the Pólya-Aeppli distribution of order k with parameters $\Lambda (s,t),\rho ,$ that is (33) $$f_m^n(t,u)=\mathbb{P}[N_{\mathit{PAk}}^n(t+u)-N_{\mathit{PAk}}^n(u)=m]=q_m(\Lambda (u,u+t)),m=0,1,2,\dots ,$$(33) where q_m are given by (6).

Then the marginal distributions of the process $N_{\mathit{PAk}}^n(t)$ are $\mathbb{P}[N_{\mathit{PAk}}^n(t)=m]=f_m^n(t,0).$ An alternative definition can be given in terms of transition probabilities.

Definition 8.

The counting process $N_{\mathit{PAk}}^n(t)$ is said to be a non-homogeneous Pólya -Aeppli process of order k with the rate function $\lambda (t)$ and parameter $\rho \in [0,1)$ if (1) $N_{\mathit{PAk}}^n(0)=0;$ (2) $N_{\mathit{PAk}}^n(t)$ has independent increments; (3) for all $t\ge 0$ (34) $$\mathbb{P}[N_{\mathit{PAk}}^n(t+h)=n|N_{\mathit{PAk}}^n(t)=m]=\{\begin{array}{c}1-\lambda (t+h)h+o(h),\hspace{1em}n=m\hfill \\ \frac{1-\rho }{1-{\rho }^k}{\rho }^{i-1}\lambda (t+h)h+o(h),\hspace{1em}n=m+i,i=1,2,\dots ,k\hfill \end{array}$$(34)

It is easy to verify that the previous two definitions are equivalent.

4.1.1. Marginal distributions of the process

The following theorem holds.

Theorem 4.1.

The functions $f_m(t,u),m=0,1,2,\dots$ satisfy the differential equation: (35) $$\frac{d}{dt}{f}_m(t,u)=-\lambda (t+u)f_m(t,u)+\lambda (t+u)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^{m\wedge k}{\rho }^{j-1}{f}_{x-j}(t,u).$$(35)

Proof.

We first consider the case $m=0.$ By fixing u and taking a small h we can write $$\begin{array}{c}{f}_0(t+h,u)=P[I(t+h)=0]=\mathbb{P}[N_{\mathit{PAk}}^n(t+u+h)-N_{\mathit{PAk}}^n(u)=0]=\\ \mathbb{P}[\text{no}\mathrm{ }\text{events}\mathrm{ }\text{in}\mathrm{ }(u,u+t]\cap \mathrm{ }\text{no}\mathrm{ }\text{events}\mathrm{ }\text{in}\mathrm{ }(u+t,u+t+h]]=\\ \mathbb{P}[\text{no}\mathrm{ }\text{events}\mathrm{ }\text{in}\mathrm{ }(u,u+t+h]]\mathbb{P}[\text{no}\mathrm{ }\text{events}\mathrm{ }\text{in}\mathrm{ }(u+t,u+t+h]]=\\ f_0(t+h)[1-\lambda (t+u)h+o(h)]\end{array}$$

Thus $$\frac{f_0(t+h,u)-f_0(t,u)}{h}=-\lambda (t+u)f_0(t,u)+\frac{o(h)}{h}.$$

Letting $h\to 0$ yields $$\frac{d}{dt}{f}_0(t,u)=-\lambda (t+u)f_0(t,u).$$

For $m\ge 1$ we have $$\begin{array}{c}{f}_m(t+h,u)=\mathbb{P}[\{\mathrm{m}\mathrm{ }\text{events}\mathrm{ }\text{in}\mathrm{ }(u,u+t+h]\}\cap \{\text{no}\mathrm{ }\text{events}\mathrm{ }\text{in}\mathrm{ }(u+t,u+t+h]\}\\ \cup \{m-1\mathrm{ }\text{events}\mathrm{ }\text{in}\mathrm{ }(u,u+t]\}\cap \{1\mathrm{ }\text{event}\mathrm{ }\text{in}(u+t,u+t+h]\}\cup \dots \\ \cup \{\mathrm{ }0\mathrm{ }\text{events}\mathrm{ }\text{in}\mathrm{ }(u,u+t]\}\cap \{\mathrm{m}\mathrm{ }\text{events}\mathrm{ }\text{in}\mathrm{ }(u+t,u+t+h]\}]=\\ f_m(t+h,u)\left[1-\lambda (t+u)h+o(h)\right]+f_{m-1}(t+h,u)\left[\frac{1-\rho }{1-{\rho }^k}\lambda (t+u)h{\rho }^{1-1}+o(h)\right]+\dots \\ +f_0(t+h,u)\left[\frac{1-\rho }{1-{\rho }^k}\lambda (t+u)h{\rho }^{m\wedge k-1}+o(h)\right]=\\ =\lambda (t+u)f_m(t+h,u)+\lambda (t+u)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^{m\wedge k}{\rho }^{j-1}{f}_{m-j}(t,u).\end{array}$$

Letting $h\to 0$ yields $$\frac{d}{dt}{f}_m(t,u)=-\lambda (t+u)f_m(t,u)+\lambda (t+u)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^{m\wedge k}{\rho }^{j-1}{f}_{m-j}(t,u),$$ which was the statement of the theorem. □

Note that in case $k\to \infty ,$ the Pólya-Aeppli process $N_{\mathit{PAk}}^n(t)$ coincides with the non-homogeneous Pólya-Aeppli process defined in Chukova and Minkova (2019), but for fixed k the Pólya-Aeppli process $N_{\mathit{PAk}}^n(t)$ is new.

4.2. Fractional Pólya-Aeppli process of order k

To the best of our knowledge, fractional versions of PAk processes have not been considered yet. We define a fractional Pólya-Aeppli process of order $k$ as a Pólya-Aeppli process of order k time-changed by the process $\{Y_{\alpha }(t);t\ge 0\},$ such that (36) $$N_{\alpha }^h(t)=N_{\mathit{PAk}}(Y_{\alpha }(t)),\hspace{1em}0<\alpha <1,$$(36) where (i) $N_1=\{N_1(t);t\ge 0\}$ is the homogeneous Poisson process with intensity $\lambda ;$ (ii) $N_{\mathit{PAk}}(t)=X_1+\dots +X_{N_1(t)};$ (iii) $\{Y_{\alpha }(t);t\ge 0\},0<\alpha <1$ is the inverse α-stable subordinator, defined in (7) and independent of $N_1(t).$

4.2.1. Marginal distributions

We shall now obtain governing equations for the marginal distributions of the fractional PAk process $$p_x^{\alpha }(t)=\mathbb{P}[N_{\mathit{PAk}}(Y_{\alpha }(t))=m]=\underset{0}{\overset{\infty }{\int }}{p}_m(u)h_{\alpha }(t,u)du,\hspace{1em}m=0,1,\dots ,$$ where $p_m(u)$ is given by (30).

Theorem 4.2.

The probabilities $p_x^{\alpha }(t),\hspace{1em}x=0,1,\dots$ satisfy the fractional differential-difference equations: (37) $$D_t^{\alpha }{p}_0^{\alpha }(t)=-\lambda p_0^{\alpha }(t)$$(37) (38) $$D_t^{\alpha }{p}_x^{\alpha }(t)=-\lambda p_x^{\alpha }(t)+\lambda \frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^{x\wedge k}{\rho }^{j-1}{p}_{x-j}^{\alpha }(t),$$(38) where $D_t^{\alpha }f(t)$ is the fractional Caputo derivative of the function f given by (17).

Proof.

We first consider the case $m\ge 1.$ By taking the fractional Caputo derivative of both sides in (29) and using the property (20), we get $$\begin{array}{c}{D}_t^{\alpha }{p}_m^{\alpha }(t)=-\underset{0}{\overset{\infty }{\int }}{p}_m(u)\frac{\partial }{\partial u}{h}_{\alpha }(t,u)du=\\ =\underset{0}{\overset{\infty }{\int }}\left[-\lambda p_m(u)+\lambda \frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^{x\wedge k}{\rho }^{j-1}{p}_{m-j}(t)\right]h_{\alpha }(t,u)du-p_m(u)h_{\alpha }(t,u){|}_0^{\infty }=\\ =-\lambda p_m^{\alpha }(t)+\lambda \frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^{x\wedge k}{\rho }^{j-1}{p}_{m-j}^{\alpha }(t).\end{array}$$

For m = 0 we have $$\begin{array}{c}{D}_t^{\alpha }{p}_0^{\alpha }(t)=-\underset{0}{\overset{\infty }{\int }}{p}_0(u)\frac{\partial }{\partial u}{h}_{\alpha }(t,u)du=\\ =\underset{0}{\overset{\infty }{\int }}[-\lambda p_0(u)]h_{\alpha }(t,u)du=-\lambda p_0^{\alpha }(t).\end{array}$$

4.3. Correlation structure and long-range dependence property

In this sub-section we shall obtain several important characteristics of the fractional Pólya-Aeppli process of order k such as its expectation, variance and covariance. After that, we are able to study the correlation structure of the process. For the fractional Pólya-Aeppli process of order k, $N_{\alpha }^h(t)=N_{\mathit{PAk}}(Y_{\alpha }(t)),$ we can use the property of the conditional expectation to write (see (Leonenko et al. 2014)) $$\begin{array}{c}\mathbb{E}[N_{\alpha }^h(t)]=\mathbb{E}\left[\mathbb{E}\right[N_{\alpha }^h(t)|Y_{\alpha }(t)]|Y_{\alpha }(t)]=\underset{0}{\overset{\infty }{\int }}\mathbb{E}[N_{\mathit{PAk}}(u)]h_{\alpha }(t,u)du=\\ =\lambda \mathbb{E}[N_{\mathit{PAk}}(1)]\frac{t^{\alpha }}{\Gamma (\alpha +1)},\\ \text{Var}[N_{\alpha }^h(t)]=\frac{t^{\alpha }\text{Var}[N_{\mathit{PAk}}(1)]}{\Gamma (\alpha +1)}+\frac{t^{2\alpha }{(\mathbb{E}[N_{\mathit{PAk}}(1)])}^2}{\alpha }\left(\frac{1}{\Gamma (2\alpha )}-\frac{1}{\alpha \Gamma {(\alpha )}^2}\right).\end{array}$$

The covariance function can be calculated via the formula: $$\text{Cov}[N_{\alpha }^h(t),N_{\alpha }^h(s)]=\text{Var}[N_{\mathit{PAk}}(1)]\frac{\min {(t,s)}^{\alpha }}{\Gamma (1+\alpha )}+{(\mathbb{E}[N_{\mathit{PAk}}(1)])}^2\text{Cov}[Y_{\alpha }(t),Y_{\alpha }(s)],$$ where the covariance of the process $Y_{\alpha }(t)$ is given by Equation (11)(11) $$\text{Cov}[Y_{\alpha }(t),Y_{\alpha }(s)]=\frac{1}{\Gamma (1+\alpha )\Gamma (\alpha )}\underset{0}{\overset{\min (t,s)}{\int }}((t-\tau {)}^{\alpha }+{(s-\tau )}^{\alpha }){\tau }^{\alpha -1}d\tau -\frac{{(st)}^{\alpha }}{{\Gamma }^2(1+\alpha )}.$$(11) .

Theorem 4.3.

The process $N_{\alpha }^h(t)$ has the LRD property.

Proof.

Using the results from (Leonenko et al. 2014) similarly to the previous section, we get $$\text{Corr}[N_{\alpha }^h(t),N_{\alpha }^h(s)]\sim t^{-\alpha }C(\alpha ,s)\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}t\to \infty ,$$ where $C(\alpha ,s)={(\frac{1}{\Gamma (2\alpha )}-\frac{1}{\alpha {(\Gamma (\alpha ))}^2})}^{-1}\left[\frac{\alpha \text{Var}[N_{\mathit{PAk}}(1)]}{\Gamma (1+\alpha ){(\mathbb{E}[N_{\mathit{PAk}}(1)])}^2}+\frac{\alpha s^{\alpha }}{\Gamma (1+2\alpha )}\right],$ and $\mathbb{E}[N_{\mathit{PAk}}(1)],\text{Var}[N_{\mathit{PAk}}(1)]$ are given by (31). Thus the correlation function of FPAk process decays at rate $t^{-\alpha },\alpha \in (0,1)$ and satisfies the LRD property. □

4.4. Non-homogeneous fractional PAk process

As we did before, we can now define a non-homogeneous fractional Pólya-Aeppli process of order k as $$N_{\alpha }^n(t)=N_{\mathit{PAk}}(\Lambda (Y_{\alpha }(t))),\hspace{1em}t\ge 0,\hspace{1em}0<\alpha <1,$$ where all the symbols have the usual meaning defined above. We assume that the inverse subordinator $Y_{\alpha }$ is independent of the process $N_{\mathit{PAk}}.$ In this sub-section, we shall derive governing equations for the probabilities $$p_m^{\ast \ast }(t,v)=\mathbb{P}[N_{\mathit{PAk}}(\Lambda (Y_{\alpha }(t)+v))-N_{\mathit{PAk}}(\Lambda (v))=m].$$

Theorem 4.4.

The marginal distributions $p_x^{\ast \ast }(t,v)$ satisfy the following fractional differential-difference integral equations (39) $$D_t^{\alpha }{p}_0^{\ast \ast }(u,v)=-\underset{0}{\overset{\infty }{\int }}\lambda (u+v)f_0^n(u,v)h_{\alpha }(t,u)du$$(39) (40) $$D_t^{\alpha }{p}_m^{\ast \ast }(u,v)=\underset{0}{\overset{\infty }{\int }}\lambda (u+v)\left[-f_m^n(u,v))+\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^{m\wedge k}{\rho }^{j-1}{f}_{m-j}^n(u,v)\right]h_{\alpha }(t,u)\mathit{\text{du\hspace{1em}m}}=1,2,\dots$$(40) with the initial condition $p_m^{\ast \ast }(0,v)={\delta }_{m,0},$ where $f_m^n(u,v)$ is given by (33).

Proof.

Using (9), the mgf of $f_m^n(u,v)$ can be written in the form: $${\widehat{f}}_s^n(u,v)=\mathbb{E}\left[s^{N^n(v+u)-N^n(v)}\right]=\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{\Lambda (v,u+v)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^k{\rho }^{j-1}(s^j-1)\right\},$$ while the Laplace transform with respect to t of $h_{\alpha }(t,u)$ is given by (9). Taking both the mgf and the Laplace transform in (24) as above, we have (41) $$\begin{array}{c}{\overline{p}}_s^{\ast \ast }(u,v)=r^{\alpha -1}\underset{0}{\overset{\infty }{\int }}{\widehat{f}}_s^n(u,v){\tilde{h}}_{\alpha }(r,u)du=\\ \underset{0}{\overset{\infty }{\int }}\left[\text{exp}\hspace{0.17em}\left\{\Lambda (v,u+v)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^k{\rho }^{j-1}(s^j-1)\right\}\right]e^{-u{r}^{\alpha }}du.\end{array}$$(41)

Note that for $U(u)=\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{\Lambda (v,u+v)\frac{1-\rho }{1-{\rho }^k}{\sum }_{j=1}^k{\rho }^{j-1}(s^j-1)\right\}$ one can take derivative in u as follows: (42) $$\frac{d}{du}U(u)=\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^k{\rho }^{j-1}(s^j-1)[\lambda (v,u+v)]\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{\Lambda (v,u+v)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^k{\rho }^{j-1}(s^j-1)\right\}.$$(42)

Thus, integrating (41) by parts with $$U=\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{\Lambda (v,u+v)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^k{\rho }^{j-1}(s^j-1)\right\},\hspace{1em}V=-\frac{1}{r^{\alpha }}{e}^{-u{r}^{\alpha }},$$ we get (43) $$\begin{array}{c}{\overline{p}}_s^{\ast \ast }(u,v)=\frac{1}{r^{\alpha }}[r^{\alpha -1}+\frac{1-\rho }{1-{\rho }^k}[(s-1)+\rho (s^2-1)+\dots +{\rho }^{k-1}(s^k-1)]\times \\ \times \underset{0}{\overset{\infty }{\int }}\lambda (v,u+v)\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{\Lambda (v,u+v)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^k{\rho }^{j-1}(s^j-1)\right\}{\mathrm{e}}^{-u{r}^{\alpha }}{r}^{\alpha -1}du\end{array}$$(43) where $p_s^{\ast \ast }(0^{+},v)=1,$ since $Y_{\alpha }(0)=0$ a.s. Hence, by (43) $$\begin{array}{c}{r}^{\alpha }{\overline{p}}_s^{\ast \ast }(r,v)-r^{\alpha -1}{\overline{p}}_s^{\ast \ast }(0,v)=L_r\{D_t^{\alpha }{\overline{p}}_s^{\ast \ast }(r,v)\}(r)=\\ =\frac{1-\rho }{1-{\rho }^k}[(s-1)+\rho (s^2-1)+\dots +{\rho }^{k-1}(s^k-1)]\times \\ \times \underset{0}{\overset{\infty }{\int }}\lambda (u+v)\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{\Lambda (v,u+v)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^k{\rho }^{j-1}(s^j-1)\right\}{r}^{\alpha -1}{\mathrm{e}}^{-u{r}^{\alpha }}du.\end{array}$$

Inverting the Laplace transform yields $$\begin{array}{c}{D}_t^{\alpha }{\widehat{p}}_s^{\ast \ast }(t,v)=\frac{1-\rho }{1-{\rho }^k}[(s-1)+\rho (s^2-1)+\dots +{\rho }^{k-1}(s^k-1)]\times \\ \times \underset{0}{\overset{\infty }{\int }}\lambda (u+v)\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{\Lambda (v,u+v)\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^k{\rho }^{j-1}(s^j-1)\right\}{h}_{\alpha }(t,u)du=\\ =\underset{0}{\overset{\infty }{\int }}\lambda (u+v)\left[\frac{1-\rho }{1-{\rho }^k}[(s-1)+\rho (s^2-1)+\dots +{\rho }^{k-1}(s^k-1)\right]{\widehat{f}}_s^n(u,v)h_{\alpha }(t,u)du,\end{array}$$ where the mgf is $${\widehat{f}}_s^n(u,v)=\sum _m{s}^m{f}_m^n(u,v).$$

Finally, by inverting the mgf $$[\frac{1-\rho }{1-{\rho }^k}[(s-1)+\rho (s^2-1)+\dots +{\rho }^{k-1}(s^k-1)]{\widehat{f}}_s^n(u,v),$$ we obtain: $$D_t^{\alpha }{p}_m^{\ast }(u,v)=\underset{0}{\overset{\infty }{\int }}\lambda (u+v)\left[-f_m^n(u,v)+\frac{1-\rho }{1-{\rho }^k}\sum _{j=1}^{m\wedge k}{\rho }^{j-1}{f}_{m-j}^n(u,v)\right]h_{\alpha }(t,u)du.$$ □

5. Final notes

The counting processes of order k that we have discussed in this paper have this general form (44) $$N(t)=\sum _{i=1}^{\mathcal{N}(t)}{X}_i,$$(44) where ${\left\{X_i\right\}}_{i=1}^{\infty }$ is a sequence of i.i.d. integer random variables assuming values in $1,\dots ,k$ and $\mathcal{N}(t)$ is a counting process independent from the sequence. One further assumes that $N(0)=0.$ A simple algorithm in R is given in arXiv:2008.09421 [math.PR] when $\mathcal{N}(t)$ is the fractional Poisson process of renewal type used above and discussed by Mainardi, Gorenflo, and Scalas (2004) and when X₁ is uniformly distributed in $1,\dots ,k.$

Acknowledgments

The authors thank Mr. Mostafizar Khandakar and Dr. Kuldeep Kumar Kataria for their thorough reading of the manuscript that highlighted some typos and a wrong formula. Enrico Scalas acknowledges partial support from the Dr Perry James (Jim) Browne Research Center at the Department of Mathematics, University of Sussex.