Splitting Algorithms for Rare Events of Semimartingale Reflecting Brownian motions

TL;DR

This paper develops a large deviation framework and particle-based simulation algorithms for estimating rare event probabilities in high-dimensional semimartingale reflecting Brownian motions, demonstrating stability and superiority over standard Monte Carlo methods.

Contribution

It introduces a novel large deviation approach and constructs subsolutions for rare event simulation in SRBMs, improving estimation efficiency and stability.

Findings

The probability of the rare event satisfies a large deviation principle.

The proposed particle-based algorithm is stable and theoretically superior to standard Monte Carlo.

The method applies broadly to positive recurrent SRBMs.

Abstract

We study rare event simulations of semimartingale reflecting Brownian motions (SRBMs) in an orthant. The rare event of interest is that a $d$-dimensional positive recurrent SRBM enters the set $B_n = \{z\in\mathbb{R}^d: \sum_{k=1}^d z_k = n\}$ before reaching a small neighborhood of the origin as $n\to\infty$. We show that under a proper scaling and some regularity conditions, the probability of interest satisfies a large deviation principle. We then construct a subsolution to the variational problem for our rare event, and based on this subsolution construct particle based simulation algorithms to estimate the probability of the rare event. It is shown that the proposed algorithm is stable and theoretically superior to standard Monte Carlo for a broad class of positive recurrent SRBMs.

Equations

$$f(t)\sim g(t)$$

$$f(t)=o(g(t))$$

$$f(t)=O(g(t))$$

$$f(t)=\Theta (g(t))$$

$$Z(t)=z+X(t)+\theta t+RY(t),$$

$$\Gamma :{\mathcal{C}}_0([0,\infty );{\mathbb{R}}^d)\to \mathcal{C}([0,\infty );{\mathbb{R}}_{+}^d);\psi \mapsto \varphi .$$

$$Z_n(t)=\frac{Z(nt)}{n},t\ge 0.$$

$$p_n(z_n)=\mathbb{P}(Z_n({\tau }_n)\in B|Z_n(0)=z_n).$$

$$I_{T,z}(\varphi )$$

$${\mathcal{A}\mathcal{C}}_z([0,T];{\mathbb{R}}^d)=\{\psi \in \mathcal{A}\mathcal{C}([0,T];{\mathbb{R}}^d):\psi (0)=z\}.$$

$$\underset{n\to \infty }{\mathrm{limsup}}\frac{1}{n}\log {\mathbb{P}}_{z_n}\left(Z_n\in F\right)\le -\underset{\varphi \in F}{\inf }{I}_{T,z}(\varphi );$$

$$\underset{n\to \infty }{\mathrm{liminf}}\frac{1}{n}\log {\mathbb{P}}_{z_n}\left(Z_n\in G\right)\ge -\underset{\varphi \in G}{\inf }{I}_{T,z}(\varphi ).$$

$$\underset{n\to \infty }{\mathrm{limsup}}\underset{z\in {\mathbb{R}}_{+}^d}{\sup }\frac{1}{n}\log {\mathbb{E}}_{z/n}\left(e^{cn{\tau }_n}\right)<\infty .$$

$$I(z)=\underset{T\ge 0}{\inf }\underset{\psi \in {\mathcal{A}\mathcal{C}}_0([0,\infty );{\mathbb{R}}^d),\Gamma (\psi )(T)=z}{\inf }\frac{1}{2}{\int }_0^T|\dot{\psi }(t)-\theta {|}_D^{2}dt,$$

$$\underset{n\to \infty }{\lim }\frac{1}{n}\log p_n(z_n)=-\underset{x\in B}{\inf }I(x).$$

$$\frac{1}{2}{\int }_0^T|\dot{\psi }(t)-\theta {|}_D^{2}dt=I(z).$$

$$R=\left[\begin{array}{cc}1& r_2\\ r_1& 1\end{array}\right].$$

$${\theta }_1+r_2{\theta }_2^{-}<0,{\theta }_2+r_1{\theta }_1^{-}<0.$$

$${\tilde{a}}^i=⟨a^i,e^i{⟩}_D{e}^i-⟨a^i,\Sigma n^i{⟩}_D\Sigma n^i.$$

$$C_i=\{\alpha e^i+\beta {\tilde{a}}^i:\alpha \ge 0,\beta \ge 0\}.$$

$$I^0(z)$$

$$I^1(z)$$

$$I^2(z)$$

$$\mathcal{I}(w,v)=\underset{\psi \in {\mathcal{A}\mathcal{C}}_w([0,\infty );{\mathbb{R}}^K)}{\inf }\underset{\varphi \in \Gamma (\psi ),{\tau }_{\varphi }(v)<\infty }{\inf }\frac{1}{2}{\int }_0^{{\tau }_{\varphi }(v)}|\dot{\psi }(t)-\theta {|}_D^{2}dt,$$

$${\tau }_{\varphi }(v)=\inf \{t\ge 0:\varphi (t)=v\}.$$

$${\mathcal{I}}_K^{\ast }(w,v)=\underset{\psi \in {\mathcal{A}\mathcal{C}}_w([0,\infty );{\mathbb{R}}^K)}{\inf }\underset{\varphi \in \Gamma (\psi ),\varphi (t)\in F_K,\forall t\in (0,{\tau }_{\varphi }(v)),{\tau }_{\varphi }(v)<\infty }{\inf }\frac{1}{2}{\int }_0^{{\tau }_{\varphi }(v)}|\dot{\psi }(t)-\theta {|}_D^{2}dt.$$

$$B^J$$

$$A^J$$

$${\alpha }^J(w,v)$$

$${\lambda }^J(w,v)$$

$${\mathcal{I}}_K^{\ast }(w,v)$$

$$=|A^{J_{\ast }}\theta {|}_D|A^{J_{\ast }}(v-w){|}_D-⟨A^{J_{\ast }}\theta ,A^{J_{\ast }}(v-w){⟩}_D,$$

$$x_K(t)=w+t{b}_K(w,v),0\le t\le {\alpha }^{J_{\ast }}(w,v),$$

$$b_K(w,v)=\{\begin{array}{ll}{\alpha }^{J_{\ast }}(w,v)(v-w)-R_{L{J}_{\ast }}{\alpha }^{J_{\ast }}(w,v)B^{J_{\ast }}(v-w)+R_{L{J}_{\ast }}{B}^{J_{\ast }}\theta ,& \text{mbox}when{J}_{\ast }\ne \varnothing ,\\ {\alpha }^{J_{\ast }}(w,v)(v-w),& \text{mbox}when{J}_{\ast }=\varnothing .\end{array}$$

$${\mathcal{I}}_{\varnothing }^{\ast }(w,v)=|\theta {|}_D|v-w{|}_D-⟨\theta ,v-w{⟩}_D,$$

Full text

Keywords: Semimartingale Brownian motions; Rare event simulations; Particle splitting algorithms; Large deviation principle; Variational problems.

MSC: 60F10, 60J60, 65C05, 60J85.

1. 1 Introduction

2. 2 Problem formulation

3. 3 Large deviation principle for the rare events

4. 4 Variational problems

5. 4.1 Explicit solution for the VP in $\mathbb{R}^{2}_{+}$

6. 4.2 Subsolution for VP in $\mathbb{R}^{d}_{+}$ when $d\geq 3$

7. 5 Splitting algorithms for SRBMs

8. 5.1 Discrete approximation for SRBMs

9. 5.2 Particle splitting methods

10. 5.3 RESTART splitting

11. 5.4 Estimation of Rare Event Probabilities

12. 6 Numerical Examples

13. 6.1 Splitting algorithm

14. 6.2 RESTART algorithm

15. 7 Discussion

16. 8 Proofs

17. 8.1 Proofs for Section 3

18. 8.2 Proofs for Section 4

19. 8.3 Proofs for Section 5

20. 8.3.1 Proofs of Propositions 5.1 and 5.2

21. 8.3.2 Proofs of Propositions 5.3 and 5.4

1 Introduction

In many scenarios one is interested in the dynamics of an entity that is constrained to a subset of its possible state space by some mechanism. For example, when studying population levels of many interacting species (e.g., chemical compounds or biological populations), a natural constraint is that all population levels are non-negative. Alternatively, one may simply be interested in the dynamics of particles in a container that they cannot escape from. The so-called reflecting diffusions are a group of stochastic processes that have been developed to study the phenomena of constrained stochastic dynamics. Roughly speaking, a reflecting diffusion process behaves as an ordinary diffusion process in the interior of its domain, and it is pushed into the domain whenever the boundary of the domain is hit.

A substantial amount of research on reflecting diffusions has been carried out going back to the seminal works of Skorokhod [24] and Kingman [22] in one dimension. Watanabe [27] and Stroock & Varadhan [25] proved the existence of reflecting diffusions in higher dimensions when the domain is sufficiently smooth. When dealing with a particle constrained to an orthant, one does not have a smooth domain and an important extension is the work of Harrison and Reiman [17] that studied semimartingale reflecting Brownian motions (SRBMs) constrained to orthants. The SRBMs in orthants arise as diffusion approximations for open queueing networks in heavy traffic; cf. [19, 16, 29]. Some important results on the existence, uniqueness and positive recurrence of SRBMs can be found in the survey paper by Williams [31] (see our Section 2). In particular, Hobson & Rogers [20] and Williams [28] established the sufficient and necessary condition for the positive recurrence of SRBMs in quadrants. A fluid-based sufficient condition for the positive recurrence of SRBMs in general dimensions was developed in Dupuis and Williams [13], which was shown to be a necessary condition in three dimension by Bramson, Dai, and Harrison [3].

Since these papers a vast literature on the topic has emerged, and of particular interest is the body of work concerned with the tail probabilities of stationary distributions of SRBMs (cf. [1, 12, 7]). Until now all results on tail probabilities have been for either two dimensional systems or special cases of higher dimensional systems. For example, [23] assumed the property of rotational symmetry for three dimensional SRBMs, while [6] studied how to ensure the stationary distribution of an SRBM has a product form, and how this can be used to identify the asymptotics of its tail probabilities. The work [14] analyzed the properties of the most likely trajectories that give rise to tail events in SRBMs. Given the challenge of developing analytical results for the tail probabilities of SRBMs, our work focuses on developing numerical algorithms to aid in the estimations of these probabilities.

Broadly speaking, there are two types of algorithms for estimating rare event probabilities – importance sampling and particle methods. The general idea of importance sampling is to simulate the system of interest under a new probability measure that makes the event of interest less rare and then adjust the estimator by the Radon-Nikodym derivative between the original measure and the sampling measure. In particle methods, one simulates multiple copies of the system of interest and uses some selection method to favor copies of the system where the rare event of interest is more likely. In the current work we will focus on a particularly simple particle method called the “splitting algorithm”. In this method, one establishes landmarks in the state space and each time a particle crosses a landmark it is replaced with a random number of offspring. The creation of landmarks where splitting events occur is of course the challenging aspect of creating efficient and stable splitting algorithms. In [8], Dean and Dupuis showed how one can use a large deviation principle to ensure that the splitting algorithm is stable and efficiently estimates the rare event of interest. There have been several works (e.g., [8, 9, 2]) that have looked at rare event simulations using particle methods, but to our knowledge, our work is the first to study rare event simulations for SRBMs using any method.

In this work we develop splitting algorithms to aid in the estimations of tail probabilities for SRBMs. We focus on the particular rare event that a positive recurrent SRBM exits the $L^{1}$ cube of radius $n$ before returning to the cube of radius $\epsilon$. We show that our splitting algorithm can stably estimate this probability for positive recurrent SRBMs and furthermore provides an improvement over standard Monte Carlo method. Given that it is not in general possible to exactly simulate an SRBM, we use the Euler method to simulate an approximate process. We establish conditions on the discretization size to ensure our approximating process satisfies the same large deviation principle as the original SRBM. As a result our splitting algorithm may be used for SRBMs or the approximating processes with identical asymptotic properties. Given the bias of using an approximating process we also analyze the relative error between the tail probabilities under the SRBM versus the approximating process.

The remainder of this work is organized as follows. In Section 2, we define the SRBMs, state the basic assumptions, and define the rare event of interest. In Section 3, we show that the rare event probability satisfies a large deviation principle. The variational problems for SRBMs are studied in Section 4, and especially, in Section 4.2, a subsolution to the variational problem in general dimensions is constructed. In Section 5, we introduce the discretizations of SRBMs, apply the splitting method to the discretized process, and state our main results on the performance of the splitting method for estimating rare events. Finally, in Section 6, we present numerical examples.

Throughout the manuscript we shall use the following notation. For a pair of $d$ dimensional vectors $u$ and $v$, we denote their standard inner product by $\langle u,v\rangle=\sum_{i=1}^{d}u_{i}v_{i}$ and the standard Euclidean norm by $|u|=\sqrt{\langle u,u\rangle}$. We write $u\leq v$ if $u_{i}\leq v_{i}$ for $i=1,\ldots,d.$ For a positive definite $d\times d$ matrix $A$ and $d$ dimensional vectors $u$ and $v$, define the associated inner product by $\langle u,v\rangle_{A}=\langle u,Av\rangle$, and the associated norm $|u|_{A}=\sqrt{\langle u,u\rangle_{A}}.$ The trace of $A$ is denoted by $tr(A)$. For a Polish space $S$, denote by $\mathcal{C}([0,\infty);S)$ the space of continuous functions from $[0,\infty)$ to $S$, $\mathcal{AC}([0,\infty);S)$ the space of absolutely continuous functions from $[0,\infty)$ to $S$, and $\mathcal{C}^{2}([0,\infty);S)$ the space of twice differentiable functions from $[0,\infty)$ to $S$. For $f\in\mathcal{AC}([0,\infty);S)$, let $\dot{f}$ denote its gradient vector, and $\ddot{f}$ its Hessian matrix. For a Borel set $B\subset\mathbb{R}^{d}$, and a constant $c\in\mathbb{R}$, $cB=\{cz\in\mathbb{R}^{d}:z\in B\}.$ For a set $K$ of finite number of elements, denote by $\textit{card}(K)$ its cardinality. Throughout the paper we will use the following notation for the asymptotic behavior of positive functions:

[TABLE]

where $C$ and $c$ are positive constants.

2 Problem formulation

Semimartingale refecting Brownian motions (SRBMs) in the nonnegative orthant were first introduced by Harrison and Reiman [17]. Roughly speaking, such a diffusion process behaves like a Brownian motion in the interior of the orthant. Upon reaching the boundary, it is reflected back into the interior in some specified direction, so as to constrain the process in the orthant. The rigorous definition of SRBMs is stated as follows.

Definition 2.1** (SRBM)**

An SRBM associated with data $(\theta,\Sigma,R)$ in $\mathbb{R}^{d}_{+}$ is a continuous $\{\mathcal{F}_{t}\}$-adapted $d$-dimensional process $Z$ defined on a filtered probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq 0},\mathbb{P})$ such that

