# Evaluating Approximations of Count Distributions and Forecasts for Poisson-
Lindley Integer Autoregressive Processes

Rachel D. Gidaro111Department of Mathematical Sciences, United States Military
Academy, West Point e-mail<EMAIL_ADDRESS>
Department of Statistical Science, Baylor University, Waco, TX 76798-7140
Jane L. Harvill  e-mail<EMAIL_ADDRESS>
Department of Statistical Science, Baylor University, Waco, TX 76798-7140
ORCID ID: #0000-0001-8585-5945

(June 2024
Word Count: 3861)

###### Abstract

Although many time series are realizations from discrete processes, it is
often that a continuous Gaussian model is implemented for modeling and
forecasting the data, resulting in incoherent forecasts. Forecasts using a
Poisson-Lindley integer autoregressive (PLINAR) model are compared to
variations of Gaussian forecasts via simulation by equating relevant moments
of the marginals of the PLINAR to the Gaussian AR. To illustrate utility, the
methods discussed are applied and compared using a discrete series with model
parameters being estimated using each of conditional least squares, Yule-
Walker, and maximum likelihood.

Keywords: Discrete time series, coherent forecasts, Poisson-Lindley integer
autoregressive process, Kullback-Leibler divergence, Kolmogorov metric.

##  1 Introduction

Time series analysis is used across a variety of areas of scientific inquiry,
from business to health sciences, physics, geology, and social sciences, to
name a few. In many circumstances, researchers look to model a process over
time and ultimately forecast future events. Traditional time series analysis
uses methods based on continuous processes, such as an autoregressive-moving
average (ARMA) with Gaussian innovations, for analysis and forecasting, even
when the process is discrete. In fact many time series textbooks commonly do
this, see Box and Jenkins, (1976), Shumway and Stoffer, (2017), and Brockwell
and Davis, (1996).

One desirable property of a time series model is that it produces forecasts
that do not violate known constraints; that is, the model produces “coherent
forecasts.” Traditional application of a continuous ARMA model for forecasting
a discrete process will result in violating coherency. Thus, coherent
forecasting is one area where discrete modeling can yield improvements over
the traditional continuous models. In a series of seminal papers (Jacobs and
Lewis, 1978a, ; Jacobs and Lewis, 1978b, ; Jacobs and Lewis,, 1983), Jacobs
and Lewis introduced the theory and analysis for discrete time series
processes. The disadvantages of these early models are well-documented; see,
for example Al-Osh and Alzaid, (1987); Alzaid and Al-Osh, (1988).

A collection of discrete time series models that addresses some of the
disadvantages of Jacobs and Lewis incorporates binomial thinning. The first of
these applies binomial thinning to a process with Poisson marginals, yielding
the first-order integer-valued autoregressive (INAR(1)) model (Al-Osh and
Alzaid,, 1987; McKenzie,, 1988). To address overdispersion, Alzaid and Al-Osh,
(1993) developed an INAR(1) process with generalized Poisson marginals. Other
approaches to addressing overdispersion apply binomial thinning to processes
with geometric or negative binomial marginals (McKenzie,, 1986). The Poisson-
Lindley distribution was introduced by Sankaran, (1970) to model count data.
This distribution is unimodal, has larger variability than the traditional
Poisson distribution, and has smaller skewness and kurtosis than the negative
binomial model (Ghitany and Al-Mutairi,, 2009). In the time series framework,
Mohammadpour et al., (2018) introduced the INAR(1) process with Poisson-
Lindley marginals. One downside to the Poisson-Lindley INAR(1) model is that
it can be complex and difficult to estimate the model parameters.

Many papers have looked at addressing the efficacy of Gaussian methods applied
to discrete time series. In Maiti et al., (2016), various INAR(1) processes
were compared to the Gaussian autoregressive process, where the forecast for
the Gaussian model was rounded to the nearest integer. They found that the
Gaussian autoregressive process produced similar results to that of the
INAR(1) process. However, in Alwan and Weiß, (2017) the marginal forecast
distribution of INAR(1) processes was approximated by the Gaussian AR(1)
conditional forecast distribution, and it was found that the approximation was
worse, especially in the range of low counts data. Additionally Bisaglia and
Gerolimetto, (2015) showed when systematically comparing a coherent point
forecast method of an INAR(1) model to the approximate point forecast obtained
from a Gaussian AR(1), the Gaussian had higher forecast mean square error and
forecast mean absolute error. Homburg, (2020) thoroughly investigated the
differences between the Poisson INAR(1) to the Gaussian AR(1) using a variety
of measures to compare the forecasts from the competing models. Her main
conclusion is the quality of the approximation is influenced by the amount of
autocorrelation present. Finally, Gidaro and Harvill, (2024) provide a
thorough investigation of forecasting discrete time series using a GINAR(1).
In their work, they compare the conditional distributions for a discretized
Gaussian AR(1) from two different moment-equated methods to the GINAR(1)
conditional distribution. They also compare the forecasts of the different
models that takes into account different parameter estimation techniques. They
find that in some circumstances, the discretized Gaussian forecasts performs
as well as the forecasts based on the GINAR model, but in most cases, the
GINAR outperforms the Gaussian AR forecasts.

In what follows, we investigate the performance of Poisson-Lindley forecasts
compared to their Gaussian counterpart. The rest of the paper is organized as
follows. In Section 2 the first order Poisson-Lindley integer autoregressive
process (PLINAR(1)), first introduced by Mohammadpour et al., (2018), is
defined. In Section 3, we present the $k$-step-ahead conditional mean
forecasts for a PLINAR(1) process and for a traditional Gaussian first order
autoregressive (AR(1)) model. In Section 4, we develop two methods for
equating moments of the PLINAR(1) to the Gaussian AR(1). The closeness of the
conditional PMFs for the forecasting methods are compared in Section 4 via
Kullback-Leibler divergence and Kolmogorov metric. To illustrate utility and
to further investigate efficacy, the methods are then applied to a discrete
time series data set. Prior to that, in Section 6, three methods for parameter
estimation are presented along with criteria for comparing the forecasts.
Concluding remarks are in Section 8.

##  2 Poisson-Lindley Integer-Valued Autoregressive Processes

The Poisson-Lindley integer-valued autoregressive model of order one, or
PLINAR(1), first proposed by Mohammadpour et al., (2018) is defined as
follows. Let $B_{j}$ be independent and identically distributed
Bernoulli($\alpha$) random variables where for each $j,~{}P(B_{j}=1)=\alpha$.
The binomial thinning operator $\circ$ is

$\alpha\circ X=\sum_{j=1}^{X}B_{j}.$

(Steutel and van Harn,, 1979). The process $\\{X_{t}\\}$ follows a PLINAR(1)
if

$X_{t}=\alpha\circ X_{t-1}+I_{t}H_{t},\qquad t=1,2,\ldots,\quad
0\leq\alpha\leq 1,$

where $I_{t}$ has $P(I_{t}=0)=1-P(I_{t}=1)=\alpha$, and for all $t,~{}H_{t}$
has probability mass function

$\displaystyle g(h)=P(H_{t}=h)$
$\displaystyle=\frac{\theta^{2}(1-\alpha)^{2}+\theta(1-\alpha^{2})+2\alpha}{\left\\{\theta(1-\alpha)+1)^{2}\right\\}}\frac{\theta}{1+\theta}\left(1-\frac{\theta}{1+\theta}\right)$
$\displaystyle\qquad+\frac{(1-\alpha)}{\theta(1-\alpha)+1}(h+1)\left(\frac{\theta}{1+\theta}\right)^{2}\left(1-\frac{\theta}{1+\theta}\right)^{h}$
$\displaystyle\qquad-\frac{\alpha}{\left\\{\theta(1-\alpha)+1\right\\}^{2}}\frac{\theta+1}{\theta+1+\alpha}\left(1-\frac{\theta+1}{\theta+1+\alpha}\right)^{h},\qquad
h=0,1,2,\ldots.$

It is further assumed that $X_{0}$ is a PL($\theta$) random variable.
Sankaran, (1970) showed that the probability mass function (PMF) can be
written as follows:

$P(X_{0}=x)=\frac{\theta^{2}(x+\theta+2)}{(1+\theta)^{x+3}},\qquad
x=0,1,2,\ldots,\quad\theta>0.$

Then, the mean and variance of the Poisson-Lindley($\theta$) distribution are
$(\theta+2)/\\{\theta(\theta+1)\\}$ and
$(\theta^{3}+4\theta^{2}+6\theta+2)/\left\\{\theta^{2}(\theta+1)^{2}\right\\}$,
respectively. Since the marginal distributions of the $X_{t}$ are Poisson-
Lindley, it follows that

${\mbox{E}}\left[X_{t}\right]={\mbox{E}}\left[H_{t}\right]=\frac{\theta+2}{\theta(\theta+1)}\quad{\text{and}}\quad{\mbox{Var}}\left(X_{t}\right)=\frac{\theta^{3}+4\theta^{2}+6\theta+2}{\theta^{2}(\theta+1)^{2}}.$

##  3 Forecasting a PLINAR Process

Traditionally, given a series of length
$n,~{}\mathbf{X}_{n}=(X_{n},\ldots,X_{1})$, the $k$-step-ahead forecast is the
conditional mean ${\mbox{E}}\left[X_{n+k}\,|\,\mathbf{X}_{n}\right]$. The
$k$-step ahead conditional PMF for the PLINAR(1), given by Mohammadpour et
al., (2018) is

$P(X_{n+k}=y\,|\,\mathbf{X}_{n},X_{n}=x)=\sum_{j=0}^{\min(y,x)}\binom{x}{j}\alpha^{kj}(1-\alpha^{k})^{x-j}P(Z_{n+k}=y-j),$

where

