\documentclass[12pt,reqno]{article}

\usepackage[usenames]{color}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{amsfonts}
\usepackage{amscd}
\usepackage{graphicx}

\usepackage[colorlinks=true,
linkcolor=webgreen,
filecolor=webbrown,
citecolor=webgreen]{hyperref}

\definecolor{webgreen}{rgb}{0,.5,0}
\definecolor{webbrown}{rgb}{.6,0,0}

\usepackage{color}
\usepackage{fullpage}
\usepackage{float}

\usepackage{psfig}
\usepackage{graphics}
\usepackage{latexsym}
\usepackage{epsf}
\usepackage{breakurl}

\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{.1in}
\setlength{\evensidemargin}{.1in}
\setlength{\topmargin}{-.1in}
\setlength{\textheight}{8.4in}

\newcommand{\seqnum}[1]{\href{https://oeis.org/#1}{\rm \underline{#1}}}



\begin{document}
\begin{center}
\epsfxsize=4in
\leavevmode\epsffile{logo129.eps}
\end{center}

\theoremstyle{plain}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}

\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{conjecture}[theorem]{Conjecture}

\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{question}[theorem]{Question}

\begin{center}
\vskip 1cm{\Large\bf 
Polynomial Identities Involving Binomial Coefficients \\
\vskip .01in
and Double and Rising Factorials \\
\vskip .1in
via Probabilistic Interpretations
}
\vskip 1cm
\large
Pawe\l \ J. Szab\l owski\footnote{Retired.} \\
Department of Mathematics and Information Sciences \\
Warsaw University of Technology\\
ul. Koszykowa 75 \\
00-662 Warsaw \\
Poland \\
\href{mailto:pawel.szablowski@gmail.com}{\tt pawel.szablowski@gmail.com}
\end{center}

\vskip .2 in

\begin{abstract}
We formulate several polynomial identities. One side of these identities has a
nice simple form, whereas the other has the form of a polynomial whose
coefficients contain binomial coefficients, double factorials, or rising
factorials. The origins and the proofs of these identities are probabilistic.
However, their form suggests universal applications in simplifying
expressions. Many useful simplifying formulae are presented in the sequel.
\end{abstract}

\section{Introduction}

In this note, we will present some identities involving factorials, binomial
coefficients, and the so-called rising factorials (sometimes called Pochhammer
symbols). Many of them have the forms of polynomials whose coefficients often
have the form of binomial coefficients or rising factorials of some
additional variables. These variables appear on both sides of the identity.
The domains of these variables can be extended to all complex numbers. All
these identities have a common origin, except for the sometimes similar form.
Namely, the origins of them are a few, rather deep probabilistic
interpretations and sometimes following non-trivial computations.

The basic
idea behind the method of getting these identities is to calculate moments of
a probability distributions in two ways and compare the results. One of the
ways is straightforward and direct. The other way involves the calculation of
some conditional expectation and then calculating the proper moment. This
method was already successfully applied to obtain nontrivial relationships
between Catalan numbers and the so-called Catalan triangles. In this case, the
so-called Kesten distribution was used. For details, see \cite{SzabKes}. In
this paper, we use the bivariate Normal and Gamma distributions to obtain
nontrivial polynomial identities (see Theorem \ref{Main}). These identities
can be the source of an infinite number of relationships between various
number sequences important for example in combinatorics. As a corollary, we
present several identities involving binomial coefficients, factorials, double
factorials, and raising factorials.

The paper is organized as follows. The next section is dedicated to the
presentation of the identities and some of the particular, interesting
particular cases. The next section is devoted to the presentation of the
probabilistic background of the results presented in the previous section and
then, finally, the presentation of the calculations leading to the identities.

\section{Identities}

In this paper, the symbol $n!!$ denotes the so-called double-factorial, i.e.,
$$n!!=\prod_{j=0}^{\left\lfloor (n-1)/2\right\rfloor }(n-2j).$$ 
We set $(-1)!!=(0)!!=1$.
The symbol $\left\lfloor x\right\rfloor $ denotes the largest integer not
exceeding $x$, i.e., the so-called \textquotedblleft floor
function\textquotedblright. Further, let the symbol $\binom{n}{k}$ denote the
binomial coefficient and we set $\binom{n}{k}=0$ when
$n<k$. The symbol $i$ always denotes the imaginary unit, i.e., $i
=\sqrt{-1}=\exp(i\pi/2)$. In order to
simplify the notation, let us introduce also the so-called \textquotedblleft
rising factorial\textquotedblright\ (sometimes called also Pochhammer symbol),
which is the following function:
\[
(x)^{(n)}=x(x+1)\cdots(x+n-1),
\]
defined for all complex $x$. Notice that we have for all $x\neq0$:
\[
(x)^{(n)}=\frac{\Gamma(x+n)}{\Gamma(x)}\text{,}
\]
where $\Gamma(x)$ denotes Euler's gamma function. To learn more about
Pochhammer symbols and binomial coefficients, see, e.g., \cite{Grad14}. Notice
only that, e.g., $\left(  1\right)  ^{\left(  n\right)  }
= n!$ and $\left(  1/2\right)  ^{\left(  n\right)  }
=(2n-1)!!/2^{n}$. We will use these values below.

\begin{theorem}
\label{Main}For all nonegative integers $k,n,m$ and all complex $\rho$ and
$\beta$ we have

\begin{enumerate}
\item[(i)]
\begin{equation}
\frac{(2k)!}{k!}(1+\rho)^{k}=\sum_{j=0}^{k}\binom{2k}{2j}(1+\rho
)^{2k-2j}(1-\rho^{2})^{j}(2j-1)!!(2k-2j-1)!!. \label{G1}
\end{equation}


\item[(ii)]
\begin{equation}
(1+\rho)^{k}=\frac{1}{2^{k}}\sum_{j=0}^{2k}\sum_{m=0}^{\left\lfloor
j/2\right\rfloor }(2\rho)^{j-2m}\binom{k}{j-2m}\binom{k-j+2m}{m}. \label{G2}
\end{equation}


\item[(iii)]
\begin{equation}
\sum_{j=0}^{n}(-1)^{j}\binom{n}{j}\frac{\left(  \beta\right)  ^{\left(
j+m\right)  }}{\left(  \beta\right)  ^{\left(  j\right)  }}=(-1)^{n}
n!\binom{m}{n}\frac{\left(  \beta\right)  ^{\left(  m\right)  }}{\left(
\beta\right)  ^{\left(  n\right)  }}.\label{Ex1}
\end{equation}


\item[(iv)]
\begin{gather}
\sum_{m=0}^{n}\left(  -1\right)  ^{m}\binom{n}{m}\left(  \beta\right)
^{\left(  m\right)  }\left(  \beta\right)  ^{\left(  n-m\right)  }\sum
_{k=0}^{m}\binom{m}{k}\binom{n-m}{k}\frac{\rho^{k}k!}{\left(  \beta\right)
^{\left(  k\right)  }}\label{Ex2}\\
=
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{n!}{(n/2)!}\left(  \beta\right)  ^{\left(  n/2\right)  }(1-\rho
)^{n/2}, & \text{if }n\text{ is even.}
\end{cases}
\nonumber
\end{gather}


\item[(v)]
\begin{align}
&  \sum_{m=0}^{n}(-1)^{n-m}\binom{n}{m}\sum_{j=0}^{m}\binom{m}{j}(1-\rho
)^{j}\rho^{m-j}\left(  \beta\right)  ^{\left(  n-j\right)  }\left(
\beta+m-j\right)  ^{\left(  j\right)  }\label{Ex3}\\
&  = \quad\quad
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{n!}{(n/2)!}\left(  \beta\right)  ^{\left(  n/2\right)  }(1-\rho
)^{n/2}, & \text{if }n\text{ is even.}
\end{cases}
\nonumber
\end{align}

\end{enumerate}
\end{theorem}

Below we present some particular cases, remarks and corollaries.

\begin{remark}
First, taking $\rho=\frac{2}{3}$, then $\rho
=\frac{1}{3}$,  and finally $\rho
=\frac{4}{3}$ in (\ref{G1}),  we get
\begin{align*}
\sum_{j=0}^{k}\binom{2k}{2j}5^{j}(2j-1)!!(2k-2j-1)!!  &  =\frac{(2k)!3^{k}
}{k!},\\
\sum_{j=0}^{k}\binom{2k}{2j}2^{j}(2j-1)!!(2k-2j-1)!!  &  =\frac{(2k)!3^{k}
}{k!2^{k}},\\
\sum_{j=0}^{k}\binom{2k}{2j}(-7)^{j}(2j-1)!!(2k-2j-1)!!  &  =(-1)^{k}
\frac{(2k)!3^{k}}{k!}.
\end{align*}

\end{remark}

\begin{remark}
Let us take $\rho=\sqrt{5}$ in (\ref{G1}) and recall
that $(1+\sqrt{5})^{n}=2^{n}(\frac{1+\sqrt{5}}{2}
)^{n}=2^{n-1}(L_{n}+\sqrt{5}F_{n})$, where $L_{n}$ and
$F_{n}$ are, respectively, the $n$-th Lucas and Fibonacci numbers (also compare
Corollary \ref{Fib}, below). Now, equating integer and irrational parts, we end
up with the following two identities valid for $k\geq0$:
\begin{align*}
\frac{(2k)!}{k!}L_{k}  &  =2^{k}\sum_{j=0}^{k}\binom{2k}{2j}(-1)^{j}
L_{2k-2j}(2j-1)!!(2k-2j-1)!!,\\
\frac{(2k)!}{k!}F_{k}  &  =2^{k}\sum_{j=0}^{k}\binom{2k}{2j}(-1)^{j}
F_{2k-2j}(2j-1)!!(2k-2j-1)!!
\end{align*}

\end{remark}

\begin{remark}
Take $\rho= i$ in (\ref{G1}). Further notice that
$(1+i)^{k}=2^{k/2}\exp(ik\pi/4)$. Now first set
$k=4n$ and cancel $(-4)^{n}$ on both sides. Then for
all $n\geq0$ we get
\[
\frac{(8n)!}{(4n)!}=\sum_{j=0}^{4n}\binom{8n}{2j}(2j-1)!!(8n-2j-1)!!.
\]
Next, set $k=4n+1$ and equate the real and imaginary
parts. For all $n\geq0$ we get
\begin{align*}
\frac{(8n+2)!}{(4n+1)!}(-4)^{n}  &  =2^{4n+1}\sum_{m=0}^{2n}(-1)^{m}
\binom{8n+2}{4m+2}(4m+1)!!(8n-4m-1)!!,\\
\frac{(8n+2)!}{(4n+1)!}(-4)^{n}  &  =2^{4n+1}\sum_{m=0}^{2n}(-1)^{m}
\binom{8n+2}{4m}(4m-1)!!(8n-4m+1)!!.
\end{align*}

\end{remark}

\begin{remark}
Take $\rho=1/2$ in (\ref{G2}). We then get the following
identity:
\[
1=\sum_{j=0}^{2k}(-1)^{j}\sum_{m=0}^{\left\lfloor j/2\right\rfloor }\binom
{k}{j-2m}\binom{k-j+2m}{m}.
\]

\end{remark}

\begin{remark}
Take $\rho= i$ and $k=4n$ in
(\ref{G2}). Further, notice that $(1+i)^{4n}=(-4)^{n}$
and $(2i)^{j-2m}= i^{j}(-1)^{m}$ and finally split the
right-hand side into real and imaginary parts. We get 
\begin{align*}
4^{n}  &  =\frac{(-1)^{n}}{16^{n}}\sum_{j=0}^{4n}(-1)^{j}\sum_{m=0}
^{j}(-1)^{m}4^{j-m}\binom{4n}{2j-2m}\binom{4n-2j+2m}{m},\\
0  &  =\sum_{j=0}^{4n}(-1)^{j}\sum_{m=0}^{j}(-1)^{m}4^{j-m}\binom{4n}
{2j-2m+1}\binom{4n-2j-1+2m}{m},
\end{align*}
for all $n\geq 0$.

\end{remark}

\begin{remark}
First set $\beta=1$ and then $\beta
=1/2$ in (\ref{Ex1}). Then for all integers $k$ and $n$, we get
\begin{align*}
\sum_{j=0}^{n}(-1)^{j}\binom{n}{j}\binom{j+k}{j}  &  =(-1)^{n}\binom{k}{n},\\
\sum_{j=0}^{n}(-1)^{j}\frac{(2n-1)!!(2j+2k-1)!!}{j!(n-j)!(2j-1)!!}  &
=(-1)^{n}2^{n}\binom{k}{n}(2k-1)!! .
\end{align*}

\end{remark}

\begin{remark}
First set $\beta=1/2$ and then $\beta
=1$ in (\ref{Ex2}). Then for all integers $k$ and $n$ and all
complex $\rho$, we get
\begin{align}
\sum_{j=0}^{n}(-1)^{n-j}\binom{n}{j}\frac{(2n+2k+2j-1)!!}{(2j-1)!!}  &
=2^{n}\frac{(n+k)!(2n+2k-1)!!}{k!(2n-1)!!},\nonumber\\
\sum_{m=0}^{n}\left(  -1\right)  ^{m}\sum_{k=0}^{m}\binom{m}{k}\binom{n-m}
{k}\rho^{k}  &  =
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
(1-\rho)^{n/2}, & \text{if }n\text{ is even.}
\end{cases}
\label{explr}
\end{align}

\end{remark}

\begin{remark}
Now let us change the order of summation in (\ref{explr}) and compare the
coefficients in expansions in powers of $\rho$. We get
\[
\sum_{m=k}^{n}(-1)^{m-k}\binom{m}{k}\binom{n-m}{k}=
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
\binom{n/2}{k}, & \text{if }n\text{ is even.}
\end{cases}
\]

\end{remark}

\begin{remark}
Let us set additionally $\rho=1$ in (\ref{explr}). Then
for all $n\geq0$, we get
\[
\sum_{m=0}^{n}\left(  -1\right)  ^{m}\sum_{k=0}^{m}\binom{m}{k}\binom{n-m}
{k}=0.
\]

\end{remark}

\begin{remark}
Let us take $\beta=1/2$ in (\ref{Ex2}). Next, let us set
$2\rho= x$ and cancel $n!$. We then get 
\begin{gather}
\sum_{m=0}^{n}\left(  -1\right)  ^{m}(2m-1)!!(2n-2m-1)!!\sum_{k=0}^{m}\frac
{1}{k!(m-k)!(n-m-k)!}\frac{x^{k}}{(2k-1)!!}\label{Ex21}\\
=
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{(n-1)!!}{(n/2)!}(2-x)^{n/2}, & \text{if }n\text{ is even.}
\end{cases}
\nonumber
\end{gather}
for all $n\geq0$.
Additionally, setting $x=0$ in (\ref{Ex21}), we get
\[
\sum_{m=0}^{n}\left(  -1\right)  ^{m}\binom{n}{m}(2m-1)!!(2n-2m-1)!!=
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{n!(n-1)!!}{(n/2)!}2^{n/2}, & \text{if }n\text{ is even.}
\end{cases}
\]
\end{remark}

\begin{remark}
Let us set $\rho=1$ in (\ref{Ex2}). For all complex
$\beta\neq0$ and $n\geq0$ we get
\[
\sum_{m=0}^{n}\left(  -1\right)  ^{m}\binom{n}{m}\left(  \beta\right)
^{\left(  m\right)  }\left(  \beta\right)  ^{\left(  n-m\right)  }\sum
_{k=0}^{m}\binom{m}{k}\binom{n-m}{k}\frac{k!}{\left(  \beta\right)  ^{\left(
k\right)  }}=0.
\]

\end{remark}

\begin{remark}
Let us set $\rho=0$ in (\ref{Ex3}). We get 
\[
\sum_{m=0}^{n}(-1)^{n-m}\binom{n}{m}\left(  \beta\right)  ^{\left(
n-m\right)  }\left(  \beta\right)  ^{\left(  m\right)  }=
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{n!}{(n/2)!}\left(  \beta\right)  ^{\left(  n/2\right)  }, & \text{if
}n\text{ is even,}
\end{cases}
\]
for all $n\geq0$.
This is related to, but not a direct generalization of, the equality 
\[
\sum_{m=0}^{n}\binom{n}{m}\left(  \beta\right)  ^{\left(  n-m\right)  }\left(
\alpha\right)  ^{\left(  m\right)  }=\left(  \alpha+\beta\right)  ^{\left(
n\right)  }.
\]

\end{remark}

\begin{remark}
Let us set $\rho=1/2$ in (\ref{Ex3}). Then for all
complex $\beta$ and $n\geq0$, we get
\begin{gather*}
\sum_{m=0}^{n}(-1)^{n-m}\binom{n}{m}2^{-m}\sum_{j=0}^{m}\binom{m}{j}\left(
\beta\right)  ^{\left(  n-j\right)  }\left(  \beta+m-j\right)  ^{\left(
j\right)  }\\
=
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{n!}{2^{n/2}(n/2)!}\left(  \beta\right)  ^{\left(  n/2\right)  }, &
\text{if }n\text{ is even.}
\end{cases}
\end{gather*}

\end{remark}

\begin{remark}
Let us change the order of summation in (\ref{Ex3}), then compare the
coefficients in expansions in powers of $\rho$, then multiply both sides by
$\left(  \beta\right)  ^{\left(  k\right)  }(n-k)!$ and finally cancel 
$n!$. Then for all complex $\beta$ and $0\leq k\leq n$ we get
\begin{align*}
&  \sum_{m=k}^{n}(-1)^{m-k}\frac{(n-k)!}{(m-k)!(n-m-k)!}\left(  \beta\right)
^{\left(  m\right)  }\left(  \beta\right)  ^{\left(  n-m\right)  }\\
&  =
\begin{cases}
0, & \text{if }n\text{ is odd; }\\
(n/2)!\binom{n-k}{n/2}\left(  \beta\right)  ^{\left(  k\right)  }\left(
\beta\right)  ^{\left(  n/2\right)  }, & \text{if }n\text{ is even.}
\end{cases}
\end{align*}
\end{remark}

We also have the following relationships between Fibonacci $\left\{
F_{n}\right\}  _{n\geq0}$ and Lucas $\left\{  L_{n}\right\}  _{n\geq0}$
numbers. For Fibonacci and Lucas numbers, see, e.g.,
\cite{Ball03,Cull05} and \cite[pp.~153--256]{Knuth98}. 
For other relationships between Fibonacci and Lucas
numbers and the so-called Catalan triangles (expressed by binomial
coefficients), see, e.g., \cite{Szabl21}.

\begin{corollary}
\label{Fib}Let $\phi$ denote the number $(\sqrt
{5}-1)  /2$, which is the reciprocal of the golden
ratio. Take $\rho=-\phi$ in
(\ref{G1}). We get
\begin{align*}
2\frac{\left(  2k\right)  !}{k!}L_{k} &  =\sum_{j=0}^{k}(-1)^{j}\binom{2k}
{2j}(2j-1)!!(2k-2j-1)!!(L_{2k-2j}L_{j}-5F_{2k-2j}F_{j}),\\
2\frac{\left(  2k\right)  !}{k!}F_{k} &  =\sum_{j=0}^{k}(-1)^{j}\binom{2k}
{2j}(2j-1)!!(2k-2j-1)!!(F_{2k-2j}L_{j}-L_{2k-2j}F_{j}).
\end{align*}

\end{corollary}

\begin{proof}
We know that for all $n\geq0$ we have
\[
(1+\phi)^{n}=((1+\sqrt{5})/2)^{n}=L_{n}/2+F_{n}\sqrt{5}/2.
\]
Now, also notice that we have $\phi=1/(1+\phi)$ and
that
\[
L_{n}^{2}-5F_{n}^{2}=(-1)^{n}4.
\]
Hence, we have
\[
\phi^{n}=(-1)^{n}(L_{n}/2-F_{n}\sqrt{5}/2).
\]
Now we take $\rho=-\phi$ and substitute it in (\ref{G1}).
Notice that $1-\rho=(1+\phi)$, and hence $(1-\rho
)^{2k-2j}= L_{2k-2j}/2+
F_{2k-2j}\sqrt{5}/2$ and $1-\rho^{2}=\phi$, consequently
$(1-\rho^{2})^{j}=(-1)^{j}(L_{j}/2
- F_{j}\sqrt{5}/2)$. It remains to multiply and
separate terms with and without $\sqrt{5}$.
\end{proof}

\section{Probabilistic background and the proofs}

\subsection{Probabilistic background}

All identities presented in Theorem \ref{Main} stem from calculating the
moments $E(X-Y)^{n}$, $n\geq0$, where $X$ and $Y$ are two normalized
random variables. The joint distribution of these random variables has one
parameter and is either the bivariate Normal (Gaussian) or bivariate Gamma
distribution. More precisely, in the first case, the joint distribution of
$(X,Y)$ has the following well-known density, valid for all $x,y\in\mathbb{R}$
and $\left\vert \rho\right\vert <1$:
\begin{equation}
f_{N}(x,y|\rho)=\frac{1}{2\pi(1-\rho^{2})}\exp\left(  -\frac{x^{2}-2\rho
xy+y^{2}}{2(1-\rho^{2})}\right)  , \label{2xN}
\end{equation}
while 
in the second case 
\[
f_{G}(x,y|\rho)=f_{g}(x|\beta)f_{g}(y|\beta)\frac{\exp\left(  -\frac
{\rho\left(  x+y\right)  }{1-\rho}\right)  }{\left(  1-\rho\right)  \left(
xy\rho\right)  ^{(\beta-1)/2}}I_{\beta-1}\left(  \frac{2\sqrt{xy\rho}}{1-\rho
}\right)  ,
\]
valid for all $x,y,\beta>0$ and $\left\vert \rho\right\vert <1$,
where
\[
f_{g}(x|\beta)= x^{\beta-1}\exp(-x)/\Gamma(\beta),
\]
for $x>0$ and $0$ otherwise. We let $I_{\alpha}$ denote, here and below, the
modified Bessel function of the first kind. The most important property of
these joint densities is that they allow the so-called Lancaster expansions.
To learn more about Lancaster expansions, their probabilistic interpretations,
and convergence problems associated with them, see, e.g.,
\cite{Alexits61,Lancaster58,Lancaster75,Szeg75,Szabl21,SzabStac}. The
orthogonal polynomials that we are using are well presented, for example, in
two monographs \cite{Chih79,Szeg75}.

In the first case, we have the so-called Poisson-Mehler expansion
\begin{equation}
f_{N}(x,y|\rho)=f_{N}(x)f_{N}(x)\sum_{n\geq0}\rho^{n}h_{n}(x)h_{n}(y),
\label{MN}
\end{equation}
where we write
\[
f_{N}(x)=\exp(-x^{2}/2)/\sqrt{2\pi}
\]
for simplicity, and
\[
h_{n}(x)= H_{n}(x)/\sqrt{n!}.
\]
$h_{n}(x)$ is the orthonormal modification of the so-called probabilistic
Hermite polynomials $\left\{  H_{n}(x)\right\}  $, i.e., polynomials defined by
the following three-term recurrence:
\[
H_{n+1}(x)= xH_{n}(x)-nH_{n-1}(x),
\]
with $H_{-1}(x)=0$ and $H_{0}(x)
=1$. It is also known that
\begin{equation}
\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}h_{n}(x)h_{m}(x)\exp
(-x^{2}/2)dx=\delta_{nm}, \label{ortH}
\end{equation}
where $\delta_{nm}$ denotes Kronecker's delta.