[TABLE]

where $z\in\mathbb{R}^{d}_{+}$, $X$ is a $d$-dimensional driftless Brownian motion with covariance matrix $\Sigma$ and $X(0)=0$, and $Y$ is a $d$-dimensional stochastic process such that for $i=1,\ldots,d$,

• (i)

$Y_{i}(0)=0;$**

• (ii)

$Y_{i}$* is continuous and non-decreasing;*

• (iii)

$Y_{i}$* can increase only when $Z_{i}$ is on $F_{i}=\{x\in\mathbb{R}^{d}_{+}:x_{i}=0\}$, i.e., $\int_{0}^{\infty}Z_{i}(t)dY_{i}(t)=0.$*

The parameters $\theta,\Sigma$, and $R$ are known as the drift, covariance matrix, and reflection matrix of the SRBM. We assume that $\Sigma$ is invertible, and for ease of notation we let $D=\Sigma^{-1}$. In [26], the authors consider the weak solution of the SRBM equation (2.1), that is the existence of the triplet $(X,Y,Z)$ which satisfies (2.1). It is shown in [26] that there exists an SRBM $Z$ associated with $(\theta,\Sigma,R)$ if and only if $R$ is completely-$\mathcal{S}$, and in this case, $Z$ is unique in law and is a Feller continuous strong Markov process. A matrix is said to be completely-$\mathcal{S}$ if for its every $k\times k$ principle submatrix $G$, there exists a $k$-dimensional vector $v_{G}$ such that $v_{G}\geq 0$ and $Gv_{G}>0.$

An important tool in the study of SRBMs is the Skorokhod problem (SP). Let $\mathcal{C}_{0}([0,\infty);\mathbb{R}^{d})$ $=\{f\in\mathcal{C}([0,\infty);\mathbb{R}^{d}):f(0)\in\mathbb{R}^{d}_{+}\}.$

Definition 2.2** (SP)**

Let $\psi\in\mathcal{C}_{0}([0,\infty);\mathbb{R}^{d})$. A pair $(\phi,\eta)\in\mathcal{C}([0,\infty);\mathbb{R}^{d}_{+})\times\mathcal{C}([0,\infty);\mathbb{R}^{d}_{+})$ is called a solution of the SP for $\psi$ associated with $R$ if the following hold.

• (i)

For all $t\geq 0,\phi(t)=\psi(t)+R\eta(t)\geq 0$.

• (ii)

$\eta$* satisfies the following: (a) $\eta_{i}(0)=0$, (b) $\eta_{i}(\cdot)$ is nondecreasing, and (c) $\eta_{i}(\cdot)$ increases only when $\phi_{i}(\cdot)=0$, that is, $\int_{0}^{\infty}\phi_{i}(t)d\eta_{i}(t)=0.$*

Theorem 2.2 in [32] guarantees that there is a solution to the SP for $\psi$ if and only if $R$ is completely-$\mathcal{S}$. Define

[TABLE]

If $\Gamma(\psi)$ is unique for each $\psi\in\mathcal{C}_{0}([0,\infty),\mathbb{R}^{d})$, then $\Gamma$ is called the Skorokhod map (SM) associated with the matrix $R$, and we will call $\Gamma(\psi)$ the $R$-regulation for $\psi$. Throughout this work we will make the assumption below to insure that our SRBM satisfies a large deviation principle.

Assumption 2.1

The SM $\Gamma$ is well defined on $\mathcal{C}_{0}([0,\infty);\mathbb{R}^{d})$, and is Lipschitz continuous with respect to the topology of uniform convergence on compact sets.

Note that Assumption 2.1 implies that $R$ is completely-$\mathcal{S}$. The work [11] develops conditions that guarantee the Lipschitz continuity of $\Gamma$.

We require that the SRBM $Z$ is positive recurrent such that the expected time that $Z$ reaches any closed set with positive Lebesgue measure will be finite.

Definition 2.3** (Positive recurrence)**

The SRBM $Z$ is positive recurrent if for every closed set $F$ in $\mathbb{R}^{d}_{+}$ with positive Lebesgue measure, $\mathbb{E}_{z}[\tau_{F}]<\infty$, where $\tau_{F}=\inf\{t\geq 0:Z(t)\in F\}$.

Based on Theorem 2.6 of [13], we make the following assumption which provides a sufficient condition for the positive recurrence of our SRBM.

Assumption 2.2

Let $\psi\in\mathcal{C}_{0}([0,\infty);\mathbb{R}^{d})$ be of the form $\psi(t)=z+\theta t$, $t\geq 0$. Then $\lim_{t\to\infty}\Gamma(\psi)(t)=0$.

For a positive recurrent SRBM $Z$, we introduce a large deviation scaling with the scaling parameter $n\in\mathbb{N}$ as follows: For $n\in\mathbb{N},$ define

[TABLE]

Let $\epsilon\in(0,1)$. Consider two sets $A_{n}=\{z\in\mathbb{R}^{d}_{+}:\sum_{k=1}^{d}z_{k}\leq\epsilon/n\}$ and $B=\{z\in\mathbb{R}^{d}_{+}:\sum_{k=1}^{d}z_{k}\geq 1\}$. We are interested in the probability that the process $Z_{n}$ first enters $B$ before reaching $A_{n}$ with a starting point $z_{n}\notin A_{n}\cup B$. Letting $\tau_{n}=\inf\{t\geq 0:Z_{n}(t)\in A_{n}\cup B\}$, the probability of interest can be defined as

[TABLE]

The goal of this work is to build algorithms that can stably estimate the probability $p_{n}(z_{n})$ and provide an improvement over standard Monte Carlo. Given the positive recurrence assumption on the SRBM $Z$ it is intuitively obvious that $p_{n}$ will go to 0 as $n\to\infty$. Furthermore, the exponential decay rate of the probability will clearly be related to the large deviation properties of the sequence $\{Z_{n}\}$. One minor point that we will have to deal with is that the sets $\{A_{n}\}_{n\geq 1}$ are dependent on the scaling factor $n$. Ideally we would like to replace the sequence of sets $\{A_{n}\}_{n\geq 1}$ with the set $\{0\}$; however for positive recurrent SRBMs it is possible that the expected time to hit the point [math] is infinite. In the next section, we develop the necessary large deviation theory to characterize (in terms of a variational problem) the exponential decay rate of $p_{n}(z_{n})$.

