Hostname: page-component-cb9f654ff-5jtmz Total loading time: 0 Render date: 2025-08-15T10:27:07.313Z Has data issue: false hasContentIssue false

Hitting probabilities of constrained random walks representing tandem networks

Published online by Cambridge University Press:  13 August 2025

Ali Devin Sezer*
Affiliation:
Institute of Applied Mathematics, Middle East Technical University, Ankara 06800, Turkey
Rights & Permissions [Opens in a new window]

Abstract

We develop anapproximation for the buffer overflow probability of a stable tandem network in dimensions three or more. The overflow event in terms of the constrained random walk representing the network is the following: the sum of the components of the process hits n before hitting 0. This is one of the most commonly studied rare events in the context of queueing systems and the constrained processes representing them. The approximation is valid for almost all initial points of the process and its relative error decays exponentially in n. The analysis is based on an affine transformation of the process and the problem; as $n\rightarrow \infty$ the transformed process converges to an unstable constrained random walk. The approximation formula consists of the probability of the limit unstable process hitting a limit boundary in finite time. We give an explicit formula for this probability in terms of the utilization rates of the nodes of the network.

Information

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2025. Published by Cambridge University Press.

1. Introduction

Tandem networks, where each customer must visit a sequence of servers, are one of the simplest forms of Jackson networks (see Figure 1). The embedded random walk X corresponding to a Jackson network represents the number of customers at each node of the network waiting for service right after arrivals and service completions. Each step of X corresponds to either an arrival to or a service completion of the network. Since the number of customers is always nonnegative, X is constrained to remain in ${\mathbb Z}_+^d$ where ${\mathbb Z}_+ =\{0,1,2,3,...\}$ and dimension d is the number of service nodes (or queues) in the system. This paper concerns embedded constrained random walks of tandem networks (i.e., tandem walks) with $d\ge 3$. We assume the network (equivalently, the corresponding random walk X) to be stable, i.e., we assume that the arrival rate to the system is less than the service rate at each node (see (2) and (5)). We are interested in the probability pn that the sum of the components of X equals n before X hits the origin given that X starts from an initial point $x \in {\mathbb Z}_+^d$.

Figure 1. d tandem queues.

The probability pn has the following physical interpretation: suppose that the customers/packets in the system are stored in a system-wide joint buffer of size n − 1. Then the first time the components of X sum to n is also the first time the system experiences a buffer overflow event and pn is the probability that a buffer overflow occurs before the system empties (i.e., before X hits the origin) given that the process starts from the initial position $X_0 = x$ (note that pn is a function of the initial position x; when no emphasis is needed on the initial position, we will simply write pn). If we divide the trajectory of X into independent cycles consisting of times spent between X’s hitting times to the origin, pn is the probability that the buffer overflows in the current cycle. The probability pn has a geometric interpretation as well: it is the probability that X hits the exit boundary consisting of all points in ${\mathbb Z}_+^d$ whose coordinates sum to n before hitting the origin. Stability of X implies that pn decays exponentially in n, hence our event of interest (the probability of hitting the exit boundary before hitting the origin) is rare for n large.

Computation/approximation of pn turns out to be a nontrivial problem (see the literature review below); the difficulty arises from the multidimensionality of the problem and the constraining boundaries of the process. A hallmark of Jackson networks is the product form formulas for their stationary distributions. The work [Reference Sezer65] develops an explicit approximation formula similar to these product form formulas for pn in the case of the two-dimensional tandem walk (the exact formula is (18) with $y=(n-x(1),x(2))$ and x(i) are the components of $x \in {\mathbb Z}_+^2$) and proves that the relative error of the approximation decays exponentially in n. The analysis in [Reference Sezer65] is based on an affine transformation of X and pn; the goal of the present work is to extend these results and analysis to dimensions $d \ge 3$.

The affine transformation approach consists of changing the coordinate system by an affine transformation so that the origin of the new coordinate system is on the exit boundary; letting $n\rightarrow \infty$ gives a limit process Y and a limit probability (see Figure 2 for an example in two dimensions). The limit process Y is unconstrained in its first coordinate which makes it unstable; the limit probability is the probability of the process Y hitting the limit exit boundary in finite time. The main results of the present work are:

  1. (1) Theorem 2.1 which proves that the limit probability approximates pn with decaying exponential error for almost all initial points in scaled coordinates,

  2. (2) Definition 4.7 of a harmonic system (consisting of a regular graph and a system of equations associated with its nodes and edges) and Theorem 4.8 which constructs a harmonic function given the solution of a harmonic system.

  3. (3) Definitions and results in Section 5 culminating with Theorem 5.7 giving an explicit formula for the limit probability. Definitions and results in Section 5 preceding Theorem 5.7 define a specific class of harmonic systems and their solutions for the d dimensional tandem walk. By Theorem 4.8 these solutions give harmonic functions for the limit process. Theorem 5.7 proves that the right linear combination of these functions gives us the limit probability.

  4. (4) Corollary 6.1 which derives the exponential decay rate of $p_n(x_n)$ for any sequence xn such that $x_n/n \rightarrow x$ with $\sum_{j=1}^d x(j) \le 1$ (i.e., the large deviations analysis of pn).

Figure 2. Derivation of the limit problem via an affine transformation in two dimensions.

More detailed statement of our results and a summary of our analysis are provided in subsection 2.1 below.

In addition to [Reference Sezer65], the works [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5] employ the affine transformation approach to derive approximation results similar to those derived in the present work: [Reference Ünlü and Sezer1] treats the case of the constrained simple random walk in two dimensions and [Reference Başoğlu Kabran and Sezer5] treats the Markov modulated tandem random walk in two dimensions. We next discuss the novelties of the present work as compared to [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65]. As already noted, all of these works treat only two dimensions; they use the two-dimensional structure of the problem at every step of their analysis. Let us begin with the derivation of the formula for the limit probability. A function is said to be harmonic with respect to a Markov process if its value at a given point equals the average of its values at the neighboring points where the average is taken with respect to the jump probabilities of the process (see [Reference Revuz55, Definition 1.1, page 40] or (20) in our setting); sub and super harmonic functions are defined similarly. In the two-dimensional nonmodulated case (treated in [Reference Ünlü and Sezer1, Reference Sezer65]) the formula for the limit probability is constructed from the roots of two quadratic equations in two variables (upon fixing one of the variables, the determination of the solutions reduces to the quadratic formula), one has to check if the resulting functions are harmonic with respect to a limit process and whether an appropriate linear combination of them has the right value on the limit boundary. Since the resulting functions have at most two terms, only two functions are needed in the linear combination and the limit boundary is a line (a single dimension). These make the above steps straightforward in [Reference Ünlü and Sezer1, Reference Sezer65]. In the modulated case treated in [Reference Başoğlu Kabran and Sezer5], the situation is similar, i.e., the arguments crucially make use of the two-dimensionality of the problem.

In d dimensions treated in the present work, the computation of the limit probability proceeds in two stages. In the first stage we introduce harmonic systems (Definition 4.7) for general Jackson networks: these are regular graphs whose nodes and edges define equations to be solved. We first prove that solutions of these systems give harmonic functions (Theorem 4.8); the corresponding stage and results are elementary in two dimensions and are not presented as a separate stage. In the second stage we introduce a specific class of harmonic systems for the tandem case whose graphs have up to $2^{d-1}$ nodes (see Section 5 and the graphs $G_{d,{\boldsymbol d}}$ defined in (127)-(129)). That 1) the solutions of the equations represented by these graphs exist 2) they have a linear combination that has the correct value on the limit boundary are nontrivial to prove (see Propositions 5.5 and Theorem 5.7). The two-dimensional case provides almost no insight into whether these arguments, structures and solutions exist in higher dimensions. Furthermore, the dimension d is now a variable itself; so a case by case analysis in fixed number of steps (as is possible in a fixed dimension) is not possible.

A difficulty in d dimensions is the lower dimensional boundaries of the state space. In the computation of the limit probability these boundaries come up when checking the harmonicity of the candidate limit function in d dimensions: the arguments must treat all $2^d-2$ boundaries of different dimensions whereas in 2 dimensions there is only a single one-dimensional boundary. The increasing complexity of the geometry of the problem is also reflected in the analysis of the relative error: a priori it is not clear that the relative error analysis for the two-dimensional case provided in [Reference Sezer65] would generalize to d dimensions. And indeed, we have not been able to produce a direct generalization. One of the reasons is as follows: in [Reference Sezer65] (as well as in [Reference Ünlü and Sezer1] and [Reference Başoğlu Kabran and Sezer5]) the harmonic functions that make up the limit approximating function are also used in the relative error analysis. These works are able to use these harmonic functions in the relative error analysis in two dimensions because in two dimensions they consist of a few terms. The number of terms grows exponentially in d resulting in functions with complex formulas (see (17); in comparison the same formula in two dimensions is (18)) and it is not clear how they can be used in the error analysis. For this reason the present work develops superharmonic functions with simpler structures that allow us to establish upper bounds on the relative error (see subsections 3.1, 3.2 and 3.3). We further comment on how dimension impacts the convergence argument in subsection 3.4.

The use of superharmonic functions in the convergence analysis leads to a third novelty in the current paper compared to [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65]: in their relative error analysis, all of these previous works make an assumption on the inequality of the jump rates. Such an assumption is no longer needed in the relative error analysis of the present work thanks to the use of the superharmonic functions mentioned above whose construction requires only stability and not the inequality assumption. We further comment on this improvement at the end of subsection 2.1.

Another novelty compared to [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65] is the more precise nature of the relative error analysis: these earlier works simply state that the relative error decays exponentially in n. In the present work, we give precise rates of convergence as a function of the initial position. Finally, as a corollary of our results we derive the large deviation decay rate of the buffer overflow probability for any initial point in the scaled domain (see the paragraph below or Section 6); this was not done in [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65]. These and further challenges and new features of the analysis compared to [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65] are discussed in detail using mathematical notation in subsection 2.1 below.

Next, we give a summary of other results and approaches from the current literature on the approximation of pn and similar quantities. As far as analytical approximations go, the currently available literature (excepting [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65] discussed above) focuses on large deviations (LD) analysis, i.e., on the computation of the decay rate of pn in n. This rate is computed for general Jackson networks in [Reference Glasserman and Kou32, Reference Ignatiouk-Robert35] as

(1)\begin{equation} \lim_{n\rightarrow \infty} -\frac{1}{n}\log p_n( x_n) =-\log\rho, \end{equation}

for $x=x_n$, $x_n/n \rightarrow 0$. As already noted above, a corollary of our approximation formula is the extension of this result to any sequence $x_n/n\rightarrow x\in A:=\{x\in\mathbb{R}_+^d,\sum_{j=1}^dx(j)\;\leq1\}$ (Corollary 6.1 in Section 6 and (23) below). We are not aware of other results in the currently available literature that gives LD decay rates for exit probabilities of constrained random walks in arbitrary dimension where the initial point is arbitrary (as is the case in Corollary 6.1).

To the best of our knowledge, beyond LD analysis, most of the currently available literature on the approximation of pn focuses on simulation techniques. Since pn is small for typical values of n, its efficient simulation requires variance reduction algorithms such as importance sampling (IS). The article [Reference Parekh and Walrand53] appears to be the first work on IS simulation of pn for two tandem queues and the exit boundary $\{x \in {\mathbb Z}_+^d: x(1) + x(2)= n\}$ (the same as the one studied in the present work); this work notes that the constrained dynamics of X causes static changes of measure implied by large deviations analysis to have poor performance. As a remedy, [Reference Parekh and Walrand53] introduces dynamic IS measures which depend on the position of X. The work [Reference Glasserman and Kou32] notes that static changes of measure perform poorly for the exit boundary $\{x: x(1)+ x(2) = n \}$ for a range of parameter values. An IS change of measure to estimate pn is said to be asymptotically optimal if the variance of the corresponding IS estimator decays at twice the rate of the large deviation decay rate of $p_n.$ The work [Reference Dupuis, Sezer and Wang23] develops asymptotically optimal and dynamic IS changes of measure for this problem using subsolutions of the limit Hamilton Jacobi Bellman (HJB) equation arising in the LD analysis. Works [Reference Dupuis, Leder and Wang22, Reference Dupuis and Wang26, Reference Sezer60, Reference Sezer61] generalize these results to higher dimensions, more general dynamics and exit boundaries. A recent article on the construction of IS algorithms based on subsolutions of the limit HJB equations is [Reference Buijsrogge, de Boer and Scheinhardt10], which treats the approximation of pn for the tandem walk. All of these works assume that the initial position of X satisfies $x=x_n$ where $x_n/n \rightarrow 0$.

A general Jackson network is said to be stable if the average arrival rate to each node is less than the service rate of that node. A well known fact about the embedded random walk of such a network is that its stationary distribution is multivariate geometric where the parameter of each component is equal to the utilization rate of that node (this is what we refer to above as the “product form formulas” for the stationary distribution). Since the embedded random walk has a stationary distribution, it is reversible, i.e., when it is observed backward in time it is again a Markov chain and its transition probabilities are $K'(x,y) = K(y,x) \pi(y)/\pi(x)$ where K is the transition matrix of the forward process and π is the stationary distribution. A vein of research in the effective simulation of the hitting probability ${\mathbb{P}}_x(\tau_n \lt \tau_0)$ makes use of the reversed process. Heuristic simulation algorithms based on reversibility are developed in [Reference Nicola and Zaburnenko52]. The work [Reference Blanchet6] develops strongly efficient simulation algorithms using the reversed process; the analysis and the simulation algorithm are based on an exact representation of the hitting probability in terms of the reversed process with initial distribution taken to be the stationary distribution conditioned on the exit set; [Reference Blanchet6] assumes that the initial point x remains bounded as $n\rightarrow \infty$ (in particular $x_n/n \rightarrow 0$ as $n\rightarrow \infty$).

An approach to the approximation of $p_n(x_n)$, $x_n/n\rightarrow 0$, combining asymptotic analysis and simulation is given in [Reference McDonald47]; see Section 8 for a review of this work. There is a vast literature on the analysis and simulation of rare events similar to pn of constrained random walks in particular, and on the analysis of constrained random walks in general; see, e.g., [Reference Anantharam, Heidelberger and Tsoucas2, Reference Asmussen and Glynn3, Reference Blanchet6, Reference Blanchet, Glynn and Leder7, Reference Borovkov and Al’fredovich Mogul’skii9, Reference Collingwood, Foley and McDonald11Reference Dean and Dupuis19, Reference Dupuis and Ellis21, Reference Dupuis, Leder and Wang22, Reference Dupuis and Wang24, Reference Dupuis and Wang26Reference Frater, Lennon and Anderson30, Reference Ignatiouk, Malyshev and Scherbakov34Reference Juneja and Shahabuddin38, Reference Kobayashi and Miyazawa40Reference Louchard and Schott44, Reference Maier46, Reference Miretskiy, Scheinhardt and Mandjes48, Reference Miyazawa49, Reference Nicola and Zaburnenko52, Reference Randhawa and Juneja54, Reference Ridder56Reference Rubino and Tuffin58, Reference Sezer61, Reference Sezer62, Reference Yao66]. [Reference Sezer65, Section 6] gives a detailed review of a number of these works that are most relevant to the analysis in the present paper. For the reader’s convenience we reproduce an updated summary in Section 8. For further reviews we refer the reader to [Reference Ünlü and Sezer1] and [Reference Başoğlu Kabran and Sezer5]; the review in [Reference Ünlü and Sezer1] has a greater focus on simple random walks (i.e., random walks with increments ei and $-e_i$) and the one in [Reference Başoğlu Kabran and Sezer5] focuses on Markov modulated dynamics. An extensive review of the importance sampling literature that is related to the problem studied in this paper can be found in [Reference Blanchet and Lam8].

The next section provides the mathematical framework of the paper (i.e., definitions of the domains, boundaries, the processes, the stopping times and probabilities of interest) and lists the main results of the paper. The same section also provides a detailed mathematical outline of the ideas involved in the proofs. The error analysis linking the probability of interest to the limit probability is presented in Section 3; since the proofs can be involved due to the general dimension d of the problem, to give the reader some idea in a simpler setting, subsection 3.4 provides some of the arguments for the fixed dimension d = 2 and comments on how they change with $d\ge 3$.

Sections 4 and 5 are devoted to the computation of the limit probability: Section 4 introduces a framework that reduces the construction of harmonic functions with respect to the limit process to the construction of systems of equations (which we call harmonic systems) represented by graphs and their solutions. This reduction can be done for general Jackson networks and for complex valued harmonic functions without additional effort. For this reason Section 4 is presented in this generality. Note that the results in Section 4 do not construct any harmonic functions, they give one way of constructing them by constructing harmonic systems and solving them. Section 5 explicitly constructs harmonic systems for the tandem case and provides solutions for them. By the results of Section 4, these solutions define a class of harmonic functions for the limit process. Theorem 5.7, the final result in Section 5, presents an explicit formula for the limit probability as a linear combination of these harmonic functions.

Section 6 derives the large deviations decay rate of pn as a corollary of Theorems 2.1 and 5.7, for any initial point in the scaled domain. Section 7 provides two numerical examples on the approximation results obtained in Theorems 2.1 and 5.7. In the first example we take d = 4 and in the second example d = 14. In both examples we take n = 60. In the first example the numbers d and n are small enough to allow a precise numerical calculation of pn by iterating the harmonic equation that pn satisfies; the approximation of pn given by Theorems 2.1 and 5.7 are compared with the results of this calculation. For d = 14 and n = 60 the exact computation based on the iteration of the harmonic equation satisfied by pn is no longer feasible. For the comparison we use the large deviations approximation and an IS estimation based on the change of measure implied by the approximate formula given in Theorem 5.7.

Section 8 provides a detailed and comparative review of some of the works cited above. In Section 9 we further comment on our results and on future work. A list of symbols and notation used throughout the paper can be found after Section 9.

2. Definitions and statement of results

For an integer $d \ge 3$ let X be a random walk with independent and identically distributed increments $\{I_1,I_2,I_3,...\}$, $I_k \in {\mathcal V} \subset {\mathbb Z}^d$, constrained to remain in ${\mathbb Z}_+^d$, i.e., the sequence $\{I_k\}$ is independent of $X_0 \in {\mathbb Z}_+^d$ and

(2)\begin{equation} X_{k+1} = X_k + \pi(X_k,I_{k+1}),~~ k=0,1,2,3,..., \end{equation}

where

\begin{equation*} \pi: {\mathbb Z}^{d}_+ \times {\mathbb Z}^d \mapsto {\mathbb Z}^d_+,~~~ \pi(x,v) := \begin{cases} v, &\ \text{if } x +v \in {\mathbb Z}_+^d, \\ 0 , &\ \text{otherwise.} \end{cases} \end{equation*}

The random variables X 0, I 1, I 2,.... are assumed to be defined on a measurable pair $(\Omega, {\mathscr F})$ equipped with a family of probability measures, ${\mathbb{P}}_x$, $x \in {\mathbb Z}_+^d$, such that ${\mathbb{P}}_x( X_0 = x) =1.$ When the initial point plays no role we will simply write ${\mathbb{P}}.$ Define

\begin{equation*}{\mathscr F}_0:=\{\Omega,\varnothing\},\text{ }\text{ }{\mathscr F}_k:=\sigma(\{I_1,I_2,...,I_k\}),\text{ }\text{ } k \in \{1,2,3,...\}.\end{equation*}

The process X is adapted to the filtration ${\mathbb F}= \{{\mathscr F}_k, k=1,2,3,...\}.$ All of the processes appearing below are adapted to the same filtration.

We will use the function notation to indicate the components of a vector, e.g., $x=(x(1),x(2),...,x(d))$ for $x \in {\mathbb Z}_+^d.$ The constraining boundaries of X are

(3)\begin{equation} \partial_j = \{x \in {\mathbb Z}^d: x(j) = 0 \}, j \in \{1,2,3,...d\}. \end{equation}

We focus on the case when X represents d queues in tandem (a tandem network). This means that the set of possible increments of X are

(4)\begin{align} {\mathcal V} &= \{e_1, -e_1 + e_2,..., -e_j+e_{j+1},...,-e_{d-1}+e_d, -e_d \},\\ e_i(j) &= \begin{cases} 1, &\ \text{if } i = j,\\ \nonumber 0, &\ \text{otherwise,} \end{cases} \end{align}

$i,j=1,2,3,...,d$; $\{e_i,i=1,2,3,...,d\}$ are the unit vectors in ${\mathbb Z}^d$. The distribution of the increments is given as follows:

\begin{align*} {\mathbb{P}}(I_k = e_1) &= \lambda, \\ {\mathbb{P}}(I_k = e_{i+1}-e_{i}) &= \mu_i,~~~ i=1,2,3,...,d-1,\\ {\mathbb{P}}(I_k = -e_d) &= \mu_d, \end{align*}

where $\lambda + \sum_{i=1}^d \mu_i =1$ and $\lambda, \mu_1,\mu_2,...\mu_d \gt 0$ (if λ, µi are given as jump rates of a continuous time tandem network, we renormalize them so that they sum to 1 and we use the renormalized jump rates as the jump probabilities of the embedded random walk). We assume X to be stable:

(5)\begin{equation} \lambda \lt \min_{i=1}^d \mu_i \end{equation}

which implies

(6)\begin{equation}\rho_i:= \lambda/\mu_i,\;\;\rho:=\max_i\rho_i\;\leq1.\end{equation}

For

(7)\begin{equation} A_n := \left\{x \in {\mathbb Z}_+^d: \sum_{i=1}^d x(i) \le n \right\}, \end{equation}

and $\partial A_n := \{x \in {\mathbb Z}_+^d: \sum_{i=1}^d x(i) = n \}$ define the stopping time

(8)\begin{equation} \tau_n := \inf \{k \ge 0: X_k \in \partial A_n\}; \end{equation}

τn is the first time the sum of components of X equals n, in particular, τ 0 is the first time X hits the origin.

Our probability of interest is

(9)\begin{equation} p_n(x) := {\mathbb{P}}_x( \tau_n \lt \tau_0), \end{equation}

$x \in A_n$.

2.1. Summary of results and analysis

As already discussed in the introduction, we will approximate pn using a limit process Y, a limit boundary $\partial B$, the associated hitting time τ and the limit probability that τ is finite. The process Y differs from X only in the following ways: the 1) the first components of its jumps are reversed (i.e., the jump e 1 is replaced with $-e_1$ and the jump $-e_1 +e_2$ is replaced with $e_1 + e_2$) and 2) it is not constrained in its first component (see Figure 2 for an example in two dimensions). The formal definition is as follows. Let ${\mathcal I}_1 \in {\mathbb R}^{d\times d}$ be the diagonal matrix with diagonal entries ${\mathcal I}_1(1,1) =-1$, ${\mathcal I}_1(j,j) =1$, $j\in \{2,3,...,d\} $, and

(10)\begin{equation} T_n:{\mathbb Z}^d \mapsto {\mathbb Z}^d,~~~ T_n(x)= n e_1 + {\mathcal I}_1x. \end{equation}

Recall the definition (2) of the original walk X, in particular, Ik are the unconstrained increments of X; define $\Omega_Y := {\mathbb Z} \times {\mathbb Z}_+^{d-1}$,

(11)\begin{equation} J_k := {\mathcal I}_1 I_k,~~~ Y_{k+1} := Y_k + \pi_1(Y_k,J_{k+1}), \end{equation}

where

\begin{equation*} \pi_1(y,v) := \begin{cases} v, &\ \text{if } y +v \in \Omega_Y,\\ 0 , &\ \text{otherwise.} \end{cases} \end{equation*}

The initial point Y 0 is assumed to satisfy $Y_0 = y$ almost surely under ${\mathbb{P}}_y.$ Define

\begin{equation*} B := \left\{y \in \Omega_Y, y(1) \ge \sum_{j=2}^d y(j) \right\}; \end{equation*}

the boundary of B is

(12)\begin{equation} \partial B = \left\{y \in \Omega_Y, y(1) = \sum_{j=2 }^d y(j) \right\}. \end{equation}

The limit stopping time

(13)\begin{equation} \tau := \inf \{k \ge 0: Y_k \in \partial B\}, \end{equation}

is the first time Y hits $\partial B$. As in [Reference Sezer65], our goal is to approximate ${\mathbb{P}}_x(\tau_n \lt \tau_0)$ by the limit probability ${\mathbb{P}}_{T_n(x)}(\tau \lt \infty)$. Define the process

\begin{equation*} Y^n =T_n(X), \end{equation*}

Yn is the same process as X except in a new coordinate system whose origin is the point $n e_1 \in \partial A_n$ in the original coordinate system (i.e., it is defined as Y above except that it has an additional constraining boundary on the line $\{y \in {\mathbb Z}_+^2: y(1) = n \}$). The limit process Y and the limit probability ${\mathbb{P}}_y(\tau \lt \infty)$ are obtained by letting $n\rightarrow \infty$, see Figure 2 for an illustration in two dimensions.

Define

(14)\begin{align} R_\rho & := \bigcap_{i =1 }^d \left\{x \in {\mathbb R}^d_+: \sum_{j=1}^{i} x(j) \le \left(1 - \frac{\log \rho}{\log \rho_i}\right) \right\},\notag\\ \bar{R}_{\rho,n} &:= \bigcup_{i =1}^d \left\{x \in {\mathbb Z}^d_+: \sum_{j=1}^{i} x(j) \ge 1+ n\left(1 - \frac{\log \rho}{\log \rho_i}\right) \right\},\notag\\ A &= \left \{x \in {\mathbb R}_+^d: \sum_{j=1}^d x(j) \le 1 \right \},\\ \nonumber g:{\mathbb R}_+^d \mapsto {\mathbb R},~~~ g(x) &:= \min_{i \in \{1,2,..,d\}} \end{align}
\begin{align*} \left( 1-\sum_{j=1}^i x(j)\right) \log_{\rho}\rho_i. \end{align*}

Recall that $\rho = \max_{i=1}^d \rho_j$. Let $i^* = {\ \text{argmax}}_{i=1}^d \rho_i$; then $R_\rho \subset\left \{x \in {\mathbb R}_+^d: \sum_{j=1}^{i^*} x(j) = 0 \right\}$; in particular 1) Rρ is always a measure zero subset of A and 2) $R_\rho = \{0 \in {\mathbb R}_+^d \}$ when $i^* = d.$ The main convergence result of the paper is the following:

Theorem 2.1. For ϵ > 0 there exists $n_0 \gt 0$ such that

(15)\begin{equation} \frac{|{\mathbb{P}}_{x}(\tau_n \lt \tau_0) - {\mathbb{P}}_{T_n(x)}(\tau \lt \infty)|}{{\mathbb{P}}_{x}(\tau_n \lt \tau_0)} \le \rho^{n(1 - g(x/n)-\epsilon)}, \end{equation}

for all $n \gt n_0$ and for any $x \in \bar{R}_{\rho,n}\subset{\mathbb Z}_+^d$. In particular, for $x_n/n \rightarrow x \in A - R_\rho\subset{\mathbb R}_+^d$ the relative error decays exponentially with rate $-\log(\rho)(1-g(x)) \gt 0$, i.e.,

(16)\begin{equation} \liminf_n -\frac{1}{n} \log \left( \frac{|{\mathbb{P}}_{x_n}(\tau_n \lt \tau_0) - {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty)|}{{\mathbb{P}}_{x_n}(\tau_n \lt \tau_0)} \right) \ge -\log(\rho)(1-g(x)). \end{equation}

For a more precise characterization of the sets R and $\bar{R}_{\rho,n}$ using an order relation on $\{1,2,3,...,d\}$, see (80).

For a finite set a let $|a|$ denote its cardinality. Theorem 5.7 gives the following formula for ${\mathbb{P}}_y(\tau \lt \infty)$ under the further assumption $\mu_i \neq \mu_j$ for ij:

(17)\begin{align} {\mathbb{P}}_y(\tau \lt \infty) = \sum_{m=1}^{d} \left( \prod_{l=m+1}^{d} \frac{\mu_l -\lambda}{\mu_l - \mu_m} \right) &\rho_m^{y(1)-\sum_{j=2}^m y(j)}\\ &\times \left( \sum_{a \subset \{1,2,3,..m-1\}} (-1)^{|a|} \prod_{j=1}^{|a|} \prod_{l = a(j) + 1}^{a(j+1)} \frac{\mu_l - \lambda}{\mu_l-\mu_{a(j)}} \rho_{a(j)}^{y(l)}\right),\notag \end{align}

where a(j) denotes the j th smallest element of a, $a=\varnothing$ is always included in the last sum, and by convention, $a(|a|+1) =m.$ For d = 2 (17) reduces to