In the second case we have
\begin{equation}
f_{G}(x,y|\rho)=f_{g}(x|\beta)f_{g}(y|\beta)\sum_{j\geq0}\rho^{n}l_{n}
(x|\beta)l_{n}(y|\beta), \label{MG}
\end{equation}
where $l_{n}(x|\beta)=\sqrt{\frac{n!}{\left(
\beta\right)  ^{\left(  n\right)  }}}L_{n}(x|\beta)$ is the orthonormal
modification of the so-called generalized Laguerre polynomials $\left\{
L_{n}(x|\beta)\right\}  $, defined by the following expansions:
\begin{equation}
L_{n}(x|\beta)=\sum_{k=0}^{n}(-1)^{k}\frac{(\beta)^{(n)}}{(n-k)!(\beta)^{(k)}
}x^{k}/k!. \label{defL}
\end{equation}


Let us remark that the expansion (\ref{MG}) is known under the name
Hardy-Hille formula (see, e.g., \cite{Szeg75}, p.~102). One knows that
\begin{equation}
\int_{0}^{\infty}l_{n}(x|\beta)l_{m}(x|\beta)f_{g}(x|\beta)dx=\delta_{n,m}.
\label{ortL}
\end{equation}
Our aim is to calculate $\left\{  E(X-Y)^{n}\right\}  _{n\geq0}$.  We will do
it in two ways.

