This is a generalization of the question Solving another non-trivial recurrence relation. Let $\lambda^C \ge 0$, $\lambda^M \ge 0$ and $\Lambda \ge 0$ and $q\in (0,1)$..Without loss of generality we can set $\lambda^C = 1$. Now, we consider a following recurrence relation: \begin{eqnarray} 0 &=& {\mathbb P}(n+1) \cdot \left( \lambda^{M} +{\mathcal C}^{(\theta)}(n+1) \cdot \lambda^{C}\right) + \\ &&{\mathbb P}(n) \cdot \left(-\lambda^{M} 1_{n \ge 1} - {\mathcal C}^{(\theta)}(n) \cdot \lambda^{C} - \Lambda \right) +\\ && \Lambda \cdot \sum\limits_{i=1}^n {\mathbb P}(n-i) \cdot q \cdot (1-q)^{i-1} \quad (i) \end{eqnarray} Again, as considered in 1 the equation above describes a steady state of a queuing system with orders conforming to three Poisson processes (cancellations, limit sell orders and market buy orders) and limit sell order sizes conforming to a negative binomial distribution with parameter $q$. Here the likelihood of a cancellation is given by ${\mathcal C}^{(\theta)}(n) \cdot \lambda^C$ were ${\mathcal C}^{(\theta)}(n)$ is a polynomial in $n$ of order $\theta$.
Now, we have found the solutions to the equations above. They read as follows: \begin{eqnarray} {\mathfrak N} \cdot {\mathbb P}(n) = (1-q)^n \cdot \prod\limits_{\xi=0}^{\theta-1} \frac{\left(0-\zeta_\xi(\lambda^M + \frac{\Lambda}{1-q})\right)^{(n)}}{\left(1-\zeta_\xi(\lambda^M )\right)^{(n)}} \cdot \left[ 1_{n=0} + \frac{\Lambda}{\Lambda + \lambda^M \cdot (1-q)} \cdot 1_{n\ge 1} \right] \quad (ii) \end{eqnarray} where $\left\{ \zeta_\xi(x) \right\}_{\xi=0}^{\theta-1}$ are roots of a polynomial equation, i.e.: \begin{eqnarray} {\mathcal C}( \zeta_\xi(x) ) + x = 0 \quad \mbox{for $\xi=0,\cdots,\theta-1$} \end{eqnarray} and the normalization constant ${\mathfrak N}$ reads: \begin{eqnarray} {\mathfrak N} := q \cdot F_{\theta+1,\theta} \left[ \begin{array}{rr} 1 & \left( 1- \zeta_\xi(\lambda^M + \frac{\Lambda}{1-q}) \right)_{\xi=0}^{\theta-1} \\ & \left( 1- \zeta_\xi(\lambda^M) \right)_{\xi=0}^{\theta-1} \end{array} ;1-q \right] \end{eqnarray} and $F_{\theta+1,\theta}\left[\right]$ is the generalized hypergeometric function.
Below is a Mathematica code that verifies the solution.
In[1597]:= (*As above where the likelihhod of cancellations is a \
polynomial of order th.*)
{lM, L, q} = RandomReal[{0, 1}, 3, WorkingPrecision -> 50];
M = 100; XX = Table[0, {M}]; XX[[1]] = 1; th = 5; n =.;
coeffs = RandomInteger[{1, 10}, th];
cLHd[n_] :=
n^Range[1,
th] .coeffs;(*Cancellation likelihood as a function of pending \
orders in the book.*)
Print[cLHd[n]];
XX = Table[
If[n == 0,
1, (1 - q)^n Product[L/(1 - q) + lM + cLHd[j], {j, 0, n - 1}]/
Product[lM + cLHd[j], {j, 1, n}] L/( L + lM (1 - q)) ], {n, 0,
M - 1}];
x =.;
rts0 = x /. NSolve[cLHd[x] + lM == 0, x];
rts1 = x /. NSolve[cLHd[x] + L/(1 - q) + lM == 0, x];
XX0 = Table[
If[n == 0,
1, (1 - q)^
n Product[Pochhammer[0 - rts1[[xi + 1]], n], {xi, 0, th - 1}]/
Product[Pochhammer[1 - rts0[[xi + 1]], n], {xi, 0, th - 1}] L/(
L + lM (1 - q)) ], {n, 0, M - 1}];
Total[(XX0 - XX)^2]
(*Check if the master equations are satisfied.*)
Table[ (lM + cLHd[n + 1] lC) XX[[2 + n]] -
XX[[1 + n]] (If[n == 0, 0, lM] + cLHd[n] lC + L + EE) +
L Sum[ XX[[1 + n - i]] q (1 - q)^(i - 1), {i, 1, n}], {n, 0,
Length[XX] - 2}]
NN = q HypergeometricPFQ[
Join[{1}, Table[1 - rts1[[xi + 1]], {xi, 0, th - 1}]],
Table[1 - rts0[[xi + 1]], {xi, 0, th - 1}], 1 - q];
Total[XX]/NN
During evaluation of In[1597]:= 4 n+8 n^2+6 n^3+2 n^4+8 n^5
Out[1606]= 0.*10^-105 + 0.*10^-105 I
Out[1607]= {0.*10^-52, 0.*10^-52, 0.*10^-52, 0.*10^-53, 0.*10^-53,
0.*10^-54, 0.*10^-55, 0.*10^-56, 0.*10^-56, 0.*10^-57, 0.*10^-58,
0.*10^-58, 0.*10^-59, 0.*10^-60, 0.*10^-61, 0.*10^-61, 0.*10^-62,
0.*10^-63, 0.*10^-64, 0.*10^-64, 0.*10^-65, 0.*10^-66, 0.*10^-67,
0.*10^-67, 0.*10^-68, 0.*10^-69, 0.*10^-70, 0.*10^-70, 0.*10^-71,
0.*10^-72, 0.*10^-73, 0.*10^-73, 0.*10^-74, 0.*10^-75, 0.*10^-76,
0.*10^-77, 0.*10^-77, 0.*10^-78, 0.*10^-79, 0.*10^-80, 0.*10^-80,
0.*10^-81, 0.*10^-82, 0.*10^-83, 0.*10^-83, 0.*10^-84, 0.*10^-85,
0.*10^-86, 0.*10^-86, 0.*10^-87, 0.*10^-88, 0.*10^-89, 0.*10^-90,
0.*10^-90, 0.*10^-91, 0.*10^-92, 0.*10^-93, 0.*10^-93, 0.*10^-94,
0.*10^-95, 0.*10^-96, 0.*10^-96, 0.*10^-97, 0.*10^-98, 0.*10^-99,
0.*10^-99, 0.*10^-100, 0.*10^-101, 0.*10^-102, 0.*10^-103,
0.*10^-103, 0.*10^-104, 0.*10^-105, 0.*10^-106, 0.*10^-106,
0.*10^-107, 0.*10^-108, 0.*10^-109, 0.*10^-109, 0.*10^-110,
0.*10^-111, 0.*10^-112, 0.*10^-113, 0.*10^-113, 0.*10^-114,
0.*10^-115, 0.*10^-116, 0.*10^-116, 0.*10^-117, 0.*10^-118,
0.*10^-119, 0.*10^-119, 0.*10^-120, 0.*10^-121, 0.*10^-122,
0.*10^-123, 0.*10^-123, 0.*10^-124, 0.*10^-125}
Out[1609]= 1.000000000000000000000000000000000000000000000000 +
0.*10^-50 I
Below we also plot the distributions in question as a function along with their 99.9 percentiles. Here we took ${\mathcal C}^{(\theta)}(n) = n^\theta$ and $(\lambda^M,q,\Lambda) = (0.1100,0.893540,4.84097)$. We have:
On the left hand side we have a family of distributions corresponding to $\theta = 0.5 + i/19$ for $i=0,\cdots,19$ (from Violet to Red respectively) whereas on the right hand side we plotted the 99.9 percentile of those distributions (Blue) as a function of the parameter $\theta$ along with a power-law fit $A\cdot \theta^b$ with parameters $(A,b)$ given in the label.
Now, having said all this my question would be the following. Clearly the quantity $\left( {\mathbb P}(n) \right)_{n=0}^\infty $ is a probability distribution (the steady state of the Markov chain in question). As such this probability distribution has moments (mean, variance, skewness, kurtosis). Can we compute those in closed form ?
1 "The Order Book as a Queueing System" in: F Abergel et al, Limit Order Books, Physics of Society: Econophysics and Sociophysics, Cambridge University Press 2016
Now we proceed to computing the moments of the distribution in question. Assume for the time being that $\theta \in {\mathbb N}$ and $\theta \ge 1$ however we will waive this assumption at the end. From $(ii)$ we can immediately write dwon the Z transform $X(z) := \sum\limits_{n=0}^\infty {\mathbb P}(n) \cdot z^n$ of the distribution. We have: \begin{eqnarray} X(z)&=& \frac{\lambda^M \cdot (1-q)}{\Lambda + \lambda^M \cdot (1-q)} + \frac{\Lambda}{\Lambda + \lambda^M \cdot (1-q)} \cdot F_{\theta+1,\theta} \left[ \begin{array}{rr} 1 & - \zeta_\xi (\lambda^M + \Lambda/(1-q))^{1/\theta} \\ & 1- \zeta_\xi (\lambda^M)^{1/\theta} \end{array}; (1-q) \cdot z \right] \\ &=& \left(1-(1-q) \cdot z\right) \cdot F_{\theta+1,\theta} \left[ \begin{array}{rr} 1 & 1 - \zeta_\xi (\lambda^M + \Lambda/(1-q))^{1/\theta} \\ & 1- \zeta_\xi (\lambda^M)^{1/\theta} \end{array}; (1-q) \cdot z \right] \end{eqnarray} Here $\zeta_\xi := \exp(\imath \pi/\theta (1+2 \xi))$ is the $\theta$-th root of minus unity and $F_{\theta+1,\theta}\left[\right]$ is the generalized hypergeometric function.
Now, clearly the moments are computed by differentiating the Z transform in question at $z=1$. We have: \begin{eqnarray} E\left[ A^n \right] &=& \frac{1}{X(1)} \cdot \left. \frac{d^n}{d z^n} X(z) \right|_{z=1} \\ &=& \frac{\Lambda}{\Lambda + \lambda^M \cdot (1-q)} \cdot \frac{1}{X(1)} \cdot \frac{d^n}{d z^n} \left. F_{\theta+1,\theta} \left[ \begin{array}{rr} 1 & - \zeta_\xi (\lambda^M + \Lambda/(1-q))^{1/\theta} \\ & 1- \zeta_\xi (\lambda^M)^{1/\theta} \end{array}; (1-q) \cdot z \right] \right|_{z=1} \\ &=& \frac{n!}{X(1)} \cdot \frac{\Lambda}{1+\lambda^M} \cdot (1-q)^{n-1} \cdot \prod\limits_{p=2}^n \left(\frac{(p-1)^\theta +\lambda^M + \Lambda/(1-q)}{p^\theta + \lambda^M}\right) \cdot F_{\theta+1,\theta} \left[ \begin{array}{rr} n+1 & n - \zeta_\xi (\lambda^M + \Lambda/(1-q))^{1/\theta} \\ & n+1- \zeta_\xi (\lambda^M)^{1/\theta} \end{array}; (1-q) \right] \end{eqnarray} Here $n=1,2,3,\cdots$. Now even though we assumed $\theta $ to be an integer the expression above can be generalized to any real value of $\theta$ by replacing the hypergeometric function above by the following sum: \begin{equation} F_{\theta+1,\theta} \left[ \right] = \sum\limits_{m=0}^\infty \frac{\prod\limits_{j=0}^{n-1}\left( (j+n)^\theta +\lambda^M + \Lambda/(1-q)\right)}{\prod\limits_{j=0}^{n-1}\left( (j+n+1)^\theta + \lambda^M\right) } \cdot \binom{n+m}{m} \cdot (1-q)^m \end{equation}