---
title: Multiparameter Skellam Process
url: https://www.emergentmind.com/topics/multiparameter-skellam-process
type: topic
---

# Multiparameter Skellam Process

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 [2003.09471], [1707.00523], [2509.10870], [2504.08374], [2407.19227], [2509.12729]. 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$ [2003.09471]. 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 [2003.09471]. 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)$ [1707.00523].

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 $M_1,M_2$ are generalized counting processes with jump amplitudes $1,2,\ldots,k$ and rate vectors $\{\lambda_j\}_{j=1}^k$ and $\{\mu_j\}_{j=1}^k$ [2107.08307]. The corresponding pgf is
$$
G_S(u,t)=\exp\left\{-t\sum_{j=1}^k\left[\lambda_j(1-u^j)+\mu_j(1-u^{-j})\right]\right\},
$$
and the first two cumulants are determined by
$$
m_1=\sum_{j=1}^k j(\lambda_j-\mu_j),\qquad
m_2=\sum_{j=1}^k j^2(\lambda_j+\mu_j),
$$
so that $\mathbb{E}[S(t)]=m_1t$ and ${\rm Var}[S(t)]=m_2t$ [2107.08307].

A second notion of multiparameter is set-indexed or field-indexed. For $M\ge 1$, a Skellam random field on $\mathbb{R}_+^M$ is defined by
$$
S(B)=N_1(B)-N_2(B),
$$
where $N_1,N_2$ are independent Poisson random fields and $B$ is a Borel set [2509.10870]. For $M=2$, the two-parameter Skellam sheet $X(s,t)=N_1(s,t)-N_2(s,t)$ has rectangular increments that are independent and stationary, and
$$
P\{X(s,t)=k\}
=e^{-(\lambda_1+\lambda_2)st}\left(\frac{\lambda_1}{\lambda_2}\right)^{k/2}
I_{|k|}(2st\sqrt{\lambda_1\lambda_2})
$$
[2509.10870].

A third notion of multiparameter arises from time change. The process may be subordinated by a general subordinator $D_f$, by stable or tempered stable subordinators, by inverse stable subordinators, or by compound Poisson–Gamma clocks, thereby adding parameters such as $(\alpha_1,\alpha_2)$, $(\mu_1,\mu_2)$, or $(\nu,\alpha,\beta)$ to the underlying Skellam structure [2003.09471], [1707.00523], [2504.08374]. 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-$k$ and generalized jump-size constructions

One major extension replaces the Poisson components by Poisson processes of order $k$. A Poisson process of order $k$, denoted $N^k(t)$, allows jumps of sizes $1,2,\ldots,k$ per arrival epoch, has Lévy measure
$$
\nu_{N^k}(x)=\lambda\sum_{j=1}^k\delta_j(x),
$$
and satisfies
$$
\phi_{N^{k}(t)}(u)=e^{-\lambda t\left(k-\sum_{j=1}^{k}e^{iuj}\right)},\qquad
G^{N^{k}}(s,t)=e^{-\lambda t\left(k-\sum_{j=1}^{k} s^{j}\right)}
$$
[2003.09471]. Its mean and variance are
$$
\mathbb{E}[N^{k}(t)] = \frac{k(k+1)}{2}\lambda t,\qquad
{\rm Var}[N^{k}(t)]= \frac{k(k+1)(2k+1)}{6}\lambda t
$$
[2003.09471].

The Skellam process of order $k$ is then
$$
S^{k}(t)=N_1^k(t)-N_2^k(t),
$$
with independent $N_1^k,N_2^k$ of intensities $\lambda_1,\lambda_2$ [2003.09471]. Its pgf and characteristic function are
$$
G^{S^{k}}(s,t)=e^{-t\left(k(\lambda_1 + \lambda_2) -\lambda_1\sum_{j=1}^{k}s^{j} -\lambda_2 \sum_{j=1}^{k}s^{-j}\right)},
$$
$$
\phi_{S^{k}(t)}(u) = e^{-t\left[k(\lambda_1 + \lambda_2) -\lambda_1\sum_{j=1}^{k}e^{iju} -\lambda_2 \sum_{j=1}^{k}e^{-iju}\right]}.
$$
The Lévy measure is supported on $\{\pm1,\ldots,\pm k\}$:
$$
\nu_{S^k}(x)=\lambda_1\sum_{j=1}^k\delta_j(x)+\lambda_2\sum_{j=1}^k\delta_{-j}(x),
$$
so the process has jumps of sizes $+1,\ldots,+k$ at rates $\lambda_1$ and $-1,\ldots,-k$ at rates $\lambda_2$ [2003.09471].

