# Fluctuation theory of continuous-time skip-free downward Markov chains with
applications to branching processes with immigration

R. Loeffen  University of Manchester, School of Mathematics, Manchester, M13
9PL, UK<EMAIL_ADDRESS>,  P. Patie  School of Operations
Research and Information Engineering, Cornell University, Ithaca, NY 14853.
<EMAIL_ADDRESS>and  J. Wang  School of Operations Research and
Information Engineering, Cornell University, Ithaca, NY 14853.
<EMAIL_ADDRESS>
###### Abstract.

We develop a new methodology for the fluctuation theory of continuous-time
skip-free Markov chains, extending the recent work of Choi and Patie [6] for
discrete-time skip-free Markov chains. As the main application we use it to
derive a full set of fluctuation identities regarding exiting a finite or
infinite interval for Markov branching processes with immigration, thereby
uncovering many new results for this classical family of continuous-time
Markov chains. The theory also allows us to recover in a simple manner
fluctuation identities for skip-free downward compound Poisson processes.

###### Key words and phrases:

Keywords: Skip-free Markov chains, First passage time problems, fluctuation
theory, Branching processes with immigration, Potential theory

##  1\. Introduction

First passage times are ubiquitous in probability theory and, more broadly, in
many fields of sciences, such as epidemiology, physics, neurology, insurance
and financial mathematics. In spite of the wide applicability of first-passage
phenomena, finer characterizations of the law of these random variables beyond
as solutions of a boundary value problem associated with a linear operator,
have been obtained only for very few and specific classes of Markov processes,
including most notably one-dimensional diffusion processes and their discrete-
state space analogue birth-death chains on the one hand and processes with
stationary and independent increments on the other hand, both for which
tailor-made approaches have been developed. Unfortunately, attempts to carry
over techniques from these cases have had limited success, providing some
evidence towards the need to find a novel way to tackle this issue. In this
vein, in the recent paper [6], Choi and Patie develop an original and
comprehensive approach based on a combination of techniques such as the theory
of Martin boundary and potential theory to provide an expression for the
resolvent of a discrete-time skip-free Markov chain in terms of the so-called
fundamental excessive functions of three different chains, namely the given
skip-free downward (upward) Markov chain, its dual and the original one killed
when going above (below) a given level. This allowed them to derive various
fluctuation identities involving the first hitting time of a given point and
the first exit times out of a given finite or semi-infinite interval.

In this paper we will follow the approach of Choi and Patie [6] to extend
their results to the case of continuous-time skip-free Markov chains, while
also improving on their methodology as explained in the beginning of Section
3. This is our first main contribution. Besides our time homogeneous chain
being skip-free (downward) we further assume that all its states are stable
(i.e. no instantaneous states) and transient, though we allow for state-
dependent killing and for explosion with the precise setting being laid out in
Section 2. Though we already mention that for our applications covered in
Sections 4 and 5 we can also deal with the cases where not all states are
transient by employing the usual trick of including a constant killing rate,
say $q>0$, which makes the chain transient, and then letting $q\downarrow 0$.
Our second main contribution is to use the general theory presented in Section
3 to derive a collection of fluctuation identities for Markov branching
processes with immigration (MBIs), a well-studied class of continuous-time
Markov chains. Most of these identities are new as explained at the end of
Section 5, which shows the power of our methodology. As a quick illustration
of the general theory and in order to contrast it with the case of MBIs we
will also derive the same type of fluctuation identities for continuous-time
skip-free Markov chains with stationary and independent increments, i.e. skip-
free compound Poisson processes (CPPs), though for this class the resulting
identities are (essentially) known.

Of course for deriving certain probabilities involving first passage times
there is the analytical approach. To illustrate this, let $p(x)$ be the
probability that our skip-free downward continuous-time Markov chain with
state space $E\subset\mathbb{Z}$ and starting at $x\in E$ exits a given
interval $(a,b)\subset E$ by hitting $a$. Then, with some of the terminology
explained later on in Section 2, $p(x),x\in[a-1,b+1]$, is harmonic for the
Markov chain killed at exiting $(a,b)$ and as a consequence, with $Q$ the
Q-matrix of the chain, $p(x)$, $x\geq a$, solves $\sum_{y\in E}Q(x,y)p(y)$,
$x\in[a+1,b-1]$, with the boundary conditions $p(a)=1$ and $p(x)=0$ for $x\geq
b$. Consequently, $p(x)$ can be computed by solving a linear system of $b-a-1$
equations. Note that this approach becomes more difficult when $E$ is
unbounded and $a=-\infty$ or $b=\infty$ as one then has to deal with an
infinite system of equations. Our main theorems in Section 3 provide in
particular a structure for the solution of these systems. Besides yielding
more insight into the solution, even in the case of a finite system, our
results for MBIs and CPPs allow one to compute probabilities like $p(x)$ for
these classes in a more efficient way than by numerically solving a system of
linear equations, especially when one wants to compute $p(x)$ for various
values of $a$ and $b$.

The rest of the paper is organized as follows. In Section 2 we introduce
various notations and concepts that will be needed later on and in particular
define the class of Markov chains that we will be working with. Section 3
contains the main results where we, following Choi and Patie [6], express, for
a general Markov chain, the resolvent and derive various fluctuation
identities in terms of the three fundamental excessive functions, whereas in
Sections 4 and 5 we apply these results to the class of skip-free CPPs and
MBIs respectively.

##  2\. Preliminaries

Let $E=[\mathfrak{l},\mathfrak{h}]\cap\mathbb{Z}$ where
$\mathfrak{l},\mathfrak{h}\in\mathbb{Z}\cup\\{-\infty,\infty\\}$ with
$\mathfrak{l}<\mathfrak{h}$. We set $E_{\partial}=E\cup\\{\partial\\}$ where
$\partial$ denotes the cemetery state. We always work in a canonical setup
where (i) $\Omega$ is the space of functions $\omega:[0,\infty)\to
E_{\partial}$ with the property that there exists $\zeta(\omega)\in[0,\infty]$
such that $\omega$ is $E$-valued and right-continuous on $[0,\zeta(\omega))$
and $\omega(t)=\partial$ for $t\geq\zeta(\omega)$, (ii) $X=(X_{t})_{t\geq 0}$
is the process defined by $X_{t}(\omega)=\omega(t)$ for $\omega\in\Omega$,
(iii) $(\mathcal{F}_{t})_{t\geq 0}$ is the natural filtration of $X$ and (iv)
$\mathcal{F}$ is the smallest sigma-algebra containing $\cup_{t\geq
0}\mathcal{F}_{t}$. The random variable $\zeta$ then represents the lifetime
$\inf\\{t\geq 0:X_{t}=\partial\\}$ of $X$. With $J_{0}=0$ and
$J_{n}=\inf\\{t\geq J_{n-1}:X_{t}\neq X_{J_{n-1}}\\}$ for $n\geq 1$ , the
random variables $J_{1},J_{2},\ldots$ are the consecutive jump times of $X$.
We denote by $J_{\infty}=\lim_{n\to\infty}J_{n}$ the explosion time.

Let $Q=(Q(x,y))_{x,y\in E}$ be a given $Q$-matrix on $E$, i.e. a matrix
satisfying $0\leq Q(x,y)<\infty$ for $x\neq y$ and $\sum_{y\neq x}Q(x,y)\leq
Q(x)<\infty$ where $Q(x):=-Q(x,x)$. We set $Q(x,\partial)=Q(x)-\sum_{y\neq
x}Q(x,y)\geq 0$ for $x\in E$. Note that $Q$ is not assumed to be conservative,
i.e. $Q(x,\partial)$ is not necessarily equal to $0$. Given $Q$ we next
specify a continuous-time, time homogeneous Markov chain on
$(\Omega,\mathcal{F})$ via the usual jump chain/holding time definition.
Namely, we assume $\mathbb{P}=(\mathbb{P}_{x})_{x\in E_{\partial}}$ is a
family of probability measures on $(\Omega,\mathcal{F})$ such that (i)
$\mathbb{P}_{x}(X_{0}=x)=1$ for $x\in E_{\partial}$, (ii) for all $n\geq 1$
and $x\in E$, given $J_{n-1}<\infty$ and $X_{J_{n-1}}=x$ the $n$-th holding
time $H_{n}:=J_{n}-J_{n-1}$ is exponentially distributed with parameter $Q(x)$
if $Q(x)>0$ and equal to $\infty$ otherwise and is further conditionally
independent from $\mathcal{F}_{J_{n-1}}$ (iii) for all $n\geq 1$, $x\in E$ and
$y\in E_{\partial}$, given $J_{n-1}<\infty$ and $X_{J_{n-1}}=x$, $X_{J_{n}}=y$
with probability $\frac{Q(x,y)}{Q(x)}$ if $Q(x)>0$ and $X_{J_{n}}$ is
conditionally independent from $\mathcal{F}_{J_{n-1}}$ and $H_{n}$, (iv) $X$
is sent to the cemetery state if explosion occurs, i.e.
$X_{J_{\infty}}=\partial$ if $J_{\infty}<\infty$. We call $(X,\mathbb{P})$ the
Markov chain generated by $Q$ with state space $E$. Note that the Markov chain
can end up in the cemetery state either due to explosion or due to a
transition to $\partial$ by the jump chain $(X_{J_{n}})_{n\geq 0}$.

We write $P_{t}(x,y)=\mathbb{P}_{x}(X_{t}=y)$, $x,y\in E$, $t\geq 0$, for the
transition kernel of $(X,\mathbb{P})$. From the theory of continuous time
Markov chains (see e.g. Section 2.8 in [13]111In this reference it is assumed
that $Q$ is conservative but the arguments go through if $Q$ is not
conservative.), the transition kernel $P_{t}(x,y)$, $x,y\in E$, $t\geq 0$, is
the minimal nonnegative solution of both the forward equations,
$\frac{d}{dt}P_{t}(x,y)=\sum_{z\in E}P_{t}(x,z)Q(z,y)$,
$P_{0}(x,y)=\delta_{xy}$ with $\delta_{xy}$ the Kronecker delta, as well as
the backward equations, $\frac{d}{dt}P_{t}(x,y)=\sum_{z\in
E}Q(x,z)P_{t}(z,y)$, $P_{0}(x,y)=\delta_{xy}$. We will use the notation
$P_{t}f(x):=\sum_{y\in
E}P_{t}(x,y)f(y)=\mathbb{E}_{x}[f(X_{t})\mathds{1}_{\\{t<\zeta\\}}]$ and
$Qf(x):=\sum_{y\in E}Q(x,y)f(y)$ for a function $f:E\to\mathbb{R}$ and where
$\mathbb{E}_{x}$ is the expectation operator associated with $\mathbb{P}_{x}$.

Throughout the paper we assume the Q-matrix $Q$ is such that for $x,y\in E$,

(2.1) $Q(x,x-1)>0\quad\text{if $x>\mathfrak{l}$}\quad\text{and}\quad
Q(x,y)=0\quad\text{if $y\leq x-2$}.$

This assumption means that $(X,\mathbb{P})$ is skip-free downwards: from any
state $(X,\mathbb{P})$ can reach a lower state with positive probability but
only by visiting all states in between. We further assume, for the rest of
this section, that all states in $E$ are transient. The resolvent of
$(X,\mathbb{P})$ is defined by $G(x,y)=\int_{0}^{\infty}P_{t}(x,y)dt$. By the
transience assumption $G(x,y)<\infty$ for all $x,y\in E$ (see e.g. Theorem
3.4.2 in [13]) and by the skip-free downward assumption in combination with
the right-continuity of the sample paths $G(x,y)>0$ for all $x,y\in E$ such
that $y\leq x$.

We say that a function $h:E\to[0,\infty)$ is excessive, respectively harmonic,
for $(X,\mathbb{P})$, if $P_{t}h(x)\leq h(x)$, respectively $P_{t}h(x)=h(x)$,
for all $x\in E$ and $t\geq 0$. Note that by the semi-group property of
$P_{t}(x,y)$, $P_{t}h(x)$ is nonincreasing in $t$ when $h$ is excessive and so
$\lim_{t\downarrow 0}P_{t}h(x)=h(x)$ by the monotone convergence theorem and
the right-continuity of $P_{t}$ at $t=0$. Similarly, we call a measure
$\pi:E\to[0,\infty)$ excessive, respectively stationary, if $\sum_{x\in
E}\pi(x)P_{t}(x,y)\leq\pi(y)$, respectively $\sum_{x\in
E}\pi(x)P_{t}(x,y)=\pi(y)$, for all $y\in E$ and $t\geq 0$. Under our
assumptions there exists an excessive measure $\pi$ for $(X,\mathbb{P})$ that
is positive (i.e. $\pi(y)>0$ for all $y\in E$), e.g. $\pi(y)=\sum_{x\in
E}2^{-\sigma(x)}G(x,y)$ for some injective map $\sigma:E\to\\{1,2,\ldots\\}$
will do, see also Section 5.2 in [10]. We now fix such a positive, excessive
measure $\pi$ as our reference measure and we denote by
$\mathtt{g}(x,y)=G(x,y)/\pi(y)$ the resolvent density.

