On the Ergodic Properties and Invariant Measure of a Two-Dimensional Reflected Ornstein–Uhlenbeck Process

1. Introduction

Reflected stochastic processes have attracted significant attention in recent years because of their close connections to queueing systems that arise in computer networks, telecommunications, mathematical biology, and transportation problems. Among the various areas of research, investigations into their ergodic properties have made striking progress. For instance, the recurrence of reflected Brownian motions and reflected diffusions has been widely explored in the literature, including, but not limited to, the references Williams (1985), Hobson and Rogers (1993), Dupuis and Williams (1994), Chen (1996), Budhiraja and Dupuis (1999), and Atar et al. (2001). Additionally, the stationary distributions of two-dimensional reflected Brownian motions have been thoroughly investigated, with explicit analytical expressions available. Relevant references include Harrison et al. (1985), Williams (1985), Dupuis and Ramanan (2002), Dieker and Moriarty (2009), Dai and Miyazawa (2011), and Franceschi and Raschel (2019). However, in higher dimensions, results on the stationary distribution are limited. Harrison and Williams (1987b) showed that under a certain skew symmetry condition, the stationary densities of reflected Brownian motions take an exponential product form. Despite this, Kang and Ramanan (2014) provided a characterization of stationary distributions of reflected diffusions.

However, when focusing specifically on the reflected Ornstein–Uhlenbeck process, the existing literature concentrates on the one-dimensional case. For example, Ward and Glynn (2003b) conducted a comprehensive study of the one-dimensional reflected Ornstein–Uhlenbeck (1-d ROU) process, including its explicit stationary density, a perturbation expansion for the transition density, and approximations for level crossing times. Furthermore, stationary distributions have been explored for certain generalized 1-d ROU processes. Specifically, Xing et al. (2009) investigated the 1-d ROU process with jumps and the 1-d Markov-modulated ROU process, whereas Zhang and Jiang (2009) studied the 1-d ROU process with two-sided barriers. Therefore, the main purpose of this paper is to extend the scope of research by studying a special case of the two-dimensional reflected Ornstein–Uhlenbeck process, as introduced below.

The two-dimensional reflected Ornstein–Uhlenbeck (2-d ROU) process $(X,Y)$ in ${\mathbb{R}}_{+}^2$ studied in this paper is defined as follows:

with initiation condition $(X_0,Y_0)=(x_0,y_0)\in {\mathbb{R}}_{+}^2$, where $({\mu }_1,{\mu }_2)\in {\mathbb{R}}^2$, $({\theta }_1,{\theta }_2)>0$, $({\sigma }_1,{\sigma }_2)>0$, $(W,B)$ is a two-dimensional Brownian motion with covariance matrix $\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right)$ with $0<\rho <1$, and $(L_t^X,L_t^Y)$ denotes the local times at the boundaries and corresponds to the coefficients in front of the normal reflections. Specifically, $(L_0^X,L_0^Y)=(0,0)$, $L_t^X$, and $L_t^Y$ are nondecreasing, and ${\int }_0^{\infty }{\mathbf{1}}_{(0,\infty )}(X_s)d{L}_s^X=0$ and ${\int }_0^{\infty }{\mathbf{1}}_{(0,\infty )}(Y_s)d{L}_s^Y=0$. The existence and uniqueness of a strong solution to the 2-d ROU process are guaranteed by a careful extension of the results of Lions and Sznitman (1984).

Before studying our 2-d process, we first provide an overview of the emergence of ROU processes. These processes arise as the diffusion scaling limit of a sequence of single-server queues with balking or reneging under heavy-traffic conditions (See Ward and Glynn 2003a, 2005 for details). Such queueing models have broad applications across various domains, including telecommunications, supply chain management, and computing and cloud services.

Beyond their derivation from such queueing models, ROU processes have also been shown to serve as effective approximations in various settings. Borovkov (1984) and Srikant and Whitt (1996) demonstrate that ROU processes approximate the number-in-system process in a G/M/s/0 queueing model when both the number of servers and the arrival rate are large. Additionally, Boxma et al. (2016) establish that the ROU process provides a heavy-traffic approximation for a sequence of reflected AR(1) processes.

However, these results focus primarily on the one-dimensional case. In contrast, our 2-d ROU process arises as the diffusion scaling limit for parallel single-server queues with abandonment in the critically loaded regime, which have a randomly split Hawkes process or a bivariate Hawkes process as the arrival processes (Li and Pang 2024). See Figure 1 for an illustration of these queueing models with split Hawkes arrivals. The correlation between these parallel queues emerges from dependencies in the arrival processes, which, in turn, translate into the correlated Brownian motions that drive the dynamics of the ROU process. Our 2-d ROU process is, of course, not the most general form of reflected OU processes; particularly, the drift is linear, the covariance coefficient is constant, and the reflection is normal. Furthermore, because there are no other interactions among the queues other than the arrivals, the reflection terms will be exactly like the reflection for each separate single-server queue with abandonment, as studied in Ward and Glynn (2005).

The analysis of our 2-d ROU process contributes to the understanding of stationary distributions and transient performance measures of these parallel server queueing models, where explicit expressions are often intractable. Furthermore, we believe our results may provide valuable insights into the stationary workload distribution of bivariate M/G/1 systems with coupled input and parallel service, as posed in Mandjes (2022).

For the one-dimensional ROU process, the ergodic properties are well understood, as discussed in Ward and Glynn (2003b), where the invariant measure is explicitly shown to have a truncated Gaussian distribution. Additionally, the exponential rate of convergence in total variation distance is studied in Lund et al. (1996), and the rate in Wasserstein distance is investigated in Sarantsev (2020). In this paper, our goal is to understand such ergodic properties for the simple 2-d ROU process arising from the parallel server queueing models mentioned above.

In pursuit of this goal, one potential method involves extending the methodologies used for 1-d ROU processes. However, this strategy appears to be intractable. Indeed, in both Ward and Glynn (2003b) and Xing et al. (2009), the study of the stationary distributions and the positive recurrence relies extensively on the explicit solutions to the corresponding ordinary differential equations (ODEs). When transitioning from the 1-d to the 2-d case, these ODEs transform into partial differential equations (PDEs), whose explicit solutions seem to be impossible to find. Moreover, their proofs use the result that the 1-d ROU processes are regenerative. However, in the two-dimensional case, this property appears no longer to hold, as 2-d processes are unlikely to revisit the same point infinitely often.

Furthermore, it might be tempting to consider the 2-d ROU process in (1.1) as a special case in the general class of constrained diffusion processes, which has been extensively studied in literature (Budhiraja and Dupuis 1999, Atar et al. 2001, Budhiraja et al. 2014, Kang and Ramanan 2014). However, these results cannot be directly applied to our 2-d ROU processes. For example, whereas the positive recurrence of constrained diffusions is well studied in Atar et al. (2001), their outcomes require strict constraints on the drift. As a result, they can only demonstrate the positive recurrence of our 2-d ROU processes when both ${\mu }_1$ and ${\mu }_2$ are $\le 0$ (see Atar et al. 2001, condition 2.3). Nonetheless, a 2-d ROU process with positive ${\mu }_1,{\mu }_2$ does arise in practical applications. Specifically, if a sequence of parallel single-server queues has traffic intensities $({\rho }_1^{(n)},{\rho }_2^{(n)})$ such that the limits of $\sqrt{n}({\rho }_1^{(n)}-1)$ and $\sqrt{n}({\rho }_2^{(n)}-1)$ exist and are positive; then, such a 2-d ROU process emerges as the diffusion scaling limit. Although this result was not explicitly stated in Li and Pang (2024), the derivation of the 2-d ROU process does not impose any constraints on these limits. On the other hand, it is important to explore the positive recurrence of 2-d ROU processes when ${\mu }_1,{\mu }_2\ge 0$. Indeed, when ${\mu }_1,{\mu }_2\ge 0$, the 2-d ROU process $(X_t^1,Y_t^1)$ “dominates” the 2-d ROU process $(X_t^2,Y_t^2)$ when ${\mu }_1,{\mu }_2\le 0$, meaning $X_t^1\ge X_t^2$ and $Y_t^1\ge Y_t^2$ for all $t\ge 0$. Therefore, the positive recurrence of $(X_t^1,Y_t^1)$ immediately implies that of $(X_t^2,Y_t^2)$. Hence, it is sufficient to prove the positive recurrence of $(X_t^1,Y_t^1)$ with ${\mu }_1,{\mu }_2\ge 0$. A detailed explanation will be provided in Section 2. In this paper, we will establish the positive recurrence of the 2-d ROU process in all cases, thereby slightly extending the corresponding results in Atar et al. (2001) as applied to this process.

A common tool for studying the positive recurrence of reflected diffusion processes involves the use of the so-called “fluid path” or “fluid model,” which is a solution to the ODEs obtained by eliminating the randomness from the stochastic differential equations (SDEs) that define the reflected diffusion processes. It has been shown that the “fluid path” tending to the origin generally implies the recurrence of the reflected diffusion process. For more details, refer to Dupuis and Williams (1994), Atar et al. (2001), and Bramson et al. (2010). However, this tool is not applicable to our 2-d ROU processes when ${\mu }_1$ and ${\mu }_2$ are both positive. In fact, solving the following ODEs:

subject to the initial conditions $X_0=x_0\ge 0$ and $Y_0=y_0\ge 0$, yields the “fluid path”:which tends to $({\mu }_1,{\mu }_2)\ne (0,0)$ as $t\to \infty$. As a result, we need a new approach to study the stability of our processes. Using a method that shares the same flavor with those in Hobson and Rogers (1993) and Ernst et al. (2021), we have succeeded in demonstrating the positive recurrence of the 2-d ROU processes, regardless of the values of ${\mu }_1$ and ${\mu }_2$. Indeed, we will prove that the process reaches an arbitrary neighborhood of the origin within a finite mean time. The strategy is summarized as follows:

• We first construct a discrete-time jump process $(X_{S_n},Y_{S_n})$ from the original 2-d ROU process by restricting it to a sequence of elaborately selected stopping times such that one of $X_{S_n}$ and $Y_{S_n}$ is zero for $n\ge 1$.

• We then demonstrate that the above jump process reaches a given neighborhood of the origin within a finite mean time when starting from a position outside that neighborhood. Consequently, the original process exhibits the same property.

• Finally, by focusing on the distance between $(X_t,Y_t)$ and the origin and using the aforementioned property, we derive a similar regenerative process, enabling us to establish the positive recurrence.

It is worth mentioning that we believe our method, in conjunction with the Skorokhod maps for 1-d ROU processes, may be extended to address the case of oblique reflections. We have a strong inclination that 2-d ROU processes are always positive recurrent, irrespective of the reflection directions. We leave this question for future investigation.

The second contribution of this paper lies in studying the exponential rates of convergence to the stationary distribution. The existence and the uniqueness of the stationary distribution, which are the prerequisites for our study, have been previously explored in Kinnally and Williams (2010). Furthermore, our results shall provide an alternative proof of these properties. Following a similar argument as in (Khasminskii 2011, section 4.4), the existence can be established through the positive recurrence, whereas our subsequent results on convergence rates ensure the uniqueness.

Inspired by Lund et al. (1996) and Sarantsev (2020), we establish the convergence rates of the transition probability measures to the invariant measure in both total variation distance and Wasserstein distance. The derivation of the convergence rate in the total variation distance relies on the tail estimate of the first hitting time for the Ornstein–Uhlenbeck process (see Lemma 13), which affects the exponent involved in the convergence rate. On the other hand, for the Wasserstein distance, we establish a concise bound for the Wasserstein distance between two transition probability measures, resulting in an explicit exponent being $\min \{{\theta }_1,{\theta }_2\}$, governing the rate of convergence to the invariant measure.

We next investigate the stationary distribution of the 2-d ROU process by a numerical scheme based on a Monte Carlo method. For diffusions, two common approaches for investigating the stationary distribution are the PDE method (refer to, e.g., Kurtz 1991, Dai and Kurtz 1994, Williams 1995) and the Monte Carlo method (see, e.g., Lamberton and Pages 2002, Talay 2002, Graham and Talay 2013). However, explicitly solving the PDEs derived from our process appears to be impossible. Therefore, we only focus on the Monte Carlo method. For reflected diffusions in polyhedral domains, Budhiraja et al. (2014) studied the numerical computation of their stationary distributions using an Euler scheme. However, their results cannot be directly applied to our case, as Budhiraja et al. (2014, conditions 1.3 and 1.4) are not satisfied. Therefore, utilizing the unique characteristics of our process, we develop our own numerical scheme for approximating the stationary distribution and also analyze the convergence of our scheme to the stationary distribution. This study is also in a similar spirit to the recent work in Jin et al. (2024a, b) on numerical methods using Euler–Maruyama schemes to compute the stationary distributions of limiting diffusions (without reflections) arising in approximations of queueing networks in heavy traffic. Similar to Budhiraja et al. (2014) and Jin et al. (2024a), our numerical scheme takes a decreasing step size.

Finally, it is worth mentioning that unlike the extensively studied ergodicity properties and stationary distributions of two-dimensional reflected Brownian motions (see, for instance, the recent work in El Kharroubp et al. 2000 and Franceschi and Raschel 2019), such studies on 2-d ROU processes seem to be lacking. Whereas the 2-d ROU processes investigated in this paper are not the most general, we contribute to this area by examining those properties for a specific case. In the future, we aim to extend our results to the most general ROU processes.

1.1. Overview of Our Contributions

Our contributions are as follows:

• Prove the positive recurrence of the 2-d ROU process defined in (1.1); see Theorem 12.

• Establish the exponential rates of convergence to its stationary distribution; see Theorems 15 and 17.

• Provide a numerical scheme to approximate its stationary distribution; see Theorem 18.

1.2. Organization of the Paper

We first summarize the notation used in this paper at the end of this section. Section 2 is devoted to investigating the positive recurrence of the 2-d ROU process defined in (1.1), utilizing the aforementioned strategy. In Section 3, we examine the exponential rates of convergence to the stationary distribution in both total variation distance and Wasserstein distance. Finally, Section 4 focuses on analyzing a numerical scheme based on the Monte Carlo method to approximate the stationary distribution of the 2-d ROU process.

1.3. Notation

Let $(\mathrm{\Omega },\mathcal{F},{\{{\mathcal{F}}_t\}}_{t\ge 0},P)$ be the natural filtration space generated by the two-dimensional Brownian motion $(W_t,B_t)$. Because $(X_t,Y_t)$ are the solutions to the SDE (1.1), they are also processes on $(\mathrm{\Omega },\mathcal{F},{\{{\mathcal{F}}_t\}}_{t\ge 0},P)$. Throughout the paper, ${\mathbb{R}}_{+}^2≔\{(x,y):x,y\ge 0\}$ denotes the two-dimensional positive orthant, and $\partial {\mathbb{R}}_{+}^2$ represents its boundary. Additionally, we use ${\mathcal{O}}_B$ to denote the set $\{(x,y)\in {\mathbb{R}}_{+}^2:x+y\le B\}$. For $(x,y)\in {\mathbb{R}}_{+}^2$, we use $P^{(x,y)}$ and $E^{(x,y)}$ to denote the probability and the expectation, respectively, when the 2-d ROU process starts from $(x,y)$. Similarly, for a 1-d process starting from $z\in \mathbb{R}$, $P^z$ and $E^z$ are the corresponding probability and expectation. Moreover, for $x,y\in \mathbb{R}$, we define $x\wedge y≔\min \{x,y\}$, $x\vee y≔\max \{x,y\}$, and $x^{+}≔\max \{x,0\}$. Additionally, ${\delta }_{(x,y)}$ is the Dirac measure on ${\mathbb{R}}^2$ concentrated at $(x,y)$. For a set A, we use ${\mathrm{𝟙}}_A$ to represent the indicator function of A.

When studying the convergence rates to the stationary distribution, we will need the following notation. We denote the stationary distribution as $\pi$. Furthermore, for $(x,y)\in {\mathbb{R}}_{+}^2$, $P_t((x,y),·)$ represents the transition probability measure for the 2-d ROU process. Finally, we introduce total variation distance and Wasserstein distance. Given two measures $m_1$ and $m_2$ on ${\mathbb{R}}_{+}^2$, the total variation distance is defined by

where ${\mathcal{B}}_{{\mathbb{R}}_{+}^2}$ is the family of all Borel sets on ${\mathbb{R}}_{+}^2$. Further, if both $\int \Vert x\Vert m_1(dx)$ and $\int \Vert x\Vert m_2(dx)$ are finite, we define the Wasserstein distance between $m_1$ and $m_2$ bywhere $\text{Lip}(1)$ is the set of all Lipschitz functions with Lipschitz constant $\le 1$.

Additional notation is introduced in the paper as needed.

2. Positive Recurrence

The aim of this section is to study the positive recurrence of the 2-d ROU process defined in (1.1), following the strategy outlined in Section 1.

2.1. Preliminary

We start by citing a lemma from Ward and Glynn (2003b), which demonstrates that a 1-d ROU process is “controlled” or “dominated” by a 1-d reflected Brownian motion (RBM) with the same parameters. This lemma allows us to exploit the properties of the 1-d RBM to study the ROU process. Subsequently, we extend this lemma to compare between two 1-d ROU processes, and the result is summarized in Lemma 2.

Lemma 1.

Let X be a 1-d ROU process satisfying

$$d{X}_t=\theta (\mu -X_t)dt+\sigma d{W}_t+d{L}_t^X$$

starting from $x_0\ge 0$. Also, let $\tilde{X}$ be a 1-d RBM associated with X, satisfying

$$d{\tilde{X}}_t=\theta \mu \hspace{0.17em}dt+\sigma d{W}_t+d{L}_t^{\tilde{X}}$$