The paper also gives the marginal probabilities
$$
R_{m}(t)=e^{-kt(\lambda_1+\lambda_2)}{\left(\frac{\lambda_1}{\lambda_2}\right)}^{m/2}I_{|m|}(2tk\sqrt{\lambda_1 \lambda_2}),\qquad m\in\mathbb{Z},
$$
together with governing equations
$$
\frac{d}{dt} R_{m}(t) = -k(\lambda_{1}+\lambda_2)R_{m}(t)+\lambda_{1} \sum_{j=1}^{k} R_{m-j}(t) + \lambda_{2} \sum_{j=1}^{k} R_{m+j}(t)
$$
[2003.09471]. For $k=1$ the model reduces to the classical Skellam process, while $\lambda_1=\lambda_2$ gives symmetry [2003.09471].

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_j\}$ and $\{\mu_j\}$ on $j=1,\ldots,k$ [2107.08307]. The non-homogeneous extension replaces these by deterministic time-dependent intensities $\lambda_j(t)$ and $\gamma_j(t)$ with cumulative rates $\Lambda_j(t)$ and $T_j(t)$, and defines
$$
\mathcal{S}(t)=M_1(t)-M_2(t)
$$
[2407.19227]. The resulting pgf is
$$
G_{\mathcal{S}}(u,t)=\exp\Big\{\sum_{j=1}^k \big(\Lambda_j(t)(u^{j}-1)+T_j(t)(u^{-j}-1)\big)\Big\},
$$
and the mean and variance are
$$
\mathbb{E}\mathcal{S}(t)=\sum_{j=1}^k j\big(\Lambda_j(t)-T_j(t)\big),\qquad
\mathbb{V}\mathcal{S}(t)=\sum_{j=1}^k j^2\big(\Lambda_j(t)+T_j(t)\big)
$$
[2407.19227].

An even broader formulation is the Poisson–Skellam family
$$
S(t)=\sum_{i\in I} iN_i(t),
$$
with a countable jump set $I\subset\mathbb{R}\setminus\{0\}$ and independent Poisson processes $N_i$ [2504.07672]. In the homogeneous case it is a Lévy process with exponent
$$
\psi(u)=\sum_{i\in I}\lambda_i(e^{iui}-1)
$$
and Lévy measure
$$
\nu(dx)=\sum_{i\in I}\lambda_i\delta_i(dx)
$$
[2504.07672]. This formulation includes classical, order-$k$, 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 $D_f(t)$ with Laplace transform
$$
\mathbb{E}[e^{-sD_f(t)}]=e^{-tf(s)},
$$
the time-changed order-$k$ process
$$
Z_f(t)=S^k(D_f(t))
$$
has moment generating function
$$
\mathbb{E}[e^{\theta Z_f(t)}]
=e^{-t f\left(k(\lambda_1 + \lambda_2) -\lambda_1\sum_{j=1}^{k}e^{\theta j} -\lambda_2 \sum_{j=1}^{k}e^{-\theta j}\right)}
$$
[2003.09471]. The stable choice $f(s)=s^\alpha$ yields a space-fractional model, and the tempered stable choice $f(s)=(s+\mu)^\alpha-\mu^\alpha$ yields a tempered space-fractional model [2003.09471].

The space-fractional Skellam process is defined by
$$
S_{\alpha_1,\alpha_2}(t)=N_1(D_{\alpha_1}(t))-N_2(D_{\alpha_2}(t)),
$$
with independent stable subordinators $D_{\alpha_1},D_{\alpha_2}$ [2003.09471]. Its mgf is
$$
\mathbb{E}[e^{\theta S_{\alpha_1,\alpha_2}(t)}]
= e^{-t\left[\lambda_{1}^{\alpha_{1}}(1-e^{\theta})^{\alpha_{1}}+\lambda_{2}^{\alpha_{2}}(1-e^{-\theta})^{\alpha_{2}}\right]},
$$
and the governing equations are fractional difference–differential equations involving $(1-B)^{\alpha_1}$ and $(1-F)^{\alpha_2}$ [2003.09471]. Its Lévy measure has infinite support,
$$
\nu_{S_{\alpha_1,\alpha_2}}(x)
={\lambda_1}^{\alpha_1}\sum_{n_1=1}^{\infty}(-1)^{n_1+1} {\alpha_1 \choose n_1} \delta_{n_1}(x)
+\lambda_2^{\alpha_2}\sum_{n_2=1}^{\infty}(-1)^{n_2+1} {\alpha_2 \choose n_2} \delta_{-n_2}(x),
$$
so fractionalization replaces finitely many jump sizes by infinitely many combinatorially weighted sizes [2003.09471].