The first way is to calculate the generating function of the set of these
numbers, i.e., the function
\[
g(t)=\sum_{n\geq0}t^{n}E(X-Y)^{n}/n!.
\]
Since it is known that all moments exist,
it can be calculated by exchanging integration and summation. Namely, we have
\[
g(t)=\int\int\exp\left(  tx-ty\right)  f(x,y|\rho)dxdy.
\]
Obviously here and below the integral depends on the case of
integration over all of $\mathbb{R}^{2}$,
or just over $\mathbb{R}^{+}\times\mathbb{R}^{+}$.

The second way is to calculate the number using an expansion:
\[
E(-X+Y)^{n}=\sum_{m=0}^{n}(-1)^{m}EX^{m}Y^{n-m}.
\]
Now
\[
EX^{m}Y^{n-m}=\int\int x^{m}y^{n-m}f(x,y|\rho)dxdy.
\]
where $f$ is either $f_{N}$ or $f_{G}$ presented above. Using one of the
expansions (\ref{MN}) or (\ref{MG}) we will find these moments using the
numbers
\[
H_{m,j}=EX^{m}k_{j}(X),
\]
where $k_{i}(x)$ is either $h_{j}(x)$, if we consider the Normal case, or
$l_{j}(x)$, if we consider the Gamma case. Since both expansions (\ref{MN}) and
(\ref{MG}) have a similar structure, we have
\begin{equation}
EX^{m}Y^{n-m}=\sum_{j=0}^{\min(m,n-m)}\rho^{j}H_{m,j}H_{n-m,j}.\label{ExpH}
\end{equation}
This is so since $X$ and $Y$ have the same distributions and since
$H_{m,j}=0$ if $j>m$ because of the orthogonality of
polynomials $h$ or $l$. We will calculate this function assuming either
expansion (\ref{MN}) or (\ref{MG}).