(18)\begin{equation} {\mathbb{P}}_y(\tau \lt \infty) = \frac{\mu_2 - \lambda}{\mu_2 - \mu_1} \rho_1^{y(1)} + \rho_2^{y(1)-y(2)} \left( 1 - \frac{\mu_2-\lambda}{\mu_2 -\mu_1}\rho_1^{y(2)}\right), \end{equation}

which equals the formula derived in [Reference Sezer65] for this probability for the case d = 2.

A precise statement corresponding to Figure 2 is [Reference Sezer65, Proposition 1] which states

(19)\begin{equation} \lim_n {\mathbb{P}}_{T_n(y)}(\tau_n \lt \tau_0) = {\mathbb{P}}_y(\tau \lt \infty), y \in B, \end{equation}

for any stable Jackson network in any dimension. In (19) the initial point of the process is specified and fixed in y-coordinates; in (16) it is specified in scaled x coordinates (as is done in large deviations analysis). For fixed $y \in B$, the probability ${\mathbb{P}}_{T_n(y)}(\tau_n \lt \tau_0)$ does not decay to 0 in n but converges to the nonzero probability ${\mathbb{P}}_y(\tau \lt \infty).$ In (16), where the initial position is fixed in scaled x coordinates, both of the probabilities ${\mathbb{P}}_{x_n}(\tau_n \lt \tau_0)$ and ${\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty)$ decay to 0 exponentially. The limit (16) expresses that the difference between them decays exponentially faster than ${\mathbb{P}}_{x_n}(\tau_n \lt \tau_0)$.

A function $h:\Omega_Y \mapsto {\mathbb R}$ is said to be Y-(sub/super)harmonic (or (sub/super) harmonic with respect to Y) if

(20)\begin{equation} h(y) = (\le ~/ \ge) {\mathbb E}_y[h(Y_1)], \end{equation}

for all $y \in \Omega_Y$. Theorem 2.1 reduces the approximation of ${\mathbb{P}}_x(\tau_n \lt \tau_0)$ to the computation of ${\mathbb{P}}_y(\tau \lt \infty).$ Considered as a function $y\mapsto {\mathbb{P}}_y(\tau \lt \infty)$ of y, this probability is the unique $\partial B$-determined Y-harmonic function taking the value 1 on $\partial B$ (see (122)). Section 4 introduces a method to construct Y-harmonic functions from solutions to a system of equations defined by a graph G with labeled edges (a “harmonic system,” see Definition 4.7; the graph G is a variable itself, each harmonic system has its own graph G). Harmonic systems can be defined using complex numbers resulting in complex valued Y-harmonic functions and they can be defined for constrained random walks arising from any Jackson network.Footnote 1 Since this generality comes with no additional effort, in Section 4 we will work in this generality. For $\beta \in {\mathbb C}$ and $\alpha \in {\mathbb C}^{\{2,3,4...,d\}}$ define the function $[(\beta,\alpha),\cdot]:{\mathbb Z}^d \mapsto {\mathbb C}$ as

(21)\begin{equation} [ (\beta,\alpha) , y ] = \beta^{y(1)-\sum_{j=2}^n y(j)} \prod_{j=2}^d \alpha(j)^{y(j)}, \end{equation}

since $y \mapsto \log([(\beta,\alpha),y])$ is linear in y, we call $[(\beta,\alpha),\cdot]$ log-linear.Footnote 2 Each node i of the graph G of a harmonic system has associated with it a variable $(\beta_i,\alpha_i) \in {\mathbb C}^d$ and a coefficient $c_i \in {\mathbb C}$; each node represents a constraint of the form $(\beta_i,\alpha_i) \in {\mathcal H}$ where ${\mathcal H}$ is the characteristic surface of Y (see (94)). The edges between nodes correspond to (1) conjugacy relations between the points $(\beta_i,\alpha_i)$ and (2) relations between the coefficients ci. The first main result of Section 4 is Theorem 4.8; which states that any solution to a harmonic system with graph G gives a Y-harmonic function of the form $\sum_{i \in G} c_i [(\beta_i,\alpha_i),\cdot]$. Subsection 4.1 introduces the concept of a simple extension of a constrained random walk associated with a Jackson network (defined in terms of its jump probability matrix) and the corresponding simple extension of regular graphs. The main idea here is the following: if a lower dimensional process is extended in a “simple” way (in the sense of Definition 4.9, an example is shown in Figure 3) to a higher dimensional process then any harmonic system associated with the lower dimensional process can also be extended to a harmonic system associated with the higher dimensional process. Furthermore, if the lower dimensional system has a solution, that solution can be extended to a solution of the higher dimensional system (Proposition 4.11).

Figure 3. Networks corresponding to p 1 and p 2 of (110), second is a simple extension of the first.

The limit probability ${\mathbb{P}}_y(\tau \lt \infty)$ is a Y-harmonic function. It has another important property: it is $\partial B$-determined, i.e., it is completely determined by the value it takes on the boundary $\partial B$ (namely, ${\mathbb{P}}_y(\tau \lt \infty) = {\mathbb E}[1_{\{\tau \lt \infty\}}])$. The formal definition for an arbitrary Y-harmonic function is given in subsection 4.2. We want the Y-harmonic functions we construct from harmonic systems to have this property as well, i.e., we want them to be completely determined by their values on $\partial B$. Proposition 4.13 gives conditions on the solution of a harmonic system so that the harmonic function defined from it is guaranteed to be $\partial B$-determined.

The results in Section 4 reduce the task of construction of harmonic functions to the construction of harmonic systems and their solutions, i.e., they show that if a given harmonic system has a solution, then the solution can be used to construct a harmonic function. Section 5 presents a particular class of harmonic systems for tandem networks and provides solutions for them. In this section, the dimension of the system is denoted by ${\boldsymbol d}$ and $d \le {\boldsymbol d}$ is used as a variable. The harmonic systems defined in Section 5 are based on a sequence of increasing regular graphs $G_{d,{\boldsymbol d}}$, $d=1,2,3,...,{\boldsymbol d}$; the nodes of $G_{d,{\boldsymbol d}}$ are the sets $a \cup \{d\}$, $a \subset \{1,2,3,...,d-1\}$, two nodes are connected with an edge labeled $d \ge j \ge 2$ if j lies in their intersection and the set difference between the nodes equals $\{j-1\}$ (see (128), see Figure 5 for an example). The following properties directly follow from their definitions: 1) $G_{d,{\boldsymbol d}}$ is a simple extension of $G_{d,d}$ and 2) $G_{d+1,{\boldsymbol d}}$ can be written as a disjoint union of the graphs $G_{k,{\boldsymbol d}}$, $k=1,2,...,d$; both of these properties are used in the proofs of the results that follow. Proposition 5.5 provides a solution to the harmonic systems defined by the graphs $G_{d,{\boldsymbol d}}$. The results of Section 4 imply that these solutions define $\partial B$-determined Y-harmonic functions. Theorem 5.7 then shows that the right linear combination of these functions equals $y\mapsto {\mathbb{P}}_y(\tau \lt \infty)$. In addition to the stability assumption, the results in Section 4 (in particular Theorem 5.7) assume

(22)\begin{equation} \mu_i \neq \mu_j\ \ \text{for }\ i \neq j. \end{equation}

The relative error analysis (Theorem 2.1) does not require (22); this is in contrast to previous works [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65] using the affine limit approach, all of which do make an assumption analogous to (22) for both the error analysis and the computation/approximation of ${\mathbb{P}}_y(\tau \lt \infty).$ Last paragraph of this section further comments on how assumption (22) was removed from the relative error analysis in the present work. The last result of the paper is Corollary 6.1 of Theorems 2.1 and 5.7 which generalizes (1) to

(23)\begin{equation}\lim_{n\rightarrow\infty}-\frac1n\log p_n(x_n)=-\log(\rho)g(x),\end{equation}

for any $x_n/n \rightarrow x \in A \subset{\mathbb R}_+^d.$ Since this result is based on Theorem 5.7 it also uses the assumption 22.

Results similar to Theorems 2.1 and 5.7 reported in the prior literature are as follows: [Reference Sezer65, Proposition 8], [Reference Ünlü and Sezer1, Theorem 6.1] and [Reference Başoğlu Kabran and Sezer5, Theorem 6.1] are convergence results similar to Theorem 2.1; the first concerns the tandem walk with d = 2, the second the constrained simple random walk in two dimensions (i.e., the constrained random walk with increments $\pm e_i$, $i=1,2$) and the third treats again the two-dimensional tandem walk but the dynamics are assumed to be Markov modulated (i.e., in addition to X there is an additional finite state Markov chain M that determines the jump distributions of X). These prior results state that the relative error on the left side of (16) converges to 0 exponentially at a rate depending on x; as already noted above, the precise formulation of the decay rate as in the right side of (16) is new.

The prelimit statement (15) that is specified in terms of an unscaled initial point x is also new. As already noted, Theorem 2.1 uses only the stability assumption on model parameters; all of the prior relative error results use an additional assumption corresponding to (22). On the computational side, for d = 2 [Reference Sezer65, Proposition 8] proves ${\mathbb{P}}_y(\tau \lt \infty) = f(y)$, where f is the function given in (18). For the two-dimensional constrained simple random walk and the Markov modulated tandem walk exact formulas for ${\mathbb{P}}_y(\tau \lt \infty)$ turn out to be not possible in general. Both of these works develop formulas similar to (18) that approximate ${\mathbb{P}}_y(\tau \lt \infty)$ with bounded relative error (see [Reference Ünlü and Sezer1, Propositions 7.2-7.6] and [Reference Başoğlu Kabran and Sezer5, Propositions 7.2, 7.3 and 8.3]).

As in [Reference Sezer65] we will use the following idea in our analysis of the relative error: because the dynamics of X and Y differ only on $\partial_1$, the events $\{\tau_n \lt \tau_0\}$ and $\{\tau \lt \infty\}$ mostly overlap. This is proved as follows: (1) find an event containing the difference of the events $\{\tau_n \lt \tau_0\}$ and $\{\tau \lt \infty\}$ and (2) prove that the upper bound event has a small probability. Let $\bar{\tau}_0 = \inf\left\{k \ge 0: \sum_{j=1}^d Y_k(j) = n\right\} = \inf\left\{k \ge 0: \sum_{j=1}^n T_n(Y_k(j)) = 0 \right\}$ and let $\sigma_{j,j+1}$ be the first time X hits $\partial_{j+1}$ after hitting $\partial_j$ (see (24) for the precise definition). Lemmas 3.1, 3.2 and 3.3 show that the difference between $\{\tau_n \lt \tau_0\}$ and $\{\tau \lt \infty\}$ is contained in the union of $\{\bar{\tau}_0 \lt \tau \lt \infty\}$ and $\{\sigma_{d-1,d} \lt \tau_n \lt \tau_0\}$. These results are either evident or much simpler to prove in two dimensions. When d > 2 the constraining boundaries have lower dimensional subsets. The treatment of these and a general dimension d requires nontrivial inductive and case by case and arguments (see the proofs of Lemmas 3.1-3.3). An upper bound on the probability ${\mathbb{P}}_y(\bar{\tau}_0 \lt \tau \lt \infty)$ follows from the Markov property of Y and an upper bound on ${\mathbb{P}}_y(\tau \lt \infty )$; in subsection 3.1 we construct an upper bound for this probability. In previous works treating two dimensions the upper bound follows directly from the computation/approximation of ${\mathbb{P}}_y(\tau \lt \infty)$: in [Reference Sezer65] there is a simple explicit formula for ${\mathbb{P}}_y(\tau \lt \infty)$ and in [Reference Başoğlu Kabran and Sezer5] an upper bound can be constructed in terms the Y-harmonic functions used in the computation of ${\mathbb{P}}_y(\tau \lt \infty).$ In the present setup the Y-harmonic functions that make up ${\mathbb{P}}_y(\tau \lt \infty)$ are more complex (see Theorem 5.7); therefore, in subsection 3.1 we construct simpler Y-superharmonic functions and corresponding supermartingales that imply the bound we seek on ${\mathbb{P}}_y(\tau \lt \infty)$. Subsection 3.2 derives an upper bound on the probability ${\mathbb{P}}_x(\sigma_{d-1,d} \lt \tau_n \lt \tau_0)$ (Proposition 3.9). For the proof, we construct a supermartingale corresponding to the event $\{\sigma_{d-1,d} \lt \tau_n \lt \tau_0\}$. The event happens in d stages (the process moves from stage j to j + 1 upon hitting $\partial_{j+1}$); the supermartingale is obtained by applying one of the Y-superharmonic functions of subsection 3.1 to X at each stage. Because Y has unconstrained dynamics on $\partial_1$, the process resulting from the application of these functions to X is not a supermartingale on $\partial_1$; to compensate for this we add a strictly decreasing term to the resulting process (see (62)). As in previous works [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65], we truncate time to manage this additional term (see (71)).

As already noted, a key difference between the present work and [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65] is that we no longer need to assume $\mu_i \neq \mu_j$, ij in our error analysis (Theorem 2.1). The construction of Y-harmonic functions used in the computation of $P_y(\tau \lt \infty)$ requires $\mu_i \neq \mu_j$ (assumption (22)) or analogous assumptions; this is the case in the present work (see Section 5) and in all of the previous works [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65]. Previous works [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65] use the same Y-harmonic functions in their error analysis as well, therefore their error analysis also requires an assumption analogous to (22). The graphs associated with the Y-harmonic functions constructed in Section 5 grow exponentially with dimension d; this means that the corresponding Y-harmonic functions get more complex as d increases. Their evaluation presents no difficulty unless d is large (a numerical example with d = 14 is given in Section 7); nonetheless their use in the error analysis is no longer straightforward and we have not been able to carry out an error analysis based on them. The Y-harmonic functions of Section 5 are linear combinations of log-linear functions of the form (21) where at least some of the parameters β and α take values in $\{1,\rho_1,\rho_2,...\rho_d\}.$ When we allow the β, α parameters to take values that are strictly greater than ρ, it turns out to be possible to construct Y-superharmonic functions as linear combinations of d or less log-linear functions, see Propositions 3.4 and 3.6. These superharmonic functions have a much simpler structure and we base our relative error analysis on them. Their construction requires only the stability assumption. This is why we no longer need the assumption (22) in our error analysis.

We now move on to the proof of Theorem 2.1, our main convergence result.

3. Error analysis

The goal of this section is to prove Theorem 2.1. This theorem generalizes [Reference Sezer65, Proposition 8], which treats d = 2, to an arbitrary positive dimension d > 0. The proof is based on the stopping times

(24)\begin{align} \sigma_{0,1} &:= \inf \{k \ge 0: X_k \in \partial_1\},\nonumber\\ \sigma_{j-1,j} &:= \inf \{k \gt \sigma_{j-2,j-1}: X_k \in \partial_j\}, j=2,3,...,d. \end{align}

Time τ 0 is the first time the set of all components of X equal 0; this definition and the dynamics of X imply $\tau_0 \ge \sigma_{d-1,d}.$ We will use these stopping times to show that the events $\{\tau_n \lt \tau_0\}$ and $\{\tau \lt \infty \}$ mostly overlap. Define

(25)\begin{equation} \bar{X}_{k+1} = \bar{X}_k + \pi_1(\bar{X}_k, I_{k+1}). \end{equation}

$\bar{X}$ has the same dynamics as X except on $\partial_1$ where it is not constrained. We assume $\bar{X}$ and X processes start from the same point:

(26)\begin{equation} X_0 = \bar{X}_0 = x_n. \end{equation}

For the limit analysis we will assume that the Y process has initial point $Y_0 = T_n(x_n).$ The definitions (25) and (11) of $\bar{X}$ and Y and $Y_0 = T_n(x_n)$ imply $Y = T_n(\bar{X})$ and $\bar{X} = T_n(Y).$ Define

(27)\begin{equation} \bar{\tau}_n := \inf\left \{k: \sum_{j=1}^d \bar{X}_k(j) = n \right\}. \end{equation}

We note that $Y=T_n(\bar{X})$ implies $\bar{\tau}_n = \tau$, therefore:

(28)\begin{equation} {\mathbb{P}}_{x_n}(\bar{\tau}_n \lt \infty) = {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty). \end{equation}

Define ${\boldsymbol S} : {\mathbb Z}^d \rightarrow {\mathbb Z}$ as

(29)\begin{equation} {\boldsymbol S}(x) = \sum_{j=1}^d x(j). \end{equation}

Why the stopping time $\sigma_{d-1,d}$ plays a key role in our analysis is encoded in the next lemma:

Lemma 3.1.

(30)\begin{align} \bar{X}_k(l) &\ge X_k(l), l=2,3,...,j+1, \end{align}
(31)\begin{align} \bar{X}_k(l) &= X_k(l), l= j+2,j+3,...,d, \end{align}
(32)\begin{align} {\boldsymbol S}( X_k) &= {\boldsymbol S}(\bar{X}_k) \end{align}

for $k \le \sigma_{j,j+1}$, $j \in \{0,1,2,3,...,d-1\}$;

(33)\begin{equation} {\boldsymbol S}( X_k) \ge {\boldsymbol S}(\bar{X}_k) \end{equation}

for $k \gt \sigma_{d-1,d}.$

Note that (32) holds for all j, i.e., ${\boldsymbol S}( X_k) = {\boldsymbol S}(\bar{X}_k)$, $k \le \sigma_{d-1,d}.$

Before we give a proof let us outline the main idea: our analysis concerns an event defined in terms of ${\boldsymbol S}(X)$; X and $\bar{X}$ may start to differ from each other once X hits $\partial_1$, but the tandem dynamics of X imply that for ${\boldsymbol S}(X)$ and ${\boldsymbol S}(\bar{X})$ to differ, the process must visit sequentially all of the boundaries $\partial_i$, $i=2,3,...,d$.

Proof. The processes X and $\bar{X}$ have the same dynamics except on $\partial_1$ where only X is constrained and by assumption (26) they start from the same point. Therefore, until they hit $\partial_1$ they move together, i.e.,

\begin{equation*} X_k = \bar{X}_k, \end{equation*}

for $k \le \sigma_{0,1}.$ These prove (30), (31) and (32) for j = 0. For $j \ge 1$ we will use induction. Assume (30), (31) and (32) hold for $j=j_0 \lt d-1$; let us prove that they will also hold for $j=j_0 + 1.$ We will do this by another induction on k, $ \sigma_{j_0, j_0+1} \le k \le \sigma_{j_0 + 1, j_0 + 2}.$ Note that there is a nested induction here, one induction on j another on k- we will refer to the induction on j as the outer induction and to the one on k as the inner induction. For $k = \sigma_{j_0,j_0 + 1}$ the statements hold by the outer induction hypothesis. Now assume that (30), (31) and (32), $j=j_0 +1$, hold for $k \le k_0 $ for some $\sigma_{j_0,j_0+1} \le k_0 \lt \sigma_{j_0+1,j_0+2}.$ We want to show that they must also hold for $k=k_0 + 1.$ We argue based on the possible positions of X and $\bar{X}$ at time k 0:

  1. (1) $X_{k_0} \in {\mathbb Z}_+^d - \bigcup_{j=1}^d \partial_j$: this and the inner induction hypothesis imply $\bar{X}_{k_0}(l) \gt 0$ for $l=2,3,...,d$; furthermore $\bar{X}$ is not constrained on $\partial_1$. These imply

    \begin{equation*} X_{k+1} = X_k + I_{k+1},~~ \bar{X}_{k+1} = \bar{X}_k + I_{k+1}, \end{equation*}

    i.e., both X and $\bar{X}$ change by the same increment. Therefore, all of the relations (30), (31) and (32) are preserved from time $k=k_0$ to $k=k_0 + 1.$

  2. (2) $X_{k_0}$ can also be on the boundary of ${\mathbb Z}_+^d$; recall that $\sigma_{j_0 + 1,j_0 + 2}$ is the first time X hits $\partial_{j_0 + 2}$ after time $\sigma_{j_0, j_0+1}$. Therefore, $X_{k_0} \notin \partial_{j_0 + 2}$, since $\sigma_{j_0, j_0 +1} \le k_0 \lt \sigma_{j_0 + 1, j_0 + 2}$. Then if $X_{k_0}$ is on the boundary of ${\mathbb Z}_+^d$ it must be on one of the following:

    \begin{equation*} X_{k_0} \in \partial_{\mathcal M} := \left(\bigcap_{m \in {\mathcal M}} \partial_m \right) \bigcap\left( \bigcap_{m \in {\mathcal M}^c} \partial_m^c\right), \end{equation*}

    for some ${\mathcal M} \subset\{1,2,3,...,j_0 + 1, j_0+3,...,d\}$:

    1. (a) if $I_{k_0 + 1} = e_1$, or $I_{k_0 + 1} = -e_m + e_{m+1}$ for some $m \in {\mathcal M}^c$: the increment e 1 is not constrained for X and $\bar{X}$ regardless of their position. For the case $I_{k_0 +1} = -e_m + e_{m+1}$: $X_{k_0} \in \partial_m^c$ means $X_{k_0}(m) \gt 0$. This and the inner induction hypothesis ((30) and (31)) imply $\bar{X}_{k_0}(m) \gt 0$ if m > 1; furthermore $\bar{X}$ is not constrained on $\partial_1$. These imply

      \begin{equation*} X_{k_0+1} = X_{k_0} + I_{k_0+1},~~ \bar{X}_{k_0+1} = \bar{X}_{k_0} + I_{k_0+1}. \end{equation*}

      Once again this implies that the relations (30) and (31) are preserved from time $k=k_0$ to $k= k_0 + 1.$

    2. (b) If $I_{k_0 + 1} = -e_m + e_{m+1}$ for $m \ge (j_0 +1)+2$, $m \in {\mathcal M}$: $X_{k_0} \in \partial_{\mathcal M}$ implies $X_{k_0}(m) = 0.$ By the inner induction hypothesis $\bar{X}_{k_0}(m) = X_{k_0}(m)$ for $m \ge (j_0 + 1) + 2$. Therefore, $\bar{X}_{k_0}(m) = 0$ as well. These imply that the increment $-e_m + e_{m+1}$ is constrained both for X and $\bar{X}$:

      (34)\begin{equation} X_{k_0+1} = X_{k_0},~~ \bar{X}_{k_0+1} = \bar{X}_{k_0}, \end{equation}

      and the relations (30), (31) and (32) are trivially preserved from time $k=k_0$ to $k= k_0 + 1.$

    3. (c) If $I_{k_0 + 1} = -e_m + e_{m+1}$, $2 \le m \le j_0 + 1$, $m \in {\mathcal M}$: we know by the induction hypothesis that $\bar{X}_{k_0}(m) \ge X_{k_0}(m)$. If $\bar{X}_{k_0}(m) = X_{k_0}(m)$ then the increment $-e_m + e_{m+1}$ is constrained both for X and $\bar{X}$, (34) holds and the relations (30), (31) and (32) are trivially preserved from time $k=k_0$ to $k= k_0 + 1.$ If $\bar{X}_{k_0}(m) \gt X_{k_0}(m)$ then the increment $-e_m + e_{m+1}$ is unconstrained for $\bar{X}$ while it is constrained for X:

      \begin{equation*} X_{k_0+1} = X_{k_0},~~ \bar{X}_{k_0+1} = \bar{X}_{k_0} - e_m + e_{m+1}. \end{equation*}

      The linearity of ${\boldsymbol S}$ and ${\boldsymbol S}(-e_m + e_{m+1}) = 0$ imply that (32) is preserved at time $k_0 + 1.$ All of the components $\bar{X}(l)$, $l \neq m, m+1$ remain unchanged from k 0 to $k_0 + 1$. Therefore, the relations (30) and (31) are trivially preserved for these components; in particular, this shows that (31) holds at time $k_0 + 1$ with $j=j_0 + 1$ because $m, m+1 \le (j_0 + 1) + 2.$ To complete the proof it suffices to show that (30) holds for $j=j_0 + 1$ and $k= k_0 +1$ for components l = m and $l=m+1$. For $l = m+1$,

      \begin{equation*} \bar{X}_{k_0+1}(m+1) = \bar{X}_{k_0}(m+1) + 1 \ge X_{k_0}(m+1) = X_{k_0 + 1}(m+1). \end{equation*}

      For l = m: recall that we are treating the case $\bar{X}_{k_0}(m) \gt X_{k_0}(m)$, i.e., $\bar{X}_{k_0}(m) \ge X_{k_0}(m) + 1$. Therefore:

      \begin{equation*} \bar{X}_{k_0+1}(m+1) = \bar{X}_{k_0}(m+1) - 1 \ge X_{k_0}(m+1) = X_{k_0 + 1}(m+1); \end{equation*}

      these prove that (30) holds at time $k_0 + 1$ with $j=j_0 + 1.$

    4. (d) Finally, it may happen that $1 \in {\mathcal M}$ and $I_{k_0 + 1} = -e_1 + e_2$. In this case, $X_{k_0} \in \partial_1$ and therefore the increment $I_{k_0 + 1}$ is canceled by the constraining map π for X; $\bar{X}$ is unconstrained on $\partial_1$, therefore, the increment $I_{k_0+1}$ is not constrained for $\bar{X}$. Therefore,

      \begin{equation*} \bar{X}_{k_0 + 1}(l) = \bar{X}_{k_0}(l),~~ {X}_{k_0 + 1}(l) = X_{k_0}(l), l =3,4,...,d, \end{equation*}

      and

      \begin{equation*} \bar{X}_{k_0 + 1}(2) = \bar{X}_{k_0}(2) + 1,~~ {X}_{k_0 + 1}(2) = X_{k_0}(2). \end{equation*}

      These imply that the relation (30) and (31) for $j = j_0 +1$ are preserved from time k 0 to $k_0 + 1.$ The preservation of (32) follows from the linearity of ${\boldsymbol S}$ and ${\boldsymbol S}(-e_1 + e_2) =0$ as in the last part.

This case by case analysis completes the inner induction step and hence the outer induction step.

Finally, let us consider the case $k \ge \sigma_{d-1,d}.$ Note that $X_k \in \partial_d$ for $k=\sigma_{d-1,d}$. If $I_{k+1} = -e_d$ and $\bar{X}_k \notin \partial_d$ we have:

(35)\begin{equation} X_{k+1} = X_k,~~ \bar{X}_{k+1} = \bar{X_k} - e_d, \end{equation}

an application of ${\boldsymbol S}$ to both sides of the above equations and (32) imply ${\boldsymbol S}(X_{k+1}) = {\boldsymbol S}(\bar{X}_{k+1}) +1$; thus ${\boldsymbol S}(X_{k+1}) \gt {\boldsymbol S}(\bar{X}_{k+1})$ can happen after time $\sigma_{d-1,d}$. A case by case analysis parallel to the one given above shows that (33) is preserved at all times after $\sigma_{d-1,d}.$

The previous lemma implies

Lemma 3.2. Suppose $X_0 =x \in {\mathbb Z}_+^d$, x0. The stopping times τm, $\bar{\tau}_m$, for any $m \ge 0$, and $\sigma_{d-1,d}$ satisfy:

  1. (1) $\sigma_{d-1,d} \ge \tau_m $ if and only if $\sigma_{d-1,d} \ge \bar{\tau}_m.$

  2. (2)

    (36)\begin{equation} \tau_m = \bar{\tau}_m \end{equation}

    over the event $\{\sigma_{d,d+1} \ge \tau_m\}=\{\sigma_{d,d+1} \ge \bar{\tau}_m\}.$

  3. (3)

    (37)\begin{equation} \tau_m \ge \bar{\tau}_m \end{equation}

    if $m \lt {\boldsymbol S}(x)$ and

    (38)\begin{equation} \tau_m \le \bar{\tau}_m \end{equation}

    if $ m \gt {\boldsymbol S}(x).$

Proof. By definition $\tau_m \le \sigma_{d-1,d}$ if and only if

\begin{equation*} {\boldsymbol S}(X_k) =m, \end{equation*}

for some $k \le \sigma_{d-1,d}$ and $\bar{\tau}_m \le \sigma_{d-1,d}$ if and only if

\begin{equation*} {\boldsymbol S}(\bar{X}_k) =m, \end{equation*}

for some $k \le \sigma_{d-1,d}$. By the previous lemma ${\boldsymbol S}(X_k) = {\boldsymbol S}({\bar X}_k)$ for $k \le \sigma_{d-1,d}.$ These imply the first two parts of the current lemma. Similarly,

\begin{equation*} \tau_m = \inf \{k \ge 0: {\boldsymbol S}(X_k) = m \},~~ \bar{\tau}_m = \inf \{k \ge 0: {\boldsymbol S}(\bar{X}_k) = m \}; \end{equation*}

${\boldsymbol S}(\bar{X}_k) = {\boldsymbol S}(X_k)$ for $k \le \sigma_{d-1,d}$ by Lemma 3.1, (32). Therefore, for $\tau_m \le \sigma_{d,d+1}$