3 Large deviation principle for the rare events

For $T\geq 0$, using Schilder’s theorem and the contraction principle, $\{Z_{n}(t);t\in[0,T]\}$ with initial value $Z_{n}(0)=z$ satisfies the large deviation principle (LDP) in $\mathcal{C}([0,T];\mathbb{R}^{d}_{+})$ with the good rate function

[TABLE]

where $\Gamma_{T}(\psi)$ is the $R$-regulation for $\{\psi(t);t\in[0,T]\}$, and

[TABLE]

In particular we have the following result for any $T>0$. Its proof is a simple application of the contraction principle and is thus omitted.

Lemma 3.1

Assume that $z_{n}\to z\in\mathbb{R}^{d}_{+}$ and let $I_{T,z}$ be as in (3.1). Then

• (i)

for any closed $F\subset C([0,T];\mathbb{R}^{d}_{+}),$

[TABLE]

• (ii)

for any open $G\subset C([0,T];\mathbb{R}^{d}_{+}),$

[TABLE]

However, the probability $p_{n}(z_{n})$ is time independent and we need to extend the previous result to handle infinite time horizon events. A standard technique is to use the stability properties of the stochastic process to establish an upper bound on the size of the stopping time $\tau_{n}$. The result is present in the following proposition.

Proposition 3.1

There exists $c\in(0,\infty)$ such that

[TABLE]

where $\mathbb{E}_{z/n}$ denotes the expectation conditioning on $Z_{n}(0)=z/n$.

We now consider an infinite horizon variational problem which minimizes over all paths that travel from 0 to a given point $z$:

[TABLE]

where $\Gamma(\psi)(T)=z$ if there is a $\phi\in\Gamma(\psi)$ such that $\phi(T)=z$. With Proposition 3.1 the finite horizon LDP can be extended to infinite time horizon events.

Theorem 3.1

Assume $z_{n}=z/n\notin A_{n}\cup B$. Then

[TABLE]

4 Variational problems

This section focuses on the study of the variational problem (VP) defined in (3.2). A triplet $(\psi,\eta,\phi)$ is called a solution to the VP (3.2) if $(\phi,\eta)$ is a solution of the SP for $\psi$, $\phi(T)=z$ for some $T\geq 0$, and

[TABLE]

The solution $(\psi,\eta,\phi)$ is also referred to as an optimal triple, and $\psi$ and $\phi$ are referred to as optimal path and optimal reflected path, respectively. It is known from [1] that the VP in $\mathbb{R}_{+}^{2}$ admits an explicit analytical solution. We summarize this result in Section 4.1. For the VP in $\mathbb{R}^{3}_{+}$, a special case on rotationally symmetric VP is studied in [23], and general paths properties are studied in [14]. In [15], the value of the VP in $\mathbb{R}^{d}_{+}$ is characterized as a solution to a partial differential equation. However, the complete analysis of the VP in three or higher dimensional space is still open. In Section 4.2, we apply the path properties from [14] to construct a subsolution for the VP in $\mathbb{R}^{d}_{+}$ when $d\geq 3$.

4.1 Explicit solution for the VP in $\mathbb{R}^{2}_{+}$

The VP in $\mathbb{R}^{2}_{+}$ is analytically solved in [1], in which the solution is characterized by defining “cones of boundary influence”. In this section, we summarize the result in [1] and adopt its notation. Without loss of generality, we let

[TABLE]

The parameters $\theta$ and $R$ are required to satisfy the following assumption, which is sufficient and necessary for the positive recurrence of $Z$.

Assumption 4.1

The drift $\theta$ and reflection matrix $R$ satify that

[TABLE]

For $i=1,2$, we choose a vector $p^{i}$ which is orthogonal (under the usual Euclidean inner product) to the $i$th column of $R$, and $|\Sigma p^{i}|_{D}=1.$ We define the face $F_{i}=\{x\in\mathbb{R}^{2}_{+}:x_{i}=0\},$ and the associated exit velocity $a^{i}=\theta-2(\theta^{\prime}p^{i})\Sigma p^{i}.$ The face $F_{i}$ is said to be reflective if the $i$th component of $a^{i}$ is negative, i.e., $a^{i}_{i}<0$. Denote $e^{1}=[0\;\,1]^{T}$, $e^{2}=[1\;\,0]^{T}$. Let $n^{i}$ be a vector that is normal to $F_{i}$, pointing to the interior of the state space, and normalized such that $|\Sigma n^{i}|_{D}=1.$ We introduce the key definition from [1] as follows:

Definition 4.1** (Cone of influence)**

• (i)

When $F_{i}$ is not reflective, the cone of influence $C_{i}$ is defined to be an empty set. In this case, the face $F_{i}$ has no boundary influence on solutions to the VP for any $z\in\mathbb{R}^{2}_{+}$.

• (ii)

When $F_{i}$ is reflective, the cone of influence $C_{i}$ is defined to be the cone generated by $e^{i}$ and $\tilde{a}^{i}$, where

[TABLE]

That is

[TABLE]

For $z\in\mathbb{R}^{2}_{+},$ the cone of influence $C_{i}$ is said to be active if $z\in C_{i}$. We define three cost functions:

[TABLE]

where $\tilde{a}^{0}(z)=z|\theta|_{D}/|z|_{D}$. With Assumption 4.1, the authors in [1] showed that an optimal path can always be found among these three paths. The search for such path depends on the location of $z$ and its active cone of influence. The following is the main theorem from [1].

Theorem 4.1** (Theorem 3.1 of [1])**

Under Assumptions 2.1 and 4.1, the following hold.

• (i)

If $z\notin C_{1}\cup C_{2},$ the optimal value of the VP is given by $I^{0}(z)$.

• (ii)

If $z\in C_{1}\backslash C_{2},$ the optimal value of the VP is given by $I^{1}(z)$.

• (iii)

If $z\in C_{2}\backslash C_{1},$ the optimal value of the VP is given by $I^{2}(z)$.

• (iv)

If $z\in C_{1}\cap C_{2},$ the optimal value of the VP is given by $\min\{I^{1}(z),I^{2}(z)\}$.

4.2 Subsolution for VP in $\mathbb{R}^{d}_{+}$ when $d\geq 3$

In this section we construct a sub-solution to the VP (3.2) in ${\mathbb{R}}^{d}_{+}$, $d\geq 3$. We first define $L=\{1,2,...,d\}$. For $J,K\subset L$, define $R_{JK}$ to be a $\textit{card}(J)\times\textit{card}(K)$ matrix by deleting rows and columns with indices in $L\backslash J$ and $L\backslash K$, respectively. Let $R_{j}$ denote the $j$th column of $R$. Finally, define the face associated with $K\subset L$ to be $F_{K}=\{x\in\mathbb{R}^{d}_{+}:x_{i}=0,i\in K,x_{j}>0,j\in L\backslash K\}.$

A $d\times d$ matrix is called a nonsingular $\mathcal{M}$-matrix if it can be expressed as $\mu\mathbb{I}_{d}-A$, where $\mathbb{I}_{d}$ is a $d\times d$ identity matrix, $A$ is a $d\times d$ matrix with nonnegative entries, and $\mu$ is greater than the spectral radius of $A$. We also note that $R^{-1}>0$, which in particular implies that $R^{-1}\theta<0$ under Assumption 4.2. In this section, we make the following assumption on $R$

Assumption 4.2

The reflection matrix $R$ is a nonsingular $\mathcal{M}-$matrix, and $\theta<0.$

This implies that the SRBM $Z$ is positive recurrent (see [18]), and also implies Assumptions 2.1 and 2.2 (see [5, 18]).

For any $w,v\in\mathbb{R}^{d}_{+},w\neq v,$ define the globally optimal cost from $w$ to $v$ as follows:

[TABLE]

where

[TABLE]

A path $\psi^{*}$ which leads to the globally optimal cost is called a globally optimal path from $w$ to $v$. We next define the locally optimal cost between points $w$ and $v$, where we constrain the path to lie on a face $F_{K}$. In particular we will consider the local cost between points satisfying the following condition with respect to $K\subset L$.

Condition 4.1

$w\neq v;w_{i}=v_{i}=0,\forall i\in K$*, and for each $j\notin K$ $w_{j}\neq 0$ or $v_{j}\neq 0$. *

For the pair $v$ and $w$ satisfying Condition 4.1 with respect to $K\subset L$, we define the locally optimal cost from $w$ to $v$ with path constrained in $F_{K}$ as follows:

[TABLE]

We now introduce the definition of a sub-solution to the VP (3.2).

Definition 4.2

A function $\bar{T}:\mathbb{R}^{d}_{+}\rightarrow\mathbb{R}$ is called a sub-solution to the VP (3.2) if

• (i)

$\bar{T}$* is continuous;*

• (ii)

$\bar{T}(x)\leq 0$* for all $x\in B$;*

• (iii)

$\bar{T}(w)-\bar{T}(v)\leq\mathcal{I}(w,v)$* for all $w,v\in\mathbb{R}^{d}_{+}\setminus(\{0\}\cup B).$*

The following lemma presents the main idea of constructing a sub-solution.

Lemma 4.1

Let $\bar{T}$ be a function that satisfies conditions (i) and (ii) in Definition 4.2. If in addition, for all $K\subseteq L$, $\bar{T}(w)-\bar{T}(v)\leq\mathcal{I}^{*}_{K}(w,v)$ for all $w,v$ satisfying Condition 4.1 with respect to $K$, then it is a sub-solution to the VP (3.2).

In [14], the authors have investigated a local variational problem which is very similar to the VP (4.4). Specifically, two types of locally optimal paths were investigated: interior escape paths and single segment boundary escape paths. An interior escape path is the one for which no reflection is applied. Let $K$ and $K^{\prime}$ be two subsets of $L$ such that $K\subset K^{\prime}$ and $0<\textit{card}(K)<\textit{card}(K^{\prime})\leq d$. For $w\in F_{K^{\prime}}$, $v\in F_{K}$, an optimal single segment boundary escape path is a path constrained in $F_{K}$ which has the cost equal to $\mathcal{I}^{*}_{K}$. The result in [14] can be applied to the VP (4.4) directly as the proofs are still valid.

In order to state the results in [14] on $\mathcal{I}^{*}_{K}$, we need some further notation. Given $K\subseteq L$, and two points $v,w\in\mathbb{R}^{d}_{+}$ satisfying Condition 4.1 with respect to $K$, for each $J\subseteq K$, define