\label{Gaussian}So let us start with the expansion (\ref{MN}). First let us
calculate the auxiliary numbers $H_{j,n}$. We have
\begin{equation}
H_{j,n}=EX^{j}h_{n}(X)=\frac{1}{\sqrt{n!}}EX^{j}H_{n}(X)=
\begin{cases}
0, & \text{if }n>j\text{ or }j-n\text{ odd;}\\
\frac{j!}{2^{(j-n)/2}((j-n)/2)!\sqrt{n!}}, & \text{if }j-n\text{ is even.}
\end{cases}
\label{HjnN}
\end{equation}
This is so, since polynomials $h_{n}$ are orthogonal and also since it is
common knowledge that 
\[
x^{j}= j!\sum_{m=0}^{\left\lfloor j/2\right\rfloor
}\frac{1}{2^{m}m!(j-2m)!}H_{j-2m}(x)
\]
for all $n\geq0$.
Secondly, let us calculate the following auxiliary function that will simplify
many further calculations.
\begin{align*}
m_{n}(t) &  = Eh_{n}(X)\exp(tX)
=\frac{1}{\sqrt{2\pi n!}}\int_{-\infty}^{\infty}H_{n}
(x)\exp(tx)\exp(-x^{2}/2)dx\\
&  =\frac{\exp(t^{2}/2)}{\sqrt{2\pi n!}}\int_{-\infty}^{\infty}H_{n}
(x)\exp(-(x-t)^{2}/2)dx.
\end{align*}
Now let us change the variable under the integral by setting $y
=x-t$. Then we get
\[
m_{n}(t)=\frac{\exp(t^{2}/2)}{\sqrt{2\pi n!}}
\int_{-\infty}^{\infty}H_{n}(y+t)\exp(-y^{2}/2)dy.
\]
Next, we utilize the following the expansion:
\[
H_{n}(x+y)=\sum_{j=0}^{n}\binom{n}{j}H_{j}(x)y^{n-j}.
\]
Now, since we have (\ref{ortH}), we see that $m_{n}(t)=
t^{n}\frac{\exp(t^{2}/2)}{\sqrt{n!}}$ and we get
\begin{align*}
E\exp\left(  tX-tY\right)  &=\sum_{n\geq0}E\exp\left(
tX-tY\right) \\
&=\frac{1}{2\pi}\sum_{n\geq0}\rho^{n}\int_{-\infty}^{\infty}\int_{-\infty
}^{\infty}\exp\left(  tx-ty\right)  \exp\left(  -\left(  x^{2}+y^{2}\right)
/2\right)  h_{n}(x)h_{n}(y)dxdy 
\end{align*}
and
$$
\sum_{n\geq0}\rho^{n}m_{n}(t)m_{n}(-t)=\exp(t^{2})\sum_{n\geq0}(-\rho
)^{n}t^{2n}/n!=\exp(t^{2}(1-\rho)) \sum_{n\geq0}t^{2n}(1-\rho)^{n}/n!.
$$
Thus, we deduce that
\begin{equation}
E(X-Y)^{n}=
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{n!}{(n/2)!}(1-\rho)^{n/2}, & \text{if }n\text{ is even.}
\end{cases}
\label{nG}
\end{equation}

