2. Preliminaries

Let $(\Omega, \mathbb{F}, \{\mathbb{F}_t\}, \mathbb{P})$ be a complete probability space with filtration $\{\mathbb{F}_t\}_{t\geq0}$ satisfying the common conditions (i.e. it is right continuous and increasing while $\mathbb{F}_0$ contains all $\mathbb{P}$-null sets). Let $(r_t)_{t\geq0}$ be a continuous-time Markov chain on the state space $\mathbb{S}=\{1, 2, \ldots, N\}$, with $N<\infty$, satisfying

\begin{eqnarray*}\mathbb{P}\{r_{t+\Delta}=j \mid r_t=i\}=\begin{cases}q_{ij}\Delta+o(\Delta) ,& \text{if}\ i\neq j,\\[5pt] 1+q_{ii}+o(\Delta) ,& \text{if}\ i=j,\end{cases}\end{eqnarray*}

where $\Delta>0$. The Q-matrix $Q=(q_{ij})_{N\times N}$ is assumed to be irreducible and conservative. Here, $q_{ij}\geq0$ is the transition rate from i to j if $i\neq j$, and $q_{ii}=-\sum_{i\neq j}q_{ij}$. Since the Markov chain $(r_t)$ is irreducible and conservative, it has a unique stationary distribution $\pi=(\pi_1, \pi_2, \ldots, \pi_N)$. The stationary distribution $\pi$ is the solution of the following equation:

\begin{equation*}\pi Q=0 \quad \text{and} \quad \sum_{i=1}^N\pi_i=1.\end{equation*}

Next, we consider the solution process of the SDE (1). Although $(X_t)$ is not a Markov process, the process $(X_t, r_t)$ is a Markov process. The infinitesimal generator $\mathscr{A}$ of $(X_t,r_t)_{t\geq0}$ corresponding to the the equation (1) is given by

(3) \begin{align} \nonumber \mathscr{A}V(x,i) & = x(a_{i}-b_{i}x) V_x(x,i)+QV(x,\cdot)(i) \\ & \quad + \int_{\mathbb{R}\setminus\{0\}}\big(V(x+\sigma_i x z,i) - V(x)-V_x(x,i) \sigma_i x z 1_{\{|z|\leq 1\}}\big)\nu(\mathrm{d} z,i),\end{align}

where $V_x(x,i)=\frac{\partial}{\partial x}V(x, i)$ and $QV(x, \cdot)(i)=\sum_{j\neq i}q_{ij}\left(V(x, j)-V(x, i)\right)$. Here, $V(\cdot, \cdot)\in C^2(\mathbb{R}\times \mathbb{S})$ means, for all $i\in \mathbb{S}$, $V(x,\cdot)\in C^2(\mathbb{R})$.

For any Borel set A of $\mathbb{R}^+$, $x\in\mathbb{R}^+$, $i, j\in\mathbb{S}$, set

\begin{equation*}\mathbb{P}(t, (x,i), (A,j))=\mathbb{P}(X_t\in A, r_t=j \mid X_0=x, r_0=i).\end{equation*}

Lemma 1. Under Assumptions 1–4, for any given initial data $(x_0, r_0)\in\mathbb{R}^+\times\mathbb{S}$, the SDE (1) has a unique positive solution $(X_t,r_t)$:

\begin{equation*} \mathbb{P}\{X_t\in\mathbb{R}^+\ \text{for all}\ t\geq0\}=1. \end{equation*}

Proof. The proof of the existence and uniqueness of the solution is almost the same as those of [20, Theorem 2.1] and [19, Lemma 3.1]. The sketch proof is given as follows. Rewrite the Markov chain as $r_t=r_0+\sum_{n=1}^\infty M_n I(\tau_n\leq t)$, where $\tau_n$, $M_n$ have the following conditional distributions. Given $r(\tau_k)=i$, $\tau_{k+1}-\tau_k$ follows the exponential distribution with parameters $q_{kk}$, the jump $M_k=r(\tau_{k+1})-r(\tau_k)$ is independent of the past, and $r_{\tau_{k+1}}$ hits the state j with probability $-\frac{q_{kj}}{q_{kk}}$. According to [20, Theorem 2.1], there exists a unique solution to the equation

\begin{equation*} \mathrm{d} Y_t^i = Y_t^i[(a_{i}-b_iX_t) \, \mathrm{d} t + \sigma_i \, \mathrm{d} Z_t^i], \qquad Y_0^i=y. \end{equation*}

For each i, denote the unique solution by $\{Y_t^{y,i}\}_{t\geq 0}$. With a little abuse of notation, define, on $t\in [0,\tau_1)$, $X_t=Y_t^{X_0,r_0}$. Iterating for n times, for $\tau_{n-1}\leq t<\tau_n$, we define $X_t=Y_t^{X_{\tau_{n-1}},r_{\tau_{n-1}}}$. It is not hard to see that, for any $t\in \mathbb{R}^+$,

(4) \begin{align} X_t & = X_0+\sum_{k=1}^\infty Y_t^{Y_{\tau_{n-1}},r_{\tau_{n-1}}}I(\tau_{k-1}\leq t<\tau_k) \nonumber \\ & = X_0+\int_0^tX_s(a_{r_s}-b_{r_s}X_s) \, \mathrm{d} s + \int_0^t \sigma_{r_s}X_s \, \mathrm{d} Z_s^{r_s}. \end{align}

Then, the construction solution $\{X_t\}_{t\geq 0}$ is the unique solution to the SDE (1).

To prove the non-negativity of the solution, let $t_0$ be any jump moment of the stable process $(Z_t)$. First, we prove that no jump causes the state of $X_t$ to be negative. Rewriting (1), we have, for any $r_{t_0^-}=j$, $\smash{X_{t_0}=X_{t_0^-}+\sigma_j X_{t_0^-}\triangle Z_{t_0}^j}$, where $\smash{\triangle Z_{t_0}^j=Z_{t_0}^j-Z_{t_0^-}^j}$ represents a jump of the stable process $\smash{(Z_t^j)}$ at time $t_0^-$. Using Assumption 4, we have $\smash{\triangle Z_{t_{0}}^j>-\frac{1}{\sigma_j}}$. If $\smash{X_{t_0^-}>0}$, then, almost surely (a.s.), $\sigma_j X_{t_0^-}\triangle Z_{t_{0}}^j>-X_{t_0^-}$. It follows from this that $X_{t_0}>0$ a.s.

Next, we prove that if there is no jump for the stable noise between $[t, t+\Delta t)$, then the continuous motion corresponding to the drift term cannot cause the path of $X_{t+\triangle t}$ to be negative if $X_{t}>0$. To prove the assertion, note that the continuous motion satisfies the SDE