[TABLE]

Theorem 4.2

Given $K\subseteq L$, and two points $v,w\in\mathbb{R}^{d}_{+}$ which satisfy Condition 4.1 with respect to $K$, there exists a unique $J_{*}\subseteq K$ such that $\langle A^{J_{*}}R_{j},\alpha^{J_{*}}(w,v)(v-w)-\theta\rangle_{D}\leq 0,\ j\in K\backslash J_{*}$, and when $J_{*}\neq\emptyset$, each component of ${\lambda}^{J_{*}}(w,v)$ is positive. Furthermore,

[TABLE]

and an optimal path is

[TABLE]

where

[TABLE]

From Theorem 4.2, we can solve the VP (4.4) by identifying the set $J_{*}$ for the pair of points $v$ and $w$. In particular, when $|K|=0$, the problem (4.4) only concerns paths in the interior of $\mathbb{R}^{d}_{+}$.

Lemma 4.2

Given two points $v,w\in\mathbb{R}^{d}_{+}$ such that at least one of them is in the interior of $\mathbb{R}^{d}_{+}$, we have

[TABLE]

and the optimal path is

[TABLE]

where $b_{\emptyset}(w,v)=\frac{|\theta|_{D}}{|v-w|_{D}}(v-w)$.

When $|K|\geq 1$, we can find $J_{*}$ by testing all subsets of $K$ through the conditions from Theorem 4.2.

In the rest of section, we construct a sub-solution for the VP (4.4), which by Lemma 4.1 is sufficient for establishing a subsolution to (3.2). We start with $T(\cdot)$ being a norm on $\mathbb{R}^{d}$. It is required that $T(\cdot)$ is monotone on $\mathbb{R}^{d}_{+}$, that is if $v,w\in\mathbb{R}^{d}_{+}$ and $v\leq w$, then $T(v)\leq T(w)$. (In Section 5.2, we let $T(x)=\sum_{i=1}^{d}|x_{i}|$ for $x\in\mathbb{R}^{d}$.) We will then rescale this norm in an appropriate fashion by the functions $\mathcal{I}_{K}^{*}(\cdot,\cdot),K\subset L,$ and shift it to allow for the largest possible value at the origin.

Consider a nonempty subset $K\subset L$, and two points $v,w$ satisfying Condition 4.1. We note that, from Theorem 4.2, if $w^{\prime},v^{\prime}\in\mathbb{R}^{d}_{+}$ also satisfy Condition 4.1, and $v^{\prime}-w^{\prime}=\alpha(v-w)$ for some $\alpha\in(0,\infty)$, then $\mathcal{I}^{*}_{K}(w^{\prime},v^{\prime})=\alpha\mathcal{I}^{*}_{K}(w,v)$. Therefore, it is proper to define $\mathcal{I}^{*}_{K}(v-w)=\mathcal{I}^{*}_{K}(w,v)$ for any $v,w\in\mathbb{R}^{d}_{+}$ that satisfy Condition 4.1. One note that $\mathcal{I}^{*}_{K}(\cdot)$ is defined on $\mathbb{R}^{d}.$

Lemma 4.3

Given $K\subseteq L$, and two points $v,w\in\mathbb{R}^{d}_{+}$ which satisfy Condition 4.1 with respect to $K$, we have

• (i)

If $\mathcal{I}^{*}_{K}(w,v)=0$, then $J_{*}=K$, $v_{i}-w_{i}<0,\forall i\in L\backslash K$, and $T(v)-T(w)\leq 0$.

• (ii)

If $\mathcal{I}^{*}_{K}(v)=\mathcal{I}^{*}_{K}(w)=0$, then $v=cw$ for some $c>0$.

For a nonempty subset $K\subset L$, define a direction set $\mathbb{D}_{K}\subset\mathbb{R}^{d}$ such that a vector $v\in\mathbb{D}_{K}$ if and only if $v$ satisfies the following conditions:

• (1)

$v_{i}=0,\forall i\in K$,

• (2)

$\sum_{i=1}^{d}|v_{i}|=1$,

• (3)

there exists $j\in L\backslash K$ such that $v_{j}\geq 0$.

Notice that $\mathbb{D}_{K}$ is a closed set and thus compact due to the continuity of the norm $|\cdot|$. By the definition of $\mathcal{I}^{*}_{K}(v)$, it is easy to see that $\mathcal{I}^{*}_{K}(v)$ is continuous. We now show that $\mathcal{I}^{*}_{K}$ is bounded from below on the set $\mathbb{D}_{K}$.

Proposition 4.1

There exists $\varepsilon>0$ such that for all $K\subset L$ and $v\in\mathbb{D}_{K}$, $\mathcal{I}^{*}_{K}(v)>\varepsilon$.

The above proposition allows us to consider the maximum of the ratio $T$ to $\mathcal{I}^{*}_{K}$ over the set $\mathbb{D}_{K}$. Let

[TABLE]

and define our subsolution to be

[TABLE]

Theorem 4.3

$\bar{T}$* is a subsolution to the VP (3.2).*

5 Splitting algorithms for SRBMs

Particle splitting methods and importance sampling methods are the two most widely used methodologies to obtain numerical estimations of probabilities of rare events. Oftentimes, particle splitting methods are used to simulate a class of rare events called first entrance probabilities. The main idea of particle splitting is to partition the state space of a process into a series of nested subsets. Then the rare event of interest can be considered as a nested sequence of events. When a given subset is entered by a simulated particle for the first time, a number of offspring will be generated at the entry point according to the assigned splitting mechanism. All the offspring will follow the same law of the original process independently. More refined versions of splitting algorithms such as RESTART have been introduced in the last decades (see [9] for details).

In Section 5.1, we construct the Euler discretization of the SRBM, and show that the discretized SRBM and the SRBM are exponentially equivalent. In Sections 5.2 and 5.3, we use the solution/sub-solution developed in Section 4 to develop algorithms to estimate the probability for our rare event of interest. Finally in Section 5.4 we analyze the performance of our estimation algorithms.

5.1 Discrete approximation for SRBMs

In this subsection, we construct a discrete approximation to the SRBM $Z$. Recall that the $d$-dimensional SRBM $Z$ in $\mathbb{R}^{d}_{+}$ is given by

[TABLE]

where $X$ is a $d$-dimensional driftless Brownian motion with covariance matrix $\Sigma$, and $Y$ is the corresponding reflecting process. Fix $\Delta>0$. Consider the discretized process $X^{\Delta}$ and $\tilde{X}^{\Delta}$ defined by $X^{\Delta}(t)=X((k-1)\Delta)$, and $\tilde{X}^{\Delta}(t)=X^{\Delta}(t)+\theta(k-1)\Delta$ for $t\in[(k-1)\Delta,k\Delta)$, and $k\in\mathbb{N}$. Define $(Z^{\Delta},Y^{\Delta})$ as the solution of the Skorokhod problem for $\tilde{X}^{\Delta}$ with the initial $z$. Then we have

[TABLE]

We solve the Skorokhod problem for $\tilde{X}^{\Delta}$ by solving a sequence of Linear Complementarity Problems as outlined in [30].

Introducing the large deviation scaling to $Z^{\Delta}$, we define

[TABLE]

Let $z_{n}=z/n$. The first entrance probability associated with $Z^{\Delta}$ is defined to be

[TABLE]

where $\tau^{\Delta}_{n}=\inf\{t\geq 0:Z^{\Delta}_{n}(t)\in A_{n}\cup B\}.$ In the following, the first proposition establishes the exponential equivalence between $Z_{n}$ and $Z^{\Delta}_{n}$, and the second one characterize the relative error between $p^{\Delta}_{n}(z_{n})$ and $p_{n}(z_{n})$.

Proposition 5.1

Under Assumption 4.2, if $\Delta\equiv\Delta(n)\to 0$ as $n\to\infty$, then for any $T\geq 0$ and $\varepsilon>0,$

[TABLE]

and consequently, for $z_{n}={z}/{n}\notin A_{n}\cup B$,

[TABLE]

Proposition 5.2

Let $\eta\in(0,1/2)$ and $\Delta\equiv\Delta(n)$ satisfy $\lim\limits_{n\rightarrow\infty}{n\Delta}^{\eta}=0$, assume that $z_{n}=z/n\notin A_{n}\cup B$, and denote $z_{n}=(z_{n,1},\ldots,z_{n,d})$. Then under Assumption 4.2,

[TABLE]

where $\mathcal{C}_{n}=\{x\in\mathbb{R}^{d}_{+}:\sum_{k=1}^{d}x_{k}=\sum_{k=1}^{d}z_{n,k}\}.$

Ideally Proposition 5.2 would be a statement of the form

[TABLE]

but unfortunately we were not able to prove such a statement. Since $p_{n}(x)$ will have the same exponential decay rate for all $x\in\mathcal{C}_{n}$ we see that considering a maximum over the set $\mathcal{C}_{n}$ introduces at most a subexponential effect.

5.2 Particle splitting methods

Following [8], we first briefly describe the splitting algorithm to estimate the probabilities $\{p^{\Delta}_{n}(z_{n})\}_{n\geq 1}$. Our state space $\mathbb{R}^{d}_{+}$ is partitioned according to a nested collection of sets $B=C_{0}^{n}\subset C_{1}^{n}\subset\cdots\subset C_{K}^{n}\cdots$. These sets are often constructed as the level sets of a particular function $V$, which is called the importance function in [8], and will be defined in terms of the solution/sub-solution of the VP (3.2).

Given the sequence of nested sets $\{C^{n}_{j}\}_{j\geq 0}$ and a splitting number $r$, the most simple splitting algorithm works in the following way. Define $l^{n}(x)=\min\{j\geq 0:x\in C^{n}_{j}\}$. The algorithm generates a time inhomogeneous branching process with generations indexed by $i\in\{0,1,\ldots,l^{n}(z_{n})\}.$

• •

For generation [math], a single particle starts at $z_{n}$, and it evolves according to the law of $Z_{n}^{\Delta}$ with initial condition $Z_{n}^{\Delta}(0)=z_{n}$.

• •

Let $N^{n}_{i-1}$ denote the number of particles in generation $i-1$, and for $j=1,\ldots,N^{n}_{i-1}$, let $X^{n}_{i-1,j}$ denote the position of the $j$-th particles in generation $i-1$. The $j$-th particle behaves according to the law of $Z_{n}^{\Delta}$ with initial condition $Z_{n}^{\Delta}(0)=X^{n}_{i-1,j}$, and it stops moving when reaching either the set $A_{n}$ or the next level set $C^{n}_{l^{n}(X^{n}_{i-1,j})-1}=C^{n}_{l^{n}(z_{n})-i}.$