The tempered space-fractional Skellam process,
$$
S^{\mu_1,\mu_2}_{\alpha_1,\alpha_2}(t)
=N_1(D_{\alpha_1,\mu_1}(t))-N_2(D_{\alpha_2,\mu_2}(t)),
$$
has mgf
$$
\mathbb{E}[e^{\theta S^{\mu_1,\mu_2}_{\alpha_1,\alpha_2}(t)}]
= e^{-t\left[\{(\lambda_1(1-e^{\theta})+\mu_1)^{\alpha_1}-\mu_1^{\alpha_1}\}+\{(\lambda_2(1-e^{-\theta})+\mu_2)^{\alpha_2}-\mu_2^{\alpha_2}\}\right]}
$$
[2003.09471]. As $\mu_i\to0$, the stable case is recovered [2003.09471].

Another generalized space-time fractional construction is the generalized space-time fractional Skellam process
$$
S_\beta^\alpha(t)=S(D_\beta(Y_\alpha(t))),
$$
where $S$ is a generalized Skellam process, $D_\beta$ is a stable subordinator, and $Y_\alpha$ is an inverse stable subordinator [2504.08374]. Its p.g.f. is
$$
G_{S_\beta^\alpha}(u,t)
=E_{\alpha,1}\!\left( -\Big(\sum_{j=1}^k (1-u^j)\lambda_j\Big)^\beta t^\alpha \right)
E_{\alpha,1}\!\left( -\Big(\sum_{j=1}^k (1-u^{-j})\mu_j\Big)^\beta t^\alpha \right),
$$
and the closed-form p.m.f. is expressed through derivatives of Mittag–Leffler functions with respect to $\Lambda=\sum_{j=1}^k\lambda_j$ and $T=\sum_{j=1}^k\mu_j$ [2504.08374]. The paper states that for $\beta\in(0,1)$, $D_\beta(t)$ has infinite moments, so mean and variance may be infinite in space-fractional cases [2504.08374].

Time change by inverse stable subordinators produces non-Lévy fractional Skellam models. The generalized fractional Skellam process
$$
S^*(t)=S(E_t),
$$
with inverse $\alpha$-stable subordinator $E_t$, has pgf
$$
G_{S^*}(u,t)=E_{\alpha,1}\!\left(t^\alpha\sum_{j=1}^k[\lambda_j(u^j-1)+\mu_j(u^{-j}-1)]\right)
$$
[2107.08307]. Its one-dimensional distributions are not infinitely divisible [2107.08307]. The non-homogeneous generalized fractional Skellam process
$$
\mathcal{S}^\alpha(t)=\mathcal{S}(Y_\alpha(t))
$$
is defined analogously, with integral representation
$$
p^\alpha(n,t)=\int_0^\infty p(n,u)h_\alpha(u,t)\,du,
$$
where $h_\alpha$ is the density of the inverse stable subordinator [2407.19227].

Compound Poisson–Gamma time changes provide another multiparameter family. If $Z_t=G(N_t)$ is a compound Poisson–Gamma subordinator with parameters $(\nu,\alpha,\beta)$, then
$$
\widetilde S_t=S_{Z_t}
$$
has characteristic function
$$
\phi_{\widetilde S}(\theta,t)
=\exp\left\{-t\nu\left[1-\left(\frac{\beta}{\beta+\lambda_1+\lambda_2-\lambda_1e^{i\theta}-\lambda_2e^{-i\theta}}\right)^\alpha\right]\right\},
$$
mean
$$
E[\widetilde S_t]=(\lambda_1-\lambda_2)\nu t(\alpha/\beta),
$$
and variance
$$
{\rm Var}(\widetilde S_t)
=(\lambda_1+\lambda_2)\nu t(\alpha/\beta)+(\lambda_1-\lambda_2)^2\nu t[\alpha(1+\alpha)/\beta^2]
$$
[1707.00523]. In contrast, inverse-time changes $S_{L_t}$ are not Lévy and have non-stationary dependent increments [1707.00523].