(5) \begin{equation}\mathrm{d} X_t = X_t(a_{r_t}-b_{r_t}X_t) \, \mathrm{d} t, \qquad X_0=x>0 . \end{equation}

According to the construction of the solution (4), to prove the non-negativity of the solution to the SDE (5), we only need to prove that, for each $i\in \mathbb{S}$, the solution to the SDE

(6) \begin{equation}\mathrm{d} X_t^i = X_t^i(a_{i}-b_{i}X_t^i) \, \mathrm{d} t, \qquad X_0^i=x>0\end{equation}

is non-negative. Observe that (6) is the famous logistic differential equation. Let $Y_t=\frac{1}{X_t^i}$. Then, (6) can be rewritten as $\mathrm{d} Y_t = (b_i-a_iY_t) \, \mathrm{d} t$.

For any $a_i>0$, $b_i>0$, $X_0=x>0$, $t\in \mathbb{R}^+$, $Y_t=\frac{b_i}{a_i}+\big(\frac{1}{x}-\frac{b_i}{a_i}\big)\mathrm{e}^{-a_it}=\frac{b_i}{a_i}(1-\mathrm{e}^{-a_it})+\frac{1}{x}\mathrm{e}^{-a_it}>0$. Thus, $X_t^i=\frac{1}{Y_t}>0$ for all $t\in \mathbb{R}^+$. This, together with (4), implies that the solution $\{X_t\}$ to the SDE (5) is non-negative. This completes the proof.

Finally, we provide a modified definition of ergodicity.

3. Ergodicity

In this section we will consider the ergodicity of the population dynamics model (1). Our methods depend heavily on a modification of Khasminskii’s results.

Proof. The proof is similar to that of [17, Theorem 4.1] and [1, Theorem 5.1]. From condition (ii), the process $(X_t,r_t)$ is positive recurrent with respect to the domain $U\times \mathbb{S}$, and recurrent in the domain $\mathbb{R}^+\times \mathbb{S}$. By the maximum principle associated with the operator in (3) (see Appendix A), we get the existence of a stationary distribution. The proof of uniqueness is presented in Appendix A.

Theorem 1. For $\alpha\in(1,2)$, let Assumptions 1–4 hold. If

\[\sum_{i=1}^{N}\pi_{i}\left(a_{i}-\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}\right)>0,\]

then the SDE (1) is ergodic and has a unique stationary distribution.

Proof. We will prove this theorem according to Lemma 2. We will verify that the conditions (i) and (ii) hold in our setting. Let $N_1$ be a sufficiently large number and set $U=\big\{x\in \mathbb{R}^+\,:\,\frac{1}{N_1}<x<N_1\big\}$. Since $\sigma_i>0$, $i\in\mathbb{S}$, $\inf_{x\in{U}}\sigma_{i}^2 x^2=\frac{\sigma_\mathrm{min}^2}{N^2}>0$, where $\sigma_\mathrm{min}^2=\min\big\{\sigma_{1}^2, \sigma_{2}^2, \ldots, \sigma_{N}^2 \big\}$. This verifies condition (i) in Lemma 2.

Next, we will verify condition (ii) in Lemma 2. For simplicity, the column vector $\mu=(\mu_1, \mu_2, \ldots, \mu_N)^\top\in\mathbb{R}^N$ is defined as

\[\mu_i=\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}-a_{i}.\]

Also, let

\[\beta=-\pi\mu=\sum_{i=1}^N\pi_i\left(a_{i}-\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}\right).\]

Note that $\beta>0$. Let $\textbf{1}=(1,\ldots,1)^\top$; then $\pi(\mu+\beta \textbf{1})=-\beta+\beta=0$. From [16, Lemma A.12], this implies that the linear equation $Qc=\mu+\beta \textbf{1}$ has a solution $c=(c_1, c_2, \ldots, c_N)^\top\in\mathbb{R}^N$. Thus, we have

\begin{equation*} \beta=\sum_{j=1}^{N}q_{ij}c_{j}-\mu_i,\qquad i\in\mathbb{S}. \end{equation*}

For each $i\in\mathbb{S}$, consider the Lyapunov function $V\,:\,\mathbb{R}^+\times \mathbb{S}\rightarrow\mathbb{R}$, $V(x,i)=(1-\gamma c_i)(x^{\gamma}+x^{-\gamma})$, where $\gamma\in(0,\alpha-1)$ is sufficiently small that $1-\gamma c_i>0$ for each $i\in\mathbb{S}$. Hence, it follows that

\begin{align*} \mathscr{A}V(x,i) & = (a_{i}-b_{i}x)(1-\gamma c_i)\gamma (x^{\gamma}-x^{-\gamma})-\sum_{j\neq i}q_{ij}\gamma(c_j-c_i)(x^{\gamma}+x^{-\gamma}) \\ & \quad + (1-\gamma c_i)\int_{|z|\leq\frac{1}{\sigma_i}}\Big({(x+\sigma_i x z)^{\gamma}+(x+\sigma_i x z)^{-\gamma}-x^{\gamma}-x^{-\gamma}}\\ & \qquad \qquad \qquad \qquad \qquad - \gamma \sigma_i x zI_{\{|z|\leq K_i\}} \big(x^{\gamma-1}-x^{-\gamma-1}\big)\Big)\frac{C_{\alpha}}{|z|^{1+\alpha}} \, \mathrm{d} z, \end{align*}