starting from $x_0\ge 0$. Then, for each $t\ge 0$,

$$X_t\le {\tilde{X}}_t.$$

Proof.

See Ward and Glynn (2003b, proof of roposition 2). □

Lemma 2.

Let X be a 1-d ROU process satisfying

$$d{X}_t=\theta (\mu -X_t)dt+\sigma d{W}_t+d{L}_t^X$$

starting from $x_0\ge 0$, where $L^X$ represents the local time of X at zero. Furthermore, let $\tilde{X}$ be another 1-d ROU process driven by the same Brownian motion $W_t$, satisfying

$$d{\tilde{X}}_t=\tilde{\theta }(\tilde{\mu }-{\tilde{X}}_t)dt+\sigma d{W}_t+d{L}_t^{\tilde{X}}$$

starting from ${\tilde{x}}_0\ge 0$, where $L^{\tilde{X}}$ represents the local time of $\tilde{X}$ at zero. Suppose that ${\tilde{x}}_0\ge x_0$, $\tilde{\theta }\le \theta$, and $\tilde{\theta }\tilde{\mu }\ge \theta \mu$. Then, for each $t\ge 0$,

$$X_t\le {\tilde{X}}_t.$$

Proof.

Suppose that there exists $t>0$ for which $X_t>{\tilde{X}}_t$. By assumptions, $X_0\le {\tilde{X}}_0$. Furthermore, noting that $X-\tilde{X}$ has continuous sample paths, there exists $s\in [0,t)$ such that $X_s={\tilde{X}}_s$, but $X_u>{\tilde{X}}_u$ for $s<u\le t$. Hence, $\theta X_u\ge \tilde{\theta }{\tilde{X}}_u$ for $s<u\le t$. Note that

$$\begin{array}{l}{X}_t-{\tilde{X}}_t=\left(X_s+{\displaystyle {\int }_s^t\theta }(\mu -X_u)\hspace{0.17em}du+{\displaystyle {\int }_s^t\sigma }\hspace{0.17em}d{W}_u+{\displaystyle {\int }_s^{t}d}{L}_u^X\right)-\left({\tilde{X}}_s+{\displaystyle {\int }_s^t\tilde{\theta }}\left(\tilde{\mu }-{\tilde{X}}_u\right)\hspace{0.17em}du+{\displaystyle {\int }_s^t\sigma }\hspace{0.17em}d{W}_u+{\displaystyle {\int }_s^{t}d}{L}_u^{\tilde{X}}\right)\hfill \\ =\hspace{0.17em}{X}_s-{\tilde{X}}_s-\left(\tilde{\theta }\tilde{\mu }-\theta \mu \right)(t-s)-{\displaystyle {\int }_s^t\left(\theta X_u-\tilde{\theta }{\tilde{X}}_u\right)}\hspace{0.17em}du+\left(L_t^X-L_s^X\right)-\left(L_t^{\tilde{X}}-L_s^{\tilde{X}}\right)\hfill \\ \le \hspace{0.17em}\left(L_t^X-L_s^X\right)-\left(L_t^{\tilde{X}}-L_s^{\tilde{X}}\right).\hfill \end{array}$$

Because $L^{\tilde{X}}$ is nondecreasing, it follows

$$X_t-{\tilde{X}}_t\le L_t^X-L_s^X.$$

Noting that $X_u>{\tilde{X}}_u\ge 0$ for $u\in (s,t]$ and $L^X$ is a continuous process that increases only when X equals zero, we conclude $L_t^X=L_s^X$. Therefore, $X_t\le {\tilde{X}}_t$, which is a contradiction. □

Remark 3.

By using a similar argument, the result in Lemma 2 still holds if the 1-d ROU processes are replaces by two 1-d Ornstein–Uhlenbeck processes. Furthermore, we can also prove a 1-d Ornstein–Uhlenbeck process is bounded above by a 1-d ROU process with the same parameters.

With the above lemma in hand, we are ready to show that it is sufficient to prove the positive recurrence when both ${\mu }_1$ and ${\mu }_2$ are nonnegative, as promised in Section 1. Indeed, if at least one of ${\mu }_1$ and ${\mu }_2$ is negative, we can construct a new 2-d ROU process $(\overline{X},\overline{Y})$ by the following SDEs:

subject to ${\overline{X}}_0=x_0$ and ${\overline{Y}}_0=y_0$. Suppose $(X,Y)$ is the 2-d ROU derived from (1.1). Applying Lemma 2, we obtain that $(X,Y)$ is “controlled” by $(\overline{X},\overline{Y})$; that is, for any $t\ge 0$, $X_t\le {\overline{X}}_t$, and $Y_t\le {\overline{Y}}_t$. Then, for any neighborhood of the origin, the hitting time of $(X,Y)$ to that neighborhood is always bounded by the hitting time of $(\overline{X},\overline{Y})$. Hence, the positive recurrence of $(\overline{X},\overline{Y})$ implies the positive recurrence of $(X,Y)$. As a result, in the rest of the section, we will always assume both ${\mu }_1$ and ${\mu }_2$ are nonnegative.

With the above preparation, we now return to the proof of the positive recurrence. For any neighborhood $\mathcal{N}$ of the origin, we define

In what follows, we will prove that $E^{(x,y)}[T]<\infty$ for any $(x,y)\in {\mathbb{R}}_{+}^2$. Consequently, the positive recurrence follows. Here, the symbol $E^{(x,y)}$ represents the expectation when the 2-d ROU process starts from $(x,y)$. The proof will be divided into three subsections, aligning with the three steps described in Section 1.

2.2. Stopping Times

In this subsection, we will first introduce certain stopping times and estimate their expectations. These will serve as the basis for defining a sequence ${\{S_n\}}_{n=1}^{\infty }$ of stopping times, enabling us to construct the jump process. Moreover, we will provide estimates for the expectation of $S_1$, as well as the position of $(X,Y)$ at $S_1$.

Define

Indeed, ${\tau }^X$ and ${\tau }^Y$ are the first hitting times of X and Y to zero, respectively. If X starts at zero, we define

and

Thus, ${\delta }^X$ is the first time that X returns to zero after it hits one. Additionally, ${\delta }^X$ can also be written as

Here, $\theta (S)$ is a shift operator, the effect of which on a path $\omega$ is to cut off the part of the path before $S(\omega )$ and to shift the remaining part in time. Similarly, if Y starts at zero, we define

and

Again, ${\delta }^Y$ is the first time that Y returns to zero after it hits one, and ${\delta }^Y$ can be written as

Once these stopping times are defined, we can proceed to provide estimates for their expectations, which are summarized in Lemma 4, Lemma 5, and Lemma 6.

Lemma 4.

For any $(x,y)\in {\mathbb{R}}_{+}^2$, we have

$$E^{(x,y)}\left[{\tau }^X\right]=\frac{2}{{\sigma }_1^2}{\displaystyle {\int }_0^x{e}^{\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}}{\displaystyle {\int }_v^{\infty }{e}^{-\frac{{\theta }_1{(u-{\mu }_1)}^2}{{\sigma }_1^2}}}\hspace{0.17em}\mathit{\text{dudv}},$$(2.1)

$$E^{(x,y)}\left[{\tau }^Y\right]=\frac{2}{{\sigma }_2^2}{\displaystyle {\int }_0^y{e}^{\frac{{\theta }_2{(v-{\mu }_2)}^2}{{\sigma }_2^2}}}{\displaystyle {\int }_v^{\infty }{e}^{-\frac{{\theta }_2{(u-{\mu }_2)}^2}{{\sigma }_2^2}}}\hspace{0.17em}\mathit{\text{dudv}}.$$(2.2)

Furthermore,

$$E^{(x,y)}\left[{\tau }^X\right]\le C_1+\frac{1}{{\theta }_1}\log (x+1),$$(2.3)

$$E^{(x,y)}\left[{\tau }^Y\right]\le C_2+\frac{1}{{\theta }_2}\log (y+1),$$(2.4)

where

$$\begin{array}{l}{C}_1≔\frac{2}{{\sigma }_1}\sqrt{\frac{\pi }{{\theta }_1}}{\displaystyle {\int }_0^{{\mu }_1+1}{e}^{\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}}\hspace{0.17em}dv,\hfill \\ C_2≔\frac{2}{{\sigma }_2}\sqrt{\frac{\pi }{{\theta }_2}}{\displaystyle {\int }_0^{{\mu }_2+1}{e}^{\frac{{\theta }_2{(v-{\mu }_2)}^2}{{\sigma }_2^2}}}\hspace{0.17em}dv.\hfill \end{array}$$

Proof.

The first statement is an immediate consequence of Ward and Glynn (2003b, proposition 4). We then turn to the second statement. By the symmetry, we only need to derive (2.3) from (2.1).

By substituting the variable $z=\sqrt{2{\theta }_1/{\sigma }_1^2}(u-{\mu }_1)$, we have

$${\displaystyle {\int }_v^{\infty }{e}^{-\frac{{\theta }_1{(u-{\mu }_1)}^2}{{\sigma }_1^2}}}\hspace{0.17em}du=\sqrt{\frac{{\sigma }_1^2}{2{\theta }_1}}{\displaystyle {\int }_{\sqrt{\frac{2{\theta }_1}{{\sigma }_1^2}}(v-{\mu }_1)}^{\infty }{e}^{-\frac{z^2}{2}}}\hspace{0.17em}dz.$$

If $v\le {\mu }_1+1$,

$${\displaystyle {\int }_v^{\infty }{e}^{-\frac{{\theta }_1{(u-{\mu }_1)}^2}{{\sigma }_1^2}}}\hspace{0.17em}du=\sqrt{\frac{{\sigma }_1^2}{2{\theta }_1}}{\displaystyle {\int }_{\sqrt{\frac{2{\theta }_1}{{\sigma }_1^2}}(v-{\mu }_1)}^{\infty }{e}^{-\frac{z^2}{2}}}\hspace{0.17em}dz\le \sqrt{\frac{{\sigma }_1^2}{2{\theta }_1}}{\displaystyle {\int }_{-\infty }^{\infty }{e}^{-\frac{z^2}{2}}}\hspace{0.17em}dz=\sqrt{\frac{\pi {\sigma }_1^2}{{\theta }_1}}.$$

If $v>{\mu }_1+1$, using the upper bound for the tail of the standard normal distribution,

$$\begin{array}{ll}\hfill & {\displaystyle {\int }_v^{\infty }{e}^{-\frac{{\theta }_1{(u-{\mu }_1)}^2}{{\sigma }_1^2}}}\hspace{0.17em}du=\sqrt{\frac{{\sigma }_1^2}{2{\theta }_1}}{\displaystyle {\int }_{\sqrt{\frac{2{\theta }_1}{{\sigma }_1^2}}(v-{\mu }_1)}^{\infty }{e}^{-\frac{z^2}{2}}}\hspace{0.17em}dz\le \sqrt{\frac{{\sigma }_1^2}{2{\theta }_1}}·\frac{1}{\sqrt{\frac{2{\theta }_1}{{\sigma }_1^2}}(v-{\mu }_1)}{e}^{-\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}=\frac{{\sigma }_1^2}{2{\theta }_1}\frac{1}{v-{\mu }_1}\hspace{0.17em}{e}^{-\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}.\hfill \end{array}$$

When $x\ge {\mu }_1+1$, combining the above two displays together with (2.1), we have

$$\begin{array}{l}{E}^{(x,y)}\left[{\tau }^X\right]=\frac{2}{{\sigma }_1^2}{\displaystyle {\int }_0^x{e}^{\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}}{\displaystyle {\int }_v^{\infty }{e}^{-\frac{{\theta }_1{(u-{\mu }_1)}^2}{{\sigma }_1^2}}}\hspace{0.17em}\mathit{\text{dudv}}\hfill \\ =\frac{2}{{\sigma }_1^2}{\displaystyle {\int }_0^{{\mu }_1+1}{e}^{\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}}{\displaystyle {\int }_v^{\infty }{e}^{-\frac{{\theta }_1{(u-{\mu }_1)}^2}{{\sigma }_1^2}}}\hspace{0.17em}\mathit{\text{dudv}}+\hspace{0.17em}\frac{2}{{\sigma }_1^2}{\displaystyle {\int }_{{\mu }_1+1}^x{e}^{\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}}{\displaystyle {\int }_v^{\infty }{e}^{-\frac{{\theta }_1{(u-{\mu }_1)}^2}{{\sigma }_1^2}}}\hspace{0.17em}\mathit{\text{dudv}}\hfill \\ \le \frac{2}{{\sigma }_1^2}{\displaystyle {\int }_0^{{\mu }_1+1}{e}^{\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}}\sqrt{\frac{\pi {\sigma }_1^2}{{\theta }_1}}\hspace{0.17em}dv+\hspace{0.17em}\frac{2}{{\sigma }_1^2}{\displaystyle {\int }_{{\mu }_1+1}^x{e}^{\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}}\frac{{\sigma }_1^2}{2{\theta }_1}\frac{1}{v-{\mu }_1}\hspace{0.17em}{e}^{-\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}\hspace{0.17em}dv\hfill \\ =C_1+\frac{1}{{\theta }_1}{\displaystyle {\int }_{{\mu }_1+1}^x\frac{1}{v-{\mu }_1}}\hspace{0.17em}dv=C_1+\frac{1}{{\theta }_1}\log (x-{\mu }_1)\hfill \\ \le C_1+\frac{1}{{\theta }_1}\log (x+1).\hfill \end{array}$$

When $x<{\mu }_1+1$, it is immediate that $E^{(x,y)}\left[{\tau }^X\right]\le C_1\le C_1+\frac{1}{{\theta }_1}\log (x+1)$. Combining these two cases, (2.3) follows. □

Lemma 5.

For any $x,y\ge 0$, we have

$$\begin{array}{l}{E}^{(0,y)}\left[{\zeta }^X\right]=\frac{2}{{\sigma }_1^2}{\displaystyle {\int }_0^1{e}^{\frac{{\theta }_1{(v-{\mu }_1)}^2}{{\sigma }_1^2}}}{\displaystyle {\int }_0^v{e}^{-\frac{{\theta }_1{(u-{\mu }_1)}^2}{{\sigma }_1^2}}}\hspace{0.17em}\mathit{\text{dudv}}<\infty ,\hfill \\ E^{(x,0)}\left[{\zeta }^Y\right]=\frac{2}{{\sigma }_2^2}{\displaystyle {\int }_0^1{e}^{\frac{{\theta }_2{(v-{\mu }_2)}^2}{{\sigma }_2^2}}}{\displaystyle {\int }_0^v{e}^{-\frac{{\theta }_2{(u-{\mu }_2)}^2}{{\sigma }_2^2}}}\hspace{0.17em}\mathit{\text{dudv}}<\infty .\hfill \end{array}$$

Thus, we use constants $C_3$ and $C_4$ to denote $E^{(0,y)}[{\zeta }^X]$ and $E^{(x,0)}[{\zeta }^Y]$, respectively.

Proof.

See Ward and Glynn (2003b, theorem 2). □

With the above two lemmas and the fact that ${\delta }^X={\zeta }^X+{\tau }^X°\theta ({\zeta }^X)$ and ${\delta }^Y={\zeta }^Y+{\tau }^Y°\theta ({\zeta }^Y)$, we are ready to establish Lemma 6.

Lemma 6.

For any $x,y\ge 0$, we have

$$E^{(0,y)}\left[{\delta }^X\right]\le C_5,$$(2.5)

$$E^{(x,0)}\left[{\delta }^Y\right]\le C_6,$$(2.6)

where

$$\begin{array}{l}{C}_5≔C_1+C_3+\frac{1}{{\theta }_1}\log (2),\hfill \\ C_6≔C_2+C_4+\frac{1}{{\theta }_2}\log (2).\hfill \end{array}$$

Proof.

We only prove (2.5), and (2.6) follows similarly.

By the strong Markov property,

$$\begin{array}{l}{E}^{(0,y)}\left[{\delta }^X\right]=E^{(0,y)}\left[{\zeta }^X\right]+E^{(0,y)}\left[{\tau }^X°\theta ({\zeta }^X)\right]\hfill \\ =E^{(0,y)}\left[{\zeta }^X\right]+E^{(0,y)}[E^{(1,Y_{{\zeta }^X})}\left[{\tau }^X\right]].\hfill \end{array}$$

By Lemma 5, we have

$$E^{(0,y)}\left[{\zeta }^X\right]\le C_3.$$

By Lemma 4, we have

$$E^{(1,Y_{{\zeta }^X})}\left[{\tau }^X\right]\le C_1+\frac{1}{{\theta }_1}\log (2).$$

The desired result follows by combining the above three displays. □

The three lemmas above provide bounds for the expectations of the stopping times. To define and analyze the jump process mentioned in Section 1, we also need to evaluate the position of $(X,Y)$ at a given stopping time. The following lemma, which provides a bound for the expectation of a one-dimensional Ornstein–Uhlenbeck process at any stopping time, will be a valuable tool for this purpose.

Lemma 7.

Let Z be a one-dimensional Ornstein–Uhlenbeck process satisfying

$$d{Z}_t=\theta (\mu -Z_t)dt+\sigma d{B}_t$$

and starting at z, where $\theta >0$, $\sigma >0$, $\mu \ge 0$ and B is a standard Brownian motion. If R is a stopping time with a finite mean, then

$$E^z\left[|Z_R|\right]\le |z|E^z[e^{-\theta R}]+\mu +\sigma ·{\left\{E^z\left[R\right]\right\}}^{\frac{1}{2}}.$$

Proof.

It is immediate that Z has the expression