## 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 $N_1,N_2$ on $\mathbb{R}^M_+$, the Skellam field is
$$
S(B)=N_1(B)-N_2(B),\qquad B\subset\mathbb{R}^M,
$$
with
$$
E[S(B)]=(\lambda_1-\lambda_2)|B|,\qquad
{\rm Var}(S(B))=(\lambda_1+\lambda_2)|B|
$$
[2509.10870]. For $M=2$, the Skellam sheet $X(s,t)=N_1(s,t)-N_2(s,t)$ has covariance
$$
{\rm Cov}(X(s,t),X(s',t'))=(\lambda_1+\lambda_2)(s\wedge s')(t\wedge t')
$$
[2509.10870].

The point probabilities and generating functions retain the classical Skellam form with area or volume replacing time. For rectangles in $\mathbb{R}_+^2$,
$$
P\{X(s,t)=k\}
=e^{-(\lambda_1+\lambda_2)st}
\left(\frac{\lambda_1}{\lambda_2}\right)^{k/2}
I_{|k|}(2st\sqrt{\lambda_1\lambda_2}),
$$
and
$$
G(z;s,t)=\exp(st(\lambda_1(z-1)+\lambda_2(z^{-1}-1))).
$$
The forward equations are
$$
\partial_s p_k(s,t)=t[\lambda_1p_{k-1}(s,t)-(\lambda_1+\lambda_2)p_k(s,t)+\lambda_2p_{k+1}(s,t)],
$$
$$
\partial_t p_k(s,t)=s[\lambda_1p_{k-1}(s,t)-(\lambda_1+\lambda_2)p_k(s,t)+\lambda_2p_{k+1}(s,t)]
$$
[2509.10870]. The generator acting on $f:\mathbb{Z}\to\mathbb{R}$ is
$$
\mathcal{L}f(k)=\lambda_1(f(k+1)-f(k))+\lambda_2(f(k-1)-f(k))
$$
[2509.10870].

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
$$
N(t)\sim {\rm Poisson}(\Lambda\cdot t),\qquad \Lambda\cdot t=\sum_{k=1}^M\lambda_kt_k,
$$
and representation
$$
N(t)\overset{d}{=}\sum_{k=1}^M N_k(t_k),
$$
where the $N_k$ are independent one-parameter Poisson processes [2509.12729]. The multiparameter Skellam process is then
$$
\widetilde S(t)=N_1(t)-N_2(t),\qquad t\in\mathbb{R}_+^M,
$$
with independent multiparameter Poisson processes $N_1,N_2$ of rate vectors $\Lambda_1,\Lambda_2$ [2509.12729]. Its marginal law is
$$
P(\widetilde S(t)=k)
=\exp(-(\Lambda_1\cdot t+\Lambda_2\cdot t))
\left(\frac{\Lambda_1\cdot t}{\Lambda_2\cdot t}\right)^{k/2}
I_{|k|}\!\left(2\sqrt{(\Lambda_1\cdot t)(\Lambda_2\cdot t)}\right),
$$
with
$$
E[\widetilde S(t)]=(\Lambda_1-\Lambda_2)\cdot t,\qquad
{\rm Var}[\widetilde S(t)]=(\Lambda_1+\Lambda_2)\cdot t
$$
[2509.12729].

The generalized multiparameter Skellam process in the same framework is
$$
S(t)=\sum_{j\in\mathcal{J}} jN_j(t),
$$
for finite $\mathcal{J}\subset\mathbb{R}\setminus\{0\}$ and independent multiparameter Poisson processes $N_j$ [2509.12729]. Its pgf is
$$
E[u^{S(t)}]=\exp\left(\sum_{j\in\mathcal{J}}(\Lambda_j\cdot t)(u^j-1)\right),
$$
which is the field-indexed analogue of the one-parameter generalized Skellam pgf [2509.12729].

A related field perspective also appears in the additive Lévy framework of the generalized Poisson–Skellam family,
$$
X_t=\sum_{k=1}^m a_kN_k(t),
$$
with multiparameter Poisson drivers and fixed jumps $a_k$ [2504.07672]. The resulting random field has independent increments over disjoint rectangles and characteristic function
$$
\phi_{X_t}(u)=\exp\left(\sum_{k=1}^m\Lambda_k((0,t])(e^{iua_k}-1)\right)
$$
[2504.07672].