where $K_i=\left\{ \begin{array}{cc} \frac{1}{\sigma_i}, & \sigma_i\geq 1, \\[4pt] 1, & \sigma_i<1. \\ \end{array} \right.$

By the variable substitution $y=\sigma_i z$, set

(7) \begin{equation} \mathscr{A}V(x,i) \,:\!=\, A_1(x,i)+A_2(x,i)+A_3(x,i), \end{equation}

where

\begin{align*} A_1(x,i) & = (a_{i}-b_{i}x)(1-\gamma c_i)\gamma (x^{\gamma}-x^{-\gamma})-\sum_{j\neq i}q_{ij}\gamma(c_j-c_i)(x^{\gamma}+x^{-\gamma}), \\[3pt] A_2(x,i) & = (1-\gamma c_i)x^{\gamma}\sigma^\alpha_i\int_{|y|\leq1}\left((1+y)^{\gamma}-1-\gamma yI_{\{|y|\leq \sigma_iK_i\}}\right)\frac{C_{\alpha}}{|y|^{1+\alpha}} \, \mathrm{d} y, \\[6pt] A_3(x,i) & = (1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\int_{|y|\leq1}\left((1+y)^{-\gamma}-1+\gamma yI_{\{|y|\leq \sigma_iK_i\}}\right)\frac{C_{\alpha}}{|y|^{1+\alpha}} \, \mathrm{d} y. \\ \end{align*}

First, we calculate the term $A_1(x,i)$:

\begin{align*} A_1(x,i) & = (a_{i}-b_{i}x)(1-\gamma c_i)\gamma (x^{\gamma}-x^{-\gamma})-\sum_{i\neq j}q_{ij}\gamma(c_j-c_i)(x^{\gamma}+x^{-\gamma}) \\ & = \gamma (1-\gamma c_i)\left[\bigg(a_{i}-\sum^N_{j=1}q_{ij}c_j\bigg)x^{\gamma}-b_{i}x^{\gamma+1} +b_{i}x^{1-\gamma}-\bigg(a_{i}+\sum^N_{j=1}q_{ij}c_j\bigg)x^{-\gamma}+o(\gamma)\right]. \end{align*}

To calculate $A_2(x,i)$, the following approximation is needed:

\begin{align*} \sum_{i\neq j}q_{ij}\frac{c_j-c_i}{1-\gamma c_i} & = \sum_{i\neq j}q_{ij}c_j+\sum_{i\neq j}q_{ij}\frac{\gamma c_ic_j-c_i}{1-\gamma c_i} \\ & = \sum_{j=1}^{N}q_{ij}c_j+\sum_{i\neq j}q_{ij}\frac{\gamma c_j(c_j-c_i)}{1-\gamma c_i}\gamma \\ & = \sum_{j=1}^{N}q_{ij}c_j+o(\gamma). \end{align*}

Then, $A_2(x,i)$ is computed as

\begin{align*} A_2(x,i) & = (1-\gamma c_i)x^{\gamma}\sigma^\alpha_i\int_{-1}^{0}\left((1+y)^{\gamma}-1-\gamma yI_{\{|y|\leq \sigma_iK_i\}})\right)\frac{C_{\alpha}}{|y|^{1+\alpha}} \, \mathrm{d} y \\[3pt] & \quad + (1-\gamma c_i)x^{\gamma}\sigma^\alpha_i\int_{0}^{1}\left((1+y)^{\gamma}-1-\gamma yI_{\{|y|\leq \sigma_iK_i\}})\right)\frac{C_{\alpha}}{|y|^{1+\alpha}} \, \mathrm{d} y \\[3pt] & = (1-\gamma c_i)x^{\gamma}\sigma^\alpha_i\int_{0}^{1}\left((1-y)^{\gamma}+(1+y)^{\gamma}-2\right)\frac{C_{\alpha}}{y^{1+\alpha}} \, \mathrm{d} y. \end{align*}

Define $F_1(y)=(1-y)^{\gamma}+(1+y)^{\gamma}-2$, where $y\in(0,1)$. By the Taylor expansion of the function $F_1(y)$ in the interval (0,1), $F_1(y)=2\sum^{+\infty}_{n=1}C^{2n}_{\gamma}y^{2n}.$ Thus, it follows that

(8) \begin{align} \nonumber A_2(x,i) & = (1-\gamma c_i)x^{\gamma}\sigma^\alpha_i\int_{0}^{1}2\sum^{+\infty}_{n=1}C^{2n}_{\gamma} y^{2n}\frac{C_{\alpha}}{y^{1+\alpha}} \, \mathrm{d} y \\[5pt] & = (1-\gamma c_i)x^{\gamma}\sigma^\alpha_i\frac{2\sum^{+\infty}_{n=1}C^{2n}_{\gamma}C_{\alpha}}{2n-\alpha}. \end{align}

As for the term $A_{3}(x,i)$, splitting the integral into two terms, we have

\begin{align*} A_3(x,i) & = (1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\int_{-1}^{1}\left((1+y)^{-\gamma}-1+\gamma yI_{\{|y|\leq \sigma_iK_i\}})\right)\frac{C_{\alpha}}{|y|^{1+\alpha}} \, \mathrm{d} y \\[5pt] & \,=\!:\, A_{31}(x,i)+A_{32}(x,i), \end{align*}

where