Next we introduce some auxiliary Markov chains associated with
$(X,\mathbb{P})$. First, writing
$\mathbb{\widehat{P}}=(\widehat{\mathbb{P}}_{y})_{y\in E_{\partial}}$, the
dual $(X,\mathbb{\widehat{P}})$ of $(X,\mathbb{P})$ with respect to $\pi$ is
the Markov chain with state space $E$ generated by the Q-matrix $\widehat{Q}$
given by

$\widehat{Q}(y,x)=\frac{\pi(x)Q(x,y)}{\pi(y)},\quad y,x\in E.$

Note that $\widehat{Q}$ is a well-defined Q-matrix since $\pi$ is a positive
measure and, for each $y\in E$,

$\begin{split}\pi(y)\sum_{x\in E}\widehat{Q}(y,x)=\sum_{x\in
E}\pi(x)Q(x,y)\leq&\liminf_{t\downarrow 0}\sum_{x\neq
y}\pi(x)\frac{P_{t}(x,y)}{t}+\pi(y)\lim_{t\downarrow
0}\frac{P_{t}(y,y)-1}{t}\\\ =&\liminf_{t\downarrow 0}\frac{\sum_{x\in
E}\pi(x)P_{t}(x,y)-\pi(y)}{t}\\\ \leq&0,\end{split}$

where the first inequality is due to Fatou’s lemma and, from the backward
equations, $\frac{d}{dt}P_{t}(x,y)|_{t=0}=Q(x,y)$ and the last inequality is
because $\pi$ is an excessive measure for $(X,\mathbb{P})$. Denoting by
$P_{t}(y,x):=\widehat{\mathbb{P}}_{y}(X_{t}=x)$ the transition kernel of
$(X,\mathbb{\widehat{P}})$, we easily see that
$\frac{\pi(x)P_{t}(x,y)}{\pi(y)}$, $y,x\in E$, $t\geq 0$, is the minimal
nonnegative solution of the backward (or forward) equations associated with
$\widehat{Q}$ and therefore
$\widehat{P}_{t}(y,x)=\frac{\pi(x)P_{t}(x,y)}{\pi(y)}$. Denoting by
$\widehat{G}$ and $\widehat{\mathtt{g}}$ respectively the resolvent and the
resolvent density (with respect to $\pi$) of the dual, we immediately observe
that $\widehat{\mathtt{g}}(y,x)=\mathtt{g}(x,y)$, $x,y\in E$. Now denote by
$T_{A}$ the first hitting time of the set $A\subseteq\mathbb{R}$, that is

$T_{A}=\inf\\{t\geq 0;X_{t}\in A\\},$