## 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
$$
\frac{d}{dt}s_k(t)=\lambda_1(s_{k-1}(t)-s_k(t))-\lambda_2(s_k(t)-s_{k+1}(t)),
$$
with initial condition $s_0(0)=1$, $s_k(0)=0$ for $k\ne0$ [2003.09471]. For the generalized Skellam process built from generalized counting processes, the pmf satisfies
$$
\frac{d}{dt}q(n,t)=A[q(n-1,t)-q(n,t)]-\bar A[q(n,t)-q(n+1,t)],
$$
where $A=\sum_{j=1}^k\lambda_j$ and $\bar A=\sum_{j=1}^k\mu_j$ [2107.08307].

In the order-$k$ setting, the operator structure becomes
$$
\frac{d}{dt} R_m(t)
= -\lambda_1\sum_{j=1}^{k}(1-B^{j}) R_m(t)-\lambda_2 \sum_{j=1}^{k}(1-F^{j})R_m(t),
$$
where $B$ and $F$ are backward and forward shift operators [2003.09471]. In the space-fractional setting, fractional powers of these difference operators appear:
$$
\frac{d}{dt}H_k(t)
= -\lambda_{1}^{\alpha_{1}} (1-B)^{\alpha_{1}} H_k(t)-\lambda_{2}^{\alpha_{2}} (1-F)^{\alpha_{2}} H_k(t)
$$
[2003.09471]. Tempering replaces these by $(\mu_1+\lambda_1(1-B))^{\alpha_1}-\mu_1^{\alpha_1}$ and $(\mu_2+\lambda_2(1-F))^{\alpha_2}-\mu_2^{\alpha_2}$ [2003.09471].

Inverse-stable time changes replace ordinary derivatives by Caputo derivatives. For the generalized fractional Skellam process,
$$
S^*(t)=S(E_t),
$$
the state probabilities solve a fractional equation with Caputo derivative, and the pgf is a Mittag–Leffler eigenfunction of the generator [2107.08307]. In the non-homogeneous generalized fractional Skellam process,
$$
\frac{d^\alpha}{dt^\alpha}p^\alpha(n,t)=\int_0^\infty h_\alpha(u,t)\,\frac{\partial}{\partial u}p(n,u)\,du,
$$
where $h_\alpha$ is the inverse stable density [2407.19227]. The generalized space-time fractional Skellam process has a p.g.f. governing equation involving a Caputo derivative and series coefficients $C_i,D_l$ built from $\sum_{j=1}^k (1-u^j)\lambda_j$ and $\sum_{j=1}^k (1-u^{-j})\mu_j$ [2504.08374].

Random-field variants have coupled partial difference–differential equations. For the Skellam sheet,
$$
\partial_sE[f(X(s,t))]=tE[\mathcal{L}f(X(s,t))],\qquad
\partial_tE[f(X(s,t))]=sE[\mathcal{L}f(X(s,t))]
$$
[2509.10870]. For the fractional Skellam random field of Type I,
$$
D_C^\beta D_C^\alpha G^{\alpha,\beta}(u,s,t)
=[\lambda_1(u-1)+\lambda_2(u^{-1}-1)]G^{\alpha,\beta}(u,s,t)
+\frac{(\lambda_1(u-1)u+\lambda_2(1-u))^2}{\lambda_1u^2-\lambda_2}\,\partial_u G^{\alpha,\beta}(u,s,t)
$$
[2509.10870].

A general operator-level formulation is provided by Bernstein-function subordination. If $\mathcal{A}_S$ is the Skellam generator and $\Phi(\theta)=\nu(1-(\beta/(\beta+\theta))^\alpha)$ is the Bernstein function of a compound Poisson–Gamma clock, then the subordinated semigroup has generator $\Phi(\mathcal{A}_S)$ [1707.00523]. More generally, the homogeneous Bernstein-fractional Poisson–Skellam family has transform
$$
E[e^{iuS_f(t)}]=\exp\left(-t\sum_{i\in I}f_i(\lambda_i(1-e^{iui}))\right),
$$
while inverse-Bernstein time changes lead to convolution-type derivatives $\mathcal{D}_t^f$ in the governing equations [2504.07672], [2510.12531].

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