\begin{align*} A_{31}(x,i) & = (1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\int_{-1}^{1}\left((1+y)^{-\gamma}-1\right)\frac{C_{\alpha}}{|y|^{1+\alpha}} \, \mathrm{d} y, \\[5pt] A_{32}(x,i) & = (1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\int_{-1}^{1}\gamma yI_{\{|y|\leq \sigma_iK_i\}}\frac{C_{\alpha}}{|y|^{1+\alpha}} \, \mathrm{d} y. \end{align*}

To estimate $A_{31}(x,i)$, by the Taylor expansion for the function $(1+y)^{-\gamma}-1$ for $y\in (-1, 1)$, it follows that

\begin{align*} A_{31}(x,i) & = (1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\sum_{n=1}^{\infty} C_{-\gamma}^{n}\int_{-1}^{1}\frac{y^{n}C_{\alpha}}{|y|^{1+\alpha}} \, \mathrm{d} y \\[5pt] & = (1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\sum_{n=1}^{\infty} C_{-\gamma}^{2n}\int_{-1}^{1}\frac{y^{2n-1}C_{\alpha}\text{sgn}(y)} {|y|^{\alpha}} \, \mathrm{d} y \\[5pt] & \quad + (1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\sum_{n=1}^{\infty} C_{-\gamma}^{2n-1}\int_{-1}^{1}\frac{y^{2n-1}C_{\alpha}} {|y|^{1+\alpha}} \, \mathrm{d} y. \end{align*}

Note that

\begin{equation*} \int_{-1}^{1}\frac{y^{2n-1}C_{\alpha}}{|y|^{1+\alpha}}dy=0,\qquad \text{for\ all}\ n\in \mathbb{N}^+. \end{equation*}

Hence, we have

\begin{align*} A_{31}(x,i) & = (1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\sum_{n=1}^{\infty}C_{-\gamma}^{2n}\int_{-1}^{1}\frac{y^{2n-1}C_{\alpha}\mathrm{sgn}(y)} {|y|^{\alpha}} \, \mathrm{d} y \\ & = 2(1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\sum_{n=1}^{\infty}C_{-\gamma}^{2n}\int_{0}^{1}y^{2n}\frac{C_{\alpha}}{y^{1+\alpha}} \, \mathrm{d} y \\ & = 2(1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\sum_{n=1}^{\infty}C_{-\gamma}^{2n}\frac{C_{\alpha}}{2n-\alpha}. \end{align*}

Observing that $A_{32}(x,i)=(1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\int_{-1}^{1}\gamma yI_{\{|y|\leq \sigma_iK_i\}}\frac{C_{\alpha}}{|y|^{1+\alpha}} \, \mathrm{d} y=0$, we get

(9) \begin{equation} A_3(x,i) = 2(1-\gamma c_i)x^{-\gamma}\sigma^\alpha_i\sum_{n=1}^{\infty}C_{-\gamma}^{2n}\frac{C_{\alpha}}{2n-\alpha}. \end{equation}

Substituting the equalities (8) and (9) into (7) gives

\begin{align*} \mathscr{A}V(x,i) & = \gamma(1-\gamma c_i)\left[\left(a_{i}-\sum_{j=1}^{N}q_{ij}c_j+\sum^{+\infty}_{n=1}C^{2n}_{\gamma}\frac{2\sigma^\alpha_iC_{\alpha}}{\gamma(2n-\alpha)}\right)x^{\gamma} -b_{i}x^{\gamma+1}\right. \\[5pt] & \left. \quad + b_{i}x^{1-\gamma}-\left(a_{i}+\sum_{j=1}^{N}q_{ij}c_j-\sum_{n=1}^{\infty}C_{-\gamma}^{2n}\frac{2C_{\alpha}\sigma^\alpha_i}{\gamma(2n-\alpha)}\right)x^{-\gamma}\right]. \end{align*}

Letting $\gamma\rightarrow0$, we obtain

\begin{align*} \lim\limits_{\gamma\rightarrow0}\sum^{+\infty}_{n=1}C^{2n}_{\gamma}\frac{2\sigma^\alpha_iC_{\alpha}}{\gamma(2n-\alpha)} & = -\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}, \\[5pt] \lim\limits_{\gamma\rightarrow0}\sum^{+\infty}_{n=1}C^{2n}_{-\gamma}\frac{2C_{\alpha}\sigma^\alpha_i}{\gamma(2n-\alpha)} & = \sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}. \end{align*}

It follows from these limits that

(10) \begin{align} \nonumber \mathscr{A}V(x,i) & = \gamma(1-\gamma c_i)\left(a_{i}-\sum_{j=1}^{N}q_{ij}c_j-\sum_{n=1}^{N}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}+o(\gamma)\right)x^{\gamma} \\[3pt] & \quad - (1-\gamma c_i) b_{i}x^{\gamma+1}+(1-\gamma c_i)b_{i}x^{1-\gamma} \nonumber \\[3pt] & \quad - (1-\gamma c_i)\gamma\left(a_{i}-\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}+\sum_{j=1}^{N}q_{ij}c_j+o(\gamma)\right)x^{-\gamma} \nonumber \\[3pt] & = \gamma(1-\gamma c_i)\left(a_{i}-\sum_{j=1}^{N}q_{ij}c_j-\sum_{n=1}^{N}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}+o(\gamma)\right)x^{\gamma} \nonumber \\[3pt] & \quad - (1-\gamma c_i) b_{i}x^{\gamma+1}+(1-\gamma c_i)b_{i}x^{1-\gamma} \nonumber\\[3pt] & \quad - (1-\gamma c_i)\gamma\left(\beta+o(\gamma)\right)x^{-\gamma}. \end{align}

Therefore, there exists a sufficiently large $N_1>0$ such that, for all $x>N_1$ and $x<\frac{1}{N_1}$, $i\in \mathbb{S}$, it follows that

(11) \begin{equation} \mathscr{A}V(x,i)\leq-1. \end{equation}

Again, by equality (10), for any $x\in \mathbb{R}^+$, $i\in \mathbb{S}$, there exists a positive constant C such that $\mathscr{A}V(x,i)\leq C$. By Itô’s formula [Reference Applebaum2], we get

\begin{align*} V(X_{t},r_{t})=V(X_{0},r_{0})+\int_{0}^{t}\mathscr{A}V(X_{s},r_{s}) \, \mathrm{d} s+M_t, \end{align*}

where

\begin{equation*} M_t=\int_0^t\left(V(X_s+\sigma_{r_s}X_sz,r_s)-V(X_s,r_s)- V_x(X_s,r_s)\sigma_{r_s}X_szI_{\{|z|<1\}}\right)\widetilde{N}^{r_s}(\mathrm{d} s,\mathrm{d} z). \end{equation*}

Note that the process $(M_t)$ is a local martingale. Let $T_n$ be a localizing sequence, such that $M_{t\wedge T_n}$ is a martingale in t for any fixed n. By the Itô formula again, we have

\begin{equation*} V(X_{t\wedge T_n\wedge \tau_U},r_{t\wedge T_n\wedge \tau_U})=V(X_{0},r_{0})+\int_{0}^{t\wedge T_n \wedge \tau_U}\mathscr{A}V(X_{s},r_{s}) \, \mathrm{d} s+M_{t\wedge T_n\wedge \tau_U}. \end{equation*}

For any $(x_0,r_0)\in U^c\times \mathbb{S}$, it follows that

(12) \begin{equation}\mathbb{E} V(X_{t\wedge T_n\wedge \tau_U},r_{t\wedge T_n\wedge \tau_U})=V(X_{0},r_{0})+ \mathbb{E} \int_{0}^{t\wedge T_n \wedge \tau_U}\mathscr{A}V(X_{s},r_{s}) \, \mathrm{d} s.\end{equation}

Thus, applying the dominated convergence theorem yields

(13) \begin{eqnarray}\lim\limits_{n\rightarrow +\infty} \mathbb{E} \int_{0}^{t\wedge T_n \wedge \tau_U}\mathscr{A}V(X_{s},r_{s}) \, \mathrm{d} s=\mathbb{E} \int_{0}^{t\wedge \wedge \tau_U}\mathscr{A}V(X_{s},r_{s}) \, \mathrm{d} s.\end{eqnarray}

Substituting (13) into (12) gives

\begin{equation*} \lim\limits_{n\rightarrow +\infty} \mathbb{E} V(X_{t\wedge T_n\wedge \tau_U},r_{t\wedge T_n\wedge \tau_U})=V(X_{0},r_{0})+ \mathbb{E} \int_{0}^{t \wedge \tau_U}\mathscr{A}V(X_{s},r_{s}) \, \mathrm{d} s.\end{equation*}

It follows from Fatou’s lemma that

\begin{align*}\mathbb{E} V(X_{t\wedge\tau_{U}},r_{t\wedge\tau_{U}}) & = \mathbb{E} \lim\limits_{n\rightarrow +\infty} V(X_{t\wedge T_n\wedge\tau_{U}}, r_{t\wedge T_n\wedge\tau_{U}}) \\ & = \mathbb{E} \liminf\limits_{n\rightarrow +\infty} V(X_{t\wedge T_n\wedge\tau_{U}}, r_{t\wedge T_n\wedge\tau_{U}}) \\ & \leq \liminf\limits_{n\rightarrow +\infty} \mathbb{E} V(X_{t\wedge T_n\wedge\tau_{U}}, r_{t\wedge T_n\wedge\tau_{U}}) \\ & = \lim\limits_{n\rightarrow +\infty} \mathbb{E} V(X_{t\wedge T_n\wedge\tau_{U}}, r_{t\wedge T_n\wedge\tau_{U}}) \\ & = V(x,i)+ \mathbb{E} \int_0^{t\wedge\tau_{U}}\mathscr{A}V(X_s,r_s) \, \mathrm{d} s, \end{align*}

which together with (11) implies that, for all $(X_{0}=x,r_{0}=i)\in U^c\times\mathbb{S}$, $\mathbb{E} (\tau_{U} \mid X_0=x_0, r_0=i_0)< V(x,i)<+\infty$. Thus, condition (ii) of Lemma 2 is verified. It follows from Lemma 2 that the desired result holds.

Theorem 1 provides a sufficient condition to check the ergodicity of the model (1) for $1<\alpha<2$. It is a natural question whether, when $\alpha$ tends to 2, this condition simplifies the condition established in [Reference Mao13] for the Brownian motion case. The following corollary will answer this question.

Corollary 1. Let Assumptions 1–4 hold. If $\alpha$ tends to 2, then the condition

\[\sum_{i=1}^{N}\pi_i\left(a_i-\sum_{n=1}^{\infty}\frac{C_{\alpha}\sigma^\alpha_i}{n(2n-\alpha)}\right)>0\]

reduces to $\sum_{i=1}^{N}\pi_{i}(a_{i}-\sigma_{i}^{2})>0$.

Proof. Using the fact that $\lim\limits_{x\rightarrow 0^+}x\Gamma(x) = 1$ yields $\lim\limits_{\alpha\rightarrow 2^-}C_{\alpha}=0$, $\lim\limits_{\alpha\rightarrow 2^-}\frac{C_{\alpha}}{2-\alpha} = 1$. Therefore, it follows that, for any $i\in \mathbb{S}$,

\begin{equation*} \lim\limits_{\alpha\rightarrow2}\sum_{n=1}^{N}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)} = \sigma^2_i.\end{equation*}

This implies the desired result.

Although for truncated stable processes corresponding to the Lévy jump measure given by Assumption 4 a sufficient condition is given in Theorem 1, it is worth considering what happens if the jump measure in Assumption 4 is replaced by Assumption 5. Next, we consider the jump process $\Delta Z_t\,=\!:\,Z_t-Z_{t^-}$ under Assumption 5.

Theorem 2. Let Assumptions 1–3 and 5 hold. If

\[\sum_{i=1}^N\pi_i\left(a_{i} -\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}+\sigma_i^\alpha \frac{C_\alpha}{\alpha}\right)>0,\]

then the SDE model (1) is ergodic and has a unique stationary distribution.

Proof. The proof is similar that of Theorem 1. For simplicity, we only provide the part where it differs. Using the Taylor formula we have, for $y\in (1,+\infty)$, $\gamma>0$,

\begin{equation*} (1+y)^{\gamma} =\sum\limits_{j=0}^{+\infty}C_\alpha C_\gamma^j{y^{-j+\gamma}}, \qquad (1+y)^{-\gamma} =\sum\limits_{j=0}^{+\infty}C_\alpha C_\gamma^j{y^{-j-\gamma}}. \end{equation*}

Hence, it follows that

(14) \begin{align} I_1\,:\!=\,\int_1^{+\infty} \left((1+y)^\gamma -1\right)\frac{C_\alpha \,\mathrm{d} y}{y^{1+\alpha}} & = \sum\limits_{j=0}^{+\infty}C_\gamma^j \frac{C_\alpha}{j+\alpha-\gamma} -\frac{C_\alpha}{\alpha} , \end{align}

(15) \begin{align} I_2\,:\!=\,\int_1^{+\infty} \left((1+y)^{-\gamma} -1\right)\frac{C_\alpha \,\mathrm{d} y}{y^{1+\alpha}} & = \sum\limits_{j=0}^{+\infty}C_{-\gamma}^j \frac{C_\alpha}{j+\alpha+\gamma} -\frac{C_\alpha}{\alpha}. \end{align}

Rewriting the limits in (14) and (15), we get $I_1=I_2=-\frac{C_\alpha}{\alpha}+o(\gamma)$.

Set the column vector $\overline{\mu}=(\overline{\mu}_1, \ldots, \overline{\mu}_N)^\top\in\mathbb{R}^N$ with

\[\overline{\mu}_i=\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}-\sigma_i^\alpha \frac{C_\alpha}{\alpha}-a_{i}.\]