with the usual convention that $\inf\\{\emptyset\\}=\infty$. We will write
$T_{a}:=T_{\\{a\\}}$ for $a\in E$. For $A\subseteq\mathbb{R}$ and $B\subseteq
E_{\partial}$ we write $B^{A}:=B\backslash A$. We will also use the notation
$[b$ for the set $[b,\infty)$. For $A\subseteq\mathbb{R}$ and writing
$\mathbb{P}^{A}=\left(\mathbb{P}^{A}_{x}\right)_{x\in E_{\partial}^{A}}$, the
second auxiliary Markov chain $(X,\mathbb{P}^{A})$ is the one with state space
$E^{A}$ and generated by the Q-matrix $Q^{A}$ given by

$Q^{A}(x,y)=Q(x,y),\quad x,y\in E^{A}.$

Let $P^{A}_{t}(x,y):=\mathbb{P}^{A}_{x}(X_{t}=y)$ be the transition kernel of
$(X,\mathbb{P}^{A})$. By mimicking Steps 1 and 2 of the proof of Theorems
2.8.3 and 2.8.4 in [13] one can show that $\mathbb{P}_{x}(X_{t}=y,t<T^{A})$,
$x,y\in E^{A}$, $t\geq 0$, is the minimal nonnegative solution of the backward
equations associated with $Q^{A}$ and thus
$P^{A}_{t}(x,y)=\mathbb{P}_{x}(X_{t}=y,t<T_{A})$. Therefore we call
$(X,\mathbb{P}^{A})$ the Markov chain obtained by killing $(X,\mathbb{P})$ at
$T_{A}$. Using the semi-group property one can show that, for all $t\geq 0$,
$x\in E$ and $\Gamma\in\mathcal{F}_{t}$,
$\mathbb{P}_{x}^{A}(\Gamma)=\mathbb{P}_{x}(\Gamma\cap\\{t<T_{A}\\})$. Then,
for any stopping time $T$ one can show by first considering stopping times
with a discrete range and then approximating a general stopping time by a non-
increasing sequence of stopping times with a discrete range,

(2.2)
$\mathbb{P}_{x}^{A}(\Gamma\cap\\{T<\infty\\})=\mathbb{P}_{x}(\Gamma\cap\\{T<T_{A}\\}),\quad\Gamma\in\mathcal{F}_{T},\
x\in E.$

The resolvent of $(X,\mathbb{P}^{A})$ is denoted by
$G^{A}(x,y)=\int_{0}^{\infty}P^{A}_{t}(x,y)dt$, $x,y\in E^{A}$ and we extend
it to the whole of $E\times E$ by setting $G^{A}(x,y)=0$ for $(x,y)\notin
E^{A}\times E^{A}$. The corresponding resolvent density is denoted by
$\mathtt{g}^{A}(x,y):=G^{A}(x,y)/\pi(y)$, $x,y\in E$. We can consider the dual
of $(X,\mathbb{P}^{A})$ or kill the dual at $T_{A}$. Clearly by looking at the
Q-matrices involved, the operations of taking the dual (with respect to $\pi$)
and killing at $T_{A}$ commute and the following switching identity holds:

(2.3) $\widehat{\mathtt{g}}^{A}(y,x)=\mathtt{g}^{A}(x,y),\quad x,y\in E,$

where $\widehat{G}^{A}(y,x)$ is the resolvent and
$\widehat{\mathtt{g}}^{A}(y,x)$ its density of the dual killed at $T_{A}$.
Very closely related to $(X,\mathbb{P}^{A})$, is the third auxiliary Markov
chain that we introduce. For $A\subseteq\mathbb{R}$ and writing
$\overline{\mathbb{P}}^{A}=\left(\overline{\mathbb{P}}^{A}_{x}\right)_{x\in
E_{\partial}}$, let $(X,\overline{\mathbb{P}}^{A})$ be the Markov chain with
state space $E$ and generated by the Q-matrix $\overline{Q}^{A}$ given by

$\overline{Q}^{A}(x,y)=\begin{cases}Q(x,y)&\text{if $x\in E^{A}$ and $y\in
E$},\\\ 0&\text{if $x\in A$ and $y\in E$}.\end{cases}$

Similarly as for the previous auxiliary Markov chain, one can show that the
transition kernel
$\overline{P}^{A}_{t}(x,y):=\overline{\mathbb{P}}^{A}_{x}(X_{t}=y)$ of
$(X,\overline{\mathbb{P}}^{A})$ admits the expression
$\overline{P}^{A}_{t}(x,y)=\mathbb{P}_{x}(X_{t\wedge T_{A}}=y)$, $x,y\in E$,
$t\geq 0$. Therefore we call $(X,\overline{\mathbb{P}}^{A})$ the Markov chain
obtained by stopping $(X,\mathbb{P})$ at $T_{A}$. Finally, for $a\in E$, with
slight abuse of notation, we write
$E^{a}_{\partial}:=E_{\partial}\cap[a,\infty)$ and
$\overline{\mathbb{P}}^{a}=\left(\overline{\mathbb{P}}^{a}_{x}\right)_{x\in
E^{a}_{\partial}}$ and denote by $(X,\overline{\mathbb{P}}^{a})$ the Markov
chain obtained by stopping $(X,\mathbb{P})$ at $T_{a}$. Note that the state
space of $(X,\overline{\mathbb{P}}^{a})$ is chosen to be
$E^{a}:=E\cap[a,\infty)$ and not $E$, which we can do because $(X,\mathbb{P})$
is skip-free downward. Its transition kernel is denoted by
$\overline{P}^{a}_{t}$.

##  3\. Potential theory and fluctuation identities of continuous-time skip-
free Markov Process

The following two theorems are, to some extent, the analogue of the results in
Section 3 and the first part of Section 4 in Choi and Patie [6] who deal with
discrete time Markov chains that are skip-free upwards. Though there are some
differences with the most notable ones being the following. First, translating
their approach to our setting of continuous-time skip-free downward Markov
chains, Choi and Patie [6] distinguish between the cases (i)
$\mathfrak{l}=-\infty$ or $\mathfrak{l}>-\infty$ with
$\mathbb{P}_{\mathfrak{l}}(X_{J_{1}}\in E\backslash\\{\mathfrak{l}\\})=0$ on
the one hand and (ii) $\mathfrak{l}>-\infty$ with
$\mathbb{P}_{\mathfrak{l}}(X_{J_{1}}\in E\backslash\\{\mathfrak{l}\\})>0$ on
the other hand, with their results being less clean in the latter case. By
using a different definition for the fundamental excessive function for the
killed process than the one in Theorem 3.1(2) in [6] (though both definitions
are equivalent in case (i)), we are able to unify the two cases and at the
same time give a more elementary proof that avoids the Martin boundary theory
employed in [6]. Second, we allow for state-dependent killing whereas Choi and
Patie [6] worked with a constant killing rate. Recall that notation and
terminology used in the following two theorems can be found in Section 2.

######  Theorem 3.1.

Let $(X,\mathbb{P})$ be a Markov chain with state space $E$ generated by the
Q-matrix $Q$ satisfying (2.1). Assume further all states in $E$ are transient
and fix a reference point $\mathfrak{o}\in E$.

  1. (1)

Define, for $x\in E$,
$H(x)=\lim_{y\downarrow\mathfrak{l}}\frac{G(x,y)}{G(\mathfrak{o},y)}:=\begin{cases}\lim_{y\downarrow-\infty}\frac{G(x,y)}{G(\mathfrak{o},y)}&\text{if
$\mathfrak{l}=-\infty$},\\\
\frac{G(x,\mathfrak{l})}{G(\mathfrak{o},\mathfrak{l})}&\text{if
$\mathfrak{l}>-\infty$}.\end{cases}$ Then for $x,a\in E$ such that $x\geq a$,

(3.1) $\mathbb{P}_{x}(T_{a}<\zeta)=\frac{H(x)}{H(a)}.$

  2. (2)

The function $H$ is excessive for $(X,\mathbb{P})$ and it is the unique
positive-valued function on $E$ satisfying the following four properties: (i)
nonincreasing, (ii) $H(\mathfrak{o})=1$, (iii) for any $a\in E$, $H$
restricted to $E^{a}$ is harmonic for $(X,\overline{\mathbb{P}}^{a})$, (iv) if
$\mathfrak{h}=\infty$ and if there exists $x>a$ such that
$\mathbb{P}_{x}(T_{a}=\zeta=\infty)$, then $\lim_{x\to\infty}H(x)=0$.

  3. (3)

Define
$\widehat{H}(y)=\lim_{x\uparrow\mathfrak{h}}\frac{\widehat{G}(y,x)}{\widehat{G}(\mathfrak{o},x)}$,
$y\in E$, and
$\mathbf{H}^{[b}(x)=\lim_{y\downarrow\mathfrak{l}}\frac{G^{[b}(x,y)}{G(\mathfrak{o},y)}$,
$x\in E$. Then for $x,y\in E$,

(3.2)
$\mathtt{g}(x,y)=\mathfrak{C}\widehat{H}(y)\left(H(x)-\mathbf{H}^{[y}(x)\right)\quad\text{with
\ $\mathfrak{C}=\mathtt{g}(\mathfrak{o},\mathfrak{o})$}.$

###### Proof.

In the proof $x,y,z,a,b$ will always denote elements in $E$.

(1). We first prove the following claim:

(3.3)
$\frac{G(x,y)}{G(\mathfrak{o},y)}=\frac{G(x,x\wedge\mathfrak{o})}{G(o,x\wedge\mathfrak{o})}\quad
y\leq x\wedge\mathfrak{o}.$

By the strong Markov property, $G(z,a)=\mathbb{P}_{z}(T_{a}<\zeta)G(a,a)$ for
any $z,a\in E$. Hence,

$\begin{split}\frac{G(b,y)}{G(y,y)}=\mathbb{P}_{b}(T_{y}<\zeta)=&\mathbb{P}_{b}(T_{z}<\zeta)\mathbb{P}_{z}(T_{y}<\zeta)=\frac{G(b,z)}{G(z,z)}\frac{G(z,y)}{G(y,y)},\quad
y\leq z\leq b,\end{split}$

where the middle equality follows by the strong Markov property and the skip-
free downward property of $(X,\mathbb{P})$. When $x\leq\mathfrak{o}$, take
$z=x$ and $b=\mathfrak{o}$, to get for $y\leq x\leq\mathfrak{o}$,
$G(\mathfrak{o},y)=\frac{G(\mathfrak{o},x)}{G(x,x)}G(x,y)$ which shows (3.3)
in the case where $x\leq\mathfrak{o}$. Similarly, if $\mathfrak{o}\leq x$,
take $z=\mathfrak{o}$ and $b=x$ instead to get for $y\leq\mathfrak{o}\leq x$,
$G(x,y)=\frac{G(x,\mathfrak{o})}{G(\mathfrak{o},\mathfrak{o})}G(\mathfrak{o},y)$,
which shows (3.3) in the case where $x\geq\mathfrak{o}$. From (3.3) it follows
that $H(x)=\lim_{y\downarrow\mathfrak{l}}\frac{G(x,y)}{G(\mathfrak{o},y)}$ is
well-defined. By the strong Markov property and the skip-free downward
property, $G(x,y)=\mathbb{P}_{x}(T_{a}<\zeta)G(a,y)$ for $y\leq a\leq x$.
Hence

$\mathbb{P}_{x}(T_{a}<\zeta)=\frac{G(x,y)}{G(a,y)}=\frac{G(x,y)/G(\mathfrak{o},y)}{G(a,y)/G(\mathfrak{o},y)},\quad
y\leq a\leq x.$

Now taking limits as $y\downarrow\mathfrak{l}$ yields (3.1).

(2). By the semi-group property of $P_{t}$, for $y\leq x$,

$\sum_{z\in E}P_{t}(x,z)G(z,y)=\sum_{z\in
E}P_{t}(x,z)\int_{0}^{\infty}P_{s}(z,y)ds=\int_{0}^{\infty}P_{t+s}(x,y)ds=\int_{t}^{\infty}P_{u}(x,y)du<G(x,y).$

So if $\mathfrak{l}>-\infty$, $P_{t}H(x)<H(x)$ for all $x\in E$ and so $H$ is
excessive (but not harmonic) for $(X,\mathbb{P})$. On the other hand if
$\mathfrak{l}=-\infty$, then $H$ is the limit as $n\to\infty$ of the excessive
functions $H_{n}(x):=\frac{G(x,-n)}{G(\mathfrak{o},-n)}$, $-n\in E$, and thus
is excessive since by Fatou’s lemma,

$P_{t}H=P_{t}\lim H_{n}\leq\liminf P_{t}H_{n}\leq\liminf H_{n}=H.$

The function $H$ is nonincreasing by (3.1) and $(X,\mathbb{P})$ being skip-
free downward. Obviously, $H(\mathfrak{o})=1$ by definition. Fix $a\in E$ and
let $h_{a}(x):=\mathbb{P}_{x}(T_{a}<\zeta)$, $x\in E^{a}$. Then for any $x\in
E^{a}$ and $t\geq 0$,

$\begin{split}\overline{P}_{t}^{a}h_{a}(x)=&\mathbb{E}_{x}\left[h_{a}(X_{t\wedge
T_{a}})\mathds{1}_{\\{t\wedge T_{a}<\zeta\\}}\right]\\\
=&\mathbb{E}_{x}\left[\mathds{1}_{\\{t\geq
T_{a},T_{a}<\zeta\\}}\mathbb{P}_{X_{T_{a}}}(T_{a}<\zeta)\right]+\mathbb{E}_{x}\left[\mathds{1}_{\\{t<T_{a}\wedge\zeta\\}}\mathbb{P}_{x}(T_{a}<\zeta|\mathcal{F}_{t})\right]\\\
=&\mathbb{P}_{x}(t\geq
T_{a},T_{a}<\zeta)+\mathbb{P}_{x}(t<T_{a}\wedge\zeta,T_{a}<\zeta)\\\
=&h_{a}(x),\end{split}$

where we have used the Markov property in the second equality. So $h_{a}$ is
harmonic for $(X,\overline{\mathbb{P}}^{a})$. Then by (3.1) $H$ restricted to
$E^{a}$ is harmonic for $(X,\overline{\mathbb{P}}^{a})$. Let
$h:E\to(0,\infty)$ satisfy (i) and (iii). Then for $x\in E^{a}$,

(3.4)
$h(x){=}\lim_{t\to\infty}\overline{P}_{t}^{a}h(x){=}\mathbb{E}_{x}\left[\lim_{t\to\infty}h(X_{t\wedge
T_{a}})\mathds{1}_{\\{t\wedge
T_{a}<\zeta\\}}\right]=h(a)\mathbb{P}_{x}(T_{a}<\zeta)+\lim_{y\to\infty}h(y)\mathbb{P}_{x}(T_{a}=\zeta=\infty),$

where the second equality is due to the dominated convergence theorem which is
applicable because $h$ is nonincreasing and positive, $x\geq a$ and
$(X,\mathbb{P})$ is skip-free downward and for the last equality we used that
by transience $X_{t}\to\infty$ as $t\to\infty$ if $T_{a}=\zeta=\infty$ (which
can only happen if $\mathfrak{h}=\infty$). So if $h$ satisfies (i)-(iv), then
we conclude $h(x)=h(a)\mathbb{P}_{x}(T_{a}<\zeta)$ for all $x\geq a$, which
yields that $h$ is proportional to $H$ on $E^{a}$ by (3.1). Since $a\in E$ was
chosen arbitrarily, $h$ is positive-valued and
$h(\mathfrak{o})=H(\mathfrak{o})$, it follows that $h=H$ on $E$ if $h$
satisfies (i)-(iv). Since we have showed that $H$ satisfy (i) and (iii), (3.4)
holds for $H$ and thus $H$ satisfies property (iv) by (3.1).

(3). By the strong Markov property and the skip-free upward property of the
dual, for $x,y,b\in E$ such that $y\leq b$,

$\widehat{G}(y,x)=\widehat{G}^{[b}(y,x)+\widehat{\mathbb{P}}_{y}(T_{[b}<\zeta)\widehat{G}(b,x)=\widehat{G}^{[b}(y,x)+\frac{\widehat{H}(y)}{\widehat{H}(b)}\widehat{G}(b,x),$

where for the second equality we used the analogue of (3.1) for the dual which
is justified because $(X,\widehat{\mathbb{P}})$ is skip-free upward and all
its states are transient. Then by the switching identity (2.3),

(3.5)
$G^{[b}(x,y)=G(x,y)-G(x,b)\frac{\widehat{H}(y)\pi(y)}{\widehat{H}(b)\pi(b)},\quad
y\leq b.$

Since by the same arguments
$\frac{G(\mathfrak{o},y)}{\pi(y)}=\mathtt{g}(\mathfrak{o},y)=\widehat{\mathtt{g}}(y,\mathfrak{o})=\widehat{\mathbb{P}}_{y}(T_{\mathfrak{o}}\leq\zeta)\mathtt{g}(\mathfrak{o},\mathfrak{o})=\mathfrak{C}\frac{\widehat{H}(y)}{\widehat{H}(\mathfrak{o})}$
for $y\leq\mathfrak{o}$,

$\frac{G^{[b}(x,y)}{G(\mathfrak{o},y)}=\frac{G(x,y)}{G(\mathfrak{o},y)}-\frac{G(x,b)}{\mathfrak{C}\widehat{H}(b)\pi(b)},\quad
y\leq b\wedge\mathfrak{o}.$

Taking limits as $y\downarrow-\mathfrak{l}$ shows that $\mathbf{H}^{[b}$ is
well-defined and
$\mathbf{H}^{[b}(x)=H(x)-\frac{G(x,b)}{\mathfrak{C}\widehat{H}(b)\pi(b)}$. The
expression for the resolvent density of $G$ then follows. ∎

######  Remark 3.2.

  1. (1)

The function $H$ is not necessarily harmonic for $(X,\mathbb{P})$ if
$\mathfrak{l}=-\infty$, which is in contrast to the discrete time case, see
[6]. For instance when $(X,\mathbb{P})$ is the negative of a birth process
which explodes, then $H\equiv 1$ but $H$ is not harmonic since the lifetime is
finite a.s..

  2. (2)

If $h:E\to(0,\infty)$ is harmonic for $(X,\mathbb{P})$ then $h$ satisfies
property (iii) in Theorem 3.1(2). Indeed, if $h$ is harmonic, then by the
Markov property for any $s\leq t$ and $x\in E$,
$\mathbb{E}_{x}[h(X_{t})\mathds{1}_{\\{t<\zeta\\}}|\mathcal{F}_{s}]=\mathds{1}_{\\{s<\zeta\\}}\mathbb{E}_{X_{s}}[h(X_{t-s})\mathds{1}_{\\{t-s<\zeta\\}}]=h(X_{s})\mathds{1}_{\\{s<\zeta\\}}$.
So $(h(X_{t})\mathds{1}_{\\{t<\zeta\\}})_{t\geq 0}$ is a martingale and
therefore so is $(h(X_{t\wedge T_{a}})\mathds{1}_{\\{{t\wedge
T_{a}}<\zeta\\}})_{t\geq 0}$ for any $a\in E$. This yields, for any $a\in E$,
$x\in E^{a}$ and $t\geq 0$,
$\overline{P}_{t}^{a}h(x)=\overline{P}_{0}^{a}h(x)=h(x)$.

  3. (3)

Property (iv) in Theorem 3.1(2) is needed to characterise $H$. Indeed, if
$\mathfrak{h}=\infty$ and $(X,\mathbb{P})$ is such that
$\lim_{t\to\infty}X_{t}=\infty$ $\mathbb{P}_{x}$-a.s., then $h\equiv 1$
satisfies properties (i)-(iii) but $H\neq h$.

######  Theorem 3.3.

Under the assumptions of Theorem 3.1 and using its notation we have the
following.

  1. (1)

For any $x,y\in E$,

(3.6) $\mathbb{P}_{x}(T_{y}<\zeta)=\frac{H(x)-\mathbf{H}^{[y}(x)}{H(y)}.$

  2. (2)

For any $x,a,b\in E$ such that $x\geq a$ and $b>a$,

(3.7)
$\mathbb{P}_{x}(T_{a}<T_{[b}\wedge\zeta)=\frac{\mathbf{H}^{[b}(x)}{\mathbf{H}^{[b}(a)}.$

  3. (3)

Let $a,b\in E$ such that $a<b$. Assume $f:E\to\mathbb{R}$ is such that
$Q|f|(x)<\infty$ for all $x\in E^{(a,b)^{\complement}}$ where
$|f|(x):=|f(x)|$. Then for $x\in E$,

(3.8)
$\mathbb{E}_{x}\left[f\left(X_{T_{(a,b)^{\complement}}}\right)\mathds{1}_{\\{T_{(a,b)^{\complement}}<\zeta\\}}\right]=f(x)+\sum_{z\in
E^{(a,b)^{\complement}}}Qf(z)G^{{(a,b)^{\complement}}}(x,z),$

where the density of $G^{{(a,b)^{\complement}}}$ (with respect to $\pi$)
admits the expression,

(3.9)
$\mathtt{g}^{{(a,b)^{\complement}}}(x,y)=\mathfrak{C}\widehat{H}(y)\left(\frac{\mathbf{H}^{[b}(x)}{\mathbf{H}^{[b}(a)}\mathbf{H}^{[y}(a)-\mathbf{H}^{[y}(x)\right),\quad
x,y\in E^{(a,b)^{\complement}}.$

In particular for $x\in E$,

(3.10) $\mathbb{P}_{x}(T_{(a,b)^{\complement}}<\zeta)=1-\sum_{z\in
E^{(a,b)^{\complement}}}Q(z,\partial)G^{{(a,b)^{\complement}}}(x,z).$

  4. (4)

If $\mathfrak{l}=-\infty$, denote
$\mathbf{H}^{[b}(-\infty):=\lim_{x\to-\infty}\mathbf{H}^{[b}(x)\in(0,\infty]$.
Then for $x,b\in E$,

(3.11) $\mathbb{P}_{x}(T_{[b}<\zeta)=\begin{cases}1-\sum_{z\in
E^{[b}}Q(z,\partial)G^{[b}(x,z)-\frac{\mathbf{H}^{[b}(x)}{\mathbf{H}^{[b}(-\infty)}&\text{if
$\mathfrak{l}=-\infty$},\\\ 1-\sum_{z\in
E^{[b}}Q(z,\partial)G^{[b}(x,z)&\text{if $\mathfrak{l}>-\infty$},\end{cases}$

where the density of $G^{[b}$ admits the expression,

(3.12)
$\mathtt{g}^{[b}(x,y)=\mathfrak{C}\widehat{H}(y)\left(\mathbf{H}^{[b}(x)-\mathbf{H}^{[y}(x)\right),\quad
x,y\in E^{[b}.$

######  Remark 3.4.

Although somewhat hidden, the identities in Theorem 3.3 provide in particular
expressions for the Laplace transform of the various first hitting/passage
times by including a constant term in the killing rates $Q(x,\partial)$, $x\in
E$. Namely, if, for $q>0$, $(X,\mathbb{P}^{(q)})$ is the Markov chain with
state space $E$ generated by the Q-matrix $Q^{(q)}$ given by
$Q^{(q)}(x,y)=Q(x,y)-q\delta_{xy}$, then its transition kernel is given by
$P_{t}^{(q)}(x,y)=e^{-qt}P_{t}(x,y)$ and so, following the arguments leading
up to (2.2), one has for any stopping time $T$,

$\mathbb{P}^{(q)}_{x}(\Gamma\cap\\{T<\infty\\})=\mathbb{E}_{x}\left[e^{-qT}\mathds{1}_{\Gamma\cap\\{T<\infty\\}}\right],\quad\Gamma\in\mathcal{F}_{T},\
x\in E.$

So using the identities in Theorem 3.3 for $(X,\mathbb{P}^{(q)})$ yields
identities involving the Laplace transform of the corresponding first
hitting/passage time with Laplace parameter $q$ for $(X,\mathbb{P})$, see also
Remark 5.5 below.

###### Proof.

(1). By the strong Markov property $G(x,y)=\mathbb{P}_{x}(T_{y}<\zeta)G(y,y)$
for $x,y\in E$. The identity (3.6) then follows from (3.2) since
$\mathbf{H}^{[y}(x)=0$ for $x\geq y$.

(2). Since $(X,\mathbb{P}^{[b})$ is skip-free downward with all states
transient (3.1) is applicable for this Markov chain. Hence by (2.2) and with
$H^{[b}(x):=\lim_{y\downarrow\mathfrak{l}}\frac{G^{[b}(x,y)}{G^{[b}(\mathfrak{o}^{[b},y)}$
where $\mathfrak{o}^{[b}\in E^{[b}$ is a reference point,

$\mathbb{P}_{x}(T_{a}<T_{[b}\wedge\zeta)=\mathbb{P}^{[b}(T_{a}<\zeta)=\frac{H^{[b}(x)}{H^{[b}(a)},\quad
x\geq a,$

The required identity then follows since
$\mathbf{H}^{[b}(x)=\lim_{y\downarrow\mathfrak{l}}\frac{G^{[b}(x,y)}{G^{[b}(\mathfrak{o}^{[b},y)}\frac{G^{[b}(\mathfrak{o}^{[b},y)}{G(\mathfrak{o},y)}=H^{[b}(x)\mathbf{H}^{[b}(o^{[b})$,
$x\in E$.

(3). We first prove (3.9) and (3.12). The latter follows directly from (3.5)
and (3.2). Since by the strong Markov property and $(X,\mathbb{P})$ being
skip-free downward,

$\mathtt{g}^{[b}(x,y)=\mathtt{g}^{{(a,b)^{\complement}}}(x,y)+\mathbb{P}_{x}^{[b}(T_{a}<\zeta)\mathtt{g}^{[b}(a,y),\quad
x\geq a,y\in E,$

(3.9) follows from (3.7) and (3.12). Now we prove (3.8). Note that without
loss of generality we can assume that $f$ is bounded because the general case
can be proved by approximation from the case where $f$ is bounded. By the
forward equations,
$\overline{P}^{(a,b)^{\complement}}_{t}(x,y)=\delta_{xy}+\int_{0}^{t}\sum_{z\in
E}\overline{P}^{(a,b)^{\complement}}_{s}(x,z)\overline{Q}^{(a,b)^{\complement}}(z,y)ds$,
$x,y\in E$, $t\geq 0$. Therefore, for $f:E\to\mathbb{R}$ bounded, $x\in E$ and
$t\geq 0$, by definition of $\overline{Q}^{(a,b)^{\complement}}$,

$\begin{split}\left(\overline{P}^{(a,b)^{\complement}}_{t}f\right)(x)=&f(x)+\sum_{z\in
E^{(a,b)^{\complement}}}Qf(z)\int_{0}^{t}\overline{P}^{(a,b)^{\complement}}_{s}(x,z)ds,\end{split}$

where the interchanging of the two sums is justified because
$E^{(a,b)^{\complement}}$ is a finite set. By the transience assumption,
$T_{(a,b)^{\complement}}\wedge\zeta<\infty$ $\mathbb{P}_{x}$-a.s. for all
$x\in E^{(a,b)^{\complement}}$. As $f$ is assumed to be bounded, the dominated
convergence theorem yields,

$\begin{split}\mathbb{E}_{x}\left[f(X_{T_{(a,b)^{\complement}}})\mathbf{1}_{\\{T_{(a,b)^{\complement}}<\zeta\\}}\right]=&\lim_{t\to\infty}\mathbb{E}_{x}\left[f(X_{t\wedge
T_{(a,b)^{\complement}}})\mathbf{1}_{\\{t\wedge
T_{(a,b)^{\complement}}<\zeta\\}}\right]\\\
=&\lim_{t\to\infty}\left(\overline{P}^{(a,b)^{\complement}}_{t}f\right)(x)\\\
=&f(x)+\sum_{z\in
E^{(a,b)^{\complement}}}Qf(z)\int_{0}^{\infty}\overline{P}^{(a,b)^{\complement}}_{s}(x,z)ds\\\
=&f(x)+\sum_{z\in
E^{(a,b)^{\complement}}}Qf(z)G^{(a,b)^{\complement}}(x,z).\end{split}$

The identity (3.10) then simply follows by taking $f\equiv 1$ and the
definition of $Q(x,\partial)$.

(4). When $\mathfrak{l}>-\infty$ we can use the same arguments for (3.10) to
prove (3.11) since $E^{[b}$ is a finite set and $T_{[b}\wedge\zeta<\infty$
$\mathbb{P}_{x}$-a.s. for all $x\in E^{[b}$ when $\mathfrak{l}>-\infty$. Now
assume $\mathfrak{l}=-\infty$. By (3.10) and (3.7), for $x,a,b\in E$ such that
$a<b,x$,

$\mathbb{P}_{x}(T_{[b}<T_{a}\wedge\zeta)=\mathbb{P}_{x}(T_{(a,b)^{\complement}}<\zeta)-\mathbb{P}_{x}(T_{a}<T_{[b}\wedge\zeta)=1-\sum_{z\in
E^{(a,b)^{\complement}}}Q(z,\partial)G^{{(a,b)^{\complement}}}(x,z)-\frac{\mathbf{H}^{[b}(x)}{\mathbf{H}^{[b}(a)}.$

Since $\mathfrak{l}=-\infty$ we have for any $x\in E$, $\mathbb{P}_{x}$-a.s.,
$T_{a}\wedge\zeta\to\zeta$ as $a\to-\infty$ since all states are assumed to be
transient and the Markov chain is sent to the cemetery state when explosion
occurs. So (3.11) for $\mathfrak{l}=-\infty$ follows by taking $a\to-\infty$
in the last equation and using the monotone convergence theorem. ∎

##  4\. Downward skip-free compound Poisson processes

In this section, we apply the methodology described in Section 3 to recover
various fluctuation identities for downward skip-free $\mathbb{Z}$-valued
compound Poisson processes. We recall that for the continuous-state space
analogue of this class of Markov processes, namely spectrally positive Lévy
processes, this problem has been well studied and has found an impressive
range of applications, such as insurance mathematics, epidemiology, financial
mathematics and queuing theory, see [11, 7]. For the discrete-time analogue,
i.e. for random walks, Spitzer [15] solved this problem by means of the
celebrated Wiener-Hopf factorization and alternative interesting proofs, based
either on excursion theory or martingales devices have been proposed, see
again [11, 7] and the references therein. However, all of these approaches
rely on the stationarity and independent increments property of the random
walks/Lévy processes, which makes it difficult to extend them to a wider
context. Note that the continuous-time, discrete-state space case considered
here has been studied in [16].

In this section we assume that $E=\mathbb{Z}$ and the Q-matrix generating
$(X,\mathbb{P})$ is given by, for $x,y\in\mathbb{Z}$,

(4.1) $Q(x,y)=\begin{cases}\mu(y-x+1)\alpha&\text{if $y\geq x+1$ or
$y=x-1$},\\\ -(\alpha+p)&\text{if $y=x$},\\\ 0&\text{if $y\leq
x-2$},\end{cases}$

where $\alpha>0$, $\mu=(\mu(j))_{j\geq 0}$ is a probability measure with
$\mu(0)>0$ and $\mu(1)=0$ and $p\geq 0$. Note that $Q(x,\partial)=p$ for
$x\in\mathbb{Z}$. This means that the holding times are i.i.d. and
exponentially distributed and the jump sizes are i.i.d. and independent from
the holding times. So $(X,\mathbb{P})$ corresponds to a $\mathbb{Z}$-valued
skip-free downward compound Poisson process with intensity $\alpha$, jump
distribution $\mathbb{P}_{0}(X_{J_{1}}=j)=\mu(j+1)$, $j\geq-1$, and killed at
the constant rate $p$. For $f:E\to\mathbb{R}$ we denote its generating
function (GF) by $f^{*}[s]:=\sum_{y\in E}s^{y}f(y)$. We have for
$x\in\mathbb{Z}$ and $s\in(0,1]$,

(4.2) $Q(x,\cdot)^{*}[s]=\psi(s)s^{x-1},$

where, setting $\mu(j)=0$ for $j\leq-1$,

$\psi(s)=\alpha\left(\mu^{*}[s]-s\right)-ps.$

Since $\psi$ is convex, $\psi(0)=\alpha\mu_{0}>0$ and $\psi(1)=-p\leq 0$, we
have $s_{0}\in(0,1]$ and $\psi(s)>0$ for $s\in[0,s_{0})$. Note that
$P_{t}(x,\cdot)^{*}[s]$ is well-defined for $s=e^{iu}$ for $u\in\mathbb{R}$.
Using the forward equations, Fubini and (4.2), we get for $t\geq 0$,
$x\in\mathbb{Z}$ and $s=e^{iu}$ with $u\in\mathbb{R}$,

$P_{t}(x,\cdot)^{*}[s]=1+\int_{0}^{t}\sum_{z\in\mathbb{Z}}P_{v}(x,z)Q(z,\cdot)^{*}[s]dv=1+\frac{\psi(s)}{s}\int_{0}^{t}P_{v}(x,\cdot)^{*}[s]dv.$

So by taking derivatives in $t$, we get
$\frac{d}{dt}P_{t}(x,\cdot)^{*}[s]=\frac{\psi(s)}{s}P_{t}(x,\cdot)^{*}[s]$.
Solving this simple ODE with the boundary condition
$P_{0}(x,\cdot)^{*}[s]=s^{x}$ yields, for $t\geq 0$, $x\in\mathbb{Z}$ and
$s=e^{iu}$ with $u\in\mathbb{R}$,

(4.3)
$\mathbb{E}_{x}\left[s^{X_{t}}\mathds{1}_{\\{t<\zeta\\}}\right]=P_{t}(x,\cdot)^{*}[s]=s^{x}\exp\left(\frac{\psi(s)}{s}t\right).$

By analytic extension (4.3) also holds for $s\in(0,1]$. Also (4.3) implies the
spatial homogeneity property $P_{t}(x,y)=P_{t}(z,y-x+z)$ for
$x,y,z\in\mathbb{Z}$. Clearly all states are transient if $p>0$. Further, as
can be deduced from the theory of random walks, if $p=0$ then all states are
transient if and only if $\mathbb{E}_{0}[X_{1}]\neq 0$ and by differentiating
in $s$ both sides of (4.3) we see that
$\mathbb{E}_{0}[X_{1}]=\psi^{\prime}(1)=\alpha(m-1)$ where $m:=\sum_{j\geq
2}j\mu(j)$. Hence all states are transient for $(X,\mathbb{P})$ if and only if
$p>0$ or $m\neq 1$. Otherwise all states are recurrent. In order to present
the fluctuation identities we introduce the function
$W:\mathbb{Z}\to[0,\infty)$ defined by $W(x)=0$ for $x<0$,
$W(0)=\frac{1}{\alpha\mu(0)}$ and for $x\geq 0$,

(4.4)
$W(x+1)=\frac{1}{\alpha\mu(0)}+\frac{1}{\mu(0)}\sum_{j=0}^{x}W(j)\left(\overline{\mu}(x-j+1)+\frac{p}{\alpha}\right),$

where $\overline{\mu}(k)=1-\sum_{j\leq k}\mu(j)$. It is easy to verify that
the function $W$ on $\mathbb{Z}\cap[0,\infty)$ is identical to the one defined
in Lemma 5.1 below and so from that lemma we get that the GF of $W$ equals
$W^{*}[s]=\frac{1}{\psi(s)}$ for $s\in[0,s_{0})$. We now provide a suite of
fluctuation identities for skip-free downward compound Poisson processes
(killed at a constant rate). Recall that some of these identities can be found
in [16] whereas all of them for the closely related class of spectrally
negative Lévy processes can be found in Chapter 8 of [11].

######  Theorem 4.1.

We have the following.

  1. (1)

If $p>0$ or $m\neq 1$, then the resolvent of $(X,\mathbb{P})$ is given by, for
$x,y\in\mathbb{Z}$,

$\begin{split}G(x,y)=-\frac{s_{0}^{x-y}}{\psi^{\prime}(s_{0})}-W(y-x-1).\end{split}$

  2. (2)

In any case, for $x,y\in\mathbb{Z}$,

$\mathbb{P}_{x}(T_{y}<\zeta)=s_{0}^{x-y}+\psi^{\prime}(s_{0})W(y-x-1).$

  3. (3)

In any case, for $a\leq x\leq b-1$,

$\begin{split}\mathbb{P}_{x}(T_{a}<T_{[b}\wedge\zeta)=&\frac{W(b-x-1)}{W(b-a-1)},\\\
\mathbb{P}_{x}(T_{(a,b)^{\complement}}<\zeta)=&1+p\sum_{z=0}^{b-x-2}W(z)-\frac{W(b-x-1)}{W(b-a-1)}p\sum_{z=0}^{b-a-2}W(z),\\\
\mathbb{P}_{x}(T_{[b}<\zeta)=&\begin{cases}1+p\sum_{z=0}^{b-x-2}W(z)-\frac{ps_{0}}{1-s_{0}}W(b-x-1)&\text{if
$p>0$ or $m>1$},\\\ 1+\psi^{\prime}(1)W(b-x-1)&\text{if $p=0$ and $m\leq
1$}.\end{cases}\end{split}$

###### Proof.

We first assume $p>0$. Then all states are transient and so we can use the
results from Section 3. We choose for the reference point $\mathfrak{o}=0$ and
we choose for the reference measure $\pi\equiv 1$. Note that by spatial
homogeneity
$\sum_{y\in\mathbb{Z}}P_{t}(x,y)=\sum_{y\in\mathbb{Z}}P_{t}(0,y-x)\leq 1$ for
any $x\in\mathbb{Z}$ and so $\pi\equiv 1$ is excessive for $(X,\mathbb{P})$.
We have by (4.3), for any $x\in\mathbb{Z}$,

$\sum_{y\in
E}P_{t}(x,y)s_{0}^{y}=s_{0}^{x}\exp\left(\frac{\psi(s_{0})}{s}t\right)=s_{0}^{x}.$

So $x\mapsto s_{0}^{x}$ is harmonic for $(X,\mathbb{P})$ and so, via Remark
(3.2) and since $s_{0}\in(0,1)$ as $p>0$, satisfies properties (i)-(iv) in
Theorem 3.1(2). Hence $H(x)=s_{0}^{x}$. Further, as
$\widehat{G}(y,x)=G(x,y)=G(-y,-x)$ for all $y\in\mathbb{Z}$, we have by the
definitions of $H$ and $\widehat{H}$ in Theorem 3.1, for all $y\in\mathbb{Z}$,

$\widehat{H}(y)=\lim_{x\uparrow\infty}\frac{\widehat{G}(y,x)}{\widehat{G}(0,x)}=\lim_{x\uparrow\infty}\frac{G(-y,-x)}{G(0,-x)}=\lim_{z\downarrow-\infty}\frac{G(-y,z)}{G(0,z)}=H(-y)=s_{0}^{-y}.$

Since $\psi(s)<0$ for $s\in(s_{0},1]$ we have by (4.3), for $s\in(s_{0},1]$,

$G(0,\cdot)^{*}[s]=\int_{0}^{\infty}\exp\left(\frac{\psi(s)}{s}t\right)dt=-\frac{s}{\psi(s)}.$

By (3.2) and recalling $\mathbf{H}^{[y}(0)=0$ for $y\leq 0$, we have for
$s\in(s_{0},1]$,

(4.5) $\begin{split}\sum_{y\geq 1}s^{y}G(0,y)=G(0,\cdot)^{*}[s]-\sum_{y\leq
0}\mathfrak{C}(s/s_{0})^{y}=-\frac{s}{\psi(s)}-\frac{\mathfrak{C}}{1-s_{0}/s}=&-\frac{(s-s_{0})s+\mathfrak{C}\psi(s)s}{\psi(s)(s-s_{0})}.\end{split}$

Since $G(0,y)\leq 1/p$ for all $y\in\mathbb{Z}$, we must have $\sum_{y\geq
1}s^{y}G(0,y)<\infty$ for all $s\in[0,1]$. Hence by continuity of a power
series

$-\lim_{s\downarrow
s_{0}}\frac{(s-s_{0})s+\mathfrak{C}\psi(s)s}{\psi(s)(s-s_{0})}=\lim_{s\downarrow
s_{0}}\sum_{y\geq 1}s^{y}G(0,y)\in[0,\infty).$

But by L’Hôpital, if $\mathfrak{C}\neq-\frac{1}{\psi^{\prime}(s_{0})}$, then

$\left|\lim_{s\downarrow
s_{0}}\frac{(s-s_{0})s+\mathfrak{C}\psi(s)s}{\psi(s)(s-s_{0})}\right|=\left|\lim_{s\downarrow
s_{0}}\frac{s-s_{0}+s(1+\mathfrak{C}\psi^{\prime}(s))+\mathfrak{C}\psi(s)}{\psi^{\prime}(s)(s-s_{0})+\psi(s)}\right|=\infty.$

So we must have $\mathfrak{C}=-\frac{1}{\psi^{\prime}(s_{0})}$. Then by
analytic continuation, (4.5) holds for all $s\in[0,1]$. Therefore, since
$G(0,y)=\mathfrak{C}s_{0}^{-y}(1-\mathbf{H}^{[y}(0))$ for $y\in\mathbb{Z}$ by
(3.2), we have for $s\in[0,s_{0})$,

$\sum_{y\geq
1}(s/s_{0})^{y}\mathbf{H}^{[y}(0)=-\frac{1}{\mathfrak{C}}\sum_{y\geq
1}s^{y}G(0,y)+\frac{s/s_{0}}{1-s/s_{0}}=-\psi^{\prime}(s_{0})\frac{s}{\psi(s)}.$

So $\mathbf{H}^{[y}(0)=-\psi^{\prime}(s_{0})s_{0}^{y}W(y-1)$ (recall
$W^{*}[s]=\frac{1}{\psi(s)}$ for $s\in[0,s_{0})$) and so
$G(0,y)=\mathfrak{C}s_{0}^{-y}(1-\mathbf{H}^{[y}(0))=-\frac{s_{0}^{-y}}{\psi^{\prime}(s_{0})}-W(y-1)$.
By spatial homogeneity $G(x,y)=G(0,y-x)$ which proves item (1) for the case
$p>0$ and further implies by (3.2), for $x,y\in\mathbb{Z}$,

$\mathbf{H}^{[y}(x)=-\psi^{\prime}(s_{0})s_{0}^{y}W(y-x-1).$

If $p=0$, we consider for $q>0$ the process $(X,\mathbb{P}^{(q)})$ which is
defined to be the process $(X,\mathbb{P})$ but killed at rate $q$, i.e. its
Q-matrix is given by (4.1) but with $p=q$. Then clearly from (4.3), the
transition kernel of $(X,\mathbb{P}^{(q)})$ is given by
$P^{(q)}_{t}(x,y)=e^{-qt}P_{t}(x,y)$. By the same reasoning as for (2.2), we
further have for any stopping time $T$,

(4.6)
$\mathbb{P}_{x}^{(q)}(\Gamma\cap\\{T<\infty\\})=\mathbb{E}_{x}\left[e^{-qT}\mathds{1}_{\Gamma\cap\\{T<\infty\\}}\right],\quad\Gamma\in\mathcal{F}_{T},\
x\in E.$

Since all states are transient for $(X,\mathbb{P}^{(q)})$ we have with the
obvious notation,

$G^{(q)}(x,y)=\int_{0}^{\infty}e^{-qt}P_{t}(x,y)dy=-\frac{s_{0}(q)^{x-y}}{\psi^{\prime}(s_{0}(q))-q}-W^{(q)}(y-x-1),\quad
x,y\in\mathbb{Z}.$

By the monotone convergence theorem $G^{(q)}(x,y)\to G(x,y)$ as $q\downarrow
0$ and it is easy to see by convexity of $\psi$ that $s_{0}(q)\to s_{0}$ and
thus also $\psi^{\prime}(s_{0}(q))\to\psi^{\prime}(s_{0})$ as $q\downarrow 0$.
Further, via (4.4) it is easily seen that $W^{(q)}(x)\to W(x)$ as $q\downarrow
0$ for any $x\in\mathbb{Z}$. Since $\psi^{\prime}(s_{0})\neq 0$ when $p=0$ and
$m\neq 1$, we can conclude that item (1) also holds when $m\neq 1$. Note that
the above expressions found for $H(x)$, $\widehat{H}(y)$ (given the same
reference measure $\pi\equiv 1$) and $\mathbf{H}^{[y}(x)$ still hold for the
case $p=0$ and $m\neq 1$ (for which we know all states are transient as well).
Then the identities in items (2) and (3) in the case $p>0$ or $m\neq 1$
immediately follow from Theorem 3.3 in combination with, relevant for the last
identity in item (3),

$\lim_{x\to\infty}W(x)=\begin{cases}\infty&\text{if $p>0$ or $m>1$},\\\
-\frac{1}{\psi^{\prime}(1)}&\text{if $p=0$ and $m<1$}.\end{cases}$

In order to see this, note that $\lim_{x\to\infty}W(x)=\lim_{s\uparrow
1}\Delta W^{*}[s]$ where $\Delta W(x)=W(x)-W(x-1)$, $x\in\mathbb{Z}$, and that
(i) $s_{0}<1$ if $p>0$ or $m>1$ and so $\lim_{s\uparrow 1}\Delta
W^{*}[s]\geq\lim_{s\uparrow s_{0}}\Delta W^{*}[s]=\lim_{s\uparrow
s_{0}}\frac{1-s}{\psi(s)}=\infty$ and (ii) $s_{0}=1$ and $\psi(1)=0$ if $p=0$
and $m<1$ and so $\lim_{s\uparrow 1}\Delta W^{*}[s]=\lim_{s\uparrow
1}\frac{1-s}{\psi(s)}=-\frac{1}{\psi^{\prime}(1)}$. If on the other hand,
$p=0$ and $m=1$, then $s_{0}=1$, $\psi^{\prime}(1)=0$ and all states are
recurrent. Hence item (2) and the last identity in item (3) still hold when
$p=0$ and $m=1$. The other two identities in item (3), in the case where $p=0$
and $m=1$, can be proved by considering the process $(X,\mathbb{P}^{(q)})$ for
$q>0$ and taking limits as $q\downarrow 0$ while using (4.6). ∎

##  5\. Markov Branching Processes with Immigration

In this section, we apply the methodology developed in Section 3 to establish
the fluctuation theory of Markov branching processes with immigration, which
are downward skip-free continuous-time Markov chains. We emphasize that
although this class of processes have been intensively studied, only the law
of the downward first passage time, that is the one with continuous crossing,
has been characterized through its Laplace transform, see [5, 17] and in the
continuous state space setting [8]. Further remarks about the literature can
be found near the end of the section.

We assume that $E=\mathbb{Z}\cap[0,\infty)$ is the set of nonnegative integers
and the Q-matrix generating $(X,\mathbb{P})$ is given by, for $x,y\in E$,

(5.1) $Q(x,y)=\begin{cases}\mu(y-x+1)\alpha x+\nu(y-x)\beta&\text{if $y\geq
x+1$},\\\ -(\alpha+p)x-(\beta+q)&\text{if $y=x$},\\\ \mu(0)\alpha x&\text{if
$y=x-1$},\\\ 0&\text{if $y\leq x-2$},\end{cases}$

where $\alpha>0$, $\mu=(\mu(j))_{j\geq 0}$ is a probability measure with
$\mu(0)>0$ and $\mu(1)=0$, $p\geq 0$, $\beta\geq 0$, $\nu=(\nu(j))_{j\geq 1}$
is a probability measure and $q\geq 0$. Note that $Q(x,\partial)=px+q$ for
$x\in E$. We call $(X,\mathbb{P})$ a Markov branching process with immigration
(MBI) with parameters $(\alpha,\mu,p,\beta,\nu,q)$. If $\beta=q=0$ then we
call the MBI a Markov branching process (MBP) with parameters
$(\alpha,\beta,p)$. MBIs are skip-free downward (note $\alpha\mu(0)>0$)
continuous-time Markov chains with transition rates that are affine in the
state variable and were introduced by Sevast’yanov [14]. For more background
information on MBIs/MBPs we refer to Section 1 of Li et al. [12] and the
references therein or Chapter V of Harris [9] or Chapter III of Athreya and
Ney [4]. For $f:E\to\mathbb{R}$ we denote its generating function (GF) by
$f^{*}[s]:=\sum_{y\in E}s^{y}f(y)$ and for $f,g:E\to\mathbb{R}$ we denote by
$f\star g:E\to\mathbb{R}$ the convolution of $f$ and $g$, i.e. $(f\star
g)(x)=\sum_{y=0}^{x}f(x-y)g(y)$. The branching and immigration mechanisms
$\psi$ and $\phi$ of an MBI are defined by, for $s\in[0,1]$,

$\begin{split}\psi(s)=&\alpha\left(\mu^{*}[s]-s\right)-ps,\\\
\phi(s)=&\beta\left(1-\nu^{*}[s]\right)+q,\end{split}$

where we have set $\nu(0)=0$. Further, we define

$s_{0}=\inf\\{s>0:\psi(s)=0\\}.$

Since $\psi$ is convex, $\psi(0)=\alpha\mu_{0}>0$ and $\psi(1)=-p\leq 0$, we
have $s_{0}\in(0,1]$ and $\psi(s)>0$ for $s\in[0,s_{0})$. From (5.1),

(5.2) $Q(x,\cdot)^{*}[s]=\psi(s)xs^{x-1}-\phi(s)s^{x},\quad x\geq 0,\
s\in[0,1],$

where we understand that $xs^{x-1}=0$ when $x=0$ and $s=0$. Before we derive
the fluctuation identities for the MBI using the theory from Section 3, we
need a few preparatory lemmas.

######  Lemma 5.1.

Let $\overline{\mu}=(\overline{\mu}(k))_{k\geq 0}$ and
$\overline{\nu}=(\overline{\nu}(k))_{k\geq 0}$ be the tail measures of $\mu$
and $\nu$, i.e. for $k\geq 0$,

$\overline{\mu}(k)=1-\sum_{0\leq j\leq
k}\mu(j),\quad\overline{\nu}(k)=1-\sum_{0\leq j\leq k}\nu(j),$

where we have set $\nu(0):=0$. Define the function $\Delta W:E\to[0,\infty)$
recursively by $\Delta W(0)=\frac{1}{\alpha\mu(0)}$ and for $k\geq 0$,

$\Delta W(k+1)=\frac{1}{\mu(0)}\sum_{j=0}^{k}\Delta
W(j)\left(\overline{\mu}(k-j+1)+\frac{p}{\alpha}\right).$

Define the functions $W,\kappa:E\to[0,\infty)$ by

$W(k)=\sum_{j=0}^{k}\Delta W(j)\quad\text{and}\quad\kappa(k)=\beta(\Delta
W\star\overline{\nu})(k)+qW(k),\quad k\geq 0.$

Define the nonnegative measure $\pi=(\pi(k))_{k\geq 0}$ and the signed measure
$\varpi=(\varpi(k))_{k\geq 0}$ recursively by $\pi(0)=\varpi(0)=1$ and for
$k\geq 0$,

$\pi(k+1)=\frac{1}{k+1}(\pi\star\kappa)(k)\quad\text{and}\quad\varpi(k+1)=-\frac{1}{k+1}(\varpi\star\kappa)(k).$

Then the GFs of $\Delta W$, $W$, $\kappa$, $\pi$ and $\varpi$ are given by,
for $s\in[0,s_{0})$,

$\displaystyle\Delta W^{*}[s]=$ $\displaystyle\frac{1-s}{\psi(s)}$
$\displaystyle W^{*}[s]=$ $\displaystyle\frac{1}{\psi(s)}$
$\displaystyle\kappa^{*}[s]=$ $\displaystyle\frac{\phi(s)}{\psi(s)}$
$\displaystyle\pi^{*}[s]=$
$\displaystyle\exp\left(\int_{0}^{s}\frac{\phi(u)}{\psi(u)}du\right)$
$\displaystyle\varpi^{*}[s]=$
$\displaystyle\exp\left(-\int_{0}^{s}\frac{\phi(u)}{\psi(u)}du\right).$

###### Proof.

Let $f:E\to[0,\infty)$ with $f(0)=0$ be given and let $f^{\star n}$ denote the
$n$-fold convolution of $f$ with $f^{\star 0}(x)=\delta_{0x}$. It is well-
known that, given $u_{0}\geq 0$, $u(x)=u_{0}\sum_{n\geq 0}f^{\star n}(x)$, is
well-defined and satisfies the renewal sequence

(5.3) $u(x+1)=\sum_{j=0}^{x}u(j)f(x+1-j),\quad x\geq 0,\quad u(0)=u_{0},$

see e.g. Theorem V.2.4 in [3] for the arguments. So if $f^{*}[s]<1$ for some
$s>0$, then $u^{*}[s]=\frac{u_{0}}{1-f^{*}[s]}<\infty$. Let
$\eta(x)=\frac{1}{\mu(0)}(\overline{\mu}(x)+\frac{p}{\alpha})$ for $x\geq 1$
and set $\eta(0)=0$. Since, for $s\in[0,1)$,
$\overline{\mu}^{*}[s]=\frac{1-\mu^{*}[s]}{1-s}$ we have

$\eta^{*}[s]=\frac{1}{\mu(0)}\left(\frac{1-\mu^{*}[s]}{1-s}-(1-\mu(0))+\frac{p}{\alpha(1-s)}\right)=1-\frac{\psi(s)}{\alpha\mu(0)(1-s)}.$

So $\eta^{*}[s]<1$ for $s\in[0,s_{0})$ and so we conclude that for
$s\in[0,s_{0})$, $\Delta W^{*}[s]=\frac{\Delta
W(0)}{1-\eta^{*}[s]}=\frac{1-s}{\psi(s)}$. Consequently,
$W^{*}[s]=\frac{\Delta W^{*}[s]}{1-s}=\frac{1}{\psi(s)}$ for $s\in[0,s_{0})$.
Since $\overline{\nu}^{*}[s]=\frac{1-\nu^{*}[s]}{1-s}$, we have
$\kappa^{*}[s]=\beta\left(\frac{1-s}{\psi(s)}\frac{1-\nu^{*}[s]}{1-s}\right)+\frac{q}{\psi(s)}=\frac{\phi(s)}{\psi(s)}$
for $s\in[0,s_{0})$. From the definitions of $\pi$ and $\varpi$ it follows
that $|\varpi(x)|\leq\pi(x)\leq u(x)$ where $u$ is given by the renewal
sequence (5.3) with $f(j)=\kappa(j-1)$, $j\geq 1$. Since
$f^{*}[s]=\kappa^{*}[s]s$ there exists $\epsilon>0$ such that $f^{*}[s]<1$ and
consequently $u^{*}[s]<\infty$ for $s\in[0,\epsilon]$. Hence $\pi^{*}[s]$ and
$\varpi^{*}[s]$ are well-defined for $s\in[0,\epsilon]$. By the recursive
definitions of $\pi$ and $\varpi$,
$\frac{d}{ds}\pi^{*}[s]=\pi^{*}[s]\frac{\phi(s)}{\psi(s)}$ with $\pi^{*}[0]=1$
and $\frac{d}{ds}\varpi^{*}[s]=-\varpi^{*}[s]\frac{\phi(s)}{\psi(s)}$ with
$\varpi^{*}[0]=1$. By solving these ODEs we get
$\pi^{*}[s]=\exp\left(\int_{0}^{s}\frac{\phi(u)}{\psi(u)}du\right)=\frac{1}{\varpi^{*}[s]}$
for $s\in[0,\epsilon]$. Since the middle part of the last equation is well-
defined for all $s\in[0,s_{0})$ it follows by properties of power series that
this last equation holds for all $s\in[0,s_{0})$. ∎

The next two lemmas are more of less covered in Section 2 of Li et al. [12],
though we provide a proof since [12] has a different starting definition of
MBIs and uses different notation.

######  Lemma 5.2.

Define $\Lambda:[0,s_{0})\to[0,\infty)$ by
$\Lambda(s)=\int_{0}^{s}\frac{1}{\psi(u)}\mathrm{d}u$. The function $\Lambda$
is bijective and so the inverse $\Lambda^{-1}:[0,\infty)\to[0,s_{0})$ is well-
defined. Then, for $x\in E$, $t\geq 0$ and $s\in[0,s_{0})$, the GF of the MBI
$(X,\mathbb{P})$ is given by,

(5.4)
$\mathbb{E}_{x}\left[s^{X_{t}}\mathds{1}_{\\{t<\zeta\\}}\right]=\Psi_{t}(s)^{x}\Phi_{t}(s),\quad$

where

(5.5) $\displaystyle\Psi_{t}(s)=$ $\displaystyle\Lambda^{-1}(t+\Lambda(s))$
(5.6) $\displaystyle\Phi_{t}(s)=$
$\displaystyle{\pi^{*}[s]}{\varpi^{*}\left[\Psi_{t}(s)\right]}.$

###### Proof.

Note that $\Lambda$ is bijective since it is increasing, continuous,
$\Lambda(0)=0$ and, as $\psi$ is a power series with $\psi(s_{0})=0$,
$\lim_{s\uparrow s_{0}}\Lambda(s)=\infty$. The first part of the proof of
Theorem 2.2 in [12] in combination with p.83 in [2] shows that, for $x,y\in E$
and $t\geq 0$, $P_{t}(x,y)$ is the unique nonnnegative solution to the forward
equations satisfying $\sum_{y\in E}P_{t}(x,y)\leq 1$. Let
$\widetilde{P}_{t}(x,y)=\left(P_{t}(0,\cdot)\star
P_{t}^{\circ}(x,\cdot)\right)(y)$ where $P_{t}^{\circ}(x,y)$ is the transition
kernel associated with the MBP with parameters $(\alpha,\mu,p)$ whose Q-matrix
we denote by $Q^{\circ}$. Clearly, $\widetilde{P}_{t}(x,y)\geq 0$ and
$\sum_{y\in E}\widetilde{P}_{t}(x,y)\leq 1$. Since $P_{t}$, respectively
$P^{\circ}_{t}$, satisfies the forward equations associated with $Q$,
respectively $Q^{\circ}$ (see e.g. Section 2.8 in [13]), we have

$\begin{split}\frac{d}{dt}\widetilde{P}_{t}(x,y)=&\sum_{l=0}^{y+1}\sum_{k=0}^{y+1}P_{t}(0,l)P_{t}^{\circ}(x,k)\left(Q(l,y-k)+Q^{\circ}(k,y-l)\right)\\\
=&\sum_{l=0}^{y+1}\sum_{k=0}^{y+1}P_{t}(0,l)P_{t}^{\circ}(x,k)Q(l+k,y)\\\
=&\sum_{m=0}^{y+1}\sum_{l=0}^{m}P_{t}(0,l)P_{t}^{\circ}(x,m-l)Q(m,y)\\\
=&\sum_{m\in E}\widetilde{P}_{t}(x,m)Q(m,y),\end{split}$

where the second equality is due to (5.1) and where we have set
$Q(x,-1)=Q^{\circ}(x,-1)=0$. So $\widetilde{P}_{t}$ is a solution to the
forward equations and therefore $P_{t}=\widetilde{P}_{t}$. From Theorem 3.1(2)
in Anderson or Theorem V.4.1 in Harris [9]222In these references it is assumed
that $p=0$ but the arguments go through for $p>0$., we get that
$P_{t}^{\circ}$ satisfies the branching property
$P_{t}^{\circ}(x,y)=P_{t}^{\circ}(1,\cdot)^{\star x}(y)$ where $\star x$
denotes $x$-fold convolution. Combining this with
$P_{t}(x,y)=\left(P_{t}(0,\cdot)\star P_{t}^{\circ}(x,\cdot)\right)(y)$, we
see that (5.4) holds with $\Psi_{t}(s)=P_{t}^{\circ}(1,\cdot)^{*}[s]$ and
$\Phi_{t}(s)=P_{t}(0,\cdot)^{*}[s]$. Since $P_{t}$, respectively
$P^{\circ}_{t}$, also satisfies the backward equations associated with $Q$,
respectively $Q^{\circ}$, we have with the help of Fubini, for $s\in[0,1]$ and
$t\geq 0$,

$\begin{split}P_{t}^{\circ}(1,\cdot)^{*}[s]=&s+\int_{0}^{t}\sum_{z\in
E}Q^{\circ}(1,z)P^{\circ}_{u}(z,\cdot)^{*}[s]du\\\
P_{t}(0,\cdot)^{*}[s]=&1+\int_{0}^{t}\sum_{z\in
E}Q(0,z)P_{u}(z,\cdot)^{*}[s]du.\end{split}$

So by taking derivatives in $t$ and using (5.2), we get, for $s\in[0,1]$ and
$t\geq 0$,

(5.7) $\displaystyle\frac{d}{dt}\Psi_{s}(t)=$
$\displaystyle\psi(\Psi_{t}(s)),$ $\displaystyle\Psi_{s}(0)=s,$ (5.8)
$\displaystyle\frac{d}{dt}\Phi_{s}(t)=$
$\displaystyle-\phi(\Psi_{t}(s))\Phi_{s}(t),$ $\displaystyle\Phi_{s}(0)=1.$

Since $\psi$ is locally Lipschitz on $(-1,1)$ and, for any $s\in[0,s_{0})$,
the right-hand side of (5.5) solves the ODE (5.7) for all $t\geq 0$, it
follows by basic ODE theory that (5.5) holds. Further, the linear ODE (5.8)
has a unique solution and via (5.7) and Lemma 5.1 the right-hand side of (5.5)
solves (5.8) for all $t\geq 0$ for any given $s\in[0,s_{0})$, which implies
that (5.6) holds. ∎

######  Lemma 5.3.

Tbe nonnegative measure $\pi$ defined in Lemma 5.1 is a stationary measure for
$(X,\mathbb{P})$.

###### Proof.

We need to prove $\sum_{x\in E}\pi(x)P_{t}(x,y)=\pi(y)$, for all $y,t\geq 0$
and $t\geq 0$. We do this by showing the GFs in $y$ of both sides are equal.
Fix $t\geq 0$ and $s\in[0,s_{0})$. We have by respectively Tonelli, (5.4),
(5.5) (which implies $\Psi_{t}(s)\in[0,s_{0})$) and (5.6) (recalling
$\pi^{*}[s]=1/\varpi^{*}[s]$),

$\begin{split}\sum_{y\geq 0}\sum_{x\in E}\pi(x)P_{t}(x,y)s^{y}=&\sum_{x\in
E}\pi(x)\sum_{y\geq 0}P_{t}(x,y)s^{y}=\sum_{x\in
E}\pi(x)\Psi_{t}(s)^{x}\Phi_{t}(s)=\pi^{*}[\Psi_{t}(s)]\Phi_{t}(s)=\pi^{*}[s].\end{split}$

∎

Let us mention the following about transience and recurrence for MBIs. If
$\beta=q=0$, i.e. $(X,\mathbb{P})$ is an MBP, then $0$ is an absorbing state
and thus recurrent whereas all other states are transient (as $0$ is
attainable from any other state). If $(X,\mathbb{P})$ is not an MBP, then
either $q>0$ in which case all states are transient or $\beta>0$ and $q=0$ in
which case the state space is irreducible and so either all states are
transient or all states are recurrent. In the latter case results on whether
the state space is transient or recurrent can be found in Section 3 in [12].
So there are three possible cases: either all states are transient, all states
are recurrent or $(X,\mathbb{P})$ is an MBP. The following theorem provides a
full set of fluctuation identities for MBIs.

######  Theorem 5.4.

Assume $(X,\mathbb{P})$ is an MBI with parameters
$(\alpha,\mu,p,\beta,\nu,q)$. Let for $x,y\geq 0$,

$\mathscr{H}^{[y}(x)=\sum_{l=0}^{y-x-1}\pi(y-x-1-l)\frac{(W\star\varpi)(l)}{l+x+1},$

where we understand that the sum equals $0$ if $y\leq x$. Then for $x\geq 0$
and $s\in[0,s_{0})$,

$\sum_{y\geq
0}\mathscr{H}^{[y}(x)s^{y}=\int_{0}^{s}\frac{v^{x}}{\psi(v)}\exp\left(\int_{v}^{s}\frac{\phi(u)}{\psi(u)}du\right)dv.$

If $\beta=q=0$, then $\mathscr{H}^{[y}(x)=\frac{W(y-x-1)}{y}$ where we set
$W(z)=0$ for $z\leq-1$. We further have the following.

  1. (1)

If all states of $(X,\mathbb{P})$ are transient, then the resolvent of
$(X,\mathbb{P})$ is given by, for $x,y\geq 0$,

$\begin{split}G(x,y)=\pi(y)\int_{0}^{s_{0}}\frac{v^{x}}{\psi(v)}\exp\left(-\int_{0}^{v}\frac{\phi(u)}{\psi(u)}du\right)dv-\mathscr{H}^{[y}(x),\\\
\end{split}$

whereas if $(X,\mathbb{P})$ is an MBP, then, for $x\geq 0$ and $y\geq 1$,

$\begin{split}G(x,y)=\frac{1}{y}\left(s_{0}^{x}W(y-1)-W(y-x-1)\right).\end{split}$

  2. (2)

If all states are transient, then for $x,y\geq 0$,

$\mathbb{P}_{x}(T_{y}<\zeta)=\frac{\pi(y)\int_{0}^{s_{0}}\frac{v^{x}}{\psi(v)}\exp\left(-\int_{0}^{v}\frac{\phi(u)}{\psi(u)}du\right)dv-\mathscr{H}^{[y}(x)}{\pi(y)\int_{0}^{s_{0}}\frac{v^{y}}{\psi(v)}\exp\left(-\int_{0}^{v}\frac{\phi(u)}{\psi(u)}du\right)dv},$

whereas if $(X,\mathbb{P})$ is an MBP, then for $x\geq 0$,

$\mathbb{P}_{x}(T_{y}<\zeta)=\begin{cases}\frac{s_{0}^{x}W(y-1)-W(y-x-1)}{s_{0}^{y}W(y-1)}&\text{if
$y\geq 1$},\\\ s_{0}^{x}&\text{if $y=0$}.\end{cases}$

Clearly, $\mathbb{P}_{x}(T_{y}<\zeta)=1$ for all $x,y\geq 0$ if all states are
recurrent.

  3. (3)

In any case, for $0\leq a\leq x\leq b-1$,

$\begin{split}\mathbb{P}_{x}(T_{a}<T_{[b}\wedge\zeta)=&\frac{\mathscr{H}^{[b}(x)}{\mathscr{H}^{[b}(a)},\\\
\mathbb{P}_{x}(T_{(a,b)^{\complement}}<\zeta)=&1-\sum_{z=a+1}^{b-1}(pz+q)\left(\frac{\mathscr{H}^{[b}(x)}{\mathscr{H}^{[b}(a)}\mathscr{H}^{[z}(a)-\mathscr{H}^{[z}(x)\right),\\\
\mathbb{P}_{x}(T_{[b}<\zeta)=&\begin{cases}1-\sum_{z=0}^{b-1}(pz+q)\left(\frac{\pi(z)}{\pi(b)}\mathscr{H}^{[b}(x)-\mathscr{H}^{[z}(x)\right)&\text{if
$\beta+q>0$},\\\
1+p\sum_{z=0}^{b-x-2}W(z)-\frac{W(b-x-1)}{W(b-1)}\left(1+p\sum_{z=0}^{b-2}W(z)\right)&\text{if
$\beta=q=0$}.\end{cases}\end{split}$

###### Proof.

We have for $x\geq 0$ and $s\in[0,s_{0})$,

(5.9) $\begin{split}\sum_{y\geq 0}\mathscr{H}^{[y}(x)s^{y}=\sum_{y\geq
x+1}\mathscr{H}^{[y}(x)s^{y}=&\pi^{*}[s]\sum_{l\geq
0}\frac{(W\star\varpi)(l)}{l+x+1}s^{x+1+l}\\\
=&\pi^{*}[s]\int_{0}^{s}\sum_{l\geq 0}(W\star\varpi)(l)v^{x+l}dv\\\
=&\pi^{*}[s]\int_{0}^{s}v^{x}W^{*}[v]\varpi^{*}[v]dv,\end{split}$

where we used Fubini for the penultimate equality which is applicable because
clearly $|\varpi(l)|\leq\pi(l)$ for all $l\geq 0$ by their definitions. So the
first statement follows from Lemma 5.1. If $\beta=q=0$, then by Lemma 5.1 we
easily see that $\pi(0)=\varpi(0)=1$ and $\pi(l)=\varpi(l)=0$ for $l\geq 1$
and consequently $\mathscr{H}^{[y}(x)=\frac{W(y-x-1)}{y}$. Next we prove item
(1). Assume all states are transient so that $G(x,y)<\infty$ for all $x,y\geq
0$. For $x\geq 0$ and $s\in[0,s_{0})$ we have by Lemma 5.2 and the change of
variables $t=\Lambda(v)-\Lambda(s)$,

$\begin{split}\sum_{y\geq
0}G(x,y)s^{y}=\int_{0}^{\infty}P_{t}(x,y)s^{y}dt=&\int_{0}^{\infty}\left(\Lambda^{-1}(t+\Lambda(s))\right)^{x}\pi^{*}[s]\varpi^{*}\left[\Lambda^{-1}(t+\Lambda(s))\right]dt\\\
=&\pi^{*}[s]\int_{s}^{s_{0}}\frac{v^{x}\varpi^{*}[v]}{\psi(v)}\mathrm{d}v\\\
=&\pi^{*}[s]\int_{0}^{s_{0}}\frac{v^{x}\varpi^{*}[v]}{\psi(v)}\mathrm{d}v-\pi^{*}[s]\int_{0}^{s}\frac{v^{x}\varpi^{*}[v]}{\psi(v)}dv.\end{split}$

So the expression for $G(x,y)$ follows by (5.9) and Lemma 5.1 in the case
where all states are transient. Now assume $(X,\mathbb{P})$ is an MBP, i.e.
$\beta=q=0$. Let, for $r>0$, $(X,\mathbb{P}^{(r)})$ be the MBI with parameters
$(\alpha,\mu,p,0,\nu,r)$. Then via (5.4), the transition kernel of
$(X,\mathbb{P}^{(r)})$ is given by $P^{(r)}_{t}(x,y)=e^{-rt}P_{t}(x,y)$,
$x,y\in E$, $t\geq 0$. Moreover, by the same reasoning as for (2.2), for any
stopping time $T$,

(5.10)
$\mathbb{P}_{x}^{(r)}(\Gamma\cap\\{T<\infty\\})=\mathbb{E}_{x}\left[e^{-rT}\mathds{1}_{\Gamma\cap\\{T<\infty\\}}\right],\quad\Gamma\in\mathcal{F}_{T},\
x\geq 0.$

Let $G^{(r)}(x,y)=\int_{0}^{\infty}P_{t}^{(r)}(x,y)dt$ be the resolvent of the
MBI with parameters $(\alpha,\mu,p,0,\nu,r)$ and denote by $\pi_{r}$ and
$\varpi_{r}$ the measures defined by $\pi_{r}(0)=\varpi_{r}(0)=1$ and, for
$k\geq 0$, $\pi_{r}(k+1)=\frac{r}{k+1}(\pi_{r}\star W)(k)$ and
$\varpi_{r}(k+1)=-\frac{r}{k+1}(\varpi_{r}\star W)(k)$. Since for the MBI with
parameters $(\alpha,\mu,p,0,\nu,r)$ all states are transient as $r>0$, we can
use the formula for the resolvent that we just proved to deduce for $x,y\geq
0$,

$G^{(r)}(x,y)=\pi_{r}(y)\int_{0}^{s_{0}}\frac{v^{x}}{\psi(v)}\exp\left(-\int_{0}^{v}\frac{r}{\psi(u)}du\right)dv-\sum_{l=0}^{y-x-1}\pi_{r}(y-x-1-l)\frac{(W\star\varpi_{r})(l)}{l+x+1}.$

It is easy to prove by induction that $\lim_{r\downarrow
0}\pi_{r}(y)=\lim_{r\downarrow 0}\varpi_{r}(y)=0$ and $\lim_{r\downarrow
0}\frac{\pi_{r}(y)}{r}=\frac{W(y-1)}{y}$ for $y\geq 1$ and further by an
integration by parts, for $x\geq 0$,

(5.11) $\lim_{r\downarrow
0}r\int_{0}^{s_{0}}\frac{v^{x}}{\psi(v)}\exp\left(-\int_{0}^{v}\frac{r}{\psi(u)}du\right)dv\\\
=\lim_{r\downarrow
0}\left(1+\int_{0}^{s_{0}}xv^{x-1}\exp\left(-\int_{0}^{v}\frac{r}{\psi(u)}du\right)dv\right)=s_{0}^{x}.$

So by the monotone convergence theorem we deduce, for all $x\geq 0$ and $y\geq
1$,

$\begin{split}G(x,y)=\lim_{r\downarrow 0}G^{(r)}(x,y)=&\lim_{r\downarrow
0}\frac{\pi_{r}(y)}{r}r\int_{0}^{s_{0}}\frac{v^{x}}{\psi(v)}\exp\left(-\int_{0}^{v}\frac{r}{\psi(u)}du\right)dv-\frac{W(y-x-1)}{y}\\\
=&\frac{s_{0}^{x}W(y-1)}{y}-\frac{W(y-x-1)}{y},\end{split}$

which finishes the proof of item (1).

For items (2) and (3) we use the results of Section 3. First assume that all
states are transient. Then $\beta>0$ or $q>0$ and thus $\pi$ from Lemma 5.1 is
a positive (not merely nonnegative) measure and by Lemma 5.3 it is a
stationary and thus excessive measure for $(X,\mathbb{P})$. So we can take
$\pi$ as our reference measure. We further let the reference point be
$\mathfrak{o}=0$. As $\pi$ is stationary, $\mathbb{P}_{x}(t<\zeta)=\sum_{y\geq
0}\widehat{P}_{t}(x,y)=\frac{\pi(x)P_{t}(x,y)}{\pi(y)}=1$ for all $x\geq 0$
and any $t\geq 0$, so the dual Markov chain $(X,\widehat{\mathbb{P}})$ has
infinite lifetime. Because all states are transient for the dual and the dual
is skip-free upwards, it follows that
$\widehat{\mathbb{P}}_{y}(T_{b}<\zeta)=1$ for all $y\leq b$. Consequently, by
the analogue of (3.1) for the dual we deduce $\widehat{H}(y)=1$ for all $y\geq
0$. By the definition of $H$ in Theorem 3.1(1), $H(x)=\frac{G(x,0)}{G(0,0)}$.
Then by (3.2) and the expression for $G(x,y)$ from item (1),

$\mathbf{H}^{[y}(x)=H(x)-\frac{G(x,y)}{\mathfrak{C}\widehat{H}(y)\pi(y)}=\frac{1}{\mathfrak{C}}\left(G(x,0)-\frac{G(x,y)}{\pi(y)}\right)=\frac{\mathscr{H}^{[y}(x)}{\mathfrak{C}\pi(y)},\quad
x,y\geq 0.$

The identities in items (2)-(3) then immediately follow from the identities in
Theorem 3.3 in the case where all states are transient. When not all states
are transient for $(X,\mathbb{P})$, we consider for $r>0$, the MBI
$(X,\mathbb{P}^{(r)})$ for which all states in $E$ are transient and so for
which the identities in items (2)-(3) apply. Then via (5.10) and the monotone
convergence theorem we can take $r\downarrow 0$ on both sides of the
identities for $(X,\mathbb{P}^{(r)})$ to get the identities for
$(X,\mathbb{P})$ using that, with the obvious notation, $\lim_{r\downarrow
0}\pi_{r}(x)=\pi(x)$ and $\lim_{r\downarrow
0}\mathscr{H}^{[y}_{r}(x)=\mathscr{H}^{[y}(x)$ for $x\geq 0$ and using in the
case where $\beta=q=0$, $\lim_{r\downarrow
0}\frac{\pi_{r}(y)}{r}=\frac{W(y-1)}{y}$ for $y\geq 1$ and (5.11). ∎

######  Remark 5.5.

Regarding the role of the killing parameters $p$ and $q$, if $(X,\mathbb{P})$
is an MBI with parameters $(\alpha,\mu,p,\beta,\nu,q)$ and $(X,\mathbf{P})$ an
MBI with parameters $(\alpha,\mu,0,\beta,\nu,0)$, then following the arguments
leading up to (2.2) one can show for any stopping time $T$,

$\mathbb{P}_{x}(\Gamma\cap\\{T<\infty\\})=\mathbf{E}_{x}\left[e^{-qT-p\int_{0}^{T}X_{s}ds}\mathds{1}_{\Gamma\cap\\{T<\infty\\}}\right],\quad\Gamma\in\mathcal{F}_{T},\
x\geq 0,$

where $\mathbf{E}_{x}$ is the expectation operator associated with
$\mathbf{P}_{x}$. So all the identities in Theorem 5.4 can be rephrased in
terms of the MBI $(X,\mathbf{P})$ for which the killing parameters $p$ and $q$
are $0$. As an example,

$\mathbf{E}_{x}\left[e^{-qT_{a}-p\int_{0}^{T_{a}}X_{s}ds}\mathds{1}_{\\{T_{a}<T_{[b}\wedge\zeta\\}}\right]=\mathbb{P}_{x}(T_{a}<T_{[b}\wedge\zeta)=\frac{\mathscr{H}^{[b}(x)}{\mathscr{H}^{[b}(a)},\quad
0\leq a\leq x\leq b-1.$

A similar remark holds of course for the role of the killing parameter $p$ for
downward skip-free compound Poisson processes considered in Section 4.

Most of the identities in Theorem 5.4 are new including in the sense that no
analogues have appeared in the literature for the discrete-time analogue of
MBIs, i.e. Galton-Watson processes, or the continuous-state space analogue of
MBIs, i.e. continuous-state branching processes with immigration (CBIs). For
the special case of MBPs some of these identities are covered in the
literature because the Lamperti transform tells us that an MBP is a time-
changed skip-free downward compound Poisson process stopped at hitting 0 and
killed at a constant rate, see Theorem III.11.3 in [4] for the Lamperti
transform for nondecreasing MBPs. This implies that, in the MBP case, the
first two identities in Theorem 5.4(3) and the identity in Theorem 5.4(2) for
$y\leq x$ are identical to those for skip-free downward compound Poisson
processes, which is confirmed by comparing Theorems 4.1 and 5.4, see also
Theorem 12.8 and Corollary 12.9 in [11] for the continuous-state space case.
With immigration the only identity that, to our knowledge, has appeared before
in the literature is the one for the downwards hitting probability, which from
Theorem 5.4(2) reads, in the case where all states are transient,

(5.12)
$\mathbb{P}_{x}(T_{a}<\zeta)=\frac{\int_{0}^{s_{0}}\frac{v^{x}}{\psi(v)}\exp\left(-\int_{0}^{v}\frac{\phi(u)}{\psi(u)}du\right)dv}{\int_{0}^{s_{0}}\frac{v^{a}}{\psi(v)}\exp\left(-\int_{0}^{v}\frac{\phi(u)}{\psi(u)}du\right)dv},\quad
0\leq a\leq x.$

Namely, the analogue of this identity for CBIs has been derived by Duhalde et
al. in [8] and the MBI case can be found in [17], see also [5] for the case
$\beta=0$. We remark that given the input parameters
$(\alpha,\mu,p,\beta,\nu,q)$ and given $y\geq 0$ the object
$\mathscr{H}^{[y}(x)$, $0\leq x\leq y-1$, can be computed exactly in a finite
number of steps via the recursions in Lemma 5.1. On the other hand for e.g.
the downwards hitting probability (5.12) one needs to evaluate the branching
and immigration mechanisms $\psi$ and $\phi$ which might not always have a
fully explicit form depending on the input parameters. Since
$\mathbb{P}_{x}(T_{a}<\zeta)=\lim_{b\to\infty}\mathbb{P}_{x}(T_{a}<T_{[b}\wedge\zeta)$
one can alternatively approximate $\mathbb{P}_{x}(T_{a}<\zeta)$ by evaluating
$\mathbb{P}_{x}(T_{a}<T_{[b}\wedge\zeta)$ for a sufficiently large $b$ which
avoids needing to compute $\psi$ and $\phi$ or the integrals in (5.12). In
Theorem 5.4 we see that the expressions simplify for the special case of MBPs.
This is also the case for another subclass of MBIs.

######  Example 5.6.

Assume the parameters $(\alpha,\mu,p,\beta,\nu,q)$ of the MBI are such that
$0<m\leq 1+\frac{p}{\alpha}$, $\beta=\alpha m$,
$\nu(j)=\frac{(j+1)\mu(j+1)}{m}$ for $j\geq 1$ and $q=p-\alpha(m-1)$, where
$m:=\sum_{j\geq 2}j\mu(j)$. Then $\phi(s)=-\psi^{\prime}(s)$ for all $s\geq 0$
and consequently, for $s\in[0,s_{0})$ and $y\geq 0$, via Lemma 5.1 and Theorem
5.4,

$\int_{0}^{s_{0}}\frac{v^{x}}{\psi(v)}\exp\left(-\int_{0}^{v}\frac{\phi(u)}{\psi(u)}du\right)dv=\frac{1}{\psi(0)}\frac{s_{0}^{x+1}}{x+1},\quad\pi^{*}[s]=\frac{\psi(0)}{\psi(s)},\quad\sum_{y\geq
0}\mathscr{H}^{[y}(x)s^{y}=\frac{1}{\psi(s)}\frac{s^{x+1}}{x+1}.$

This yields, for $x,y\geq 0$,

$\pi(x)=\psi(0)W(x),\quad\mathscr{H}^{[y}(x)=\frac{W(y-x-1)}{x+1}.$

## References

  * [1]
  * [2] W. Anderson.  _Continuous-Time Markov Chains: An Applications-Oriented Approach_.  Springer-Verlag, New York, 1991.
  * [3] S. Asmussen.  _Applied probability and queues_ , volume 51 of _Applications of Mathematics (New York)_.  Springer-Verlag, New York, second edition, 2003.  Stochastic Modelling and Applied Probability.
  * [4] K. B. Athreya and P. E. Ney.  Branching processes, Springer-Verlag, New York-Heidelberg, 1972\.
  * [5] F. Avram, P. Patie and J. Wang,  Purely excessive functions and hitting times of continuous-time branching processes,  _Methodol. Comput. Appl. Probab._ , 21:391-399, 2019.
  * [6] M. Choi and P. Patie.  Skip-free Markov chains.  _Trans. Amer. Math. Soc._ , 371(10), 7301-7342, 2019.
  * [7] R. A. Doney.  Fluctuation Theory for Lévy Processes. Ecole d’Eté de Probabilités de Saint-Flour XXXV-2005. Springer, 2007.
  * [8] Xan Duhalde, Clément Foucart, and Chunhua Ma.  On the hitting times of continuous-state branching processes with immigration.  Stochastic Processes and their Applications, 124(12):4182–4201, 2014.
  * [9] T. E. Harris. _The Theory of Branching Processes_.  Springer-Verlag, Berlin; Prentice-Hall, Inc., Englewood Cliffs, N.J, 1963.
  * [10] J. G. Kemeny, J. L. Snell, and A. W. Knapp. _Denumerable Markov chains_.  Springer-Verlag, New York-Heidelberg-Berlin, second edition, 1976. With a chapter on Markov random fields, by David Griffeath, Graduate Texts in Mathematics, No. 40.
  * [11] A. E. Kyprianou.  Introductory lectures on fluctuations of Lévy processes with applications.  Second Edition. Springer-Verlag, Berlin, 2014.
  * [12] J. Li, A. Chen and A. G. Pakes.  Asymptotic properties of the Markov branching process with immigration _J. Theoret. Probab._ , 25:122–143, 2012.
  * [13] J. R. Norris.  Markov chains.  Cambridge University Press, Cambridge, 1997.
  * [14] B. A. Sevast’yanov.  Limit theorems for branching stochastic processes of special form. _Theory Probab. Appl._ , 2(3), 321–331, 1957.
  * [15] F. Spitzer. Principles of random walk. _The University Series in Higher Mathematics. D. Van Nostrand Co., Inc., Princeton, N.J.-Toronto-London_ , 1964.
  * [16] M. Vidmar. Fluctuation Theory for Upwards Skip-Free Lévy Chains _Risks_ , 6(3), 102, 2018.
  * [17] M. Vidmar. Some harmonic functions for killed Markov branching processes with immigration and culling, _arXiv 1908.04714v5_ , 2021.