Multiparameter Skellam Process

• Multiparameter Skellam Process is an extension of the classical Skellam process, capturing signed count dynamics with multidimensional parameterizations including diverse jump sizes and rate vectors.
• It employs various constructions such as order-k generalizations, set-indexed random fields, and time-change mechanisms (fractional and tempered) to address complex stochastic behaviors.
• Its tractable structure features explicit generating functions, closed-form state probabilities, and governing equations, offering robust tools for modeling asymmetric and heavy-tailed count data.

Searching arXiv for the cited Skellam-process papers to ground the synthesis in the latest relevant literature. A multiparameter Skellam process is a family of integer-valued stochastic models obtained by extending the classical Skellam process—the difference of two independent Poisson processes—along more than one structural axis. In the literature, “multiparameter” appears in several distinct but related senses: an order parameter controlling admissible jump sizes, multiple rate vectors governing positive and negative components, random-field indexing over $\mathbb{R}^M_+$, fractional and tempered time changes through subordinators and inverse subordinators, and interacting vector constructions with migration-type dependence (Gupta et al., 2020, Buchak et al., 2017, Vishwakarma, 13 Sep 2025, Tathe et al., 11 Apr 2025, Tathe et al., 2024, Vishwakarma, 16 Sep 2025). Across these variants, the common core is a signed count process or field built from independent Poisson-type drivers, with explicit probability generating functions, characteristic functions, state probabilities, governing difference–differential equations, and Lévy or subordination representations.

1. Classical Skellam structure and the meaning of multiparameter

The classical Skellam process is defined by

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$

where $N_1$ and $N_2$ are independent homogeneous Poisson processes with intensities $\lambda_1>0$ and $\lambda_2>0$ (Gupta et al., 2020). Its fixed-time law is the Skellam distribution,

$$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$

with modified Bessel function of the first kind (Gupta et al., 2020). The mean and variance are

$$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$

and the process is a Lévy process with Lévy measure $$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$ (Buchak et al., 2017).

A first notion of multiparameter generalization replaces the two scalar rates by richer parameter sets. In the generalized Skellam process,

$S(t)=M_1(t)-M_2(t),$

the components $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$0 are generalized counting processes with jump amplitudes $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$1 and rate vectors $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$2 and $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$3 (Kataria et al., 2021). The corresponding pgf is

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$4

and the first two cumulants are determined by

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$5

so that $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$6 and $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$7 (Kataria et al., 2021).

A second notion of multiparameter is set-indexed or field-indexed. For $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$8, a Skellam random field on $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$9 is defined by

$N_1$0

where $N_1$1 are independent Poisson random fields and $N_1$2 is a Borel set (Vishwakarma, 13 Sep 2025). For $N_1$3, the two-parameter Skellam sheet $N_1$4 has rectangular increments that are independent and stationary, and

$N_1$5

(Vishwakarma, 13 Sep 2025).

A third notion of multiparameter arises from time change. The process may be subordinated by a general subordinator $N_1$6, by stable or tempered stable subordinators, by inverse stable subordinators, or by compound Poisson–Gamma clocks, thereby adding parameters such as $N_1$7, $N_1$8, or $N_1$9 to the underlying Skellam structure (Gupta et al., 2020, Buchak et al., 2017, Tathe et al., 11 Apr 2025). This suggests that “multiparameter Skellam process” is not a single object but a class of signed Poisson-type models with several simultaneous control parameters.

2. Order-$N_2$0 and generalized jump-size constructions

One major extension replaces the Poisson components by Poisson processes of order $N_2$1. A Poisson process of order $N_2$2, denoted $N_2$3, allows jumps of sizes $N_2$4 per arrival epoch, has Lévy measure

$N_2$5

and satisfies

$N_2$6

(Gupta et al., 2020). Its mean and variance are

$N_2$7

(Gupta et al., 2020).

The Skellam process of order $N_2$8 is then

$N_2$9

with independent $\lambda_1>0$0 of intensities $\lambda_1>0$1 (Gupta et al., 2020). Its pgf and characteristic function are

$\lambda_1>0$2

$\lambda_1>0$3

The Lévy measure is supported on $\lambda_1>0$4:

$\lambda_1>0$5

so the process has jumps of sizes $\lambda_1>0$6 at rates $\lambda_1>0$7 and $\lambda_1>0$8 at rates $\lambda_1>0$9 (Gupta et al., 2020).

The paper also gives the marginal probabilities

$\lambda_2>0$0

together with governing equations

$\lambda_2>0$1

(Gupta et al., 2020). For $\lambda_2>0$2 the model reduces to the classical Skellam process, while $\lambda_2>0$3 gives symmetry (Gupta et al., 2020).

A related but more general jump-size framework appears in the generalized Skellam process of Kataria and Khandakar, where the positive and negative legs are generalized counting processes with arbitrary rate vectors $\lambda_2>0$4 and $\lambda_2>0$5 on $\lambda_2>0$6 (Kataria et al., 2021). The non-homogeneous extension replaces these by deterministic time-dependent intensities $\lambda_2>0$7 and $\lambda_2>0$8 with cumulative rates $\lambda_2>0$9 and $$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$0, and defines

$$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$1

(Tathe et al., 2024). The resulting pgf is

$$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$2

and the mean and variance are

$$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$3

(Tathe et al., 2024).

An even broader formulation is the Poisson–Skellam family

$$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$4

with a countable jump set $$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$5 and independent Poisson processes $$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$6 (Cinque et al., 10 Apr 2025). In the homogeneous case it is a Lévy process with exponent

$$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$7

and Lévy measure

$$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$8