\label{Gamma}Now let us consider expansion (\ref{MG}). We need to recall some
simple facts (see any textbook on probability). The Gamma
distribution with rate parameter zero and shape parameter $\beta>0$ is the
distribution with the density $f_{g}$ presented above. It is easy to see,
recalling the definition of the Euler's Gamma function that
\begin{equation}
EX^{n}=\int_{0}^{\infty}x^{n}f_{g}(x|\beta)dx=\left(  \beta\right)  ^{\left(
n\right)  }. \label{momG}
\end{equation}
As before, let us calculate the set of auxiliary quantities and functions. We
start with the numbers
\[
H_{j,n}=EX^{j}l_{n}(X).
\]
We have
\begin{equation}
H_{j,n}=\frac{\sqrt{n!}}{\Gamma(\beta)\sqrt{\left(  \beta\right)  ^{\left(
n\right)  }}}\int_{0}^{\infty}x^{j}L_{n}^{(\beta)}(x)x^{\beta-1}\exp(-x)dx.
\label{defH}
\end{equation}
Now, we use the following expansion of $x^{j}$ in powers of Laguerre
polynomials (see \cite{Chih79,Szeg75}):
\begin{equation}
x^{j}= j!\sum_{k=0}^{j}(-1)^{k}\frac{(\beta)^{(j)}
}{(j-k)!(\beta)^{(k)}}L_{k}^{(\beta)}(x), \label{xL}
\end{equation}
and then we use the orthogonality of the Laguerre polynomials. Hence, we have
\begin{equation}
H_{j,n}=(-1)^{n}\binom{j}{n}(\beta)^{(j)}\sqrt{\frac{n!}{\left(  \beta\right)
^{\left(  n\right)  }}}. \label{HjnG}
\end{equation}


Let us now calculate the following auxiliary functions:
\begin{align*}
m_{n}(t) &  =El_{n}(X)\exp(tX)=\sqrt{\frac{n!}{\left(  \beta\right)  ^{\left(
n\right)  }}}E\exp(tX)L_{n}(X|\beta)\\
&  =\sqrt{\frac{n!}{\left(  \beta\right)  ^{\left(  n\right)  }}
}\frac{1}{\Gamma(\beta)}\int_{0}^{\infty}\exp(xt)L_{n}(x|\beta)x^{\beta-1}
\exp(-x)dx\\
&  =\sqrt{\frac{n!}{\left(  \beta\right)  ^{\left(  n\right)  }}
}\frac{1}{\Gamma(\beta)}\int_{0}^{\infty}L_{n}(x|\beta)x^{\beta-1}
\exp(-x(1-t))dx.
\end{align*}
Now, let us change variables under the integral by considering $y
= x(1-t)$. Then we get 
\[
m_{n}(t)=\sqrt{\frac{n!}{\left(  \beta\right)  ^{\left(  n\right)  }}}\frac
{1}{1-t}\frac{1}{\Gamma(\beta)}\int_{0}^{\infty}L_{n}(\frac{y}{1-t}
|\beta)\left(  \frac{y}{1-t}\right)  ^{\beta-1}\exp(-y)dy.
\]
We will now utilize the following recurrence relation for Laguerre polynomials
(see, e.g., \cite[Chap.~22]{Abram64}):
\[
L_{n}(y|\beta)=\sum_{j=0}^{n}L_{n-j}(x|\beta+j)(y-x)^{j}/j!.
\]
So we get
$$
m_{n}(t)=\sqrt{\frac{n!}{\left(  \beta\right)  ^{\left(  n\right)
}}}\frac{1}{(1-t)^{\beta}}\frac{1}{\Gamma(\beta)}\times\\
\int_{0}^{\infty}\left(  \sum_{j=0}^{n}\left(  \frac{-ty}{1-t}\right)
^{j}L_{n-j}(x|(\beta+j))/j!\right)  y^{\beta-1}\exp(-y)dy.
$$
Further, we have
\begin{align*}
m_{n}(t)&=\sqrt{\frac{n!}{\left(  \beta\right)  ^{\left(
n\right)  }}}\frac{1}{(1-t)^{\beta}}\frac{1}{\Gamma(\beta)}\times
\sum_{j=0}^{n}\int_{0}^{\infty}(-1)^{j}\frac{t^{j}}{(1-t)^{j}}L_{n-j}
(y\left\vert \left(  \beta+j\right)  \right\vert )y^{j+\beta-1}\exp
(-y)dy\\
&=(-1)^{n}\sqrt{\frac{n!}{\left(  \beta\right)  ^{\left(  n\right)
}}}\frac{t^{n}\Gamma(\beta+n)}{\Gamma(\beta)n!(1-t)^{\beta+n}} \\
&=(-1)^{n}\sqrt{\frac{n!}{\left(  \beta\right)  ^{\left(  n\right)
}}}\frac{t^{n}\left(  \beta\right)  ^{\left(  n\right)  }}{n!(1-t)^{\beta+n}
}.
\end{align*}
Hence, we now have
\begin{align*}
 g(t)&=E\exp(tX-tY)=\sum_{n\geq0}\rho^{n}\int_{0}^{\infty}