$$Z_t=z{e}^{-\theta t}+\mu \left(1-e^{-\theta t}\right)+\sigma e^{-\theta t}{\displaystyle {\int }_0^t{e}^{\theta s}}\hspace{0.17em}d{B}_s.$$

By substituting R for t and applying the triangle inequality, it easily follows that

$$|Z_R|\le |z|e^{-\theta R}+\mu +\sigma e^{-\theta R}\left|{\displaystyle {\int }_0^R{e}^{\theta s}}\hspace{0.17em}d{B}_s\right|.$$(2.7)

We then evaluate the expectation of $e^{-\theta R}|{\int }_0^R{e}^{\theta s}\hspace{0.17em}d{B}_s|$. By Itô’s lemma,

$$d\left(e^{-2\theta t}{\left({\displaystyle {\int }_0^t{e}^{\theta s}}\hspace{0.17em}d{B}_s\right)}^2\right)=-2\theta e^{-2\theta t}{\left({\displaystyle {\int }_0^t{e}^{\theta s}}\hspace{0.17em}d{B}_s\right)}^{2}dt+2e^{-\theta t}\left({\displaystyle {\int }_0^t{e}^{\theta s}}\hspace{0.17em}d{B}_s\right)d{B}_t+dt.$$

Using the optional sampling theorem yields

$$\begin{array}{l}{E}^z\left[e^{-2\theta (t\wedge R)}{\left({\displaystyle {\int }_0^{t\wedge R}{e}^{\theta s}}\hspace{0.17em}d{B}_s\right)}^2\right]=-2\theta E^z\left[{\displaystyle {\int }_0^{t\wedge R}{e}^{-2\theta s}}{\left({\displaystyle {\int }_0^s{e}^{\theta u}}\hspace{0.17em}d{B}_u\right)}^2\hspace{0.17em}ds\right]+E^z\left[t\wedge R\right]\hfill \\ \le E^z\left[t\wedge R\right].\hfill \end{array}$$

Letting $t\to \infty$, together with Fatou’s lemma and monotone convergence theorem, we have

$$E^z\left[e^{-2\theta R}{\left({\displaystyle {\int }_0^R{e}^{\theta s}}\hspace{0.17em}d{B}_s\right)}^2\right]\le E^z\left[R\right].$$

It follows by Jensen’s inequality that

$$E^z\left[e^{-\theta R}\left|{\displaystyle {\int }_0^R{e}^{\theta s}}\hspace{0.17em}d{B}_s\right|\right]\le {\left\{E^z\left[e^{-2\theta R}{\left({\displaystyle {\int }_0^R{e}^{\theta s}}\hspace{0.17em}d{B}_s\right)}^2\right]\right\}}^{\frac{1}{2}}\le {\left\{E^z\left[R\right]\right\}}^{\frac{1}{2}}.$$

By taking expectation on both sides of (2.7) and combining it with the above display, the desired result follows. □

Next, we introduce a new stopping time, S, related to ${\tau }^X,{\tau }^Y,{\delta }^X$, and ${\delta }^Y$. This stopping time is a key ingredient in defining the jump process and is defined as follows:

• If $(X,Y)$ starts at $(x,0)$, then $S={\tau }^X\wedge {\delta }^Y$.

• If $(X,Y)$ starts at $(0,y)$, then $S={\tau }^Y\wedge {\delta }^X$.

In other words, if the process $(X,Y)$ starts from a point on $\partial {\mathbb{R}}_{+}^2\setminus \{(0,0)\}=\{(x,y)\ne (0,0):x=0 \text{or} y=0\}$, the stopping time S means the first time when either one component of the pair $(X_t,Y_t)$ reaches zero or the other that starts at zero reaches one and then returns to zero. In particular, if the process $(X,Y)$ starts from $(0,0)$, both two definitions result in $S=0$.

Before introducing the sequence of stopping times that defines the jump process, we pause to examine the properties of S. Indeed, we will estimate the expectations of S and $X_S+Y_S$, as summarized in Lemma 8 and Lemma 9, respectively.

We begin with Lemma 8.

Lemma 8.

For any $(x,y)\in \partial {\mathbb{R}}_{+}^2$, we always have

$$E^{(x,y)}\left[S\right]\le C_5\vee C_6,$$

where $C_5$ and $C_6$ are defined in Lemma 6.

Proof.

If $(X,Y)$ starts at $(x,0)$, it follows by Lemma 6 and the definition of S that

$$E^{(x,0)}[S]=E^{(x,0)}[{\tau }^X\wedge {\delta }^Y]\le E^{(x,0)}[{\delta }^Y]\le C_6.$$

Similarly, if $(X,Y)$ starts at $(0,y)$, we have $E^{(0,y)}[S]\le C_5$. The desired result follows by combining the above two cases. □

The following lemma concerns the expectation of $X_S+Y_S$.

Lemma 9.

For any $x,y\ge 0$, we have

$$E^{(x,0)}\left[X_S+Y_S\right]\le C_{7}x+{\sigma }_2{\theta }_1^{-1/2}\sqrt{x}+C_9,$$(2.8)

$$E^{(0,y)}\left[X_S+Y_S\right]\le C_{8}y+{\sigma }_1{\theta }_2^{-1/2}\sqrt{y}+C_{10},$$(2.9)

where

$$\begin{array}{l}{C}_7≔E^{(x,0)}\left[e^{-{\theta }_1{\delta }^Y}\right]<1,\hfill \\ C_8≔E^{(0,y)}\left[e^{-{\theta }_2{\delta }^X}\right]<1,\hfill \end{array}$$

and

$$\begin{array}{l}{C}_9≔2+{\mu }_1+{\mu }_2+{\sigma }_1\sqrt{C_6}+{\sigma }_2\sqrt{C_1}+{\sigma }_2{\theta }_1^{-1/2}\sqrt{{\mu }_1+{\sigma }_1\sqrt{C_4}},\hfill \\ C_{10}≔2+{\mu }_1+{\mu }_2+{\sigma }_2\sqrt{C_5}+{\sigma }_1\sqrt{C_2}+{\sigma }_1{\theta }_2^{-1/2}\sqrt{{\mu }_2+{\sigma }_2\sqrt{C_3}}.\hfill \end{array}$$

Proof.

We only prove (2.8), and (2.9) follows similarly.

By the definition of S, when $(X,Y)$ starts at $(x,0)$,

$$X_S=X_{{\tau }^X\wedge {\delta }^Y}=X_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\tau }^X\le {\delta }^Y\}}+X_{{\delta }^Y}·{\mathrm{𝟙}}_{\{{\delta }^Y<{\tau }^X\}}=X_{{\delta }^Y}·{\mathrm{𝟙}}_{\{{\delta }^Y<{\tau }^X\}},$$

where the last equality follows from the fact that $X_{{\tau }^X}=0$, as ${\tau }^X$ is the first time X hits zero. On the event $\{t<{\tau }^X\}$, the 1-d ROU process X is exactly an Ornstein–Uhlenbeck process and, hence, has the representation

$$X_t=x{e}^{-{\theta }_{1}t}+{\mu }_1\left(1-e^{-{\theta }_{1}t}\right)+{\sigma }_1{e}^{-{\theta }_{1}t}{\displaystyle {\int }_0^t{e}^{\theta s}}\hspace{0.17em}d{W}_s.$$

Applying Lemma 7 yields

$$\begin{array}{l}{E}^{(x,0)}\left[X_S\right]=E^{(x,0)}\left[X_{{\delta }^Y}·{\mathrm{𝟙}}_{\{{\delta }^Y<{\tau }^X\}}\right]\\ =E^{(x,0)}\left[\left(x{e}^{-{\theta }_1{\delta }^Y}+{\mu }_1\left(1-e^{-{\theta }_1{\delta }^Y}\right)+{\sigma }_1{e}^{-{\theta }_1{\delta }^Y}{\displaystyle {\int }_0^{{\delta }^Y}{e}^{\theta s}}\hspace{0.17em}d{W}_s\right){\mathrm{𝟙}}_{\{{\delta }^Y<{\tau }^X\}}\right]\\ \le E^{(x,0)}\left[\left|x{e}^{-{\theta }_1{\delta }^Y}+{\mu }_1\left(1-e^{-{\theta }_1{\delta }^Y}\right)+{\sigma }_1{e}^{-{\theta }_1{\delta }^Y}{\displaystyle {\int }_0^{{\delta }^Y}{e}^{\theta s}}\hspace{0.17em}d{W}_s\right|\right]\hfill \\ \le x{E}^{(x,0)}\left[e^{-{\theta }_1{\delta }^Y}\right]+{\mu }_1+{\sigma }_1{\left\{E^{(x,0)}\left[{\delta }^Y\right]\right\}}^{\frac{1}{2}}\hfill \\ \le x{E}^{(x,0)}\left[e^{-{\theta }_1{\delta }^Y}\right]+{\mu }_1+{\sigma }_1\sqrt{C_6},\hfill \end{array}$$(2.10)

where, in the last inequality, we have invoked Lemma 6.

Similarly, when $(X,Y)$ starts at $(x,0)$, noting that $Y_{{\delta }^Y}=0$, we have

$$\begin{array}{ll}{Y}_S\hfill & =Y_{{\tau }^X\wedge {\delta }^Y}=Y_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\tau }^X<{\delta }^Y\}}+Y_{{\delta }^Y}·{\mathrm{𝟙}}_{\{{\delta }^Y\le {\tau }^X\}}=Y_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\tau }^X<{\delta }^Y\}}\hfill \\ \hfill & =Y_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\tau }^X\le {\zeta }^Y\}}+Y_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\zeta }^Y<{\tau }^X<{\delta }^Y\}}.\hfill \end{array}$$

Because $Y_t\le 1$ for $t\le {\zeta }^Y$, we have $Y_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\tau }^X\le {\zeta }^Y\}}\le 1$; thus,

$$Y_S\le 1+Y_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\zeta }^Y<{\tau }^X<{\delta }^Y\}}.$$

Taking expectation yields

$$\begin{array}{ll}{E}^{(x,0)}\left[Y_S\right]\hfill & \le 1+E^{(x,0)}\left[Y_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\zeta }^Y<{\tau }^X<{\delta }^Y\}}\right]\hfill \\ \hfill & =1+E^{(x,0)}\left[{\mathrm{𝟙}}_{\{{\zeta }^Y<{\tau }^X\}}{Y}_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\tau }^X<{\delta }^Y\}}\right]\hfill \\ \hfill & =1+E^{(x,0)}\left[{\mathrm{𝟙}}_{\{{\zeta }^Y<{\tau }^X\}}{E}^{(X_{{\zeta }^Y},1)}\left[Y_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\tau }^X<{\tau }^Y\}}\right]\right],\hfill \end{array}$$(2.11)

where the last equality follows by the strong Markov property. On the event $\{t<{\tau }^Y\}$, the 1-d ROU process Y is exactly an Ornstein–Uhlenbeck process and, hence, has the form

$$Y_t=e^{-{\theta }_{2}t}+{\mu }_2\left(1-e^{-{\theta }_{2}t}\right)+{\sigma }_2{e}^{-{\theta }_{2}t}{\displaystyle {\int }_0^t{e}^{\theta s}}\hspace{0.17em}d{B}_s.$$

Applying Lemma 7 yields

$$\begin{array}{l}{E}^{(X_{{\zeta }^Y},1)}\left[Y_{{\tau }^X}·{\mathrm{𝟙}}_{\{{\tau }^X<{\tau }^Y\}}\right]=E^{(X_{{\zeta }^Y},1)}\left[\left(e^{-{\theta }_2{\tau }^X}+{\mu }_2\left(1-e^{-{\theta }_2{\tau }^X}\right)+{\sigma }_2{e}^{-{\theta }_2{\tau }^X}{\displaystyle {\int }_0^{{\tau }^X}{e}^{\theta s}}\hspace{0.17em}d{B}_s\right){\mathrm{𝟙}}_{\{{\tau }^X<{\tau }^Y\}}\right]\hfill \\ \le E^{(X_{{\zeta }^Y},1)}\left[\left|e^{-{\theta }_2{\tau }^X}+{\mu }_2\left(1-e^{-{\theta }_2{\tau }^X}\right)+{\sigma }_2{e}^{-{\theta }_2{\tau }^X}{\displaystyle {\int }_0^{{\tau }^X}{e}^{\theta s}}\hspace{0.17em}d{B}_s\right|\right]\hfill \\ \le E^{(X_{{\zeta }^Y},1)}\left[e^{-{\theta }_2{\tau }^X}\right]+{\mu }_2+{\sigma }_2{\left\{E^{(X_{{\zeta }^Y},1)}\left[{\tau }^X\right]\right\}}^{\frac{1}{2}}\hfill \\ \le 1+{\mu }_2+{\sigma }_2·{\left(C_1+\frac{1}{{\theta }_1}\log \left(X_{{\zeta }^Y}+1\right)\right)}^{\frac{1}{2}}\hfill \\ \le 1+{\mu }_2+{\sigma }_2\sqrt{C_1}+{\sigma }_2{\theta }_1^{-1/2}{(\log \left(X_{{\zeta }^Y}+1\right))}^{\frac{1}{2}}\hfill \\ \le 1+{\mu }_2+{\sigma }_2\sqrt{C_1}+{\sigma }_2{\theta }_1^{-1/2}·X_{{\zeta }^Y}^{\frac{1}{2}},\hfill \end{array}$$

where the third last inequality is due to Lemma 4, the second last inequality follows by the inequality $\sqrt{a+b}\le \sqrt{a}+\sqrt{b}$, and in the last inequality, the inequality $\log (x+1)\le x$ is applied. Combining the last display with (2.11), after a routine calculation, we have

$$\begin{array}{ll}{E}^{(x,0)}\left[Y_S\right]\hfill & \le 2+{\mu }_2+{\sigma }_2\sqrt{C_1}+{\sigma }_2{\theta }_1^{-1/2}·E^{(x,0)}[{\mathrm{𝟙}}_{\{{\zeta }^Y<{\tau }^X\}}{X}_{{\zeta }^Y}^{\frac{1}{2}}]\hfill \\ \hfill & \le 2+{\mu }_2+{\sigma }_2\sqrt{C_1}+{\sigma }_2{\theta }_1^{-1/2}·{\left\{E^{(x,0)}\left[X_{{\zeta }^Y}{\mathrm{𝟙}}_{\{{\zeta }^Y<{\tau }^X\}}\right]\right\}}^{\frac{1}{2}},\hfill \end{array}$$(2.12)

where in the last inequality, we have invoked Jensen’s inequality. Again, using the argument that $X_t$ is an Ornstein–Uhlenbeck process when $t\le {\tau }^X$ and Lemma 7, we have

$$\begin{array}{l}{E}^{(x,0)}\left[X_{{\zeta }^Y}{\mathrm{𝟙}}_{\{{\zeta }^Y<{\tau }^X\}}\right]\le x{E}^{(x,0)}\left[e^{-{\theta }_1{\zeta }^Y}\right]+{\mu }_1+{\sigma }_1{\left\{E^{(x,0)}\left[{\zeta }^Y\right]\right\}}^{\frac{1}{2}}\hfill \\ \le x+{\mu }_1+{\sigma }_1\sqrt{C_4},\hfill \end{array}$$

where, in the last inequality, Lemma 5 is applied. Plugging the last display into (2.12) yields

$$\begin{array}{ll}{E}^{(x,0)}\left[Y_S\right]\hfill & \le 2+{\mu }_2+{\sigma }_2\sqrt{C_1}+{\sigma }_2{\theta }_1^{-1/2}·{\left(x+{\mu }_1+{\sigma }_1\sqrt{C_4}\right)}^{\frac{1}{2}}\hfill \\ \le 2+{\mu }_2+{\sigma }_2\sqrt{C_1}+{\sigma }_2{\theta }_1^{-1/2}\sqrt{{\mu }_1+{\sigma }_1\sqrt{C_4}}+{\sigma }_2{\theta }_1^{-1/2}\sqrt{x}.\hfill \end{array}$$

Combining the last display with (2.10), the desired result follows. □

With Lemma 9 in hand, We immediately obtain Corollary 10 as follows:

Corollary 10.

There are constants $C_{11}<1$, $C_{12}$ and $C_{13}$ such that for any $(x,y)\in \partial {\mathbb{R}}_{+}^2$,

$$E^{(x,y)}\left[X_S+Y_S\right]\le C_{11}(x+y)+C_{12}\sqrt{x+y}+C_{13}.$$

Consequently, there exist sufficiently large B and $\alpha \in [0,1)$ such that if $(x,y)\in \partial {\mathbb{R}}_{+}^2$ and $x+y\ge B$, then

$$\begin{array}{ll}{E}^{(x,y)}\left[X_S+Y_S\right]\hfill & \le C_{11}(x+y)+C_{12}\sqrt{x+y}+C_{13}\hfill \\ \hfill & \le \alpha (x+y).\hfill \end{array}$$

In the rest of the section, we will use ${\mathcal{O}}_B$ to denote the set $\{(x,y):x+y\le B\}$.

With the above preparation, we are ready to introduce the stopping times and the jump process as mentioned in Section 1.

Define $S_0=0$, $S_1=S$, and

By the strong Markov property and Lemma 8, we have, for $(x,y)\in \partial {\mathbb{R}}_{+}^2$,

Additionally, we define the jump process as follows:

for $n=0,1,2,\dots$. By the definition of $S_n$, it is immediate that $({\xi }_n,{\eta }_n)\in \partial {\mathbb{R}}_{+}^2$. Therefore, $\{({\xi }_n,{\eta }_n):n=1,2,\dots \}$ is a jump process on $\partial {\mathbb{R}}_{+}^2$.

2.3. Reachability of the Jump Process to the Neighborhood of the Origin