• •

If the $j$-th particle in generation $i-1$ reaches $A_{n}$ first, do nothing.

• •

If the $j$-th particle in generation $i-1$ reaches $C^{n}_{l^{n}(z_{n})-i}$ first, denoting by $X^{n}_{i,j}$ its location of entrance to $C^{n}_{l^{n}(z_{n})-i}$, generate $r$ number of new particles in generation $i$ with position $X^{n}_{i,j}$.

Once all the generations have been calculated, define the estimator

[TABLE]

Now consider a sub-solution $\bar{T}$. If the problem is two dimensional, we let $\bar{T}(v)=-I(v)+\inf_{x\in B}I(x)$, where $I(v)$ is the explicit solution to the VP discussed in Section 4.1. It is easy to check that this is indeed a sub-solution, and $\bar{T}(0)=\inf_{x\in B}I(x)$. If the problem is three or higher dimensional, we use the sub-solution constructed in (4.5) of Section 4.2, where the norm $T$ on $\mathbb{R}^{d}$ is defined to be $T(x)=\sum_{i=1}^{d}|x_{i}|$. Define the importance function $V(z)=\delta\bar{T}(z)/\log(r).$ The level sets of the importance function $V$ are defined to be

[TABLE]

Now define the sets $C^{n}_{j}$ to be the level sets of $V$. More precisely, for some $\delta>0$, let

[TABLE]

Here $\delta$ is referred to as the level of the splitting algorithm. In the proposition below we assume that $z_{n}=z/n\notin A_{n}\cup B$.

Proposition 5.3

The splitting algorithm for $Z^{\Delta}$ is stable, i.e., the total number of particles ever used grows subexponentially as $n\to\infty$. Letting $\{s_{n}^{(k)}\}_{k\in\mathbb{N}}$ be an i.i.d sequence of copies of $s_{n}$, then for each $N\in\mathbb{N},$

[TABLE]

Furthermore, the variance of $s_{n}$ can be measured using the following rate of decay of the second moment.

[TABLE]

Proposition 5.3 is a combination of Proposition 7 and Theorem 8 from [8]. It should be noted that [8] assumes the rare event of interest is governed by a large deviations principle with a specific form that is different from what we consider in the current paper. However, the proofs of Proposition 7 and Theorem 8 do not depend on that assumption and still apply to the current setting.

5.3 RESTART splitting

The particle splitting algorithm presented in Section 5.2 is very efficient with regard to the mean and second moment of the number of particles generated if a sub-solution to the associated variational problem can be obtained. However, it can still require a lot of computational effort. The source of inefficiency comes from the fact that most of the particles generated will not end up in the rare event set (in our case, set $B$) [9] . For those particles generated near set $B$ by splitting, most of them will still reach $A_{n}$ before $B$, and it takes much time to simulate those trajectories.

Here, we briefly introduce the RESTART splitting algorithm we implemented specifically for our problem. More details about the RESTART algorithm for rare event simulations can be found in [9]. We first define a sequence of importance functions $\{V^{n}\}$ as follows.

[TABLE]

where $\Delta=\log(r)$ with $r$ being a fixed positive integer, $z_{n}$ is the starting point, and $U(\cdot)$ is a subsolution satisfying the conditions in definition 4.2. For each $n$, we define a sequence of nested sets $C^{n}_{0}\supset C^{n}_{1}\supset\cdots\supset C^{n}_{J_{n}}\supset C^{n}_{J_{n}+1}$, where $C^{n}_{0}=\mathbb{R}^{d}_{+}$, and $C^{n}_{J_{n}+1}=\emptyset$, such that $y\in C^{n}_{j}\setminus C^{n}_{j+1}$ if and only if $V^{n}(y)=j\Delta$. For simplicity, we denote $C^{n}_{j}\setminus C^{n}_{j+1}$ by $D_{j}$.

The algorithm starts from an initial particle $y_{1}$ located at $z_{n}$. Each particle will be assigned a killing threshold $l$ when it is generated. The killing threshold for the initial particle is set to be [math]. A particle is killed whenever it reaches any $D_{j}$ such that $j<l$. A killed particle will not be simulated further. During the simulation, whenever a particle moves from $D_{j}$ to $D_{k}$, $k>j$, in one step, a total number of $r^{k-j}-1$ (excluding the original particle) offspring will be generated. For each integer $j<\alpha\leq k$, $(r-1)r^{\alpha-j-1}$ of the offspring will be assigned the killing threshold $\alpha$. A particle is said to be stopped whenever it reaches set $A_{n}$ or $B$. Notice that splitting can happen at the moment of a particle being stopped. In such case, all the offspring will also be stopped. A stopped particle will not be simulated further. The algorithm terminates when all the particles are stopped or killed.

Let $\mathbb{I}$ be the index of those particles which reach set $B$. For a single trial of simulation, $p_{n}(z_{n})$ is then approximated by $r_{n}=\sum_{i\in I}\exp(V^{n}(y_{i}(\tau)))$, where $y_{i}(\tau)$ is the position of $y_{i}$ when it is stopped. The following proposition summarizes the main results on the RESTART algorithm for $Z^{\Delta}_{n}.$

Proposition 5.4

• (i)

The estimator constructed from the RESTART algorithm is unbiased.

• (ii)

Suppose that $z_{n}\rightarrow 0$ with each $z_{n}\notin(A_{n}\cup B)$. Then

[TABLE]

• (iii)

Let $w^{n}$ denotes the sum of the lifetimes of all the particles simulated in a single trial. Then

[TABLE]

This proposition is a combination of results from [9]. The RESTART algorithm does not necessarily have a faster decay rate of its second moment than standard splitting. However, due to the lower computational effort per replication it often has superior performance to standard splitting. In particular, there might be a slight increase in variance of the estimator due to the use of RESTART instead of splitting, but this increase is often made up for by not simulating particles that are very unlikely to hit the target set.

5.4 Estimation of Rare Event Probabilities

Our goal is to estimate the rare event probability $p_{n}(z_{n})$, and we do this with the estimators $s_{n}$ or $r_{n}$. We will focus our discussion on $s_{n}$, but a similar analysis can be carried out for $r_{n}$. A standard metric for the quality of the estimator is given by the relative expected mean square error

[TABLE]

The previous display shows that $rMSE(s_{n})$ can be decomposed into a relative variance term and a relative bias term. The relative variance term can be rewritten as

[TABLE]

From Proposition 5.2 we know that the second term in the previous display grows at most at subexponential rate. From Theorem 3.1 and Proposition 5.3 we know that

[TABLE]

As mentioned earlier for the 2-d problem we have that $\bar{T}(0)=\inf_{x\in B}I(x)$ and thus the relative variance will grow subexponentially in $n$. An unbiased estimator whose relative variance grows subexponentially is said to be efficient, i.e., in two dimensions $s_{n}$ is an efficient estimator for $p_{n}^{\Delta}(z_{n})$. For three and higher dimensions we have that

[TABLE]

and therefore the relative variance term grows with exponential rate

[TABLE]

which is guaranteed to be non-negative due to the subsolution property.

As a point of comparison suppose one were to use standard Monte Carlo to estimate $p_{n}(z_{n})$ and denote the estimator by $mc_{n}(z_{n})$. Note that the bias terms would remain the same but the relative variance term would grow with exponential rate of $\inf_{x\in B}I(x)$. Therefore the estimator $s_{n}$ has theoretically superior mean square error to $mc_{n}(z_{n})$. Another point to mention, is that analysis of $rMSE(s_{n})$ leaves out the computational work for generating a single copy of $s_{n}$. However, Proposition 5.3 guarantees that the expected work per copy of $s_{n}$ grows subexponentially and thus when normalized for expected work per replicate $s_{n}$ will still be a superior estimator to $mc_{n}(z_{n})$.

Proposition 5.2 is unfortunately not strong enough to provide precise control over the relative bias term, but we do know that if it does grow it does so in a subexponential fashion. It is our hope that a stronger statement of the form (5.3) can be achieved.

6 Numerical Examples

6.1 Splitting algorithm

Before we run the splitting algorithm for a specific example, there are some initial steps needed to be done. We first obtain a sub-solution $\bar{T}$ to the associated variation problem. If the problem is two dimensional, we let $T(v)=I(v)$, $\forall v\in\mathbb{R}_{+}^{2}$, where $I(v)$ is the explicit solution to the VP discussed in Section 4.1. Then we let $\bar{T}(v)=-T(v)+\inf_{x\in B}T(x)$. If the problem is three or higher dimensional, we can find the scaling factor $\mathbf{r}$ by solving the corresponding optimization problem with the help of any optimization solver, and then construct the sub-solution as defined in (4.5) of Section 4.2. Then we choose a positive integer $r>1$ as the number of new particles when a splitting occurs, and choose $\delta\in(0,1]$ as the level size. The importance function is then obtained as $V(z)=\delta\bar{T}(z)/\log(r)$. We then define a collection of level sets $\{C_{j}^{n}:j\in\mathbb{N}\}$, where $C_{0}^{n}=B,C^{n}_{j}=L_{(j-1)\delta/n},j\in\mathbb{N}$, where $L_{a}=\{z:V(z)\leq a\}$. Finally, for a point $x$, its level can be calculated through the level function $l_{n}(x)=\min\{j\geq 0;x\in C^{n}_{j}\}$.

We conduct numerical experiments for one 2-dimensional SRBM, and one 3-dimensional SRBM. The first splitting simulation used 1,000 runs with step size $1/(1000*n)$. The second splitting simulation used 200,000 runs with step size $1/(1000*n)$. Note that our step size does not comply with the requirement in Proposition 5.2 to achieve the asymptotic error bound. However, for every value of $n$ used in the following simulations, $1/(1000*n)\leq 1/(n^{2.5})$, where the latter step size satisfies the requirement in Proposition 5.2. In this way, we increase the accuracy of simulation results when $n$ is relatively small, and the polynomial factor is not too large as $n$ increases. Simulations are run in Matlab with computer having Intel i7-7500U CPU @ 2.70GHz, and RAM 8.00 GB(7.88 GB useable).

For a $2$-dimensional SRBM, we consider the initial position $(0.1,0.1)$, $R=[1,0;-1,1]$, $\theta=[-2;1]$, $\Sigma=[1,0;0,1]$, $A_{n}=\{z\in\mathbb{R}^{2}_{+}:\sum_{k=1}^{2}z_{k}\leq 0.15/n\}$, $B=\{z\in\mathbb{R}^{2}_{+}:\sum_{k=1}^{2}z_{k}\geq 1\}$.