Define

\[\overline{\beta}=-\pi\overline{\mu}=\sum_{i=1}^N\pi_i\left(a_{i} -\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}+\sigma_i^\alpha \frac{C_\alpha}{\alpha}\right).\]

By the conditions of Theorem 2, we have $\overline{\beta}>0$. Using the same discussion as for Theorem 1 yields $\overline{\beta}=\sum_{j=1}^{N}q_{ij}c_{j}-\overline{\mu}_i$, $i\in\mathbb{S}$, where $c=(c_1,\ldots, c_N)$ is a solution of the Poisson equation $Qc=\overline{\mu}+\overline{\beta} \textbf{1}$.

Let V(x, i) be the Lyapunov function defined in Theorem 1. Similar to the discussion of Theorem 1, and together with the equalities (14) and (15), it follows that

\begin{align*} \mathscr{A}V(x,i) & = \gamma(1-\gamma c_i)\left(a_{i}-\sum_{j=1}^{N}q_{ij}c_j-\sum_{n=1}^{N}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}-\sigma_1^\alpha\frac{C_\alpha}{\alpha}+o(\gamma)\right)x^{\gamma} \\ & \quad -(1-\gamma c_i) b_{i}x^{\gamma+1}+(1-\gamma c_i)b_{i}x^{1-\gamma} \\ & \quad -(1-\gamma c_i)\gamma\left(\overline{\beta}+o(\gamma)\right)x^{-\gamma}. \end{align*}

Since $\overline{\beta}>0$ and $b_i>0$, it follows from Lemma 1 that the desired result holds.

Next, we will rewrite Theorem 2 to compare with the case $\alpha=2$. Let $D_\alpha= \sum^{+\infty}_{n=2}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}-\sigma_i^\alpha \frac{C_\alpha}{\alpha}$.