In this subsection, we will prove that if the process $\{({\xi }_n,{\eta }_n):n=0,1,2,\dots \}$ starts from a point on $\partial {\mathbb{R}}_{+}^2$, then it will eventually reach $\partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$ within a finite mean time. Consequently, the same property holds for the original 2-d ROU process.

The preceding statements are summarized in the following lemma.

Lemma 11.

Define $N≔\inf \{n:({\xi }_n,{\eta }_n)\in {\mathcal{O}}_B\}$. We have $N<\infty$ almost surely. Furthermore,

$$E^{(x,y)}\left[S_N\right]\le C_5\vee C_6·\frac{x+y}{(1-\alpha )B}.$$(2.13)

Therefore, the 2-d ROU process $(X,Y)$, starting from a point $(x,y)$ on $\partial {\mathbb{R}}_{+}^2$, will eventually reach the region $\partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$ within an expected time of at most $(C_5\vee C_6)(x+y)/((1-\alpha )B)$.

Proof.

We claim $\{{\alpha }^{-(n\wedge N)}({\xi }_{n\wedge N}+{\eta }_{n\wedge N})\}$ is a nonnegative supermartingale with respect to the filtration ${\{{\mathcal{F}}_{S_n}\}}_{n\in N}$. Recall that ${({\mathcal{F}}_t)}_{t\ge 0}$ is the natural filtration generated by $(W,B)$. Therefore, ${\mathcal{F}}_{S_n}$ represents the corresponding $\sigma$-field stopped at $S_n$. It is immediate that ${\xi }_n$ and ${\eta }_n$ are ${\mathcal{F}}_{S_n}$-measurable because ${\xi }_n=X_{S_n}$ and ${\eta }_n=Y_{S_n}$.

We now turn to verify the above claim.

$$\begin{array}{l}{E}^{(x,y)}\left[{\alpha }^{-((n+1)\wedge N)}({\xi }_{(n+1)\wedge N}+{\eta }_{(n+1)\wedge N})|{\mathcal{F}}_{S_n}\right]\hfill \\ =E^{(x,y)}\left[{\alpha }^{-((n+1)\wedge N)}({\xi }_{(n+1)\wedge N}+{\eta }_{(n+1)\wedge N}){\mathrm{𝟙}}_{\{N\le n\}}|{\mathcal{F}}_{S_n}\right]+E^{(x,y)}\left[{\alpha }^{-((n+1)\wedge N)}({\xi }_{(n+1)\wedge N}+{\eta }_{(n+1)\wedge N}){\mathrm{𝟙}}_{\{N>n\}}|{\mathcal{F}}_{S_n}\right]\hfill \\ =E^{(x,y)}\left[{\alpha }^{-(n\wedge N)}\left({\xi }_{n\wedge N}+{\eta }_{n\wedge N}\right){\mathrm{𝟙}}_{\{N\le n\}}|{\mathcal{F}}_{S_n}\right]+\hspace{0.17em}{E}^{(x,y)}\left[{\alpha }^{-(n+1)}\left({\xi }_{n+1}+{\eta }_{n+1}\right){\mathrm{𝟙}}_{\{N>n\}}|{\mathcal{F}}_{S_n}\right]\hfill \\ ={\alpha }^{-(n\wedge N)}\left({\xi }_{n\wedge N}+{\eta }_{n\wedge N}\right){\mathrm{𝟙}}_{\{N\le n\}}+{\mathrm{𝟙}}_{\{N>n\}}·E^{({\xi }_n,{\eta }_n)}\left[{\alpha }^{-(n+1)}\left(X_S+Y_S\right)\right]\hfill \\ \le {\alpha }^{-(n\wedge N)}\left({\xi }_{n\wedge N}+{\eta }_{n\wedge N}\right){\mathrm{𝟙}}_{\{N\le n\}}+{\mathrm{𝟙}}_{\{N>n\}}·{\alpha }^{-n}\left({\xi }_n+{\eta }_n\right)\hfill \\ ={\alpha }^{-(n\wedge N)}\left({\xi }_{n\wedge N}+{\eta }_{n\wedge N}\right),\hfill \end{array}$$

where the second last equality is due to the strong Markov property, and in the last inequality, Corollary 10 is applied. By the standard properties of a supermartingale,

$$\begin{array}{ll}x+y\hfill & =E^{(x,y)}\left[{\xi }_0+{\eta }_0\right]\ge E^{(x,y)}\left[{\alpha }^{-(n\wedge N)}\left({\xi }_{n\wedge N}+{\eta }_{n\wedge N}\right)\right]\hfill \\ \hfill & \ge E^{(x,y)}\left[{\alpha }^{-(n\wedge N)}\left({\xi }_{n\wedge N}+{\eta }_{n\wedge N}\right){\mathrm{𝟙}}_{\{N>n\}}\right]\hfill \\ \hfill & =E^{(x,y)}\left[{\alpha }^{-n}\left({\xi }_n+{\eta }_n\right){\mathrm{𝟙}}_{\{N>n\}}\right]\hfill \\ \hfill & \ge E^{(x,y)}\left[{\alpha }^{-n}·B·{\mathrm{𝟙}}_{\{N>n\}}\right]\hfill \\ \hfill & ={\alpha }^{-n}B·P^{(x,y)}\left(N>n\right),\hfill \end{array}$$

where the last inequality follows by the fact that ${\xi }_n+{\eta }_n>B$ on the event $\{N>n\}$. Thus,

$$P^{(x,y)}\left(N>n\right)\le \frac{x+y}{B}{\alpha }^n.$$

Given that $\alpha \in [0,1)$, it follows by letting $n\to \infty$ that $N<\infty$ almost surely. Therefore, the jump process $({\xi }_n,{\eta }_n)$ inevitably hits the region $\partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$. If we look back at the 2-d ROU process $(X,Y)$, we can conclude that $(X,Y)$ must hit $\partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$ if it starts from a point on $\partial {\mathbb{R}}_{+}^2$.

Note that the first time when $(X,Y)$ hits $\partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$ is less than or equal to $S_N$. In the next step, we estimate the expectation of $S_N$ and, consequently, the expected time for the process $(X,Y)$ to visit $\partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$. Indeed, by Lemma 8 and the strong Markov property,

$$\begin{array}{ll}{E}^{(x,y)}\left[S_N\right]\hfill & =E^{(x,y)}\left[{\displaystyle \sum _{n=0}^{N-1}\left(S_{n+1}-S_n\right)}\right]\hfill \\ \hfill & =E^{(x,y)}\left[{\displaystyle \sum _{n=0}^{\infty }\left(S_{n+1}-S_n\right)} {\mathrm{𝟙}}_{\{N>n\}}\right]\hfill \\ \hfill & ={\displaystyle \sum _{n=0}^{\infty }{E}^{(x,y)}}\left[\left(S_{n+1}-S_n\right){\mathrm{𝟙}}_{\{N>n\}}\right]\hfill \\ \hfill & ={\displaystyle \sum _{n=0}^{\infty }{E}^{(x,y)}}\left[{\mathrm{𝟙}}_{\{N>n\}}·E^{(X_{S_n},Y_{S_n})}\left[S\right]\right]\hfill \\ \hfill & \le {\displaystyle \sum _{n=0}^{\infty }{E}^{(x,y)}}\left[{\mathrm{𝟙}}_{\{N>n\}}·C_5\vee C_6\right]\hfill \\ \hfill & =C_5\vee C_6·{\displaystyle \sum _{n=0}^{\infty }{P}^{(x,y)}}\left(N>n\right)\hfill \\ \hfill & \le C_5\vee C_6·{\displaystyle \sum _{n=0}^{\infty }\frac{x+y}{B}}{\alpha }^n\hfill \\ \hfill & =C_5\vee C_6·\frac{x+y}{(1-\alpha )B}.\hfill \end{array}$$

Therefore, for the 2-d ROU process $(X,Y)$ starting from a point $(x,y)$ on $\partial {\mathbb{R}}_{+}^2$, it eventually reaches the region $\partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$ with an expected time no greater than $(C_5\vee C_6)(x+y)/((1-\alpha )B)$. □

2.4. Positive Recurrence

The purpose of this subsection is to provide the main (theorem Theorem 12) of this section, which focuses on the positive recurrence of the 2-d ROU process $(X,Y)$.

Recall T is the first time the process hits the neighborhood $\mathcal{N}$ of the origin.

Theorem 12.

There exists a constant $C_{14}$ (see (2.19)) such that for any $(x,y)\in {\mathbb{R}}_{+}^2$

$$E^{(x,y)}\left[T\right]\le C_{14}\left(1+x+y\right).$$

Proof.

Define

$$T_0≔\inf \{t:(X_t,Y_t)\in \partial {\mathbb{R}}_{+}^2\}$$

and

$$H≔\inf \{t\ge T_0:(X_t,Y_t)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B\}.$$

H is the first time that the process $(X,Y)$ hits the area $\partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$. We first prove

$$E^{(x,y)}\left[H\right]\le C_{15}(1+x+y),$$(2.14)

where

$$C_{15}≔\frac{C_5\vee C_6}{(1-\alpha )B}·(1+{\mu }_1+{\mu }_2+\left({\sigma }_1+{\sigma }_2\right)\left(\sqrt{C_1}+{\theta }_1^{-1/2}\right))+C_1+\frac{1}{{\theta }_1}.$$

Note that $T_0\le {\tau }^X\wedge {\tau }^Y\le {\tau }^X$; thus, by Lemma 4,

$$E^{(x,y)}\left[T_0\right]\le E^{(x,y)}\left[{\tau }^X\right]\le C_1+\frac{1}{{\theta }_1}\log (x+1).$$(2.15)

To bound $E^{(x,y)}[H]$, it remains to provide an upper bound for $E^{(x,y)}[H-T_0]$, which depends on the estimate for $E^{(x,y)}[X_{T_0}+Y_{T_0}]$. To achieve this, we note that for $t\le T_0$, $(X_t,Y_t)$ behaves like a 2-d Ornstein–Uhlenbeck process. Using a similar argument as in the proof of Lemma 9, we conclude

$$\begin{array}{ll}{E}^{(x,y)}\left[X_{T_0}+Y_{T_0}\right]\hfill & \le x{E}^{(x,y)}\left[e^{-{\theta }_1{T}_0}\right]+{\mu }_1+{\sigma }_1{\left\{E\left[T_0\right]\right\}}^{\frac{1}{2}}+y{E}^{(x,y)}\left[e^{-{\theta }_1{T}_0}\right]+{\mu }_2+{\sigma }_2{\left\{E\left[T_0\right]\right\}}^{\frac{1}{2}}\hfill \\ \hfill & \le x+y+{\mu }_1+{\mu }_2+\left({\sigma }_1+{\sigma }_2\right){\left(C_1+\frac{1}{{\theta }_1}\log (x+1)\right)}^{\frac{1}{2}}.\hfill \end{array}$$(2.16)

Then, by the strong Markov property,

$$\begin{array}{l}{E}^{(x,y)}\left[H-T_0\right]=E^{(x,y)}\left[E^{(X_{T_0},Y_{T_0})}\left[\inf \left\{t\ge 0:(X_t,Y_t)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B\right\}\right]\right]\hfill \\ \le E^{(x,y)}\left[E^{((X_{T_0},Y_{T_0}))}\left[S_N\right]\right]\hfill \\ \le E^{(x,y)}\left[C_5\vee C_6·\frac{X_{T_0}+Y_{T_0}}{(1-\alpha )B}\right]\hfill \\ \le \frac{C_5\vee C_6}{(1-\alpha )B}·\left(x+y+{\mu }_1+{\mu }_2+\left({\sigma }_1+{\sigma }_2\right){\left(C_1+\frac{1}{{\theta }_1}\log (x+1)\right)}^{\frac{1}{2}}\right),\hfill \end{array}$$

where the second last inequality is due to (2.13), and the last inequality follows from (2.16). Combining the last display with (2.15), after a straightforward algebra, (2.14) follows.

We have already proven that the process $(X,Y)$, starting from any point on $\partial {\mathbb{R}}_{+}^2$, will reach $\partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$ within a time H whose mean is bounded by $1+x+y$ (up to a constant). We then define

$$R≔1\wedge T\wedge \inf \{t\ge 0:X_t+Y_t\ge 2B\}.$$

To finish the proof, we need the following result:

There is a constant $ϵ>0$ such that

$$P^{(x,y)}\left(R=T\right)\ge ϵ,$$(2.17)