[TABLE]

For a $3$-dimensional SRBM, consider the initial position $(0.1,0.1,0.1)$, $R=[3,-1,-1;$ $-1,2,-1;-1,-1,2]$, $\theta=[-2;-1;-1]$, $\Sigma=[2,1,1;1,2,1;1,1,3]$, $A_{n}=\{z\in\mathbb{R}^{3}_{+}:\sum_{k=1}^{3}z_{k}\leq 0.15/n\}$, $B=\{z\in\mathbb{R}^{3}_{+}:\sum_{k=1}^{3}z_{k}\geq 1\}$.

[TABLE]

Remark 6.1

Splitting algorithm does not always work well for 3-dimensional SRBM due to the fact that the subsolution might not be good enough in some cases.

We also ran the standard Monte Carlo simulation for the 3-dimensional example with $n=5$ to compare with the results generated by splitting algorithm. We ran splitting algorithm $100,000$ times so that $100,000*r^{l(x_{0})-1}=1,000,000$. We also ran standard Monte Carlo $1,000,000$ times.

[TABLE]

6.2 RESTART algorithm

We applied the RESTART algorithm to simulate the same $2$-dimensional example as in the previous section. As expected, the RESTART algorithm required much less computational effort. Hence, we could run more trials, and simulate the cases for high value of $n$ which would be very time consuming when splitting algorithm is applied.

Here we give numerical results for a $2$-dimensional SRBM. The simulation used 10,000 runs with step size $1/(1000*n)$. Simulations are run in Matlab with computer having Intel i7-7500U CPU @ 2.70GHz, and RAM 8.00 GB(7.88 GB useable).

For $2$-dimensional SRBM, consider the initial position $(0.1,0.1)$, $R=[1,0;-1,1]$, $\theta=[-2;1]$, $\Sigma=[1,0;0,1]$, $A_{n}=\{z\in\mathbb{R}^{2}_{+}:\sum_{k=1}^{2}z_{k}\leq 0.15/n\}$, $B=\{z\in\mathbb{R}^{2}_{+}:\sum_{k=1}^{2}z_{k}\geq 1\}$.

[TABLE]

[TABLE]

7 Discussion

In this project, we developed particle based simulation algorithms for estimating rare event probability related to SRBMs in a non-negative orthant. We found that splitting-type algorithms provide much better performance than standard Monte Carlo method. Furthermore, we found that the RESTART splitting algorithm is superior to regular splitting algorithm in terms of operating time. Following [8, 9] our algorithms are based on subsolutions to the variational problems associated with our rare events. In particular, the work [1] developed solutions to these variational problems in two dimensions. For three and higher dimensions we construct new subsolutions to the variational problem. Due to the results of [8, 9] we are able to show that our splitting algorithm has superior performance to standard Monte Carlo in two and higher dimensions.

Future directions for this work include developing subsolutions in three and higher dimension that are closer to the solution, and thus enable simulation algorithms with better performance. Another possible direction is to consider more general reflecting diffusions, i.e., reflecting diffusions with state dependent drift and variance.

8 Proofs

We provide all the proofs in this section.

8.1 Proofs for Section 3

To show Proposition 3.1, we introduce the following Lyapunov function introduced in [13]. Recall that $R_{i}$ is the $i$th column of the reflection matrix $R$, and for a $d\times d$ matrix $M$, its trace is denoted as $tr(M)$.

Theorem 8.1** ([13])**

Under Assumption 2.2, and the assumption that $R$ is completely-$\mathcal{S}$, there exists a Lyapunov function $L:\mathbb{R}^{K}_{+}\to\mathbb{R}_{+}$ such that

• (i)

$L\in\mathcal{C}^{2}(\mathbb{R}^{d}_{+}\backslash\{0\})$;

• (ii)

Given $M_{0}>0$, there exists an $M_{1}>0$ such that when $|z|\geq M_{1}$, we have $L(z)\geq M_{0}$;

• (iii)

Given $\varepsilon>0,$ there exists a $M_{2}>0$ such that when $|z|\geq M_{2}$, we have $|\ddot{L}(z)|\leq\varepsilon$;

• (iv)

There exists $c_{0}>0$ such that

[TABLE]

• (v)

$L(\alpha z)=\alpha L(z)$* for all $\alpha\geq 0$ and $z\in\mathbb{R}^{d}_{+}.$*

Some consquences of the above properties are obtained in [4], which are listed as follows.

• (vi)

For every $M_{3}>0$, there exists a $\gamma(M_{3})>0$ such that $\sup_{|z|\leq M_{3}}L(z)\leq\gamma(M_{3}).$

• (vii)

There exists an $M_{4}>0$ such that $\sup_{z\in\mathbb{R}^{d}_{+}\backslash\{0\}}|\dot{L}(z)|\leq M_{4}$.

• (viii)

There exists $0<c_{1}\leq c_{2}<\infty$ such that $c_{1}|z|\leq L(z)\leq c_{2}|z|.$

Lemma 8.1** ([4])**

Fix $x\in\mathbb{R}^{d}_{+}$ and $\Delta>0$. For $m\in\mathbb{N}$, define

[TABLE]

Then for any $\kappa>0$ and $m,n\in\mathbb{N},m\leq n$, we have

[TABLE]

where $M_{4}$ is as in Theorem 8.1 (vii) and $\gamma_{0}$ is a positive constant only depending on the norm of $\Sigma^{-1}.$

Proof of Proposition 3.1. Define $\tau_{\epsilon}=\inf\{t\geq 0:\sum_{k=1}^{d}Z_{k}(t)\leq\epsilon\}$ and then $n\tau_{n}\leq\tau_{\epsilon}$. Thus it suffices to show that there exists $c>0$ such that

[TABLE]

The proof idea is adapted from that of Theorem 4.3 of [4]. From Theorem 8.1 (iii), choose $M_{2}>0$ such that when $\sum_{k=1}^{d}z_{k}\geq M_{2}$, $tr(\ddot{L}(z)\Sigma)<c_{0}$, where $c_{0}$ is as in Theorem 8.1 (iv). Let $M_{0}>M_{2}$, and define $\kappa_{\epsilon}=M_{0}/\epsilon.$ For $n\in\mathbb{N}$, consider the set $A_{n}=\{\omega\in\Omega:\inf_{0\leq s\leq n\Delta}\sum_{k=1}^{d}Z_{k}(s)>\epsilon\}.$ For $\omega\in A_{n}$ and $m\leq n$, from Ito’s formula, we have

[TABLE]

Multiplying both sides of the above equation by $\kappa_{\epsilon}$, from Theorem 8.1 (iv), (v) and (viii), we have for $m\leq n$, on $A_{n}$,

[TABLE]

which implies that

[TABLE]

From (8.1) and using Lemma 8.1, we have for any $\kappa>0,$

[TABLE]

We can choose $\Delta$ large enough and $\kappa$ small enough such that $-\eta\equiv\log(2\sqrt{2})/\Delta+\kappa^{2}M_{4}^{2}\gamma_{0}-\kappa c_{0}/2<0.$ Therefore, we have

[TABLE]

For $t>0$, there exists $n_{0}\in\mathbb{N}$ such that $t\in[n_{0}\Delta,(n_{0}+1)\Delta]$, and we have

[TABLE]

When $z\in\mathbb{R}^{d}_{+}$ satisfies that $\sum_{k=1}^{d}z_{k}\leq\epsilon$ or $\sum_{k=1}^{d}z_{k}\geq n$, we have $\tau_{n}=0$. When $\epsilon<\sum_{k=1}^{d}z_{k}<n$, we have for $c\in(0,\eta)$,

[TABLE]

The result follows.

Proof of Theorem 3.1 We first show that

[TABLE]

From Proposition 3.1, it follows that

[TABLE]

where $\mathbb{P}_{z_{n}}$ and $\mathbb{E}_{z_{n}}$ are the probability and expectation conditioning on $Z_{n}(0)=z_{n}.$ Therefore we can take $T>0$ such that

[TABLE]

The following proof is similar to that of Lemma 5.7.21 in [10]. Define the following closed set of $\mathcal{C}([0,T];\mathbb{R}^{d}_{+})$:

[TABLE]

Then

[TABLE]

Also note that for $T>0$ we have

[TABLE]

According to Lemma 1.2.16 in [10],

[TABLE]

Now we show the lower bound

[TABLE]

We next note that

[TABLE]

where $\tau_{\psi}=\inf\{t\geq 0:\Gamma(\psi)(t)\in B\}.$ Thus it suffices to establish that

[TABLE]

for any $\psi\in\mathcal{AC}([0,\infty),\mathbb{R}^{d})$ such that $\psi(0)=0$ and $\tau_{\psi}<\infty.$ Noting that the possible optimal solutions of the VP in (3.2) never return to [math] before reaching $B$, we only need to consider $\psi$ such that $\sum_{k=1}^{d}\Gamma(\psi)_{k}(s)>0$ for all $s\in(0,\tau_{\psi}].$ Fix such a $\psi$. Let $\varphi^{\tilde{z}}(t)=\tilde{z}+\psi(t)$, with $\sum_{k=1}^{d}\tilde{z}_{k}>0.$ It is clear that

[TABLE]

Fix $\kappa>0$, and let $B_{\kappa}=\{\tilde{z}\in\mathbb{R}^{d}_{+}:\sum_{k=1}^{d}\tilde{z}_{k}=\kappa\}.$ For $\kappa,\delta>0,$ define

[TABLE]

For $\kappa\in(0,1)$ there exists $\delta_{0}>0$ and $N\in\mathbb{N}$ such that, when $0<\delta\leq\delta_{0}$,

[TABLE]

for all $t\in[0,\tau_{\psi}]$, and

[TABLE]

Assume that $n$ is large enough such that $\kappa>\sum_{k}z_{n,k}>\epsilon/n$. From the Markov property of $Z_{n}$, we have that

[TABLE]

Using the minimality property of the one-dimensional Skorokhod map, we see that $\sum_{k=1}^{d}Z_{n,k}(t)\geq\Gamma_{1}(\sum_{k=1}^{d}z_{n,k}+\sum_{k=1}^{d}\tilde{X}_{n,k})(t)$ for all $t\geq 0$, where $\tilde{X}_{n}(t)=\theta t+X(nt)/n$ and $\Gamma_{1}$ is the one-dimensional Skorokhod map. Thus we have

[TABLE]