\exp(tx-ty)l_{n}(x|\beta)l_{n}(y|\beta)f_{g}(x|\beta)f_{g}(y|\beta
)dxdy\\
&=\sum_{n\geq0}\rho^{n}m_{n}(t)m_{n}(-t)=
\sum_{n\geq0}\rho^{n}\frac{n!(-1)^{n}t^{2n}\left(  \beta\right)  ^{\left(
n\right)  }\left(  \beta\right)  ^{\left(  n\right)  }}{\left(  \beta\right)
^{\left(  n\right)  }n!n!(1-t^{2})^{n+\beta}}\\
&=\frac{1}{(1-t^{2})^{\beta}}\sum_{n\geq0}\rho^{n}\frac{(-1)^{n}t^{2n}\left(
\beta\right)  ^{\left(  n\right)  }}{n!(1-t^{2})^{n}}\\
&=\frac{1}{(1-t^{2})^{\beta}}\sum_{n\geq0}(-1)^{n}\left(  \frac{\rho t^{2}
}{1-t^{2}}\right)  ^{n}\left(  \beta\right)  ^{\left(  n\right)  }/n!,
\end{align*}
and further we have
\begin{align*}
g(t)&=\frac{1}{(1-t^{2})^{\beta}}\left(  1+\frac{\rho
t^{2}}{1-t^{2}}\right)  ^{-\beta}\\
&=\frac{1}{(1-t^{2}(1-\rho))^{\beta}} \\
&=\sum_{n\geq0}t^{2n}
(1-\rho)^{n}\left(  \beta\right)  ^{\left(  n\right)  }/n!.
\end{align*}
Here we twice used the formula for the binomial series. Thus, we
deduce
\begin{equation}
E(X-Y)^{n}=
\begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{n!}{(n/2)!}\left(  \beta\right)  ^{\left(  n/2\right)  }(1-\rho
)^{n/2}, & \text{if }n\text{ is even.}
\end{cases}
\label{EXYG}
\end{equation}

\subsection{Proofs}

The basis of all proofs is the consideration of different methods for
calculating $\left\{  E(X-Y)^{n}\right\}  _{n\geq1}$.

\begin{proof}[Proof of assertion (i) of Theorem \ref{Main}]
Let $(X,Y)$ have bivariate Gaussian distribution given by (\ref{2xN}) and let
$N(m,\sigma^{2})$ denote Normal (or Gaussian) distribution with mean value $m$
and variance $\sigma^{2}$. Then it is well known that every linear
transformation of such a variable also has Gaussian distribution. In particular
we deduce that $X-Y\sim N(0,(1-\rho)^{2})$. Hence, keeping in mind that
\[
EZ^{j}=
\begin{cases}
0, & \text{if }j\text{ is odd;}\\
1, & \text{if }j=0;\\
\sigma^{j}(2j-1)!!, & \text{if }j\text{ is even,}
\end{cases}
\]
if $Z\sim N(0,\sigma^{2})$. Thus, the left-hand side of assertion (i) is
equal to $(2k-1)!!(1-\rho)^{k}=\frac{(2k)!}{k!}
(1-\rho)^{k}$. In order to get the right-hand side, we calculate
\begin{align*}
E(X-Y)^{2k}&= E(X-\rho Y-(1-\rho)Y)^{2k} \\
&= \sum_{j=0}^{2k}\binom{2k}{j}(-1)^{j}(1-\rho)^{2k-j}E(X-\rho
Y)^{j}Y^{2k-j}\\
&= \sum_{j=0}^{2k}\binom{2k}{j}(-1)^{j}(1-\rho)^{2k-j}E(E(X-\rho
Y)^{j}|Y)Y^{2k-j}.
\end{align*}
Now, we use the known fact that random variables $(X-\rho Y)$ and $Y$ are
independent. Hence we get
\[
(E(X-\rho Y)^{j}|Y)=
\begin{cases}
0, & \text{if }j\text{ is odd;}\\
1, & \text{if }j=0\text{;}\\
(2j-1)!!(1-\rho^{2})^{j/2}, & \text{if }j\text{ is even.}
\end{cases}
\]
Thus, we get
\[
E(X-Y)^{2k}=\sum_{j=0}^{k}\binom{2k}{2j}(1-\rho)^{2k-2j}(1-\rho^{2}
)^{j}(2j-1)!!(2k-2j-1)!!.
\]
Now it is enough to change $\rho$ to $-\rho$ to get the assertion.
\end{proof}

\begin{proof}[Proof of assertion (ii) of Theorem \ref{Main}]