for any $(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$.

The proof of this result relies on the standard properties of Brownian motion. By Lemma 1, we can construct a 2-d reflected Brownian motion $(\tilde{X},\tilde{Y})$ such that $X_t\le {\tilde{X}}_t$ and $Y_t\le {\tilde{Y}}_t$ for any $t\ge 0$. Define $\tilde{T}$ and $\tilde{R}$ for $(\tilde{X},\tilde{Y})$ similarly to T and R, respectively. By a standard property of Brownian motion, there exists $ϵ>0$ such that for all $(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$, $P^{(x,y)}(\tilde{R}=\tilde{T})\ge ϵ$ (refer to Hobson and Rogers 1993). Thus,

$$P^{(x,y)}\left(R=T\right)\ge P^{(x,y)}(\tilde{R}=\tilde{T})\ge ϵ.$$

We now return to complete the proof. We define

$$\begin{array}{l}{H}_1=H,\hfill \\ R_1=H_1+R°\theta (H_1),\hfill \end{array}$$

and define recursively for $n\ge 2$,

$$\begin{array}{l}{H}_n=R_{n-1}+H°\theta (R_{n-1}),\hfill \\ R_n=H_n+R°\theta (H_n).\hfill \end{array}$$

It follows immediately from the definition of H and R that $(X_{H_n},Y_{H_n})\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$ and $X_{R_n}+Y_{R_n}\le 2B$ for $n\ge 1$. Considering the Estimate (2.14) and applying the strong Markov property, we have, for $n\ge 2$,

$$\begin{array}{ll}{E}^{(x,y)}\left[H_n-R_{n-1}\right]\hfill & =E^{(x,y)}\left[E^{(X_{R_{n-1}},Y_{R_{n-1}})}\left[H\right]\right]\hfill \\ \hfill & \le E^{(x,y)}\left[C_{15}(1+X_{R_{n-1}}+Y_{R_{n-1}})\right]\hfill \\ \hfill & \le C_{15}(1+2B).\hfill \end{array}$$

Thus, noting that $R\le 1$, we obtain, for $n\ge 2$,

$$\begin{array}{ll}{E}^{(x,y)}\left[H_n-H_{n-1}\right]\hfill & =E^{(x,y)}\left[H_n-R_{n-1}\right]+E^{(x,y)}\left[R_{n-1}-H_{n-1}\right]\hfill \\ \hfill & \le C_{15}(1+2B)+E^{(x,y)}\left[E^{(X_{H_{n-1}},Y_{H_{n-1}})}\left[R\right]\right]\hfill \\ \hfill & \le C_{15}(1+2B)+1.\hfill \end{array}$$

From the above estimate, we conclude that for $n\ge 1$,

$$H_n-H_1-(C_{15}(1+2B)+1)(n-1)$$

is a supermartingale with respect to the filtration ${\{{\mathcal{F}}_{H_n}\}}_{n\ge 1}$. Once again, ${\mathcal{F}}_{H_n}$ is the $\sigma$-field obtained by stopping ${\{{\mathcal{F}}_t\}}_{t\ge 0}$ at $H_n$ (recall ${\{{\mathcal{F}}_t\}}_{t\ge 0}$ is the natural filtration generated by $(W_t,B_t)$). Define

$$M≔\inf \{n\ge 1:(X_{R_n},Y_{R_n})\in \mathcal{N}\}.$$

Then, $M+1$ is a stopping time with respect to the filtration ${\{{\mathcal{F}}_{H_n}\}}_n$ because $R_{n-1}\le H_n$. Applying the optional sampling theorem yields

$$\begin{array}{ll}{E}^{(x,y)}\left[H_{n\wedge (M+1)}-H_1\right]\hfill & \le E^{(x,y)}\left[(C_{15}(1+2B)+1)(n\wedge (M+1)-1)\right]\hfill \\ \hfill & =(C_{15}(1+2B)+1)E^{(x,y)}\left[(n-1)\wedge M\right].\hfill \end{array}$$

Letting $n\to \infty$ and applying the monotone convergence theorem, we have

$$E^{(x,y)}\left[H_{M+1}-H_1\right]\le (C_{15}(1+2B)+1)E^{(x,y)}\left[M\right].$$

Together with (2.14), we have

$$E^{(x,y)}\left[H_{M+1}\right]\le C_{15}(1+x+y)+(C_{15}(1+2B)+1)E^{(x,y)}\left[M\right].$$

It follows immediately from the definition of M that $T\le R_M\le H_{M+1}$; thus,

$$E^{(x,y)}\left[T\right]\le C_{15}(1+x+y)+(C_{15}(1+2B)+1)E^{(x,y)}\left[M\right].$$(2.18)

What remains is to estimate the expectation of M.

For any $(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$, by the strong Markov property,

$$\begin{array}{ll}{P}^{(x,y)}(M>n)\hfill & =E^{(x,y)}\left[{\mathrm{𝟙}}_{\{M>n\}}\right]=E^{(x,y)}\left[{\mathrm{𝟙}}_{\{R_1\ne H_1+T°\theta (H_1)\}}·{\mathrm{𝟙}}_{\{M>n\}}\right]\hfill \\ \hfill & =E^{(x,y)}\left[{\mathrm{𝟙}}_{\{R_1\ne H_1+T°\theta (H_1)\}}·E^{(X_{H_2},Y_{H_2})}\left[{\mathrm{𝟙}}_{\{M>n-1\}}\right]\right]\hfill \\ \hfill & \le E^{(x,y)}\left[{\mathrm{𝟙}}_{\{R_1\ne H_1+T°\theta (H_1)\}}·\underset{(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B}{\sup }{P}^{(x,y)}(M>n-1)\right]\hfill \\ \hfill & =E^{(x,y)}\left[{\mathrm{𝟙}}_{\{R_1\ne H_1+T°\theta (H_1)\}}\right]·\underset{(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B}{\sup }{P}^{(x,y)}(M>n-1)\hfill \\ \hfill & =E^{(x,y)}\left[E^{(X_{H_1},Y_{H_1})}\left[{\mathrm{𝟙}}_{\{R\ne T\}}\right]\right]·\underset{(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B}{\sup }{P}^{(x,y)}(M>n-1)\hfill \\ \hfill & \le (1-ϵ)·\underset{(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B}{\sup }{P}^{(x,y)}(M>n-1),\hfill \end{array}$$

where the last inequality follows by (2.17). Taking the supremum over $(x,y)$ gives

$$\underset{(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B}{\sup }{P}^{(x,y)}(M>n)\le (1-ϵ)·\underset{(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B}{\sup }{P}^{(x,y)}(M>n-1).$$

Continuing the above procedure recursively, we have

$$\underset{(x,y)\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B}{\sup }{P}^{(x,y)}(M>n)\le {(1-ϵ)}^n.$$

Because $(X_{H_1},Y_{H_1})\in \partial {\mathbb{R}}_{+}^2\cap {\mathcal{O}}_B$,

$$\begin{array}{ll}{E}^{(x,y)}\left[M\right]\hfill & =E^{(x,y)}\left[E^{(X_{H_1},Y_{H_1})}[M]\right]=E^{(x,y)}\left[{\displaystyle \sum _{n=0}^{\infty }{P}^{(X_{H_1},Y_{H_1})}}(M>n)\right]\hfill \\ \hfill & \le E^{(x,y)}\left[{\displaystyle \sum _{n=0}^{\infty }{(1-ϵ)}^n}\right]=\frac{1}{ϵ}.\hfill \end{array}$$

Plugging the above result into (2.18) gives

$$E^{(x,y)}\left[T\right]\le C_{15}(1+x+y)+(C_{15}(1+2B)+1)/ϵ\le C_{14}(1+x+y),$$

where

$$C_{14}≔C_{15}+\frac{C_{15}(1+2B)+1}{ϵ}.$$(2.19)

We conclude the proof. □

3. Exponential Rate of Convergence

The section aims to investigate the rate of convergence to the stationary distribution of the 2-d ROU process. We will provide exponential rates of convergence in both total variation distance and Wasserstein distance. Further, we will explicitly specify the exponent for the latter.

3.1. A Tail Estimate for the First Hitting Time of an OU Process

Before addressing the rate of convergence, we first present a (lemma Lemma 13) regarding the tail estimate of the first hitting time for an Ornstein–Uhlenbeck process, which is the main ingredient in studying the rate of convergence in total variation distance.

Lemma 13.

Let Z be a one-dimensional Ornstein–Uhlenbeck process driven by the SDE

$$d{Z}_t=\theta (\mu -Z_t)dt+\sigma d{B}_t$$

starting at z, where $\theta >0$, $\sigma >0$, and B is a standard Brownian motion. Let

$${\tau }^Z≔\inf \{t\ge 0:Z_t=0\}.$$

Then, for any $t_{*}>0$, there exist constants $C_{16}$ and $C_{17}$ depending on $\theta ,\mu ,\sigma$, and $t_{*}$ such that for $t\ge t_{*}$,

$$P^z({\tau }^Z>t)\le C_{16}{e}^{-C_{17}t}+\frac{4}{\sqrt{\pi }}\frac{z\sigma }{{\theta }^2}{e}^{z\theta /{\sigma }^2}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }^2}{4{\sigma }^2}t}.$$

Proof.

We start by considering the case where $\mu \ge 0$. Our approach to estimating $P^z({\tau }^Z>t)$ will be divided into three cases based on the value of z: (i) $z=\mu +1$, (ii) $0\le z<\mu +1$, and (iii) $z>\mu +1$.

1. $z=\mu +1$. Using Alili et al. (2005, corollary 3.2) (after an appropriate transformation), we have that for any $t_{*}>0$, there exist constants $C_{16}$ and $C_{17}$ depending on $\theta ,\mu ,\sigma$, and $t_{*}$ such that for $t\ge t_{*}/2$,

   $$P^{(\mu +1)}({\tau }^Z>t)\le C_{16}{e}^{-2C_{17}t}.$$

2. $0\le z<\mu +1$. By applying Remark 3, Z is bounded above by the Ornstein–Uhlenbeck process with the same parameters, but starting at $\mu +1$; hence, for $t\ge t_{*}/2$,

   $$P^z({\tau }^Z>t)\le P^{(\mu +1)}({\tau }^Z>t)\le C_{16}{e}^{-2C_{17}t}.$$

3. $z>\mu +1$. We define

   $$\kappa ≔\inf \{t\ge 0:Z_t=\mu +1\}.$$

Thus, in the event $\{t\le \kappa \}$, we have $Z_t\ge \mu +1$. Let $\tilde{Z}$ be a Brownian motion with a drift driven by the SDE:

$$d{\tilde{Z}}_t=-\theta dt+\sigma d{B}_t,$$

where it starts at z. Analogously, define

$$\tilde{\kappa }≔\inf \{t\ge 0:{\tilde{Z}}_t=\mu +1\}.$$

By a similar argument to that in the proof of Lemma 1, we have $Z_{t\wedge \kappa }\le {\tilde{Z}}_{t\wedge \kappa }$ for $t\ge 0$. Thus, $\kappa \le \tilde{\kappa }$. Then, we have

$$P^z(\kappa >t)\le P^z(\tilde{\kappa }>t)={\displaystyle {\int }_t^{\infty }\frac{z-(\mu +1)}{\sigma }}\frac{\exp \left(-\frac{{(z-(\mu +1)-\theta s)}^2}{2{\sigma }^{2}s}\right)}{\sqrt{2\pi s^3}}\hspace{0.17em}ds,$$

where the last equality follows by the standard result of the first hitting time for Brownian motion with a drift (cf. Le Gall 2016, section 5.6). A routine calculation gives

$$\begin{array}{l}{\displaystyle {\int }_t^{\infty }\frac{z-(\mu +1)}{\sigma }}\frac{\exp \left(-\frac{{(z-(\mu +1)-\theta s)}^2}{2{\sigma }^{2}s}\right)}{\sqrt{2\pi s^3}}\hspace{0.17em}ds\le \frac{1}{\sqrt{2\pi }}\frac{z-(\mu +1)}{\sigma }{e}^{\frac{\theta (z-(\mu +1))}{{\sigma }^2}}{\displaystyle {\int }_t^{\infty }\frac{e^{-\frac{{\theta }^2}{2{\sigma }^2}s}}{\sqrt{s^3}}}\hspace{0.17em}ds\hfill \\ \le \frac{1}{\sqrt{2\pi }}\frac{z}{\sigma }\hspace{0.17em}{e}^{\frac{\theta z}{{\sigma }^2}}{\displaystyle {\int }_t^{\infty }\frac{e^{-\frac{{\theta }^2}{2{\sigma }^2}s}}{\sqrt{s^3}}}\hspace{0.17em}ds\le \frac{1}{\sqrt{2\pi }}\frac{z}{\sigma }\hspace{0.17em}{e}^{\frac{\theta z}{{\sigma }^2}}\frac{1}{\sqrt{t^3}}\frac{2{\sigma }^2}{{\theta }^2}\hspace{0.17em}{e}^{-\frac{{\theta }^2}{2{\sigma }^2}t}\hfill \\ =\sqrt{\frac{2}{\pi }}\frac{z\sigma }{{\theta }^2}{e}^{z\theta /{\sigma }^2}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }^2}{2{\sigma }^2}t}.\hfill \end{array}$$

Then, together with the strong Markov property, we have, for $t\ge t_{*}$,

$$\begin{array}{ll}{P}^z({\tau }^Z>t)\hfill & =P^z(\kappa +{\tau }^Z°\theta (\kappa )>t)\hfill \\ \hfill & \le P^z(\kappa >t/2)+P^z({\tau }^Z°\theta (\kappa )>t/2)\hfill \\ \hfill & \le \frac{4}{\sqrt{\pi }}\frac{z\sigma }{{\theta }^2}{e}^{z\theta /{\sigma }^2}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }^2}{4{\sigma }^2}t}+E^z\left[E^{Z_{\kappa }}\left[{\mathrm{𝟙}}_{\{{\tau }^Z>t/2\}}\right]\right]\hfill \\ \hfill & \le \frac{4}{\sqrt{\pi }}\frac{z\sigma }{{\theta }^2}{e}^{z\theta /{\sigma }^2}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }^2}{4{\sigma }^2}t}+C_{16}{e}^{-C_{17}t},\hfill \end{array}$$

where, in the last inequality, we have invoked the result from case (i).

Combining the above three cases gives, for $z\ge 0$ and $t\ge t_{*}$,

$$P^z({\tau }^Z>t)\le C_{16}{e}^{-C_{17}t}+\frac{4}{\sqrt{\pi }}\frac{z\sigma }{{\theta }^2}{e}^{z\theta /{\sigma }^2}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }^2}{4{\sigma }^2}t}.$$

What remains is the case that $\mu <0$. Indeed, the desired result follows from the fact that $Z_t$ is bounded by the Ornstein–Uhlenbeck process with $\mu =0$. We conclude the proof. □

3.2. Exponential Rate of Convergence in Total Variation Distance

This subsection is devoted to the rate of convergence in total variation distance. It begins by studying the total variation distance between two transition probability measures, followed by the total variation distance between the transition probability measure and the invariant measure.

Recall that $P_t((x,y),·)$ is the transition probability measure for the process $(X_t,Y_t)$ when it starts from $(x,y)$, and $\pi$ represents the stationary distribution.

Theorem 14.

Fix $(x_1,y_1),(x_2,y_2)\in {\mathbb{R}}_{+}^2$. For any $t_{*}>0$, there exist constants $C_{18},C_{19}$, depending on ${\theta }_1,{\mu }_1,{\sigma }_1$, and $t_{*}$, and $C_{20}$, $C_{21}$, depending on ${\theta }_2,{\mu }_2,{\sigma }_2$, and $t_{*}$ such that for $t\ge t_{*}$,

$$\begin{array}{l}{d}_{TV}(P_t((x_1,y_1),·),P_t((x_2,y_2),·))\le C_{18}{e}^{-C_{19}t}+C_{20}{e}^{-C_{21}t}+\frac{8}{\sqrt{\pi }}\frac{{\sigma }_1\max (x_1,x_2)}{{\theta }_1^2}\hspace{0.17em}{e}^{\frac{{\theta }_1\max (x_1,x_2)}{{\sigma }_1^2}}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }_1^2}{4{\sigma }_1}t}\hfill \\ \hspace{1em}+\hspace{0.17em}\frac{8}{\sqrt{\pi }}\frac{{\sigma }_2\max (y_1,y_2)}{{\theta }_2^2}\hspace{0.17em}{e}^{\frac{{\theta }_2\max (y_1,y_2)}{{\sigma }_2^2}}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }_2^2}{4{\sigma }_2}t}.\hfill \end{array}$$

Proof.

For $i=1,2$, we use $(X_t^{(i)},Y_t^{(i)})$ to denote the 2-d ROU process satisfying the SDEs

$$\begin{array}{ll}\hfill & d{X}_t^{(i)}={\theta }_1({\mu }_1-X_t^{(i)})dt+{\sigma }_{1}d{W}_t+d{L}_t^{X^{(i)}}\hfill \\ \hfill & d{Y}_t^{(i)}={\theta }_2({\mu }_2-Y_t^{(i)})dt+{\sigma }_{2}d{B}_t+d{L}_t^{Y^{(i)}}\hfill \end{array}$$

and starting from $(x_i,y_i)$. Define for $i=1,2$,

$$\begin{array}{l}{\tau }^{X^{(i)}}≔\inf \{t\ge 0:X_t^{(i)}=0\},\hfill \\ {\tau }^{Y^{(i)}}≔\inf \{t\ge 0:Y_t^{(i)}=0\}.\hfill \end{array}$$

Further, we define

$$\begin{array}{l}{\tau }_{\text{max}}^X≔{\tau }^{X^{(1)}}\vee {\tau }^{X^{(2)}},\hfill \\ {\tau }_{\text{max}}^Y≔{\tau }^{Y^{(1)}}\vee {\tau }^{Y^{(2)}}.\hfill \end{array}$$

We first claim that on the event $\{t\ge {\tau }_{\text{max}}^X\}$, $X_t^{(1)}=X_t^{(2)}$. This will be demonstrated subsequently. Without loss of generality, we assume that $x_1\le x_2$. Applying Lemma 2, we have $X_t^{(1)}\le X_t^{(2)}$ for $t\ge 0$. Thus, ${\tau }^{X^{(1)}}\le {\tau }^{X^{(2)}}$, and ${\tau }_{\text{max}}^X={\tau }^{X^{(2)}}$. Because $X_t^{(1)}$ is nonnegative, then

$$0\le X_{{\tau }^{X^{(2)}}}^{(1)}\le X_{{\tau }^{X^{(2)}}}^{(2)}=0,$$

which implies $X_{{\tau }^{X^{(2)}}}^{(1)}=X_{{\tau }^{X^{(2)}}}^{(2)}=0$. By the strong Markov property, we have $X_t^{(1)}=X_t^{(2)}$ for $t\ge {\tau }^{X^{(2)}}$.

Similarly, in the event $\{t\ge {\tau }_{\text{max}}^Y\}$, we have $Y_t^{(1)}=Y_t^{(2)}$.

With the above preparation, we proceed to study the total variation distance between $P_t((x_1,y_1),·)$ and $P_t((x_2,y_2),·)$. It follows from the definition of total variation distance that

$$\begin{array}{l}{d}_{TV}(P_t((x_1,y_1),·),P_t((x_2,y_2),·))=\underset{A\in {\mathcal{B}}_{{\mathbb{R}}_{+}^2}}{\sup }|P_t((x_1,y_1),A)-P_t((x_2,y_2),A)|\hfill \\ =\underset{A\in {\mathcal{B}}_{{\mathbb{R}}_{+}^2}}{\sup }|E\left[{\mathrm{𝟙}}_A(X_t^{(1)},Y_t^{(1)})\right]-E\left[{\mathrm{𝟙}}_A(X_t^{(2)},Y_t^{(2)})\right]|\hfill \\ \le \underset{A\in {\mathcal{B}}_{{\mathbb{R}}_{+}^2}}{\sup }E\left[|{\mathrm{𝟙}}_A(X_t^{(1)},Y_t^{(1)})-{\mathrm{𝟙}}_A(X_t^{(2)},Y_t^{(2)})|\right].\hfill \end{array}$$

Note that in the event $\{t\ge {\tau }_{\text{max}}^X\vee {\tau }_{\text{max}}^Y\}$, $(X_t^{(1)},Y_t^{(1)})=(X_t^{(2)},Y_t^{(2)})$. Thus,

$$\begin{array}{l}|{\mathrm{𝟙}}_A(X_t^{(1)},Y_t^{(1)})-{\mathrm{𝟙}}_A(X_t^{(2)},Y_t^{(2)})|=|{\mathrm{𝟙}}_A(X_t^{(1)},Y_t^{(1)})-{\mathrm{𝟙}}_A(X_t^{(2)},Y_t^{(2)})|·{\mathrm{𝟙}}_{\{t<{\tau }_{\text{max}}^X\vee {\tau }_{\text{max}}^Y\}}\hfill \\ \le 2·{\mathrm{𝟙}}_{\{t<{\tau }_{\text{max}}^X\vee {\tau }_{\text{max}}^Y\}}\hfill \\ \le 2·{\mathrm{𝟙}}_{\{t<{\tau }_{\text{max}}^X\}}+2·{\mathrm{𝟙}}_{\{t<{\tau }_{\text{max}}^Y\}}.\hfill \end{array}$$

Taking expectation on both sides of the last display and then taking supremum over A from ${\mathcal{B}}_{{\mathbb{R}}_{+}^2}$, we have

$$d_{TV}(P_t((x_1,y_1),·),P_t((x_2,y_2),·))\le 2P({\tau }_{\text{max}}^X>t)+2P({\tau }_{\text{max}}^Y>t).$$(3.1)

What remains is to evaluate $P({\tau }_{\text{max}}^X>t)$ and $P({\tau }_{\text{max}}^Y>t)$. Without loss of generality, we assume $x_1\le x_2$; then, ${\tau }_{\text{max}}^X={\tau }^{X^{(2)}}$ and

$$P({\tau }_{\text{max}}^X>t)=P({\tau }^{X^{(2)}}>t).$$(3.2)

On the event $\{t\le {\tau }^{X^{(2)}}\}$, $X^{(2)}$ behaves like an Ornstein–Uhlenbeck process; hence, ${\tau }^{X^{(2)}}$ is the first hitting time to zero for an Ornstein–Uhlenbeck process. By Lemma 13, for any $t_{*}$, there are constants $C_{18}$ and $C_{19}$ depending on ${\theta }_1,{\mu }_1,{\sigma }_1$, and $t_{*}$ such that for $t\ge t_{*}$,