where $\gamma=\frac{2|\sum\theta|}{\mathbf{1}^{\prime}\Sigma\mathbf{1}},a_{n}=\kappa-\sum_{k=1}^{d}z_{n,k},$ and $b_{n}=\sum_{k=1}^{d}z_{n,k}-\frac{\epsilon}{n}$. Thus

[TABLE]

We next note that

[TABLE]

Letting $\delta\downarrow 0$, we conclude that

[TABLE]

Finally, letting $\kappa\downarrow 0$, the result follows.

8.2 Proofs for Section 4

Proof of Lemma 4.1. For any $w,v\in\mathbb{R}^{d}_{+}\backslash(A\cup B)$, suppose we have an optimal triple $(\psi,\eta,\phi)$ with $\phi(0)=w,\phi(T)=v$. Denote by $\mathcal{P}(I)$ the power set of $L$. When $\phi$ traverses a face $F_{K}$ from time $t_{1}$ to $t_{2}$ for some $K\in\mathcal{P}(I)$, the segment of $(\psi,\eta,\phi)$ must be locally optimal. Otherwise, we can find a path which has lower cost. For some $K\in\mathcal{P}(I)$, we define the set $G_{K}=\{t\in[0,T]:\phi_{i}(t)=0,\phi_{j}(t)>0,\forall i\in K,\forall j\in L\backslash K\}^{o}$, where for a set $E$, $E^{o}$ denotes its interior. We notice that $G_{K}$ is a open set, and can be represented as the disjoint union of at most countably many open intervals. Denote the set of these open intervals by $\mathbb{O}_{K}$. Notice that the optimal path traverses the face $F_{K}$ from time $t_{1}$ to $t_{2}$ if $(t_{1},t_{2})\in\mathbb{O}_{K}$. Observe that $\cup_{K\in\mathcal{P}(I)}\mathbb{O}_{K}$ contains at most countably many disjoint open intervals, and

[TABLE]

Both of the equalities can be justified by the argument that there are at most countably many singleton points in the set $[0,T]\setminus\{\cup_{K\in\mathcal{P}(I)}G_{K}\}$, which contributes [math] to the cost.

Proof of Theorem 4.2. Theorem 4.2 summarizes some major results of [14], including Theorem 3, Theorem 4, Lemma 5, Proposition 1, and Proposition 2. Despite a small difference in problem settings, their proofs are still valid for our setting. To be more specific, given $K\subseteq L$, and two points $v,w\in\mathbb{R}^{d}_{+}$ which satisfy Condition 4.1 with respect to $K$, Theorems 3 and 4 of [14] assert that there is a unique $J_{*}\subseteq K$ such that $\langle A^{J_{*}}R_{j},\alpha^{J_{*}}(w,v)(v-w)-\theta\rangle_{D}\leq 0,\ j\in K\backslash J_{*}$, and when $J_{*}\neq\emptyset$, each component of ${\lambda}^{J_{*}}(w,v)$ is positive. Proposition 1 of [14] shows that $\mathcal{I}^{*}_{K}(w,v)=|A^{J_{*}}\theta|_{D}|A^{J_{*}}(v-w)|_{D}-\langle A^{J_{*}}\theta,A^{J_{*}}(v-w)\rangle_{D}$. Lemma 5 and Proposition 2 of [14] give the explicit form of an optimal path.

Proof of Lemma 4.2. See the proof of Lemma 1 in [14].

Proof of Lemma 4.3. (i) If $\mathcal{I}^{*}_{K}(w,v)=0$, then by Theorem 4.2, there exists an optimal path such that $\dot{\psi}(t)=\theta<0$ for all $t\in(0,\tau_{\phi}(v))$, where $\tau_{\phi}(v)=\inf\{t\geq 0:\phi(t)=v\}$ given that $\phi(0)=w$. Hence

[TABLE]

where $\eta_{i}(\tau_{\phi}(v))>0,\forall i\in{K}$, and $\eta_{i}(\tau_{\phi}(v))=0,\forall i\in L\backslash{K}$ since $\phi(t)\in F_{K}$ for $t\in(0,\tau_{\phi}(v))$. From Theorem 4.2, we have

[TABLE]

and

[TABLE]

Multiply both sides of the second equality in (8.6) by $\tau_{\phi}(v)$, and we can get

[TABLE]

Comparing the above equation with (8.5), we obtain that

[TABLE]

Since $R$ is invertible, all its columns are linearly independent. If $J_{*}\neq K$, for $i\in K\setminus J_{*}$, $\eta_{i}(\tau_{\phi}(v))=0$, which leads to a contradiction. Thus, we can conclude that $J=K$.

We next note that $v_{i}-w_{i}<0$ for $i\in L\backslash K$ from (8.5) and the fact that $\theta<0$; which yields that $v_{i}-w_{i}\leq 0,$ for all $i\in L$. Noting that both $w$ and $v$ are in $\mathbb{R}^{d}_{+}$, $T(v)\leq T(w)$ follows from the monotonicity of $T(\cdot)$.

(ii) Again from Theorem 4.2, we can find two optimal paths such that $\dot{\psi}_{1}(t)=\theta,t\in(0,\tau(v))$ and $\dot{\psi}_{2}(t)=\theta,t\in(0,\tau(w))$, where $\tau(v)=\inf\{t\geq 0:\phi_{1}(t)=v\}$ and $\tau(w)=\inf\{t\geq 0:\phi_{2}(t)=w\}$ given that $\phi_{1}(0)=\phi_{2}(0)=0$. For simplicity, we ignore the symbols of particular optimal path in the following proof. Similar to (8.5), we have

[TABLE]

where $\eta_{i}(\tau(v))>0,\eta_{i}(\tau(w))>0,\forall i\in{K}$, and $\eta_{i}(\tau(v))=\eta_{i}(\tau(w))=0,\forall i\in L\backslash{K}$. Let $\phi_{K}=(\phi_{i};i\in K),\eta_{K}=(\eta_{i};i\in K),$ and $\theta_{K}=(\theta_{i};i\in K)$. We first consider $\mathcal{I}^{*}_{K}(v)$, noting that $\eta_{i}(t)=0$ for $i\in L\backslash K$ and $t\in(0,\tau(v))$, we have for $t\in(0,\tau(v)]$,

[TABLE]

which gives that $\eta_{K}(t)=-R^{-1}_{KK}\theta t.$ Similarly, for $\mathcal{I}^{*}_{K}(w)$, we have $\eta_{K}(t)=-R^{-1}_{KK}\theta t$ for $t\in(0,\tau(w)].$ It now follows that for $i\in K,$

[TABLE]

and for $i\in L\backslash K,$

[TABLE]

The result follows.

Proof of Proposition 4.1. It suffices to show that for any $K\subset L$, there exists an $\epsilon>0$ such that for any $v\in\mathbb{D}_{K}$, $\mathcal{I}^{*}_{K}(v)>\epsilon$. From the extreme value theorem, $\mathcal{I}^{*}_{K}(v)$ attains the minimum in $\mathbb{D}_{K}$. It is also known that $\mathcal{I}^{*}_{K}(v)\neq 0$ for any $v\in\mathbb{D}_{K}$. To see this, if there exists a $v\in\mathbb{D}_{K}$ such that $\mathcal{I}^{*}_{K}(v)=0$, then from Lemma 4.3 (i), $v_{j}<0$ for all $j\in L\setminus K$, which contradicts the property (3) in the definition of $\mathbb{D}_{K}$. Hence there exists some positive number $\epsilon$ such that $\mathcal{I}^{*}_{K}(v)>\epsilon$ for any $v\in\mathbb{D}_{K}$.

Proof of Theorem 4.3. From Lemma 4.1, it suffices to show that $\bar{T}(x)\leq 0$ for all $x\in B$, $\bar{T}$ is continuous, and that for all $w,v\in\mathbb{R}^{d}_{+}\backslash B$, which satisfy Condition 4.1 for some $K\subset L$, we have

[TABLE]

Since the first two conditions are obviously satisfied we only need establish the third condition. First assume $v_{i}-w_{i}<0,\forall i\in L\backslash K$, then

[TABLE]

Now if there exists at least one $j\in L\backslash K$ such that $v_{j}-w_{j}\geq 0$, we have $(v-w)/|v-w|\in\mathbb{D}_{K}$, and

[TABLE]

where the first inequality follows since $T$ is a norm, the second inequality follows from the definition of $\mathbf{r}$ and the final equality follows from Theorem 4.2.

8.3 Proofs for Section 5

8.3.1 Proofs of Propositions 5.1 and 5.2

Proof of Proposition 5.1. Let $W$ be a standard Brownian motion, and for a fixed $\Delta>0$, define its discretization $W^{\Delta}(t)=W((k-1)\Delta)$ for $t\in[(k-1)\Delta,k\Delta)$ and $k\in\mathbb{N}.$ Let $T>0$, and consider $\mathbf{W}_{T}=\{W(t):0\leq t\leq T\}$ and $\mathbf{W}_{T}^{\Delta}=\{W^{\Delta}(t):0\leq t\leq T\}$. We are interested in the following discretization error.

[TABLE]

Define $V(\Delta)=\max_{0\leq t\leq\Delta}|W(t)|$. By the Markov property of Brownian motion, we have that

[TABLE]