We use (\ref{ExpH}) with (\ref{HjnN}), getting
\[
\frac{(2k)!}{k!}(1-\rho)^{k}=2\sum_{j=0}^{k-1}(-1)^{j}\binom{2k}{j}\sum
_{n=0}^{j}\frac{\rho^{n}}{n!}H_{j,n}H_{2k-j,n}+(-1)^{k}\binom{2k}{k}\sum
_{n=0}^{k}\frac{\rho^{n}}{n!}H_{k,n}^{2},
\]
where
\[
H_{j,n}=
\begin{cases}
0, & \text{if }n>j\text{ or }j-n\text{ odd;}\\
\frac{j!}{2^{(j-n)/2}((j-n)/2)!}, & \text{if }j-n\text{ is even.}
\end{cases}
\]
Notice that, if $j-n$ is even, then 
$$\binom{2k}{j}H_{j,n}H_{2k-j,n}
=2^{n}\binom{2k}{k}\binom{k}{n}\binom{k-n}
{(j-n)/2}/2^{k},$$ and when $k-n$ is even, we have 
$$H_{k,n}^{2}= 2^{n}\binom{k}{n}\binom{k-n}{(k-n)/2}/2^{k}.$$
Hence (\ref{G2}) can be
reformulated in the following way:
\[
\frac{1}{2^{k-1}}\sum_{j=0}^{k-1}(-1)^{j}\sum_{n=0-}^{j}(2\rho)^{n}\binom
{k}{n}mh_{j,n,k}+\frac{(-1)^{k}}{2^{k}}\sum_{n=0}^{k}(2\rho)^{n}\binom{k}
{n}mh_{k,n,k},=(1-\rho)^{k},
\]
where
\[
mh_{j,n,k}=
\begin{cases}
0, & \text{if }j>k-1\text{ or }j-n\text{ is odd;}\\
\binom{k-n}{(j-n)/2}, & \text{if }j-n\text{ is even.}
\end{cases}
\]
Now it remains to change the order of summation in the internal sums and
$\rho$ to $-\rho$.
\end{proof}

\begin{proof}[Proof of assertion (iii) of Theorem \ref{Main}]
Let us recall Eqs.~(\ref{defL}) and (\ref{defH}) and let us calculate $H_{j,n}$
directly, getting
\begin{align*}
H_{j,n}&=\frac{\sqrt{n!}}{\Gamma(\beta)\sqrt{\Gamma(n+\beta)}}\sum_{k=0}
^{n}\frac{(-1)^{k}}{k!}\frac{(\beta)^{(n)}}{(n-k)!(\beta)^{(k)}}\int
_{0}^{\infty}x^{k+j}x^{\beta-1}\exp(-x)dx\\
&=\frac{\sqrt{n!}}{\sqrt{\Gamma(n+\beta)}}\sum_{k=0}^{n}\frac{(-1)^{k}}
{k!}\frac{(\beta)^{(n)}(\beta)^{(k+j)}}{(n-k)!(\beta)^{(k)}}.
\end{align*}
Now, we compare it with the already-calculated $H_{j,n}$ for the Gamma
distribution, i.e., $(-1)^{n}\binom{j}{n}(\beta)^{(j)}\sqrt{\frac{n!}{\left(
\beta\right)  ^{\left(  n\right)  }}}$.
\end{proof}

\begin{proof}[Proof of assertion (iv) of Theorem \ref{Main}]

Recall Eq.~(\ref{HjnG}). We have
\begin{align*}
E(-X+Y)^{n} &=\sum_{m=0}^{n}(-1)^{m}\binom{n}{m}EX^{m}Y^{n-m} \\
&= \begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{n!}{(n/2)!}\left(  \beta\right)  ^{\left(  n/2\right)  }(1-\rho
)^{n/2}, & \text{if }n\text{ is even.}
\end{cases}
\end{align*}
Now we use (\ref{ExpH}) with $H_{j,m}$ given by (\ref{HjnG}). By (\ref{EXYG})
we have
\begin{align*}
E(-X+Y)^{n}&=\sum_{m=0}^{n}(-1)^{m}\binom{n}{m}\sum_{j=0}^{\min(m,n-m)}\rho
^{j}H_{m,j}H_{n-m,j}\\
&=\sum_{m=0}^{n}(-1)^{n-m}\binom{n}{m}\sum_{j=0}^{\min(m,n-m)}\rho^{j}\binom
{m}{j}\binom{n-m}{j}\left(  \beta\right)  ^{\left(  m\right)  }\left(
\beta\right)  ^{\left(  n-m\right)  }\frac{j!}{\left(  \beta\right)  ^{\left(
j\right)  }}\\
&=\sum_{m=0}^{n}(-1)^{n-m}\binom{n}{m}\left(  \beta\right)  ^{\left(  m\right)
}\left(  \beta\right)  ^{\left(  n-m\right)  }\sum_{j=0}^{m}\binom{m}{j}
\binom{n-m}{j}\frac{j!\rho^{j}}{\left(  \beta\right)  ^{\left(  j\right)  }}\\
&= \begin{cases}
0, & \text{if }n\text{ is odd;}\\
\frac{n!}{(n/2)!}\left(  \beta\right)  ^{\left(  n/2\right)  }(1-\rho
)^{n/2}, & \text{if }n\text{ is even.}
\end{cases}
\end{align*}
This is so, since $\binom{m}{j}\binom{n-m}{j}$ is zero whenever $m>j$ or
$n-m>j$.
\end{proof}

\begin{proof}[Proof of assertion (v) of Theorem \ref{Main}]