\begin{equation*} \tau_m = \inf \{\sigma_{d-1,d} \ge k \ge 0: {\boldsymbol S}(X_k) = m \} = \inf \{\sigma_{d-1,d} \ge k \ge 0: {\boldsymbol S}(\bar{X}_k) = m \} = \bar{\tau}_m, \end{equation*}

i.e, (36) holds.

The relations (32) and (33) imply that

(39)\begin{equation} {\boldsymbol S}(\bar{X}_k) \le {\boldsymbol S}(X_k). \end{equation}

for all $k \ge 0.$ We will argue the case when $m \lt {\boldsymbol S}(x)$, the opposite case is argued similarly. By definition, ${\boldsymbol S}(X_{\tau_m}) = m.$ This and (39) imply ${\boldsymbol S}(\bar{X}_{\tau_m}) \le m.$ The process ${\boldsymbol S}(\bar{\boldsymbol X})$ jumps by increments of −1 (happens when $\bar X$ jumps by $-e_d$) and 1 (happens when $\bar X$ jumps by e 1). It follows that $\bar{X}$ must take all of the values $m,m+1,m+2,...,{\boldsymbol S}(x)$ in the time interval $k \in \{0,1,2,...., \tau_m\}.$ This implies $\bar \tau_m \le \tau_m.$

We now express the difference between the events $\{\tau_n \le \tau_0\}$ and $\{\tau \lt \infty \} = \{\bar{\tau}_n \lt \infty \}$ in terms of the stopping times $\sigma_{d-1,d}$ and $\bar{\tau}_0.$ For two events E 1 and E 2 let Δ denote their symmetric difference:

(40)\begin{equation} E_1 \Delta E_2 := ( E_1 - E_2) \cup (E_2 - E_1). \end{equation}

Lemma 3.3. For $X_0 = x$, $0 \lt {\boldsymbol S}(x) \lt n$

(41)\begin{equation} \{\tau_n \lt \tau_0\} \Delta \{\bar{\tau}_n \lt \infty\} \subset \{\tau_n \lt \tau_0, \tau_n \gt \sigma_{d-1,d} \} \cup \{\bar{\tau}_0 \lt \bar{\tau}_n \lt \infty\} \end{equation}

holds.

Proof. Break down $\{\tau_n \lt \tau_0\}$ and $\{\bar{\tau}_n \lt \infty\}$ into two as

(42)\begin{align} \{\tau_n \lt \tau_0 \} &= \{\tau_n \lt \tau_0, \tau_n \le \sigma_{d-1,d} \} \cup \{\tau_n \lt \tau_0, \tau_n \gt \sigma_{d-1,d} \},\nonumber\\ \{\bar{\tau}_n \lt \infty \} &= \{\bar{\tau}_n \lt \infty, \bar{\tau}_n \le \sigma_{d-1,d} \} \cup \{\bar{\tau}_n \lt \infty, \bar{\tau}_n \gt \sigma_{d-1,d} \}. \end{align}

That $\tau_n = \bar{\tau}_n$ for $\tau_n \le \sigma_{d-1,d}$ and $\sigma_{d-1,d} \le \tau_n$ if and only if $\sigma_{d-1,d} \le \bar{\tau}_n$ imply

(43)\begin{equation} \{\tau_n \lt \tau_0, \tau_n \le \sigma_{d-1,d} \} \subset \{\bar{\tau}_n \lt \infty , \bar{\tau}_n \le \sigma_{d-1,d} \}. \end{equation}

On the other hand,

(44)\begin{equation} \{\bar{\tau}_n \lt \infty , \bar{\tau}_n \le \sigma_{d-1,d} \} = \{\bar{\tau}_n \lt \infty , \bar{\tau}_n \le \sigma_{d-1,d}, \bar{\tau}_0 \lt \bar{\tau}_n \} \cup \{\bar{\tau}_n \lt \infty , \bar{\tau}_n \le \sigma_{d-1,d}, \bar{\tau}_0 \gt \bar{\tau}_n \}; \end{equation}

$\tau_0 \ge \bar{\tau}_0$ by (38) (for the case m = 0) and $\tau_n = \bar{\tau}_n$ for $\bar{\tau}_n \le \sigma_{d-1,d}$ by (36); therefore

\begin{equation*} \{\bar{\tau}_n \lt \infty , \bar{\tau}_n \le \sigma_{d-1,d}, \bar{\tau}_0 \gt \bar{\tau}_n \} \subset \{\tau_n \le \sigma_{d-1,d}, \tau_0 \gt \tau_n \}. \end{equation*}

The last line, (42), (43) and (44) imply

(45)\begin{align} &\{\tau_n \lt \tau_0\} \Delta \{\bar{\tau}_n \lt \infty\}\nonumber\\ &~~~~\subset \{\tau_n \lt \tau_0, \tau_n \gt \sigma_{d-1,d} \} \cup \{\bar{\tau}_n \lt \infty, \bar{\tau}_n \gt \sigma_{d-1,d} \} \cup \{\bar{\tau}_n \lt \infty , \bar{\tau}_n \le \sigma_{d-1,d}, \bar{\tau}_0 \le \bar{\tau}_n \}. \end{align}

Next we decompose $\{\bar{\tau}_n \lt \infty, \bar{\tau}_n \gt \sigma_{d-1,d} \}$ into two:

\begin{equation*} \{\bar{\tau}_n \lt \infty, \bar{\tau}_n \gt \sigma_{d-1,d} \} = \{\bar{\tau}_n \lt \infty, \bar{\tau}_n \gt \sigma_{d-1,d} , \bar{\tau}_0 \lt \bar{\tau}_n \} \cup \{\bar{\tau}_n \lt \infty, \bar{\tau}_n \gt \sigma_{d-1,d} , \bar{\tau}_0 \gt \bar{\tau}_n \}. \end{equation*}

The assumption $ 0 \lt {\boldsymbol S}(x) \lt n$, (37) and (38) imply $\tau_n \le \bar{\tau}_n$ and $\tau_0 \ge \bar{\tau}_0$; furthermore by the first part of Lemma 3.2 $\tau_n \gt \sigma_{d-1,d}$ if and only if $\bar{\tau}_n \gt \sigma_{d-1,d}$; these imply

\begin{equation*} \{\bar{\tau}_n \lt \infty, \bar{\tau}_n \gt \sigma_{d-1,d} , \bar{\tau}_0 \gt \bar{\tau}_n\} \subset \{\tau_n \gt \sigma_{d-1,d} , \tau_0 \gt \tau_n\}. \end{equation*}

The last two displays give:

\begin{equation*} \{\bar{\tau}_n \lt \infty, \bar{\tau}_n \gt \sigma_{d-1,d} \} \subset \{\bar{\tau}_n \lt \infty, \bar{\tau}_n \gt \sigma_{d-1,d} , \bar{\tau}_0 \lt \bar{\tau}_n \} \cup \{\tau_n \gt \sigma_{d-1,d} , \tau_0 \gt \tau_n\}. \end{equation*}

This and (45) imply (41).

By this lemma, the numerator of the relative error

\begin{equation*} \frac{|{\mathbb{P}}_{x}(\tau_n \lt \tau_0) - {\mathbb{P}}_{T_n(x)}(\tau \lt \infty)|}{{\mathbb{P}}_{x}(\tau_n \lt \tau_0)} \end{equation*}

is bounded by

\begin{align*} | {\mathbb{P}}_{x_n}(\tau_n \lt \tau_0 ) - {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty) | &\le {\mathbb{P}}_{x_n}( \{\tau_n \lt \tau_0 \} \Delta \{\bar{\tau}_n \lt \infty\}) \\ &\le {\mathbb{P}}_{x_n}( \tau_n \lt \tau_0, \tau_n \gt \sigma_{d-1,d} ) +{\mathbb{P}}_{x_n}(\bar{\tau}_0 \lt \bar{\tau}_n \lt \infty). \end{align*}

The next two subsections derive upper bounds on the last two probabilities. Subsection 3.3 following them 1) derives a lower bound on the denominator ${\mathbb{P}}_{x_n}(\tau_n \lt \tau_0)$ of the relative error 2) combines the upper and lower bounds obtained to give a proof of Theorem 3.15.

3.1. Upper bound on ${\mathbb{P}}_x(\bar{\tau}_n$ < ∞)

Recall ${\mathbb{P}}_x(\bar{\tau}_n \lt \infty) = {\mathbb{P}}_{T_n(x)}( \tau \lt \infty)$ (see (28)). The function $y\mapsto {\mathbb{P}}_{y}( \tau \lt \infty)$ is Y-harmonic. In Sections 4 and 5 we will compute this Y-harmonic function exactly and see that $y\mapsto {\mathbb{P}}_{y}( \tau \lt \infty)$ has a rather intricate structure. It turns out to be possible to derive the upper bounds we need using much simpler Y-superharmonic functions. We will first derive an upper bound on the probability ${\mathbb{P}}_{x_n}(\bar{\tau}_n \lt \infty)={\mathbb{P}}_{T_n(x)}(\tau \lt \infty)$; the bound we seek on ${\mathbb{P}}_{x_n}(\bar{\tau}_0 \lt \bar{\tau}_n \lt \infty)$ will follow from the bound on ${\mathbb{P}}_{y}(\tau \lt \infty)$ by the Markov property of Y.

The n variable plays no role here and therefore we will derive the bound and the superharmonic functions in terms of the Y process- the bound for the $\bar{X}$-process will follow by the change of variable $x = T_n(y).$

A real valued function h is said to be Y-superharmonic on a set $A \subset \Omega_Y$ if

(46)\begin{equation} {\mathbb E}_y[h(Y_1)] \le h(y), y \in A. \end{equation}

We say h is Y-superharmonic if $A= \Omega_Y$.

Define

(47)\begin{equation} h_{k,r}(y) := r^{y(1) - \sum_{j=2}^k y(j)}, k \in \{1,2,3,...,d\}, r \gt 0. \end{equation}

Note that $\log(h_{k,r})$ is linear in y, i.e., $h_{k,r}$ is log-linear. This property of the function is compatible with the dynamics of the Y process in the sense that it reduces questions about its (sub/super) harmonicity to algebraic equations/inequalities involving r. The choice of the particular linear combination $y(1) - \sum_{j=2}^k y(k)$ has to do with exit boundary $\partial B$: note that $h_{d,r}(y) =1$ for $y \in \partial B$, which is the boundary behavior that we are interested in; $h_{k,r}$ for k < d can be thought of as lower dimensional versions of $h_{d,r}.$ We will be using log-linear functions in the next sections as well when we construct classes of Y-harmonic functions.

Proposition 3.4. The function $h_{k,r}$ satisfies

(48)\begin{equation} {\mathbb E}_y[h_{k,r}(Y_1)] -h_{k,r}(y) = \begin{cases} h_{k,r}(y) \left( \lambda \left(\frac{1}{r}-1\right) + \mu_1 (r-1)\right), &\ \text{if } k=1,\\ h_{k,r}(y) \left( \lambda \left(\frac{1}{r}-1\right) + \mu_k (r-1) \mathbf{1}_{\partial_k^c}(y)\right), &\ \text{if } k \in \{2,3,...,d\}. \end{cases} \end{equation}

In particular, for $r \in (\rho,1)$, $h_{1,r}$ is Y-superharmonic and for $k \in \{2,3,4,...,d\}$ $h_{k,r}$ is Y-superharmonic on $\Omega_Y - \partial_k.$

In the proof we will use the following basic fact:

Lemma 3.5. For $r \in (\rho,1)$ and any $k \in \{1,2,3,4,...,d\}$

(49)\begin{equation} \lambda {\big{(}}\frac{1}{r}-1{\big{)}} + \mu_k (r-1) \lt 0. \end{equation}

Proof. The function $r\mapsto \lambda (\frac{1}{r}-1) + \mu_k (r-1)$ is strictly convex for $r \in (0,\infty)$ and equals 0 for $r=\lambda/\mu_k \le \rho$ and r = 1. It follows that it is strictly below 0 on the interval $(\rho,1).$

Proof of Proposition 3.4

The distribution of Y 1 and the definition of $h_{k,r}$ imply

(50)\begin{equation} {\mathbb E}_y[h_{k,r}(Y_1)] = h_{k,r}(y) \left( \lambda \frac{1}{r} + \mu_k r{\mathbf{1}}_{\partial_k^c}(y) + \mu_k {\mathbf{1}}_{\partial_k}(y) + \sum_{j \neq k} \mu_j \right); \end{equation}

subtracting $ h_{k,r}(y) = h_{k,r}(y) \left( \lambda_1 + \sum_{j=1}^d \mu_j\right)$ from the last expression gives (48). For k = 1, Y 1 is not constrained on $\partial_1$, therefore (50) in that case reduces to

\begin{equation*} {\mathbb E}_y[h_{k,r}(Y_1)] = h_{k,r}(y) \left( \lambda \frac{1}{\rho} + \mu_1 \rho + \sum_{j \neq 1} \mu_j \right). \end{equation*}

The rest of the argument remains the same for k = 1. The inequality (49) and (48) imply

\begin{align*} {\mathbb E}_y[h_{k,r}(Y_1)] -h_{k,r}(y) & \lt 0 , y \in \Omega_Y, \ \text{if } k=1,\\ {\mathbb E}_y[h_{k,r}(Y_1)] -h_{k,r}(y) & \lt 0 , y \in \Omega_Y -\partial_k\ \text{if } k\in \{2,3,...,d\}. \end{align*}

This proves the Y-superharmonicity of $h_{k,r}$ (on $\Omega_Y$ for k = 1 and on $\Omega-\partial_k$ for k > 1).

We note above that $h_{1,r}$ is Y-superharmonic. The function $h_{k,r}$, for k > 1, on the other hand is Y-superharmonic everywhere except on $\partial_k.$ For example, $h_{2,r}$ (i.e., the case k = 2) is not Y-superharmonic on $\partial_2$. We ask the question: can we linearly combine it with $h_{1,r}$, which is Y-superharmonic on $\partial_2$, so that the linear combination is Y-superharmonic everywhere? An attempt to answer this question for general k and the differences (48) lead to the following coefficients:

(51)\begin{equation} \gamma_1 := 1, \gamma_k := \frac{1}{d}\frac{\min_{j \lt k}\gamma_j(\lambda(1 -1/r) + \mu_{j}(1-r))}{\lambda(1/r - 1) }, k=2,3,...,d. \end{equation}

By (49), $\gamma_k \gt 0$ for $r \in (\rho,1)$. Now define

(52)\begin{equation} h_{2,k,r} := \sum_{j=1}^k \gamma_j h_{j,r}. \end{equation}

Proposition 3.6. For any $k=1,2,3,...d$ and $r \in (\rho,1)$ the function $h_{2,k,r}$ is

Y-superharmonic.

Proof. We assume throughout that $r \in (\rho,1).$ The proof is by induction. By definition $h_{2,k,r} = h_{k,r}$ for k = 1 and we know that $h_{1,r}$ is Y-superharmonic by the previous proposition. Now assume that $h_{2,k,r}$ is Y-superharmonic for some k < d; we will prove that $h_{2,k+1,r}$ must also be Y-superharmonic. The function $h_{k+1,r}$ is Y-superharmonic on $\Omega - \partial_{k+1}$ by the previous proposition; the function $h_{2,k,r}$ is Y-superharmonic by the induction hypothesis. These and $\gamma_{k+1} \gt 0$ imply that $h_{2,k+1,r}$ is Y-superharmonic on $\Omega - \partial_{k+1}.$ Therefore, it suffices to prove that $h_{2,k+1,r}$ is Y-superharmonic on $\partial_{k+1}.$ Choose any $y\in \partial_{k+1}$ and let

\begin{equation*} k_0 = \max\{j \le k+1: y(j) \gt 0 \} \vee 1, \end{equation*}

where, by convention the max of the empty set is $-\infty.$ By the induction hypothesis $h_{2,k_0-1,r}$ is Y-superharmonic on $\partial_k$. Therefore it suffices to prove that

\begin{equation*} h_{2,k_0,k+1,r} := \sum_{j=k_0}^{k+1} \gamma_j h_{j,r} \end{equation*}

satisfies the Y-superharmonicity condition for the chosen y. The definition of k 0 implies $y(j) = 0$ for $k_0 \lt j \le k+1.$ This and the definition of $h_{k,r}$ imply $h_{j,r}(y) = h_{k_0,r}(y)$ for all $k_0 \le j \le k+1$. Then

(53)\begin{align} {\mathbb E}_y&[h_{2,k_0,k+1,r}(Y_1)] -h_{2,k_0,k+1,r}(y) \notag\\ &= h_{k_0,r}(y) \left( \lambda \left(\frac{1}{r}-1\right) \sum_{j=k_0+1}^{k+1} \gamma_j + \gamma_{k_0} \left( \lambda\left(\frac{1}{r}-1\right) + \mu_{k_0}(r-1)\right) \right). \end{align}

The definition (51) and $j \gt k_0$ implies

\begin{equation*} \gamma_j \le \frac{1}{d} \frac{1}{\lambda\left(\frac{1}{r}-1\right) } \gamma_{k_0} \left( \lambda\left(1-\frac{1}{r}\right) +\mu_{k_0}(1-r)\right). \end{equation*}

Substituting this in (53) gives

\begin{align*} {\mathbb E}_y&[h_{2,k_0,k+1,r}(Y_1)] -h_{2,k_0,k+1,r}(y) \\ &\le h_{k_0,r}(y) \gamma_{k_0}\frac{d-(k+1)+k_0}{d} \left( \lambda\left(\frac{1}{r}-1\right) + \mu_{k_0}(r-1)\right) \lt 0, \end{align*}

which completes the induction step.

Applying the Y-harmonic function $h_{2,d,r}$ to the process Y we obtain the supermartingale $h_{2,d,r}(Y_k)$, which gives us the bound we seek on ${\mathbb{P}}_y(\tau \lt \infty)$:

Proposition 3.7. For $y \in \Omega_Y$ and $ r \in (\rho, 1)$

(54)\begin{equation} {\mathbb{P}}_y(\tau \lt \infty) \le\frac{1}{\gamma_d} h_{2,d,r}(y). \end{equation}

Proof. Let m be a positive integer. The optional sampling theorem applied to the supermartingale $k \mapsto h_{2,d,r}(Y_k)$ at the stopping time $m \wedge \tau$ gives

\begin{align*} h_{2,d,r}(y) &\ge {\mathbb E}_{y}[ h_{2,d,r}(Y_{\tau \wedge m})]\\ &= {\mathbb E}_{y}[ h_{2,d,r}(Y_{\tau}){\mathbf{1}}_{\{\tau \le m \}}] +{\mathbb E}_{y}[ h_{2,d,r}(Y_m){\mathbf{1}}_{\{\tau \gt m \}}]. \end{align*}

By definition $h_{2,d,r} \ge 0$ and $h_{2,d,r}/\gamma_d \ge h_{d,r} = 1$ on $\partial B.$ These and the previous display imply

\begin{equation*} h_{2,d,r}(y)/\gamma_d \ge {\mathbb{P}}_{y}(\tau \le m ). \end{equation*}

Letting $m \rightarrow \infty$ gives (54).

By definition

(55)\begin{equation} {\mathbb{P}}_x(\bar{\tau}_n \lt \infty) = {\mathbb{P}}_{T_n(x)}(\tau \lt \infty). \end{equation}

The bound on ${\mathbb{P}}_x(\bar{\tau}_0 \lt \bar{\tau}_n \lt \infty)$ now follows from the previous proposition and the Markov property of $\bar{X}$:

Proposition 3.8. For $r \in (\rho,1)$, $x \in {\mathbb Z} \times {\mathbb Z}_+^{d-1}$ and $ 0 \lt {\boldsymbol S}(x) \lt n$

(56)\begin{equation} {\mathbb{P}}_x(\bar{\tau}_0 \lt \bar{\tau}_n \lt \infty) \le r^n \frac{1}{\gamma_d} \sum_{j=1}^d \gamma_j. \end{equation}

Proof. By (55), (54) in x coordinates is

\begin{equation*} {\mathbb{P}}_x(\bar{\tau}_n \lt \infty) \le \frac{1}{\gamma_d}h_{2,d,r}(T_n(x))= \frac{1}{\gamma_d}\sum_{j=1}^d \gamma_j h_{d,r}(T_n(x)) = \frac{1}{\gamma_d}\sum_{j=1}^d \gamma_j r^{n- \sum_{k=1}^{j} x(k)}. \end{equation*}

This, $x \in {\mathbb Z} \times {\mathbb Z}_+^{d-1}$ and $0 \lt r \lt 1$ imply

(57)\begin{equation} {\mathbb{P}}_x(\bar{\tau}_n \lt \infty) \le r^{n}\frac{1}{\gamma_d} \sum_{j=1}^d \gamma_j, \end{equation}

for ${\boldsymbol S}(x) = 0.$ We now condition on ${\mathscr F}_{\bar{\tau}_0}$:

\begin{align*} {\mathbb{P}}_x(\bar{\tau}_0 \lt \bar{\tau}_n \lt \infty) &= {\mathbb E}_x[1_{\{\bar{\tau}_0 \lt \infty\}} 1_{\{\bar{\tau}_0 \lt \bar{\tau}_n\}} 1_{\{\bar{\tau}_n \lt \infty\}}]\\ &={\mathbb E}_x[ {\mathbb E}[ 1_{\{\bar{\tau}_0 \lt \bar{\tau}_n\}} 1_{\{\bar{\tau}_0 \lt \infty\}} 1_{\{\bar{\tau}_n \lt \infty\}} |{\mathscr F}_{\bar{\tau}_0}]]\\ &= {\mathbb E}_x[ 1_{\{\bar{\tau}_0 \lt \infty\}} 1_{\{\bar{\tau}_0 \lt \bar{\tau}_n \}} {\mathbb E}[ 1_{\{\bar{\tau}_n \lt \infty\}}| {\mathscr F}_{\bar{\tau}_0}]], \end{align*}

The strong Markov property of $\bar{X}$ implies

\begin{align*} &={\mathbb E}_{x}\left[ 1_{\{\bar{\tau}_0 \lt \infty\}} 1_{\{\bar{\tau}_0 \lt \bar{\tau}_n \}} {\mathbb{P}}_{\bar{X}_{\bar{\tau}_0}}( \bar{\tau}_n \lt \infty)\right]. \end{align*}

This, ${\boldsymbol S}(\bar X_{\bar{\tau}_0}) = 0$ and (57) imply (56).

3.2. Upper bound on ${\mathbb{P}}_x(\sigma_{d-1,d} $ < $\tau_n $ < $\tau_0)$

The goal of this subsection is to prove the following bound:

Proposition 3.9. For any ϵ > 0 there exists $n_0 \gt 0$ such that

(58)\begin{equation} {\mathbb{P}}_{x}(\sigma_{d-1,d} \lt \tau_n \lt \tau_0) \lt \rho^{n(1-\epsilon)}, \end{equation}

for all $n \gt n_0$ and $x \in A_n.$

As in the previous section and as in two dimensions treated in [Reference Sezer65] we will construct a supermartingale to upper-bound ${\mathbb{P}}_x(\sigma_{d-1,d} \lt \tau_n \lt \tau_0)$. The event $\{\sigma_{d-1,d} \lt \tau_n \lt \tau_0\}$ consists of at most d + 1 stages: the process X starts on or away from $\partial_1$, then hits $\partial_2$, then hits $\partial_3$, etc. and finally hits $\partial A_n$ after hitting $\partial_d$ without ever hitting 0. Roughly, the supermartingale will be constructed by applying one of the functions $h_{2,k,r}$ to the process X at each of these stages. The next lemma is used to adjust the definition so that the defined process remains a supermartingale as X jumps from one stage to the next.

For $k \in \{2,3,...,d\}$ define

(59)\begin{equation} \gamma_{k-1,k} := \frac{\gamma_{k-1}}{\gamma_{k-1} + \gamma_{k}}. \end{equation}

Lemma 3.10. For $k \in \{2,3,...,d\}$ and $r \in (\rho,1)$

(60)\begin{equation} \min_{y \in \partial_{k}} \frac{h_{2,k-1,r}(y)}{h_{2,k,r}(y)} \ge \gamma_{k-1,k} \end{equation}

for $y \in \partial_{k}.$

Proof. By their definition

\begin{equation*} h_{k-1,r}(y) = h_{k,r}(y) = r^{y(1)- \sum_{j=2}^{k-1} y(j)}. \end{equation*}

for $y \in \partial_{k}.$ This and the definition of $h_{2,k-1,r}, h_{2,k,r}$ imply

\begin{equation*} \frac{h_{2,k-1,r}(y)}{h_{2,k,r}(y)} = \frac{ r^{y(1)- \sum_{j=1}^{k-1} y(j)} \left( \gamma_{k-1} + \sum_{j=1}^{k -2} \gamma_j r^{\sum_{l=2}^j y(l)} \right)} {r^{y(1)- \sum_{j=1}^{k-1} y(j)} \left( \gamma_{k} +\gamma_{k-1}+ \sum_{j=1}^{k -2} \gamma_j r^{\sum_{l=2}^j y(l)} \right)} \end{equation*}

for $y \in \partial_{k}$; this and $\gamma_j, r \gt 0$ imply (60).

Define

(61)\begin{equation} \Gamma_j := \prod_{i=2}^{j} \gamma_{i-1,i}, \end{equation}

and

\begin{equation*} S_k' := \Gamma_j h_{2,j,r}(T_n(X_k)) \ \text{for }\sigma_{j-1,j} \lt k \le \sigma_{j,j+1}, j=0,1,2,3,...,d, \end{equation*}

where, by convention, $\Gamma_{0} = \Gamma_1 = 1$, $h_{2,0,r} = r^n$, $\sigma_{-1,0} = -1$ and $\sigma_{d,d+1} = \infty$; in particular, $S_k' = r^n$ for $k \le \sigma_{0,1}$ and $S_k' = \Gamma_d h_{2,d,r}(T_n(X_k))$ for $k \gt \sigma_{d-1,d}.$ The supermartingale that we will use to upper-bound the probability ${\mathbb{P}}_{x}(\sigma_{d-1,d} \lt \tau_n \lt \tau_0)$ is

(62)\begin{equation} S_k := S_k' - k\left( \lambda \left(\frac{1}{r} -1\right) \sum_{j=1}^d \gamma_j \right) r^n. \end{equation}

Proposition 3.11. The process $\{S_k, k=0,1,2,3,...\}$ is a supermartingale.

Proof. The proof is a case by case analysis. We begin by $X_k \notin \partial_1$, i.e., $X_k(1) \gt 0.$ There are two subcases two consider: $k = \sigma_{j-1,j}$ for some $j \ge 1$ and $k \neq \sigma_{j-1,j}$ for all $j \ge 1$. For $k \neq \sigma_{j-1,j}$, $S_k' = \Gamma_j h_{2,j,r}(T_n(X_k))$ $S_{k+1}' = \Gamma_j h_{2,j,r}(T_n(X_{k+1}))$ for some $j \in \{0,1,2,3,...,d\}$; the functions $h_{2,j,r}$ are Y-superharmonic by Proposition 3.6 and therefore $h_{2,j,r}(T_n(\cdot))$ are X-superharmonic on $\partial_1^c.$ It follows from these that

(63)\begin{equation} {\mathbb E}_x[S_{k+1}' | {\mathscr F}_k] \le S_k' \end{equation}

over the event $\{X_k \notin \partial_1\} \cap \{k \neq \sigma_{j-1,j}, j=1,2,3,...,d\}\}.$ If $k = \sigma_{j-1,j}$ for some $j \in \{2,3,4,...,d\}$ we have $S_k' = \Gamma_{j-1}h_{2,j-1,r}(T_n(X_k))$ and $S_{k+1}' = \Gamma_j h_{2,j,r}(T_n(X_{k+1})$; $h_{2,j,r}$ is Y-superharmonic and therefore $h_{2,j,r}(T_n(\cdot))$is X-superharmonic on $\partial_1^c$. These imply

(64)\begin{equation} \Gamma_j h_{2,j,r}(T_n(X_{k})) \ge {\mathbb E}[\Gamma_j h_{2,j,r}(T_n(X_{k+1})) |{\mathscr F}_k] = {\mathbb E}[S_{k+1}' |{\mathscr F}_k] \end{equation}

over the event $\{k = \sigma_{j-1,j}\} \cap \{X_k \notin \partial_1\}.$ That $X_{k} \in \partial_j$ for $k=\sigma_{j-1,j}$ and Lemma 3.10 imply

\begin{equation*} S_{k}'=\Gamma_{j-1} h_{2,j-1,r}(T_n(X_{k}))\ \ge \Gamma_j h_{2,j,r}(T_n(X_k)). \end{equation*}

This and (64) imply

\begin{equation*} S_k' \ge {\mathbb E}[S_{k+1}' |{\mathscr F}_k] \end{equation*}

over the event $\{k = \sigma_{j-1,j}\} \cap \{X_k \notin \partial_1.\}.$ This and (63) imply $S_k' \ge {\mathbb E}[S_{k+1}' |{\mathscr F}_k]$; subtracting $k\left( \frac{\lambda}{r} + \sum_{j=1}^d \mu_j \right) r^n$ from the left and $(k+1)\left( \frac{\lambda}{r} + \sum_{j=1}^d \mu_j \right) r^n$ from the right gives