(Cinque et al., 10 Apr 2025). This formulation includes classical, order-$$s_{k}(t)=e^{-t(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{k/2}I_{|k|}(2t\sqrt{\lambda_1 \lambda_2}),\qquad k\in \mathbb{Z},$$9, and signed fixed-jump-size Skellam-type processes as special cases.

3. Fractional, tempered, and time-changed Skellam dynamics

A central strand of the literature builds multiparameter Skellam models by randomizing the operational time. For a general subordinator $$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$0 with Laplace transform

$$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$1

the time-changed order-$$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$2 process

$$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$3

has moment generating function

$$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$4

(Gupta et al., 2020). The stable choice $$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$5 yields a space-fractional model, and the tempered stable choice $$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$6 yields a tempered space-fractional model (Gupta et al., 2020).

The space-fractional Skellam process is defined by

$$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$7

with independent stable subordinators $$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$8 (Gupta et al., 2020). Its mgf is

$$\mathbb{E}[S(t)]=(\lambda_1-\lambda_2)t,\qquad {\rm Var}[S(t)]=(\lambda_1+\lambda_2)t,$$9

and the governing equations are fractional difference–differential equations involving $$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$0 and $$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$1 (Gupta et al., 2020). Its Lévy measure has infinite support,

$$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$2

so fractionalization replaces finitely many jump sizes by infinitely many combinatorially weighted sizes (Gupta et al., 2020).

The tempered space-fractional Skellam process,

$$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$3

has mgf

$$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$4

(Gupta et al., 2020). As $$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$5, the stable case is recovered (Gupta et al., 2020).

Another generalized space-time fractional construction is the generalized space-time fractional Skellam process

$$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$6

where $$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$7 is a generalized Skellam process, $$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$8 is a stable subordinator, and $$\nu(dy)=\lambda_1\delta_1(dy)+\lambda_2\delta_{-1}(dy)$$9 is an inverse stable subordinator (Tathe et al., 11 Apr 2025). Its p.g.f. is

$S(t)=M_1(t)-M_2(t),$0

and the closed-form p.m.f. is expressed through derivatives of Mittag–Leffler functions with respect to $S(t)=M_1(t)-M_2(t),$1 and $S(t)=M_1(t)-M_2(t),$2 (Tathe et al., 11 Apr 2025). The paper states that for $S(t)=M_1(t)-M_2(t),$3, $S(t)=M_1(t)-M_2(t),$4 has infinite moments, so mean and variance may be infinite in space-fractional cases (Tathe et al., 11 Apr 2025).

Time change by inverse stable subordinators produces non-Lévy fractional Skellam models. The generalized fractional Skellam process

$S(t)=M_1(t)-M_2(t),$5

with inverse $S(t)=M_1(t)-M_2(t),$6-stable subordinator $S(t)=M_1(t)-M_2(t),$7, has pgf

$S(t)=M_1(t)-M_2(t),$8

(Kataria et al., 2021). Its one-dimensional distributions are not infinitely divisible (Kataria et al., 2021). The non-homogeneous generalized fractional Skellam process

$S(t)=M_1(t)-M_2(t),$9

is defined analogously, with integral representation

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$00

where $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$01 is the density of the inverse stable subordinator (Tathe et al., 2024).

Compound Poisson–Gamma time changes provide another multiparameter family. If $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$02 is a compound Poisson–Gamma subordinator with parameters $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$03, then

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$04

has characteristic function

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$05

mean

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$06

and variance

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$07

(Buchak et al., 2017). In contrast, inverse-time changes $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$08 are not Lévy and have non-stationary dependent increments (Buchak et al., 2017).

4. Random fields and genuinely multiparameter indexing

A distinct branch of the theory studies Skellam random fields indexed by multidimensional positive orthants. For independent Poisson random fields $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$09 on $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$10, the Skellam field is

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$11

with

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$12

(Vishwakarma, 13 Sep 2025). For $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$13, the Skellam sheet $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$14 has covariance

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$15

(Vishwakarma, 13 Sep 2025).

The point probabilities and generating functions retain the classical Skellam form with area or volume replacing time. For rectangles in $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$16,

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$17

and

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$18

The forward equations are

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$19

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$20

(Vishwakarma, 13 Sep 2025). The generator acting on $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$21 is

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$22

(Vishwakarma, 13 Sep 2025).

The paper "Skellam Processes via Multiparameter Poisson Process" adopts a different, additive multiparameter Poisson framework. A multiparameter Poisson process with linear increments has one-dimensional marginals

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$23

and representation

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$24

where the $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$25 are independent one-parameter Poisson processes (Vishwakarma, 16 Sep 2025). The multiparameter Skellam process is then

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$26

with independent multiparameter Poisson processes $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$27 of rate vectors $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$28 (Vishwakarma, 16 Sep 2025). Its marginal law is

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$29

with

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$30

(Vishwakarma, 16 Sep 2025).

The generalized multiparameter Skellam process in the same framework is

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$31

for finite $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$32 and independent multiparameter Poisson processes $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$33 (Vishwakarma, 16 Sep 2025). Its pgf is

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$34

which is the field-indexed analogue of the one-parameter generalized Skellam pgf (Vishwakarma, 16 Sep 2025).

A related field perspective also appears in the additive Lévy framework of the generalized Poisson–Skellam family,

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$35

with multiparameter Poisson drivers and fixed jumps $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$36 (Cinque et al., 10 Apr 2025). The resulting random field has independent increments over disjoint rectangles and characteristic function

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$37

(Cinque et al., 10 Apr 2025).

5. Governing equations, transforms, and structural characterizations

One reason Skellam-type models remain tractable is the availability of explicit transform and operator formulas. In the classical one-parameter case, the forward equation is

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$38

with initial condition $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$39, $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$40 for $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$41 (Gupta et al., 2020). For the generalized Skellam process built from generalized counting processes, the pmf satisfies

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$42

where $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$43 and $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$44 (Kataria et al., 2021).

In the order-$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$45 setting, the operator structure becomes

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$46

where $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$47 and $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$48 are backward and forward shift operators (Gupta et al., 2020). In the space-fractional setting, fractional powers of these difference operators appear:

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$49

(Gupta et al., 2020). Tempering replaces these by $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$50 and $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$51 (Gupta et al., 2020).

Inverse-stable time changes replace ordinary derivatives by Caputo derivatives. For the generalized fractional Skellam process,

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$52

the state probabilities solve a fractional equation with Caputo derivative, and the pgf is a Mittag–Leffler eigenfunction of the generator (Kataria et al., 2021). In the non-homogeneous generalized fractional Skellam process,

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$53

where $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$54 is the inverse stable density (Tathe et al., 2024). The generalized space-time fractional Skellam process has a p.g.f. governing equation involving a Caputo derivative and series coefficients $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$55 built from $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$56 and $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$57 (Tathe et al., 11 Apr 2025).

Random-field variants have coupled partial difference–differential equations. For the Skellam sheet,

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$58

(Vishwakarma, 13 Sep 2025). For the fractional Skellam random field of Type I,

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$59

(Vishwakarma, 13 Sep 2025).

A general operator-level formulation is provided by Bernstein-function subordination. If $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$60 is the Skellam generator and $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$61 is the Bernstein function of a compound Poisson–Gamma clock, then the subordinated semigroup has generator $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$62 (Buchak et al., 2017). More generally, the homogeneous Bernstein-fractional Poisson–Skellam family has transform

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$63

while inverse-Bernstein time changes lead to convolution-type derivatives $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$64 in the governing equations (Cinque et al., 10 Apr 2025, Cinque et al., 14 Oct 2025).

6. Running averages, interacting vectors, asymptotics, and modeling implications

Beyond the base processes, several papers study integral and running-average functionals. For the order-$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$65 Skellam process,

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$66

has characteristic function

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$67

and admits a compound Poisson representation with a mixed double uniform jump law (Gupta et al., 2020). Its mean and variance are

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$68

(Gupta et al., 2020).

For the generalized space fractional Skellam process, the running average

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$69

also admits a compound Poisson representation,

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$70

with $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$71 an SFPP and explicitly given characteristic function of $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$72 (Tathe et al., 11 Apr 2025). In the generalized Skellam process with constant intensities, the running average

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$73

satisfies

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$74

and has compound Poisson representation with mixed double uniform jumps (Tathe et al., 2024).

Interacting-vector versions extend Skellam processes to dependent multicomponent point processes. In the two-group model of "Interacting point processes", the vector $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$75 on $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$76 has a migration-type infinitesimal structure and admits the decomposition

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$77

where $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$78 are independent non-homogeneous Skellam processes with explicitly stated rate pairs (Cinque et al., 14 Oct 2025). The covariance is

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$79

for $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$80 (Cinque et al., 14 Oct 2025). In the homogeneous case, the vector has a compound Poisson representation with jump types $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$81 and corresponding probabilities (Cinque et al., 14 Oct 2025).

Asymptotic and dependence properties depend strongly on the time-change mechanism. The generalized fractional Skellam process $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$82 has long-range dependence, and its one-dimensional distributions are not infinitely divisible (Kataria et al., 2021). The non-homogeneous generalized fractional Skellam process inherits long and short range dependence patterns through the covariance structure of the inverse stable time change (Tathe et al., 2024). For the generalized space-time fractional Skellam process,

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$83

as $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$84, while the one-dimensional marginals are not infinitely divisible; by contrast, the generalized space fractional Skellam process with $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$85 is infinitely divisible (Tathe et al., 11 Apr 2025).

The random-field literature provides weak convergence from discrete triangular arrays to Skellam fields. In the Skellam random field setting, if the array probabilities converge to the required intensity limits and maximal probabilities vanish, then the finite-dimensional distributions converge to those of the generalized Skellam random field, with the classical Skellam field obtained when the jump set is $S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$86 (Vishwakarma, 13 Sep 2025). A parallel triangular-array weak convergence result is proved for the alternative multiparameter construction

$S(t)=N_1(t)-N_2(t), \qquad t\ge 0,$87

in the additive multiparameter Poisson framework (Vishwakarma, 16 Sep 2025).

These constructions collectively show that multiparameter Skellam processes form a tractable but heterogeneous class. Order parameters control admissible jump sizes, rate vectors determine asymmetry and dispersion, fractional indices introduce heavy-tailed waiting times and non-Markovian effects, tempering parameters modify the tail behavior of the Lévy measure, and multidimensional indexing replaces one-dimensional time by area, volume, or coordinatewise additive evolution (Gupta et al., 2020, Vishwakarma, 13 Sep 2025, Tathe et al., 11 Apr 2025). A plausible implication is that model selection within this class is driven less by a single canonical definition than by which notion of “multiparameter” is required: jump-order heterogeneity, time-change heterogeneity, set-indexed Lévy structure, or multicomponent interaction.