One can easily notice, based on (\ref{MG}), that $E(l_{n}(X)|Y)
=\rho^{n}l_{n}(Y)$. Hence, let us first calculate the conditional
moments $\eta_{n}(y|\beta,\rho)= E(X^{n}|Y=y)$. 
By (\ref{xL}) we have
\begin{align*}
\eta_{j}(y|\beta,\rho) &  =E(X^{j}|Y=y)= j!\sum
_{k=0}^{j}(-1)^{k}\frac{(\beta)^{(j)}}{(j-k)!(\beta)^{(k)}}\rho^{k}
L_{k}^{(\beta)}(y)\\
&  =j!\sum_{k=0}^{j}(-1)^{k}\frac{(\beta)^{(j)}}{(j-k)!(\beta)^{(k)}}\rho
^{k}\sum_{m=0}^{k}(-1)^{m}\frac{(\beta)^{(k)}}{(k-m)!(\beta)^{(m)}}y^{m}/m!\\
&  =j!\sum_{m=0}^{j}\frac{(\beta)^{(j)}}{(\beta)^{(m)}(j-m)!m!}(\rho
y)^{m}\sum_{k=m}^{j}(-1)^{k-m}\frac{(j-m)!}{(k-m)!(j-k)!}\rho^{k-m}\\
&  =\sum_{m=0}^{j}\binom{j}{m}\frac{(\beta)^{(j)}}{(\beta)^{(m)}}(\rho
y)^{m}(1-\rho)^{j-m}=j!(1-\rho)^{j}L_{j}(-\frac{\rho y}{1-\rho}).
\end{align*}
Further, we get
\begin{align*}
E\left(  (X-Y)^{n}|Y=y\right)  &=\sum_{j=0}^{n}(-1)^{n-j}\binom{n}{j}\eta
_{j}(y|\beta,\rho)y^{n-j}\\
&=\sum_{j=0}^{n}(-1)^{n-j}y^{n-j}\binom{n}{j}\sum_{m=0}^{j}\binom{j}{m}
\frac{(\beta)^{(j)}}{(\beta)^{(m)}}(\rho y)^{m}(1-\rho)^{j-m}
\end{align*}
\begin{align*}
&=\sum_{m=0}^{n}\frac{n!}{m!(n-m)!}\frac{(\rho y)^{m}}{(\beta)^{(m)}}\sum
_{j=m}^{n}(-1)^{n-j}\frac{(n-m)!(\beta)^{(j)}y^{n-j}(1-\rho)^{j-m}
}{(n-j)!(j-m)!}\\
&\quad\quad\quad \sum_{m=0}^{n}\frac{n!}{m!(n-m)!}(\rho y)^{m}\sum_{s=0}^{n-m}(-1)^{n-m-s}
\binom{n-m}{s}y^{n-m-s}(1-\rho)^{s}(\beta+m)^{(s)}\\
&=\sum_{t=0}^{n}(-1)^{n-t}y^{n-t}\frac{n!}{t!(n-t)!}\sum_{s=0}^{t}\binom{t}
{s}(\rho y)^{t-s}(1-\rho)^{s}(\beta+t-s)^{(s)} \\
&=\sum_{s=0}^{n}\binom{n}{s}y^{n-s}(1-\rho)^{s}\sum_{t=s}^{n}(-1)^{n-t}
\frac{(n-s)!}{(t-s)!t!}\rho^{t-s}(\beta+t-s)^{(s)}\\
&=\sum_{s=0}^{n}\binom{n}{s}y^{n-s}(1-\rho)^{s}\sum_{m=0}^{n-s}(-1)^{n-m-s}
\binom{n-s}{m}\rho^{m}(\beta+m)^{(s)}.
\end{align*}
Hence, using (\ref{momG}), we get 
\begin{align*}
E(X-Y)^{n}&=\int_{0}^{\infty}(\sum_{s=0}^{n}\binom{n}{s}y^{n-s}(1-\rho
)^{s}\times 
\sum_{m=0}^{n-s}(-1)^{n-m-s}\binom{n-s}{m}\rho^{m}(\beta+m)^{(s)}
)f_{g}(x|\beta)dx\\
&=\sum_{s=0}^{n}\binom{n}{s}(\beta)^{(n-s)}(1-\rho)^{s}\sum_{m=0}
^{n-s}(-1)^{n-m-s}\binom{n-s}{m}\rho^{m}(\beta+m)^{(s)}\\
&=\sum_{t=0}^{n}(-1)^{n-t}\frac{n!(\beta)^{(n-t)}}{t!(n-t)!}\sum_{s=0}^{t}
\frac{t!}{s!(t-s)!}(1-\rho)^{s}\rho^{t-s}(\beta+n-t)^{(t-s)}(\beta+t-s)^{(s)}.
\end{align*}
In the last line, we changed the order of summation. To get the
assertion, we use the formula
\[
\left(  \beta\right)  ^{\left(  n\right)  }\left(  \beta+n\right)  ^{\left(
m\right)  }=\left(  \beta\right)  ^{\left(  n+m\right)  },
\]
true for all complex $\beta$.
\end{proof}

\begin{thebibliography}{99}                                                                                               


\bibitem {Alexits61}G. Alexits, \textit{Convergence Problems of Orthogonal
Series}, in International Series of
Monographs in Pure and Applied Mathematics, Vol. 20, Pergamon Press, 1961.

\bibitem {Abram64}M. Abramowitz and I. A. Stegun, eds.,
\textit{Handbook of Mathematical Functions with Formulas, Graphs, and
Mathematical Tables}, 
Applied Mathematics Series, Vol. 55, Dover, 1983.

\bibitem {Ball03}K. Ball, \textit{Strange Curves,
Counting Rabbits, and Other Mathematical Explorations},
Princeton University Press, 2003.  

\bibitem {Cull05}P. Cull,  M. Flahive, and R. Robson, \textit{Difference
Equations. From Rabbits to Chaos.} 
Springer, 2005.

\bibitem {Chih79}T. S. Chihara, \textit{An Introduction to Orthogonal
Polynomials}, Mathematics and its Applications, Vol. 13,  Gordon and
Breach, 1978.

\bibitem {Grad14}I. S. Gradshteyn and I. M. Ryzhik, \textit{Table of
Integrals, Series, and Products}, 8th edition,  Academic Press, 2014.

\bibitem {Knuth98}R. L. Graham, D. E. Knuth, and O. Patashnik, \textit{Concrete
Mathematics}, 2nd ed., Addison-Wesley, 1994. 

\bibitem {Lancaster58}H. O. Lancaster, The structure of bivariate
distributions, \emph{Ann. Math. Statistics}, \textbf{29} (1958),
719--736.

\bibitem {Lancaster75}H. O. Lancaster, Joint probability distributions in the
Meixner classes. \emph{J. Roy. Statist. Soc. Ser. B }\textbf{37} (1975),
434--443. 

\bibitem {SzabKes}P. J. Szab\l owski, Yet another way of calculating moments of
the Kesten's distribution and its consequences for Catalan numbers and Catalan
triangles, \emph{Discrete Math.} \textbf{345} (2022), Paper No.\ 112891.

\bibitem {SzabStac}P. J. Szab\l owski, On stationary Markov processes with
polynomial conditional moments. \emph{Stoch. Anal. Appl.} \textbf{35} (2017),
852--872.

\bibitem {Szabl21}P. J. Szab\l owski, On positivity of orthogonal series and
its applications in probability, \emph{Positivity} \textbf{26} (2022),
Article 19. 

\bibitem {Szeg75}G. Szeg\H{o}, {\it Orthogonal Polynomials}, 4th edition,
\emph{Amer. Math. Soc. Colloq. Publ.}, Vol. \textbf{23}, Amer. Math. Soc.,
1975.
\end{thebibliography}
\bigskip
\hrule
\bigskip

\noindent 2020 {\it Mathematics Subject Classification}:
Primary 05A10; Secondary 05D40, Secondary 11B65, 11B73.

\noindent \emph{Keywords: } binomial coefficient, Pochhammer symbol, bivariate normal distribution,
bivariate gamma distribution, Hermite polynomial, Laguerre polynomial,
Fibonacci number, Lucas number.

\bigskip
\hrule
\bigskip

\noindent (Concerned with sequences
\seqnum{A000032},
\seqnum{A000045},
\seqnum{A001147},
\seqnum{A007318}, and
\seqnum{A221954}.)

\bigskip
\hrule
\bigskip

\vspace*{+.1in}
\noindent
Received July 7 2023; 
revised versions received  November 7 2023; November 13 2023.
Published in {\it Journal of Integer Sequences}, January 3 2024.

\bigskip
\hrule
\bigskip

\noindent
Return to
\htmladdnormallink{Journal of Integer Sequences home page}{https://cs.uwaterloo.ca/journals/JIS/}.
\vskip .1in


\end{document}