\begin{equation*} S_k \ge {\mathbb E}[S_{k+1} |{\mathscr F}_k] \end{equation*}

over the event $\{X_k \notin \partial_1\}$.

For $x \in {\mathbb Z}_+^d$, define

(65)\begin{equation} L^*(x) := \{l \in \{2,3,4,...d\}, x(l) \neq 0 \} \end{equation}

and

\begin{equation*} d^*(x) := |L^*(x)|; \end{equation*}

$l_1(x) \lt l_2(x) \lt \cdots \lt l_{d^*}(x)$ are the members of $L^*$. Two conventions 1) $l_{d^*+1}(x) = d+1$ and 2) $l_1 = d+1$ if $d^* = 0$, i.e., if $L^*(x) = \varnothing.$ For $X_k \in \partial_1$, there are two cases to consider: 1) $k = \sigma_{j-1,j}$ for some j and 2) $k \neq \sigma_{j-1,j}$ for all j. For the latter case

(66)\begin{equation} S_k' = \Gamma_j h_{2,j,r}(T_n(X_k)), S_{k+1}' = \Gamma_j h_{2,j,r}(T_n(X_{k+1})), \end{equation}

for some j. For ease of notation, let us abbreviate $l_m(X_{k})$ to lm, and $d^*(X_{k})$ to $d^*$. Decompose $h_{2,j,r}(T_n(x))$ as

(67)\begin{equation} h_{2,j,r}(T_n(x)) = \sum_{l=1}^{(l_1-1) \wedge j}\gamma_l h_{l,r}(T_n(x)) + \sum_{m=1}^{d^*} \sum_{l=l_m}^{(l_{m+1} -1)\wedge j} \gamma_l h_{l,r}(T_n(x)), \end{equation}

where we use the conventions set above. Let us begin by considering any of the inner sums in the second sum in the last display. By the definition of lm and $l_{m+1}$, $X_k(l_m) \gt 0$ and $X_k(l) = 0$ for $l_m \lt l \lt l_{m+1}$. This implies $h_{l,r}(T_n(X_k)) = h_{l_m,r}(T_n(X_k))$ for $l_m \lt l \lt l_{m+1}$. These, $l_m \gt 1$, Proposition 3.4, the definition (51) of γl and the dynamics of X imply

\begin{align*} &{\mathbb E} \left [\sum_{l=l_m}^{(l_{m+1} -1)\wedge j} \gamma_l h_{l,r}(T_n(X_{k+1})) | {\mathscr F}_k \right] - \sum_{l=l_m}^{(l_{m+1} -1)\wedge j} \gamma_l h_{l,r}(T_n(X_k))\\ &~~~~= h_{l,r}(T_n(X_k)) \left( \lambda \left( \frac{1}{r} - 1 \right) \sum_{l=l_m+1}^{(l_{m+1} -1)\wedge j} \gamma_l + \gamma_{l_m} \left( \lambda\left(\frac{1}{r} -1\right) + \mu_{l_m}(r-1) \right)\right) \le 0. \end{align*}

Summing the last inequality over m gives

(68)\begin{equation} {\mathbb E} \left [\sum_{m=1}^{d^*}\sum_{l=l_m}^{(l_{m+1} -1)\wedge j} \gamma_l h_{l,r}(T_n(X_{k+1})) | {\mathscr F}_k \right] - \sum_{m=1}^{d^*}\sum_{l=l_m}^{(l_{m+1} -1)\wedge j} \gamma_l h_{l,r}(T_n(X_k))\\ \le 0. \end{equation}

Similarly, the definition of l 1 implies, $X_k(l) = 0$ for $l \lt l_1$, this and $h_{l,r}(T_n(x)) = r^{n-\sum_{m=1}^l x(m)}$ imply $h_{l,r}(T_n(X_k)) = r^n$ for $l \lt l_1$ over the event $\{X_k \in \partial_1\}$. These and the dynamics of X imply

\begin{align*} &{\mathbb E} \left [\sum_{l=1}^{(l_1 -1)\wedge j} \gamma_l h_{l,r}(T_n(X_{k+1})) | {\mathscr F}_k \right] - \sum_{l=1}^{(l_{1} -1)\wedge j} \gamma_l h_{l,r}(T_n(X_k))\\ &~~~~= h_{1,r}(T_n(X_k)) \left( \lambda \left( \frac{1}{r} - 1 \right) \sum_{l=1}^{(l_{1} -1)\wedge j} \gamma_l \right) = r^n \left( \lambda \left( \frac{1}{r} - 1 \right) \sum_{l=1}^{(l_{1} -1)\wedge j} \gamma_l \right) \ge 0 \end{align*}

over the event $\{X_k \in \partial_1\}.$ Putting together the last display, (68), (67) and $0 \lt \Gamma_j \le 1$ give

(69)\begin{equation} {\mathbb E}[ \Gamma_{j} h_{2,j,r}(T_n(X_{k+1}))| {\mathscr F}_k] \le \Gamma_{j} h_{2,j,r}(T_n(X_{k})) + \left( \lambda \left( \frac{1}{r} - 1 \right) \sum_{l=1}^{d} \gamma_l \right) \end{equation}

over the event $\cap_{j=1}^d \{X_k \in \partial_1, k \neq \sigma_{j-1,j}\}$; this and (66) imply

\begin{equation*} {\mathbb E}[S_{k+1}' | {\mathscr F}_k] \le S_{k}' + r^n \left( \lambda \left( \frac{1}{r} - 1 \right) \sum_{l=1}^{d} \gamma_l \right). \end{equation*}

Moving the last expression to the left of the inequality sign and subtracting $k \left( \lambda \left( \frac{1}{r} - 1 \right) \sum_{l=1}^{d} \gamma_l \right) $ from both sides give ${\mathbb E}[S_{k+1} | {\mathscr F}_k] \le S_{k}$ over the same event. It remains to show

(70)\begin{equation} {\mathbb E}[S_{k+1} | {\mathscr F}_k] \le S_{k} \end{equation}

over the event $\cup_{j=1}^d \{X_k \in \partial_1, k = \sigma_{j-1,j} \}.$ In this case

\begin{equation*} S_{k}'=\Gamma_{j-1} h_{2,j-1,r}(T_n(X_{k})), S_{k+1}' = \Gamma_j h_{2,j,r}(T_n(X_{k+1})), \end{equation*}

for some $j \in \{1,2,3,...,d\}$ and $X_k \in \partial_j.$ By Lemma 3.10

\begin{equation*} S_{k}'=\Gamma_{j-1} h_{2,j-1,r}(T_n(X_{k})) \ge \Gamma_{j} h_{2,j,r}(T_n(X_{k})); \end{equation*}

this and (69) imply (70).

The upper bound on ${\mathbb{P}}_{x}(\sigma_{d-1,d} \lt \tau_n \lt \tau_0)$ now follows from the supermartingale constructed above:

Proof of Proposition 3.9

There is no loss of generality in assuming ϵ < 1, otherwise $1-\epsilon \le 0$ and (58) holds trivially.

To use the supermartingale S to bound ${\mathbb{P}}_{x}(\sigma_{d-1,d} \lt \tau_n \lt \tau_0) \lt \rho^{n(1-\epsilon))}$ we need to truncate time by an application of the following fact (see [Reference Dupuis, Sezer and Wang23, Proposition A.1 and its proof] or [Reference Sezer59, Theorem A.1.13]): there exists $c_1 \gt 0$ and $n_0 \gt 0$ such that ${\mathbb{P}}_x(\tau_n \wedge \tau_0 \gt c_1n) \le \rho^{2n}$ for $n \gt n_0$. Although they give the same results, the truncation argument varies in [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65]; below we closely follow the one given in [Reference Başoğlu Kabran and Sezer5]. We decompose ${\mathbb{P}}_{x}(\sigma_{d-1,d} \lt \tau_n \lt \tau_0)$:

(71)\begin{equation} {\mathbb{P}}_{x}(\sigma_{d-1,d} \lt \tau_n \lt \tau_0) \le {\mathbb{P}}_{x}(\sigma_{d-1,d} \lt \tau_n \lt \tau_0 \le c_6 n) + \rho^{2n}. \end{equation}

To bound the last probability, we apply the optional sampling theorem to the supermartingale S of (62) at the bounded terminal time $\eta = c_6 n \wedge \tau_0 \wedge \tau_n$:

(72)\begin{align} r^n = S_0 &\ge {\mathbb E}_{x}[ S_{\eta}]\notag\\ &= {\mathbb E}_{x}\left[ S'_\eta - \eta \left( \lambda \left(\frac{1}{r} -1\right) \sum_{j=1}^d \gamma_j \right) r^n \right]\notag \\ &\ge {\mathbb E}_{x}\left[ S'_\eta\right] - c_6 n\left( \lambda \left(\frac{1}{r} -1\right) \sum_{j=1}^d \gamma_j \right) r^n, \end{align}

$S' \gt 0$ implies

(73)\begin{equation} {\mathbb E}_x[S'_\eta] \ge {\mathbb E}[S'_\eta {\mathbf{1}}_{\{\sigma_{d-1,d} \lt \tau_n \lt \tau_0 \le c_6n\}}]. \end{equation}

Over the event $\{\sigma_{d-1,d} \lt \tau_n \lt \tau_0 \le c_6n \}$ we have:

\begin{align*} \eta &= \tau_n, \\ S'_\eta &= \Gamma_d h_{2,d,r}(T_n(X_{\tau_n})) = \Gamma_d \sum_{j=1}^d \gamma_j h_{j,r}(T_n(X_{\tau_n})) \ge \Gamma_d \gamma_d h_{d,r}(T_n(X_{\tau_n})) = \Gamma_d \gamma_d. \end{align*}

This, (72) and (73) imply

(74)\begin{equation} \frac{1}{\gamma_d \Gamma_d} r^n \left( 1 + c_6 n\left( \lambda \left(\frac{1}{r} -1\right) \sum_{j=1}^d \gamma_j \right)\right)\ge {\mathbb{P}}_x(\sigma_{d-1,d} \lt \tau_n \lt \tau_0 \le c_6n). \end{equation}

This inequality holds for any $r \in (\rho,1)$ in particular for $r= \rho^{1-\epsilon/3}.$ Now choose n 0 so that for $n \gt n_0$ we have

\begin{equation*} \rho^{-n\epsilon/3} \ge\left( 1 + c_6 n\left( \lambda \left(\frac{1}{r} -1\right) \sum_{j=1}^d \gamma_j\right)\right)\frac{1}{\gamma_d \Gamma_d} \end{equation*}

and $\rho^{n(1-\epsilon)} \gt \rho^{n(1-2\epsilon/3)} + \rho^{2n}$. Then (74) (with $r=\rho^{1-\epsilon/3}$) and (71) imply (58).

3.3. Completion of the analysis

As the last step, we derive a lower bound on ${\mathbb{P}}_{x}(\tau_n \lt \tau_0).$ Following [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5] we will do this via subharmonic functions. Define

\begin{equation*} x\in {\mathbb Z}_+^d \mapsto g_{i,n}(x) = h_{i,\rho_i}(T_n(x))= \rho_i^{n- \sum_{j=1}^i x(j)}, k=1,2,3,...d, \end{equation*}

and

(75)\begin{equation} g_n(x) := \max_{i \in \{1,2,3,...,d\}} g_{i,n}(x). \end{equation}
Proposition 3.12.

(76)\begin{equation} g_n(x)-\rho^n \le {\mathbb{P}}_x(\tau_n \lt \tau_0). \end{equation}

Proof. That $g_{i,n}(x) = h_{i,\rho_i}(T_n(x))$ and the calculation in the proof of Proposition 3.4 give

\begin{equation*} {\mathbb E}_x[g_{i,n}(X_1)] -g_{i,n}(x) = g_{i,n}(x) \left( \lambda (\frac{1}{\rho_i}-1) + \mu_k (\rho_i-1) {\mathbf{1}}_{\partial_k^c}(x)\right). \end{equation*}

The right side of this equality is 0 for $x \in \partial_k^c$ and positive for $x \in \partial_k.$ It follows that $g_{i,n}$ is X-subharmonic on ${\mathbb Z}_+^d.$ Therefore, $k\mapsto g_{i,n}(X_k)$ is a submartingale. The stability of X implies that $\tau_0 \lt \infty$ almost surely. This, that $k\mapsto g_{i,n}(X_k)$ is a submartingale and the optional sampling theorem give

\begin{align*} g_{i,n}(x) &\le {\mathbb E}[ g_{i,n}(X_{\tau_n})1_{\{\tau_n \lt \tau_0\}}] + {\mathbb E}[g_{i,n}(0)1_{\{\tau_n \gt \tau_0\}}] \end{align*}

$g_{i,n} \le 1$ on $\partial A_n$ and $g_{i,n}(0) = \rho_i^n$ imply

\begin{align*} &\le {\mathbb{P}}_x(\tau_n \lt \tau_0) + \rho_i^n. \end{align*}

Applying $\max_{i\in \{1,2,3,...,d\}}$ to both sides gives (76).

Define the order relation $\preccurlyeq$ on the nodes $\{1,2,3,...,d\}$ as follows:

(77)\begin{equation} i \preccurlyeq j \ \text{if }\rho_i \le \rho_j\ \text{and }i \le j. \end{equation}

It follows from its definition that $\preccurlyeq$ is a partial order relation (it is reflexive, antisymmetric, transitive). Define

(78)\begin{equation} {\mathscr M} := \{i \in \{1,2,3,...,d\}: \nexists j \in \{1,2,3,..,d\} \ \text{such that } \rho_j \ge \rho_i \ \text{and } j \gt i \}. \end{equation}

The set ${\mathscr M}$ consists exactly of the maximal elements of the relation $\preccurlyeq$. d is the maximum of $\{1,2,3,...,d\}$, therefore there can be no $j \in \{1,2,3,...,d\}$ satisfying j > d, this implies that $d \in {\mathscr M}$ always holds, in particular, ${\mathscr M}$ is never empty. A similar argument implies $\rho_i \neq \rho_j$ for $i,j \in {\mathscr M}$. Let us label members of ${\mathscr M}$ by i 1, i 2,..., $i_{|{\mathscr M}|}$ so that

(79)\begin{equation} \rho_{i_1} \gt \rho_{i_2} \gt \rho_{i_3} \gt \cdots \gt \rho_{i_{|{\mathscr M}|}}; \end{equation}

$i_1 \lt i_3 \lt \cdots \lt i_{|{\mathscr M}|}$ and $i_{|{\mathscr M}|} = d$ once again follow from the definitions just given.

The point $x \in A_n$ must satisfy $g_n(x) \gt \rho^n$ for the bound (76) to be nontrivial. Note that for $g_n(x) \gt \rho^n$ we have $g_n(x)/(g_n(x) - \rho^n) \gt 0.$ In the proof of Theorem 3.15, we also need a uniform upper bound on this ratio. To identify points $ x\in A_n \subset{\mathbb Z}_+^d$ for which such an upper bound holds we need the following definitions. Let

(80)\begin{align} R_\rho &:= \bigcap_{i \in {\mathscr M} } \left\{x \in {\mathbb R}^d_+: \sum_{j=1}^{i} x(j) \le \left(1 - \frac{\log \rho}{\log \rho_i}\right) \right\},\\ R_{\rho,n} &:= \{x \in {\mathbb Z}_+^d: x/n \in R_\rho \}, \notag \\ \nonumber \bar{R}_{\rho,n} &:= \bigcup_{i \in {\mathscr M} } \left\{x \in {\mathbb Z}^d_+: \sum_{j=1}^{i} x(j) \ge 1+ n\left(1 - \frac{\log \rho}{\log \rho_i}\right) \right\}. \end{align}

$\bar{R}_{\rho,n}$ is almost the complement of $R_{\rho,n}$ (we say “almost” because of the 1 appearing in the definition of $\bar{R}_{\rho,n}$. The 1 corresponds to taking an additional step out of the boundary of $R_{\rho,n}$. This ensures a uniform upper bound on $g_n(x)/(g_n(x) - \rho^n)$, see the proof of the next proposition). The next proposition shows that $0 \lt g_n(x)/(g_n(x) - \rho^n) \lt 1/(1-\rho)$ holds for $x \in \bar{R}_{\rho,n}$. Note that if $\sum_{j=1}^{\min {\mathscr M} } x(j) \ge 1$ then $x \in \bar{R}_{\rho,n} \subset{\mathbb Z}_+^d$; in particular $\bar{R}_{\rho,n} = \{x \in{\mathbb Z}_+^d: \sum_{j=1}^{d} x(j) \ge 1\}$ if ${\mathscr M} = \{d\}$ (i.e., if $\rho_d \ge \rho_i$ for all i).

Proposition 3.13. The following hold:

(81)\begin{equation} g_n(x) = \max_{i \in {\mathscr M}} \rho_i^{n -\sum_{j=1}^i x(j)}\ \text{for } x \in {\mathbb Z}_+^d, ~~ \min_{x \in A_n} g_n(x) = \rho^n, \end{equation}
(82)\begin{equation} \{x : \in {\mathbb Z}_+^d: g_n(x) = \rho^n\} = R_{\rho,n}, \end{equation}

and

(83)\begin{equation} 0 \lt \frac{g_n(x)}{g_n(x) - \rho^n} \le \frac{1}{1-\rho}, \end{equation}

for $x \in \bar{R}_{\rho,n} \subset{\mathbb Z}_+^d$.

Proof. For ij and $x \in {\mathbb Z}_+^d$, the definitions of $\preccurlyeq$, $g_{i,n}$ and $g_{j,n}$ imply

\begin{equation*} g_{i,n}(x) \le g_{j,n}(x), \end{equation*}

if $i \preccurlyeq j$. Therefore, one can replace the index set $\{1,2,3,...,d\}$ in (75) with ${\mathscr M}$; this implies the first statement in (81).

Reversing the order of min and max gives

\begin{equation*} \min_{x \in A_n} g_n(x) = \min_{x \in A_n} \max_{i \in \{1,2,3,...,d\}}g_{i,n}(x) \ge \max_{i \in \{1,2,3,...,d\}} \min_{x \in A_n} g_{i,n}(x). \end{equation*}

By the definition of $g_{i,n}$ we have

\begin{equation*} \min_{x \in A_n} g_{i,n}(x) = \rho_i^n. \end{equation*}

The last two displays imply

\begin{equation*} \min_{x \in A_n}g_n(x) \ge \rho^n. \end{equation*}

On the other hand, $g_{n}(0) = \max_{i \in \{1,2,3,...d\}} \rho_i^n = \rho^n$; this and the last display imply the second statement in (81).

Once again, the definition of $g_{i,n}$ implies

\begin{equation*} \{x: g_{i,n}(x) \le \rho^n \} = \left\{x \in {\mathbb Z}^d_+: \sum_{j=1}^{i} x(j) \le n \left(1 - \frac{\log \rho}{\log \rho_i}\right) \right\}. \end{equation*}

This and $g_n(x) = \max_{i \in {\mathscr M}} \rho_i^{n -\sum_{j=1}^i x(j)}$ imply (82).

Finally, if $x\in {\mathbb Z}_+^d$ satisfies $\sum_{j=1}^{i} x(j) \ge 1+ n\left(1 - \frac{\log \rho}{\log \rho_i}\right)$ we have $g_{i,n}(x) \ge \rho^n/ \rho_i.$ Then

\begin{equation*} 0 \lt \frac{g_n(x)}{g_n(x) - \rho^n} \le \frac{g_{i,n}(x)}{g_{i,n}(x) - \rho^n} \le \frac{\rho^n/\rho_i}{\rho^n/\rho_i - \rho^n} \le \frac{1}{1-\rho_i} \le \frac{1}{1-\rho}, \end{equation*}

which implies (83).

Recall the convention (79); i.e., $\rho_{i_1} = \rho$ and $i_1 = \min{\mathcal M}$. Therefore, $R_{\rho} \subset \{x \in {\mathbb R}_+^d: \sum_{j=1}^{i_1} x(j) = 0 \}$; in particular Rρ has strictly lower dimension than d. If $|{\mathscr M}| =1$, i.e., if $\rho_d \gt \rho_i$ for all i < d we have $R_{\rho} =\{0\}.$

Define

\begin{equation*} g: {\mathbb R}_+^d \mapsto {\mathbb R}, g(x) = \min_{i \in \{1,2,..,d\}} (1-\sum_{j=1}^i x(j))\log_{\rho}\rho_i,~~A := \{x \in {\mathbb R}_+^d, \sum_{j=1}^d x(j) \le 1\}. \end{equation*}

The following lemma follows from the definition of g and the arguments of the previous proposition:

Lemma 3.14. $g_n(x) = \rho^{ng(x/n)}$, $\max_{x \in A} g(x) =1$, $\{x \in {\mathbb R}_+^d: g(x) = 1 \} = R_\rho.$

The upper bound on the approximation error follows from the bounds above:

Theorem 3.15 For ϵ > 0 there exists $n_0 \gt 0$ such that

(84)\begin{equation} \frac{|{\mathbb{P}}_{x}(\tau_n \lt \tau_0) - {\mathbb{P}}_{T_n(x)}(\tau \lt \infty)|}{{\mathbb{P}}_{x}(\tau_n \lt \tau_0)} \le \rho^{n(1 - g(x/n)-\epsilon)} \end{equation}

for all $n \gt n_0$ and for any $x \in \bar{R}_{\rho,n}\subset{\mathbb Z}_+^d$. In particular, for $x_n/n \rightarrow x \in A - R_\rho \subset {\mathbb R}_+^d$ the relative error decays exponentially with rate $-\log(\rho)(1-g(x)) \gt 0$, i.e.,

(85)\begin{equation} \liminf_n -\frac{1}{n} \log \left( \frac{|{\mathbb{P}}_{x_n}(\tau_n \lt \tau_0) - {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty)|}{{\mathbb{P}}_{x_n}(\tau_n \lt \tau_0)} \right) \ge -\log(\rho)(1-g(x)). \end{equation}

Proof. By Lemma 3.3

(86)\begin{equation} |{\mathbb{P}}_{x_n}(\tau_n \lt \tau_0) - {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty)| \le {\mathbb{P}}_x(\bar{\tau}_0 \lt \bar{\tau}_n \lt \infty) + {\mathbb{P}}_{x}(\sigma_{d-1,d} \lt \tau_n \lt \tau_0) \end{equation}

By Proposition 3.9 we can choose n 0 large enough so that

\begin{equation*} {\mathbb{P}}_{x}(\sigma_{d-1,d} \lt \tau_n \lt \tau_0) \lt \rho^{n(1-\epsilon/2)} \end{equation*}

for $n \gt n_0.$ This, (86) and Proposition 3.8 (with $r=\rho^{1-\epsilon/2}$) give

\begin{equation*} |{\mathbb{P}}_{x_n}(\tau_n \lt \tau_0) - {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty)| \le \rho^{n(1-\epsilon/2)} + \rho^{n(1-\epsilon/2)} \frac{1}{\gamma_d} \sum_{j=1}^d \gamma_j. \end{equation*}

for $n \gt n_0.$ This and the lower bound on ${\mathbb{P}}_x(\tau_n \lt \tau)$ given in Proposition 3.12 imply

\begin{align*} \frac{|{\mathbb{P}}_{x}(\tau_n \lt \tau_0) - {\mathbb{P}}_{T_n(x)}(\tau \lt \infty)|}{{\mathbb{P}}_{x}(\tau_n \lt \tau_0)} &\le\frac{1}{g_n(x)-\rho^n} \left( \rho^{n(1-\epsilon/2)} + \rho^{n(1-\epsilon/2)} \frac{1}{\gamma_d} \sum_{j=1}^d \gamma_j\right) \end{align*}

and by Proposition 3.13 and Lemma 3.14

\begin{align*} &\le \frac{1}{1-\rho}\rho^{-n g(x/n)} \left( \rho^{n(1-\epsilon/2)} + \rho^n \frac{1}{\gamma_d} \sum_{j=1}^d \gamma_j\right)\\ &=\frac{1}{1-\rho}\rho^{-n g(x/n)+n(1-\epsilon)} \rho^{n\epsilon/2} \left( 1+ \frac{1}{\gamma_d} \sum_{j=1}^d \gamma_j\right). \end{align*}

Finally, increase n 0 if necessary so that $ \frac{1}{1-\rho} \rho^{n\epsilon/2}\left(1 + \frac{1}{\gamma_d}\sum_{j=1}^d \gamma_j\right) \lt 1$ for $n \gt n_0$; then the last display and this choice of n 0 imply (84); (85) follows from (84), the continuity of g and from the fact that ϵ > 0 can be chosen arbitrarily small.

3.4. Example in two dimensions

One of the key steps in the above analysis is the the upper bound on ${\mathbb{P}}_x(\sigma_{d-1,d} \lt \tau_n \lt \tau_0)$ derived in Proposition 3.9. Let us go over this argument in the case d = 2. We will also comment briefly on the three dimensional case. The derivation is based on a supermartingale S constructed from the functions $h_{j,r}$ and $h_{2,j,r}$, $j=0,1,2,...,d$, whose definitions are given in (52). For d = 2 these functions are:

\begin{align*} h_{1,r}(y) &= r^{y(1)}, ~~ h_{2,r} = r^{y(1)-y(2)},\\ h_{2,1,r}(y) &= h_{1,r}(y) = r^{y(1)}, \\ h_{2,2,r}(y) &= h_{1,r}(y) + \gamma_2 h_{2,r}(y) = r^{y(1)} + \gamma_2 r^{y(1) - y(2)}, \end{align*}

where (by definition (51))

\begin{equation*} \gamma_2 = \frac{1}{2} \frac{\lambda(1-1/r) + \mu_1(1-r)}{\lambda(1/r -1)}. \end{equation*}

The process S is now defined as:

\begin{align*} S_k &= S_k' - k (\lambda (1/r -1)(1+ \gamma_2))r^n,\\ S_k' &= \begin{cases} r^n, &k \le \sigma_1,\\ h_{2,1,r}(T_n(X_k)) = h_{1,r}(T_n(X_k)) = r^{n-X_k(1)}, & \sigma_1 \lt k \le \sigma_2,\\ h_{2,2,r}(T_n(X_k)) = \Gamma_2 \left( r^{n-X_k(1)} + \gamma_2 r^{n -(X_k(1) + X_k(2))}\right), & k \gt \sigma_2, \end{cases} \end{align*}

where, by (61), $\Gamma_2 = \gamma_{1,2} = \gamma_1/(\gamma_1 + \gamma_2) = 1/(1+\gamma_2).$ Note that S is defined in terms of X. Proving S is a supermartingale requires the treatment of the following cases: $k \lt \sigma_1$, $ \sigma_1 \lt k \lt \sigma_2$, $k \gt \sigma_2$, $k=\sigma_1$ and $k=\sigma_2$. The first case is trivial since for $k \lt \sigma_1$, S is strictly and deterministically decreasing. For each of $ \sigma_1 \lt k \lt \sigma_2$ and $ k \gt \sigma_2$ we have the further three cases $X_k(1),X_k(2) \gt 0$, $X_k(1) \in \partial_2 $ and $X_k(2)\in \partial_1.$ The first two subcases follow from the Y-superharmonicity of $h_{2,1,r}$ and $h_{2,2,r}$ established in Proposition 3.6. Let us also check this directly here using the fact that we are dealing with d = 2. We begin with $h_{2,1,r}$: by definition $h_{2,1,r}(y) = h_{1,r}(y) = r^{y(1)}.$ From the dynamics of Y 1:

\begin{equation*} {\mathbb E}_y[r^{Y_1}] = \lambda r^{y(1)-1} + \mu_1 r^{y(1) + 1} = r^{y(1)} (\lambda/r + \mu_1 r). \end{equation*}

Then

\begin{equation*} {\mathbb E}_y[r^{Y_1}] - r^{y(1)} = r^{y(1)} ( \lambda (1/r -1 ) + \mu_1(r-1)); \end{equation*}