$$\begin{array}{ll}P({\tau }^{X^{(2)}}>t)\hfill & \le \frac{1}{2}\hspace{0.17em}{C}_{18}{e}^{-C_{19}t}+\frac{4}{\sqrt{\pi }}\frac{{\sigma }_1{x}_2}{{\theta }_1^2}\hspace{0.17em}{e}^{\frac{{\theta }_1{x}_2}{{\sigma }_1^2}}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }_1^2}{4{\sigma }_1}t}\hfill \\ \hfill & =\frac{1}{2}\hspace{0.17em}{C}_{18}{e}^{-C_{19}t}+\frac{4}{\sqrt{\pi }}\frac{{\sigma }_1\max (x_1,x_2)}{{\theta }_1^2}\hspace{0.17em}{e}^{\frac{{\theta }_1\max (x_1,x_2)}{{\sigma }_1^2}}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }_1^2}{4{\sigma }_1}t}.\hfill \end{array}$$

Together with (3.2), we have

$$P({\tau }_{\text{max}}^X>t)=\frac{1}{2}\hspace{0.17em}{C}_{18}{e}^{-C_{19}t}+\frac{4}{\sqrt{\pi }}\frac{{\sigma }_1\max (x_1,x_2)}{{\theta }_1^2}\hspace{0.17em}{e}^{\frac{{\theta }_1\max (x_1,x_2)}{{\sigma }_1^2}}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }_1^2}{4{\sigma }_1}t}.$$(3.3)

Similarly, there exist constants $C_{20}$ and $C_{21}$ depending on ${\theta }_2,{\mu }_2,{\sigma }_2$, and $t_{*}$ such that for $t\ge t_{*}$,

$$P({\tau }_{\text{max}}^Y>t)=\frac{1}{2}\hspace{0.17em}{C}_{20}{e}^{-C_{21}t}+\frac{4}{\sqrt{\pi }}\frac{{\sigma }_2\max (y_1,y_2)}{{\theta }_2^2}\hspace{0.17em}{e}^{\frac{{\theta }_2\max (y_1,y_2)}{{\sigma }_2^2}}·\frac{1}{\sqrt{t^3}}{e}^{-\frac{{\theta }_2^2}{4{\sigma }_2}t}.$$(3.4)

The desired result follows by combining (3.1), (3.3), and (3.4). □

With Theorem 14 in hand, we are ready to study the rate of convergence to the stationary distribution in total variation distance.

Theorem 15.