where $\{V_{k}(\Delta):k\in\{1,\ldots,\left\lfloor{T/\Delta}\right\rfloor\}$ is an independent sequence of copies of $V(\Delta)$. From [21], an upper bound on the tail probabilities of $V(\Delta)$ can be derived. In particular, for $z>0$,

[TABLE]

It follows that for $z>0$

[TABLE]

where $T(\Delta)=\left\lfloor{T/\Delta}\right\rfloor$. Next we use the inequality that $\log(z)\geq-(1-z)/z$ for $z\in(0,1)$ to conclude that

[TABLE]

Considering the exponent in the above equation, we have

[TABLE]

and thus

[TABLE]

Now consider a scaled version of $\mathbf{W}_{T}$ defined by $W_{n}(t)=\frac{1}{n}W(nt),t\in[0,T]$, which is denoted by $\mathbf{W}_{n,T}=\{W_{n}(t):0\leq t\leq T\}$. For $\Delta>0$ we consider the discretization $\mathbf{W}^{\Delta}_{n,T}=\{W_{n}^{\Delta}(t):0\leq t\leq T\}$ defined by $W_{n}^{\Delta}(t)=W_{n}((k-1)\Delta)$ for $t\in[(k-1)\Delta,k\Delta).$ Finally, we observe that

[TABLE]

Therefore we can use (8.9) to obain

[TABLE]

Under Assumption 4.2, the Skorokhod map is Lipschitz continuous, and there exists a constant $C$ which only depends on $R$, such that for $t\geq 0$,

[TABLE]

For $z>0$, we have

[TABLE]

Denote $\tilde{z}_{\Delta}=(z/C-\Delta|\theta|)/(d|\Sigma^{1/2}|).$ From (8.9), we have

[TABLE]

Now consider the large derivation scaled processes $Z_{n}(t)=Z(nt)/n,X_{n}(t)=X(nt)/n$, and $Y_{n}(t)=Y(nt)/n$. From the homogenity property of the Skorokhod map for linear scalings, $(Z_{n},Y_{n})$ is the solution of the Skorokhod problem for $\{\tilde{X}_{n}(t)=X_{n}(t)+\theta t;t\geq 0\}$. It follows from (8.10) that

[TABLE]

Hence

[TABLE]

When $\Delta=\Delta(n)=o(1)$, we have

[TABLE]

To show Proposition 5.2, we introduce the following notation. Let $\eta\in(0,1/2)$, and take $\Delta\equiv\Delta(n)$ small enough such that $\epsilon/n-{\Delta}^{\eta}>0$ and $\epsilon/n+{\Delta}^{\eta}<z_{n}$.

Lemma 8.2

If $\Delta(n)\rightarrow 0$ as $n\rightarrow\infty$, then for $\eta\in(0,1/2)$,

[TABLE]

Proof. From (8.12) in the proof of Proposition 5.1, letting $\xi=\frac{{{\Delta}^{\eta}}/C-\Delta|\theta|}{d|\Sigma^{1/2}|}$, we have that for $T\geq 0,$

[TABLE]

Next noting that

[TABLE]

it follows that

[TABLE]

We next observe that for $T>0$

[TABLE]

and therefore by Lemma 1.2.16 in [10]

[TABLE]

The result then follows from Proposition 3.1 and sending $T\to\infty$.

Define:

[TABLE]

and

[TABLE]

At last, for $\tilde{z}\in\mathbb{R}^{d}_{+}$, let $p_{n}^{-}(\tilde{z})=\mathbb{P}(\tau_{n}^{-}(B_{n}^{-})<\tau_{n}^{-}(A_{n}^{-})|Z_{n}(0)=\tilde{z}).$

Lemma 8.3

For $\eta\in(0,1/2)$, let $\Delta\equiv\Delta(n)$ satisfy $\lim\limits_{n\rightarrow\infty}n\Delta^{\eta}=0$. Then

[TABLE]

where $\mathcal{C}_{n}=\{\tilde{z}\in\mathbb{R}^{d}_{+}:\sum_{i=1}^{d}\tilde{z}_{i}=\sum_{i=1}^{d}z_{n,i}\}.$

Proof. Define the boundaries of $A^{-}_{n},B^{-}_{n}$ and $A_{n}$ inside $\mathbb{R}^{d}_{+}$ as follows:

[TABLE]

and

[TABLE]

Using the Markov property,

[TABLE]

We observe that $\{\tau^{-}_{n}(B_{n}^{-})<\tau_{n}(A_{n})|Z(0)=z_{n}\}\subset\{\tau_{n}^{-}(B^{-}_{n})<\tau_{n}^{-}(A_{n}^{-})|Z(0)=z_{n}\}$, and thus

[TABLE]

Next note that for $\tilde{z}\in\mathcal{A}_{n},$

[TABLE]

where $\tau_{n}(\mathcal{C}_{n})=\inf\{t\geq 0:Z_{n}(t)\in\mathcal{C}_{n}\}$. It follows that

[TABLE]

which yields that

[TABLE]

Taking maximum over $z_{n}\in\mathcal{C}_{n}$, we have

[TABLE]

and

[TABLE]

Similar to the proof of Theorem 3.1, using the minimality property of the one-dimensional Skorokhod map $\Gamma_{1}$, we have

[TABLE]

where $\gamma=\frac{2|\sum\mu|}{\mathbf{1}^{\prime}\Sigma\mathbf{1}}$ and $a_{n}=1-\Delta^{\eta}-\epsilon/n.$

We next use the Lyapunov function $L$ introduced in Theorem 8.1 to derive an upper bound for $\max_{\tilde{z}\in\mathcal{A}_{n}}\mathbb{P}(\tau_{n}(C_{n})<\tau_{n}^{-}(A_{n}^{-})|Z_{n}(0)=\tilde{z})$. Choose $\kappa_{1}>0$ such that when $\sum_{i=1}^{d}z_{i}\geq\kappa_{1}$, $tr(\ddot{L}(z)\Sigma)<c_{0}$. Let $\kappa_{2}=2\kappa_{1}/\epsilon$, and define $\check{Z}(t)=\kappa_{2}Z(t)$ for $t\geq 0$. Then for $\tilde{z}\in\mathcal{A}_{n}$,

[TABLE]

From Theorem 8.1 (v), for $\check{z}\in\mathbb{R}^{d}_{+},$ $L(\check{z})=(\sum_{i=1}^{d}\check{z}_{i})L(\check{z}/\sum_{i=1}^{d}\check{z}_{i})$. Let $\underline{l}=\inf\{L(z):z\in\mathbb{R}^{d}_{+},\sum_{i=1}^{d}z_{i}=1\}$. Then $\underline{l}\in(0,\infty)$, and

[TABLE]

We further define the following two sets

[TABLE]

We will see that $\tilde{A}\subset\kappa_{2}nA^{-}_{n}$ and $\kappa_{2}nC_{n}\subset\tilde{C}$, where $C_{n}=\{\check{z}\in\mathbb{R}^{d}_{+}:\sum_{i=1}^{d}\check{z}_{i}\geq\sum_{i=1}^{d}z_{n,i}\}.$ Indeed, for $\check{z}\in\tilde{A}$, from (8.19),

[TABLE]

which says $\check{z}\in\kappa_{2}nA^{-}_{n}$. Now for $\check{z}\in\kappa_{2}n\mathcal{C}_{n}$, that is $\sum_{i=1}^{d}\check{z}_{i}\geq\kappa_{2}\sum_{i=1}^{d}z_{i}$, again from (8.19), we have

[TABLE]

Fix $\tilde{z}\in A_{n}$. Hence from (8.19),

[TABLE]

We next consider the following martingale: For $\lambda>0$,

[TABLE]

where

[TABLE]

Let $\tau(\tilde{z})=\inf\{t\geq 0:\kappa_{3}L(\check{Z}(t))\in\tilde{A}\cup\tilde{C}\}.$ Now let $m=n=1$ and $\Delta=t$ in (8.1) in the proof of Proposition 3.1, we have for positive constants $\kappa_{4}$,

[TABLE]

From Theorem 8.1 (vii), $|\dot{L}(Z(t))|\leq M_{4}$, and we can choose $\lambda$ such that

[TABLE]

Using optional sampling theorem, we have

[TABLE]

Hence

[TABLE]

Applying (8.20) and (8.18) to (8.16), we have that

[TABLE]

This completes the proof.

Proof of Proposition 5.2: We note that

[TABLE]

Recall that $\tau_{n}^{-}=\inf\{t\geq 0:\sum_{k=1}^{d}Z_{n,k}(t)\geq 1-\Delta^{\eta},\ \mbox{or}\ \sum_{k=1}^{d}Z_{n,k}(t)\leq\epsilon/n-\Delta^{\eta}\}$, and $p_{n}^{-}(z_{n})=\mathbb{P}(\sum_{k=1}^{d}Z_{n,k}(\tau_{n}^{-})=1-\Delta^{\eta}|Z_{n}(0)=z_{n}).$ We note that Lemma 8.2 is still valid when $\tau_{n}$ is replaced by $\tau_{n}^{-}$. Thus, from the proof of Lemma 8.2, we have

[TABLE]

where $\xi=\frac{{{\Delta}^{\eta}}/C-\Delta|\theta|}{d|\Sigma^{1/2}|}.$ Also Theorem 3.1 is valid for $p_{n}^{-}(z_{n})$, which yields that

[TABLE]

We observe that

[TABLE]

Noting that $\frac{\xi^{2}}{2\Delta}=\Theta(\Delta^{2\eta-1}),$ for some constant $C>0$, we have

[TABLE]

At last from Lemma 8.3, we obtain

[TABLE]

Combining the above estimates, we establish that

[TABLE]

To show the other direction, we define $\tau_{n}^{+}=\inf\{t\geq 0:\sum_{k=1}^{d}Z_{n,k}(t)\geq 1+\Delta^{\eta},\ \mbox{or}\ \sum_{k=1}^{d}Z_{n,k}(t)\leq\epsilon/n+\Delta^{\eta}\}$, and $p_{n}^{+}(z_{n})=\mathbb{P}(\sum_{k=1}^{d}Z_{n,k}(\tau_{n}^{+})=1+\Delta^{\eta}|Z_{n}(0)=z_{n}).$ Similar to Lemma 8.3 (replace $p^{-}_{n}(z_{n})$ and $p_{n}(z_{n})$ by $p_{n}(z_{n})$ and $p^{+}_{n}(z_{n})$, respectively, in the proof), it can be shown that

[TABLE]

Further note that

[TABLE]

Finally, we have

[TABLE]

8.3.2 Proofs of Propositions 5.3 and 5.4

Proof of Theorem 5.3. The stability follows directly from Proposition 7 in [8], and (5.5) follows from Lemma 5 in the same reference. To prove (5.6), we need to verified that Conditions 1 and 4 in [8] hold for our problem. Condition 1 of [8] follows from Theorem 3.1. For their Condition 4, due to the positive recurrence of our defined process, $\sigma_{n}$ is finite almost surely. Then by strong Markov property

[TABLE]

Hence using the strong Markov property of $Z_{n}$, we have

[TABLE]

By the large deviation upper bound, we can get

[TABLE]

Now (5.6) follows directly from Theorem 8 in [8].

Proof of Proposition 5.4: It suffices to verify the five conditions in [9], which are Conditions 3.1, 4.3, 5.1, 5.2, 5.3, are satisfied. We note that for our specific problem, those conditions essentially require that:

• •

$A_{n}$ and $B$ are closed sets.

• •

$Z_{n}^{\Delta}(t)$ satisfies the large deviation principle.

• •

For any compact $K\subset\mathbb{R}^{d}_{+}$, all zero cost trajectories with initial positions in $K$ will enter $(A_{n}\cup B)^{\circ}$ by some fixed finite time depended on $K$.

Indeed all these conditions are satisfied, and thus the proposition follows.