Beyond the base processes, several papers study integral and running-average functionals. For the order-$k$ Skellam process,
$$
S_A^k(t)=\frac{1}{t}\int_0^t S^k(s)\,ds
$$
has characteristic function
$$
\phi_{S_A^k(t)}(u)=\exp\left[-kt\left\{\lambda_1\left(1-\frac{1}{k}\sum_{j=1}^k\frac{e^{iuj}-1}{iuj}\right)+\lambda_2\left(1-\sum_{j=1}^k\frac{1-e^{-iuj}}{iuj}\right)\right\}\right]
$$
and admits a compound Poisson representation with a mixed double uniform jump law [2003.09471]. Its mean and variance are
$$
\mathbb{E}[S_A^k(t)]=\frac{k(k+1)}{4}(\lambda_1-\lambda_2)t,\qquad
{\rm Var}[S_A^k(t)]=\frac{1}{18}[k(k+1)(2k+1)](\lambda_1+\lambda_2)t
$$
[2003.09471].

For the generalized space fractional Skellam process, the running average
$$
S_\beta^A(t)=\frac{1}{t}\int_0^t S_\beta(s)\,ds
$$
also admits a compound Poisson representation,
$$
S_\beta^A(t)\overset{d}{=}\sum_{i=1}^{N(D_\beta(t))}X_i,
$$
with $N(D_\beta(t))$ an SFPP and explicitly given characteristic function of $X_i$ [2504.08374]. In the generalized Skellam process with constant intensities, the running average
$$
\mathcal{S}_A(t)=\frac{1}{t}\int_0^t \mathcal{S}(s)\,ds
$$
satisfies
$$
\phi_{\mathcal{S}_A(t)}(u)
=\exp\Big\{t\Big[\sum_{j=1}^k \lambda_j\Big(\frac{e^{iuj}-1}{iuj}-1\Big)+\sum_{j=1}^k \mu_j\Big(\frac{1-e^{-iuj}}{iuj}-1\Big)\Big]\Big\}
$$
and has compound Poisson representation with mixed double uniform jumps [2407.19227].

Interacting-vector versions extend Skellam processes to dependent multicomponent point processes. In the two-group model of "Interacting point processes", the vector $(N_1,N_2)$ on $\mathbb{Z}^2$ has a migration-type infinitesimal structure and admits the decomposition
$$
N_1\overset{d}{=}S_1+S_3+S_4,\qquad
N_2\overset{d}{=}S_2+S_3-S_4,
$$
where $S_1,S_2,S_3,S_4$ are independent non-homogeneous Skellam processes with explicitly stated rate pairs [2510.12531]. The covariance is
$$
{\rm Cov}(N_1(s),N_2(t))
=\int_0^s\Big[(\lambda_1(r)-\mu_1(r))(\lambda_2(r)-\mu_2(r))-(\eta_{12}+\eta_{21})\Big]dr
$$
for $0\le s\le t$ [2510.12531]. In the homogeneous case, the vector has a compound Poisson representation with jump types $(1,1),(-1,-1),(1,-1),(-1,1),(1,0),(0,1),(-1,0),(0,-1)$ and corresponding probabilities [2510.12531].

Asymptotic and dependence properties depend strongly on the time-change mechanism. The generalized fractional Skellam process $S^*(t)=S(E_t)$ has long-range dependence, and its one-dimensional distributions are not infinitely divisible [2107.08307]. The non-homogeneous generalized fractional Skellam process inherits long and short range dependence patterns through the covariance structure of the inverse stable time change [2407.19227]. For the generalized space-time fractional Skellam process,
$$
\frac{S_\beta^\alpha(t)}{t^{\alpha/\beta}}
\Rightarrow (Y_\alpha(1))^{1/\beta}D_\beta(1)\sum_{j=1}^k j(\lambda_j-\mu_j)
$$
as $t\to\infty$, while the one-dimensional marginals are not infinitely divisible; by contrast, the generalized space fractional Skellam process with $\alpha=1$ is infinitely divisible [2504.08374].

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 $J=\{\pm1\}$ [2509.10870]. A parallel triangular-array weak convergence result is proved for the alternative multiparameter construction
$$
\mathcal{S}(t)=\sum_{j\in\mathcal{J}}jN_j(t_j)
$$
in the additive multiparameter Poisson framework [2509.12729].

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 [2003.09471], [2509.10870], [2504.08374]. 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.

Source: https://www.emergentmind.com/topics/multiparameter-skellam-process