$P(Z_{n+k}=z)=\left\\{\begin{array}[]{ll}\alpha^{k}+(1-\alpha^{k})\left[A_{k}\frac{\theta}{1+\theta}+B_{k}\left(\frac{\theta}{1+\theta}\right)^{2}+C_{k}\frac{1+\theta}{1+\theta+\alpha^{k}}\right]&{\mbox{if
$z=0$,}}\\\ \\\
(1-\alpha^{k})\left[A_{k}\frac{\theta}{1+\theta}\left(\frac{1}{1+\theta}\right)^{z}+B_{k}(z+1)\left(\frac{\theta}{1+\theta}\right)^{2}\left(\frac{1}{1+\theta}\right)^{z}+C_{k}\frac{1+\theta}{1+\theta+\alpha^{k}}\left(\frac{\alpha^{k}}{1+\theta+\alpha^{k}}\right)^{z}\right]&{\mbox{if
$z=1,2,\ldots,$}}\\\ \end{array}\right.$

and

$A_{k}=\frac{\theta^{2}(1-\alpha^{k})^{2}+\theta(1-\alpha^{2k})+2\alpha^{k}}{\left\\{\theta(1-\alpha^{k})+1\right\\}^{2}},\qquad
B_{k}=\frac{1-\alpha^{k}}{\theta(1-\alpha^{k})+1},\qquad
C_{k}=\frac{-\alpha^{k}}{\left\\{\theta(1-\alpha^{k})+1\right\\}^{2}}.$

The $k$-step ahead conditional mean and variance provided by Mohammadpour et
al., (2018) are, respectively

$\displaystyle{\mbox{E}}\left[X_{n+k}\,|\,\mathbf{X}_{n},X_{n}=x\right]$
$\displaystyle=\alpha^{k}x+(1-\alpha^{k})\frac{\theta+2}{\theta(\theta+1)}$
(1) $\displaystyle{\mbox{Var}}\left(X_{n+k}\,|\,\mathbf{X}_{n},X_{n}=x\right)$
$\displaystyle=\alpha^{k}(1-\alpha^{k})x+\left(\frac{1-\alpha^{2k}}{1-\alpha^{2}}\right)\sigma^{2}_{\epsilon}+\left[\frac{(1-\alpha^{k})(\alpha-\alpha^{k+1})}{1-\alpha^{2}}\right]\mu_{\epsilon},$

where $\mu_{\epsilon}$ and $\sigma^{2}_{\epsilon}$ are the marginal mean and
variance of $\epsilon_{t}$ given by

$\mu_{\epsilon}=\frac{(1-\alpha)(\theta+2)}{\theta(\theta+1)}\quad{\mbox{and}}\quad\sigma^{2}_{\epsilon}=\frac{(1-\alpha)\left[\theta^{3}+4\theta^{2}+6\theta+2+\alpha(\theta^{2}+4\theta+2)\right]}{\theta^{2}(\theta+1)^{2}}.$

As
$k\rightarrow\infty,{\mbox{E}}\left[X_{n+k}\,|\,\mathbf{X}_{n}\right]\rightarrow{\mbox{E}}\left[X_{n}\right]$
and
${\mbox{Var}}\left(X_{n+k}\,|\,\mathbf{X}_{n}\right)\rightarrow{\mbox{Var}}\left(X_{n}\right)$.

As noted earlier, it is common practice to use a Gaussian autoregressive
process (AR) for modeling and forecasting discrete-valued time series. If
$\\{W_{t}\\}$ denotes a Gaussian autoregressive process of order $p$, AR($p$),
$W_{t}$ can be written

$W_{t}=\phi_{0}+\phi_{1}W_{t-1}+\cdots+\phi_{p}W_{t-p}+\varepsilon_{t},\quad
t=1,2,\ldots,$

where $\varepsilon_{t}$ are IID normal random variables with finite mean
$\mu_{\varepsilon}$ and finite variance $\sigma^{2}_{\varepsilon}$,
independent of the $\\{W_{s}:\,s<t\\}$. If the mean of the innovations
${\mbox{E}}\left[\varepsilon_{t}\right]=0$, then $\phi_{0}=0$. Assuming the
process is second-order stationary, the autocovariance function of $W_{t}$
satisfies the Yule-Walker equations. In particular, for the AR(1) with
$\phi=\phi_{1}$, the marginal mean and variance of $W_{t}$ are

${\mbox{E}}\left[W_{t}\right]=\frac{\mu_{\varepsilon}}{1-\phi}\quad{\mbox{and}}\quad{\mbox{Var}}\left(W_{t}\right)=\frac{\sigma^{2}_{\varepsilon}}{1-\phi^{2}}.$

Let ${\mathbf{W}}_{n}=(W_{n},\ldots,W_{1})$. Then for $W_{n}=w$, the $k$-step-
ahead forecast is the $k$-step conditional mean given by

${\mbox{E}}\left[W_{n+k}\,|\,{\mathbf{W}}_{n},W_{n}=w\right]=\phi^{k}w,$

the $k$-step-ahead conditional distribution is normal,

$W_{n+k}\,|\,\mathbf{W}_{n}\sim{\mbox{N}}\left(\mu_{W_{n+k}\,|\,\mathbf{W}_{n}},\sigma^{2}_{W_{n+k}\,|\,\mathbf{W}_{n}}\right),$

where for $W_{n}=w$

$\displaystyle\mu_{W_{n+k}\,|\,\mathbf{W}_{n}}$
$\displaystyle=\phi^{k}w+\mu_{\varepsilon}\frac{1-\phi^{k}}{1-\phi}$ (2)
$\displaystyle\sigma^{2}_{W_{n+k}\,|\,\mathbf{W}_{n}}$
$\displaystyle=\left(\frac{1-\phi^{2k}}{1-\phi^{2}}\right)\sigma^{2}_{\varepsilon},$
(3)

##  4 Comparing PLINAR(1) to Gaussian AR(1) Predictions

In order to investigate how well the Gaussian AR(1) process can be used for
$k$-step-ahead forecasting of a PLINAR(1) process, two approaches were
considered for equating first- and second-moments. The first approach, which
we refer to as the “innovation method,” equates the means and the variances of
the innovation series; the second equates the means and variances of the
marginal distributions. We refer to the second method as the “marginal
method.” For both, the Gaussian AR(1) coefficient $\phi$ is set equal to the
PLINAR(1) thinning parameter $\alpha$. We note that since $0<\alpha<1$, the
Gaussian AR(1) will be second-order stationary.

To equate the mean and variance of the innovations, let

$\displaystyle{\mbox{E}}\left[\varepsilon\right]$
$\displaystyle={\mbox{E}}\left[\epsilon\right]$
$\displaystyle{\mbox{Var}}\left(\varepsilon\right)$
$\displaystyle={\mbox{Var}}\left(\epsilon\right)$
$\displaystyle\mu_{\varepsilon}$
$\displaystyle=\frac{(1-\alpha)(\theta+2)}{\theta(\theta+1)}$
$\displaystyle\sigma^{2}_{\varepsilon}$
$\displaystyle=\frac{(1-\alpha)\left[\theta^{3}+4\theta^{2}+6\theta+2+\alpha(\theta^{2}+4\theta+2)\right]}{\theta^{2}(\theta+1)^{2}}.$
(4)

Thus, for the innovations method, $\mu_{\varepsilon}$ and
$\sigma^{2}_{\varepsilon}$ in (2) and (3) are replaced by the corresponding
expressions in (4) for computing $k$-step-ahead forecasts using the Gaussian
AR(1).

For the marginal method, we first find an expression for $\mu_{\varepsilon}$
so that the marginal means of $X_{t}$ and $W_{t}$ are equal. This gives

$\displaystyle{\mbox{E}}\left[W_{t}\right]$
$\displaystyle={\mbox{E}}\left[X_{t}\right]$
$\displaystyle\frac{\mu_{\varepsilon}}{1-\alpha}$
$\displaystyle=\frac{\theta+2}{\theta(\theta+1)}$
$\displaystyle\mu_{\varepsilon}$
$\displaystyle=\frac{(1-\alpha)(\theta+2)}{\theta(\theta+1)},$ (5)

which is equal to the expression for $\mu_{\varepsilon}$ in (4). However, to
equate the marginal variances, $\sigma^{2}_{\varepsilon}$ is

$\displaystyle{\mbox{Var}}\left(W_{t}\right)$
$\displaystyle={\mbox{Var}}\left(X_{t}\right)$
$\displaystyle\frac{\sigma^{2}_{\varepsilon}}{1-\alpha^{2}}$
$\displaystyle=\frac{\theta^{3}+4\theta^{2}+6\theta+2}{\theta^{2}(\theta+1)^{2}}$
$\displaystyle\sigma^{2}_{\varepsilon}$
$\displaystyle=\frac{(1-\alpha^{2})(\theta^{3}+4\theta^{2}+6\theta+2)}{\theta^{2}(\theta+1)^{2}}$
(6)

These expressions are used in (2) and (3) to find $k$-step-ahead forecasts
using the Gaussian AR(1). Note that denominators in the expressions for the
variances in (4) and (6) are equal. It can be shown that for all $\theta>0$,
the numerator of the variance for the innovation method is less than the
numerator of the variance for the marginal method. Specifically, as $\alpha$
gets closer to one and as $\theta$ increases the marginal variance becomes
larger faster than the innovation variance. Consequently, there is less
variability in the forecasts computed when equating innovation moments than in
equating marginal moments. For values considered in our simulation study in
Section 5.3, the ratio of the variances ranges from 1.069444 when $\alpha=0.1$
and $\theta=0.5$, to 1.780168 when $\alpha=0.9$ and $\theta=5$.

##  5 Assessing Forecasts

In this section, we investigate how well the conditional PMF based on the
Gaussian AR(1) compares to the conditional PLINAR(1) PMF. Our investigation is
similar to that of Gidaro and Harvill, (2024) and Homburg, (2020). Homburg,
considers a Poisson AR(1) and sets $X_{n}=$ round($\mu$). Here, for the
Poisson-Lindley AR(1), we consider the effects of a variety of combinations of
values of $\alpha,\theta$, and $X_{n}$. Additionally, our target distribution
is from a PLINAR(1), compared to Poisson. Before presenting the results of the
study, a brief overview of the Kullback-Leibler divergence and the Kolmogorov
metric is provided.

###  5.1 Kullback-Leibler Divergence

The Kullback-Leibler divergence (KLD) is a commonly used statistical distance
that measures the dissimilarity between two different probability
distributions and is defined as

$KL_{p,q}=\sum_{y}p(y)\log\left(\frac{p(y)}{q(y)}\right)$

where $p(y)$ is the target distribution function and $q(y)$ is the
approximation of $p(y)$. The Kullback-Leibler divergence is a nonsymmetric
measure of the aggregate disagreement of two distribution functions. Exact
agreement will result in $KL_{p,q}=0$. Larger values of $KL_{p,q}$ imply the
two distributions are more dissimilar. For the purposes of this study, we note
that $p(y)$ is the conditional PMF given in (3). Since the conditional
distribution of the Gaussian AR(1) is zero for all $y=0,1,\ldots$, an
equivalent measure to (3) must be found for $q(y)$. To accomplish this, let
$\Phi(\cdot)$ represent the standard normal CDF, and define $\tilde{G}(\cdot)$
as

$\tilde{G}(y)=\Phi\left(\frac{y-\mu}{\sigma}\right),\quad{\mbox{for all
$y=0,1,\ldots$,}}$

where expressions for the mean $\mu$ and the standard deviation
$\sigma=\sqrt{\sigma^{2}}$ are provided in (4), or (5) and (6), depending upon
which approximation method is being evaluated. Thus $q(y)$ is defined as

$q(y)=\tilde{G}(y)-\tilde{G}(y-1)\quad{\text{for $y=1,2,\ldots.$}}$

At $y=0,~{}q(0)=\tilde{G}(0)$.

###  5.2 Kolmogorov Metric

Although the Kullback-Leibler divergence is an effective method for
quantifying the dissimilarity between two different probability measures, it
is not a metric; therefore, $KL_{p,q}\neq KL_{q,p}$. For this reason, we also
consider the Kolmogorov metric which is the maximum distance between any two
CDFs and is defined as

$K_{F,G}=\max_{y}\bigl{|}F(y)-G(y)\bigr{|}.$

Interpreting the Kolmogorov metric is straight-forward compared interpreting
$KL_{p.q}$. Simply put, the Kolmogorov metric is the maximum distance between
the target cumulative distribution function $F(\cdot)$ from the approximating
$G(\cdot)$ across all values of $y$. In the context of comparing forecasts,
$F(\cdot)$ is the conditional cumulative distribution function (CDF) for the
PLINAR(1) process based on the conditional PMF in (3) and $G(\cdot)$ is the
conditional CDF for the Gaussian AR(1) process given in (3) for $k=1$.

###  5.3 Simulation Study

A simulation study was conducted to compare the performances of one-step-ahead
($k=1$) forecasting using a PLINAR(1) model versus a Gaussian AR(1).
Comparisons were based on $KL_{p.q}$ and $K_{F,G}$ for all combinations of
$\alpha=0.1,0.5$ and 0.9, and $\theta=0.5,2$, and 5. Results are visually
displayed in Figures 1 through 18. The figures are organized as follows. Each
figure contains three graphs that show the Kullback-Leibler divergence or the
Kolmogorov metric for combinations of values of $\alpha,\theta$, and $X_{n}$.
On each graph, are two curves. The black curve is the $KL_{p,q}$ or $K_{F,G}$
for the marginal method; the red is for the innovation method.

Figures 1 and 2 contain the two distances for $\theta=0.5$, $X_{n}=0,15$, and
30, and $\alpha$ ranging from 0.05 to 0.95 by increments of 0.01. For
$\theta=0.5$, we note that ${\mbox{E}}\left[X_{n}\right]=3.3333$, which is why
we consider larger values for $X_{n}$. We notice that the marginal and
innovation methods perform equally well. In this case, because the prediction
variance of the innovation method is smaller than the prediction variance of
the marginal method, the innovation method would be preferred.

Figures 3 and 4 provide similar information for $\theta=2$ and $X_{n}=0,2$,
and 5, and Figures 5 and 6 for $\theta=5$ and $X_{n}=0,2$, and 5. Maximum
values for $X_{n}$ were chosen based on the mean
${\mbox{E}}\left[X_{n}\right]$. Overall, the conclusions are the same.
However, we note that values of the Kullback-Leibler divergences are larger as
$\theta$ and $\alpha$ increase. So the marginal prediction method may be
preferred in these cases even though it has a larger variance.

A prominent feature of the graphs containing the Kolmogorov distance is the
increasing and decreasing pattern. This pattern is governed by values of
$\theta,\alpha$ and $X_{n}$. Specifically, (2) shows the conditional mean of
the forecast using the Gaussian approximation includes an additive term
$\mu_{\varepsilon}=(1-\alpha)(\theta+2)/(\theta(\theta+1))$. Since $k=1$, the
factor $(1-\phi^{k})/(1-\phi)=1$. When $\theta$ is larger, the local minimums
of the Kolmogorov distances occur close to integer values of $\alpha X_{n}$.
However, as $\theta$ becomes smaller, $\mu_{\varepsilon}$ becomes more heavily
weighted in the conditional mean.

Figures 7 through 12 contain the Kullback-Leibler and Kolmogorov values for
values of $\alpha=0.1,0.5$, and 0.9 and values of $X_{n}$ determined by
${\mbox{E}}\left[X_{n}\right]$ as for the previous study. However, in the
following figures the value of $\theta$ ranges from 0.05 to 5 in increments of
0.01. When comparing the distances via $KL_{p,q}$ and $K_{F,G}$, the
innovation method is still outperformed by the marginal method for one-step-
ahead forecasting. In most cases, the improvement in these distances is slight
or negligible, and so since the innovation method has the smaller forecast
error, one could argue in favor of either method.

In the graphs of the Kullback-Leibler divergence for $\alpha=0.9$ and
$X_{n}=15$ or 30, an irregularity occurs. This can be explained by noting that
for $\alpha=0.9$, the value of $X_{n+1}$ would be expected to be similar to
that of $X_{n}$. However, as $\theta$ increases the mean would become smaller,
which contraindicates many values of $X_{n}$ being as large as 15. The erratic
behavior is even more pronounced when $X_{n}=30$. In our applications, in
Section 7, the largest estimated $\alpha$ is 0.3942. Therefore these
irregularities are noteworthy, but perhaps not too important in determining
which prediction method should be used for those data sets.

Finally, Figures 13 through 18 show $KL_{p,q}$ and $K_{F,G}$ for fixed values
of $\alpha=0.1,0.5$, and 0.9, and $\theta=0.5,2$, and 5, allowing
$X_{n}=0,1,2,\ldots$, up to 5 or 30, depending upon
${\mbox{E}}\left[X_{n}\right]$ as previously described. When $\alpha=0.1$ the
performances of the two methods are essentially the same as measured by the
two distances. Therefore the innovation method might be preferred as it has a
smaller one-step-ahead prediction error. When $\alpha=0.5$ and $\theta=2$ or
5, the marginal method has smaller distances, and so would be the preferred
method. When $\alpha=0.5$ and $\theta=0.5$ the two are close to equivalent.
Lastly, when $\alpha=0.9$ the marginal method outperforms the innovation
method expect when $\theta=0.5$, when the two are equivalent. In summary, for
combinations of small to moderate $\alpha$ and small $\theta$, the marginal
method appears to be more accurate in one-step-ahead prediction.

Figure 1: Kullback-Leibler divergences for the marginal method (black) and the
innovation method (red) when $\theta=0.5$. Figure 2: Kolmogorov Metric for the
marginal method (black) and the innovation method (red) when $\theta=0.5$.
Figure 3: Kullback-Leibler divergences for the marginal method (black) and the
innovation method (red) when $\theta=2$. Figure 4: Kolmogorov Metric for the
marginal method (black) and the innovation method (red) when $\theta=2$.
Figure 5: Kullback-Leibler divergences for the marginal method (black) and the
innovation method (red) when $\theta=5$. Figure 6: Kolmogorov Metric for the
marginal method (black) and the innovation method (red) when $\theta=5$.
Figure 7: Kullback-Leibler divergences for the marginal method (black) and the
innovation method (red) when $\alpha=0.1$. Figure 8: Kolmogorov Metric for the
marginal method (black) and the innovation method (red) when $\alpha=0.1$.
Figure 9: Kullback-Leibler divergences for the marginal method (black) and the
innovation method (red) when $\alpha=0.5$. Figure 10: Kolmogorov Metric for
the marginal method (black) and the innovation method (red) when $\alpha=0.5$.
Figure 11: Kullback-Leibler divergences for the marginal method (black) and
the innovation method (red) when $\alpha=0.9$. Figure 12: Kolmogorov Metric
for the marginal method (black) and the innovation method (red) when
$\alpha=0.9$. Figure 13: Kullback-Leibler divergences for the marginal method
(black) and the innovation method (red) when $\theta=0.5$. Figure 14:
Kolmogorov Metric for the marginal method (black) and the innovation method
(red) when $\theta=0.5$. Figure 15: Kullback-Leibler divergences for the
marginal method (black) and the innovation method (red) when $\theta=2$.
Figure 16: Kolmogorov Metric for the marginal method (black) and the
innovation method (red) when $\theta=2$. Figure 17: Kullback-Leibler
divergences for the marginal method (black) and the innovation method (red)
when $\theta=5$. Figure 18: Kolmogorov Metric for the marginal method (black)
and the innovation method (red) when $\theta=5$.

##  6 Estimation of Parameters for the PLINAR(1) Model

In what follows, we will use the forecasting algorithms described in Section 3
to find $k$-step ahead forecasts for a discrete time series. To accomplish
this, model parameters must be estimated. Therefore in our application
section, we investigate the effect of three different estimation techniques,
described in Section 6.1, on the forecasts using three different forecast
measures provided in Section 6.2.

###  6.1 Estimation Techniques

Three techniques are prevalent in the literature for estimating the parameters
of the PLINAR(1) model defined in Section 2; Yule-Walker, conditional least
squares, and maximum likelihood. We provide a brief overview of the three
techniques below. Mohammadpour et al., (2018) conducted a simulation study to
investigate and compare the perfomances of the three estimation methods and
concluded the maximum likelihood estimators are preferred.

Conditional least squares estimators are found by minimizing

$\sum^{n}_{t=2}(X_{t}-\alpha X_{t-1}-\mu(1-\alpha))^{2}$

, where $\mu=E(X_{t})=\theta+2/\theta(\theta+1)$. Thus, we have the estimates:

$\displaystyle\hat{\alpha}_{CLS}$
$\displaystyle=\frac{(n-1)\sum^{n}_{t=2}(X_{t-1}X_{t})-\sum^{n}_{t=2}X_{t}\sum^{n}_{t=2}X_{t-1}}{(n-1)\sum^{n}_{t=2}X^{2}_{t-1}-\left(\sum_{t=2}^{n}X_{t-1}\right)^{2}}$
$\displaystyle\hat{\mu}_{CLS}$
$\displaystyle=\frac{\sum^{n}_{t=2}X_{t}-\hat{\alpha}_{CLS}\sum^{n}_{t=2}X_{t-1}}{(n-1)(1-\hat{\alpha}_{CLS})}$
$\displaystyle\hat{\theta}_{CLS}$
$\displaystyle=\frac{-(\hat{\mu}_{CLS}-1)+\sqrt{(\hat{\mu}_{CLS}-1)^{2}+8\hat{\mu}_{CLS}}}{2\hat{\mu}_{CLS}}$

The Yule-Walker estimators of $\alpha$ and $\mu$ are

$\displaystyle\hat{\alpha}_{YW}$
$\displaystyle=\frac{\hat{\gamma}_{1}}{\hat{\gamma}_{0}}=\frac{\sum^{n}_{t=2}(X_{t}-\bar{X})(X_{t-1}-\bar{X})}{\sum_{t=1}^{n}\left(X_{t}-\bar{X}\right)^{2}}$
$\displaystyle\hat{\mu}_{YW}$
$\displaystyle=\bar{X}=\frac{\hat{\theta}+2}{\hat{\theta}(\hat{\theta}+1)}$
$\displaystyle\hat{\theta}_{YW}$
$\displaystyle=\frac{-(\bar{X}-1)+\sqrt{(\bar{X}-1)^{2}+8\bar{X}}}{2\bar{X}}$

Mohammadpour et al., (2018) prove that the asymptotic normal distributions for
both the conditional least squares estimators and the Yule-Walker estimators
are the same.

Finally, the maximum likelihood estimators of $\alpha$ and $\theta$ are the
values $\hat{\alpha}_{ML}$ and $\hat{\theta}_{ML}$ that maximize the log-
likelihood given by

$\displaystyle\ln
L(X_{1},\ldots,X_{n};\theta,\alpha)=2\log(\theta)-(X_{1}+3)\log(1+\theta)+\log(2+\theta+X_{1})+\sum_{t=2}^{n}\ln
p(X_{t}|X_{t-1};\theta,\alpha)$

###  6.2 Descriptive Measures of Forecasting Accuracy

To investigate the accuracy of the forecast in Section 7, we partition an
observed dataset of length $n$ into two sets. The longer set, called the
training set, contains the first $m$ observations, and is used to estimate the
parameters of the models. The remaining $n-m$ observations are used for
evaluating forecasts using each of the following three criteria. Our design is
similar to Maiti et al., (2016). They consider a variety of models, but do not
match moments, as does our marginal and innovation methods. Their comparisons
only consider PRMSE and PMAD, but not PTP, which we feel provides a clearer
picture of which models perform better. In their simulated studies, they
conclude that the discretized Gaussian forecasts are comparable to or out-
perform some discrete models. In their work, as in ours, the idea of coherency
is central. A coherent forecast is a forecast that is in the support of the
process. When the support is restricted, as is the case with integer
autoregressive processes, the conditional mean may produce forecasts that are
not coherent. In that case, the noncoherent value will be rounded in the
traditional manner to produce a coherent forecast.

The first of the three criteria is the $k$-step-ahead prediction root mean
square error, PRMSE($k$), defined as

$\displaystyle{\text{PRMSE}}(k)$
$\displaystyle=\sqrt{\frac{1}{n-m-k+1}\sum_{t=m+k}^{n}\left(X_{t}-\hat{X}^{me}_{t}\right)^{2}},$

where $\hat{X}^{me}_{t+k}$ is estimated $k$-step-ahead conditional mean given
in (1).

The second forecasting criteria we consider is the $k$-step-ahead prediction
mean absolute deviation, PMAD($k$), given by

$\displaystyle{\text{PMAD}}(k)$
$\displaystyle=\frac{1}{n-m-k+1}\sum_{t=m+k}^{n}\left|X_{t}-\hat{X}^{med}_{t}\right|,$

where $\hat{X}^{med}_{t+k}$ is the median of the estimated $k$-step-ahead
conditional distribution in (3).

The third and final criteria used to investigate the performances of the
forecasts is the $k$-step-ahead percentage of true prediction, PTP($k$), given
by

$\displaystyle{\text{PTP}}(k)$
$\displaystyle=\frac{\sum^{n}_{t=m+k}I\left(X_{t}=\hat{X}_{t}\right)}{n-m-k+1}\times
100\%,$

where $I(\cdot)$ is an indicator function and $\hat{X}_{t}$ is one of the
$k$-step-ahead conditional mean rounded to the nearest integer, the median, or
the mode of the estimated $k$-step-ahead conditional distribution given in
(3).

##  7 Application to Discrete Valued Time Series Data

In this section we consider the data given in Table 1. Each record is the
number of sex offences reported in the 21st police car beat in Pittsburgh for
one month, beginning in January 1990 through December 2001 (Ristić et al.,,
2009). The data is publicly available in the _R_ library tsinteger under the
dataset name sexoffences. In what follows, we provide a description of the
data’s general features. Following that, we determine an optimal
autoregressive order using Akaike information criteria (AIC), corrected AIC
(AICc) and the Bayesian information criteria (BIC). We then fit both a PLINAR
and Gaussian AR to the series using the suggested optimal AR order. The models
are then used to forecast the data using the methods described in Section 4.
Because the parameters of the models must be estimated, forecasts are also
compared for three estimation techniques; namely, conditional least squares,
Yule-Walker, and maximum likelihood.

###  7.1 Data Description

The 144 observations in Table 1 are counts of sex offences reported. The data
has a high percentage (62.5%) of observations equal to zero. The values in the
data range from zero to six.

| Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec
---|---|---|---|---|---|---|---|---|---|---|---|---
1990 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0
1991 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 0
1992 | 0 | 0 | 0 | 0 | 1 | 1 | 2 | 1 | 0 | 1 | 0 | 0
1993 | 1 | 2 | 0 | 0 | 0 | 0 | 1 | 0 | 2 | 0 | 0 | 0
1994 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 1 | 0 | 3 | 1 | 0
1995 | 1 | 1 | 1 | 0 | 3 | 1 | 0 | 0 | 1 | 2 | 2 | 0
1996 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0
1997 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0
1998 | 0 | 0 | 0 | 0 | 0 | 1 | 2 | 2 | 0 | 2 | 0 | 0
1999 | 1 | 1 | 0 | 3 | 2 | 0 | 0 | 2 | 0 | 0 | 0 | 0
2000 | 1 | 1 | 6 | 5 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 0
2001 | 0 | 1 | 1 | 0 | 1 | 0 | 1 | 5 | 0 | 0 | 0 | 0
Table 1: Number of sex offences reported in the 21st police car beat in Pittsburgh monthly from January 1990 to December 2001. Figure 19: Frequency histogram of the number of sex offences reported in the 21st police car beat in Pittsburgh monthly from January 1990 to December 2001. Minimum | Maximum | Median | Mode | Mean | Variance
---|---|---|---|---|---
0 | 6 | 0 | 0 | 0.5903 | 1.0268
Table 2: Summary statistics for number of sex offences reported in the 21st
police car beat in Pittsburgh monthly from January 1990 to December 2001.

The histogram in Figure 19 reveals a right-skewed distribution. The value of
zero is the mode of the distribution. For the 144 observed records, 90 (62.5%)
have a value of zero; 36 (25.00%) have a value of one; 12 (8.33%) a value of
two, and 3 (2.08%) a value of three. No records have a value of four, and only
three (2.08%) have a value of five or six. Table 2 contains the summary
statistics for the data. The median of the data is zero. The mean number of
sex offences is $\bar{X}=0.5903$, and the variance is $s_{X}^{2}=1.0268$. The
index of dispersion is $\hat{I}=s^{2}_{X}/\bar{X}=1.7395$, which is large
enough to conclude that the data is over-dispersed, which provides a good
reason for considering Poisson-Lindley marginals.

Figure 21 contains the time plot of the number of sex offences reported in the
21st police car beat in Pittsburgh monthly from January 1990 to December 2001.
The data record (Table 1) and the time plot both show that the number of sex
offences in the winter months (December - March) is much less than during the
remainder of the year. The month-to-month values tend to be close to the same.
In March of 2000 there were six sex offences and in April of 2000 and August
of 2001 there were five, both of which are unusually high values. Figure 21
contains the sample autocorrelation function (left) and sample partial
autocorrelation function (right) with large sample 95% confidence bounds for
white noise superimposed. The sample ACF indicates lag one is significant. The
sample PACF has only lag one as significant as well. The indication is that an
AR(1) would be a good fit to the data. To further aid in determining the
optimal AR order for fitting a model, we compute AIC, AICc, and BIC. Those
values are supplied in Table 3. All three measures indicate that the optimal
order for the AR model is one.

| AIC | AICc | BIC
---|---|---|---
AR(1) | 409.3019 | 409.4734 | 418.2114
AR(2) | 411.2420 | 411.5298 | 423.1212
AR(3) | 413.2324 | 413.6671 | 428.0814
Table 3: Autoregressive order comparison for number of sex offences reported
in the 21st police car beat in Pittsburgh monthly from January 1990 to
December 2001.

Figure 20: Time plot of number of sex offences reported in the 21st police car
beat in Pittsburgh monthly from January 1990 to December 2001.

Figure 21: Sample autocorrelation function (left) and sample partial
autocorrelation function (right) for the sex offences data with large sample
95% confidence bounds for white noise superimposed (blue dashed lines).

Table 4 provides estimates for $\alpha$ and $\theta$ using the conditional
least squares (CLS), Yule-Walker (YW) and maximum likelihood (ML) when fitting
the PLINAR(1) model. We note that the ML estimates for $\alpha$ is smaller
than the CLS and YW estimates. The small value of the MLE can be attributed to
the large number of zeros in the data.

| CLS | YW | ML
---|---|---|---
$\alpha$ | 0.2297 | 0.2291 | 0.1028
$\theta$ | 2.1671 | 2.1804 | 2.1900
Table 4: Conditional least squares, Yule-Walker, and maximum likelihood
estimates of $\alpha$ and $\theta$ for the sex offences data.

Tables 5, 6, and 7 provide the conditional $k$-step-ahead probability mass
functions given $X_{n}=0$ where model parameters are estimated using each of
CLS, YW, and ML based on the traditional PLINAR(1) model, the model using the
marginal means method, and the model using the innovation method,
respectively. To find the estimates of the PMFs, the first 141 values of the
data were used. In this study, we consider only $X_{n}=0$ since all four of
observation numbers 141 through 144 are zero. In subsequent application,
conditional $k$-step-ahead forecasts for values of $X_{n}$ other than zero are
considered. The conditional PMFs are used for computing the means, medians,
and modes that are used for the $k$-step-ahead forecasts in Table 8 for the
PLINAR(1), marginal, and innovation methods, respectively. As well as the
traditional Gaussian AR(1) was considered. We note that Table 5 provides the
conditional PMFs for only up to $X_{n+k}=8$, and using these values one is
able to compute the means found in Table 8.

| CLS | YW | MLE
---|---|---|---
| $k=1$ | $k=2$ | $k=3$ | $k=\infty$ | $k=1$ | $k=2$ | $k=3$ | $k=\infty$ | $k=1$ | $k=2$ | $k=3$ | $k=\infty$
$p_{k}(0|0)$ | 0.703 | 0.636 | 0.621 | 0.616 | 0.704 | 0.637 | 0.622 | 0.618 | 0.657 | 0.623 | 0.619 | 0.619
$p_{k}(1|0)$ | 0.188 | 0.229 | 0.238 | 0.241 | 0.188 | 0.229 | 0.238 | 0.241 | 0.217 | 0.238 | 0.240 | 0.240
$p_{k}(2|0)$ | 0.070 | 0.086 | 0.090 | 0.091 | 0.070 | 0.086 | 0.089 | 0.090 | 0.081 | 0.089 | 0.090 | 0.090
$p_{k}(3|0)$ | 0.025 | 0.032 | 0.033 | 0.033 | 0.025 | 0.031 | 0.033 | 0.033 | 0.029 | 0.032 | 0.033 | 0.033
$p_{k}(4|0)$ | 0.009 | 0.011 | 0.012 | 0.012 | 0.009 | 0.011 | 0.012 | 0.012 | 0.010 | 0.012 | 0.012 | 0.012
$p_{k}(5|0)$ | 0.003 | 0.004 | 0.004 | 0.004 | 0.003 | 0.004 | 0.004 | 0.004 | 0.004 | 0.004 | 0.004 | 0.004
$p_{k}(6|0)$ | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001
$p_{k}(7|0)$ | 0.000 | 0.000 | 0.001 | 0.001 | 0.000 | 0.000 | 0.000 | 0.001 | 0.000 | 0.000 | 0.000 | 0.000
$p_{k}(8|0)$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
Table 5: Estimated conditional $k$-step-ahead probability mass function for
$X_{n}=0$ for three estimation methods for $\alpha$ and $\theta$ based on
PLINAR(1) model for the sex offences data.

| CLS | YW | ML
---|---|---|---
| $k=1$ | $k=2$ | $k=3$ | $k=\infty$ | $k=1$ | $k=2$ | $k=3$ | $k=\infty$ | $k=1$ | $k=2$ | $k=3$ | $k=\infty$
$p_{k}(0|0)$ | 0.309 | 0.276 | 0.267 | 0.265 | 0.310 | 0.276 | 0.268 | 0.265 | 0.286 | 0.268 | 0.266 | 0.266
$p_{k}(1|0)$ | 0.405 | 0.395 | 0.393 | 0.393 | 0.406 | 0.396 | 0.395 | 0.395 | 0.400 | 0.396 | 0.396 | 0.396
$p_{k}(2|0)$ | 0.234 | 0.260 | 0.266 | 0.267 | 0.233 | 0.259 | 0.265 | 0.267 | 0.251 | 0.265 | 0.266 | 0.266
$p_{k}(3|0)$ | 0.048 | 0.064 | 0.067 | 0.068 | 0.047 | 0.063 | 0.066 | 0.067 | 0.058 | 0.065 | 0.066 | 0.066
$p_{k}(4|0)$ | 0.004 | 0.006 | 0.006 | 0.006 | 0.003 | 0.005 | 0.006 | 0.006 | 0.005 | 0.006 | 0.006 | 0.006
$p_{k}(5|0)$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
$p_{k}(6|0)$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
$p_{k}(7|0)$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
$p_{k}(8|0)$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
Table 6: Estimated conditional $k$-step-ahead probability mass function fo
$X_{n}=0$ for three estimation methods for $\alpha$ and $\theta$ based on
marginal method for the sex offences data.

| CLS | YW | ML
---|---|---|---
| $k=1$ | $k=2$ | $k=3$ | $k=\infty$ | $k=1$ | $k=2$ | $k=3$ | $k=\infty$ | $k=1$ | $k=2$ | $k=3$ | $k=\infty$
$p_{k}(0|0)$ | 0.298 | 0.262 | 0.254 | 0.251 | 0.298 | 0.263 | 0.254 | 0.252 | 0.280 | 0.261 | 0.259 | 0.259
$p_{k}(1|0)$ | 0.429 | 0.418 | 0.417 | 0.416 | 0.431 | 0.420 | 0.419 | 0.418 | 0.411 | 0.408 | 0.408 | 0.407
$p_{k}(2|0)$ | 0.232 | 0.262 | 0.268 | 0.270 | 0.231 | 0.261 | 0.267 | 0.269 | 0.252 | 0.266 | 0.267 | 0.267
$p_{k}(3|0)$ | 0.039 | 0.054 | 0.057 | 0.058 | 0.038 | 0.053 | 0.056 | 0.057 | 0.053 | 0.060 | 0.061 | 0.061
$p_{k}(4|0)$ | 0.002 | 0.004 | 0.004 | 0.004 | 0.002 | 0.003 | 0.004 | 0.004 | 0.004 | 0.005 | 0.005 | 0.005
$p_{k}(5|0)$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
$p_{k}(6|0)$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
$p_{k}(7|0)$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
$p_{k}(8|0)$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
Table 7: Estimated conditional $k$-step-ahead probability mass function fo
$X_{n}=0$ for three estimation methods for $\alpha$ and $\theta$ based on
innovation method for the sex offences data.

|  | CLS | YW | ML
---|---|---|---|---
|  | $k=1$ | $k=2$ | $k=3$ | $k=1$ | $k=2$ | $k=3$ | $k=1$ | $k=2$ | $k=3$
PLINAR(1) | Mean | 0.468 | 0.575 | 0.600 | 0.465 | 0.571 | 0.596 | 0.538 | 0.593 | 0.599
Med. | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
Mode | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
Marginal Method | Mean | 0.468 | 0.575 | 0.600 | 0.465 | 0.571 | 0.596 | 0.538 | 0.593 | 0.599
Med. | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Mode | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Innovation Method | Mean | 0.468 | 0.575 | 0.600 | 0.465 | 0.571 | 0.596 | 0.538 | 0.593 | 0.599
Med. | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Mode | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Gauss. |  |  |  |  |  |  |  |  |  |
AR(1) | 0.468 | 0.575 | 0.600 | 0.465 | 0.571 | 0.596 | 0.463 | 0.569 | 0.593
| Truth | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Table 8: Conditional $k$-step-ahead forecasts for the sex offences data for
each of the PLINAR(1) model, the marginal method, innovation method. For each
of the three discrete model methods, model parameters were estimated using the
three estimation methods in Section 6.1. For the Gaussian AR(1) and Gaussian
AR(2), model parameters were estimated using ordinary least squares, Yule-
Walker, or maximum likelihood via the _R_ function ar.ols, ar.yw, and ar.mle.

For this data $X_{141}=X_{142}=X_{143}=X_{144}=0$. From Table 8, the median
and mode of the PLINAR(1) all provide correct forecasts for $k=1,2$ and 3. The
rounded mean for all four models provides a correct forecast for one-step-
ahead. However, for $k=2$ and 3, the means incorrectly forecasts a value of
one. For the marginal and innovation methods, the median and modes incorrectly
forecast a value of one for all $k=1,2$, and 3.

|  | PRMSE | PMAD | PTP (Med.) | PTP (Mode)
---|---|---|---|---|---
|  | $k=1$ | $k=2$ | $k=3$ | $k=1$ | $k=2$ | $k=3$ | $k=1$ | $k=2$ | $k=3$ | $k=1$ | $k=2$ | $k=3$
CLS | PLINAR(1) | 1.597 | 1.669 | 1.700 | 0.931 | 0.929 | 0.963 | 0.517 | 0.536 | 0.519 | 0.517 | 0.536 | 0.519
Marginal | 1.597 | 1.669 | 1.700 | 0.966 | 1.000 | 1.000 | 0.345 | 0.357 | 0.370 | 0.345 | 0.357 | 0.370
Innovation | 1.597 | 1.669 | 1.700 | 0.966 | 1.000 | 1.000 | 0.345 | 0.357 | 0.370 | 0.345 | 0.357 | 0.370
Traditional | 1.597 | 1.669 | 1.700 | 0.966 | 0.929 | 0.889 | 0.414 | 0.464 | 0.519 | 0.414 | 0.464 | 0.519
YW | PLINAR(1) | 1.597 | 1.669 | 1.895 | 0.931 | 0.929 | 0.963 | 0.517 | 0.536 | 0.519 | 0.517 | 0.536 | 0.519
Marginal | 1.597 | 1.669 | 1.895 | 0.966 | 1.000 | 1.000 | 0.345 | 0.357 | 0.370 | 0.345 | 0.357 | 0.370
Innovation | 1.597 | 1.669 | 1.895 | 0.966 | 1.000 | 1.000 | 0.345 | 0.357 | 0.370 | 0.345 | 0.357 | 0.370
Traditional | 1.597 | 1.669 | 1.895 | 0.966 | 0.929 | 1.000 | 0.414 | 0.464 | 0.481 | 0.414 | 0.464 | 0.481
ML | PLINAR(1) | 1.597 | 1.852 | 1.886 | 0.931 | 0.929 | 0.963 | 0.517 | 0.536 | 0.519 | 0.517 | 0.536 | 0.519
Marginal | 1.597 | 1.852 | 1.886 | 0.966 | 1.000 | 1.000 | 0.345 | 0.357 | 0.370 | 0.345 | 0.357 | 0.370
Innovation | 1.597 | 1.852 | 1.886 | 0.966 | 1.000 | 1.000 | 0.345 | 0.357 | 0.370 | 0.345 | 0.357 | 0.370
Traditional | 1.597 | 1.669 | 1.895 | 0.966 | 0.929 | 1.000 | 0.414 | 0.464 | 0.481 | 0.414 | 0.464 | 0.481
Table 9: Accuracy measures of forecasts of sex offences data based on a
PLINAR(1) model, the marginal method, innovation method, and traditional
Gaussian AR(1) where model parameters are estimated using conditional least
squares (CLS), Yule-Walker (YW) and maximum likelihood (ML).

The measures of forecasts accuracy in Table 9 were computed training the
models on the first 80% (115 observations) of the data. From that training
set, $k$-step-ahead forecasts were computed for the remaining 20% (29
observations) for $k=1,2$, and 3. Focusing on PRMSE and PMAD, the forecasts
using CLS, YW and ML estimates for the models perform equally well. Focusing
on the PTP, the median and mode performed equally well for all four processes.

##  8 Conclusions

In this paper, we investigated forecasting discrete time series using a
PLINAR(1) model, a traditional Gaussian AR(1) model, and two modified version
of the traditional Gaussian AR(1) model. The resultant distribution functions
were compared via Kullback-Liebler divergence and Kolmogorov metric for a
variety of values of $\alpha$ and $\theta$. To determine which of the marginal
method or innovation method performs best in approximating the conditional
Poisson-Lindley PMF, the innovation method would be preferred when it is
believed that values of $\alpha$ and $\theta$ are close to zero or one, since
the innovation method has a smaller standard error compared to the marginal
method, and the two approximations perform equally well. However, for moderate
values of $\alpha$ and moderate to large values of $\theta$, the distances
between the marginal approximation and the conditional Poisson-Lindley PMF are
smaller, and so the marginal method would be preferred.

One application was considered to further evaluate and compare the forecasting
potential of the PLINAR(1), and the marginal, innovation, and traditional
Gaussian AR for discrete-valued time series. This data set had overdispersion,
which could provide justification for using the PLINAR. For this data set, the
PLINAR(1) model outperformed its Gaussian counterparts. It should also be
noted that the zero-inflated Poisson distribution may be a better choice of
distribution for modeling this data, and that work is being considered in
additional research.

The authors would like to thank Annika Homburg for helpful discussions and
allowing us to use her code based on the Poisson distribution to help develop
our own for the Poisson-Lindley.

## References

  * Al-Osh and Alzaid, (1987) Al-Osh, M. A. and Alzaid, A. A. (1987).  First-order integer-valued autoregressive (INAR(1)) process.  Journal of Time Series Analysis, 8(3):261–275.
  * Alwan and Weiß, (2017) Alwan, L. C. and Weiß, C. H. (2017).  INAR implementation of newsvendor model for serially dependent demand counts.  International Journal of Production Research, 55(4):1085–1099.
  * Alzaid and Al-Osh, (1988) Alzaid, A. A. and Al-Osh, M. A. (1988).  First-order integer-valued autoregressive (INAR (1)) process: Distributional and regression properties.  Statistica Neerlandica, 42(1):53–61.
  * Alzaid and Al-Osh, (1993) Alzaid, A. A. and Al-Osh, M. A. (1993).  Some autoregressive moving average processes with generalized poisson marginal distributions.  Annals of the Institute of Statistical Mathematics, 45(2):223–232.
  * Bisaglia and Gerolimetto, (2015) Bisaglia, L. and Gerolimetto, M. (2015).  Forecasting integer autoregressive processes of order 1: Are simple ar competitive?  Economics Bulletin, 35(3):1652–1660.
  * Box and Jenkins, (1976) Box, G. and Jenkins, G. M. (1976).  Time Series Analysis: Forecasting and Control.  Holden-Day.
  * Brockwell and Davis, (1996) Brockwell, P. J. and Davis, R. A. (1996).  Introduction to Time Series and Forecasting.  Springer Texts in Statistics. Springer, New York, NY.
  * Ghitany and Al-Mutairi, (2009) Ghitany, M. E. and Al-Mutairi, D. K. (2009).  Estimation methods for the discrete Poisson–Lindley distribution.  Journal of Statistical Computation and Simulation, 79(1):1–9.
  * Gidaro and Harvill, (2024) Gidaro, R. and Harvill, J. L. (2024).  Evaluating approximations of count distributions and forecasts for geometric autoregressive processes.  Communications in Statistics: Theory and Methods.
  * Homburg, (2020) Homburg, A. (2020).  Criteria for evaluating approximations of count distributions.  Communications in Statistics – Simulation and Computation, 49(12):3152–3170.
  * (11) Jacobs, P. A. and Lewis, P. A. W. (1978a).  Discrete time series generated by mixtures. i: Correlational and runs properties.  Journal of the Royal Statistical Society. Series B (Methodological), 40(1):94–105.
  * (12) Jacobs, P. A. and Lewis, P. A. W. (1978b).  Discrete time series generated by mixtures ii: Asymptotic properties.  Journal of the Royal Statistical Society. Series B (Methodological), 40(2):222–228.
  * Jacobs and Lewis, (1983) Jacobs, P. A. and Lewis, P. A. W. (1983).  Stationary discrete autoregressive-moving average time series generated by mixtures.  Journal of Time Series Analysis, 4(1):19–36.
  * Maiti et al., (2016) Maiti, R., Biswas, A., and Das, S. (2016).  Coherent forecasting for count time series using Box–Jenkins’s AR(p) model.  Statistica Neerlandica, 70(2):123–145.
  * McKenzie, (1986) McKenzie, E. (1986).  Autoregressive moving-average processes with negative-binomial and geometric marginal distributions.  Advances in Applied Probability, 18(3):679–705.
  * McKenzie, (1988) McKenzie, E. (1988).  Some ARMA models for dependent sequences of Poisson counts.  Advances in Applied Probability, 20(4):822–835.
  * Mohammadpour et al., (2018) Mohammadpour, M., Bakouch, H. S., and Shirozhan, M. (2018).  Poisson–lindley INAR(1) model with applications.  Brazilian Journal of Probability and Statistics, 32(2):262–280.
  * Ristić et al., (2009) Ristić, M. M., Bakouch, H. S., and Nastić, A. S. (2009).  A new geometric first-order integer-valued autoregressive (NGINAR(1)) process.  Journal of Statistical Planning and Inference, 139(7):2218–2226.
  * Sankaran, (1970) Sankaran, M. (1970).  275\. note: The discrete Poisson-Lindley distribution.  Biometrics, 26(1):145–149.
  * Shumway and Stoffer, (2017) Shumway, R. H. and Stoffer, D. S. (2017).  Time Series Analysis and Its Applications: With R Examples.  Springer Texts in Statistics. Springer International Publishing, Cham.
  * Steutel and van Harn, (1979) Steutel, F. W. and van Harn, K. (1979).  Discrete Analogues of Self-Decomposability and Stability.  The Annals of Probability, 7(5):893 – 899.


# Complete and incomplete Bell polynomials associated with Lah-Bell numbers
and polynomials

Taekyun Kim 1 1 Department of Mathematics, Kwangwoon University, Seoul
139-701, Republic of Korea<EMAIL_ADDRESS><EMAIL_ADDRESS>,  DAE
SAN KIM 2 2 Department of Mathematics, Sogang University, Seoul 121-742,
Republic of Korea<EMAIL_ADDRESS>,  Lee-Chae Jang3 3 Graduate School of
Education, Konkuk University, Seoul 05029, Republic of Korea
<EMAIL_ADDRESS>,  Hyunseok Lee 4 4 Department of Mathematics, Kwangwoon
University, Seoul 139-701, Republic of Korea<EMAIL_ADDRESS>and
Han-Young Kim 5 5 Department of Mathematics, Kwangwoon University, Seoul
139-701, Republic of Korea<EMAIL_ADDRESS>
###### Abstract.

The $n$th $r$-extended Lah-Bell number is defined as the number of ways a set
with $n+r$ elements can be partitioned into ordered blocks such that $r$
distinguished elements have to be in distinct ordered blocks. The aim of this
paper is to introduce incomplete $r$-extended Lah-Bell polynomials and
complete $r$-extended Lah-Bell polynomials respectively as multivariate
versions of $r$-Lah numbers and the $r$-extended Lah-Bell numbers and to
investigate some properties and identities for these polynomials. From these
investigations, we obtain some expressions for the $r$-Lah numbers and the
$r$-extended Lah-Bell numbers as finite sums.

###### Key words and phrases:

incomplete $r$-extended Lah-Bell polynomial; complete $r$-extended Lah-Bell
polynomial

###### 2010 Mathematics Subject Classification:

11B73; 11B83; 05A19

##  1\. Introduction

It is well known that the unsigned Lah number $L(n,k),\ (n\geq k\geq 0)$,
counts the number of ways a set with $n$ elements can be partitioned into $k$
non-empty linearly ordered subsets (see [4,7,8]). The $n$-th Lah-Bell number
$B_{n}^{L},\ (n\geq 0)$, is the number of ways a set with $n$ elements can be
partitioned into non-empty linearly ordered subsets. Thus, we note that

(1) $B_{n}^{L}=\sum_{k=0}^{n}L(n,k),\quad(n\geq 0),\quad(\mathrm{see}\
[7,8]).$

From (1), it can be shown that the generating function of Lah-Bell numbers is
given by

(2)
$e^{\frac{t}{1-t}}=\sum_{n=0}^{\infty}B_{n}^{L}\frac{t^{n}}{n!},\quad(\mathrm{see}\
[7,8]),$

where

(3)
$\frac{1}{k!}\bigg{(}\frac{t}{1-t}\bigg{)}^{k}=\sum_{n=k}^{\infty}L(n,k)\frac{t^{n}}{n!},\
(k\geq 0),\quad(\mathrm{see}\ [8,13,15,17]).$

Explicitly, we see from (3) that the Lah numbers are given by

$L(n,k)=\frac{n!}{k!}\binom{n-1}{k-1},\quad(n\geq k\geq 0),\quad(\mathrm{see}\
[7,8,9,10,17,19]).$

Let $n,k,r$ be non-negative integers with $n\geq k$. Then, the $r$-Lah number
$L_{r}(n,k)$ counts the number of partitions of a set with $n+r$ elements into
$k+r$ ordered blocks such that $r$ distinguished elements have to be in
distinct ordered blocks (see [17]). The $r$-extended Lah-Bell number
$B_{n,r}^{L}$ is defined as the number of ways a set with $n+r$ elements can
be partitioned into ordered blocks such that $r$ distinguished elements have
to be in distinct ordered blocks (see [7]). By the definitions of $r$-Lah
numbers and $r$-extended Lah-Bell numbers, we have

(4) $B_{n,r}^{L}=\sum_{k=0}^{n}L_{r}(n,k),\quad(n\geq 0),\quad(\mathrm{see}\
[7]).$

From (4), we see that the generating function of $r$-extended Lah-Bell numbers
is given by

(5)
$e^{\frac{t}{1-t}}\bigg{(}\frac{1}{1-t}\bigg{)}^{2r}=\sum_{n=0}^{\infty}B_{n,r}^{L}\frac{t^{n}}{n!},\quad(\mathrm{see}\
[7,15]),$

where

(6)
$\frac{1}{k!}\bigg{(}\frac{t}{1-t}\bigg{)}^{k}\bigg{(}\frac{1}{1-t}\bigg{)}^{2r}=\sum_{n=k}^{\infty}L_{r}(n,k)\frac{t^{n}}{n!},\quad(\mathrm{see}\
[7,17]),$

where $k$ is a non-negative integer.

Explicitly, the $r$-Lah numbers are given by

$L_{r}(n,k)=\frac{n!}{k!}\binom{n+2r-1}{k+2r-1},\quad(n\geq k\geq
0),\quad(\mathrm{see}\ [7,8,9,10,17,19]).$

In [7], the $r$-extended Lah-Bell polynomials are defined by

(7)
$e^{x(\frac{t}{1-t})}\bigg{(}\frac{1}{1-t}\bigg{)}^{2r}=\sum_{n=0}^{\infty}B_{n,r}^{L}(x)\frac{t^{n}}{n!}.$

It is well known that the complete Bell polynomials are defined by

(8)
$\exp\bigg{(}\sum_{j=1}^{\infty}x_{j}\frac{t^{j}}{j!}\bigg{)}=\sum_{n=0}^{\infty}B_{n}(x_{1},x_{2},\dots,x_{n})\frac{t^{n}}{n!},\quad(\mathrm{see}\
[2,3,4,5,11,14]).$

Then it can be shown that the complete Bell polynomials are given by

(9)
$B_{n}(x_{1},x_{2},\dots,x_{n})=\sum_{j_{1}+2j_{2}+\cdots+nj_{n}=n}\frac{n!}{j_{1}!j_{2}!\cdots
j_{n}!}\bigg{(}\frac{x_{1}}{1}\bigg{)}^{j_{1}}\bigg{(}\frac{x_{2}}{2!}\bigg{)}^{j_{2}}\cdots\bigg{(}\frac{x_{n}}{n!}\bigg{)}^{j_{n}},$

where the sum runs over all nonnegative integers $j_{1},j_{2},\dots,j_{n}$
satisfying $j_{1}+2j_{2}+\cdots+nj_{n}=n$. The incomplete Bell polynomials are
given by

(10)
$\frac{1}{k!}\bigg{(}\sum_{m=1}^{\infty}x_{m}\frac{t^{m}}{m!}\bigg{)}^{k}=\sum_{n=k}^{\infty}B_{n,k}(x_{1},x_{2},\dots,x_{n-k+1})\frac{t^{n}}{n!},\quad(n\geq
0),\quad(\mathrm{see}\ [5,11,14]).$

Thus, we note that

(11)
$B_{n,k}(x_{1},x_{2},\dots,x_{n-k+1})=\sum_{\pi(n,k)}\frac{n!}{j_{1}!j_{2}!\cdots
j_{n-k+1}!}\bigg{(}\frac{x_{1}}{1!}\bigg{)}^{j_{1}}\bigg{(}\frac{x_{2}}{2!}\bigg{)}^{j_{2}}\cdots\bigg{(}\frac{x_{n-k+1}}{(n-k+1)!}\bigg{)}^{j_{n-k+1}},$

where the sum runs over the set $\pi(n,k)$ of all nonnegative integers
$(j_{i})_{i\geq 1}$ satisfying $j_{1}+j_{2}+\cdots+j_{n-k+1}=k$, and
$1j_{1}+2j_{2}+\cdots+(n-k+1)j_{n-k+1}=n$.

Then the complete and incomplete Bell polynomials are related by

$B_{n}(x_{1},x_{2},\dots,x_{n})=\sum_{k=1}^{n}B_{n,k}(x_{1},x_{2},\dots,x_{n-k+1}),\quad(n\geq
1).$

Let $f$ be a $C^{\infty}$-function. Then we have

(12) $\displaystyle e^{f(x+t)}\ $ $\displaystyle=\
\exp\bigg{(}\sum_{j=0}^{\infty}f^{(j)}(x)\frac{t^{j}}{j!}\bigg{)}$
$\displaystyle=\
\exp\bigg{(}f(x)+\sum_{j=1}^{\infty}f^{(j)}(x)\frac{t^{j}}{j!}\bigg{)}$
$\displaystyle\ =\
e^{f(x)}\bigg{(}1+\sum_{n=1}^{\infty}B_{n}\big{(}f^{(1)}(x),f^{(2)}(x),\dots,f^{(n)}(x)\big{)}\frac{t^{n}}{n!}\bigg{)},$

where $f^{(j)}(x)$ is the $j$-th derivative of $f(x)$ and $\exp(t)=e^{t}$.

We observe that

(13) $\displaystyle\frac{d^{m}}{dx^{m}}e^{f(x)}\ $ $\displaystyle=\
\frac{\partial^{m}}{\partial x^{m}}e^{f(x+t)}\bigg{|}_{t=0}\ =\
\frac{\partial^{m}}{\partial t^{m}}e^{f(x+t)}\bigg{|}_{t=0}$ $\displaystyle=\
e^{f(x)}B_{m}\big{(}f^{(1)}(x),f^{(2)}(x),\dots,f^{(m)}(x)\big{)}.$

From (12) and (13), we obtain Kölbig and Coeffey’s equation as follows :

(14)
$\frac{d^{m}}{dx^{m}}e^{f(x)}=e^{f(x)}B_{m}\big{(}f^{(1)}(x),f^{(2)}(x),\dots,f^{(m)}(x)\big{)},\quad(m\geq
1),\quad(\mathrm{see}\ [6,12]).$

The exponential incomplete $r$-Bell polynomials are defined by the generating
function

(15)
$\frac{1}{k!}\bigg{(}\sum_{j=1}^{\infty}a_{j}\frac{t^{j}}{j!}\bigg{)}^{k}\bigg{(}\sum_{i=0}^{\infty}b_{i+1}\frac{t^{i}}{i!}\bigg{)}^{r}=\sum_{n=k}^{\infty}B_{n+r,k+r}^{(r)}(a_{1},a_{2},\dots:b_{1},b_{2},\dots)\frac{t^{n}}{n!}.$

From (15), we note that

(16) $\displaystyle\ B_{n+r,k+r}^{(r)}(a_{1},a_{2},\dots:b_{1},b_{2},\dots)$
$\displaystyle\quad=\
\sum_{\Lambda(n,k,r)}\bigg{[}\frac{n!}{k_{1}!k_{2}!k_{3}\cdots}\bigg{(}\frac{a_{1}}{1!}\bigg{)}^{k_{1}}\bigg{(}\frac{a_{2}}{2!}\bigg{)}^{k_{2}}\bigg{(}\frac{a_{3}}{3!}\bigg{)}^{k_{3}}\cdots\bigg{]}$
$\displaystyle\qquad\quad\times\bigg{[}\frac{r!}{r_{0}!r_{1}!r_{2}\cdots}\bigg{(}\frac{b_{1}}{0!}\bigg{)}^{r_{0}}\bigg{(}\frac{b_{2}}{1!}\bigg{)}^{r_{1}}\bigg{(}\frac{b_{3}}{2!}\bigg{)}^{r_{2}}\cdots\bigg{]},$

where $\Lambda(n,k,r)$ denotes the set of all non-negative integers
$(k_{i})_{i\geq 1}$ and $(r_{i})_{i\geq 0}$ such that

$\sum_{i\geq 1}k_{i}=k,\quad\sum_{i\geq
0}r_{i}=r\quad\mathrm{and}\quad\sum_{i\geq
1}i(k_{i}+r_{i})=n,\quad(\mathrm{see}\ [4,5,14]).$

Assume that $(a_{i})_{i\geq 1}$ and $(b_{i})_{i\geq 1}$ are sequences of
positive integers.

Then the number $B_{n+r,k+r}^{(r)}(a_{1},a_{2},\dots;b_{1},b_{2},\dots)$
counts the number of partitions of an $(n+r)$-set into $(k+r)$ blocks
satisfying:

  * •

The first $r$ elements belong to different blocks,

  * •

Any block of size $i$ containing no elements from the first $r$ elements can
be colored with $a_{i}$ colors,

  * •

Any block of size $i$ containing one element from the first $r$ elements can
be colored with $b_{i}$ colors.

The complete $r$-Bell polynomials are given by

(17)
$\exp\bigg{(}\sum_{i=1}^{\infty}a_{i}\frac{t^{i}}{i!}\bigg{)}\bigg{(}\sum_{j=0}^{\infty}b_{j+1}\frac{t^{j}}{j!}\bigg{)}^{r}=\sum_{n=0}^{\infty}B_{n}^{(r)}(a_{1},a_{2},\dots;b_{1},b_{2},\dots)\frac{t^{n}}{n!},$

(see [4,5,11,14]).

By (16) and (17), we get

(18)
$B_{n}^{(r)}(a_{1},a_{2},\dots:b_{1},b_{2},\dots)=\sum_{k=0}^{n}B_{n+r,k+r}^{(r)}(a_{1},a_{2},\dots:b_{1},b_{2},\dots),\quad(\mathrm{see}\
[5]).$

The incomplete and complete Bell polynomials have applications to such diverse
areas as combinatorics, probability, algebra and analysis. The number of
monomials appearing in the incomplete Bell polynomial
$B_{n,k}(x_{1},x_{2},\cdots,x_{n-k+1})$ is the number of partitioning $n$ into
$k$ parts and the coefficient of each monomial is the number of partitioning
$n$ as the corresponding $k$ parts. Also, the incomplete Bell polynomials
$B_{n,k}(x_{1},x_{2},\cdots,x_{n-k+1})$ appear in the Faà di Bruno formula
concerning higher-order derivatives of composite functions (see [5]). In
addition, the incomplete Bell polynomials can be used in constructing
sequences of binomial type (see [16]) and there are certain connections
between incomplete Bell polynomials and combinatorial Hopf algebras such as
the Hopf algebra of word symmetric functions, the Hopf algebra of symmetric
functions and the Faà di Bruno algebra, etc (see [1]). The complete Bell
polynomials $B_{n}(x_{1},x_{2},\cdots,x_{n})$ have applications to probability
theory (see [5,12,18]). Indeed, the $n$th moment $\mu_{n}=E[X^{n}]$ of the
random variable $X$ is the $n$th complete Bell polynomial in the first $n$
cumulants $\mu_{n}=B_{n}(\kappa_{1},\kappa_{2},\cdots,\kappa_{n})$. The reader
can refer to the Ph. D. thesis of Port [18] for many applications to
probability theory and combinatorics. Many special numbers, like Stirling
numbers of both kinds, Lah numbers and idempotent numbers, appear in many
combinatorial and number theoretic identities involving complete and
incomplete Bell polynomials. We let the reader refer to the Introduction in
[11] for further details on these.

The incomplete Lah-Bell polynomials (see (22)) and the complete Lah-Bell
polynomials (see (25)) are respectively multivariate versions of the unsigned
Lah numbers and the Lah-Bell numbers. We note here that the incomplete Bell
polynomials (see (10)) and the incomplete Lah-Bell polynomials are related as
given in (23), while the complete Bell polynomials (see (8)) and the complete
Lah-Bell polynomials are related as given in (26). The incomplete $r$-extended
Lah-Bell polynomials (see (30)) and the complete $r$-extended Lah-Bell
polynomials (see (32)) are respectively extended versions of the incomplete
Lah-Bell polynomials and the complete Lah-Bell polynomials. Further, they are
respectively multivariate versions of the $r$-Lah numbers and the $r$-extended
Lah-Bell numbers.

The aim of this paper is to introduce the incomplete $r$-extended Lah-Bell
polynomials and the complete $r$-extended Lah-Bell polynomials and to
investigate some properties and identities for these polynomials. From these
investigations, we obtain some expressions for the $r$-Lah numbers and the
$r$-extended Lah-Bell numbers as finite sums.

##  2\. Complete and incomplete $r$-extended Lah-Bell polynomials

Let $f(t)=\frac{t}{1-t}$. Then we have

(19) $f^{(n)}(t)=\frac{d^{n}}{dt^{n}}f(t)=\frac{n!}{(1-t)^{n+1}},\quad(n\geq
1).$

By (14), we get

(20)
$\frac{d^{n}}{dt^{n}}e^{\frac{t}{1-t}}\bigg{|}_{t=0}=B_{n}\big{(}1!,2!,\dots,n!\big{)}.$

From (2), we note that

(21)
$\frac{d^{n}}{dt^{n}}e^{\frac{t}{1-t}}\bigg{|}_{t=0}=\frac{d^{n}}{dt^{n}}\sum_{k=0}^{\infty}B_{k}^{L}\frac{t^{k}}{k!}\bigg{|}_{t=0}=B_{n}^{L}.$

Therefore, by (20) and (21), we obtain the following theorem.

######  Theorem 1.

For $n\geq 1$, we have

$B_{n}^{L}=B_{n}(1!,2!,\dots,n!)=\sum_{k_{1}+2k_{2}+\cdots+nk_{n}=n}\frac{n!}{k_{1}!k_{2}!\cdots
k_{n}!}.$

Let us consider the incomplete Lah-Bell polynomials given by

(22)
$\frac{1}{k!}\bigg{(}\sum_{m=1}^{\infty}x_{m}t^{m}\bigg{)}^{k}=\sum_{n=k}^{\infty}B_{n,k}^{L}\big{(}x_{1},x_{2},\dots,x_{n-k+1}\big{)}\frac{t^{n}}{n!},$

where $n,k\geq 0$ with $n\geq k$.

Note hat $B_{n,k}^{L}(1,1,\dots,1)=L(n,k)$, $(n\geq k\geq 0)$.

Indeed, by (10) and (22) , we get

(23)
$B_{n,k}^{L}\big{(}x_{1},x_{2},\dots,x_{n-k+1}\big{)}=B_{n,k}\big{(}1!x_{1},2!x_{2},\dots,(n-k+1)!x_{n-k+1}\big{)}.$

From (23), we note that

$\displaystyle L(n,k)\ $ $\displaystyle=\ B_{n,k}^{L}(1,1,\dots,1)\ =\
B_{n,k}\big{(}1!,2!,\dots,(n-k+1)!\big{)}$ $\displaystyle=\
\sum_{\begin{subarray}{c}j_{1}+j_{2}+\cdots+j_{n-k+1}=k\\\
j_{1}+2j_{2}+\cdots+(n-k+1)j_{n-k+1}=n\end{subarray}}\frac{n!}{j_{1}!j_{2}!\cdots
j_{n-k+1}!}$

Therefore, by (23), we obtain the following proposition.

######  Proposition 2.

For $n,k\geq 0$ with $n\geq k$, we have

$B_{n,k}^{L}\big{(}x_{1},x_{2},\dots,x_{n-k+1}\big{)}=B_{n,k}\big{(}1!x_{1},2!x_{2},\dots,(n-k+1)x_{n-k+1}\big{)}.$

In addition,

$L(n,k)=\sum_{\begin{subarray}{c}j_{1}+j_{2}+\cdots+j_{n-k+1}=k\\\
j_{1}+2j_{2}+\cdots+(n-k+1)j_{n-k+1}=n\end{subarray}}\frac{n!}{j_{1}!j_{2}!\cdots
j_{n-k+1}!}.$

From (23), we note that

(24) $\displaystyle B_{n,k}^{L}\big{(}\alpha x_{1},\alpha x_{2},\dots,\alpha
x_{n-k+1}\big{)}\ $ $\displaystyle=\ B_{n,k}\big{(}\alpha 1!x_{1},\alpha
2!x_{2},\dots,\alpha(n-k+1)!x_{n-k+1}\big{)}$ $\displaystyle=\
\alpha^{k}B_{n,k}\big{(}1!x_{1},2!x_{2},\dots,(n-k+1)!x_{n-k+1}\big{)}$
$\displaystyle=\ \alpha^{k}B_{n,k}^{L}(x_{1},x_{2},\dots,x_{n-k+1}).$

We now consider the complete Lah-Bell polynomials given by

(25)
$\exp\bigg{(}\sum_{i=1}^{\infty}x_{i}t^{i}\bigg{)}=\sum_{n=0}^{\infty}B_{n}^{L}(x_{1},x_{2},\dots,x_{n})\frac{t^{n}}{n!}.$

By (25), we get

(26) $\displaystyle B_{n}^{L}\big{(}x_{1},x_{2},\dots,x_{n}\big{)}\ $
$\displaystyle=\ B_{n}\big{(}1!x_{1},2!x_{2},\dots,n!x_{n}\big{)}$
$\displaystyle=\
\sum_{l_{1}+2l_{2}+\cdots+nl_{n-1}=n}\frac{n!}{l_{1}!l_{2}!\cdots
l_{n}!}x_{1}^{l_{1}}x_{2}^{l_{2}}\cdots x_{n}^{l_{n}},\quad(n\geq 0).$

From (22) and (25), we note that

(27) $\displaystyle
1+\sum_{n=1}^{\infty}B_{n}^{L}(x_{1},x_{2},\dots,x_{n})\frac{t^{n}}{n!}$
$\displaystyle=\exp\bigg{(}\sum_{i=1}^{\infty}x_{i}t^{i}\bigg{)}$
$\displaystyle=\
1+\sum_{k=1}^{\infty}\frac{1}{k!}\bigg{(}\sum_{i=1}^{\infty}x_{i}t^{i}\bigg{)}^{k}$
$\displaystyle=\
1+\sum_{k=1}^{\infty}\sum_{n=k}^{\infty}B_{n}^{L}(x_{1},x_{2},\dots,x_{n-k+1})\frac{t^{n}}{n!}$
$\displaystyle=\
1+\sum_{n=1}^{\infty}\bigg{(}\sum_{k=1}^{n}B_{n,k}^{L}(x_{1},x_{2},\dots,x_{n-k+1}\bigg{)}\frac{t^{n}}{n!}.$

Therefore, by (25) and (27), we obtain the following theorem.

######  Theorem 3.

For $n\geq 1$, we have

$\displaystyle B_{n}^{L}(x_{1},x_{2},\dots,x_{n})\ $ $\displaystyle=\
\sum_{k=1}^{n}B_{n,k}^{L}\big{(}x_{1},x_{2},\dots,x_{n-k+1}\big{)}$
$\displaystyle=\
\sum_{k=1}^{n}B_{n,k}\big{(}1!x_{1},2!x_{2},\dots,(n-k+1)!x_{n-k+1}\big{)}.$

In addition, for $n\geq 1$, we have

$B_{n}^{L}(x_{1},x_{2},\dots,x_{n})=\sum_{k=1}^{n}\sum_{\pi(n,k)}\frac{n!}{l_{1}!l_{2}!\cdots
l_{n-k+1}!}x_{1}^{l_{1}}x_{2}^{l_{2}}\cdots x_{n-k+1}^{l_{n-k+1}},$

where $\pi(n,k)$ denotes the set of all non-negative integers $(l_{i})_{i\geq
1}$ such that $l_{1}+l_{2}+\cdots+l_{n-k+1}=k$ and $1\cdot l_{1}+2\cdot
l_{2}+\cdots+(n-k+1)l_{n-k+1}=n$.

By (25), we easily get

(28)
$\sum_{n=0}^{\infty}B_{n}^{L}(1,1,\dots,1)\frac{t^{n}}{n!}=\exp\bigg{(}\sum_{i=1}^{\infty}t^{i}\bigg{)}=\exp\bigg{(}\frac{t}{1-t}\bigg{)}=\sum_{n=0}^{\infty}B_{n}^{L}\frac{t^{n}}{n!}.$

From (28), we note that

$B_{n}^{L}(1,1,\dots,1)=B_{n}^{L},\quad(n\geq 0).$

By Proposition 2, (24) and Theorem 3, we get

(29) $\displaystyle B_{n}^{L}(x,x,\dots,x)\ $ $\displaystyle=\
\sum_{k=0}^{n}B_{n,k}^{L}(x,x,\dots,x)\ =\
\sum_{k=0}^{n}B_{n,k}\big{(}1!x,2!x,\dots,(n-k+1)!x_{n-k+1}\big{)}$
$\displaystyle=\ \sum_{k=0}^{n}x^{k}B_{n,k}\big{(}1!,2!,\dots,(n-k+1)!\big{)}\
=\ \sum_{k=0}^{n}x^{k}L(n,k)\ =\ B_{n}^{L}(x).$

Assume that $\\{a_{i}\\}_{i\geq 1}$ and $\\{b_{i}\\}_{i\geq 1}$ are sequences
of positive integers. We define the incomplete $r$-extended Lah-Bell
polynomials by

(30)
$\frac{1}{k!}\bigg{(}\sum_{j=1}^{\infty}a_{j}t^{j}\bigg{)}^{k}\bigg{(}\sum_{i=0}^{\infty}b_{i+1}t^{i}\bigg{)}^{2r}=\sum_{n=k}^{\infty}B_{n+2r,k+2r}^{L}(a_{1},a_{2},\dots:b_{1},b_{2},\dots)\frac{t^{n}}{n!},$

where $k,r$ are nonnegative integers.

From (30), we have

(31) $\displaystyle
B_{n+2r,k+2r}^{L}(a_{1},a_{2},\dots:b_{1},b_{2},\dots)=B_{n+2r,k+2r}^{(2r)}(1!a_{1},2!a_{2},\dots:0!b_{1},1!b_{2},\dots)$
$\displaystyle\
=\sum_{\Lambda(n,k,2r)}\bigg{[}\frac{n!}{k_{1}!k_{2}!k_{3}!\cdots}a_{1}^{k_{1}}a_{2}^{k_{2}}a_{3}^{k_{3}}\cdots\bigg{]}\bigg{[}\frac{(2r)!}{r_{0}!r_{1}!r_{2}\cdots}b_{1}^{r_{0}}b_{2}^{r_{1}}b_{3}^{r_{2}}\cdots\bigg{]},$

where $\Lambda(n,k,2r)$ denotes the set of all nonnegative integers
$\\{k_{i}\\}_{i\geq 1}$ and $\\{r_{i}\\}_{i\geq 0}$ such that
$\displaystyle\sum_{i\geq 1}k_{i}=k,\ \sum_{i\geq 0}r_{i}=2r\displaystyle$ and
$\displaystyle\sum_{i\geq 1}i(k_{i}+r_{i})=n\displaystyle$.

We define the complete $r$-extended Lah-Bell polynomials
$B_{n}^{(L,2r)}(x|a_{1},a_{2},\dots:b_{1},b_{2},\dots)$, $(n\geq 0)$, which
are given by

(32)
$\exp\bigg{(}x\sum_{j=1}^{\infty}a_{j}t^{j}\bigg{)}\bigg{(}\sum_{i=0}^{\infty}b_{i+1}t^{i}\bigg{)}^{2r}=\sum_{n=0}^{\infty}B_{n}^{(L,2r)}(x|a_{1},a_{2},\dots:b_{1},b_{2},\dots)\frac{t^{n}}{n!}.$

Thus, we note that

(33)
$\displaystyle\exp\bigg{(}x\sum_{j=1}^{\infty}a_{j}t^{j}\bigg{)}\bigg{(}\sum_{i=0}^{\infty}b_{i+1}t^{i}\bigg{)}^{2r}\
=\
\sum_{k=0}^{\infty}\frac{x^{k}}{k!}\bigg{(}\sum_{j=1}^{\infty}a_{j}t^{j}\bigg{)}^{k}\bigg{(}\sum_{i=0}^{\infty}b_{i+1}t^{i}\bigg{)}^{2r}$
$\displaystyle\
=\sum_{k=0}^{\infty}x^{k}\sum_{n=k}^{\infty}B_{n+2r,k+2r}^{L}(a_{1},a_{2},\dots:b_{1},b_{2},\dots)\frac{t^{n}}{n!}\
=\
\sum_{n=0}^{\infty}\sum_{k=0}^{n}x^{k}B_{n+2r,k+2r}^{L}(a_{1},a_{2},\dots:b_{1},b_{2},\dots)\frac{t^{n}}{n!}.$

From (32) and (33), we have

(34) $B_{n}^{(L,2r)}\big{(}x\ |\
a_{1},a_{2},\dots:b_{1},b_{2},\dots)=\sum_{k=0}^{n}x^{k}B_{n+2r,k+2r}^{L}(a_{1},a_{2},\dots:b_{1},b_{2},\dots),\
(n\geq 0).$

By (18), (31), (32) and (34), we have

(35) $\displaystyle
B_{n}^{(2r)}(1!a_{1},2!a_{2},\dots:0!b_{1},1!b_{2},\dots)=\sum_{k=0}^{n}B_{n+2r,k+2r}^{(2r)}(1!a_{1},2!a_{2},\dots:0!b_{1},1!b_{2},\dots)$
$\displaystyle=\sum_{k=0}^{n}B_{n+2r,k+2r}^{L}(a_{1},a_{2},\dots:b_{1},b_{2},\dots)=B_{n}^{(L,2r)}(1|a_{1},a_{2},\dots:b_{1},b_{2},\dots).$

Therefore, by (31) and (34), we obtain the following theorem.

######  Theorem 4.

For $n\geq 0$, we have

$\displaystyle B_{n}^{(L,2r)}\big{(}x\ |\
a_{1},a_{2},\dots:b_{1},b_{2},\dots\big{)}\ $ $\displaystyle=\
\sum_{k=0}^{n}x^{k}B_{n+2r,k+2r}^{L}(a_{1},a_{2},\dots:b_{1},b_{2},\dots)$
$\displaystyle=\
\sum_{k=0}^{n}x^{k}B_{n+2r,k+2r}^{(2r)}(1!a_{1},2!a_{2},\dots:0!b_{1},1!b_{2},\dots).$

From (30), we note that

(36)
$\displaystyle\sum_{n=k}^{\infty}B_{n+2r,k+2r}^{L}\big{(}1,1,\dots;1,1,\dots\big{)}\frac{t^{n}}{n!}\
=\
\frac{1}{k!}\bigg{(}\sum_{j=1}^{\infty}t^{j}\bigg{)}^{k}\bigg{(}\sum_{i=0}^{\infty}t^{i}\bigg{)}^{2r}$
$\displaystyle\ =\
\frac{1}{k!}\bigg{(}\frac{t}{1-t}\bigg{)}^{k}\bigg{(}\frac{1}{1-t}\bigg{)}^{2r}\
=\ \sum_{n=k}^{\infty}L_{r}(n,k)\frac{t^{n}}{n!}.$

By (32) and (36), we get

(37)
$\displaystyle\sum_{n=0}^{\infty}B_{n}^{(L,2r)}(x|1,1,\dots:1,1,\dots)\frac{t^{n}}{n!}\
=\
\exp\bigg{(}x\sum_{j=1}^{\infty}t^{j}\bigg{)}\bigg{(}\sum_{i=0}^{\infty}t^{i}\bigg{)}^{2r}$
$\displaystyle=\ e^{x(\frac{t}{1-t})}\cdot\bigg{(}\frac{1}{1-t}\bigg{)}^{2r}\
=\
\sum_{n=0}^{\infty}\bigg{(}\sum_{k=0}^{n}x^{k}L_{r}(n,k)\bigg{)}\frac{t^{n}}{n!}.$

Thus, by (36) and (37), we have

(38) $\displaystyle B_{n}^{(L,2r)}\big{(}x\ |\ 1,1,\dots:1,1,\dots)\ $
$\displaystyle=\ \sum_{k=0}^{n}x^{k}B_{n+2r,k+2r}^{L}(1,1,\dots:1,1,\dots)$
$\displaystyle=\ \sum_{k=0}^{n}x^{k}L_{r}(n,k).$

Therefore, we obtain the following theorem.

######  Theorem 5.

For $n\geq k\geq 0$, we have

$B_{n+2r,k+2r}^{L}(1,1,\dots:1,1,\dots)=L_{r}(n,k)$

and

$B_{n}^{(L,2r)}(x\ |\
1,1,\dots:1,1,\dots)=\sum_{k=0}^{n}B_{n+2r,k+2r}^{L}(1,1,\dots:1,1,\dots)=\sum_{k=0}^{n}x^{k}L_{r}(n,k).$

From (36) and (31), we note that

$\displaystyle L_{r}(n,k)\ $ $\displaystyle=\
B_{n+2r,k+2r}^{L}(1,1,\dots:1,1,\dots)$ $\displaystyle=\
B_{n+2r,k+2r}^{(2r)}(1!,2!,\dots:1!,2!,\dots)$ $\displaystyle=\
\sum_{\Lambda(n,k,2r)}\frac{n!}{k_{1}!k_{2}!\cdots}\frac{(2r)!}{r_{0}!r_{1}!\cdots}.$

######  Corollary 6.

For $n,k,r\geq 0$ with $n\geq k$, we have

$L_{r}(n,k)=\sum_{\Lambda(n,k,2r)}\frac{n!}{k_{1}!k_{2}!\cdots}\cdot\frac{(2r)!}{r_{0}!r_{1}!\cdots},$

where $\Lambda(n,k,2r)$ denotes the set of all nonnegative integers
$\\{k_{i}\\}_{i\geq 1}$ and $\\{r_{i}\\}_{i\geq}$ such that

$\displaystyle\sum_{i\geq 1}k_{i}=k,\ \sum_{i\geq 0}r_{i}=2r\displaystyle$ and
$\displaystyle\sum_{i\geq 1}i(k_{i}+r_{i})=n.\displaystyle$

Now, we observe that

(39) $\displaystyle\exp\bigg{(}\sum_{i=1}^{\infty}x_{i}t^{i}\bigg{)}\ $
$\displaystyle=\
1+\sum_{k=1}^{\infty}\frac{1}{k!}\bigg{(}\sum_{i=1}^{\infty}x_{i}t^{i}\bigg{)}^{k}$
$\displaystyle=\
1+\frac{1}{1!}\sum_{i=1}^{\infty}x_{i}t^{i}+\frac{1}{2!}\bigg{(}\sum_{i=1}^{\infty}x_{i}t^{i}\bigg{)}^{2}+\frac{1}{3!}\bigg{(}\sum_{i=1}^{\infty}x_{i}t^{i}\bigg{)}^{3}+\cdots$
$\displaystyle=\
\sum_{k=0}^{\infty}\sum_{m_{1}+2m_{2}+\cdots+km_{k}=k}\frac{1}{m_{1}!m_{2}!\cdots
m_{k}!}x_{1}^{m_{1}}x_{2}^{m_{2}}\cdots x_{k}^{m_{k}}t^{k},$

and

(40)
$\bigg{(}\sum_{j=0}^{\infty}y_{j+1}t^{j}\bigg{)}^{2r}=\sum_{m=0}^{\infty}\sum_{l_{1}+\cdots+l_{2r}=m}y_{l_{1}+1}y_{l_{2}+1}\cdots
y_{l_{2r}+1}t^{m}.$

By (39) and (40), we get

(41)
$\displaystyle\exp\bigg{(}\sum_{i=1}^{\infty}x_{i}t^{i}\bigg{)}\bigg{(}\sum_{j=0}^{\infty}y_{j+1}\frac{t^{j}}{j!}\bigg{)}^{2r}$
$\displaystyle=\
\sum_{n=0}^{\infty}n!\bigg{(}\sum_{k=0}^{n}\sum_{m_{1}+2m_{2}+\cdots+km_{k}=k}\,\sum_{l_{1}+l_{2}+\cdots+l_{2r}=n-k}\frac{1}{m_{1}!m_{2}!\cdots
m_{k}!}x_{1}^{m_{1}}x_{2}^{m_{2}}\cdots x_{k}^{m_{k}}$
$\displaystyle\qquad\quad\times y_{l_{1}+1}y_{l_{2}+1}\cdots
y_{l_{2r}+1}\bigg{)}\frac{t^{n}}{n!}.$

Therefore, by (32) and (41), we obtain the following theorem.

######  Theorem 7.

For $n,r\geq 0$, we have

$\displaystyle B_{n}^{(L,2r)}(1\ |\ x_{1},x_{2},\dots:y_{1},y_{2},\dots)$
$\displaystyle=\
n!\sum_{k=0}^{n}\sum_{m_{1}+2m_{2}+\cdots+km_{k}=k}\,\sum_{l_{1}+l_{2}+\cdots+l_{2r}=n-k}\frac{1}{m_{1}!m_{2}!\cdots
m_{k}!}x_{1}^{m_{1}}\cdots x_{k}^{m_{k}}y_{l_{1}+1}y_{l_{2}+1}\cdots
y_{l_{2r}+1}.$

_Remark._ For $n\geq 0$, we have

(42)
$B_{n+2r,k+2r}^{L}(x,x,\dots:1,1,\dots)=x^{k}B_{n+2r,k+2r}^{L}(1,1,\dots:1,1,\dots).$

Thus, we note that

$\sum_{k=0}^{n}B_{n+2r,k+2r}^{L}(x,x,\dots:1,1,\dots)=B_{n,r}^{L}(x),\quad(n\geq
0).$

##  3\. Conclusion

There are various methods of studying special numbers and polynomials, for
example, generating functions, combinatorial methods, umbral calculus,
$p$-adic analysis, differential equations, probability theory, orthogonal
polynomials, and special functions. These ways of investigating special
polynomials and numbers can be applied also to degenerate versions of such
polynomials and numbers. Indeed, in recent years, many mathematicians have
drawn their attention to studies of degenerate versions of many special
polynomials and numbers by making use of the aforementioned means [9,10,14 and
references therein].

The incomplete and complete Bell polynomials arise in many different contexts
as we stated in the Introduction. For instance, many special numbers, like
Stirling numbers of both kinds, Lah numbers and idempotent numbers, appear in
many combinatorial and number theoretic identities involving complete and
incomplete Bell polynomials.

In this paper, we introduced the incomplete $r$-extended Lah-Bell polynomials
and the complete $r$-extended Lah-Bell polynomials respectively as
multivariate versions of $r$-Lah numbers and the $r$-extended Lah-Bell numbers
and investigated some properties and identities for these polynomials. As
corollaries to these results, we obtained some expressions for the $r$-Lah
numbers and the $r$-extended Lah-Bell numbers as finite sums.

It would be very interesting to explore many applications of the incomplete
and complete $r$-extended Lah-Bell polynomials just as the incomplete and
complete Bell polynomials have diverse applications.

Acknowledgements

The authors would like to thank the reviewers for their valuable comments and
suggestions and Jangjeon Institute for Mathematical Science for the support of
this research.

Funding

Not applicable.

Availability of data and materials

Not applicable.

Competing interests

The authors declare that they have no conflicts of interest.

Authors’ contributions

TK and DSK conceived of the framework and structured the whole paper; DSK and
TK wrote the paper; LCJ, HL, and HYK checked the results of the paper; DSK and
TK completed the revision of the paper. All authors have read and approved the
final version of the manuscript.

## References

  * [1] Aboud, A.; Bultel, J.-P.; Chouria, A.; Luque, J.-G.; Mallet, O. _Word Bell polynomials,_ Sém. Lothar. Combin. 75 (2017), B75h.
  * [2] Abramowitz, M.; Stegun, I. A. _Handbook of mathematical functions with formulas, graphs, and mathematical tables,_ National Bureau of Standards Applied Mathematics Series 55, For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C., 1964.
  * [3] Broder, A.Z. _The $r$-Stirling numbers,_ Discrete Math. 49 (3), 241-–259 (1984).
  * [4] Carlitz, L. _Arithmetic properties of the Bell polynomials,_ J. Math. Anal. Appl. 15 (1966), 33-–52.
  * [5] Comtet, L. _Advanced combinatorics. The art of finite and infinite expansions,_ Revised and enlarged edition, D. Reidel Publishing Co., Dordrecht, 1974.
  * [6] Coffey, M. W. _A set of identites for a class of altenating binomial sums arising in computing applications,_ arXiv:math-ph/0608049.
  * [7] Kim, D. S.; Kim, T. _$r$ -extended Lah-Bell numbers and polynomials associated with $r$-Lah numbers,_ arXiv:2008.06155.
  * [8] Kim, D. S.; Kim, T. _Lah-Bell numbers and polynomials,_ Proc. Jangjeon Math. Soc. 23 (2020), no. 4, 577–586.
  * [9] Kim, H. K. _Degenerate Lah-Bell polynomials arising from degenerate Sheffer sequences,_ Adv. Difference Equ. 2020, 687.
  * [10] Kim, H. K; Lee, D. S. _Note on extended Lah-Bell polynomials and degenerate extended Lah-Bell polynomials,_ Adv. Stud. Contemp. Math. (Kyungshang) 30 (2020), no. 4, 547–558.
  * [11] Kim, T.; Kim, D.S.; Jang, G.-W. _On central complete and incomplete Bell polynomials I,_ Symmetry 11 (2019), no. 2, Article ID 288.
  * [12] Kölbig, K. S. _The complete Bell polynomials for cetain arguments in terms of stirling numbers of the first kind,_ J. Comput. Appl. Math. 51 (1994), 113–116.
  * [13] Kucukoglu, I.: Simsek,Y. _Observations on identities and relations for interpolation fucntions and special numbers,_ Adv. Stud. Contemp. Math. (Kyungshang) 28 (2018), no. 1, 41-–56.
  * [14] Kwon, J.; Kim, T.; Kim, D. S.; Kim, H. Y. _Some identities for degenerate complete and incomplete $r$-Bell polynomials,_ J. Inequal. Appl. 2020, Paper No. 23, 9 pp.
  * [15] Ma, Y.; Kim, D. S.; Kim, T.; Kim, H.; Lee, H. _Some identities of Lah-Bell polynomials,_ Adv. Difference Equ. 2020, Paper no. 510, 10 pp.
  * [16] Mihoubi, M. _Bell polynomials and binomial type sequences,_ Discrete Math. 308 (2008), no. 12, 2450–-2459.
  * [17] Nyul, G.; Racz, G. _The $r$-Lah numbers,_ Discrete Math. 338 (2015), no. 10, 1660-–1666.
  * [18] Port, D. _Polynomial maps with applications to combinatorics and probability theory,_ Thesis (Ph.D.)–Massachusetts Institute of Technology, 1994
  * [19] Roman, S. _The umbral calculus,_ Pure and Applied Mathematics, 111. Academic Press, Inc. [Harcourt Brace Jovanovich, Publishers], New York, 1984.