Proof. From Theorem 2, we can immediately see that

\begin{equation*} a_{i} -\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}+\sigma_i^\alpha \frac{C_\alpha}{\alpha}=a_i-\frac{C_{\alpha}\sigma^\alpha_i} {2-\alpha}-D_\alpha. \end{equation*}

Note that

(16) \begin{equation} \sum_{n=2}^\infty\frac{1}{n(2n-\alpha)}=-\frac{2}{\alpha} \sum_{n=2}^{\infty}\left(\frac{1}{2n}-\frac{1}{2n-\alpha}\right). \end{equation}

For $1<\alpha<2$, it follows that

(17) \begin{equation} \sum_{n=2}^{\infty}\left(\frac{1}{2n}-\frac{1}{2n-\alpha}\right)<-\frac{1}{2}. \end{equation}

Equations (16) and (17) yield that, for any $i\in \mathbb{S}$,

(18) \begin{equation} D_\alpha=\sigma_i^\alpha C_\alpha\left(-\frac{1}{\alpha}+\sum_{n=2}^\infty\frac{1}{n(2n-\alpha)}\right)>0. \end{equation}

Consequently, Theorem 7 combined with (18) implies the desired result.

From Corollary (2), letting $\alpha$ tend to 2 gives the following result.

Proof. By the proof of Corollary 1, we have $\lim\limits_{\alpha\rightarrow 2^-}\frac{C_\alpha}{2-\alpha}=1$. It follows from this and Corollary 7 that the desired result holds.

4. Extinction

In this section we find sufficient conditions such that the solution of the SDE (1) tends to zero almost surely. In the terminology of population dynamics, we shall derive under what conditions the population dynamics will go extinct.

Proof. The basic idea is to adopt the Lyapunov method together with the Poisson equation. We construct a Lyapunov function $V(x,i)\,:\,\mathbb{R}^+\times \mathbb{S}\rightarrow \mathbb{R}$, that is,

\begin{eqnarray*} \lim\limits_{x\rightarrow 0^+} V(x,i)=0 ,\qquad \lim\limits_{x\rightarrow +\infty}V(x,i)=+\infty.\end{eqnarray*}

Then, we prove that, for any $(x,i)\in \mathbb{R }^+\times \mathbb{S}$, $\lim\limits_{t\rightarrow +\infty} \mathbb{E} [V(X_t,r_t) \mid X_0=x, r_0=i]=0$, which implies the desired result. For simplicity, we use the same notation as the proof of Theorem 1. Recall that $\mu=(\mu_1, \mu_2, \ldots, \mu_N)^\top\in\mathbb{R}^N$, with $\mu_i=\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i}{n(2n-\alpha)}-a_i$. Thus,

\[-\beta=\pi\mu=-\sum_{i=1}^N\pi_i\left(a_{i}-\sum^{+\infty}_{n=1}\frac{C_{\alpha}\sigma^\alpha_i}{n(2n-\alpha)}\right).\]

By the conditions of Theorem 3, we have $\beta<0$. Let $\textbf{1}=(1,\ldots,1)^\top$. Then, $\pi(-\mu-\beta \textbf{1})=\beta-\beta=0$. Again from [16, Lemma A.12], this implies that the linear equation $Qc=-\mu-\beta \textbf{1}$ has a solution $c=(c_1, c_2, \ldots, c_N)^\top\in\mathbb{R}^N$. Thus, we have

\begin{eqnarray*}\beta=-\mu_i-\sum_{j=1}^{N}q_{ij}c_{j},\qquad i\in\mathbb{S}.\end{eqnarray*}

For each $i\in\mathbb{S}$, consider the Lyapunov function $V\,:\,\mathbb{R}^+\times \mathbb{S}\rightarrow\mathbb{R}$ such that $V(x,i)=(1-\gamma c_i)x^{\gamma}$, where $\gamma\in(0,\alpha-1)$ is sufficiently small that $1-\gamma c_i>0$ for each $i\in\mathbb{S}$. Using the same discussion as in Theorem 1, we obtain

\begin{align*} \mathscr{A}V(x,i) & = \gamma(1-\gamma c_i)\left(a_{i}-\sum_{j=1}^{N}q_{ij}c_j-\sum_{n=1}^{N}\frac{C_{\alpha}\sigma^\alpha_i} {n(2n-\alpha)}+o(\gamma)\right)x^{\gamma} \\ & \quad -(1-\gamma c_i) b_{i}x^{\gamma+1} \\ & = \gamma(1-\gamma c_i)\left(\beta+o(\gamma)\right)x^{\gamma}-(1-\gamma c_i) b_{i}x^{\gamma+1}.\end{align*}

Note that $\beta<0$. It follows that, for any $(x,i)\in \mathbb{R}^+\times \mathbb{S}$, $\mathscr{A}V(x,i)<0$. Similar to (12) and (13), it can be shown that

\begin{eqnarray*} 0\leq \mathbb{E} V(X_{t},r_{t})\leq V(x,i)+\mathbb{E}\int_0^{t}\mathscr{A}V(X_s,r_s) \, \mathrm{d} s.\end{eqnarray*}

By the non-negativity of V(x, i), we immediately obtain the desired result.

Next, we reconsider Theorem 3 if Assumption 4 is replaced by Assumption 5.

Proof. The proof of this theorem is the same as Theorems 3 and 2, and is not detailed here.

Appendix A

To prove Lemma 2, we first prove a maximum principle for operators defined by (3). Assume for the following that $U\subset \mathbb{R}^+$ is open and bounded, with closure $\overline{U}\subset \mathbb{R}^+$. Let $U^c=\mathbb{R}^+-U$. For simplicity, assume that $U=\big(\frac{1}{N_1}, N_1\big)$, where $N_1>0$. Consider the equation