Fix $(x,y)\in {\mathbb{R}}_{+}^2$. For any $t_{*}>0$, there exist constants $C_{22}$, depending on ${\theta }_1,{\theta }_2,{\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2,t_{*},x$, and y, and $C_{23}$, depending on ${\theta }_1,{\theta }_2,{\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2$, and $t_{*}$, such that for $t\ge t_{*}$,

$$d_{TV}(P_t((x,y),·),\pi )\le C_{22}{e}^{-C_{23}t}.$$

Proof.

It can be easily verified that

$$d_{TV}(P_t((x,y),·),\pi )\le {\displaystyle {\int }_{{\mathbb{R}}_{+}^2}{d}_{TV}}(P_t((x,y),·),P_t((\tilde{x},\tilde{y}),·))\hspace{0.17em}\pi (d\tilde{x},d\tilde{y}).$$

Plugging the result of Theorem 14 into the above display, after arrangement, we have, for $t\ge t_{*}$,

$$\begin{array}{l}{d}_{TV}(P_t((x,y),·),\pi )\le C_{18}{e}^{-C_{19}t}+C_{20}{e}^{-C_{21}t}+\hspace{0.17em}\frac{8}{\sqrt{\pi \hspace{0.17em}{t}^3}}{e}^{-\frac{{\theta }_1^2}{4{\sigma }_1}t}{\displaystyle {\int }_{{\mathbb{R}}_{+}^2}\frac{{\sigma }_1\max (x,\tilde{x})}{{\theta }_1^2}}\hspace{0.17em}{e}^{\frac{{\theta }_1\max (x,\tilde{x})}{{\sigma }_1^2}}\hspace{0.17em}\pi (d\tilde{x},d\tilde{y})\hfill \\ +\hspace{0.17em}\frac{8}{\sqrt{\pi \hspace{0.17em}{t}^3}}{e}^{-\frac{{\theta }_2^2}{4{\sigma }_2}t}{\displaystyle {\int }_{{\mathbb{R}}_{+}^2}\frac{{\sigma }_2\max (y,\tilde{y})}{{\theta }_2^2}}\hspace{0.17em}{e}^{\frac{{\theta }_2\max (y,\tilde{y})}{{\sigma }_2^2}}\hspace{0.17em}\pi (d\tilde{x},d\tilde{y}).\hfill \end{array}$$

Note that the marginal measures of $\pi$ are stationary distributions for 1-d ROU processes; hence, they have truncated Gaussian distributions (cf. Ward and Glynn 2003b, proposition 1). Consequently,

$$\begin{array}{l}{\displaystyle {\int }_{{\mathbb{R}}_{+}^2}\frac{{\sigma }_1\max (x,\tilde{x})}{{\theta }_1^2}}\hspace{0.17em}{e}^{\frac{{\theta }_1\max (x,\tilde{x})}{{\sigma }_1^2}}\hspace{0.17em}\pi (d\tilde{x},d\tilde{y})<\infty \hfill \\ {\displaystyle {\int }_{{\mathbb{R}}_{+}^2}\frac{{\sigma }_2\max (y,\tilde{y})}{{\theta }_2^2}}\hspace{0.17em}{e}^{\frac{{\theta }_2\max (y,\tilde{y})}{{\sigma }_2^2}}\hspace{0.17em}\pi (d\tilde{x},d\tilde{y})<\infty .\hfill \end{array}$$

Then, the desired result follows. □

3.3. Exponential Rate of Convergence in Wasserstein Distance

This subsection aims to demonstrate that the rate of convergence in Wasserstein distance is exponential, with the exponent being explicit. Again, we begin by examining the Wasserstein distance between two transition probability measures and then turn to the Wasserstein distance between any transition probability measure and the stationary distribution.

Theorem 16.

For $(x_1,y_1),(x_2,y_2)\in {\mathbb{R}}_{+}^2$, we have, for $t\ge 0$,

$$d_W(P_t((x_1,y_1),·),P_t((x_2,y_2),·))\le |x_1-x_2|e^{-{\theta }_{1}t}+|y_1-y_2|e^{-{\theta }_{2}t}.$$

Proof.

We adopt the notation in the proof of Theorem 14. That is, $(X_t^{(i)},Y_t^{(i)})$ is the 2-d ROU process starting from $(x_i,y_i)$, and ${\tau }^{X^{(i)}}$ and ${\tau }^{Y^{(i)}}$ are the first hitting times to zero.

By the definition of Wasserstein distance,

$$\begin{array}{l}{d}_W(P_t((x_1,y_1),·),P_t((x_2,y_2),·))=\underset{f\in \text{Lip}(1)}{\sup }\left|{\displaystyle \int f}(\tilde{x},\tilde{y})\hspace{0.17em}{P}_t((x_1,y_1),(d\tilde{x},d\tilde{y}))-{\displaystyle \int f}(\tilde{x},\tilde{y})\hspace{0.17em}{P}_t((x_2,y_2),(d\tilde{x},d\tilde{y}))\right|\hfill \\ =\underset{f\in \text{Lip}(1)}{\sup }|E\left[f(X_t^{(1)},Y_t^{(1)})\right]-E\left[f(X_t^{(2)},Y_t^{(2)})\right]|\hfill \\ \le \underset{f\in \text{Lip}(1)}{\sup }E\left[|f(X_t^{(1)},Y_t^{(1)})-f(X_t^{(2)},Y_t^{(2)})|\right].\hfill \end{array}$$

Note that f is a Lipschitz function with Lipschitz constant $\le 1$. Thus,

$$\begin{array}{ll}|f(X_t^{(1)},Y_t^{(1)})-f(X_t^{(2)},Y_t^{(2)})|\hfill & \le \sqrt{{(X_t^{(1)}-X_t^{(2)})}^2+{(Y_t^{(1)}-Y_t^{(2)})}^2}\hfill \\ \hfill & \le |X_t^{(1)}-X_t^{(2)}|+|Y_t^{(1)}-Y_t^{(2)}|.\hfill \end{array}$$

Taking expectation on both sides of the last display and then taking supremum over f from $\text{Lip}(1)$, we have

$$d_W(P_t((x_1,y_1),·),P_t((x_2,y_2),·))\le E\left[|X_t^{(1)}-X_t^{(2)}|\right]+E\left[|Y_t^{(1)}-Y_t^{(2)}|\right].$$(3.5)

In what follows, we will evaluate the expectations of $|X_t^{(1)}-X_t^{(2)}|$ and $|Y_t^{(1)}-Y_t^{(2)}|$. We begin with the expectation of $|X_t^{(1)}-X_t^{(2)}|$. Without loss of generality, we assume that $x_1\le x_2$. Let ${\tilde{X}}_t^{(1)}$ be an Ornstein–Uhlenbeck process satisfying the SDE

$$d{\tilde{X}}_t^{(1)}={\theta }_1({\mu }_1-{\tilde{X}}_t^{(1)})dt+{\sigma }_{1}d{W}_t,$$

starting from $x_1$. Then, ${\tilde{X}}_t^{(1)}$ has the representation

$${\tilde{X}}_t^{(1)}=x_1{e}^{-{\theta }_{1}t}+{\mu }_1\left(1-e^{-{\theta }_{1}t}\right)+{\sigma }_1{\displaystyle {\int }_0^t{e}^{-{\theta }_1(t-s)}}\hspace{0.17em}d{W}_s.$$(3.6)

Using Lemma 2 and Remark 3, we have, for $t\ge 0$,

$${\tilde{X}}_t^{(1)}\le X_t^{(1)}\le X_t^{(2)}.$$

By a similar argument in the proof of Theorem 14, we conclude $X_t^{(1)}=X_t^{(2)}$ on the event $\{t\ge {\tau }^{X^{(2)}}\}$. Thus, combining these two results yields

$$\begin{array}{ll}|X_t^{(1)}-X_t^{(2)}|\hfill & =|X_t^{(1)}-X_t^{(2)}|·{\mathrm{𝟙}}_{\{t<{\tau }^{X^{(2)}}\}}=(X_t^{(2)}-X_t^{(1)})·{\mathrm{𝟙}}_{\{t<{\tau }^{X^{(2)}}\}}\hfill \\ \hfill & \le (X_t^{(2)}-{\tilde{X}}_t^{(1)})·{\mathrm{𝟙}}_{\{t<{\tau }^{X^{(2)}}\}}.\hfill \end{array}$$(3.7)

In the event $\{t<{\tau }^{X^{(2)}}\}$, $X_t^{(2)}$ behaves like an Ornstein–Uhlenbeck process and, hence, has the representation

$$X_t^{(2)}=x_2{e}^{-{\theta }_{1}t}+{\mu }_1\left(1-e^{-{\theta }_{1}t}\right)+{\sigma }_1{\displaystyle {\int }_0^t{e}^{-{\theta }_1(t-s)}}\hspace{0.17em}d{W}_s.$$(3.8)

Then, plugging (3.6) and (3.8) into (3.7) gives

$$|X_t^{(1)}-X_t^{(2)}|\le (x_2-x_1)e^{-{\theta }_{1}t}·{\mathrm{𝟙}}_{\{t<{\tau }^{X^{(2)}}\}}\le (x_2-x_1)e^{-{\theta }_{1}t}=|x_2-x_1|e^{-{\theta }_{1}t}.$$

Thus,

$$E\left[|X_t^{(1)}-X_t^{(2)}|\right]\le |x_2-x_1|e^{-{\theta }_{1}t}.$$(3.9)

Similarly,

$$E\left[|Y_t^{(1)}-Y_t^{(2)}|\right]\le |y_2-y_1|e^{-{\theta }_{2}t}.$$(3.10)

Combining (3.5), (3.9), and (3.10), the desired result follows. □

Theorem 17.

For $(x,y)\in {\mathbb{R}}_{+}^2$, there exist constants $C_{24}$, depending on ${\theta }_1,{\mu }_1,{\sigma }_1$, and x, and $C_{25}$, depending on ${\theta }_2,{\mu }_2,{\sigma }_2$, and y, such that for $t\ge 0$,

$$d_W(P_t((x,y),·),\pi )\le C_{24}{e}^{-{\theta }_{1}t}+C_{25}{e}^{-{\theta }_{2}t}.$$

Proof.

It can be easily verified that

$$d_W(P_t((x,y),·),\pi )\le {\displaystyle {\int }_{{\mathbb{R}}_{+}^2}{d}_W}(P_t((x,y),·),P_t((\tilde{x},\tilde{y}),·))\hspace{0.17em}\pi (d\tilde{x},d\tilde{y}).$$

Together with Theorem 16, we have

$$\begin{array}{ll}{d}_W(P_t((x,y),·),\pi )\hfill & \le {\displaystyle {\int }_{{\mathbb{R}}_{+}^2}(|x-\tilde{x}|e^{-{\theta }_{1}t}+|y-\tilde{y}|e^{-{\theta }_{2}t})}\hspace{0.17em}\pi (d\tilde{x},d\tilde{y})\hfill \\ \hfill & =e^{-{\theta }_{1}t}{\displaystyle {\int }_{{\mathbb{R}}_{+}^2}|}x-\tilde{x}|\hspace{0.17em}\pi (d\tilde{x},d\tilde{y})+e^{-{\theta }_{2}t}{\displaystyle {\int }_{{\mathbb{R}}_{+}^2}|}y-\tilde{y}|\hspace{0.17em}\pi (d\tilde{x},d\tilde{y}).\hfill \end{array}$$

Because the marginal measure ${\int }_{{\mathbb{R}}_{+}}\pi (·,dy)$ has a truncated Gaussian distribution,

$${\displaystyle {\int }_{{\mathbb{R}}_{+}^2}|}x-\tilde{x}|\pi (d\tilde{x},d\tilde{y})<\infty ,$$

and its value depends on ${\theta }_1,{\mu }_1,{\sigma }_1$, and x. Similarly,

$${\displaystyle {\int }_{{\mathbb{R}}_{+}^2}|}y-\tilde{y}|\pi (d\tilde{x},d\tilde{y})<\infty ,$$

and its value depends on ${\theta }_2,{\mu }_2,{\sigma }_2$, and y. We conclude the proof. □

4. A Numerical Scheme for the Stationary Distribution

The purpose of this section is to numerically study the stationary distribution for the 2-d ROU process. For diffusions, two common approaches for investigating the stationary distribution are the PDE method and the Monte Carlo method. We begin by presenting the PDE that the density of the stationary distribution satisfies, along with the corresponding boundary conditions, provided the existence of the density. However, because of its complexity, we are unable to solve this PDE ourselves and defer it to the experts in PDE. Therefore, we turn to the Monte Carlo method, which is partly inspired by Budhiraja et al. (2014). Additionally, the Monte Carlo method has an advantage in higher dimensions. Indeed, the approach we will elaborate on later can be automatically extended to high-dimensional ROU processes.

4.1. Stationary Distribution: A PDE Perspective

Assume the stationary distribution, denoted $\pi$, is absolutely continuous with density function $p(x,y)$. Indeed, we believe the assumption can be justified by modifying the argument in Harrison and Williams (1987a, section 7). However, as this is beyond the scope of the paper, we omit the detailed verification. In the following, we will characterize $p(x,y)$ by introducing the corresponding PDE.

For any smooth function f on $[0,\infty )\times [0,\infty )$ with compact support on $(0,\infty )\times (0,\infty )$, an application of Itô’s lemma gives

where $\mathcal{A}f$ is defined as

Noting that $\frac{\partial f}{\partial x}$, $\frac{\partial f}{\partial y}$ are bounded and $\frac{\partial f}{\partial x}(0,y)=\frac{\partial f}{\partial y}(x,0)=0$, it follows by taking expectation that

Because $\pi$ is the stationary distribution for $(X_t,Y_t)$, we have $E_{\pi }[f(X_t,Y_t)]=E_{\pi }[f(X_0,Y_0)]$. Then,

Hence,

Because we assume $p(x,y)$ is the density function of $\pi$, we have

Using integration by parts yields

Here, $\frac{\partial p}{\partial x},\frac{\partial p}{\partial y},\frac{{\partial }^{2}p}{\partial x^2},\dots$ can be understood as weak derivatives. By the arbitrariness of f, we conclude

We now turn to the boundary conditions for $p(x,y)$. Because $X_t$ and $Y_t$ are one-dimensional reflected Ornstein–Uhlenbeck processes, the marginal distributions of $\pi$ correspond to the stationary distributions of $X_t$ and $Y_t$; that is,

where $\varphi (·)$ and $\mathrm{\Phi }(·)$ are the density and distribution functions of a standard normal random variable. The explicit representations of $p_1(x)$ and $p_2(y)$ come from Ward and Glynn (2003b).

In summary, $p(x,y)$ satisfies the following PDE:

subject to the boundary conditions

We are unable to solve the above PDE, as its boundary conditions do not fall into standard categories such as Dirichlet, Neumann, or Robin. So we entrust the task to experts in PDE.

4.2. Numerical Scheme and Its Convergence

In this subsection, we will develop a convergent numerical procedure to approximate the stationary distribution of our ROU process. Consistently, the stationary distribution is denoted by $\pi$. We begin by constructing a discrete-time process that serves as an approximation to the ROU process.

Let ${\{{\lambda }_k\}}_{k=1}^{\infty }$ be a sequence of real numbers satisfying the following conditions:

1. $0<{\lambda }_k<1/\max \{{\theta }_1,{\theta }_2\}$ for all $k\in {\mathbb{N}}_{+}$;

2. ${\lambda }_k\to 0$ as $k\to \infty$;

3. ${\mathrm{\Lambda }}_n≔{\sum }_{k=1}^n{\lambda }_k\to \infty$ as $n\to \infty$.

Here, ${\mathbb{N}}_{+}$ is the set of all positive integers. Note that the above conditions are satisfied if ${\lambda }_k=1/(\max \{{\theta }_1,{\theta }_2\}{(k+1)}^{\nu })$ with $\nu \in (0,1]$. Let ${\{(G_{k,1},G_{k,2})\}}_{k=1}^{\infty }$ be a sequence of independent and identically distributed (i.i.d.) random vectors in ${\mathbb{R}}^2$ with mean $\mathbf{0}$ and covariance matrix $(\begin{array}{cc}1& \rho \\ \rho & 1\end{array})$. Furthermore, we pose the following assumption on $G_{k,i}$’s:

for some $\beta \in (0,\infty )$ and for all $k\in {\mathbb{N}}_{+}$, $i=1,2$, and $\lambda \in \mathbb{R}$. The above assumption is satisfied when $(G_{k,1},G_{k,2})$ has a joint Gaussian distribution. As usual, we define ${\mathcal{G}}_k=\sigma (G_{1,1},G_{1,2},\dots ,G_{k,1},G_{k,2})$ for $k\in {\mathbb{N}}_{+}$ and ${\mathcal{G}}_0=\{\varnothing ,\mathrm{\Omega }\}$.

With the above preparation, we are ready to iteratively construct a discrete-time process in ${\mathbb{R}}_{+}^2$. For $(x_0,y_0)\in {\mathbb{R}}_{+}^2$, define

where ${[z]}^{+}≔\max \{z,0\}$. Obviously, ${\{(U_k,V_k)\}}_{k=0}^{\infty }$ is a sequence of random vectors with values in ${\mathbb{R}}_{+}^2$. Employing these random vectors, we construct a sequence of random measures on ${\mathbb{R}}_{+}^2$ that approximate the stationary distribution $\pi$ as follows:where ${\delta }_{(x,y)}$ is the Dirac measure on ${\mathbb{R}}^2$. In particular, these random measures yield an approximation for any integral of the form ${\int }_{{\mathbb{R}}_{+}^2}f(x,y)\hspace{0.17em}d\pi (x,y)$ through the corresponding weighted averages:

With these random measures in hand, we are now prepared to introduce the main result of this section.

Theorem 18.

As $n\to \infty$, ${\pi }_n$ converges weakly to $\pi$ almost surely.

The proof of Theorem 18 is relegated to the following two subsections.

4.3. Tightness

The aim of this subsection is to establish the tightness of ${\{{\pi }_n\}}_{n=1}^{\infty }$ almost surely. To this end, we prepare the following lemmas.

Lemma 19.

For nonnegative integers n, l, we have

$$\begin{array}{l}{U}_{n+l}\le U_n{\displaystyle \prod _{j=n+1}^{n+l}(1-{\theta }_1{\lambda }_j)}+{\mu }_1^{+}\hspace{0.17em}+{\sigma }_1\underset{0\le i\le l}{\max }{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1},\end{array}$$(4.2)

and

$$\begin{array}{l}{V}_{n+l}\le V_n{\displaystyle \prod _{j=n+1}^{n+l}(1-{\theta }_2{\lambda }_j)}+{\mu }_2^{+}\hspace{0.17em}+{\sigma }_2\underset{0\le i\le l}{\max }{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_2{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,2}.\end{array}$$(4.3)

Proof.

We only prove (4.2), and (4.3) follows similarly.

Define

$$L_{k+1}=-\min \{U_k+{\theta }_1({\mu }_1-U_k){\lambda }_{k+1}+{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1},0\}.$$

Then,

$$\begin{array}{ll}{U}_{k+1}\hfill & ={[U_k+{\theta }_1({\mu }_1-U_k){\lambda }_{k+1}+{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1}]}^{+}\hfill \\ \hfill & =U_k+{\theta }_1({\mu }_1-U_k){\lambda }_{k+1}+{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1}+L_{k+1}\hfill \\ \hfill & =(1-{\theta }_1{\lambda }_{k+1})U_k+{\theta }_1{\mu }_1{\lambda }_{k+1}+{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1}+L_{k+1}.\hfill \end{array}$$

Dividing by ${\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)$ on both sides yields

$$\begin{array}{l}\frac{U_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}=\frac{U_k}{{\prod }_{j=1}^k(1-{\theta }_1{\lambda }_j)}+\frac{{\theta }_1{\mu }_1{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}+\frac{{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}+\frac{L_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}.\hfill \end{array}$$

In the last display, by letting $k=n,n+1,\dots ,n+l-1$, and summing up these l equations, we have

$$\begin{array}{l}\frac{U_{n+l}}{{\prod }_{j=1}^{n+l}(1-{\theta }_1{\lambda }_j)}=\frac{U_n}{{\prod }_{j=1}^n(1-{\theta }_1{\lambda }_j)}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\theta }_1{\mu }_1{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+\sum _{k=n}^{n+l-1}\frac{L_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}.\hfill \end{array}$$

For simplicity of notation, we define

$$F_{n:l}≔\frac{U_n}{{\prod }_{j=1}^n(1-{\theta }_1{\lambda }_j)}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\theta }_1{\mu }_1{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}},$$

and

$${\tilde{L}}_{n:l}≔{\displaystyle \sum _{k=n}^{n+l-1}\frac{L_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}.$$

Thus,

$$\frac{U_{n+l}}{{\prod }_{j=1}^{n+l}(1-{\theta }_1{\lambda }_j)}=F_{n:l}+{\tilde{L}}_{n:l}.$$

At this stage, we fix n and allow l to vary. Then, $U_{n+l}/{\prod }_{j=1}^{n+l}(1-{\theta }_1{\lambda }_j)$ is nonnegative for all $l\in {\mathbb{N}}_{+}$. ${\tilde{L}}_{n:l}$, viewed as a function of l, is nondecreasing. Furthermore, ${\tilde{L}}_{n:l}$ increases at l only if $L_{n+l}>0$, which occurs only if $U_{n+l}=0$ according to the definition of $L_{n+l}$. Combining the above observations, together with the uniqueness of the Skorokhod map, it follows that

$${\tilde{L}}_{n:l}=\max \{0,-F_{n:1},-F_{n:2},\dots ,-F_{n:l}\}.$$

Thus,

$$\begin{array}{ll}\frac{U_{n+l}}{{\prod }_{j=1}^{n+l}(1-{\theta }_1{\lambda }_j)}\hfill & =F_{n:l}+\max \{0,-F_{n:1},-F_{n:2},\dots ,-F_{n:l}\}\hfill \\ \hfill & =\max \{F_{n:l},F_{n:l}-F_{n:1},F_{n:l}-F_{n:2},\dots ,F_{n:l}-F_{n:(l-1)},0\}.\hfill \end{array}$$

Note that for $1\le i\le l-1$,

$$\begin{array}{ll}{F}_{n:l}-F_{n:i}\hfill & ={\displaystyle \sum _{k=n+i}^{n+l-1}\frac{{\theta }_1{\mu }_1{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+{\displaystyle \sum _{k=n+i}^{n+l-1}\frac{{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}\hfill \\ \hfill & \le {\displaystyle \sum _{k=n+i}^{n+l-1}\frac{{\theta }_1{\mu }_1^{+}{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+{\displaystyle \sum _{k=n+i}^{n+l-1}\frac{{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}\hfill \\ \hfill & ={\displaystyle \sum _{k=n+i}^{n+l-1}\frac{{\theta }_1{\mu }_1^{+}{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+{\displaystyle \sum _{m=1}^{l-i}\frac{{\sigma }_1\sqrt{{\lambda }_{n+l+1-m}}{G}_{n+l+1-m,1}}{{\prod }_{j=1}^{n+l+1-m}(1-{\theta }_1{\lambda }_j)}},\hfill \end{array}$$

where, in the last equality, we substitute m for $n+l-k$. Furthermore,

$$\begin{array}{ll}{F}_{n:l}\hfill & \le \frac{U_n}{{\prod }_{j=1}^n(1-{\theta }_1{\lambda }_j)}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\theta }_1{\mu }_1^{+}{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\sigma }_1\sqrt{{\lambda }_{k+1}}{G}_{k+1,1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}\hfill \\ \hfill & =\frac{U_n}{{\prod }_{j=1}^n(1-{\theta }_1{\lambda }_j)}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\theta }_1{\mu }_1^{+}{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+{\displaystyle \sum _{m=1}^l\frac{{\sigma }_1\sqrt{{\lambda }_{n+l+1-m}}{G}_{n+l+1-m,1}}{{\prod }_{j=1}^{n+l+1-m}(1-{\theta }_1{\lambda }_j)}}.\hfill \end{array}$$

Combining the last three displays, together with the nonnegativity of ${\theta }_1{\mu }_1^{+}{\lambda }_{k+1}$, we have

$$\begin{array}{l}\frac{U_{n+l}}{{\prod }_{j=1}^{n+l}(1-{\theta }_1{\lambda }_j)}\le \frac{U_n}{{\prod }_{j=1}^n(1-{\theta }_1{\lambda }_j)}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\theta }_1{\mu }_1^{+}{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+\hspace{0.17em}0\vee \underset{0\le i\le l-1}{\max }{\displaystyle \sum _{m=1}^{l-i}\frac{{\sigma }_1\sqrt{{\lambda }_{n+l+1-m}}{G}_{n+l+1-m,1}}{{\prod }_{j=1}^{n+l+1-m}(1-{\theta }_1{\lambda }_j)}}\hfill \\ =\frac{U_n}{{\prod }_{j=1}^n(1-{\theta }_1{\lambda }_j)}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\theta }_1{\mu }_1^{+}{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+\hspace{0.17em}0\vee \underset{1\le i\le l}{\max }{\displaystyle \sum _{m=1}^i\frac{{\sigma }_1\sqrt{{\lambda }_{n+l+1-m}}{G}_{n+l+1-m,1}}{{\prod }_{j=1}^{n+l+1-m}(1-{\theta }_1{\lambda }_j)}}\hfill \\ =\frac{U_n}{{\prod }_{j=1}^n(1-{\theta }_1{\lambda }_j)}+{\displaystyle \sum _{k=n}^{n+l-1}\frac{{\theta }_1{\mu }_1^{+}{\lambda }_{k+1}}{{\prod }_{j=1}^{k+1}(1-{\theta }_1{\lambda }_j)}}+\hspace{0.17em}\underset{0\le i\le l}{\max }{\displaystyle \sum _{m=1}^i\frac{{\sigma }_1\sqrt{{\lambda }_{n+l+1-m}}{G}_{n+l+1-m,1}}{{\prod }_{j=1}^{n+l+1-m}(1-{\theta }_1{\lambda }_j)}},\hfill \end{array}$$

where, in the second last equality, we substitute i for $l-i$, and the last equality follows by the convention that ${\sum }_{m=1}^0{x}_m=0$. In the last display, multiplying by ${\prod }_{j=1}^{n+l}(1-{\theta }_1{\lambda }_j)$ on both sides, after changing the dummy variables, we have

$$\begin{array}{l}{U}_{n+l}\le U_n{\displaystyle \prod _{j=n+1}^{n+l}(1-{\theta }_1{\lambda }_j)}+{\displaystyle \sum _{k=n+1}^{n+l}{\theta }_1}{\mu }_1^{+}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}{\lambda }_k+\hspace{0.17em}\underset{0\le i\le l}{\max }{\displaystyle \sum _{m=1}^i{\sigma }_1}{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}\sqrt{{\lambda }_{n+l+1-m}}{G}_{n+l+1-m,1}.\hfill \end{array}$$

To prove (4.2), it remains to prove

$${\displaystyle \sum _{k=n+1}^{n+l}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_k\le \frac{1}{{\theta }_1}.$$(4.4)

Indeed, noting that ${\lambda }_n\le 1/{\theta }_1$ for $\forall n\in {\mathbb{N}}_{+}$, it follows that

$$\begin{array}{l}{\displaystyle \sum _{k=n+1}^{n+l}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_k={\displaystyle \sum _{k=n+2}^{n+l}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_k+{\displaystyle \prod _{j=n+2}^{n+l}(1-{\theta }_1{\lambda }_j)}{\lambda }_{n+1}\hfill \\ \le {\displaystyle \sum _{k=n+2}^{n+l}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_k+{\displaystyle \prod _{j=n+2}^{n+l}(1-{\theta }_1{\lambda }_j)}\hspace{0.17em}\frac{1}{{\theta }_1}\hfill \end{array}$$(4.5)

$$\begin{array}{l}={\displaystyle \sum _{k=n+3}^{n+l}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_k+{\displaystyle \prod _{j=n+3}^{n+l}(1-{\theta }_1{\lambda }_j)}{\lambda }_{n+2}+{\displaystyle \prod _{j=n+3}^{n+l}(1-{\theta }_1{\lambda }_j)}\left(\frac{1}{{\theta }_1}-{\lambda }_{n+2}\right)\hfill \\ ={\displaystyle \sum _{k=n+3}^{n+l}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_k+{\displaystyle \prod _{j=n+3}^{n+l}(1-{\theta }_1{\lambda }_j)}\hspace{0.17em}\frac{1}{{\theta }_1}.\hfill \end{array}$$(4.6)

Comparing (4.3) and (4.6) and repeating this procedure, we have

$$\begin{array}{l}{\displaystyle \sum _{k=n+2}^{n+l}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_k+{\displaystyle \prod _{j=n+2}^{n+l}(1-{\theta }_1{\lambda }_j)}\hspace{0.17em}\frac{1}{{\theta }_1}={\displaystyle \sum _{k=n+l}^{n+l}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_k+{\displaystyle \prod _{j=n+l}^{n+l}(1-{\theta }_1{\lambda }_j)}\hspace{0.17em}\frac{1}{{\theta }_1}\hfill \\ ={\lambda }_{n+l}+(1-{\theta }_1{\lambda }_{n+l})\hspace{0.17em}\frac{1}{{\theta }_1}\hfill \\ =\frac{1}{{\theta }_1}.\hfill \end{array}$$

Combining the last two displays yields (4.4). We conclude the proof. □

Lemma 20.

For nonnegative integers n, l, and any $r\in {\mathbb{R}}_{+}$, we have

$$E\left[\exp \left(r{\sigma }_1\underset{0\le i\le l}{\max }{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1}\right)\right]\le 4e^{\beta r^2{\sigma }_1^2/{\theta }_1}$$(4.7)

and

$$E\left[\exp \left(r{\sigma }_2\underset{0\le i\le l}{\max }{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_2{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,2}\right)\right]\le 4e^{\beta r^2{\sigma }_2^2/{\theta }_2}.$$

Consequently, it follows by Hölder’s inequality that

$$\begin{array}{l}E[\exp \left(r{\sigma }_1\underset{0\le i\le l}{\max }{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1}\right)\hfill \\ \hspace{1em}\times \exp \left(r{\sigma }_2\underset{0\le i\le l}{\max }{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_2{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,2}\right)]\hfill \\ \le 4e^{2\beta r^2{\sigma }_1^2/{\theta }_1+2\beta r^2{\sigma }_2^2/{\theta }_2}.\hfill \end{array}$$

Proof.

By symmetry, we only prove (4.7).

Note that

$$\begin{array}{l}E\left[\exp \left(r{\sigma }_1\underset{0\le i\le l}{\max }{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1}\right)\right]\hfill \\ =E\left[\underset{0\le i\le l}{\max }\exp \left(r{\sigma }_1{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1}\right)\right].\hfill \end{array}$$

For fixed n and l, because $G_{n+1,1},G_{n+2,1},\dots ,G_{n+l,1}$ are independent with a zero mean,

$${\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1}$$

forms a martingale with respect to $\{{\mathcal{G}}_i\}$. Thus,

$$\exp \left(r{\sigma }_1{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1}\right)$$

is a submartingale. From Doob’s maximal inequality for submartingales, we have

$$\begin{array}{l}E\left[\underset{0\le i\le l}{\max }\exp \left(r{\sigma }_1{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1}\right)\right]\hfill \\ \le 4E\left[\exp \left(r{\sigma }_1{\displaystyle \sum _{m=1}^l{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1}\right)\right]\hfill \\ =4{\displaystyle \prod _{m=1}^{l}E}\left[\exp \left(r{\sigma }_1{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}\sqrt{{\lambda }_{n+l+1-m}}\hspace{0.17em}{G}_{n+l+1-m,1}\right)\right]\hfill \\ \le 4{\displaystyle \prod _{m=1}^l\exp }\left(\beta r^2{\sigma }_1^2{\displaystyle \prod _{j=n+l+2-m}^{n+l}{(1-{\theta }_1{\lambda }_j)}^2}{\lambda }_{n+l+1-m}\right)\hfill \\ \le 4{\displaystyle \prod _{m=1}^l\exp }\left(\beta r^2{\sigma }_1^2{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}{\lambda }_{n+l+1-m}\right)\hfill \\ =4\exp \left(\beta r^2{\sigma }_1^2{\displaystyle \sum _{m=1}^l{\displaystyle \prod _{j=n+l+2-m}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_{n+l+1-m}\right)\hfill \\ =4\exp \left(\beta r^2{\sigma }_1^2{\displaystyle \sum _{k=n+1}^{n+l}{\displaystyle \prod _{j=k+1}^{n+l}(1-{\theta }_1{\lambda }_j)}}{\lambda }_k\right)\hfill \\ \le 4e^{\beta r^2{\sigma }_1^2/{\theta }_1},\hfill \end{array}$$

where the second inequality follows by the Assumption (4.1) on $G_{k,i}$; in the last equality, we have substituted k for $n+l+1-m$; and in the last inequality, we have invoked (4.4). We complete the proof. □

With these two lemmas in hand, we can establish a uniform upper bound for the expectations of $\exp (r{U}_l+r{V}_l)$. This result is summarized in the following lemma.

Lemma 21.

For any $r\in {\mathbb{R}}_{+}$, we have

$$\underset{0\le l<\infty }{\sup }E\left[e^{r{U}_l+r{V}_l}\right]\le 4e^{r(x_0+y_0+{\mu }_1^{+}+{\mu }_2^{+})+2\beta r^2{\sigma }_1^2/{\theta }_1+2\beta r^2{\sigma }_2^2/{\theta }_2}.$$

Proof.

Letting $n=0$ in Lemma 19 and noting that $U_0=x_0$ yield

$$\begin{array}{l}{U}_l\le U_0\prod _{j=1}^l(1-{\theta }_1{\lambda }_j)+{\mu }_1^{+}+{\sigma }_1\underset{0\le i\le l}{\max }\sum _{m=1}^i\prod _{j=l+2-m}^l(1-{\theta }_1{\lambda }_j)\sqrt{{\lambda }_{l+1-m}}\hspace{0.17em}{G}_{l+1-m,1}\\ \le x_0+{\mu }_1^{+}+{\sigma }_1\underset{0\le i\le l}{\max }\sum _{m=1}^i\prod _{j=l+2-m}^l(1-{\theta }_1{\lambda }_j)\sqrt{{\lambda }_{l+1-m}}\hspace{0.17em}{G}_{l+1-m,1}.\end{array}$$

Similarly,

$$V_l\le y_0+{\mu }_2^{+}+{\sigma }_2\underset{0\le i\le l}{\max }\sum _{m=1}^i\prod _{j=l+2-m}^l(1-{\theta }_2{\lambda }_j)\sqrt{{\lambda }_{l+1-m}}\hspace{0.17em}{G}_{l+1-m,2}.$$

Combining the last two displays with Lemma 20, the desired result follows. □

Before presenting the next result, we introduce some necessary notation. Define $\lambda :[0,\infty )\to [0,\infty )$ and $\iota :[0,\infty )\to N$ as

where we define ${\mathrm{\Lambda }}_0=0$. Furthermore, define

For any $t>0$, it can be easily verified that

Lemma 22.

For any $r\in {\mathbb{R}}_{+}$, define

$$\begin{array}{ll}M(r)\hfill & ≔r{\mu }_1^{+}+r{\mu }_2^{+}+2\beta r^2{\sigma }_1^2/{\theta }_1+2\beta r^2{\sigma }_2^2/{\theta }_2+\ln 4,\hfill \\ \mathrm{\Delta }\hfill & ≔{\lambda }_0+\max \left\{\ln 2/{\theta }_1,\ln 2/{\theta }_2\right\}.\hfill \end{array}$$

Then, for any $t\ge 0$, we have

$$E\left[e^{r{U}_{\iota (t+\mathrm{\Delta })}+r{V}_{\iota (t+\mathrm{\Delta })}}|{\mathcal{G}}_{\iota (t)}\right]\le e^{-1}{e}^{r{U}_{\iota (t)}+r{V}_{\iota (t)}}+e^{2M(r)+1}.$$

Proof.

Applying Lemma 19 and Lemma 20, together with the independence of the sequence ${\{(G_{n,1},G_{n,2})\}}_{n=1}^{\infty }$, we have

$$\begin{array}{l}E\left[e^{r{U}_{\iota (t+\mathrm{\Delta })}+r{V}_{\iota (t+\mathrm{\Delta })}}|{\mathcal{G}}_{\iota (t)}\right]\le \exp (r{U}_{\iota (t)}{\displaystyle \prod _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}(1-{\theta }_1{\lambda }_j)}+r{V}_{\iota (t)}{\displaystyle \prod _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}(}1-{\theta }_2{\lambda }_j))\times e^{r{\mu }_1^{+}+r{\mu }_2^{+}}\\ \times E\left[\exp \right(r{\sigma }_1\underset{0\le i\le \iota (t+\mathrm{\Delta })-\iota (t)}{\max }{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=\iota (t+\mathrm{\Delta })+2-m}^{\iota (t+\mathrm{\Delta })}(1-{\theta }_1{\lambda }_j)}}\sqrt{{\lambda }_{\iota (t+\mathrm{\Delta })+1-m}}\hspace{0.17em}{G}_{\iota (t+\mathrm{\Delta })+1-m,1})\hfill \\ \hspace{1em}\hspace{1em}\hspace{0.17em}\times \exp (r{\sigma }_2\underset{0\le i\le \iota (t+\mathrm{\Delta })-\iota (t)}{\max }{\displaystyle \sum _{m=1}^i{\displaystyle \prod _{j=\iota (t+\mathrm{\Delta })+2-m}^{\iota (t+\mathrm{\Delta })}(1-{\theta }_2{\lambda }_j)}}\sqrt{{\lambda }_{\iota (t+\mathrm{\Delta })+1-m}}\hspace{0.17em}{G}_{\iota (t+\mathrm{\Delta })+1-m,2})]\hfill \end{array}$$

$$\begin{array}{l}\le \exp (r{U}_{\iota (t)}{\displaystyle \prod _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}(1-{\theta }_1{\lambda }_j)}+r{V}_{\iota (t)}{\displaystyle \prod _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}(}1-{\theta }_2{\lambda }_j))\times 4e^{r{\mu }_1^{+}+r{\mu }_2^{+}+2\beta r^2{\sigma }_1^2/{\theta }_1+2\beta r^2{\sigma }_2^2/{\theta }_2}\\ =e^{M(r)}·\exp (r{U}_{\iota (t)}{\displaystyle \prod _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}(1-{\theta }_1{\lambda }_j)}+r{V}_{\iota (t)}{\displaystyle \prod _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}(}1-{\theta }_2{\lambda }_j)),\hfill \end{array}$$(4.8)

where the last equality follows by the definition of $M(r)$. We then claim

$${\displaystyle \sum _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}{\lambda }_j}=\lambda (t+\mathrm{\Delta })-\lambda (t).$$

In fact, if $\iota (t)=k$ and $\iota (t+\mathrm{\Delta })=l$, then $\lambda (t)={\mathrm{\Lambda }}_k$ and $\lambda (t+\mathrm{\Delta })={\mathrm{\Lambda }}_l$. Therefore,

$$\lambda (t+\mathrm{\Delta })-\lambda (t)={\mathrm{\Lambda }}_l-{\mathrm{\Lambda }}_k={\displaystyle \sum _{j=k+1}^l{\lambda }_j}={\displaystyle \sum _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}{\lambda }_j}.$$

Note that $\lambda (t+\mathrm{\Delta })-\lambda (t)\ge \mathrm{\Delta }-{\lambda }_0$. It follows that

$$\begin{array}{l}{\displaystyle \prod _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}(1-{\theta }_1{\lambda }_j)}\le {\displaystyle \prod _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}{e}^{-{\theta }_1{\lambda }_j}}=\exp (-{\theta }_1{\displaystyle \sum _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}{\lambda }_j})\hfill \\ =e^{-{\theta }_1(\lambda (t+\mathrm{\Delta })-\lambda (t))}\le e^{-{\theta }_1(\mathrm{\Delta }-{\lambda }_0)}\le e^{-{\theta }_1(\ln 2/{\theta }_1+{\lambda }_0-{\lambda }_0)}=\frac{1}{2}.\hfill \end{array}$$

Similarly,

$${\displaystyle \prod _{j=\iota (t)+1}^{\iota (t+\mathrm{\Delta })}(1-{\theta }_2{\lambda }_j)}\le \frac{1}{2}.$$

Combining the last two displays with (4.8), we have

$$\begin{array}{l}E\left[e^{r{U}_{\iota (t+\mathrm{\Delta })}+r{V}_{\iota (t+\mathrm{\Delta })}}|{\mathcal{G}}_{\iota (t)}\right]\le e^{M(r)}·e^{r{U}_{\iota (t)}/2+r{V}_{\iota (t)}/2}\hfill \\ =e^{M(r)}{e}^{r{U}_{\iota (t)}/2+r{V}_{\iota (t)}/2}{\mathrm{𝟙}}_{\{U_{\iota (t)}+V_{\iota (t)}>2(M(r)+1)/r\}}+\hspace{0.17em}{e}^{M(r)}{e}^{r{U}_{\iota (t)}/2+r{V}_{\iota (t)}/2}{\mathrm{𝟙}}_{\{U_{\iota (t)}+V_{\iota (t)}\le 2(M(r)+1)/r\}}\hfill \\ \le e^{M(r)}{e}^{r{U}_{\iota (t)}+r{V}_{\iota (t)}-\frac{r}{2}\times 2(M(r)+1)/r}{\mathrm{𝟙}}_{\{U_{\iota (t)}+V_{\iota (t)}>2(M(r)+1)/r\}}+\hspace{0.17em}{e}^{M(r)}{e}^{\frac{r}{2}\times 2(M(r)+1)/r}{\mathrm{𝟙}}_{\{U_{\iota (t)}+V_{\iota (t)}\le 2(M(r)+1)/r\}}\hfill \\ \le e^{-1}{e}^{r{U}_{\iota (t)}+r{V}_{\iota (t)}}+e^{2M(r)+1}.\hfill \end{array}$$

We conclude the proof. □

With the above preparation, we arrive at our primary result in this subsection, indicating the almost sure tightness of the sequence of random measures ${\{{\pi }_n\}}_{n=1}^{\infty }$.

Proposition 23.

For any $r\in {\mathbb{R}}_{+}$, we have

$$\underset{1\le n<\infty }{\sup }\frac{1}{{\mathrm{\Lambda }}_n}{\displaystyle \sum _{k=1}^n{\lambda }_k}{e}^{r{U}_{k-1}+r{V}_{k-1}}<\infty$$

almost surely. In other words,

$$\underset{1\le n<\infty }{\sup }{\displaystyle {\int }_{{\mathbb{R}}_{+}^2}{e}^{rx+ry}}\hspace{0.17em}{\pi }_n(dx,dy)<\infty$$

almost surely. Consequently, the sequence ${\{{\pi }_n\}}_{n=1}^{\infty }$ is tight almost surely.

Proof.

The proof is similar to that of Budhiraja et al. (2014, lemma 6); therefore, we omit the specific details. □

4.4. Identification of the Limit

To complete the proof of Theorem 18, we must demonstrate that for almost every $\omega$, any weak limit of ${\pi }_n(\omega )$ is $\pi$. To this end, we follow the method in Budhiraja et al. (2014, section 2.2). It is worth noting that all the conditions necessary for the deduction in Budhiraja et al. (2014, section 2.2) are satisfied, except for the requirement of boundedness in $|b(x)|$, which corresponds to the boundedness of $|{\mu }_1-x|$ and $|{\mu }_2-y|$ in our case. However, $|{\mu }_1-x|\le |{\mu }_1|+|x|$ and $|{\mu }_2-y|\le |{\mu }_2|+|y|$. When modifying the proof in Budhiraja et al. (2014, section 2.2), it only requires the following results:

• (R1) $\frac{1}{{\mathrm{\Lambda }}_n}{\sum }_{k=1}^n{\lambda }_k^2{(U_{k-1}+V_{k-1})}^2\to 0$ almost surely as $n\to \infty$;

• (R2) For any nonnegative integer j,

  $$\underset{1\le n<\infty }{\sup }\frac{1}{{\mathrm{\Lambda }}_n}{\displaystyle \sum _{k=1}^n{\lambda }_k}{(U_{k-1}+V_{k-1})}^j{\mathrm{𝟙}}_{\{U_{k-1}+V_{k-1}>M\}}\to 0$$
  almost surely as $M\to \infty$;

• (R3)

  $${\displaystyle \sum _{k=1}^{\infty }\frac{{\lambda }_{k+1}^{2}E\left[{(U_k+V_k)}^2\right]}{{\mathrm{\Lambda }}_k^2}}<\infty .$$

All the results can be derived from Lemma 21 and Proposition 23. Indeed, note that

Therefore, (R1) follows immediately from Proposition 23 and the fact that ${\lambda }_n\to 0$ and ${\mathrm{\Lambda }}_n\to \infty$. Furthermore, (R2) follows similarly because

Finally, (R3) follows from the uniform boundedness of $E[{(U_k+V_k)}^2]$, which is a direct result of Lemma 21. This completes the convergence analysis of the numerical scheme.

4.5. Numerical Results

In this subsection, we apply our numerical scheme to approximate the stationary distribution of the 2-d ROU process under the parameter setting ${\theta }_1={\theta }_2={\mu }_1={\mu }_2={\sigma }_1={\sigma }_2=1$, $\rho =0.5$, with initial conditions $x_0=1$ and $y_0=1$. The resulting numerical joint and marginal distributions are computed using three different schemes: (i) ${\lambda }_k={(1+k)}^{-0.1}$, (ii) ${\lambda }_k={(1+k)}^{-0.7}$, and (iii) ${\lambda }_k={(1+k)}^{-1}$. The number of iterations is set to $n=200,000$. Additionally, the true marginal distributions, with density functions given by $\sqrt{2}\varphi (\sqrt{2}(x-1))/(1-\mathrm{\Phi }(-\sqrt{2}))$ and $\sqrt{2}\varphi (\sqrt{2}(y-1))/(1-\mathrm{\Phi }(-\sqrt{2}))$, are plotted as red curves to enable comparison and evaluate the accuracy of the numerical approximations.

The selection of the step size presents a challenge. When the number of iterations is relatively small, the step size should be chosen carefully to balance accuracy and stability. If the step size is too large, such as ${\lambda }_k={(1+k)}^{-0.1}$, the discretized process $(U_n,V_n)$, as defined in Section 4.2, may frequently visit $(0,0)$, resulting in an approximating distribution that is overly concentrated at $(0,0)$. Moreover, as illustrated in Figure 2, the marginal distributions exhibit a good fit except at zero. Conversely, if the step size is too small, for instance, ${\lambda }_k={(1+k)}^{-1}$, the discretized process $(U_n,V_n)$ tends to remain near the initial position $(1,1)$, leading to an approximating distribution concentrated at $(1,1)$ and thus producing an inaccurate result (see Figure 4). In contrast, when the step size is moderate, the approximating distribution is more dispersed. As suggested by the marginal distribution diagrams in Figure 3, we believe this choice provides a closer approximation to the true stationary distribution.

We proceed to numerically estimate the correlation of the 2-d ROU processes under different parameter settings, using the scheme ${\lambda }_k={(1+k)}^{-0.7}$. Specifically, we use the empirical correlation

to approximate the true correlation of the 2-d ROU process. As before, the initial position is set to $(x_0,y_0)=(1,1)$, and the number of iterations is $n=200,000$. With extensive simulation experiments, we observe that the empirical correlations converge for the choice $n=200,000$. Furthermore, because the means and variances of the marginal distributions can be explicitly computed from their density functions, we omit their numerical results. The numerical results for the correlation are summarized in Tables 1–4, as well as in Figure 5. These results strongly suggest that the correlation of the 2-d ROU process is nondecreasing with respect to the correlation of the Brownian motion terms. A formal investigation of this property is left for future work.

Table 1. Numerical Results for the Correlation of the 2-d ROU Process for Different Values of $\rho$ When ${\theta }_1={\theta }_2=1$, ${\mu }_1={\mu }_2=1$, ${\sigma }_1={\sigma }_2=1$, ${\lambda }_k={(1+k)}^{-0.7}$, and $n=200,000$

Table 2. Numerical Results for the Correlation of the 2-d ROU Process for Different Values of $\rho$ When ${\theta }_1={\theta }_2=1$, ${\mu }_1={\mu }_2=-1$, ${\sigma }_1={\sigma }_2=1$, ${\lambda }_k={(1+k)}^{-0.7}$, and $n=200,000$

Table 3. Numerical Results for the Correlation of the 2-d ROU Process for Different Values of $\rho$ When ${\theta }_1={\theta }_2=1$, ${\mu }_1=1$, ${\mu }_2=-1$, ${\sigma }_1=1$, ${\sigma }_2=\sqrt{2}$, ${\lambda }_k={(1+k)}^{-0.7}$, and $n=200,000$

Table 4. Numerical Results for the Correlation of the 2-d ROU Process for Different Values of ρ When ${\theta }_1={\theta }_2=1$, ${\mu }_1={\mu }_2=-1$, ${\sigma }_1=1$, ${\sigma }_2=\sqrt{2}$, ${\lambda }_k={(1+k)}^{-0.7}$, and $n=200,000$

We close this subsection by approximating the Laplace transform of the stationary distribution $\pi$, namely, $\int e^{-sx-ty}\hspace{0.17em}\pi (dx,dy)$, employing the empirical Laplace transform

We adopt the parameter setting: ${\theta }_1={\theta }_2=1$, ${\mu }_1={\mu }_2=-1$, ${\sigma }_1=1$, ${\sigma }_2=\sqrt{2}$, $\rho =0.5$, with initial conditions $x_0=1$ and $y_0=1$. We employ the scheme ${\lambda }_k={(1+k)}^{-0.7}$ and set the number of iterations to $n=200,000$. For the parameters s and t, we consider the values $\{0.1,0.5,1,2,10\}$ and compute the corresponding results, which are summarized in Table 5.

Table 5. Numerical Results for the Laplace Transform When ${\theta }_1={\theta }_2=1$, ${\mu }_1={\mu }_2=-1$, ${\sigma }_1=1$, ${\sigma }_2=\sqrt{2}$, $\rho =0.5$, ${\lambda }_k={(1+k)}^{-0.7}$, and $n=200,000$