by Lemma 3.5 this last expression is strictly negative, this proves that $h_{1,r}$ is Y-superharmonic (note that the argument does not depend on whether $y(2)=0$ or $y(2) \gt 0$). The argument for $h_{2,2,r}(y)= h_{1,r}(y) + \gamma_2 h_{2,r}(y) = r^{y(1)} + \gamma_2 r^{y(1)-y(2)}$ looks at two cases $y(2) \gt 0$ and $y(2) =0$ separately. In the first case, proceeding as above, we see:

\begin{equation*} {\mathbb E}_y[h_{2,r}(Y_1)] - h_{2,r}(y) = h_{2,r}(y) ( \lambda (1/r -1 ) + \mu_2(r-1)) \lt 0 \end{equation*}

where the last inequality again follows from Lemma 3.5. The last two displays and $\gamma_2 \gt 0$ establish that $h_{2,2,r} = h_{1,r} + \gamma_2 h_{2,r}$ is Y- superharmonic for $y(2) \gt 0.$ Note that for $y(2) \gt 0$ we are able to argue Y-harmonicity for $h_{1,r}$ and $h_{2,r}$ separately. For $y(2) =0$, $h_{2,r}$ is in fact Y-subharmonic. Therefore, the linear combination $h_{2,2,r} = h_{1,r} + \gamma_2 h_{2,r}$ must be considered as a whole. A calculation similar to the above for $y(2) = 0$ gives:

\begin{align*} {\mathbb E}_y[h_{2,2,r}(Y_1)] - h_{2,2,r}(y) &= h_{1,r}(y) (\lambda (1/r -1 ) + \mu_1(r-1)) + \gamma_2 h_{2,r}(y)( \lambda(1/r -1)). \end{align*}

Note that for $y(2)=0$ we have $h_{1,r}(y) = h_{2,r}(y) = r^{y(1)}$; this and the definition of γ 2 above yield

\begin{align*} &= r^{y(1)} ( \lambda(1/r -1) + \mu_1(r-1) + \frac{1}{2}( \lambda(1-1/r) + \mu_1(1-r)))\\ &= \frac{1}{2} r^{y(1)} ( \lambda(1/r -1) + \mu_1(r-1)) \end{align*}

which is again strictly negative by Lemma 3.5. These arguments establish that $h_{2,1,r}$ and $h_{2,2,r}$ are Y-superharmonic which in turn establishes the supermartingale property of S for the cases $\sigma_1 \lt k \lt \sigma_2$ and $k \gt \sigma_2$ and $X_k \notin \partial_1.$

Now let’s consider the case $\sigma_1 \lt k \lt \sigma_2$ and $X_k \in \partial_1$, for which

(87)\begin{align} S_k &= h_{1,r}(T_n(X_k)) - k \lambda (1/r-1)(1+\gamma_2)r^n \notag \\ & = r^{n - X_k(1)} - k \lambda (1/r-1)(1+\gamma_2) r^n. \end{align}

From the dynamics of X on $\partial_1$:

(88)\begin{equation} {\mathbb E}[S_{k+1} | {\mathscr F}_k] - S_k = r^n\lambda(1/r -1) - \lambda(1/r -1)(1 + \gamma_2) \lt 0, \end{equation}

where we used $X_k(1) = 0$ for $X_k \in \partial_1.$ Note that $x \mapsto h_{1,r}(T_n(x))$ is X-subharmonic on $\partial_1$ which leads to the first and positive term on the right side of the equal sign in (88). The term $- k \lambda (1/r-1)(1+\gamma_2)r^n $ in the definition of Sk compensates for this, leading to the second and negative term on the right side of the equal sign in (88) making the whole difference negative and thus ensuring the supermartingale property of S for this case.

For $\sigma_2 \lt k \le \tau_n $ and $X_k \in \partial_1$ we have

\begin{align*} S_k &= \Gamma_2(h_{1,r}(T_n(X_k)) + \gamma_2 h_{2,r}(T_n(X_k))) - k \lambda (1/r-1)(1+\gamma_2)r^n \\ &= \Gamma_2 \left( r^{n - X_k(1)} + \gamma_2 r^{n-(X_k(1)+X_k(2))}\right) - k \lambda (1/r-1)(1+\gamma_2) r^n\\ &= \Gamma_2 r^{n - X_k(1)} - k \lambda (1/r-1)(1+\gamma_2) + \Gamma_2 \gamma_2 r^{n-(X_k(1)+X_k(2))}. \end{align*}

Compared to (87) there are two new terms here: Γ2 coefficient in front of $r^{n-X_k(1)}$ and the last term. By its definition $0 \lt \Gamma_2 \lt 1$, therefore, the supermartingale property for the first two terms can proceed as in (88). Since we are working over an event that takes place before τ 0 (first hitting time to $(0,0)$) we have $X_k \neq (0,0)$ for $ \sigma_2 \lt k \le \tau_n$. Therefore, we only have to deal with $ x \in \partial_1 - \{(0,0)\} = \{x \in {\mathbb Z}_+^2: x(1)=0, x(2) \gt 0\}.$ Over this last set $ x \mapsto r^{n-(x(1) + x(2))}$ is X-superharmonic. This and $\gamma_2, \Gamma_2 \gt 0$ establish the supermartingale property for the $ \Gamma_2 \gamma_2 r^{n-(X_k(1)+X_k(2))}$ term.

Note how working in d = 2 leads to the simple case by case argument above. The second case is particularly simple: the function $x \mapsto r^{n-(x(1)+x(2))}$ is X-superharmonic on $\partial_1$ (except at $x=(0,0)$ which can be omitted because we are working over a time interval before X hits the origin). In higher dimensions, even when dimension is fixed, this argument does not work. For example, in d = 3, for $x =(0,0,1) \in \partial_1$, the term $r^{n-(X_k(1)+X_k(2))}$ cannot be handled by itself, it must be considered in a linear combination with the $-k \left(\lambda((1/r)-1)\sum_{j=1}^3 \gamma_j\right) r^n$ term; or for $x=(0,1,0) \in \partial_1$, the term $r^{n-(X_k(1)+X_k(2)+X_k(3))}$ must be considered with the $r^{n-(X_k(1)+X_k(2))}$ term. In higher dimensions, points on different lower dimensional boundary points require different linear combinations. The argument in the proof of Proposition 3.11 handles all of the possible cases in arbitrary dimension by partitioning $x \in \partial_1$ according to its nonzero elements (see (65) and (67)). Similar issues arise in proving that $h_{2,j,r}$ is Y-superharmonic. In the proof of Proposition 3.6 these issues are handled systematically in arbitrary dimension with an inductive argument.

Let us now consider the case $k=\sigma_2$, for which we have

\begin{align*} S_k &= h_{1,r}(T_n(X_k))- k ( \lambda(1/r -1) (1+\gamma_2))r^n\\ S_{k+1} &=\Gamma_2 h_{2,2,r}(T_n(X_{k+1})) - (k+1)(\lambda(1/r-1) (1+\gamma_2))r^n. \end{align*}

Note that this is one of the transitions that S goes through: S is defined by the function $h_{2,1,r} = h_{1,r}$ for $k=\sigma_2$ and by the function $h_{2,2,r}$ for k + 1. The term Γ2 is chosen so that the supermartingale property is preserved during this transition. The details are as follows: for $k=\sigma_2$ we have:

(89)\begin{align} {\mathbb E}[S_{k+1} | {\mathscr F}_k] - S_k &= {\mathbb E}[ \Gamma_2 h_{2,2,r}(T_n(X_{k+1})| {\mathscr F}_k] - \Gamma_2 h_{2,2,r}(T_n(X_k)) -(\lambda(1/r-1) (1+\gamma_2))r^n \notag \\ &+ \Gamma_2 h_{2,2,r}(T_n(X_k)) - h_{1,r}(T_n(X_k)). \end{align}

The first line is negative by the Y-superharmonicity of $h_{2,2,r}$ and the positivity of $\lambda(1/r-1) (1+\gamma_2)$. For the second line we note $k=\sigma_2$ implies $X_k \in \partial_2$, i.e., $X_k(2) = 0.$ Then (recalling $\Gamma_2 = 1/(1+\gamma_2)$)

\begin{align*} \Gamma_2 h_{2,2,r}(T_n(X_k)) &= \Gamma_2 (1+\gamma_2)r^{n-X_k(1)} = r^{n-X_k(1)}\\ h_{1,r}(T_n(X_k)) &= r^{n-X_k(1)} \end{align*}

which implies that the difference in (89) is zero. This establishes the supermartingale property of S for this case. The case $k=\sigma_1$ is handled similarly.

4. Y-harmonic functions from harmonic systems

We now build a framework for the construction of Y-harmonic functions. In this section, we assume the underlying network to be an arbitrary Jackson network (not necessarily tandem) and we also allow ${\mathbb C}$ valued harmonic functions. This generality does not lead to any significant complications in the definitions and arguments of the present section compared to the tandem case with only ${\mathbb R}$ valued harmonic functions and it can be useful in the generalization of our results to general Jackson networks. The main element of the framework is the reduction of the construction of Y-harmonic functions to the solution of certain equations represented by graphs with labeled edges, which we call harmonic systems (see Definitions 4.6 and 4.7 below). In the next section, we will define particular harmonic systems for the tandem case and provide solutions to these.

In this section, we allow ${\mathcal V}$ to be the set of possible increments of any constrained random walk arising from a Jackson network, i.e.:

\begin{equation*} {\mathcal V} = \{- e_i + e_j, i,j \in \{0,1,2,...,d\}, i \neq j\}, \end{equation*}

where $e_0 = 0 \in {\mathbb Z}^d$; the unconstrained increments of X takes values in ${\mathcal V}$ with probabilities ${\mathbb{P}}(I_k = -e_i +e_j) = p(i,j)$ where $p \in {\mathbb R}_+^{(d+1) \times (d+1)}$, $p(i,i) =0$, $i \in \{0,1,2,3,...,d\}$ and $\sum_{i,j=0}^d p(i,j) = 1.$ With this update to the set of possible increments the definition of X remains unchanged. The increment $-e_i +e_j$ represents a customer leaving node i and joining node j where node 0 represents outside of the system. For a general Jackson network the total service rates are defined as

\begin{equation*} \mu_i = \sum_{j=0}^d p(i,j), i \in \{1,2,3,...,d\}. \end{equation*}

The Y process is defined as in (11) on $\partial_1$ with possible increments

\begin{align*} {\mathcal V}_Y & = \{{\mathcal I}_1 v, v \in {\mathcal V}\}\\ & = \{v_{1,j} := e_1 +e_j, v_{i,1} := -e_i -e_1, v_{i,j} := -e_i + e_j, i,j \in \{0,1,2,3...,d\},i \neq j\} \end{align*}
(90)\begin{equation} Y_{k+1} = Y_k + \pi_1(Y_k, J_k). \end{equation}

For $\alpha \in {\mathbb C}^{d-1}$ we will index the components of the vector α with the set $\{2,3,4,...,d\}$, i.e., $\alpha = (\alpha(2),\alpha(3),...,\alpha(d))$ (so, more precisely, $\alpha \in {\mathbb C}^{\{2,3,4,...d\}}$). Y-harmonic functions resulting from solutions of harmonic systems will be linear combinations of functions of the form

(91)\begin{align} y \in {\mathbb Z}^d &\mapsto [(\beta,\alpha),y],\\ [(\beta,\alpha),y] &:= \beta^{y(1)- \sum_{j=2}^d y(j) } \prod_{j=2}^d \alpha(j)^{y(j)}.\notag \end{align}

$y\mapsto [(\beta,\alpha),y]$ is log-linear in y, i.e., $y\mapsto \log([(\beta,\alpha),y])$ is linear in y. This also means that $y\mapsto [(\beta,\alpha),y]$ itself is multiplicative, i.e., $ [(\beta,\alpha),y+v] = [(\beta,\alpha),y] [(\beta,\alpha),v]$.

For $a \subset \{2,3,...,d\}$ and $a^c = \{0,1,2,3,...,d\} - a$, define the characteristic polynomial

(92)\begin{equation} {\mathbf p}_a(\beta,\alpha) := \left( \sum_{i \in a^c, j=0 }^d p(i,j) [(\beta,\alpha) ,v_{i,j} ] + \sum_{i \in a} \mu_i\right) \end{equation}

the characteristic equation

(93)\begin{equation} {\mathbf p}_a(\beta,\alpha) = 1, \end{equation}

and the characteristic surface

(94)\begin{equation} {\mathcal H}_a := \{(\beta,\alpha) \in {\mathbb C}^d: {\mathbf p}_a(\beta,\alpha) = 0\} \end{equation}

of the boundary

\begin{equation*} \partial_a := \bigcap_{i \in a } \partial_i, a \subset \{2,3,4,...,d\}. \end{equation*}

We will write ${\mathbf p}$ instead of ${\mathbf p}_\varnothing$. ${\mathbf p}_a$ is not a polynomial but a rational function; to make it a polynomial one must multiply it by $\beta \prod_{j=2}^d \alpha(j)$; nonetheless, to keep our language simple we will refer to the rational (92) as the “characteristic polynomial.”

Conditioning Y on its first step gives

Lemma 4.1. Suppose $(\beta,\alpha) \in {\mathcal H}$. Then $[ (\beta,\alpha),\cdot]$ is Y-harmonic on $\Omega_Y - \bigcup_{j=2}^d \partial_j$.

Proof. For $y \in \Omega_Y-\bigcup_{j=2}^d \partial_j$ we have

\begin{equation*} {\mathbb E}_y[ [ (\beta,\alpha), Y_1]] - [(\beta,\alpha),y] = [(\beta,\alpha),y] ( {\mathbf p}(\beta,\alpha) -1 ) = 0, \end{equation*}

where the last equality follows from $(\beta,\alpha) \in {\mathcal H}.$

Define the operator Da acting on functions on ${\mathbb Z}^d$ and giving functions on $\partial_a$:

(95)\begin{align} &D_aV = g,~~~~ V: {\mathbb Z}^d \rightarrow {\mathbb C},\\ & g( y)\notag := \left(\sum_{i \in a} \mu_i V(y) + \sum_{i\in a^c, j=0}^d p(i,j) V( y + v_{i,j}) \right) - V(y);\notag \end{align}

Lemma 4.2. $D_a V = 0$ if and only if V is Y-harmonic on $\partial_a$.

The proof follows from the definitions. Define

(96)\begin{equation} C(i,\beta,\alpha) := \mu_i -\sum_{j=0}^d p(i,j) [ (\beta,\alpha),v_{i,j} ]. \end{equation}

Lemma 4.3. For $y \in \partial_i$ and $(\beta,\alpha) \in {\mathcal H}$:

(97)\begin{equation} D_i( [(\beta,\alpha),\cdot])(y) =C(i,\beta,\alpha) [(\beta,\alpha),y]. \end{equation}
Proof.

\begin{align*} D_i( [(\beta,\alpha),\cdot])(y) &= \left( \mu_i +\sum_{i'\neq i, j=0}^d p(i',j) [ (\beta,\alpha),v_{i',j} ] -1\right) [ (\beta,\alpha), y]\\ &= \left( \mu_i -\sum_{j=0}^d p(i,j) [ (\beta,\alpha),v_{i,j} ] \right) [ (\beta,\alpha), y], \end{align*}

where we used $1={\mathbf p}(\beta,\alpha)$ and the multiplicative property of $[(\beta,\alpha),\cdot].$

Lemma 4.4. Suppose $(\beta,\alpha) \in {\mathcal H} \cap {\mathcal H}_i$. Then $[(\beta,\alpha),\cdot]$ is Y-harmonic on $\Omega_Y-\bigcup_{j \in \{2,3,...,d\} - \{i\}} \partial_j$.

Proof. $(\beta,\alpha) \in {\mathcal H}$ and Lemma 4.1 imply that $[(\beta,\alpha),\cdot]$ is Y-harmonic on $\Omega_Y-\bigcup_{j=2}^d \partial_j.$ That ${\mathbf p}(\beta,\alpha) = 1$ and ${\mathbf p}_i(\beta,\alpha) =1$ imply

\begin{equation*} C(i,\beta,\alpha) = {\mathbf p}_i(\beta,\alpha)-{\mathbf p}(\beta,\alpha) =0. \end{equation*}

This and the last two lemmas imply that $[(\beta,\alpha),\cdot ]$ is Y-harmonic on $\partial_i.$

For $\alpha \in {\mathbb C}^{\{2,3,...,d\}}$ and $j \in \{2,3,4,...,d\}$ define $\alpha\{j\} \in {\mathbb C}^{\{2,3,...d\}}$ as follows:

(98)\begin{equation} \alpha\{j\}(i) = \begin{cases} 1,&\ \text{if } i = j\\ \alpha(i),&\ \text{otherwise.} \end{cases} \end{equation}

For example, for d = 4, j = 4 and $\alpha = (0.2,0.3,0.4)$, $\alpha\{4\} = (0.2,0.3,1).$

For $i \in \{2,3,4,...,d\}$, multiplying both sides of the characteristic equation ${\bf p}(\beta,\alpha) =0$ by $\alpha(i)$ gives a second order polynomial equation in $\alpha(i)$: denote the roots by r 1 and r 2. From the coefficients of the second order polynomial we read

(99)\begin{equation} r_1 r_2 = \frac{ \sum_{j=0}^d p(i,j) [ (\beta, \alpha\{i\}), v_{i,j}]} {\sum_{j=0}^d p(j,i) [ (\beta, \alpha\{i\}), v_{j,i}]}. \end{equation}

From these two roots, we get two points $(\beta,\alpha_1), (\beta,\alpha_2)$ on ${\mathcal H}$ whose components are

\begin{equation*} \alpha_1(i) =r_1, \alpha_2(i)=r_2, \end{equation*}

and

(100)\begin{equation} \alpha_1(j) =\alpha_2(j) = \alpha(j), j\neq i. \end{equation}

By (99)

(101)\begin{equation} \alpha_1(i)\alpha_2(i) = \frac{ \sum_{j=0}^d p(i,j) [ (\beta, \alpha\{i\}), v_{i,j}]} {\sum_{j=0}^d p(j,i) [ (\beta, \alpha\{i\}), v_{j,i}]}. \end{equation}

If $\alpha_1(i)\neq \alpha_2(i)$ we call $(\beta,\alpha_1) \neq (\beta,\alpha_2) \in {\mathcal H}$ i-conjugate. Note that $\alpha_1\{i\} = \alpha_2\{i\} = \alpha\{i\}$; therefore (101) can also be written as

(102)\begin{equation} \alpha_1(i)\alpha_2(i) = \frac{ \sum_{j=0}^d p(i,j) [ (\beta, \alpha_1\{i\}), v_{i,j}]} {\sum_{j=0}^d p(j,i) [ (\beta, \alpha_1\{i\}), v_{j,i}]}. = \frac{ \sum_{j=0}^d p(i,j) [ (\beta, \alpha_2\{i\}), v_{i,j}]} {\sum_{j=0}^d p(j,i) [ (\beta, \alpha_2\{i\}), v_{j,i}]}. \end{equation}

Next proposition generalizes [Reference Sezer65, Proposition 4] to the current setup.

Proposition 4.5. Suppose that $(\beta,\alpha_1)$ and $(\beta,\alpha_2)$ are i-conjugate and $C(i,\beta,\alpha_j)$, $j=1,2$ are well defined. Then

\begin{equation*} h_\beta := C(i,\beta,\alpha_2) [(\beta,\alpha_1),\cdot] - C(i,\beta,\alpha_1) [(\beta,\alpha_2),\cdot] \end{equation*}

is Y-harmonic on $\partial_i$.

Proof. The definition (96) of C, (97) and linearity of Di imply

\begin{align*} D_i(h_\beta) &= C(i,\beta,\alpha_2) C(i,\beta,\alpha_1) [(\beta,\alpha_1),\cdot] - C(i,\beta,\alpha_2) C(i,\beta,\alpha_1) [(\beta,\alpha_2),\cdot] \end{align*}

(100) implies $[(\beta,\alpha_1),z] = [(\beta,\alpha_2),z]$ for $z \in \partial_i$ and therefore the last line reduces to

\begin{align*} &=0. \end{align*}

Lemma 4.2 now implies that hβ is Y-harmonic on $\partial_i$.

The class of Y-harmonic functions we identify in this section is based on graphs with labeled edges; let us now give a precise definition of these. We denote any graph by its adjacency matrix G; the structure of G is as follows. Let VG, a finite set, denote the set of vertices of G; let L denote the set of labels. For two vertices ij, $G(i,j) = 0$ if they are disconnected, and $G(i,j) = l$ if an edge with label $l \in L$ connects them; such an edge will be called an l-edge. As usual, an edge from a vertex to itself is called a loop. For a vertex $j \in V_G$, $G(j,j)$ is the set of the labels of the loops on j. Thus $G(j,j)\subset L$ is set valued.

In graph theory a graph is said to be k-regular if all of its vertices have the same degree (number of edges) k [Reference Diestel20, page 5]. We generalize this definition as follows:

Definition 4.6. Let G and L be as above. If each vertex $j \in V_G$ has a unique l-edge (perhaps an l-loop) for all $l \in L$ we will call G L-regular.

Definition 4.7. A Y-harmonic system consists of a $\{2,3,4,...,d\}$-regular graph G, the variables $ (\beta,\alpha_j) \in {\mathbb C}^d$, ${\boldsymbol c}_j \in {\mathbb C}$, $j \in V_G$ and these equations/constraints:

  1. (1) $(\beta,\alpha_j) \in {\mathcal H}, {\boldsymbol c}_j \in {\mathbb C}-\{0\}, j \in V_G$,

  2. (2) $\alpha_{i} \neq \alpha_{j}$, if $i \neq j, i,j \in V_G$,

  3. (3) $\alpha_{i}, \alpha_{j}$ are $G(i,j)$-conjugate if $ G(i,j) \neq 0$, $i \neq j, i,j \in V_G$,

  4. (4)

    (103)\begin{equation} {\boldsymbol c}_{i}/{\boldsymbol c}_{j} = -\frac{C(G(i,j),\beta,\alpha_j)}{C(G(i,j),\beta,\alpha_i)}, \ \text{if } G(i,j) \neq 0, \end{equation}
  5. (5) $(\beta,\alpha_j) \in {\mathcal H}_l \ \text{for all } l \in G(j,j), j \in V_G.$

Theorem 4.8 Suppose that a Y-harmonic system with graph G has a solution $({\boldsymbol c}_j, (\beta,\alpha_j), j \in V_G)$. Then

(104)\begin{equation} h_G := \sum_{j \in V_G} {\boldsymbol c}_j [(\beta,\alpha_j),\cdot] \end{equation}

is Y-harmonic.

In the proof the following decomposition is useful: for $y \in \partial_a$ and $(\beta,\alpha) \in {\mathcal H}$:

(105)\begin{align} D_a( [(\beta,\alpha,\cdot)])(y) &= \left( \sum_{i \in a} \mu_i +\sum_{i \in a^c, j} p(i,j) [ (\beta,\alpha),v_{i,j} ] -1\right) [ (\beta,\alpha), y]\notag \\ &= \left( \sum_{i \in a} \mu_i -\sum_{i \in a, j=0}^d p(i,j) [ (\beta,\alpha),v_{i,j} ] \right) [ (\beta,\alpha), y]\notag \\ &= \sum_{i \in a} \left( \mu_i - \sum_{j=0}^d p(i,j) [ (\beta,\alpha),v_{i,j} ] \right) [ (\beta,\alpha), y]\notag \\ &= \sum_{i \in a} D_i( [(\beta,\alpha,\cdot)])(y). \end{align}

Proof of Theorem 4.8.

By Lemma 4.1, all summands of hG are Y-harmonic on $\Omega_Y-\bigcup_{j=2}^d \partial_j$ because $(\beta,\alpha_j)$, $ j \in V_G$ are all on the characteristic surface ${\mathcal H}$. It remains to show that hG is Y-harmonic on all $\partial_a \cap \Omega_Y$ $a \subset \{2,3,4,...,d\}$ and $a \neq \varnothing.$ We will do this by induction on $|a|$. Let us start with $|a|=1$, i.e., $a = \{l\}$, for some $l \in \{2,3,4,..,d\}$ Take any vertex $i \in V_G$; if $l \in G(i,i)$ then $(\beta,\alpha_i) \in {\mathcal H}_l$ and by Lemma 4.4 $[(\beta,\alpha_i),\cdot]$ is Y-harmonic on $\partial_l$. Otherwise, the definition of a harmonic system implies that there exists a unique vertex j of G such that $G(i,j) = l$. This implies, by definition, that $(\beta,\alpha_i)$ and $(\beta,\alpha_{j})$ are l-conjugate and by Proposition 4.5 and (103)

\begin{equation*} {\boldsymbol c}_i [(\beta,\alpha_i),\cdot] + {\boldsymbol c}_{j} [ (\beta,\alpha_{j}),\cdot] \end{equation*}

is Y-harmonic on $\partial_l$. Thus, all summands of hG are either Y-harmonic on $\partial_l$ or form pairs which are so; this implies that the sum hG is Y-harmonic on $\partial_l$.

Now assume hG is Y-harmonic for all aʹ with $|a'| = k-1$; fix an $a \subset \{2,3,4,...,d\}$ such that $|a| = k$ and a $i \in a$; by (105)

\begin{equation*} D_a( h_G) = D_{a-\{i\}}(h_G) + D_i(h_G). \end{equation*}

The induction assumption and Lemma 4.2 imply that the first term on the right is zero; the same lemma and the previous paragraph imply the same for the second term. Then $D_a(h_G) = 0$; this and Lemma 4.2 finish the proof of the induction step.

4.1. Simple extensions

In this subsection, we show how the solution of a harmonic system for a lower dimensional process can provide solutions for a related harmonic system of a higher dimensional process provided that the higher dimensional process is a “simple extension” (defined below) of the lower dimensional one.

For two integers $d_2 \gt d_1 \gt 0 $ let $p_i \in {\mathbb R}^{(d_i + 1) \times (d_i +1)}$, $i=1,2$, be two transition matrices. Define $p' \in {\mathbb R}^{(d_1 +1) \times (d_1+1)}$ as

(106)\begin{align} p'(i,j) &= p_2(i,j) \end{align}

if $i\in \{0,1,2,3,...,d_1\}$, $j \in \{1,2,3,...,d_1\}$ and

(107)\begin{align} p'(i,0) &= p_2(i,0) + \sum_{j=d_1+1}^{d_2} p_2(i,j), ~~~i\in \{1,2,3,...,d_1\}. \end{align}

Definition 4.9. We say that p 2 is a simple extension of p 1 if

(108)\begin{align} &p' = \left(\sum_{i,j=0}^{d_1} p'(i,j)\right) p_1 , p' \neq 0, \end{align}
(109)\begin{align} &p_2(i,j) = 0 \ \text{if }i \in \{d_1 +1,...,d_2\}, j \in \{1,2,3,...,d_1\}. \end{align}

An example:

(110)\begin{equation} p_1 = \left( \begin{matrix} 0 & 1/7 & 0 \\ 0 & 0 & 4/7 \\ 2/7 & 0 & 0 \end{matrix} \right), p_2 = \left( \begin{matrix} 0 & 0.05 & 0 & 0 & 0.02\\ 0 & 0 & 0.2 & 0 & 0 \\ 0.05 & 0 & 0 & 0.05 & 0 \\ 0 & 0 & 0 & 0 & 0.25\\ 0.2 & 0 & 0 & 0.18 & 0 \end{matrix} \right). \end{equation}

Figure 3 shows the topologies of the networks corresponding to p 1 and $p_2.$

Definition 4.10. Let G be a L-regular. Let ${L}_1 \supset L$ be another set of labels. G’s simple extension G1 to an L1-regular graph is defined as follows: $V_{G_1} = V_G$ and

(111)\begin{align} G_1(i,j) &= G(i,j), i \neq j, i,j\in V_G\\ G_1(j,j) &= G(j,j) \cup (L_1 - L), j \in V_G.\notag \end{align}

To get G 1 from G one adds to each vertex of G an l-loop for each $l \in L_1 - L$. G is L-regular implies that G 1 is L 1-regular. Figure 4 gives an example.

Figure 4. A $\{2\}$-regular graph and its simple extension to a $\{2,3,4,5\}$-regular graph.

If Y 2 is a simple extension of Y 1, any solution to a Y 1-harmonic system implies a related solution to a related Y 2-harmonic system:

Proposition 4.11. For $d_2 \gt d_1 \gt 1$ let $p_i \in {\mathbb R}_+^{(d_i + 1) \times (d_i + 1)}$, $i=1,2$ be transition matrices such that p 2 is a simple extension of $p_1.$ Let Yi be defined through (90) with $d=d_i$, $i=1,2$ and $p = p_i$, $i=1,2.$ Let Gi, $i=1,2$ be $\{2,3,...,d_i\}$-regular graphs for Y 2 and Y 1 such that G 2 is a simple extension of G 1 (in the sense of Definition 4.10). Suppose $(\beta,\alpha_k), {\boldsymbol c}_k, k \in V_{G_1}$ solve the harmonic system associated with G 1. For $k \in V_{G_2} = V_{G_1}$ define $\alpha^2_k \in {\mathbb C}^{d_2 + 1}$ as follows

(112)\begin{align} \alpha^2_k(j) &= \alpha_k^1(j),~ j\in \{2,3,4,...,d_1\} \end{align}
(113)\begin{align} \alpha^2_k(j)&= \beta, j \in \{d_1+1,d_1+2,...,d_2\}. \end{align}

Then $(\beta,\alpha^2_k), {\boldsymbol c}_k, k \in V_{G_2}$ solves the harmonic system defined by G 2 and $p_2.$

The definition (113) extends $\alpha^1_k \in {\mathbb C}^{d_1-1}$ to $\alpha^2_k \in {\mathbb C}^{d_2-1}$ by assigning the value β to the additional dimensions of $\alpha^2_k$. This, (108) and (109) imply that, when $\alpha^2_k$ is defined as above, the harmonic system defined by G 2 reduces to that defined by G 1; the details are as follows:

Proof. By assumption, $(\beta,\alpha_k^1)$, ${\boldsymbol c}_k$, $k \in V_{G_0}$ satisfy the five conditions listed under Definition 4.7 for $G=G_1$ and $p=p_1$. We want to show that this implies that the same holds for $(\beta,\alpha^2_k), {\boldsymbol c}_k, k \in V_{G_2}$ for $G=G_2$ and $p=p_2.$

Fix any $k \in V_{G_1}$; (112) and (113) imply

(114)\begin{equation} [(\beta,\alpha^2_k),v_{i,j}^2 ] = \begin{cases} [(\beta,\alpha_k^1),v_{i,j}^1 ], &\ \text{if }j \le d_1,\\ [(\beta,\alpha_k^1),v_{i,0}^1 ], &\ \text{if }j \gt d_1, \end{cases} \end{equation}

for all $i \in \{0,1,2,...,d_1\}$, $j \in \{0,1,2,...,d_2\}$, $i\neq j.$ Similarly, (113) implies

(115)\begin{equation} [(\beta,\alpha^2_k), v^2_{i,j} ] = 1 \end{equation}

for all $i,j \in \{0,d_1 + 1,d_1 +2,...,d_2\}$, $i\neq j.$

Let ${\mathbf p}^2$ denote the characteristic polynomial of Y 2 and let ${\mathcal H}^2$ denote its characteristic surface; we would like to show $(\beta,\alpha_k^2) \in {\mathcal H}^2$, i.e., ${\mathbf p}^2(\beta,\alpha_k^2) = 1.$

By (109)

(116)\begin{align} {\mathbf p}^{2}(\beta,\alpha^2_k) &= \sum_{i=0,j=1}^{d_1} p_2(i,j) [(\beta,\alpha^2_k), v^2_{i,j} ] + \sum_{i=1,j \in\{0,d_1+1,...,d_2\}}^{d_1} p_2(i,j) [(\beta,\alpha^2_k), v^2_{i,j} ]\\ &~~~+ \sum_{i,j \in\{0,d_1+1,...d_2\}} p_2(i,j) [(\beta,\alpha^2_k), v^2_{i,j} ].\notag \end{align}

(106), (114) and (115) imply

(117)\begin{align} &= \sum_{i=0,j=1}^{d_1} p'(i,j) [(\beta,\alpha_k^1), v^1_{i,j} ] + \sum_{i=1, j \in \{0,d_1 +1,...,d_2\} }^{d_1} p_2(i,j) [(\beta,\alpha_k^1), v^1_{i,0} ]\\ &~~~+ \sum_{i,j \in \{0,d_1+1,...,d_2\}} p_2(i,j) \notag \end{align}

(107) implies that the second sum above equals $\sum_{i=1}^{d_1} p'(i,0) [(\beta,\alpha_k^1), v^1_{i,0}].$ Substitute this back in (117) to get

\begin{align*} {\mathbf p}^{2}(\beta,\alpha^2_k) &= \sum_{i,j=0}^{d_1} p'(i,j) [(\beta,\alpha_k^1), v_{i,j} ] + \sum_{i,j \in\{0,d_1+1,...,d_2 } p_2(i,j) \end{align*}

which, by (108), equals

\begin{align*} & = \left(\sum_{i,j =0 }^{d_1} p'(i,j) \right) \sum_{i,j=0}^{d_1}p_1(i,j)[(\beta,\alpha_k^1), v^{1}_{i,j} ] + \sum_{i,j \in \{0,d_1+1,...,d_2\}} p_2(i,j) \end{align*}

$(\beta,\alpha_k) \in {\mathcal H}$, (106), (107) and (109) now give

\begin{align*} & = \sum_{i,j=0}^{d_2} p_2(i,j) =1 \end{align*}

i.e., $(\beta,\alpha_k^2) \in {\mathcal H}^2.$ This proves $(\beta,\alpha^2_k) \in {\mathcal H}^{2}$, $k \in V_{G_2}$,i.e., the first part of Definition 4.7 is satisfied by $(\beta,\alpha^2_k), {\boldsymbol c}_k, k \in V_{G_2}$ for $G=G_2$ and $p=p_2.$

By definition $\alpha_{i}^1 \neq \alpha_{j}^1$ for ij, this and (112) imply $\alpha_{i}^2 \neq \alpha_{j}^2$, i.e, the second part of Definition 4.7 also holds for $(\beta,\alpha^2_k), {\boldsymbol c}_k, k \in V_{G_2}$ for $G=G_2$ and $p=p_2.$

Let us now show that the third part of the same definition is also satisfied. Fix any ij with $G_2(i,j) =l \in \{2,3,4,....,d_1\}$ (that G 2 is a simple extension of G 1 means that $G_2(i,j) \in \{2,3,4,..,d_1\}$; see (111)). We want to show that $(\beta,\alpha_{i}^2)$ and $(\beta,\alpha_{j}^2)$ are l-conjugate, , i.e., that they satisfy (100) and (102):

(118)\begin{equation} \alpha_{i}^2(k) = \alpha_j^2(k), k \neq l, \end{equation}
(119)\begin{equation} \alpha_{i}^2(l) \alpha_j^2(l) =\frac{\sum_{k=0}^{d_2} p_2(l,k) [(\beta,\alpha_i^2\{l\}),v_{l,k}^2]} {\sum_{k=0}^{d_2} p_2(k,l) [(\beta, \alpha_i^2\{l\}), v_{k,l}^2]}. \end{equation}

By definition $G_2(i,j)=l$ when $G_1(i,j) =l$; $G_1(i,j)=l$ implies that $\alpha_i^1$ and $\alpha_j^1$ are l-conjugate; in particular, they satisfy (100). (118) follows from this, (112) and (113).

We next prove (119). For $l \in \{2,3,4,...d_1\}$, $\alpha_j^2(l) = \alpha^1_j(l)$ and $\alpha_i^2(l) = \alpha^2_i(l)$; therefore $\alpha_i^2(l) \alpha_j^2(l) = \alpha^1_i(l)\alpha^1_j(l)$. $\alpha^1_i$ and $\alpha^1_j$ are l-conjugate, in particular, they satisfy (102):

\begin{equation*} \alpha_{i}^1(l) \alpha_j^1(l) =\frac{\sum_{k=0}^{d_1} p_1(l,k) [(\beta,\alpha_i^1\{l\}),v_{l,k}^1]} {\sum_{k=0}^{d_1} p_1(k,l) [(\beta, \alpha_i^1\{l\}), v_{k,l}^1]}. \end{equation*}

Then to prove (119) it suffices to prove

(120)\begin{equation} \frac{\sum_{k=0}^{d_2} p_2(l,k) [(\beta,\alpha_i^2\{l\}),v^2_{l,k}]} {\sum_{k=0}^{d_2} p_2(k,l) [(\beta, \alpha_i^2\{l\}), v^2_{k,l}]} =\frac{\sum_{k=0}^{d_1} p_1(l,k) [(\beta,\alpha_i^1\{l\}),v^1_{l,k}]} {\sum_{k=0}^{d_1} p_1(k,l) [(\beta, \alpha_i^1\{l\}), v^1_{k,l}]}. \end{equation}

This follows from a decomposition parallel to the one given in (116); let us first apply it to the numerator:

\begin{align*} \sum_{k=0}^{d_2}& p_2(l,k) [(\beta,\alpha_i^2\{l\}),v_{l,k}^2] =\sum_{k=0}^{d_1} p_2(l,k) [(\beta,\alpha_i^2\{l\}),v_{l,k}^2] + \sum_{k=d_1+1}^{d_2} p_2(l,k) [(\beta,\alpha_i^2\{l\}),v_{l,k}^2] \notag \end{align*}

(106), (107), (108) and (114) imply

(121)\begin{align} &~~~=\sum_{k=1}^{d_1} p'(l,k) [(\beta,\alpha_i^1\{l\}),v_{l,k}^1] + p_2(l,0) [(\beta,\alpha_i^1\{l\}),v_{l,0}^1] + \sum_{k={d_1+1}}^{d_2} p_2(l,k) [(\beta,\alpha_i^1\{l\}),v_{l,0}^1] \notag \\ &~~~= \left(\sum_{i,j=0 }^{d_1} p'(i,j) \right) \sum_{k=0}^{d_1} p_1(k,l) [(\beta,\alpha_i^1\{l\}),v_{k,l}^1]. \end{align}

A parallel argument for the denominator gives (this time also using (109))

\begin{equation*} \sum_{k=0}^{d_2} p_2(k,l) [(\beta,\alpha_i^2\{l\}),v^2_{k,l}] = \left(\sum_{i,j=0}^{d_1} p'(i,j) \right) \sum_{k=0}^{d_1} p_1(k,l)[(\beta,\alpha_i^1\{l\}),v_{k,l}^1]. \end{equation*}

Dividing (121) by the last equality gives (120).

The proof that parts 4-5 of Definition 4.7 hold for $(\beta,\alpha^2_k), {\boldsymbol c}_k, k \in V_{G_2}$ for $G=G_2$ and $p=p_2$ is parallel to the arguments just given and is omitted.

In the following remark, we note several facts that we don’t need directly in our arguments. Their proofs are very similar to the arguments given above and are left to the reader:

Remark 4.12. Let Y 1 and Y 2 be as above, i.e, Yi is di dimensional, $d_1 \lt d_2$ and Y 2 is a simple extension of Y 1; for $y \in {\mathbb Z}^{d_2}$, let $y^{1,d_1}$ be denote the projection of y onto its first d 1 coordinates. If h is Y 1-harmonic then, $y\mapsto h(y^{1,d_1})$ is Y 2-harmonic. Similarly, let G 1 and G 2, ${\boldsymbol c}_k$, $\alpha_k^i$, $k \in V_{G_1}$ be as in the previous proposition; then $h_{G_2}(y) = h_{G_1}(y^{1,d_1})$.

4.2. $\partial B$-determined Y-harmonic functions

A Y-harmonic function h is said to be $\partial B$-determined if

(122)\begin{equation} h(y) = {\mathbb E}_y[h(Y_{\tau}) 1_{\{\tau \lt \infty \}}], y \in \Omega_Y. \end{equation}

$y \mapsto {\mathbb{P}}_y(\tau \lt \infty)$ is the unique $\partial B$-determined Y-harmonic function with the value 1 on $\partial B$. The next proposition identifies simple conditions under which a Y-harmonic function defined by a harmonic system is $\partial B$-determined.

Proposition 4.13. Let $(\beta,\alpha_j), {\boldsymbol c}_j$ be the solutions of a Y-harmonic system with its graph G and let hG be defined as in (104). If

\begin{equation*} |\beta| \lt 1,~~~|\alpha_j(i)| \le 1, i=2,3,...,d, j \in V_G, \end{equation*}

then hG is $\partial B$-determined.

The proof is identical to that of [Reference Sezer65, Proposition 5]; for ease of reference we give an outline below:

Proof. Define $\xi_n = \inf\{k: Y_k(1) = \sum_{j=2}^d Y_k(j) + n \}.$ The optional sampling theorem and the fact that hG is Y-harmonic imply

\begin{equation*} h_G(y) = {\mathbb E}_y[ h_G(Y_\tau) 1_{\{\tau \le \xi_n\}}] + {\mathbb E}_y[h_G(Y_{\xi_n}) 1_{\{\xi_n \le \tau\}}]. \end{equation*}

$|\alpha_i| \le 1$ implies $|{\mathbb E}_y[h_G(Y_{\xi_n}) 1_{\{\xi_n \le \tau\}}]| \le \beta^n |V_G|\max_{j\in V_G} |{\boldsymbol c}_j|$. This, $|\beta| \lt 1$ and letting $n\rightarrow \infty$ in the last display give

\begin{equation*} h_G(y) = {\mathbb E}_y[ h_G(Y_\tau) 1_{\{\tau \lt \infty \}}]. \end{equation*}

5. Harmonic systems for constrained random walks representing tandem networks and the computation of ${\mathbb{P}}_y(\tau $ < $ \infty)$

Throughout this section, we will denote the dimension of the system with ${\boldsymbol d}$; the arguments below for ${\boldsymbol d}$ dimensions require the consideration of all walks with dimension $d \le {\boldsymbol d}$.

We will now define a specific sequence of regular graphs for tandem walks and construct a particular solution to the harmonic system defined by these graphs. These particular solutions will give us an exact formula for ${\mathbb{P}}_y(\tau \lt \infty)$ in terms of the superposition of a finite number of log-linear Y-harmonic functions.

We will assume

(123)\begin{equation} \mu_i \neq \mu_j, i \neq j; \end{equation}

this generalizes $\mu_1 \neq \mu_2$ assumed in [Reference Sezer65]. One can treat parameter values which violate (123) by taking limits of the results of the present section, we give an example in Section 9.

The characteristic polynomials for the tandem walk are:

(124)\begin{align} {\mathbf p}(\beta,\alpha) &= \lambda \frac{1}{\beta} + \mu_1 \alpha(2) + \sum_{j=2}^{{\boldsymbol d}} \mu_j\frac{\alpha(j+1)}{\alpha(j)}\\ {\mathbf p}_i(\beta,\alpha) &= \lambda \frac{1}{\beta} + \mu_1 \alpha(2) + \mu_i + \sum_{j=2, j \neq i}^{{\boldsymbol d}} \mu_j\frac{\alpha(j+1)}{\alpha(j)}, \notag \end{align}

where by convention $\alpha({\boldsymbol d}+1) = \beta$ (this convention will be used throughout this section, and in particular, in Lemma 5.1, (125), and (126)).

(124) implies

Lemma 5.1. For $j \in \{2,3,4,...,{\boldsymbol d} \}$, $(\beta,\alpha) \in {\mathcal H} \cap {\mathcal H}_j$ $\iff$ $\mu_j \frac{\alpha(j+1)}{\alpha(j)} = \mu_j$ $\iff$ $\alpha(j+1) = \alpha(j)$.

For the tandem walk, the conjugacy relation (101) reduces to

(125)\begin{equation} \alpha_1(i)\alpha_2(i)= \begin{cases} \frac{\alpha(3)\mu_2}{\mu_1}, & i =2, \\ \frac{\alpha(i-1)\alpha(i+1)\mu_i}{\mu_{i-1}}, &i=2,3,...,{\boldsymbol d}. \end{cases} \end{equation}

For tandem walks the functions $C(j, \beta,\alpha)$ of (96) reduce to

(126)\begin{equation} C(j,\beta,\alpha) = \mu_j - \mu_j \frac{\alpha(j+1)}{\alpha(j)}, \end{equation}

We define $\{2,3,...,{\boldsymbol d}\}$-regular graphs $G_{d,{\boldsymbol d}}$, $ d \in \{1,2,3,...,{\boldsymbol d}\}$ as follows:

(127)\begin{equation} V_{G_{d,{\boldsymbol d}}} = \{a\cup\{d\}, a \subset \{1,2,3,...,d-1\}\}; \end{equation}

for $j \in (a \cup \{d\})$, j ≠ 1, define $G_{d,{\boldsymbol d}}$ by

(128)\begin{align} G_{d,{\boldsymbol d}}( a\cup\{d\}, a\cup \{d\} \cup \{j-1\}) &= j \ \text{if } j-1 \notin a \end{align}

and

(129)\begin{equation} G_{d,{\boldsymbol d}}( a\cup \{d\}, a \cup \{d\}) = \{2,3,4,...,d\}- a \cup \{d\}; \end{equation}

these and its symmetry determine $G_{d,{\boldsymbol d}}$ completely. We note that vertices of $G_{d,{\boldsymbol d}}$ are subsets of $\{1,2,3,...,{\boldsymbol d}\}$; we will assume these sets to be sorted, for $a \subset \{1,2,3,...,{\boldsymbol d}\}$, a(1) denotes the smallest element of a, $|a|$ the number of elements in a and $a(|a|)$ the greatest element of a. Figure 5 shows the graph $G_{4,4}$.

Figure 5. $G_{d,{\boldsymbol d}}$ for $d={\boldsymbol d}=4$.

The next two propositions follow directly from the above definition:

Proposition 5.2. $G_{d,{\boldsymbol d}}$ is the simple extension of $G_{d,d}$ to a $\{2,3,...,{\boldsymbol d}\}$-regular graph.

Let $G_{d+1,{\boldsymbol d}}^{k}$ denote the subgraph of $G_{d+1,{\boldsymbol d}}$ consisting of the vertices $\{a, k,d+1\}$, $a\subset \{1,2,3,...,k-1\}$.

Proposition 5.3. One can represent $G_{d+1,{\boldsymbol d}}$ as a disjoint union of the graphs $G_{k,{\boldsymbol d}}, k=1,2,..,d,$ and the vertex $\{d+1\}$ as follows: for $a \subset \{1,2,3,...,k-1\}$ map the vertex $a \cup \{k\}$ of $G_{d+1,{\boldsymbol d}}$ to $a \cup \{k, d+1\}$. This maps $G_{k,{\boldsymbol d}}$ to the subgraph $G_{d+1,{\boldsymbol d}}^{k}$ of $G_{d+1,{\boldsymbol d}}$ consisting of the vertices $a \cup\{k,d+1\}$, $a\subset \{1,2,3,...,k-1\}$. The same map preserves the edge structure of $G_{k,{\boldsymbol d}}$ as well except for the d + 1-loops. These loops on $G_{k,{\boldsymbol d}}$ are broken and are mapped to d + 1-edges between $G^{k}_{d+1,{\boldsymbol d}}$ and $G^{d-1}_{d+1,{\boldsymbol d}}$.

Figure 5 shows an example of the decomposition described in Proposition 5.3.

For $a \subset \{2,3,4,...,{\boldsymbol d}\}$ define

(130)\begin{align} {\boldsymbol c}^*_{a} &:= (-1)^{|a|-1} \prod_{j=1}^{|a|-1} \prod_{l = a(j) + 1}^{a(j+1)} \frac{\mu_l - \lambda}{\mu_l-\mu_{a(j)}}\notag\\ \alpha^*_{a}(l) &:= \begin{cases} 1 &\ \text{if } l \le a(1)\\ \rho_{a(j)}, &\ \text{if } a(j) \lt l \le a(j+1),\\ \rho_{a(|a|)} & \ \text{if } l \gt a(|a|), \end{cases}\\ \beta^*_{a}&:= \rho_{a(|a|)}, \notag \end{align}

$l \in \{2,3,...,{\boldsymbol d}\}$ (remember that we assume sets ordered and $a(|a|)$ denotes the largest element in the set). Let us give several examples to these definitions for ${\boldsymbol d} = 8$:

(131)\begin{align} {\boldsymbol c}^*_{\{5\}}&= 1\notag \\ \alpha^*_{\{5\}} &= (1,1,1,1,\rho_5,\rho_5,\rho_5)\notag \\ {\boldsymbol c}^*_{\{3,6\}} &= - \frac{\mu_4 -\lambda}{\mu_4 - \mu_3} \frac{\mu_5 - \lambda}{\mu_5-\mu_3} \frac{\mu_6-\lambda}{\mu_6-\mu_3}.\notag \\ {\boldsymbol c}^*_{\{3,5,7\}} & = (-1)^2 \frac{\mu_4 - \lambda}{\mu_4 -\mu_3} \frac{\mu_5-\lambda}{\mu_5-\mu_3} \frac{\mu_6-\lambda}{\mu_6-\mu_5} \frac{\mu_7-\lambda}{\mu_7-\mu_5}\notag \\ \alpha^*_{\{3\}} &= (1,1,\rho_3,\rho_3,\rho_3,\rho_3,\rho_3)\\ \alpha^*_{\{3,6\}} &= (1,1,\rho_3,\rho_3,\rho_3,\rho_6,\rho_6)\notag \\ \alpha^*_{\{3,5,7\}} &= (1,1,\rho_3,\rho_3,\rho_5,\rho_5,\rho_7)\notag \\ \alpha^*_{\{8\}} &= (1,1,1,1,1,1,1);\notag \end{align}

remember that we index the components of $\alpha^*$ with $\{2,3,4,...,{\boldsymbol d}\}$; therefore, e.g., the first 1 on the right side of the last line is $\alpha^*_{\{8\}}(2).$

It follows from (130) and (130) that

\begin{align*} {\boldsymbol c}^*_{a \cup \{d_1,d_2\}} &=-{\boldsymbol c}^*_{a \cup \{d_1\}} \prod_{l=d_1+1}^{d_2} \frac{\mu_l - \lambda}{\mu_l - \mu_{d_1}} \\ \alpha^*_{a \cup \{d_1\}} &= \alpha^*_{a \cup \{d_1,{\boldsymbol d}\}} \end{align*}

for any $1 \lt a(|a|) \lt d_1 \lt d_2 \le {\boldsymbol d}$ and $a \subset \{2,3,4,...,{\boldsymbol d}\}$; These and Proposition 5.3 imply

Proposition 5.4. For $d \lt {\boldsymbol d}$ and $y \in \partial B$

(132)\begin{equation} -\left(\prod_{l=d+1}^{{\boldsymbol d}} \frac{\mu_l - \lambda}{\mu_l - \mu_d}\right) \sum_{a \in V_{G_{d,{\boldsymbol d}} }} {\boldsymbol c}_a^* [ (\beta_a^*, \alpha_a^*), y ] = \sum_{a \in V_{G_{{\boldsymbol d},{\boldsymbol d}}^d }} {\boldsymbol c}_a^* [ (\beta_a^*, \alpha_a^*), y ] \end{equation}

Proposition 5.5. For $d \le \boldsymbol d$, let $G_{d,{\boldsymbol d}}$ be as in (127) and (128). Then $(\beta^*_{a\cup\{d\}}, \alpha^*_{a\cup\{d\}})$, ${\boldsymbol c}^*_{a\cup\{d\}}$, $a \subset \{1,2,3,...,d-1\}$, defined in (130), solve the harmonic system defined by $G_{d,{\boldsymbol d}}$.

Proof. A dʹ tandem walk is a simple extension of the tandem walk defined by its first $d'-1$ dimensions. This, Propositions 5.2, 4.11 and the definitions of $\beta^*$ and $\alpha^*$ above imply that it suffices to prove the current proposition only for $d=\boldsymbol d$.

The vertices of $G_{{\boldsymbol d},{\boldsymbol d}}$ are $a \cup \{{\boldsymbol d}\}$, $a \subset \{1,2,3,....,{\boldsymbol d}-1\}$, and for all of them we have $\beta^*_{a\cup\{{\boldsymbol d}\}} = \rho_{\boldsymbol d}$ by definition (130). Let us begin by showing $\left(\rho_{\boldsymbol d}, \alpha^*_{a\cup\{\boldsymbol d\}}\right)$, $a \subset \{1,2,3,....,{\boldsymbol d}-1\}$ is on the characteristic surface ${\mathcal H}$ of the tandem walk. We will write $\alpha^*$ instead of $\alpha^*_{a\cup\{\boldsymbol d\}}$, the set a will be clear from context.

Let us first consider the case when $a(1) \gt 1$, i.e., when $1 \notin a$; the opposite case is treated similarly and is left to the reader. Then $\alpha^*(l) = 1$ for $2 \le l \le a(1)$. By definition $\alpha^*(i) = \alpha^*(i+1)$ if $a(j) \lt i \lt a(j+1)$; these and $\beta^*_{a\cup\{{\boldsymbol d}\}} = \rho_{\boldsymbol d}$ give

\begin{align*} {\mathbf p}(\rho_{\boldsymbol d},\alpha^*) &= \mu_{\boldsymbol d} + \sum_{j=1}^{a(1)-1} \mu_j + \mu_{a(1)} \rho_{a(1)} + \sum_{j \in (a^c -\{1\cdots a(1)-1 \})} \mu_j\\ &~~~~~+ \sum_{j \in (a-\{a(1)\}) } \mu_j \frac{\alpha^*(j+1)}{\alpha^*(j)} + \rho_{\boldsymbol d} \frac{\mu_{\boldsymbol d}}{\alpha^*(\boldsymbol d)} \end{align*}

(where $a^c= \{1,2,3,...,{\boldsymbol d}-1\} - a$) and in the last expression we have used the convention $\alpha^*({\boldsymbol d} + 1) =\beta^*$; by definition (130) $\alpha^*(a(j+1)) = \rho_{a(j)}$, $\alpha^*(a(j)) = \rho_{a(j-1)}$ and therefore

\begin{align*} &= \mu_{\boldsymbol d} + \sum_{j=1}^{a(1)-1} \mu_j + \lambda + \sum_{j\in (a^c -\{1\cdots a(1)-1 \})} \mu_j\\ &~~~~~+ \sum_{j=2}^{|a|} \mu_{a(j)} \frac{\rho_{a(j)}}{\rho_{a(j-1)}} + \mu_{a(|a|)} \end{align*}

$\mu_{a(j)} \rho_{a(j)}/\rho_{a(j-1)} = \mu_{a(j-1)}$ implies

\begin{align*} &= \mu_{\boldsymbol d} + \sum_{j=1}^{a(1)-1} \mu_j + \lambda + \sum_{j\in (a^c -\{1\cdots a(1)-1\})} \mu_j\\ &~~~~~+ \sum_{j=2}^{|a|} \mu_{a(j-1)} + \mu_{a(|a|)}\\ &= \mu_{\boldsymbol d} + \sum_{j=1}^{a(1)-1} \mu_j + \lambda + \sum_{j\in (a^c -\{1\cdots a(1)-1 \})} \mu_j+ \sum_{j \in a }\mu_j = 1; \end{align*}

i.e., $(\rho_{\boldsymbol d}, \alpha^*) \in {\mathcal H}.$

If $a_1 \neq a_2$ take any $ i \in a_1 -a_2$ (relabel the sets if necessary so that $a_1 -a_2 \neq \varnothing$). Let j be the index of i in a 1, i.e., $a_1(j) = i$. Then by definition, $\alpha^*_{a_1\cup\{\boldsymbol d\}}(j+1) = \rho_i$; but $i \notin a_2$ and (123) imply that no component of $\alpha^*_{a_1\cup \{\boldsymbol d\}}$ equals ρi, and therefore $\alpha^*_{a_1\cup\{\boldsymbol d\}} \neq \alpha^*_{a_2\cup\{\boldsymbol d\}}$. This shows that $\alpha^*_{a \cup \{\boldsymbol d\}}$, $a \subset \{1,2,...,{\boldsymbol d}-1\}$ satisfy the second part of Definition 4.7.

Fix a vertex $a \cup \{\boldsymbol d\}$ of $G_{{\boldsymbol d},{\boldsymbol d}}$. By definition, for each of its elements l, this vertex is connected to $ a \cup \{\boldsymbol d\} \cup \{l-1\}$ if $l-1 \notin a$ or to $ a \cup \{\boldsymbol d\} - \{l-1\}$ if $l-1 \in a$. Then to show that $(\beta^*_{a \cup \{{\boldsymbol d} \}}, \alpha^*_{a \cup \{\boldsymbol d\}})$, $a \subset \{1,2,3,....,{\boldsymbol d}-1\}$ satisfy the third part of Definition 4.7 it suffices to prove that for each $a \subset \{1,2,3,...,{\boldsymbol d}-1\}$, and each $l \in a \cup \{\boldsymbol d\}$ such that $l-1 \notin a \cup \{\boldsymbol d\}$ $\alpha^*_{a \cup \{\boldsymbol d\}}$ and $\alpha^*_{a \cup \{\boldsymbol d\} \cup \{l-1\}}$ are l-conjugate. For ease of notation let us denote $a \cup \{l-1\}$ by a 1, $\alpha^*_{a\cup \{\boldsymbol d\}}$ by $\alpha^*$, $\alpha^*_{a_1 \cup \{\boldsymbol d\}}$ by $\alpha^*_1$ and $\beta^*_{a_1} =\beta^*_a$ by $\beta^*$ (because we have assumed $d={\boldsymbol d}$, both $\beta^*$ and $\beta^*_1$ are equal to $\rho_{\boldsymbol d}$). We want to show that $(\beta^*, \alpha^*)$ and $(\beta^*, \alpha^*_1)$ are l-conjugate. Let us assume $2 \lt l \lt {\boldsymbol d}$, the cases $l=2,{\boldsymbol d}$ are treated almost the same way and are left to the reader. By assumption $l \in \alpha^*$ but $l-1 \notin \alpha^*$. If l is the j th element of a, i.e., $l = a(j)$; then $a(k) =a_1(k)$ for k < j, $a_1(j)=l-1$, $a(k-1) = a_1(k)$ for k > j. This and the definition (130) of $\alpha^*$ imply

(133)\begin{equation} \alpha^*(i) = \alpha^*_1(i), i \in \{2,3,4,...,{\boldsymbol d} \}, i \neq l, \end{equation}

i.e., $\alpha^*$ and $\alpha^*_1$ satisfy (100) (for example, for $\boldsymbol d=8$, $\alpha^*_{\{3,6\}}$ is given in (131); on the other hand $\alpha^*_{\{3,5,6\}} = (1,1,\rho_3,\rho_3,\rho_5,\rho_6,\rho_6)$ and indeed $ \alpha^*_{\{3,6\}}(i) = \alpha^*_{\{3,5,6\}}(i), i \neq 6$). Definition (130) also implies

(134)\begin{equation} \alpha^*_1(l) = \rho_{a_1(j)} = \rho_{l-1},~ \alpha^*_1(l+1) = \rho_{a_1(j+1)} = \rho_l, \end{equation}

On the other hand, again by (130), and by $l-1 \notin a$, we have

\begin{equation*} \alpha^*(l) = \alpha^*(l-1) = \rho_{a(j-1)} \ \text{and } \alpha^*(l+1) = \rho_l. \end{equation*}

Then

\begin{equation*} \frac{1}{\alpha^*(l)} \frac{\alpha^*(l-1) \alpha^*(l+1) \mu_l}{\mu_{l-1}} = \rho_{l-1} \end{equation*}

and, by (134) this equals $\alpha^*_1(l)$, i.e., $\alpha^*_1$ and $\alpha^*$ satisfy (125). This and (133) mean that $(\beta^*,\alpha_1^*)$ and $(\beta^*,\alpha)$ are l-conjugate.

Now we will prove that the ${\boldsymbol c}^*_{a\cup\{\boldsymbol d\}}$, $a \subset \{2,3,4,...,{\boldsymbol d}-1\}$ defined in (130) satisfy the fourth part of Definition 4.7. The structure of $G_{{\boldsymbol d},{\boldsymbol d}}$ implies that it suffices to check that

(135)\begin{equation} \frac{{\boldsymbol c}^*_{a}}{{\boldsymbol c}^*_{a_1} } = -\frac{C(l',\rho_{\boldsymbol d}, \alpha^*_{a_1})} {C(l',\rho_{\boldsymbol d},\alpha^*_a)} \end{equation}

holds for any $ l' \in a$ such that $l'-1 \notin a$ and $a_1 = a \cup \{l'-1\}.$ There are three cases to consider: $l'=2$, $l'=\boldsymbol d$ and $2 \lt l' \lt \boldsymbol d$; we will only treat the last as the rest are similar and simpler. For $2 \lt l' \lt \boldsymbol d$ one needs to further consider the cases $a(1) = l'$ and $a(1) \lt l'$. For $b \subset \{2,3,4,...,{\boldsymbol d}-1\}$, ${\boldsymbol c}^*_{b\cup\{\boldsymbol d\}}$ of (130) is the product of a parity term and a running product of $d-b(1)$ ratios of the form $(\mu_l - \lambda)/(\mu_l - \mu_{b(j)}).$ The ratio of the parity terms of a and a 1 is −1 because a 1 has one additional term. If $a(1) = l'$ then $a_1(1) = l'-1$ and the only difference between the running products in the definitions of ${\boldsymbol c}^*$ and ${\boldsymbol c}^*_1$ is that the latter has an additional initial term $(\mu_{l'}-\lambda)/(\mu_{l'}-\mu_{l'-1})$ and therefore

\begin{equation*} \frac{{\boldsymbol c}^*_{a}}{{\boldsymbol c}^*_{a_1}} = -\frac{\mu_{l'}-\mu_{l'-1}}{\mu_{l'}-\lambda}. \end{equation*}

Because $l' \gt 2$ and $l'-1 \ge 2$, definition (130) implies $\alpha^*(l') =1$, $\alpha^*(l'+1) = \rho_l$, $\alpha^*_1(l) = \rho_{l-1}$ and $\alpha^*_1(l+1) = \rho_l$. These and (126) imply

\begin{equation*} \frac{C(l,\rho_{\boldsymbol d}, \alpha^*_{a_1})}{C(l,\rho_{\boldsymbol d},\alpha^*_a)} = \frac{\mu_{l'} - \mu_{l'-1}}{\mu_{l'} - \lambda}. \end{equation*}

The last two display imply (135) for $a(1) = l'$.

If $l' \gt a(1)$, let j > 1 be the position of l in a, i.e., $l=a(j)$. In this case, the definition (130) implies that the running products in the definitions of ${\boldsymbol c}^*_{a}$ and ${\boldsymbol c}^*_{a_1}$ are a product of the same ratios except for the $(l')^{th}$ terms, which is $(\mu_{l'} - \lambda)/(\mu_{l'} -\mu_{a(j-1)})$ for a and $(\mu_{l'} - \lambda)/(\mu_{l'} -\mu_{l'-1})$ for a 1. a 1 has one more element than a, therefore, the ratio of the parity terms is again −1; these imply

\begin{equation*} \frac{{\boldsymbol c}^*_{a}}{{\boldsymbol c}^*_{a_1}} = -\frac{\mu_{l'}-\mu_{l'-1}}{\mu_{l'}-\mu_{a(j-1)}}. \end{equation*}

On the other hand, $l' \in a$, j > 1, $a_1 = a \cup \{l'-1\}$ and the definition (130) imply $\alpha^*(l') = \rho^*(a(j-1))$, $\alpha^*(l'+1) = \rho_{l'}$, $\alpha_1^*(l') = \rho_{l'-1}$, and $\alpha_1^*(l'+1) = \rho_{l'}$ and therefore

\begin{equation*} \frac{C(l',\rho_{\boldsymbol d}, \alpha^*_{a_1})}{C(l',\rho_{\boldsymbol d},\alpha^*_a)} = \frac{\mu_{l'} - \mu_{l'-1}}{\mu_{l'} - \mu_{a(j-1)}}. \end{equation*}

The last two displays once again imply (135) for $l' \gt a(1).$

Consider a vertex $a \cup \{{\boldsymbol d} \}$ of $G_{{\boldsymbol d}, {\boldsymbol d}}$; by definition (129), the loops on this vertex are $\{2,3,...,{\boldsymbol d}\} - a \cup \{{\boldsymbol d}\}.$ For $l \in \{2,3,...,{\boldsymbol d}\} - a \cup \{{\boldsymbol d}\}$ the definition (130) implies

\begin{equation*} \alpha^*_{a\cup\{{\boldsymbol d}\}}(l) = \alpha^*_{a\cup\{{\boldsymbol d}\}}(l+1); \end{equation*}

we have already shown $\alpha^*_{a\cup\{d\}} \in {\mathcal H}$, then, Lemma 5.1 and the last display imply $\alpha^*_{a\cup \{\boldsymbol d\}} \in {\mathcal H}_l$ for $l \in \{2,3,...,{\boldsymbol d}\} - a \cup \{{\boldsymbol d}\}$; i.e., the last part of Definition 4.7 is also satisfied. This finishes the proof of the proposition.

Proposition 5.6.

(136)\begin{equation} h^*_d := \sum_{a \subset \{1,2,3,...,d-1\}} {\boldsymbol c}^*_{a\cup \{d\}} [(\rho_d, \alpha^*_{a\cup\{d\}}),\cdot], \end{equation}

$d=1,2,3,...,{\boldsymbol d}$, are $\partial B$-determined Y-harmonic functions.

Proof. That $h^*_d$ is Y-harmonic follows from Proposition 5.5 and Theorem 4.8. The components of $\alpha^*_{a\cup\{d\}}$, $a\subset \{1,2,3,...,d-1\}$ and $\beta^*_d = \rho_d$ are all between 0 and 1. This and Proposition 4.13 imply that $h^*_d$ are all $\partial B$-determined.

With definition (136) we can rewrite (132) as

(137)\begin{equation} -\left(\prod_{l=d+1}^{{\boldsymbol d}} \frac{\mu_l - \lambda}{\mu_l - \mu_d}\right) h_d^*(y) = \sum_{a \in V_{G_{\boldsymbol d}^d }} {\boldsymbol c}_a^* [ (\beta_a^*, \alpha_a^*), y ] \end{equation}

for $y \in {\partial B}.$

Theorem 5.7

(138)\begin{equation} {\mathbb{P}}_y(\tau \lt \infty) = \sum_{d=1}^{{\boldsymbol d}} \left( \prod_{l=d+1}^{\boldsymbol d} \frac{\mu_l -\lambda}{\mu_l - \mu_d} \right) h^*_d(y) \end{equation}

for $y \in B.$

For ${\boldsymbol d}=2$ (138) reduces to

\begin{equation*} {\mathbb{P}}_y(\tau \lt \infty) =f(y), \end{equation*}

where f is given in (18). To express (138) as a similar formula, we expand $h_d^*$ using (130) and (21) (i.e., $[ (\beta,\alpha) , y ] = \beta^{y(d)-\sum_{j=2}^n y(j)} \prod_{j=2}^d \alpha(j)^{y(j)}$):

(139)\begin{equation} h^*_d(y) =\rho_d^{y(1)-\sum_{j=2}^d y(j)} \left( \sum_{a \subset \{1,2,3,..d-1\}} (-1)^{|a|} \prod_{j=1}^{|a|} \prod_{l = a(j) + 1}^{a(j+1)} \frac{\mu_l - \lambda}{\mu_l-\mu_{a(j)}} \rho_{a(j)}^{y(l)}\right), \end{equation}

where, by convention, we set $a(|a|+1) = d$. Substituting this in (138) we get

(140)\begin{align} {\mathbb{P}}_y(\tau \lt \infty) = \sum_{d=1}^{\boldsymbol d} \left( \prod_{l=d+1}^{\boldsymbol d} \frac{\mu_l -\lambda}{\mu_l - \mu_d} \right) &\rho_d^{y(1)-\sum_{j=2}^d y(j)}\\ &\times \left( \sum_{a \subset \{1,2,3,..d-1\}} (-1)^{|a|} \prod_{j=1}^{|a|} \prod_{l = a(j) + 1}^{a(j+1)} \frac{\mu_l - \lambda}{\mu_l-\mu_{a(j)}} \rho_{a(j)}^{y(l)}\right).\notag \end{align}

Proof of Theorem 5.7.

Let ${\boldsymbol 1} \in {\mathbb C}^{\{2,3,..,{\boldsymbol d}\}}$ denote the vector with all components equal to 1. The decomposition of $G_{\boldsymbol d}$ into the single vertex ${\{\boldsymbol d}\}$ and $G_{\boldsymbol d}^d$, $d \lt {\boldsymbol d}$ implies that the right side of (138) equals

\begin{align*} &[(\rho_{\boldsymbol d}, {\boldsymbol 1}),y] + \sum_{d=1}^{{\boldsymbol d}-1} \sum_{a \in V_{G_{\boldsymbol d}^d}} {\boldsymbol c}^*_a [ ( \beta_a^*, \alpha^*_a), y ] +\sum_{d=1}^{{\boldsymbol d}-1} \left( \prod_{l=d+1}^{\boldsymbol d} \frac{\mu_l -\lambda}{\mu_l - \mu_d} \right) h^*_d(y) \end{align*}

for $y \in \partial B$; (137) implies

\begin{align*} &~~= [(\rho_{\boldsymbol d}, {\boldsymbol 1}),y], \end{align*}

which, for $y\in \partial B$, equals 1. Thus, we see that the right side of (138) equals 1 on $\partial B$. Proposition 5.6 says that the same function is $\partial B$-determined and is Y-harmonic. Then its restriction to B must be indeed equal to $y\rightarrow {\mathbb{P}}_y(\tau \lt \infty)$, $y \in B$, which is the unique function with those properties.

6. Large deviation rate of pn for an arbitrary initial point

In the next corollary to Theorems 2.1 and 5.7 we derive the exponential decay rate of pn in n for any scaled initial condition $x \in A = \{x \in {\mathbb R}_+^d: \sum_{j=1}^{\boldsymbol d} x(j) \le 1 \}$. Since the argument depends on Theorem 5.7, we assume (123) ( $\mu_i \neq \mu_j$ if ij) as well as the stability assumption (5); we also continue to denote the dimension by ${\boldsymbol d}.$

Recall the function g appearing in Theorem 2.1: $g: {\mathbb R}_+^{\boldsymbol d} \mapsto {\mathbb R}, g(x) = \min_{i=1}^{\boldsymbol d} (1-\sum_{j=1}^i x(j))\log_{\rho}\rho_i$.

Corollary 6.1. Suppose $x_n \in {\mathbb Z}_+^{\boldsymbol d}$ satisfies $x_n/n \rightarrow x \in A \subset{\mathbb R}_+^d.$ Then

(141)\begin{equation} \lim_{n \rightarrow \infty}-\frac{1}{n} \log{\mathbb{P}}_{x_n}(\tau_n \lt \tau_0) = -\log(\rho)g(x). \end{equation}

Remark 6.2. Note that the function g arises in the derivation of the lower bound (76) on pn (see also Lemma 3.14). The lower bound (76) implies that g provides us with the large deviations upper bound for pn (this is (143) below). Surprisingly, as computed in the next proof, Theorems 2.1 and 5.7 imply that g also determines the large deviations lower bound for pn.

Proof. To prove (141) it suffices to prove

(142)\begin{align} \limsup_{n\rightarrow \infty} -\frac{1}{n}\log {\mathbb{P}}_{x_n}(\tau_n \lt \tau_0) &\le -\log(\rho)g(x) \end{align}

and

(143)\begin{align} \liminf_{n\rightarrow \infty} -\frac{1}{n} \log {\mathbb{P}}_{x_n}(\tau_n \lt \tau_0) &\ge -\log(\rho)g(x). \end{align}

Let us first assume that $x \in A\subset{\mathbb R}_+^d$ satisfies $x(j) \gt 0$ for all j. From Proposition 3.12 we have the following lower bound:

(144)\begin{equation} \max_{i=1}^{\boldsymbol d} \rho_i^{n- \sum_{j=1}^i x_n(j)}-\rho^n \le {\mathbb{P}}_{x_n}(\tau_n \lt \tau_0). \end{equation}

Recall that $\rho = \max_{i=1}^{\boldsymbol d} \rho_i.$ Let $i^*$ be so that $\rho= \rho_{i^{*}}.$ The assumptions $x(i) \gt 0$ for all i, $\rho \in (0,1)$, $x_n/n \rightarrow x$ imply that there exists N such that for n > N we have

\begin{equation*} 0.5 \rho_{i^{*}}^{n-\sum_{j=1}^{i^{*}} x_n(j)} \lt \rho_{i^{*}}^{n-\sum_{j=1}^{i^{*}} x_n(j)} -\rho^n. \end{equation*}

Let $j^*$ be the maximizer in (144); that the function $z \mapsto 0.5 z - \rho^n$, $z \in {\mathbb R}$, is increasing in z implies that the last inequality continues to hold if we replace $i^*$ with $j^*$. These imply

(145)\begin{equation} \max_{i=1}^{\boldsymbol d} 0.5 \rho_i^{n- \sum_{j=1}^i x_n(j)}\le {\mathbb{P}}_{x_n}(\tau_n \lt \tau_0) \end{equation}

for $n \gt N.$ Applying $\limsup -\frac{1}{n}\log$ to both sides gives (142).

The assumption $x(j) \gt 0$ for all j means $x \in A - R_\rho$. For such x, Theorem 2.1 implies

\begin{equation*} \liminf_{n\rightarrow \infty} -\frac{1}{n} \log {\mathbb{P}}_{x_n}(\tau_n \lt \tau_0)= \liminf_{n\rightarrow \infty} -\frac{1}{n}\log {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty). \end{equation*}

Therefore, to prove (143) it suffices to prove

(146)\begin{equation} \liminf_{n\rightarrow \infty} -\frac{1}{n}\log {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty) \ge -\log(\rho)g(x). \end{equation}

By Theorem 5.7 we know

(147)\begin{equation} {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty) = \sum_{d=1}^{\boldsymbol d} \left( \prod_{l=d+1}^{\boldsymbol d} \frac{\mu_l -\lambda}{\mu_l - \mu_d} \right)h^*_d(T_n(x_n)). \end{equation}

By (139) and the definition (10) of Tn:

\begin{align*} h^*_d(T_n(x_n)) &= \rho_d^{n -\sum_{j=1}^d x_n(j)} \left( \sum_{a \subset \{1,2,3,..d-1\}} (-1)^{|a|} \prod_{j=1}^{|a|} \prod_{l = a(j) + 1}^{a(j+1)} \frac{\mu_l - \lambda}{\mu_l-\mu_{a(j)}} \rho_{a(j)}^{x_n(l)}\right). \end{align*}

write the $a=\varnothing$ term separately in the last display:

\begin{align*} &= \rho_d^{n -\sum_{j=1}^d x_n(j)} \left( 1 + \sum_{a \subset \{1,2,3,..d-1\}, a\neq \varnothing} (-1)^{|a|} \prod_{j=1}^{|a|} \prod_{l = a(j) + 1}^{a(j+1)} \frac{\mu_l - \lambda}{\mu_l-\mu_{a(j)}} \rho_{a(j)}^{x_n(l)}\right). \end{align*}

The assumptions $x(j) \gt 0$ for all j, $x_n/n \rightarrow x$ and $\rho_j \lt 1$ for all j imply that the last sum decays exponentially in n. This and the last display imply

(148)\begin{equation} 0.5 \rho_d^{n -\sum_{j=1}^d x_n(j)} \le h^*_d(T_n(x_n)) \le 2 \rho_d^{n -\sum_{j=1}^d x_n(j)}, \end{equation}

for n large. Define

\begin{equation*} {\mathscr M}' = \left \{d \le {\boldsymbol d}: \prod_{l=d+1}^{\boldsymbol d} \frac{\mu_l -\lambda}{\mu_l - \mu_d} \gt 0\right\}; \end{equation*}

${\mathscr M}'$ are the indices d for which the coefficients in the sum (147) are positive. This definition of ${\mathscr M}'$, (147) and (148) imply, for n large,

\begin{align*} {\mathbb{P}}_{T_n(x_n)}(\tau \lt \infty) &\le 2C_0\sum_{d \in {\mathscr M}'} \rho_d^{n -\sum_{j=1}^d x_n(j)}, \end{align*}

where $C_0 = \max_{d \in {\mathscr M}' } \prod_{l=d+1}^{\boldsymbol d} \frac{\mu_l -\lambda}{\mu_l - \mu_d} \gt 0.$ We bound the last term above by

\begin{align*} &\le 2C_0\sum_{d =1}^{\boldsymbol{d} } \rho_d^{n -\sum_{j=1}^d x_n(j)}. \end{align*}

Applying $\liminf_{n} -\frac{1}{n}\log$ to both sides of this inequality gives (146). This concludes the proof of (141) under the assumption $x(j) \gt 0$ for all $j \in \{1,2,3,...,{\boldsymbol d}\}.$

Let us now extend (141) to $\{x \in A: x(j) \gt 0, j \neq {\boldsymbol d} \}$, i.e., we allow $x({\boldsymbol d}) = 0.$ Suppose $x \in A$ satisfies $x(j) \gt 0$ for $j \neq {\boldsymbol d}$ and $x({\boldsymbol d}) = 0.$ Choose δ small enough so that $\sum_{j=1}^{\boldsymbol d} x(j) + \delta \lt 1.$ Consider the sequence $x_n' = x_n + \lfloor n \delta\rfloor e_n$. The dynamics of X and the Markov property of X imply:

(149)\begin{equation} \left(\lambda \prod_{j=1}^{\boldsymbol{d} -1} \mu_j \right)^{\lfloor n \delta \rfloor } {\mathbb{P}}_{x_n'}(\tau_n \lt \tau )\le {\mathbb{P}}_{x_n}(\tau_n \lt \tau ) \le {\mathbb{P}}_{x_n'}(\tau_n \lt \tau) \mu_{\boldsymbol d}^{-\lfloor n \delta \rfloor }. \end{equation}

The sequence $x_n'$ satisfies $x_n'/n \rightarrow x' = x + \delta e_n$. The limit $x' \in A$ satisfies $x'(j) \gt 0$ for all j; therefore, we know (141) holds for the sequence $x_n'.$ Then applying $\lim_{n\rightarrow \infty} -\frac{1}{n} \log $ to the upper and lower bound in (149) and letting $\delta \rightarrow 0$ we get (141) for the sequence xn and the limit x. The extension to the cases when all other components are allowed to be 0 can be proved by the same argument and induction.

7. Numerical example

Take a four dimensional tandem system with rates, for example,

\begin{equation*} \lambda = 1/18, \mu_1 = 3/18, \mu_2 = 7/18, \mu_3 = 2/18, \mu_4 = 5/18. \end{equation*}

For n = 60, and in four dimensions, the probability ${\mathbb{P}}_x(\tau_n \lt \tau_0)$ can be computed numerically by iterating the harmonic equation ${\mathbb{P}}_x(\tau_n \lt \tau_0)= {\mathbb E}_x[{\mathbb{P}}_{X_1}(\tau_n \lt \tau_0))]$. Let f(y) denote the right side of (138). Define $V_n = -\log({\mathbb{P}}_x(\tau_n \lt \tau_0))/n$ and $W_n = -\log f(T_n(x))/n$. The level curves of Vn and Wn and the graph of the relative error $(V-W)/V$ for $x= (i,j,0,0,0)$ and $x=(0,i,0,j,0)$, $i,j \le n=60$ are shown in Figure 6; qualitatively these graphs show results similar to those reported in [Reference Sezer65]: almost zero relative error across the domain selected, except for a boundary layer along the x(4)-axis, where the relative error is bounded by 0.05. The size of the boundary layer is determined by the set Rρ of (80) and Theorem 3.15.

Figure 6. Level curves and relative error in four dimensions.

Next we consider the 14-tandem queues with parameter values shown in Figure 7.

Figure 7. The service rates (blue) and the arrival rate (red) for a 14-dimensional tandem Jackson network.

For n = 60, An contains $60^{14}/14!= 8.99 \times 10^{13}$ states which makes impractical an exact calculation via iterating the harmonic equation satisfied by ${\mathbb{P}}_y(\tau_n \lt \tau_0).$ On the other hand, (138) has $2^{15}=32,768$ summands and can be quickly calculated. Define Wn as before. Its graph over $\{x: x(4) + x(14) \le 60, x(j) = 0, j \neq 4,14\}$ is depicted in Figure 8.

Figure 8. The graph of Wn over $\{x: x(4) + x(14)= 60, x(j) = 0, j \neq 4,14\}$.

For a finer approximation of ${\mathbb{P}}_{(1,0,\cdots,0)} (\tau_n \lt \tau_0)$ we use importance sampling based on Wn. With 12,000 samples IS gives the estimate $7.53 \times 10^{-20}$ with an estimated $95\%$ confidence interval $[6.57, 8.48] \times 10^{-20}$ (rounded to two significant figures). The value given by our approximation (138) for the same probability is $f((1,0,\cdots,0)) = 1.77 \times 10^{-20}$ which is approximately $1/4^{th}$ of the estimate given by IS. The large deviation estimate of the same probability is $(\lambda/\min_{i=1}^{14}(\mu_i))^{60} = 4.15\times10^{-23}$. The discrepancy between IS and (138) quickly disappears as x(1) increases. For example, for $x(1) = 4$, IS gives $2.47\times 10^{-19}$ and (138) gives $2.32 \times 10^{-19}$.

8. Literature review

There is a direct correspondence between the structures used in the LD analysis and the subsolution approach to IS estimation of pn of [Reference Dupuis, Sezer and Wang23, Reference Dupuis and Wang25, Reference Dupuis and Wang26, Reference Sezer59, Reference Sezer61] and those appearing in [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65] and in the present work. The characteristic polynomials of the present work correspond to the Hamiltonians in these works and the subsolutions are constructed from points lying on the zero sets of the Hamiltonians. Both approaches use the functions they construct in their asymptotic analysis and algorithm construction. We note several important differences: a key difference is the nature of the results: the results in the present work give very precise (exponentially decaying relative error) deterministic approximation formulas for the probability of interest whereas the works cited above on IS algorithms construct simulation algorithms with bounds on the exponential decay rate of the absolute (not relative) variance of the estimator. Secondly, the works on IS cover only the initial point x = 0, which significantly simplifies the construction of the subsolutions and the asymptotic analysis. Thirdly, the actual construction of the functions are different. The Y-harmonic functions in Section 4 and 5 are constructed from regular graphs and from conjugate points lying on the characteristic surfaces. Such structures (functions from graphs, conjugate points) do not come up in the above cited previous works. Finally, the asymptotic analysis in the present work and in [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5, Reference Sezer65] is based on an affine transformation of the problem and involves no scaling, whereas IS and classical LD analysis is based on a law of large number of scaling of the problem. An important strength of the IS algorithms developed in these works is that they can be applied to any stable Jackson network in any dimension (see [Reference Dupuis, Sezer and Wang23, Reference Dupuis and Wang26, Reference Sezer63]) with a fairly general exit boundary (this generality crucially depends on the initial point x = 0). Such a generalization for the results in the present work is not clear at all.

The work [Reference McDonald47] studies the buffer overflow of one of the nodes in a stable network. Let W denote the stable process representing the network; W is assumed r + m dimensional: the first dimension represents the node whose overflow event is to be studied, the dimensions $2,3,..,r$ represent those nodes that become unstable when the first node overflows. For n > 0, let τn be the first time the first component of W hits n let τ 0 denote the first time W hits the origin ${\boldsymbol{0}}$. Finally, let $\tau_\triangle$ denote the first time after time 0, one of the nodes from 1 to r hits 0, i.e., $\tau_\triangle = \inf\{k: k \gt 0, W \in {\small \triangle} \}$, where $\triangle = \{x: x_j = 0, \ \text{for some }, j \in \{1,2,3,...,r\}\}$; the main approximation result in [Reference McDonald47] is the following: let $\pi_{\triangle}$ be the stationary measure conditioned on $\triangle$ and ${\mathbb E}_{\pi_\triangle}$ be the expectation conditioned on W(0) having initial distribution $\pi_{\triangle}$. Let ${\boldsymbol{\tau}}_{0}$ be the first return time to 0, i.e., ${\boldsymbol{\tau}}_{0} = \inf\{k \gt 0: W_k = {\boldsymbol{0}}\}.$ [Reference McDonald47, Lemma 1.8] states, under the assumptions made in the paper,

\begin{equation*} \lim_{n\rightarrow \infty} \frac{| \pi({\boldsymbol{0}})P_{{\boldsymbol{0}}}(\tau_n \lt {\boldsymbol{\tau}}_{\boldsymbol{0}}) - \pi(\Delta) P_{\pi_{\Delta}}(\tau_n \lt \tau_\triangle )|}{ \pi({\boldsymbol{0}})P_{{\boldsymbol{0}}}(\tau_n \lt {\boldsymbol{\tau}}_{\boldsymbol{0}})} = 0. \end{equation*}

The analysis that leads to this result is based on the h-transform of W where h is a harmonic function of W away from the set Δ that is of the following form: $h(x) = e^{\alpha x_1}a(x)$; [Reference McDonald47] gives conditions under which such an h function exists based on results from [Reference Ney and Nummelin51]. [Reference McDonald47] develops the following representation for $\pi(\Delta) P_{\pi_{\Delta}}(\tau_n \lt \tau_\triangle )$:

(150)\begin{align} \pi(\Delta) P_{\pi_{\Delta}}(\tau_n \lt \tau_\triangle ) &= e^{-\alpha n} {\mathbb E}_{\pi_\Delta}[h(W(1))\Psi(W(1))]\notag\\ \Psi(x) &= {\mathbb E}_x [a^{-1}({\mathscr W}(\tau_n)) e^{-\alpha ({\mathscr W}(\tau_n)-n)} 1_{\{\tau_n \lt \tau_{\Delta}\}}], \end{align}

where ${\mathscr W}$ is the h-transform of the process W (if W is not a nearest neighbor random walk on ${\mathbb Z}^{r+m}_+$ the formula for Ψ needs to be slightly modified, for details we refer the reader to [Reference McDonald47]). For the computation of the expectation appearing in (150), [Reference McDonald47] suggests simulation. [Reference McDonald47, Section 3] treats the two-dimensional constrained random walk on ${\mathbb Z}_+^2$ with increments $(-1,0)$, $(1,0)$, $ (0,-1)$, $(0,1)$, $(1,1)$; for this process [Reference McDonald47] constructs explicitly an h function of the form $h(x) =a_1^{x_1} a_2^{x_2}$, where $(a_1,a_2) \in {\mathbb R}^2$ is a point on a curve whose definition is analogous to the definition of the characteristic surface ${\mathcal H}$ in two dimensions.

The work [Reference Miyazawa49] employs the ideas of removing constraints on one of the boundaries and using points on curves associated with the resulting process to study the tail asymptotics of the stationary distribution of a two-dimensional nearest neighbor random walk ${\boldsymbol{L}}$ constrained to remain in ${\mathbb Z}_+^2.$ To study the asymptotic decay rate of ${\boldsymbol{\nu}}(n,k)$ in n for a fixed k, [Reference Miyazawa49] considers the random walk ${\boldsymbol{L}}^{(1)}$, which has the same dynamics as ${\boldsymbol{L}}$ except that it is not constrained on the vertical axis. Associated with this process, [Reference Miyazawa49] defines two curves, whose definitions are similar to the definition of the two-dimensional version of the characteristic surface of the present work (see the definition of ${\mathscr D}_1$ on [Reference Miyazawa49, page 554]) and uses points on and inside these curves to define solutions to an eigenvalue/eigenvector problem associated with the problem (see [Reference Miyazawa49, Theorem 3.1]); for the study of tail asymptotics along the vertical axis, [Reference Miyazawa49] uses the same analysis but this time removing the constraint on the horizontal axis. For further works along this line of research we refer the reader to [Reference Dai and Miyazawa15, Reference Kobayashi and Miyazawa40, Reference Miyazawa50].

The work [Reference Ignatiouk-Robert35] develops an explicit formula for the large deviation local rate function $L(x,v)$ of a general Jackson network, starting from representations of these rates as limits derived in [Reference Atar and Dupuis4, Reference Dupuis and Ellis21]. For this, [Reference Ignatiouk-Robert35] employs “free processes;” these are versions of the original process obtained by removing those constraints from the original process that are not involved in a given direction v at a given point $x \in {\mathbb R}_+^d.$ The proofs in [Reference Ignatiouk-Robert35] use fluid limits for the free process under a change of measure (i.e., a twisted/h-transformed version of the free process); the changes of measures used here correspond to using h-functions of the form $e^{\langle \theta, x \rangle}$ where θ is a point on a characteristic surface (analogous to ${\mathcal H}$ in this work or H in [Reference Dupuis, Sezer and Wang23]) associated with the process being transformed (see [Reference Ignatiouk-Robert35, Section 6]). As an application of its results, [Reference Ignatiouk-Robert35] computes the limit $\lim_{n\rightarrow \infty} \frac{1}{n}\log {\mathbb E}_0[\tau_n]$ by noting from [Reference Parekh and Walrand53] that this limit equals

\begin{equation*} \lim_{n\rightarrow \infty} -\frac{1}{n}\log P_0( \tau_n \lt \tau_0), \end{equation*}

which is the LD decay rate of the probability we have studied in this paper for general stable Jackson networks; [Reference Ignatiouk-Robert35] derives the explicit formula

$\min_{1 \le i \le d} -\log(\rho_i)$ for the above LD rate using the explicit local rate functions developed in the same work and the explicit formulas available for the stationary distribution of the underlying process. Note that the result concerns the initial point x = 0.

The Martin boundary of an unstable process is a characterization of the directions through which the process may diverge to $\infty.$ The idea of using points on characteristic surfaces, and the idea of removing constraints from the process to simplify analysis, appear also in works devoted to identifying Martin boundaries of constrained or stopped processes. An example is [Reference Ignatiouk-Robert and Loree36], which identifies the Martin boundary of two-dimensional random walks in ${\mathbb Z}_+^2$ and which are stopped as soon as they hit the boundary of ${\mathbb Z}_+^2$. This work breaks up its analysis into three cases: 1)the directions $q \in {\mathbb R}_+^2$, where both components of q are nonzero, 2) the directions q such that $q(1) = 0$, and 3) directions such that $q(2) = 0$. For each of these cases, [Reference Ignatiouk-Robert and Loree36] work with what it calls local processes; the local process for the first case is a completely unconstrained random walk, the local process for the second case is a process keeping the horizontal axis (i.e., the vertical boundary is removed) and the third case is the reverse of the last. [Reference Ignatiouk-Robert and Loree36] uses LD analysis of the local processes, harmonic functions of the form

\begin{equation*} h_a(x) = \begin{cases} x_1 e^{\langle a, x\rangle} - {\mathbb E}_x[S_1(\tau) e^{\langle a,x\rangle}1_{\{\tau \lt \infty\}}], &\ \text{if } q(a) = (0,1),\\ x_2 e^{\langle a, x\rangle} - {\mathbb E}_x[S_2(\tau) e^{\langle a,x\rangle}1_{\{\tau \lt \infty\}}], &\ \text{if } q(a) = (1,0),\\ e^{\langle a, x\rangle} - {\mathbb E}_x[e^{\langle a,x\rangle}1_{\{\tau \lt \infty\}}], &\ \text{otherwise.} \end{cases} \end{equation*}

where S is the underlying process, τ is the first hitting time to the boundary of ${\mathbb Z}_+^2$, a is a given point on a surface associated with S (defined analogous to ${\mathcal H}$), q(a) is the mean direction of S under an exponential change of measure defined by a (see [Reference Ignatiouk-Robert and Loree36, page 1108]). In this connection let us also cite [Reference Kurkova and Malyshev42], which uses geometry and complex analysis to identify the Martin boundary of random walks on ${\mathbb Z}^2$, ${\mathbb Z} \times {\mathbb Z}_+$ and ${\mathbb Z}_+^2.$

Let X be the constrained random walk in ${\mathbb Z}_+^2$ with increments $(1,0)$, $(-1,0)$, $(0,1)$, and $(0,-1)$ and let τn be as in (8). A classical problem in computer science going back to [Reference Knuth39, section 2.2.2, exercise 13] is the analysis of the following expectation:

(151)\begin{equation} {\mathbb E} \left[ \max(X_1(\tau_n),X_2(\tau_n)) \right], \end{equation}

i.e., the expected size of the longest queue at the time of buffer overflow. This expectation is computed in [Reference Knuth39] for the case $P(I_k=(1,0)) =P(I_k=(0,1)) = 1/2$, $P(I_k=(-1,0)) =P(I_k=(0,-1)) = 0$. Various versions of this problem were treated in [Reference Flajolet27, Reference Guillotin-Plantard and Schott33, Reference Louchard, Schott, Tolley and Zimmermann45, Reference Maier46, Reference Yao66]. [Reference Maier46] treats a generalization of this problem where the dynamics of the random walk depend on its position; the approach of [Reference Maier46] uses large deviations techniques from [Reference Freidlin and Wentzell31]. [Reference Yao66] treats the approximation of (151) for the case when the increments have a symmetric distribution as follows: $P(I_k =(1,0)) = P(I_k(0,1)) = (1-p)/2$ and $P(I_k =(-1,0)) = P(I_k(0,-1)) = p/2$; furthermore $p \lt 1/2$ is assumed, i.e., the process is assumed unstable. Under these assumptions, [Reference Yao66] develops an approximation for the expectation in (151) as $n\rightarrow \infty.$ The main idea in [Reference Yao66] is the following: under the assumptions of the paper one can ignore both of the constraining boundaries of the process, to prove this the author uses LD bounds on iid Bernoulli sequences (see [Reference Yao66, Lemma 3]). Then an explicit computation for the unconstrained process using elementary techniques gives the desired approximation.

9. Conclusion

In Section 5, we computed ${\mathbb{P}}_y(\tau \lt \infty)$ under the assumption $\mu_i \neq \mu_j$ for $i\neq j.$ One can obtain formulas for ${\mathbb{P}}_y(\tau \lt \infty)$ when this assumption is violated by computing limits of (138) as $\mu_i \rightarrow \mu_j$; this limiting process introduces polynomial terms to the formula. For example, for d = 3 and $\mu_1 = \mu_2 = \mu_3=\mu$ we get

\begin{equation*} {\mathbb{P}}_y(\tau \lt \infty) = \rho^{\bar{y}(1)} \left( \frac{1}{2}c_0^2 (\bar{y}(1))^2 \rho^{y(2)+y(3)} + \rho^{y(3)} \left( \left( \frac{c_0^2}{2} + y(3) c_0^2 \right) \rho^{y(2)} + c_0\right) \bar{y}(1) +1 \right), \end{equation*}

where $c_0 = (\mu - \lambda)/\mu$ and $\bar{y}(1) = y(1) -( y(2) + y(3))$. Similar limits can be computed explicitly for the cases $\mu_1 = \mu_2 \neq \mu_3$, $\mu_1=\mu_3 \neq \mu_2$ and $\mu_1 \neq \mu_2 = \mu_3$. A systematic study of these cases in three and higher dimensions remain for future work.

Recall that the computation of ${\mathbb{P}}_y(\tau \lt \infty)$ (Theorem 5.7) is based on the sequence of regular graphs $G_{d,{\boldsymbol d}}$. We came up with the definitions of these graphs and the solutions of the harmonic systems defined by them through trial and error. Proposition 5.5 (which proves that the definitions (130) provide a solution to the harmonic systems defined by $G_{d,{\boldsymbol d}}$) consists of a verification argument, i.e., it consists of the proof that a candidate/proposed solution is really a solution. This type of argument is common in the solution of differential equations and in stochastic optimal control where one first guesses the solution and then verifies (using the particular structure of the guess) that it is correct. The structure of $G_{d,{\boldsymbol d}}$ and the solution of the harmonic system associated with it depend on the tandem structure of the underlying process since they are directly based on the characteristic equations, conjugacy relations and the C function associated with the tandem case (see (124), (102) and (126)). These are significantly simpler in the tandem case than they are in the general case (see (92), (101) and (96)). Therefore, we do not expect the constructions in Section 5 to easily generalize to arbitrary Jackson networks. Furthermore, at this point we do not have a systematic way of generating harmonic systems and their solutions for arbitrary Jackson networks. This appears to be a challenging problem for future research.

Our convergence analysis (Section 3) and the computations in Section 5 also depend nontrivially on the exit boundary $\partial A_n$ (7). Generalization to other exit boundaries is another important direction that can be explored in the future.

Symbols and notation

  1. (1) tandem walk X, its increments Ik and the map π constraining X to ${\mathbb Z}_+^d$: (2)

  2. (2) constraining boundaries $\partial_j$: (3)

  3. (3) possible increments ${\mathcal V}$ of the tandem walk: (4)

  4. (4) utilization rates ρj and ρ: (6)

  5. (5) domain An: (7), first hitting time τn to the boundary of An: (8), overflow probability pn: (9)

  6. (6) linear map ${\mathcal I}_1$ and affine transformation map Tn: (10)

  7. (7) domain $\Omega_Y$, limit process Y, its increments Jk and the map π 1 constraining Y to $\Omega_Y = {\mathbb Z} \times {\mathbb Z}_+^{d-1}$: (11)

  8. (8) domain B: (12), limit hitting time τ: (13)

  9. (9) Sets Rρ, A and the function g: (14)

  10. (10) $|a|$ where a is a finite set: the cardinality of a (before (17))

  11. (11) Hitting times $\sigma_{j-1,j}$: (24)

  12. (12) process $\bar{X}$: (25)

  13. (13) Hitting time $\bar{\tau}_n$: (27)

  14. (14) The summation function ${\boldsymbol S}$: (29)

  15. (15) Symmetric difference Δ: (40)

  16. (16) Y-superharmonicity: (46)

  17. (17) Function $h_{k,r}$: (47), coefficient γk: (51), function $h_{2,k,r}$: (52)

  18. (18) Coefficient $\gamma_{k-1,k}$: (59), coefficient $\Gamma_k$: (61)

  19. (19) Process S: (62)

  20. (20) Function gn: (75)

  21. (21) Order relation $\preccurlyeq$: (77), set ${\mathscr M}$: (78)

  22. (22) Log-linear functions $[(\beta,\alpha),\cdot]$ (91)

  23. (23) Characteristic polynomials ${\mathbf p}_a$: (92), characteristic equations (93), characteristic surfaces ${\mathcal H}_a$: (94)

  24. (24) Operator Da: (95)

  25. (25) Function C:(96)

  26. (26) $\alpha\{j\}$: (98)

  27. (27) Definition of an L-regular graph: Definition 4.6

  28. (28) Definition of a harmonic system: Definition 4.7

  29. (29) Simple extensions: Definition 4.9 and Definition 4.10

  30. (30) $\partial B$-determined Y-harmonic functions: (122)

  31. (31) Regular graphs $G_{d,{\boldsymbol d}}$: (128)

  32. (32) ${\boldsymbol{c}}^*$, $\alpha^*$, $\beta^*$: (130)

  33. (33) $h^*_d$: (136)

Funding statement

Sections 4, 5 and 7 of the present work are revisions of Sections 5, 6 and subsection 8.3 of the preprint [Reference Sezer64] whose underlying research was supported by the Rbuce-UP COFUND Programme (https://cordis.europa.eu/project/id/246556).

Conflict of interest

The author declares no conflict of interest.

Footnotes

1 The works [Reference Ünlü and Sezer1, Reference Başoğlu Kabran and Sezer5] use complex valued harmonic functions in their approximation formulas.

2 Recall that for two sets D and R, RD denotes the set of functions from D to R. Setting $\alpha \in {\mathbb C}^{\{2,3,...,d\}}$ corresponds to indexing the components of α by $\{2,3,...,d\}$, see (21).

References

Ünlü, K.D. & Sezer, A.D. (2019). Excessive backlog probabilities of two parallel queues. Annals of Operations Research 134.Google Scholar
Anantharam, J., Heidelberger, P. & Tsoucas, P. (1990). Analysis of rare events in continuous time Markov chains via time reversal and fluid approximation. Tech Rep, IBM Research.Google Scholar
Asmussen, S. & Glynn, P. (2007). Stochastic simulation: Algorithms and analysis, Vol. 57. Springer Science & Business Media.10.1007/978-0-387-69033-9CrossRefGoogle Scholar
Atar, R. & Dupuis, P. (1999). Large deviations and queueing networks: methods for rate function identification. Stochastic Processes and Their Applications 84(2): 255296.10.1016/S0304-4149(99)00051-4CrossRefGoogle Scholar
Başoğlu Kabran, F. & Sezer, A.D. (2022). Approximation of the exit probability of a stable Markov modulated constrained random walk. Annals of Operations Research 310: 431475.10.1007/s10479-020-03693-7CrossRefGoogle Scholar
Blanchet, J. (2013). Optimal sampling of overflow paths in Jackson networks. Mathematics of Operations Research 38(4): 698719.10.1287/moor.2013.0586CrossRefGoogle Scholar
Blanchet, J., Glynn, P. & Leder, K. (2008). Efficient simulation of light-tailed sums: an old folk song sung to a faster new tune. Monte Carlo and Quasi-Monte Carlo Methods 2008 227258.Google Scholar
Blanchet, J. & Lam, H. (2012). State-dependent importance sampling for rare-event simulation: An overview and recent advances. Surveys in Operations Research and Management Science 17(1): 3859.10.1016/j.sorms.2011.09.002CrossRefGoogle Scholar
Borovkov, A.A. & Al’fredovich Mogul’skii, A. (2001). Large deviations for markov chains in the positive quadrant. Russian Mathematical Surveys 56(5): 803916.10.1070/RM2001v056n05ABEH000398CrossRefGoogle Scholar
Buijsrogge, A., de Boer, P.-T. & Scheinhardt, W.R.W. (2021). Importance sampling for markovian tandem queues using subsolutions: exploring the possibilities. Simulation 97(12): 849866.10.1177/00375497211041351CrossRefGoogle Scholar
Collingwood, J., Foley, R.D. & McDonald, D.R. (2011). Networks with cascading overloads. Proceedings of the 6th International Conference on Queueing Theory and Network Applications, ACM. pp. 3337.10.1145/2021216.2021221CrossRefGoogle Scholar
Comets, F., Delarue, F. & Schott, R. (2007). Distributed algorithms in an ergodic markovian environment. Random Structures & Algorithms 30(1-2): 131167.10.1002/rsa.20154CrossRefGoogle Scholar
Comets, F., Delarue, F. & Schott, R. (2009). Large deviations analysis for distributed algorithms in an ergodic markovian environment. Applied Mathematics and Optimization 60(3): 341396.10.1007/s00245-009-9079-8CrossRefGoogle Scholar
Crane, M.A. & Iglehart, D.L. (1974). Simulating stable stochastic systems, I: General multiserver queues. Journal of the Association for Computing Machinery 21(1): 103113.10.1145/321796.321805CrossRefGoogle Scholar
Dai, J.G. & Miyazawa, M. et al. (2011). Reflecting Brownian motion in two dimensions: Exact asymptotics for the stationary distribution. Stochastic Systems 1(1): 146208.10.1287/10-SSY022CrossRefGoogle Scholar
De Boer, P.-T. (2006). Analysis of state-independent importance-sampling measures for the two-node tandem queue. ACM Transactions on Modeling and Computer Simulation (TOMACS) 16(3): 225250.10.1145/1147224.1147226CrossRefGoogle Scholar
De Boer, P.-T., Kroese, D.P. & Rubenstein, R.Y. (2004). A fast cross-entropy method for estimating buffer overflows in queueing networks. Management Science 50: 883895.10.1287/mnsc.1030.0139CrossRefGoogle Scholar
De Boer, P.-T. & Nicola, V.F. (2001). Adaptive state-dependent importance sampling simulation of Markovian queueing networks. European Transactions on Telecommunications 13: 303315.10.1002/ett.4460130403CrossRefGoogle Scholar
Dean, T. & Dupuis, P. (2009). Splitting for rare event simulation: A large deviation approach to design and analysis. Stochastic Processes and Their Applications 119(2): 562587.10.1016/j.spa.2008.02.017CrossRefGoogle Scholar
Diestel, R. (fifth ed.) (2017). Graph theory. Graduate Texts in Mathematics, Vol. 173: Springer, Heidelberg; New York.Google Scholar
Dupuis, P. & Ellis, R.S. (1995). The large deviation principle for a general class of queueing systems. i. Transactions of the American Mathematical Society 347(8): 26892751.Google Scholar
Dupuis, P., Leder, K. & Wang, H. (2007). Importance sampling for sums of random variables with regularly varying tails. ACM Trans. Model. Comput. Simul. 17(3): 14.10.1145/1243991.1243995CrossRefGoogle Scholar
Dupuis, P., Sezer, A.D. & Wang, H. (2007). Dynamic importance sampling for queueing networks. Annals of Applied Probability 17(4): 13061346.10.1214/105051607000000122CrossRefGoogle Scholar
Dupuis, P. & Wang, H. (2004). Importance sampling, large deviations and differential games. Stochastics and Stochastic Reports 76(6): 481508.10.1080/10451120410001733845CrossRefGoogle Scholar
Dupuis, P. & Wang, H. (2007). Subsolutions of an Isaacs equation and efficient schemes for importance sampling. Mathematics of Operations Research 32(3): 723.10.1287/moor.1070.0266CrossRefGoogle Scholar
Dupuis, P. & Wang, H. (2009). Importance Sampling for Jackson Networks. Queueing Systems 62: 113157.10.1007/s11134-009-9124-yCrossRefGoogle Scholar
Flajolet, P. (1986). The Evolution of two Stacks in Bounded Space and Random Walks in a triangle, Springer.10.1007/BFb0016257CrossRefGoogle Scholar
Foley, R.D. & McDonald, D.R. et al. (2005). Large deviations of a modified Jackson network: Stability and rough asymptotics. The Annals of Applied Probability 15(1B): 519541.10.1214/105051604000000666CrossRefGoogle Scholar
Foley, R.D. & McDonald, D.R. (2012). Constructing a harmonic function for an irreducible nonnegative matrix with convergence parameter r > 1. Bulletin of the London Mathematical Society bdr115.Google Scholar
Frater, M.R., Lennon, T.M. & Anderson, B.D.O. (1991). Optimally efficient estimation of the statistics of rare events in queueing networks. IEEE Transactions on Automatic Control 36(12): 13951405.10.1109/9.106155CrossRefGoogle Scholar
Freidlin, M.I & Wentzell, A.D. (1998). 2nd edition Random Perturbation of Dynamical systems, Springer-Verlag Telos.10.1007/978-1-4612-0611-8CrossRefGoogle Scholar
Glasserman, P. & Kou, S.-G. (1995). Analysis of an importance sampling estimator for tandem queues. ACM Transactions on Modeling and Computer Simulation 5: 2242.10.1145/203091.203093CrossRefGoogle Scholar
Guillotin-Plantard, N. & Schott, R. (2006). Dynamic Random walks: Theory and applications, Elsevier.10.1016/B978-044452735-6/50038-5CrossRefGoogle Scholar
Ignatiouk, I., Malyshev, V.A. & Scherbakov, V. (1994). Boundary effects in large deviation problems. Russian Mathematical Surveys 49(2): 4199.10.1070/RM1994v049n02ABEH002204CrossRefGoogle Scholar
Ignatiouk-Robert, I. (2000). Large deviations of Jackson networks. Annals of Applied Probability 9621001.Google Scholar
Ignatiouk-Robert, I. & Loree, C. (2010). Martin boundary of a killed random walk on a quadrant. The Annals of Probability 11061142.Google Scholar
Juneja, S. & Nicola, V. (2005). Efficient simulation of buffer overflow probabilities in Jackson networks with feedback. ACM Transactions on Modeling and Computer Simulation 15: 281315.10.1145/1113316.1113317CrossRefGoogle Scholar
Juneja, S. & Shahabuddin, P. (2006). Rare-event simulation techniques: An introduction and recent advances. Handbooks in Operations Research and Management Science 13: 291350.10.1016/S0927-0507(06)13011-XCrossRefGoogle Scholar
Knuth, D.E. (1972). Art of Computer Programming Volume 1: Fundamental algorithms, Addison-Wesley Publishing Company.Google Scholar
Kobayashi, M. & Miyazawa, M. (2013). Revisiting the tail asymptotics of the double QBD process: Refinement and complete solutions for the coordinate and diagonal directions, Matrix-analytic methods in stochastic models, Springer. pp. 145185.Google Scholar
Kroese, D.P. & Nicola, V. (2002). Efficient simulation of Jackson networks. ACM Transactions on Modeling and Computer Simulation 12: 119141.10.1145/566392.566395CrossRefGoogle Scholar
Kurkova, I. & Malyshev, V.A. (1998). Martin boundary and elliptic curves. Markov Process. Related Fields 4(2): 203272.Google Scholar
Lee, J. (2000). Fast simulation for excessive backlogs in tandem networks. Communications for Statistical Applications and Methods 7(2): 499511.Google Scholar
Louchard, G. & Schott, R. (1991). Probabilistic analysis of some distributed algorithms. Random Structures & Algorithms 2: 151186.10.1002/rsa.3240020203CrossRefGoogle Scholar
Louchard, G., Schott, R., Tolley, M. & Zimmermann, P. (1994). Random walks, heat equation and distributed algorithms. Journal of Computational and Applied Mathematics 53(2): 243274.10.1016/0377-0427(94)90048-5CrossRefGoogle Scholar
Maier, R.S. (1991). Colliding stacks: A large deviations analysis. Random Structures & Algorithms 2(4): 379420.10.1002/rsa.3240020404CrossRefGoogle Scholar
McDonald, D.R. (1999). Asymptotics of first passage times for random walk in an orthant. Annals of Applied Probability 110145.Google Scholar
Miretskiy, D., Scheinhardt, W. & Mandjes, M.R.H. (2008). State-dependent importance sampling for a Jackson tandem network.Google Scholar
Miyazawa, M. (2009). Tail decay rates in double QBD processes and related reflected random walks. Mathematics of Operations Research 34(3): 547575.10.1287/moor.1090.0375CrossRefGoogle Scholar
Miyazawa, M. (2011). Light tail asymptotics in multidimensional reflecting processes for queueing networks. Top 19(2): 233299.10.1007/s11750-011-0179-7CrossRefGoogle Scholar
Ney, P. & Nummelin, E. (1987). Markov additive processes i. eigenvalue properties and limit theorems. The Annals of Probability 561592.Google Scholar
Nicola, V. & Zaburnenko, T. (2007). Efficient importance sampling heuristics for the simulation of population overflow in Jackson networks. ACM Transactions on Modeling and Computer Simulation (TOMACS) 17(2): 10.10.1145/1225275.1225281CrossRefGoogle Scholar
Parekh, S. & Walrand, J. (1989). A quick simulation method for excessive backlogs in networks of queues. IEEE Transactions on Automatic Control 34(1): 5466.10.1109/9.8649CrossRefGoogle Scholar
Randhawa, R.S. & Juneja, S. (2004). Combining importance sampling and temporal difference control variates to simulate markov chains. ACM Transactions on Modeling and Computer Simulation 14(1): 130.10.1145/974734.974735CrossRefGoogle Scholar
Revuz, D. (1984). Markov chains, North-Holland.Google Scholar
Ridder, A. (2009). Importance sampling algorithms for first passage time probabilities in the infinite server queue. European Journal of Operational Research 199(1): 176186.10.1016/j.ejor.2008.10.028CrossRefGoogle Scholar
Robert, P. (2003). Stochastic networks and queues, Stochastic modelling and applied probability series, Vol. 52. New York: Springer.Google Scholar
Rubino, G. & Tuffin, B. (2009). Rare event simulation using Monte Carlo methods, John Wiley & Sons.10.1002/9780470745403CrossRefGoogle Scholar
Sezer, A.D. (2005). Dynamic importance sampling for queueing networks, Ph.D. thesis. Brown University Division of Applied Mathematics.Google Scholar
Sezer, A.D. (2007). Asymptotically optimal importance sampling for Jackson networks with a tree topology. arXiv preprint arXiv:0708.3260.Google Scholar
Sezer, A.D. (2007). Asymptotically optimal importance sampling for Jackson networks with a tree topology. Queueing Systems 64(2): 103117.10.1007/s11134-009-9139-4CrossRefGoogle Scholar
Sezer, A.D. (2009). Importance sampling for a Markov modulated queuing network. Stochastic Processes and Their Applications 119 (2): 491517.10.1016/j.spa.2008.02.009CrossRefGoogle Scholar
Sezer, A.D. (2010). Asymptotically optimal importance sampling for Jackson networks with a tree topology. Queueing Systems 64(2): 103117.10.1007/s11134-009-9139-4CrossRefGoogle Scholar
Sezer, A.D. (2015). Exit probabilities and Balayage of constrained random walks. arXiv preprint arXiv:1506.08674.Google Scholar
Sezer, A.D. (2018). Approximation of excessive backlog probabilities of two tandem queues. Journal of Applied Probability 55(3): 968997.10.1017/jpr.2018.60CrossRefGoogle Scholar
Yao, A.C. (1981). An analysis of a memory allocation scheme for implementing stacks. SIAM Journal on Computing 10(2): 398403.10.1137/0210029CrossRefGoogle Scholar
Figure 0

Figure 1. d tandem queues.

Figure 1

Figure 2. Derivation of the limit problem via an affine transformation in two dimensions.

Figure 2

Figure 3. Networks corresponding to p1 and p2 of (110), second is a simple extension of the first.

Figure 3

Figure 4. A $\{2\}$-regular graph and its simple extension to a $\{2,3,4,5\}$-regular graph.

Figure 4

Figure 5. $G_{d,{\boldsymbol d}}$ for $d={\boldsymbol d}=4$.

Figure 5

Figure 6. Level curves and relative error in four dimensions.

Figure 6

Figure 7. The service rates (blue) and the arrival rate (red) for a 14-dimensional tandem Jackson network.

Figure 7

Figure 8. The graph of Wn over $\{x: x(4) + x(14)= 60, x(j) = 0, j \neq 4,14\}$.