(19) \begin{eqnarray} \left\{ \begin{array}{cl} \mathscr{A}u(x,i)=0, & \mbox{in}\ \ (x,i)\in U\times \mathbb{S} , \\[5pt] u(x,i)=c_1, & \mbox{in}\ \ (x,i)\in (N_1,+\infty)\times \mathbb{S} , \\[5pt] u(x,i)=c_2, & \mbox{in}\ \ (x,i)\in (0,\frac{1}{N_1})\times \mathbb{S} , \\ \end{array} \right.\end{eqnarray}

where $c_1$, $c_2$ are two non-negative constants and $c_1\neq c_2$.

Proof. For simplicity, we divide the proof into two steps.

Step 1. Assume that $\mathscr{A}u(x,i)>0$ for all $(x,i)\in U\times \mathbb{S}$. Let us prove that the maximum of u cannot be attained in the interior. If $u(\cdot,\cdot)$ attains a maximum at $(x_0,i_0)\in U\times \mathbb{S}$, then

(20) \begin{eqnarray}u_x(x_0,i_0)=0.\end{eqnarray}

Since the measure $\nu(\cdot,i)$ is symmetric, we can rewrite the integral term on the right-hand side of (3) as

\begin{eqnarray*}I(x,i)\,:\!=\,\mbox{p.v.} \int_{\mathbb{R}-\{0\}}\left(u(x+\sigma_ixz,i)-u(x,i)\right)\nu(\mathrm{d} z,i),\end{eqnarray*}

where $\mbox{p.v.}$ stands for the Cauchy principal value.

Noting that $u(x_0,i_0)$ is the maximum value on the domain $U\times \mathbb{S}$, we have

(21) \begin{eqnarray}I(x_0,i_0)= \int_{U}\left(u(x_0+\sigma_ix_0z,i)-u(x_0,i_0)\right)\nu(\mathrm{d} z,i_0)\leq 0.\end{eqnarray}

For the term, we have

(22) \begin{align} Qu(x_0,\cdot)(i_0) & = q_{i_0i_0}u(x_0,i_0)+\sum\limits_{j\neq i_0}q_{i_0j}u(x_0,j) \nonumber \\ & \leq q_{i_0i_0}u(x_0,i_0)+\sum\limits_{j\neq i_0}q_{i_0j} u(x_0,i_0)=0.\end{align}

It follows from (3) and (20)–(22) that $\mathscr{A}u(x_0,i_0)\leq 0$. This is a contradiction with $\mathscr{A}u(x_0,i_0)>0$.

Step 2. Now we consider u such that $\mathscr{A}u(x,i)\geq 0$ for all $x\in U$, $i\in \mathbb{S}$. To proceed, we consider two cases.

Case 1: Assumptions 1–4 hold. We define the function $u_\varepsilon(x,i)=u(x,i)+\varphi(x,i)$, where $\varphi(\cdot,i)\in C^2(\mathbb{R}^+)$ and $\varphi(x,i)=K x^2$, where K is a positive constant to be determined. Since $\overline{U}$ is a bounded and closed set and $b(\cdot,i)$ is continuous on $\overline{U}$, there exists a constant $C>0$ such that, for any $x\in \overline{U}, i\in \mathbb{S}$, $|2xb(x,i)|<C$. Without loss of generality, we assume that $U=(\varepsilon_0, 1/\varepsilon_0)$, where $\varepsilon_0$ is any positive real number. Choosing $K=\frac{(C+1)(2-\alpha)}{2\sigma^2_i\varepsilon_0^2C_\alpha \sigma_i^{\alpha-2}}$, we get, for any $x\in U, i\in \mathbb{S}$,

\begin{align*} \mathscr{A}u_\varepsilon(x,i) & = \mathscr{A}u(x,i)+ b(x,i)\varphi_x(x,i) \\ & \quad +\int_{|z|<\frac{1}{\sigma_i}}(\varphi(x+\sigma_i xz,i)-\varphi(x,i)-\varphi_x(x,i)\sigma_ixz)\frac{C_\alpha \, \mathrm{d} z}{|z|^{1+\alpha}} \\ & = \mathscr{A}u(x,i)+2xb(x,i)+2K\sigma_i^2 x^2\frac{C_\alpha\sigma_i^{\alpha-2}}{2-\alpha} \\ & \geq \mathscr{A}u(x,i)-C+2K\sigma_i^2 x^2 \frac{C_\alpha\sigma_i^{\alpha-2}}{2-\alpha}\geq 1.\end{align*}

Now, the desired result follows from Step 1.

Case 2: Assumptions 1–3 and 5 hold. For any $\varepsilon>0$, define the function $u_\varepsilon(x,i)=u(x,i)+\varphi(x,i)$, where $\varphi(\cdot,i)\in C^2(\mathbb{R}^+)$ and

\begin{eqnarray*} \varphi(x,i)=\left\{ \begin{array}{l@{\quad}l} x^2, & \mbox{if} \ 0\leq x<2R, \\[3pt] 8R^2, &\mbox{if}\ x>2R+1,\\ \end{array} \right.\end{eqnarray*}

where R is a constant to be determined. Hence, we have, for any $x\in U$, $i\in \mathbb{S}$,

\begin{align*} \mathscr{A}u_\varepsilon(x,i) & = \mathscr{A}u(x,i)+b(x,i)\varphi_x(x,i) \\[2pt] & \quad + C_\alpha\sigma_i^\alpha x^\alpha\int_{B_R(0)\bigcap (-1,+\infty)}(\varphi(x+z,i)-\varphi(x,i)-\varphi_x(x,i)z) \frac{ \mathrm{d} z}{|z|^{1+\alpha}} \\[6pt] & \quad + C_\alpha\sigma_i^\alpha x^\alpha\int_{1\leq z\leq R}\!\!\!\!\!\! \varphi_x (x,i)\frac{ \mathrm{d} z}{|z|^{1+\alpha}} +\int_{B^c_{R}(0)\bigcap (-1,+\infty)}\!\!\!\!\!\!\!\!\!\!\!\!\left(\varphi(x+z)-\varphi(x)\right)\frac{ \mathrm{d} z}{|z|^{1+\alpha}}>0. \end{align*}

Note that, for $1<\alpha<2$, $x\in U$. For any $U=\big(\varepsilon_0, \frac{1}{\varepsilon_0}\big)$, the term

\begin{align*} I & \,:\!=\, C_\alpha\sigma_i^\alpha x^\alpha\int_{z>-1}(\varphi(x+z,i)-\varphi(x,i)-\varphi_x(x,i)z) \frac{ \mathrm{d} z}{|z|^{1+\alpha}} \\[5pt] & \geq 2\inf_{x\in U}C_\alpha \sigma_i^\alpha x^\alpha\int_{z>-1}z^2 \frac{\mathrm{d} z}{|z|^{1+\alpha}}=+\infty.\end{align*}

Set $K=\inf\limits_{x\in U, i\in \mathbb{S}}b(x,i)\varphi_x(x,i)$. Thus, there exists a sufficiently large R such that $I>|K|+1$. Since sufficiently large R, $x\in U\subset \mathbb{R}^+$, we have $\varphi(x+z,i)-\varphi(x,i)>0$ for all $z\in B_{R^c} \bigcap (-1,+\infty)$, and $\varphi_x>0$. Hence, it follows that $\mathscr{A}u_\varepsilon(x,i)>1$. Applying the result of Step 1, the maximum of $u_\varepsilon$ is attained on $\partial U \times \mathbb{S}$, say at $(x_\varepsilon, i)$. Since $\partial U$ is compact, choosing a convergent subsequence of $x_\varepsilon$ as $\varepsilon$ tends to 0 through a sequence, it follows that the maximum of u is attained on $\partial U$. Finally, if $\mathscr{A}u(x,i)=0$ for $x\in U$, $i\in \mathbb{S}$, then applying the above results to u and $-u$, we obtain that u attains both maximum and minimum on $\partial U \times \mathbb{S}$.

Next, we will prove the uniqueness of the stationary distribution of the process $(X_t,r_t)$. Lemma 2 can be viewed as a degenerate noise coefficient matrix of [Reference Chen, Chen, Tran and Yin5]. By the maximum principle, following almost the same argument as [17, Theorem 4.10] together with the solution to the Dirichlet problem (cf. [Reference Arapostathis, Biswas and Caffarelli1]), we can prove the existence of the stationary distribution of the process $(X_t,r_t)$. To proceed, we introduce some notation. Let (x, i) be the initial point of the process $(X_t,r_t)$, and $y\in \mathbb{R}^+$ be any point $(x\neq y)$. Set $k>\max\{|x|, |y|+1\}$. Without loss of generality, assume that $0<\varepsilon<1$. Choose c such that $0<c<k$ and $\overline{B}_c(0)\cap B_{\varepsilon}(y)=\emptyset$. Define the stopping times $\{\tau_{n}\}_{n\geq1}$ inductively as follows. Let $\tau_0=0$. For $n\in N^+$, set

\begin{align*}\tau_{2n+1} & = \inf\{s>\tau_{2n}\,:\, X_s \notin B_k(0) \} , \\\tau_{2n+2} & = \inf\{s>\tau_{2n-1}\,:\, (X_s,r_s) \ \mbox{hits\ the\ set}\ \overline{B}_c(0)\times \{i\} \}.\end{align*}

Proof of Lemma 2. We need to show the uniqueness of the invariant measure for the process $(X_t,r_t)$. This result is a modification of [17, Theorems 4.1, 4.2; Corollary 4.4]. To proceed, the positive recurrence gives $\mathbb{E}_{x,i}(\tau_{2n}-\tau_{2n-2})=\mathbb{E}_{x,i}\tau_2< +\infty$. Hence, by the classical strong law of large numbers,

\begin{eqnarray*}\frac{\tau_{2n}}{n}=\frac{\sum_{j=1}^n(\tau_{2n}-\tau_{2(n-1)})}{n}\rightarrow\mathbb{E}_{x,i}\tau_2,\end{eqnarray*}

almost surely.

Let f be a bounded real-valued measurable function on $\mathbb{R}^+\times \mathbb{S}$. For $n\in \mathbb{N}^+$, set $Y_n=\int_{\tau_{2n}}^{\tau_{2n+2}}f(X_s,r_s) \, \mathrm{d} s$. Similar to [4, Proposition 12.1], we can prove that the sequences of the random variables $\{Y_n\}$ are independent and identically distributed. Then, we can see that, applying the strong law of large numbers to the sequence $\{Y_n\}_{n=1}^\infty$, we have

\begin{eqnarray*}\lim_{n\rightarrow +\infty}\frac{1}{n}\sum_{j=1}^n Y_j = \mathbb{E}_{x, i} Y_0= \mathbb{E}_{x,i}\int_0^{\tau_2}f(X_s,r_s) \, \mathrm{d} s.\end{eqnarray*}

Similar to [4, Theorem 8.14], we have

\begin{eqnarray*}\lim\limits_{n\rightarrow +\infty}\frac{1}{\tau_{2n}}\int_0^{\tau_{2n}}f(X_s,r_s) \, \mathrm{d} s = \frac{1}{t}\int_0^t f(X_s,r_s) \, \mathrm{d} s.\end{eqnarray*}

Note that

\begin{eqnarray*}\frac{1}{\tau_{2n}}\int_0^{\tau_{2n}}f(X_s,r_s) \, \mathrm{d} s=\frac{n}{\tau_{2n}}\frac{1}{n}\sum_{j=1}^n Z_j\rightarrow \frac{1}{\mathbb{E}_{x, i}\tau_2} \mathbb{E}_{x,i}\int_0^{\tau_2} f(X_s,r_s) \, \mathrm{d} s\end{eqnarray*}

with probability 1. Hence, we have

(23) \begin{eqnarray}\mathbb{P}\left\{\lim\limits_{t\rightarrow +\infty}\frac{1}{t}\int_0^t f(X_s,r_s) \, \mathrm{d} s=\sum_{j\in \mathbb{S}}\int_{\mathbb{R}^+}f(y,i)\mu(\mathrm{d} y,j)\right\}=1.\end{eqnarray}

If $f(\cdot,i)$ is a bounded function with respect to x, then

\begin{eqnarray*}\frac{1}{t}\int_0^t \mathbb{E}_{x,i}f(X_s,r_s) \, \mathrm{d} s = \mathbb{E}_{x,i}\left(\frac{1}{t}\int_0^tf(X_s,r_s) \, \mathrm{d} s\right).\end{eqnarray*}

By (23) and the dominated convergence theorem, we have

\begin{eqnarray*}\lim\limits_{t\rightarrow +\infty} \frac{1}{t}\int_0^t \mathbb{E}_{x,i} f(X_s,r_s) \, \mathrm{d} s = \sum_{j\in \mathbb{S}}\int_{\mathbb{R}^+}f(y,j)\mu(\mathrm{d} y,j).\end{eqnarray*}

In particular, for any $A\in \mathbb{B}(\mathbb{R}^+)$, $j\in \mathbb{S}$, taking $f(x,j)=I_{A\times j}$, we have

\begin{equation*} \frac{1}{t}\int_0^t\mathbb{P}(s, (x,i), A\times j) \, \mathrm{d} s\rightarrow \mu(A, j).\end{equation*}

Last, we prove the uniqueness. Let $\mu_1(\cdot,\cdot)$ be another stationary distribution. Then

(24) \begin{eqnarray}\sum_{i\in \mathbb{S}}\int_{\mathbb{R}^+}\mathbb{P}(s,(x,i), (A,j))\mu_1(\mathrm{d} x,i)=\mu_1(A,j).\end{eqnarray}

Integrating (24) with respect to s from 0 to t, according to the limit (23), we have, for any $j\in \mathbb{S}$, $\mu_1(A,j)=\mu(A,j)$, which completes the proof.