Skew Ornstein–Uhlenbeck processes with sticky reflection and their applications to bond pricing

Abstract

We study a skew Ornstein–Uhlenbeck process with zero being a sticky reflecting boundary, which is defined as the weak solution to a stochastic differential equation (SDE) system involving local time. The main results obtained include: (i) the existence and uniqueness of solutions to the SDE system, (ii) the scale function and speed measure, and (iii) the distributional properties regarding the transition density and the first hitting times. On the application side, we apply the process to interest rate modeling and obtain the explicit pricing formula for zero-coupon bonds. Numerical examples illustrate the impacts on bond yields of skewness and stickiness parameters.

1. Introduction

Since the seminal works of [18] and [41], the skew Brownian motion and its diffusion extensions have long been the focus of much academic attention. Serving as a paradigm for developing the related theory, [17] showed that the skew Brownian motion is the unique solution to a stochastic differential equation (SDE) involving a local time term. Lots of techniques have been advanced for the in-depth study of skew Brownian motions; see, among others, [2, 3, 8, 11, 15]. There are also a number of papers devoted to the investigation of more general skew diffusions. For example, [5] introduced the skew Bessel process by determining the jumps of the local time in a certain way. By means of moving domains, [39] constructed solutions of the squared Bessel and Cox–Ingersoll–Ross (CIR) processes skewed on a time-dependent curve, and the corresponding pathwise uniqueness results were offered in [40]. The distributions of the skew Ornstein–Uhlenbeck (OU) process were probed in [42], and a series of explicit formulae obtained. However, compared to the abundant literature dealing with skew diffusions whose state space is the whole real line, there seem to be few results on skew processes with sticky reflecting boundaries.

The concept of the sticky reflecting boundary goes back to [16]. Different from other boundary behaviors, the set of times spent at such a boundary is nowhere dense and forms a Cantor set with positive Lebesgue measure. We refer the reader to [1, 13, 30, 32, 38] for more about the probabilistic properties of sticky reflecting diffusions (even with jumps). This kind of process has found a variety of applications, among which we mention its close connection with interest rate modeling in financial markets. As evidenced by [24], use of the sticky reflecting boundary is consistent with empirical observations on treasury-bill yields when near zero. In view of that finding, [27] checked the empirical performance of a two-dimensional sticky interest rate model, and [28] made a theoretical study of the Vasicek interest rate dynamics with a sticky zero lower bound. Using novel computational methods based on continuous-time Markov chain approximations, [25, 26] obtained bond prices under single- and multi-factor sticky interest rate models, respectively. On the other hand, [12] claimed that skew processes are suitable to characterize interest rates, where the skew parameters are related to uncertain monetary policies. The skew phenomenon was recently uncovered in [4] from interest rate data in the framework of skew OU diffusions. These in turn motivate us to utilize the skewed version of a sticky reflecting process for describing the interest rate.

In the coming sections, we seek to explore some basic facts about the skew OU process with sticky reflection and apply it as an interest rate model to bond pricing. Since the SDE system satisfied by this process has not been investigated before, we begin by establishing the existence and uniqueness of weak solutions in Section 2. We follow the ideas of [28, 30] for sticky reflecting OU and Bessel diffusions, while noting that some key points of their analysis cannot be extended to our case directly due to the existence of ‘skewness’. The trick to address this issue is to transform the skew processes in question into more tractable diffusions with discontinuous coefficients. Section 3 provides the infinitesimal generator, the scale function, and the speed measure, where the latter two can be derived by relating them to those of an auxiliary process with an instantaneous reflecting boundary and discontinuous coefficients. As a by-product, the explicit expression for the stationary distribution is obtained. Important distributional properties, including the Green function and the first-hitting-time problems, are studied in Section 4. In particular, inspired by [6, 29], we consider the first hitting time at a random jump boundary. The martingale approach is adopted to compute the expectation of the local time term, which plays a vital role in arriving at the final conclusion. In Section 5, we compute the price of a zero-coupon bond by using the spectral expansion approach that was proposed by [23] to solve derivative-pricing problems. Our interest rate model contains the one used by [28] as a special case and thus shows more flexibility. Numerical illustrations show that increase in the value of the skewness or stickiness parameter exerts a positive influence on the bond yield and the slope of the yield curve.

2. Existence and uniqueness of solutions

Let us take as given a filtered probability space $({\Omega},{\mathcal {F}},\{\mathcal {F}_{t},t\geq0\},{{\mathbb {P}}})$ which carries a standard Brownian motion $\{B_{t},t\geq 0\}$ with the filtration $\{\mathcal {F}_{t},t\geq0\}$ satisfying the usual conditions. For any continuous semimartingale $\{I_{t},t\geq 0\}$ defined on this space, denote by $\{\widehat{L}^{I}_{t}(x), t\geq 0\}$ (resp. $\{{L}^{I}_{t}(x), t\geq 0\}$ ) the symmetric (resp. right) semimartingale local time of I at the level x in the sense of [31, p. 212]. We consider on $[0,\infty)$ a diffusion process governed by the SDE system

with the constants $\theta\in\mathbb {R}$ , $\kappa>0$ , $\sigma>0$ , $a>0$ , $\beta>0$ , and $0<p<1$ . Above 0 the process X behaves like a skew OU process [42, 43] with the interesting feature that when hitting a it reflects upwards with probability p and downwards with probability $1-p$ , while after X touches the boundary 0 for the first time, it will spend positive amount of time at 0 but not be absorbed by 0 (see [44] for more details on this point), hence the name ‘sticky reflection’. The parameter $\beta$ dictates the intensity of stickiness; we may consider two extreme cases with $\beta\rightarrow 0$ corresponding to the absorbing boundary and $\beta\rightarrow \infty$ corresponding to the instantaneously reflecting boundary.

Before examining the probabilistic properties associated with X, we first present an existence and uniqueness result in the following theorem. See [20] for the standard definitions of weak solutions and uniqueness in law.

Theorem 2.1. The SDE system of (2.1) and (2.2) has a unique weak solution (X, B) for every arbitrary initial value $x_{0}\geq0$ . Also, the process X is strong Markov.

Proof. First, we show the existence of a weak solution to the system. Let us first consider on $[0,\infty)$ the skew OU process with zero as an instantaneously reflecting boundary:

(2.3) \begin{equation} {\textrm{d}}Y_{t} = \kappa(\theta-Y_{t})\,{\textrm{d}}t + \sigma\,{\textrm{d}}W_{t} + (2p-1)\,{\textrm{d}}\widehat{L}^{Y}_{t}(a) + \frac{1}{2}\,{\textrm{d}}L^{Y}_{t}(0), \end{equation}

with $\{W_{t},t\geq 0\}$ being a standard Brownian motion. Introduce the function

(2.4)

whose inverse function $N(\!\cdot\!)$ and symmetric derivative $g(\!\cdot\!)$ are given by

\begin{align*}N(x)=\left\{\begin{array}{l@{\quad}l} \frac{1}{p}(x-a)+a, & (1-p)a\leq x\leq a,\\[5pt] \frac{1}{1-p}(x-a)+a, & x>a,\end{array}\right.\end{align*}

\begin{align*} g(x)=\left\{\begin{array}{l@{\quad}l} p, & 0\leq x< a,\\[5pt] \frac{1}{2}, & x=a,\\[5pt] 1-p, & x>a,\end{array}\right.\end{align*}

Defining a new process $\{Z_{t}\triangleq G(Y_{t}),t\geq 0\}$ and putting

\begin{align*} \widetilde{{\mu}}(x) & = \left\{\begin{array}{l@{\quad}l} \kappa[p\theta+(1-p)a-x], & (1-p)a\leq x\leq a, \\[5pt] \kappa[(1-p)\theta+pa-x], & x>a, \end{array} \right. \end{align*}

\begin{align*} \widetilde{{\sigma}}(x) & = \left\{\begin{array}{l@{\quad}l} p\sigma, & (1-p)a\leq x\leq a, \\[5pt] (1-p)\sigma, & x>a, \end{array} \right. \end{align*}

we have, by the generalized Itô formula [14, p. 150],

\begin{align*} {\textrm{d}}Z_{t} & = \kappa(\theta-Y_{t})g(Y_{t})\,{\textrm{d}}t + \sigma g(Y_{t})\,{\textrm{d}}W_{t} + \frac{p}{2}\,{\textrm{d}}L^{Y}_{t}(0) \\[5pt] & = \kappa(\theta-N(Z_{t}))g(N(Z_{t}))\,{\textrm{d}}t + \sigma g(N(Z_{t}))\,{\textrm{d}}W_{t} + \frac{1}{2}\,{\textrm{d}}L^{Z}_{t}((1-p)a) \\[5pt] & = \widetilde{\mu}(Z_{t})\,{\textrm{d}}t + \widetilde{\sigma}(Z_{t})\,{\textrm{d}}W_{t} + \frac{1}{2}\,{\textrm{d}}L^{Z}_{t}((1-p)a), \end{align*}

where for the second equality we used the relation ${p} L^{Y}_{t}(0)= L^{Z}_{t}((1-p)a)$ due to [5, Lemma A.5], and the last equality holds almost surely because the occupation times formula [35, p. 224, Corollary 1.6] implies that Z has no occupation time in a. It is straightforward to see that the drift and diffusion coefficients of Z meet the following three conditions: [(i)]

1. (i) There exists a positive constant K such that $\widetilde{{\mu}}^{2}(x)+\widetilde{{\sigma}}^{2}(x)\leq K(1+x^{2})$ for any $x\geq (1-p)a$ .

2. (ii) $[\widetilde{{\sigma}}(x)-\widetilde{{\sigma}}(y)]^{2}\leq |f(x)-f(y)|$ for any $x,y\geq (1-p)a$ in which $f(x)\triangleq(1-2p)^{2}\sigma^{2}\textbf{1}_{\{x>a\}}$ is a bounded increasing function.

3. (iii) $|\widetilde{{\sigma}}(x)|\geq \min\{p,1-p\}\sigma$ for any $x\geq (1-p)a$ .

Hence, it follows from [37, Theorem 4.2] and [36, Theorem 1.3] that the SDE with discontinuous coefficients and instantaneous reflection satisfied by Z allows a unique strong solution, which further suggests that Y is the unique strong solution to the SDE (2.3) due to the one-to-one correspondence between Y and Z.

Now consider the random time

(2.5) \begin{equation} \eta_{t} = t + \frac{1}{2\beta}L^{Y}_{t}(0), \end{equation}

which is a strictly increasing and continuous stochastic process admitting a proper inverse $\xi_{t}=\eta_{t}^{-1}$ . We define a time-changed version of Y by $\{X_{t}\triangleq Y_{\xi_{t}},t\geq 0\}$ . It can be readily inferred from the definitions of semimartingale local times as occupation time densities [31, p. 225, Corollary 3] that $\widehat{L}^{X}_{t}(a)=\widehat{L}^{Y}_{\xi_{t}}(a)$ and $L^{X}_{t}(0)=L^{Y}_{\xi_{t}}(0)$ , and thus

(2.6) \begin{equation} \xi_{t}=t-\frac{1}{2\beta}L^{X}_{t}(0) \end{equation}

due to (2.5). This in turn implies the SDE (2.2) since

\begin{align*} \int_{0}^{t}\textbf{1}_{\{X_{s}=0\}}\,{\textrm{d}}s & = \int_{0}^{t}\textbf{1}_{\{X_{s}=0\}}\,{\textrm{d}}\xi_{s} + \frac{1}{2\beta}\int_{0}^{t}\textbf{1}_{\{X_{s}=0\}}\,{\textrm{d}}L^{X}_{s}(0) \\[5pt] & = \int_{0}^{\xi_{t}}\textbf{1}_{\{Y_{s}=0\}}\,{\textrm{d}}s + \frac{1}{2\beta}L^{X}_{t}(0) = \frac{1}{2\beta}L^{X}_{t}(0). \end{align*}

On the other hand, in light of (2.2), (2.6), and the definition of X, we obtain

(2.7) \begin{align} X_{t} & = Y_{0} + \int_{0}^{t}\kappa(\theta-Y_{\xi_{s}})\,{\textrm{d}}\xi_{s} + \sigma W_{\xi_{t}} + (2p-1)\widehat{L}^{Y}_{\xi_{t}}(a) + \frac{1}{2}L^{Y}_{\xi_{t}}(0) \nonumber \\[5pt] & = Y_{0} + \int_{0}^{t}\kappa(\theta-X_{s})\,{\textrm{d}}\bigg[s-\frac{1}{2\beta}L^{X}_{s}(0)\bigg] + \sigma W_{\xi_{t}} + (2p-1)\widehat{L}^{Y}_{\xi_{t}}(a) + \frac{1}{2}L^{Y}_{\xi_{t}}(0) \nonumber \\[5pt] & = X_{0} + \int_{0}^{t}\kappa(\theta-X_{s})\textbf{1}_{\{X_{s}>0\}}\,{\textrm{d}}s + \sigma W_{\xi_{t}} + (2p-1)\widehat{L}^{X}_{{t}}(a) + \frac{1}{2}L^{X}_{{t}}(0). \end{align}

Noting that $\{W_{\xi_{t}},t\geq 0\}$ is a continuous martingale with the quadratic variation process $\{\langle W_{\xi}\rangle_{t}=\int_{0}^{t}\textbf{1}_{\{X_{s}>0\}}\,{\textrm{d}}s,t\geq 0\}$ , we may use the same logic as in the proof of [20, Theorem 5.5.4] to express $W_{\xi_{t}}$ as a Brownian integral. In fact, the martingale representation theorem indicates that there exist a Brownian motion $\{\widetilde{W}_{t},t\geq 0\}$ and an adapted process $\{A_{t},t\geq 0\}$ defined on a possibly extended probability space $(\widetilde{\Omega},\widetilde{\mathcal {F}},\widetilde{{\mathbb {P}}})$ such that $W_{\xi_{t}} = \int_{0}^{t}A_{s}\,{\textrm{d}}\widetilde{W}_{s}$ , $\langle W_{\xi}\rangle_{t} = \int_{0}^{t}A^{2}_{s}\,{\textrm{d}}s$ . Consequently,

(2.8) \begin{align} W_{\xi_{t}} = \int_{0}^{t}\textbf{1}_{\{X_{s}>0\}}\,{\textrm{d}}B_{s}, \end{align}

where $B_{t} = \int_{0}^{t}{\textrm{sgn}}(\textbf{1}_{\{X_{s}>0\}}A_{s})\,{\textrm{d}}\widetilde{W}_{s}$ is a Brownian motion by Lévy’s theorem. Here. sgn(x) stands for the value of the sign function at x, which is equal to 1 if $x>0$ and $-1$ otherwise. Substituting (2.8) back into (2.7), we arrive at the SDE (2.1).

Next, we demonstrate the uniqueness in law of solutions to the SDE system. To start with, we refresh all the notation used earlier, and let X and B solve the SDEs (2.1) and (2.2). Define $\xi_{t} = \int_{0}^{t}\textbf{1}_{\{X_{s}>0\}}\,{\textrm{d}}s$ and its right continuous inverse $\eta_{t} = \inf\{s\geq0\colon\xi_{s}>t\}$ . Obviously we have $\eta_{\xi_{\infty}}=\infty$ .

Consider the process $\{{W}_{t}=\int_{0}^{\eta_{t}}\textbf{1}_{\{X_{s}>0\}}\,{\textrm{d}}B_{s},t\geq 0\}$ . It is a Brownian motion with the observations that $\langle W\rangle_{t}=t$ and that [20, Problem 3.4.5, assertion (iv)] guarantees the path continuity of W. As a result, the time-changed process $\{{Y}_{t}=X_{\eta_{t}},t\geq 0\}$ satisfies

(2.9) \begin{align} Y_{t} & = X_{0} + \int_{0}^{\eta_{t}}\kappa(\theta-X_{s})\,{\textrm{d}}\xi_{s} + \sigma\int_{0}^{\eta_{t}}\textbf{1}_{\{X_{s}>0\}}\,{\textrm{d}}B_{s} + (2p-1)\widehat{L}^{X}_{\eta_{t}}(a) + \frac{1}{2}L^{X}_{\eta_{t}}(0) \nonumber \\[5pt] & = Y_{0} + \int_{0}^{t}\kappa(\theta-Y_{s})\,{\textrm{d}}{s} + \sigma W_{t} + (2p-1)\widehat{L}^{Y}_{{t}}(a) + \frac{1}{2}L^{Y}_{{t}}(0) \end{align}

on the set $\{t<\xi_{\infty}\}$ . Also, by (2.2) and the definition of $\eta$ ,

(2.10) \begin{equation} \eta_{t} = \int_{0}^{\eta_{t}}\textbf{1}_{\{X_{s}>0\}}\,{\textrm{d}}s + \int_{0}^{\eta_{t}}\textbf{1}_{\{X_{s}=0\}}\,{\textrm{d}}s = t + \frac{1}{2\beta}L^{X}_{\eta_{t}}(0) = t + \frac{1}{2\beta}L^{Y}_{{t}}(0), \end{equation}

from which we know that $\eta$ is strictly increasing and continuous, and that $\eta_{\xi_{\infty}}<\infty$ on $\{\xi_{\infty}<\infty\}$ , contradicting the fact $\eta_{\xi_{\infty}}=\infty$ unless the event $\{\xi_{\infty}=\infty\}$ is of probability 1. Hence, $\xi$ is the proper inverse of $\eta$ , and $Y_{\xi_{t}}=X_{t}$ accordingly. This, together with (2.10), confirms that the law of X is unique since Y and $\xi$ are uniquely determined by the standard Brownian motion W.

Finally, to establish the strong Markov property of X, we first transform it into an SDE without skewness. Specifically, set $\overline{{X}}_{t}=G(X_{t})$ , recalling (2.4), and an application of the generalized Itô formula [14, p. 150] produces

(2.11) \begin{equation} {\textrm{d}}\overline{X}_{t} = \widetilde{\mu}(\overline{X}_{t})\textbf{1}_{\{\overline{X}_{t}>(1-p)a\}}\,{\textrm{d}}t + \widetilde{{\sigma}}(\overline{X}_{t})\textbf{1}_{\{\overline{X}_{t}>(1-p)a\}}\,{\textrm{d}}B_{t} + \frac{1}{2}\,{\textrm{d}}L^{\overline{X}}_{t}((1-p)a). \end{equation}

Moreover, by (2.2) and [5, Lemma A.5],

(2.12) \begin{equation} \textbf{1}_{\{\overline{X}_{t}=(1-p)a\}}\,{\textrm{d}}t = \textbf{1}_{\{{X}_{t}=0\}}\,{\textrm{d}}t = \frac{1}{2\beta}\,{\textrm{d}}L^{X}_{t}(0) = \frac{1}{2p\beta}\,{\textrm{d}}L^{\overline{X}}_{t}((1-p)a). \end{equation}

Therefore, the pair $(\overline{X},B)$ constitutes a weak solution to the SDE system of (2.11) and (2.12) that is unique in the law sense due to the one-to-one correspondence between $\overline{X}$ and X. Going one step further, we find that $\overline{X}$ solves

(2.13) \begin{equation} {\textrm{d}}\overline{X}_{t} = \widetilde{\mu}(\overline{X}_{t})\textbf{1}_{\{\overline{X}_{t}>(1-p)a\}}\,{\textrm{d}}t + \widetilde{{\sigma}}(\overline{X}_{t})\textbf{1}_{\{\overline{X}_{t}>(1-p)a\}}\,{\textrm{d}}B_{t} + p\beta\textbf{1}_{\{\overline{X}_{t}=(1-p)a\}}\,{\textrm{d}}t \end{equation}

by plugging (2.12) into (2.11). Conversely, if $\overline{X}$ solves (2.13), then [31, Corollary 1, p. 224] suggests that

\begin{align*} L^{\overline{X}}_{t}((1-p)a) - 2p\beta\int_{0}^{t}\textbf{1}_{\{\overline{X}_{s}=(1-p)a\}}\,{\textrm{d}}s & = L^{\overline{X}}_{t}((1-p)a-\!) \\[5pt] & = \lim_{\varepsilon\rightarrow 0}\frac{1}{\varepsilon} \int_{0}^{t}\textbf{1}_{\{(1-p)a-\varepsilon\leq\overline{X}_{s}\leq(1-p)a\}}\,{\textrm{d}}\langle\overline{X}\rangle_{s} \\[5pt] & = \lim_{\varepsilon\rightarrow 0}\frac{1}{\varepsilon} \int_{0}^{t}\textbf{1}_{\{(1-p)a-\varepsilon\leq\overline{X}_{s}\leq(1-p)a\}}\widetilde{{\sigma}}^{2}(\overline{X}_{s}) \textbf{1}_{\{\overline{X}_{s}>(1-p)a\}}\,{\textrm{d}}s \\[5pt] & = 0. \end{align*}

The relation in (2.12) thus holds, and (2.11) follows by substitution. In other words, the SDE system consisting of (2.11) and (2.12) is equivalent to the SDE (2.13). Thanks to [20, Theorem 5.4.20], we conclude that $\overline{X}$ solving (2.13) is strong Markov, and so is the process X. The proof is complete.

In view of the above theorem, we call the process X the skew OU process with sticky reflection (or with sticky reflecting boundary) hereafter.

3. Fundamental characterizations

Based on the results derived in the last section and according to the formulation by [7] (see also [18]), the skew OU process with sticky reflection is a linear diffusion that can be fully characterized through its infinitesimal generator and the corresponding domain of definition. Specifically, the infinitesimal generator $\mathscr{A}$ of X governed by (2.1) and (2.2) is given by

(3.14) \begin{equation} \mathscr{A} = \frac{\sigma^{2}}{2}\frac{{\textrm{d}}^{2}}{{\textrm{d}}x^{2}} + \kappa(\theta-x)\frac{{\textrm{d}}}{{\textrm{d}}x}\end{equation}

whose domain $\mathfrak{D}(\mathscr{A})$ is restricted to $\mathcal{C}_\textrm{b}([0,\infty))$ , the set of continuous bounded functions defined on $[0,\infty)$ . As implied by [21, 28], we have

\begin{align*} \mathfrak{D}(\mathscr{A}) = \{f\in\mathcal{C}_\textrm{b}([0,\infty))\cap\mathcal{C}^{2}((0,\infty)\backslash\{a\})\colon \\[5pt] \mathscr{A}f\in\mathcal{C}_\textrm{b}([0,\infty)), pf'(a+\!) = (1-p)f'(a-\!),(\mathscr{A}f)(0) = \beta f'(0+)\}.\end{align*}

Every linear diffusion (without killing) has two basic characteristics: scale function and speed measure. As will be seen, among others, the scale density (the derivative of the scale function) of X is discontinuous at the point a if $p\neq \frac12$ , and the speed measure puts positive mass at the boundary 0.

Proposition 3.1. Let

(3.15)

(3.16)

Then, for the process X, the scale function and the speed measure , where $\delta_{0}({\textrm{d}}x)$ is the Dirac measure with unit mass concentrated at 0.

Proof. Consider the auxiliary processes Y and Z as defined in the proof of Theorem 2.1. Using standard results [35, p. 311, Exercise 3.20], we obtain the scale and speed densities of Z given by

Next, we try to make clear the relationship between $\mathfrak{s}_{Z}(\!\cdot\!)$ and $\mathfrak{s}_{Y}(\!\cdot\!)$ , the scale density of Y. To this end, let $\tau^{I}_{b} = \inf\{t\geq0\colon I_{t}=b\}$ be the first hitting time at level b for a process $\{I_{t},t\geq 0\}$ . On one hand, we have, by the definition of the scale function [35, p. 301, Proposition 3.2], for $0\leq b_{1}<y<b_{2}$ ,

(3.17)

On the other hand, it follows from the relation $\tau^{Y}_{b}=\tau^{Z}_{G(b)}$ that

(3.18)

As a consequence, the combination of (3.17) and (3.18) leads to

(3.19) \begin{equation} \mathfrak{s}_{Y}(x)=G'(x)\mathfrak{s}_{Z}(G(x)). \end{equation}

To figure out the relationship between $\mathfrak{m}_{Z}(\!\cdot\!)$ and the speed density $\mathfrak{m}_{Y}(\!\cdot\!)$ of Y, we need the function defined by

\begin{equation*} J^{I}_{(z_{1},z_{2})}(x,y) = \frac{[\mathbb m{s}_{I}(x\wedge y)-\mathbb m{s}_{I}(z_{1})][\mathbb m{s}_{I}(z_{2})-\mathbb m{s}_{I}(x\vee y)]} {\mathbb m{s}_{I}(z_{2})-\mathbb m{s}_{I}(z_{1})}, \end{equation*}

with I being the diffusion Y or Z. Evidently, $J^{Y}_{(z_{1},z_{2})}(x,y)=J^{Z}_{(G(z_{1}),G(z_{2}))}(G(x),G(y))$ , which together with [35, Theorem 7.3.6] suggests that

\begin{align*} {\mathbb {E}}_{y}\big[\tau^{Y}_{b_{1}}\wedge\tau^{Y}_{b_{2}}\big] & = {\mathbb {E}}_{G(y)}\big[\tau^{Z}_{G(b_{1})}\wedge\tau^{Z}_{G(b_{2})}\big] \\[5pt] & = \int_{G(b_{1})}^{G(b_{2})}J^{Z}_{(G(b_{1}),G(b_{2}))}(G(y),v)\mathfrak{m}_{Z}(v)\,{\textrm{d}}v \\[5pt] & = \int_{b_{1}}^{b_{2}}J^{Z}_{(G(b_{1}),G(b_{2}))}(G(y),G(u))\mathfrak{m}_{Z}(G(u))G'(u)\,{\textrm{d}}u \\[5pt] & = \int_{b_{1}}^{b_{2}}J^{Y}_{(b_{1},b_{2})}(y,u)\mathfrak{m}_{Z}(G(u))G'(u)\,{\textrm{d}}u \end{align*}

for $0\leq b_{1}<y<b_{2}$ . This gives

(3.20) \begin{equation} \mathfrak{m}_{Y}(x)=G'(x)\mathfrak{m}_{Z}(G(x)), \end{equation}

in that is the unique Radon measure such that

For the process X, recall from the proof of Theorem 2.1 that it can be constructed as a time-changed process $Y_{\xi_{t}}$ . Then, the fact that $\tau^{X}_{b}=\eta_{\tau^{Y}_{b}}$ implies that for $0\leq b_{1}<x<b_{2}$ , thereby indicating that $\mathfrak{s}_{X}(x)=\mathfrak{s}_{Y}(x)$ , followed by the relation $J^{X}_{(b_{1},b_{2})}(x,y)=J^{Y}_{(b_{1},b_{2})}(x,y)$ . Furthermore, noting that

\begin{equation*} {\mathbb {E}}_{x}\big[\tau^{X}_{b_{1}}\wedge\tau^{X}_{b_{2}}\big] = {\mathbb {E}}_{x}\bigg[\bigg(\tau^{Y}_{b_{1}}+\frac{1}{2\beta}L^{Y}_{\tau^{Y}_{b_{1}}}(0)\bigg) \wedge \bigg(\tau^{Y}_{b_{2}}+\frac{1}{2\beta}L^{Y}_{\tau^{Y}_{b_{2}}}(0)\bigg)\bigg] = {\mathbb {E}}_{x}\big[\tau^{Y}_{b_{1}}\wedge\tau^{Y}_{b_{2}}\big], \end{equation*}

we arrive at $\mathfrak{m}_{X}(x)=\mathfrak{m}_{Y}(x)$ . The desired results (3.15) and (3.16) thus follow by direct calculation from (3.19) and (3.20). Finally, a comparison between [7, formula (c), p. 16] and the boundary condition $(\mathscr{A}f)(0)=\beta f'(0+)$ for $f(\!\cdot\!)\in \mathfrak{D}(\mathscr{A})$ yields $\mathbb m{m}_{X}(\{0\})=1/(p\beta)$ , which completes the proof.

Corollary 3.1. The process X is positively recurrent with the unique stationary distribution . Here, $\gamma$ is a constant defined by

\begin{equation*} \gamma = \frac{(1-p)\sigma\sqrt{\kappa}} {(1-p)\sigma\sqrt{\kappa} + 2\beta\sqrt{\pi}\textrm{e}^{{\kappa\theta^{2}}/{\sigma^{2}}} [p+(1-2p)\Phi({(a-\theta)\sqrt{2\kappa}}/{\sigma}) -(1-p)\Phi(\!-{\theta\sqrt{2\kappa}}/{\sigma})]}, \end{equation*}

where $\Phi(\!\cdot\!)$ denotes the standard normal cumulative distribution function, and $\mathfrak{r}(x)=p\beta\gamma \mathfrak{m}_X(x)$ for all $x>0$ , where $\mathfrak{m}_X(\!\cdot\!)$ is given in (3.16).

Proof. From Proposition 3.1 it is clear that $\mathbb m{s}_{X}(\infty)=\infty$ , which signifies that X is a recurrent diffusion with sticky reflection. In fact, when $x>b\geq0$ , for any $b_{0}>x$ ,

\begin{equation*} {\mathbb {P}}_{x}\big(\tau^{X}_{b}<\infty\big) \geq {\mathbb {P}}_{x}\big(\tau^{X}_{b}<\tau^{X}_{b_{0}}\big) = \frac{\mathbb m{s}_{X}(b_{0})-\mathbb m{s}_{X}(x)}{\mathbb m{s}_{X}(b_{0})-\mathbb m{s}_{X}(b)}, \end{equation*}

suggesting that ${\mathbb {P}}_{x}\big(\tau^{X}_{b}<\infty\big)=1$ by letting $b_{0}\rightarrow \infty$ . On the other hand, arguing as in the proof of [19, Theorem 23.15, p. 465], we obtain that ${\mathbb {P}}_{0}\big(\tau^{X}_{y}<\infty\big)=1$ for all $y>0$ . As a result, when $0\leq x< b$ , the strong Markov property of X generates

\begin{align*} {\mathbb {P}}_{x}\big(\tau^{X}_{b}<\infty\big) & = {\mathbb {P}}_{x}\big(\tau^{X}_{b}<\tau^{X}_{0}\big)+{\mathbb {P}}_{x}\big(\tau^{X}_{0}<\tau^{X}_{b}<\infty\big) \\[5pt] & ={\mathbb {P}}_{x}\big(\tau^{X}_{b}<\tau^{X}_{0}\big) + {\mathbb {P}}_{x}\big(\tau^{X}_{0}<\tau^{X}_{b}\big){\mathbb {P}}_{0}\big(\tau^{X}_{b}<\infty\big) =1. \end{align*}

Furthermore, since

\begin{align*} \mathbb m{m}_X([0,\infty)) & = \int_{0}^{\infty}\mathfrak{m}_{X}(x)\,{\textrm{d}}x + \mathbb m{m}_{X}(\{0\}) \\[5pt] & = \frac{1}{p\beta} + \frac{2}{p\sigma^{2}}\int_{0}^{a}{\textrm{e}}^{-({\kappa}/{\sigma^{2}})[(x-\theta)^{2}-\theta^{2}]}\,{\textrm{d}}x + \frac{2}{(1-p)\sigma^{2}}\int_{a}^{\infty}{\textrm{e}}^{-({\kappa}/{\sigma^{2}})[(x-\theta)^{2}-\theta^{2}]}\,{\textrm{d}}x \\[5pt] & = \frac{1}{p\beta\gamma}<\infty, \end{align*}

X is positively recurrent [7, p. 20], and the unique stationary distribution is given by $\mathbb m{r}({\textrm{d}}x) = [\mathfrak{m}_X(x)\,{\textrm{d}}x + \mathbb m{m}_{X}(\{0\})\delta_{0}({\textrm{d}}x)]/{\mathbb m{m}_X([0,\infty))}$ . The proof is thus complete.

4. Some distributional properties

So far we have worked out the scale density and speed measure of the skew OU process with sticky reflection. In this section, we turn to studying its distributional properties associated with the first hitting times and the transition density. More precisely, set $p(t;\;x,y)$ to be the transition density of X governed by the SDE system of (2.1) and (2.2) in the sense that ${\mathbb {P}}_{x}(X_{t}\in {\textrm{d}} y)=p(t;\;x,y)\mathbb m{m}_X({\textrm{d}} y)$ . We begin by computing the Laplace transform of the first hitting time $\tau^{X}_{y}$ to a constant level $y\geq 0$ , denoted by

(4.21) \begin{equation} \zeta_\lambda(x,y)={\mathbb {E}}_{x}\big[{\textrm{e}}^{-\lambda \tau^{X}_{y} }\big],\end{equation}

and deriving the explicit expression for the Green function

which is the Laplace transform of the transition density in the time domain. Actually, both of them can be expressed in terms of the functions $\phi_{\lambda}(\!\cdot\!)$ and $\psi_{\lambda}(\!\cdot\!)$ , as displayed in the following theorem. Let $H_{\nu}(\!\cdot\!)$ represent the Hermite function of degree $\nu\triangleq-\lambda/\kappa$ [22, p. 285], and $\Gamma(\!\cdot\!)$ the gamma function [22, p. 1]. It is well known that $H_{\nu}(\!\cdot\!)$ (resp. $H_{\nu}(\!-\cdot)$ ) is a decreasing (resp. increasing) and positive function on the real line which solves the ordinary differential equation (ODE) (4.31) with the asymptotic properties $H_{\nu}(\!-\infty)=\infty$ and $H_{\nu}(\infty)=0$ .

Theorem 4.1. Introduce the two functions

(4.22) \begin{align} \phi_{\lambda}(x) = \left\{\begin{array}{l@{\quad}l} H_{\nu}(z(x))+c_{1}H_{\nu}(\!-z(x)), & 0\leq x\leq a, \\[5pt] c_{2}H_{\nu}(z(x))+c_{3}H_{\nu}(\!-z(x)), & x>a, \\[5pt] \end{array}\right. \end{align}

(4.23) \begin{align} \psi_{\lambda}(x) = \left\{\begin{array}{l@{\quad}l} d_{1}H_{\nu}(z(x))+d_{2}H_{\nu}(\!-z(x)), & 0\leq x\leq a, \\[5pt] H_{\nu}(\!-z(x)), & x>a, \\[5pt] \end{array}\right. \end{align}

and the constant

\begin{equation*} w_{\lambda}=\frac{c_{2}\lambda 2^{\nu+1}\sqrt{\pi}{\textrm{e}}^{z^{2}(0)}}{(1-p)\sigma\Gamma(1-\nu)\sqrt{\kappa}\,}. \end{equation*}

Here, $z(x)=-\sqrt{\kappa}(x-\theta)/\sigma$ , and the coefficients are given by

(4.24) \begin{align} c_{1} = \frac{2\beta H_{\nu-1}(z(0))-\sigma\sqrt{\kappa}H_{\nu}(z(0))} {2\beta H_{\nu-1}(\!-z(0))+\sigma\sqrt{\kappa}H_{\nu}(\!-z(0))}, \end{align}

(4.25) \begin{align} c_{2} = \frac{\Gamma(1-\nu){\textrm{e}}^{-z^{2}(a)}}{p2^{\nu}\sqrt{\pi}} [pK_{1}H_{\nu-1}(\!-z(a))+(1-p)K_{2}H_{\nu}(\!-z(a))], \end{align}

(4.26) \begin{align} c_{3} = \frac{\Gamma(1-\nu){\textrm{e}}^{-z^{2}(a)}}{p2^{\nu}\sqrt{\pi}} [pK_{1}H_{\nu-1}(z(a))-(1-p)K_{2}H_{\nu}(z(a))], \end{align}

(4.27) \begin{align} d_{1} = \frac{(1-2p)\Gamma(1-\nu){\textrm{e}}^{-z^{2}(a)}}{(1-p)2^{\nu}\sqrt{\pi}}[H_{\nu-1}(\!-z(a))H_{\nu}(\!-z(a))], \end{align}

(4.28) \begin{align} d_{2} = \frac{\Gamma(1-\nu){\textrm{e}}^{-z^{2}(a)}}{(1-p)2^{\nu}\sqrt{\pi}} [(1-p)H_{\nu-1}(z(a))H_{\nu}(\!-z(a))+pH_{\nu-1}(\!-z(a))H_{\nu}(z(a))], \end{align}

where $K_{1} = H_{\nu}(z(a)) + c_{1}H_{\nu}(\!-z(a))$ , $K_{2} = H_{\nu-1}(z(a)) - c_{1}H_{\nu-1}(\!-z(a))$ . Then we have, for $x,y\geq 0$ and $\lambda>0$ ,

(4.29) \begin{align} \zeta_\lambda(x,y) = \left\{\begin{array}{l@{\quad}l} \dfrac{\phi_{\lambda}(x)}{\phi_{\lambda}(y)}, & x\leq y, \\[9pt] \dfrac{\psi_{\lambda}(x)}{\psi_{\lambda}(y)}, & x>y, \end{array} \right.\end{align}

(4.30)

Proof. The general theory of linear diffusions [7, pp. 18–19] tells us that the expression in (4.29) holds, where $\phi_{\lambda}(\!\cdot\!)$ and $\psi_{\lambda}(\!\cdot\!)$ are respectively increasing and decreasing continuous solutions of the ODE

(4.31) \begin{equation} \frac{\sigma^{2}}{2}f''(x)+\kappa(\theta-x)f'(x)=\lambda f(x) , \end{equation}

on $(0,\infty)$ . What is more, both the scale derivatives ${\textrm{d}}\phi_{\lambda}(\!\cdot\!)/{\textrm{d}} \mathbb m{s}_{X}$ and ${\textrm{d}}\psi_{\lambda}(\!\cdot\!)/{\textrm{d}} \mathbb m{s}_{X}$ are continuous at a, implying that

(4.32) \begin{equation} p\phi'_{\!\!\lambda}(a+\!)=(1-p)\phi'_{\!\!\lambda}(a-\!), \qquad p\psi'_{\!\!\lambda}(a+\!)=(1-p)\psi'_{\!\!\lambda}(a-\!), \end{equation}

and the following boundary conditions should be satisfied:

(4.33) \begin{equation} \phi'_{\!\!\lambda}(0+\!) = \frac{\lambda}{\beta}\phi_{\lambda}(0), \qquad \psi_{\lambda}(\infty) = 0, \qquad \frac{{\textrm{d}} \psi_{\lambda}}{{\textrm{d}} \mathbb m{s}_{X}}(\infty)=0. \end{equation}

Also, given the Wronskian defined by

(4.34) \begin{equation} w_{\lambda} = \psi_{\lambda}(x)\frac{{\textrm{d}} \phi_{\lambda}}{{\textrm{d}} \mathbb m{s}_{X}}(x) - \phi_{\lambda}(x)\frac{{\textrm{d}} \psi_{\lambda}}{{\textrm{d}} \mathbb m{s}_{X}}(x), \end{equation}

a constant independent of x, the Green function is of the form in (4.30). So it remains to make explicit $\phi_{\lambda}(\!\cdot\!)$ , $\psi_{\lambda}(\!\cdot\!)$ , and $w_{\lambda}$ .

First, consider the derivation of $\phi_{\lambda}(\!\cdot\!)$ . We conjecture that its expression takes the form (4.22) with the coefficients $c_{k}$ , $k=1,2,3$ , to be determined. From the continuity of $\phi_{\lambda}(\!\cdot\!)$ at a and the first equalities in both (4.32) and (4.33), and by virtue of the differential property

(4.35) \begin{equation} H'_{\nu}(x)=2\nu H_{\nu-1}(x), \end{equation}

we get three equations:

\begin{equation*} \left\{ \begin{aligned} & c_{1}H_{\nu}(\!-z(a))-c_{2}H_{\nu}(z(a))-c_{3}H_{\nu}(\!-z(a))=-H_{\nu}(z(a)),\\[5pt] & c_{1}(1-p)H_{\nu-1}(\!-z(a))+c_{2}pH_{\nu-1}(z(a))-c_{3}pH_{\nu-1}(\!-z(a))=(1-p)H_{\nu-1}(z(a)),\\[5pt] & c_{1}\left[2\beta H_{\nu-1}(\!-z(0))+\sigma\sqrt{\kappa}H_{\nu}(\!-z(0))\right]=2\beta H_{\nu-1}(z(0))-\sigma\sqrt{\kappa}H_{\nu}(z(0)). \end{aligned} \right. \end{equation*}

Solving for $c_{k}$ , $k=1,2,3$ , and applying the formula

(4.36) \begin{equation} H_{\nu-1}(\!-x)H_{\nu}(x)+H_{\nu-1}(x)H_{\nu}(\!-x)=\frac{2^{\nu}\sqrt{\pi}}{\Gamma(1-\nu)}{\textrm{e}}^{x^{2}} \end{equation}

yields the results (4.24)–(4.26).

Next we show that $\phi_{\lambda}(\!\cdot\!)$ is an increasing function on $[0,\infty)$ . Examine first the monotonicity on the interval [0, a]. Noting that $H_{\nu}(z(x))$ (resp. $H_{\nu}(\!-z(x))$ ) increases (resp. decreases) as x increases, we assert that $\phi_{\lambda}(\!\cdot\!)$ is increasing when $c_{1}\leq 0$ . When $c_{1}> 0$ , since

\begin{equation*} [H_{\nu}(z(x))+c_{1}H_{\nu}(\!-z(x))]'' = \frac{4\lambda(\lambda+\kappa)}{\sigma^{2}\kappa} [H_{\nu-2}(z(x))+c_{1}H_{\nu-2}(\!-z(x))]>0 \end{equation*}

by (4.35), $\phi_{\lambda}(\!\cdot\!)$ is strictly convex on [0, a]. This, combined with the observation that

\begin{equation*} \phi'_{\!\!\lambda}(0+) = \frac{\lambda}{\beta}\phi_{\lambda}(0) = \frac{2^{\nu+1}\lambda\sqrt{\pi}{\textrm{e}}^{z^{2}(0)}} {\Gamma(1-\nu)[2\beta H_{\nu-1}(\!-z(0))+\sigma\sqrt{\kappa}H_{\nu}(\!-z(0))]}>0, \end{equation*}

ensures the monotonic increase of $\phi_{\lambda}(\!\cdot\!)$ as well as its positiveness.

The increase of $\phi_{\lambda}(\!\cdot\!)$ on $(a,\infty)$ can be substantiated in a similar manner as long as we notice that $c_{2}>0$ . Indeed, simple calculations produce

\begin{align*} \phi_{\lambda}(a) & = H_{\nu}(z(a))+c_{1}H_{\nu}(\!-z(a))>0, \\[5pt] \phi'_{\!\!\lambda}(a-\!) & = \frac{2\lambda}{\sigma\sqrt{\kappa}\,}[H_{\nu-1}(z(a))-c_{1}H_{\nu-1}(\!-z(a))]>0, \end{align*}

which manifest that $K_{1},K_{2}>0$ , and therefore $c_{2}>0$ . So it is plain to see that $\phi_{\lambda}(\!\cdot\!)$ is increasing if $c_{3}\leq 0$ . If $c_{3}> 0$ ,

(4.37) \begin{equation} [c_{2}H_{\nu}(z(x))+c_{3}H_{\nu}(\!-z(x))]'' = \frac{4\lambda(\lambda+\kappa)}{\sigma^{2}\kappa}[c_{2}H_{\nu-2}(z(x))+c_{3}H_{\nu-2}(\!-z(x))]>0. \end{equation}

In addition, from the first equality in (4.32) and the fact that $\phi'_{\!\!\lambda}(a-\!)>0$ , we realize that $\phi'_{\!\!\lambda}(a+\!)>0$ , which further indicates that $\phi_{\lambda}(\!\cdot\!)$ is increasing on $(a,\infty)$ . The monotonicity property on $[0,\infty)$ follows immediately from the continuity of $\phi_{\lambda}(\!\cdot\!)$ at a.

Next, consider the derivation of $\psi_{\lambda}(\!\cdot\!)$ . We conjecture that its expression takes the form (4.23), with the coefficients $d_{k}$ , $k=1,2$ , to be determined. Using the second equality in (4.32) and the continuity of $\psi_{\lambda}(\!\cdot\!)$ at a, we obtain

\begin{equation*} \left\{\begin{aligned} & d_{1}H_{\nu}(z(a))+d_{2}H_{\nu}(\!-z(a))=H_{\nu}(\!-z(a)),\\[5pt] & d_{1}(1-p)H_{\nu-1}(z(a))-d_{2}(1-p)H_{\nu-1}(\!-z(a))=-pH_{\nu-1}(\!-z(a)), \end{aligned} \right. \end{equation*}

which produces (4.27) and (4.28) with the help of (4.36).

The decrease of $\psi_{\lambda}(\!\cdot\!)$ on $[0,\infty)$ can be demonstrated as follows. It is obvious that $\psi_{\lambda}(\!\cdot\!)$ is decreasing on $(a,\infty)$ . On [0, a], since $d_{2}>0$ , the desired conclusion holds if $d_{1}\leq 0$ (i.e. $p\geq1/2$ ). If $d_{1}> 0$ (i.e. $p<1/2$ ), the substitution of $c_{2}$ and $c_{3}$ by $d_{1}$ and $d_{2}$ respectively in (4.37) results in the strict convexity of $\psi_{\lambda}(\!\cdot\!)$ . This, together with the fact that

\begin{equation*} \psi'_{\!\!\lambda}(a-\!)=\frac{p}{1-p}\psi'_{\!\!\lambda}(a+\!)<0, \end{equation*}

indicates that $\psi_{\lambda}(\!\cdot\!)$ is decreasing on [0, a]. The monotonicity on $[0,\infty)$ follows immediately from the continuity of $\psi_{\lambda}(\!\cdot\!)$ at a.

Finally, for the Wronskian $w_{\lambda}$ , it follows from (4.34) that

\begin{align*} w_{\lambda} &= \psi_{\lambda}(a)\frac{\phi'_{\!\!\lambda}(a-\!)}{\mathfrak{s}_{X}(a-\!)} - \phi_{\lambda}(a)\frac{\psi'_{\!\!\lambda}(a+\!)}{\mathfrak{s}_{X}(a+\!)} \\[5pt] &= \frac{2\lambda}{\sigma\sqrt{\kappa}\,}{\textrm{e}}^{z^{2}(0)-z^{2}(a)} \bigg[\frac{K_{2}}{p}H_{\nu}(\!-z(a))+\frac{K_{1}}{1-p}H_{\nu-1}(\!-z(a))\bigg], \end{align*}

which is equal to the desired result in light of the explicit expression (4.25) for $c_{2}$ .

It is worth mentioning that a sticky reflecting diffusion has a probability mass at its sticky reflecting boundary. The following corollary insists on this point by offering the analytical form of

\begin{equation*} U_{\lambda}(x) = \int_{0}^{\infty}{\textrm{e}}^{-\lambda t}{\mathbb {P}}_{x}(X_{t}=0)\,{\textrm{d}}t,\end{equation*}

which is the Laplace transform of ${\mathbb {P}}_{x}(X_{t}=0)$ .

Corollary 4.1. For any $x\geq 0$ ,

\begin{equation*} U_{\lambda}(x) = \frac{H_{\nu}(z(0))+c_{1}H_{\nu}(\!-z(0))}{p\beta w_{\lambda}}\psi_{\lambda}(x), \end{equation*}

where $w_{\lambda}$ , $c_{1}$ , and $\psi_{\lambda}(\!\cdot\!)$ are defined in Theorem 4.1.

Proof. From the definition of the Green function and the results in Theorem 4.1, we see that

\begin{equation*} \int_{0}^{\infty}{\textrm{e}}^{-\lambda t}{\mathbb {P}}_{x}(X_{t}=0)\,{\textrm{d}}t = \frac{1}{p\beta}\int_{0}^{\infty}{\textrm{e}}^{-\lambda t}p(t;\;x,0)\,{\textrm{d}}t = \frac{1}{p\beta}\mathbb m{G}_{\lambda}(x,0) = \frac{\phi_{\lambda}(0)}{p\beta w_{\lambda}}\psi_{\lambda}(x), \end{equation*}

which completes the proof.

In the remainder of this section, we are concerned with the first hitting time at a random jump boundary defined by $O_{t}=b+\overline{O}\textbf{1}_{\{\varrho\leq t\}}$ , where $b>0$ is a constant, $\overline{O}$ is a positive random variable with distribution $F({\textrm{d}}y)$ , and $\varrho$ is exponentially distributed with parameter $q>0$ . Assume that $\overline{O}$ , $\varrho$ , and the process X are mutually independent, and set $\widetilde{\tau}_{b}=\inf\{t\geq0\colon X_{t}=O_{t}\}$ . We are interested in the following joint Laplace transform:

(4.38) \begin{equation} {\mathbb {E}}_{x}\big[{\textrm{e}}^{-r_{1}X_{\widetilde{\tau}_{b}}-r_{2}\widetilde{\tau}_{b}}\big], \qquad 0\leq x\leq b.\end{equation}

In applications, this kind of problem may translate into ones that emerge in storage system or insurance theory [29]. As pointed out in [6], the expectation (4.38), recalling the notation $\zeta_\lambda(x,y)$ defined in (4.21), is equal to

\begin{equation*} {\textrm{e}}^{-r_{1}b}\zeta_{q+r_{2}}(x,b) + {\textrm{e}}^{-r_{1}b}{\mathbb {E}}_{x}\big[{\textrm{e}}^{-r_{1}\overline{O}-r_{2}\varrho} \textbf{1}_{\{\widetilde{\tau}_{b}>\varrho\}}\zeta_{r_{2}}(X_{\varrho},b+\overline{O})\big],\end{equation*}

where the second term is determined by

\begin{equation*} V_{r_{1},r_{2}}^{(1)}(x)={\mathbb {E}}_{x}\big[{\textrm{e}}^{-r_{1}X_{\varrho}-r_{2}\varrho}\big], \qquad V_{r_{1},r_{2}}^{(2)}(x)={\mathbb {E}}_{x}\big[{\textrm{e}}^{-r_{1}X_{\varrho}-r_{2}\varrho}\textbf{1}_{\{\widetilde{\tau}_{b}<\varrho\}}\big],\end{equation*}

since they can jointly determine the distribution of $(X_{\varrho},\varrho,\textbf{1}_{\{\widetilde{\tau}_{b}>\varrho\}})$ by Laplace inversion. Thus, in what follows it suffices to compute $V_{r_{1},r_{2}}^{(1)}(\!\cdot\!)$ and $V_{r_{1},r_{2}}^{(2)}(\!\cdot\!)$ .

Lemma 4.1. ${\mathbb {E}}_{x}\big[\int_{0}^{\varrho}{\textrm{e}}^{-r_{2}s}\textbf{1}_{\{X_{s}=0\}}\,{\textrm{d}}s\big] = U_{q+r_{2}}(x)$ , where the explicit expression for $U_{\lambda}(\!\cdot\!)$ is provided in Corollary 4.1.

Proof. Applying the law of total expectation, we have

\begin{align*} {\mathbb {E}}_{x}\bigg[\int_{0}^{\varrho}{\textrm{e}}^{-r_{2}s}\textbf{1}_{\{X_{s}=0\}}\,{\textrm{d}}s\bigg] & = \int_{0}^{\infty}q{\textrm{e}}^{-qt}\,{\textrm{d}}t \int_{0}^{t}{\textrm{e}}^{-r_{2}s}{\mathbb {P}}_{x}(X_{s}=0)\,{\textrm{d}}s \\[5pt] & = \int_{0}^{\infty}{\textrm{e}}^{-r_{2}s}{\mathbb {P}}_{x}(X_{s}=0)\,{\textrm{d}}s \int_{s}^{\infty}q{\textrm{e}}^{-qt}\,{\textrm{d}}t \\[5pt] & = \int_{0}^{\infty}{\textrm{e}}^{-(q+r_{2})s}{\mathbb {P}}_{x}(X_{s}=0)\,{\textrm{d}}s = U_{q+r_{2}}(x). \end{align*}

Lemma 4.2. Recall the definition of the function $z(\!\cdot\!)$ provided in Theorem 4.1, and let $\bar{\nu}=-(q+r_{2})/{\kappa}$ and

\begin{equation*} \widetilde{K} = \frac{2(q+r_{2})}{\sigma\sqrt{\kappa}} [m_{1}(1-p)H_{\bar{\nu}-1}(z(a)) + (p-(1-p)m_{2})H_{\bar{\nu}-1}(\!-z(a))]. \end{equation*}

Then

(4.39) \begin{equation} \widetilde{U}_{r_{2}}(x) \triangleq {\mathbb {E}}_{x} \bigg[\int_{0}^{\varrho}{\textrm{e}}^{-r_{2}s}\,{\textrm{d}}\widehat{L}^{X}_{s}(a)\bigg] = \left\{\begin{array}{l@{\quad}l} \dfrac{1}{\widetilde{K}}[m_{1}H_{\bar{\nu}}(z(x))+m_{2}H_{\bar{\nu}}(\!-z(x))], & 0\leq x\leq a, \\[9pt] \dfrac{1}{\widetilde{K}}H_{\bar{\nu}}(\!-z(x)), & x>a, \end{array}\right. \end{equation}

where

(4.40) \begin{align} m_{1} & = \frac{H_{\bar{\nu}}(\!-z(a))[2\beta H_{\bar{\nu}-1}(\!-z(0))+\sigma\sqrt{\kappa}H_{\bar{\nu}}(\!-z(0))]} {\widehat{K}}, \end{align}

(4.41) \begin{align} m_{2} = \frac{H_{\bar{\nu}}(\!-z(a))[2\beta H_{\bar{\nu}-1}(z(0))-\sigma\sqrt{\kappa}H_{\bar{\nu}}(z(0))]} {\widehat{K}}, \end{align}

\begin{align*} \widehat{K} & = H_{\bar{\nu}}(\!-z(a))[2\beta H_{\bar{\nu}-1}(z(0))-\sigma\sqrt{\kappa}H_{\bar{\nu}}(z(0))] \nonumber \\[5pt] & \quad + H_{\bar{\nu}}(z(a))[2\beta H_{\bar{\nu}-1}(\!-z(0))+\sigma\sqrt{\kappa}H_{\bar{\nu}}(\!-z(0))]. \nonumber \end{align*}

Proof. Suppose that the process X starts at the level x. Then, for any nonzero function $f(\!\cdot\!)\in\mathcal{C}_\textrm{b}([0,\infty)) \cap \mathcal{C}_\textrm{b}^{1}([0,\infty)\backslash\{a\}) \cap \mathcal{C}^{2}([0,\infty)\backslash\{a\})$ with finite left and right derivatives of orders one and two at a, we obtain from the generalized Itô formula [14, p. 150] that

\begin{align*} f(X_{t}) & = f(x) + \int_{0}^{t}(\mathscr{A}f)(X_{s})\textbf{1}_{\{X_{s}> 0\}}\,{\textrm{d}}s + \sigma\int_{0}^{t}f'(X_{s})\textbf{1}_{\{X_{s}> 0\}}\,{\textrm{d}}B_{s} \\[5pt] & \quad + [pf'(a+\!)-(1-p)f'(a-\!)]\widehat{L}^{X}_{t}(a) + \frac{f'(0)}{2}L^{X}_{t}(0) \\[5pt] & = f(x) + \int_{0}^{t}(\mathscr{A}f)(X_{s})\,{\textrm{d}}s + \sigma\int_{0}^{t}f'(X_{s})\textbf{1}_{\{X_{s}> 0\}}\,{\textrm{d}}B_{s} \\[5pt] & \quad + [pf'(a+\!)-(1-p)f'(a-\!)]\widehat{L}^{X}_{t}(a) + \bigg[\frac{1}{2}f'(0)-\frac{1}{2\beta}(\mathscr{A}f)(0)\bigg]L^{X}_{t}(0), \end{align*}

where the second equality comes from (2.2), and $\mathscr{A}$ is the infinitesimal generator of X given in (3.14). Consequently, the Itô product rule gives

\begin{align*} {\textrm{e}}^{-(q+r_{2})t}f(X_{t}) & = f(x) + \int_{0}^{t}{\textrm{e}}^{-(q+r_{2})s}[(\mathscr{A}f)(X_{s})-(q+r_{2})f(X_{s})]\,{\textrm{d}}s \\[5pt] & \quad + \sigma\int_{0}^{t}{\textrm{e}}^{-(q+r_{2})s}f'(X_{s})\textbf{1}_{\{X_{s}> 0\}}\,{\textrm{d}}B_{s} + [pf'(a+\!)-(1-p)f'(a-\!)] \\[5pt] & \qquad \times \int_{0}^{t}{\textrm{e}}^{-(q+r_{2})s}\,{\textrm{d}}\widehat{L}^{X}_{s}(a) + \bigg[\frac{1}{2}f'(0)-\frac{1}{2\beta}(\mathscr{A}f)(0)\bigg] \int_{0}^{t}{\textrm{e}}^{-(q+r_{2})s}\,{\textrm{d}}L^{X}_{s}(0). \end{align*}

This indicates that, if $f(\!\cdot\!)$ satisfies

(4.42) \begin{equation} (\mathscr{A}f)(x)=(q+r_{2}) f(x) \end{equation}

on $[0,\infty)$ with the initial condition

(4.43) \begin{equation} f'(0)=\frac{1}{\beta}(\mathscr{A}f)(0)=\frac{q+r_{2}}{\beta}f(0), \end{equation}

we have

\begin{equation*} [(1-p)f'(a-\!)-pf'(a+\!)]{\mathbb {E}}_{x}\bigg[\int_{0}^{t}{\textrm{e}}^{-(q+r_{2})s}\,{\textrm{d}}\widehat{L}^{X}_{s}(a)\bigg] = f(x)-{\mathbb {E}}_{x}[{\textrm{e}}^{-(q+r_{2})t}f(X_{t})] \end{equation*}

as the process $\big\{\int_{0}^{t}{\textrm{e}}^{-(q+r_{2})s}f'(X_{s})\textbf{1}_{\{X_{s}> 0\}}\,{\textrm{d}}B_{s}, t\geq 0 \big\}$ is a martingale. Moreover, it is evident that $(1-p)f'(a-\!)-pf'(a+\!)\neq 0$ , since otherwise the use of the dominated convergence theorem would lead to $f\equiv 0$ on $[0,\infty)$ which violates the nonzero assumption on $f(\!\cdot\!)$ . So, we have

\begin{align*} \widetilde{U}_{r_{2}}(x) = {\mathbb {E}}_{x}\bigg[\int_{0}^{\infty}q{\textrm{e}}^{-qt}\,{\textrm{d}}t \int_{0}^{t}{\textrm{e}}^{-r_{2}s}\,{\textrm{d}}\widehat{L}^{X}_{s}(a)\bigg] & = {\mathbb {E}}_{x}\bigg[\int_{0}^{\infty}{\textrm{e}}^{-(q+r_{2})s}\,{\textrm{d}}\widehat{L}^{X}_{s}(a)\bigg] \\[5pt] & = \lim_{t\rightarrow\infty}{\mathbb {E}}_{x} \bigg[\int_{0}^{t}{\textrm{e}}^{-(q+r_{2})s}\,{\textrm{d}}\widehat{L}^{X}_{s}(a)\bigg] \\[5pt] & = \lim_{t\rightarrow\infty}\frac{f(x)-{\mathbb {E}}_{x}[{\textrm{e}}^{-(q+r_{2})t}f(X_{t})]}{(1-p)f'(a-\!)-pf'(a+\!)} \\[5pt] & = \frac{f(x)}{(1-p)f'(a-\!)-pf'(a+\!)} \end{align*}

by using the monotone convergence theorem for the third equality and the dominated convergence theorem for the last equality. Next, to work out the explicit expression for $f(\!\cdot\!)$ fulfilling (4.42) and (4.43), we construct $f(\!\cdot\!)$ via

(4.44)

where the expressions in (4.40) and (4.41) for the coefficients $m_{1}$ and $m_{2}$ can be obtained from the continuity of $f(\!\cdot\!)$ at a and the condition in (4.43). Furthermore, in light of

\begin{equation*} \lim_{x\rightarrow\infty}H_{\bar{\nu}}(\!-z(x)) = 0, \qquad \lim_{x\rightarrow\infty}H'_{\bar{\nu}}(\!-z(x)) = -\frac{2(q+r_{2})}{\sigma\sqrt{\kappa}} \lim_{x\rightarrow\infty}H_{\bar{\nu}-1}(\!-z(x)) = 0, \end{equation*}

we assert that $f(\!\cdot\!)$ of the form (4.44) obeys all the assumptions made on itself at the outset of the proof. The desired result (4.39) is thus achieved by noting that $\widetilde{K}=(1-p)f'(a-\!)-pf'(a+\!)$ .

We are now in a position to state the main result.

Theorem 4.2.

(4.45) \begin{align} V_{r_{1},r_{2}}^{(1)}(x) & = \frac{q}{\kappa}{\textrm{e}}^{\chi(r_{1})}\bigg[\int_{0}^{r_{1}}\frac{{\textrm{e}}^{-xy-\chi(y)}}{y}\,{\textrm{d}}y - U_{q+r_{2}}(x)\int_{0}^{r_{1}}\bigg(\frac{\sigma^{2}}{2}y+\beta-\kappa\theta\bigg) {\textrm{e}}^{-\chi(y)}\,{\textrm{d}}y \nonumber \\[5pt] & \qquad\qquad\quad - (2p-1)\widetilde{U}_{r_{2}}(x) \int_{0}^{r_{1}}{\textrm{e}}^{-ay-\chi(y)}\,{\textrm{d}}y\bigg], \end{align}

(4.46) \begin{align} V_{r_{1},r_{2}}^{(2)}(x) & = \zeta_{q+r_{2}}(x,b)V_{r_{1},r_{2}}^{(1)}(b), \end{align}

where the explicit expressions for $U_{\lambda}(\!\cdot\!)$ , $\widetilde{U}_{r_{2}}(\!\cdot\!)$ , and $\zeta_\lambda(\cdot,\cdot)$ are provided in Corollary 4.1, Lemma 4.2, and Theorem 4.1, respectively, and

\begin{equation*} \chi(x) = \frac{\sigma^{2}}{4\kappa}x^{2}-\theta x-\frac{q+r_{2}}{\kappa}\ln x. \end{equation*}

Proof. Conditional on $X_{0}=x$ , we have by the Itô formula [20, Theorem 3.3.6] and (2.2) that

\begin{align*} {\textrm{e}}^{-r_{1}X_{t}-r_{2}t} & = {\textrm{e}}^{-r_{1}x} + \int_{0}^{t}{\textrm{e}}^{-r_{1}X_{s}-r_{2}s} \bigg[\frac{1}{2}r_{1}^{2}\sigma^{2}-r_{1}\kappa(\theta-X_{s})-r_{2}\bigg]\,{\textrm{d}}s \\[5pt] & \quad - \bigg[\frac{1}{2}r_{1}^{2}\sigma^{2}+r_{1}(\beta-\kappa\theta)\bigg] \int_{0}^{t}{\textrm{e}}^{-r_{2}s}\textbf{1}_{\{X_{s}=0\}}\,{\textrm{d}}s - r_{1}\sigma\int_{0}^{t}{\textrm{e}}^{-r_{1}X_{s}-r_{2}s}\textbf{1}_{\{X_{s}>0\}}\,{\textrm{d}}B_s \\[5pt] & \quad - (2p-1)r_{1}{\textrm{e}}^{-r_{1}a}\int_{0}^{t}{\textrm{e}}^{-r_{2}s}\,{\textrm{d}}\widehat{L}^{X}_{s}(a), \end{align*}

which implies that

\begin{align*} V_{r_{1},r_{2}}^{(1)}(x) & = {\textrm{e}}^{-r_{1}x} + \bigg(\frac{1}{2}r_{1}^{2}\sigma^{2}-r_{1}\kappa\theta-r_{2}\bigg) {\mathbb {E}}_{x}\bigg[\int_{0}^{\varrho}{\textrm{e}}^{-r_{1}X_{s}-r_{2}s}\,{\textrm{d}}s\bigg] \\[5pt] & \quad + r_{1}\kappa{\mathbb {E}}_{x}\bigg[\int_{0}^{\varrho}X_{s}{\textrm{e}}^{-r_{1}X_{s}-r_{2}s}\,{\textrm{d}}s\bigg] - \bigg[\frac{1}{2}r_{1}^{2}\sigma^{2}+r_{1}(\beta-\kappa\theta)\bigg] {\mathbb {E}}_{x}\bigg[\int_{0}^{\varrho}{\textrm{e}}^{-r_{2}s}\textbf{1}_{\{X_{s}=0\}}\,{\textrm{d}}s\bigg] \\[5pt] & \quad - (2p-1)r_{1}{\textrm{e}}^{-r_{1}a} {\mathbb {E}}_{x}\bigg[\int_{0}^{\varrho}{\textrm{e}}^{-r_{2}s}\,{\textrm{d}}\widehat{L}^{X}_{s}(a)\bigg]. \end{align*}

Following a similar argument to the proof of Lemma 4.1 with suitable modifications, we obtain

\begin{equation*} {\mathbb {E}}_{x}\bigg[\int_{0}^{\varrho}{\textrm{e}}^{-r_{1}X_{s}-r_{2}s}\,{\textrm{d}}s\bigg] = \frac{1}{q}V_{r_{1},r_{2}}^{(1)}(x), \qquad {\mathbb {E}}_{x}\bigg[\int_{0}^{\varrho}X_{s}{\textrm{e}}^{-r_{1}X_{s}-r_{2}s}\,{\textrm{d}}s\bigg] = -\frac{1}{q}\frac{\partial V_{r_{1},r_{2}}^{(1)}(x)}{\partial r_{1}}. \end{equation*}

Thanks to Lemmas 4.1 and 4.2, we conclude that $V_{r_{1},r_{2}}^{(1)}(\!\cdot\!)$ satisfies

\begin{align*} \frac{\partial V_{r_{1},r_{2}}^{(1)}(x)}{\partial r_{1}} & = \frac{q}{r_{1}\kappa}{\textrm{e}}^{-r_{1}x} + \frac{\frac{1}{2}r_{1}^{2}\sigma^{2}-r_{1}\kappa\theta-r_{2}-q}{r_{1}\kappa}V_{r_{1},r_{2}}^{(1)}(x) - \frac{\frac{1}{2}\sigma^{2}r_{1}q+q(\beta-\kappa\theta)}{\kappa}U_{q+r_{2}}(x) \\[5pt] & \quad - \frac{q(2p-1)}{\kappa}{\textrm{e}}^{-r_{1}a}\widetilde{U}_{r_{2}}(x), \end{align*}

with $V_{0,r_{2}}^{(1)}(x)={q}/({q+r_{2}})$ . Solving this ODE gives the result in (4.45).

The derivation of the expression in (4.46) for $V_{r_{1},r_{2}}^{(2)}(\!\cdot\!)$ relies on the standard argument using the memoryless property of $\varrho$ and the strong Markov property of X (see the proof of [29, Lemma 5.2]) and is thus omitted.

5. Applications to bond pricing

In this section we present the application of the skew OU process with sticky reflection to interest rate modeling. It is assumed that the dynamics of the short-term interest rate follows the stochastic process $\{X_{t},t\geq 0\}$ which evolves according to the SDE system of (2.1) and (2.2) under the risk-neutral probability measure ${\mathbb {P}}$ . A special case is that when $p=0.5$ , our proposed model is reduced to the Vasicek interest rate model with sticky zero lower bound advocated by [28]. We concentrate on the valuation problem of the zero-coupon bond expired at time T whose initial price is given by

\begin{equation*} P(x,T) = {\mathbb {E}}_{x}\bigg[\exp\bigg\{{-}\int_{0}^{T}X_{t}\,{\textrm{d}}t\bigg\}\bigg].\end{equation*}

For ease of exposition, we consider in a more general setting the expectation

\begin{equation*} (\mathscr{P}_{T}f)(x) = {\mathbb {E}}_{x}\bigg[\exp\bigg\{{-}\int_{0}^{T}X_{t}\,{\textrm{d}}t\bigg\}f(X_{T})\bigg].\end{equation*}

Let $\varrho_{1}$ be an exponential random variable with parameter 1 independent of X, and put $\widehat{\tau} = \inf\big\{t\geq0\colon\int_{0}^{t}X_{s}\,{\textrm{d}}s \geq \varrho_{1}\big\}$ . Then the process $\{\widehat{X}_{t},t\geq 0\}$ defined by

\begin{equation*} \widehat{X}_{t} = \left\{\begin{array}{l@{\quad}l} X_{t}, & t< \widehat{\tau}, \\[5pt] \varpi, & t\geq \widehat{\tau} , \end{array}\right.\end{equation*}

where $\varpi$ denotes a cemetery state, is a diffusion process with lifetime $\widehat{\tau}$ having the infinitesimal generator $\mathscr{G}$ of the form

\begin{equation*} \mathscr{G}f = \frac{\sigma^{2}}{2}\frac{{\textrm{d}}^{2}f}{{\textrm{d}} x^{2}}+\kappa(\theta-x)\frac{{\textrm{d}} f}{{\textrm{d}} x}-xf.\end{equation*}

Thus, with the convention that $f(\varpi)=0$ , we have

\begin{equation*} (\mathscr{P}_{T}f)(x) = {\mathbb {E}}_{x}\big[f(\widehat{X}_{T})\textbf{1}_{\{t<\widehat{\tau}\}}\big] = {\mathbb {E}}_{x}[f(\widehat{X}_{T})].\end{equation*}

A feasible way to compute this quantity is to consider its Laplace transform in the time domain; we could derive the Green function $\widehat{\mathbb m{G}}_{\lambda}(\cdot,\cdot)$ of $\widehat{X}$ and get

\begin{equation*} \int_{0}^{\infty}{\textrm{e}}^{-\lambda T}{\mathbb {E}}_{x}[f(\widehat{X}_{T})]\,{\textrm{d}}T = \int_{0}^{\infty}f(y)\widehat{\mathbb m{G}}_{\lambda}(x,y)\,\mathbb m{m}_X({\textrm{d}}y)\end{equation*}

accordingly. Nevertheless, we expect to find the explicit formula for $(\mathscr{P}_{T}f)(\!\cdot\!)$ here, especially for the case when $f\equiv 1$ . To achieve this, the spectral expansion approach (see [23] for a lucid account of the theory) will be adopted. In fact, let $\mathcal {L}([0,\infty),\mathbb m{m}_X)$ denote the Hilbert space of functions on $[0,\infty)$ square-integrable with respect to $\mathbb m{m}_X$ equipped with the inner product $(f_{1},f_{2}) = \int_{0}^{\infty}f_{1}(x)f_{2}(x)\,\mathbb m{m}_X({\textrm{d}}x)$ . Then the Feller semigroup $\{\mathscr{P}_{t},t\geq0\}$ restricted to $\mathcal{C}_\textrm{b}([0,\infty)) \cap \mathcal{L}([0,\infty),\mathbb m{m}_X)$ can be extended to a strongly continuous semigroup of self-adjoint contractions in $\mathcal{L}([0,\infty),\mathbb m{m}_X)$ with the infinitesimal generator $\mathscr{G}$ . Since both 0 and $\infty$ are nonoscillatory boundaries [28], it follows from the spectral theorem that the spectrum of $\mathscr{G}$ is simple and purely discrete, and thus, for any $f(\!\cdot\!)\in\mathcal{L}([0,\infty),\mathbb m{m}_X)$ ,

(5.47) \begin{equation} (\mathscr{P}_{T}f)(x) = \sum_{n=0}^{\infty}l_{n}{\textrm{e}}^{-\lambda_{n}T}\varphi_{n}(x),\end{equation}

where $\{\lambda_{n}\}_{n=0}^{\infty}$ , $0<\lambda_{0}<\lambda_{1}<\cdots<\lambda_{n}\uparrow\infty$ , are the eigenvalues of $-\mathscr{G}$ , $\{\varphi_{n}(\!\cdot\!)\}_{n=1}^{\infty}$ are the corresponding normalized eigenfunctions, and $l_{n}=(f,\varphi_{n})$ .

Now introduce

\begin{equation*} \alpha = \sigma\sqrt{\frac{2}{\kappa^{3}}}, \qquad \varsigma_{\lambda} = \frac{\sigma^{2}}{2\kappa^{3}} + \frac{\lambda-\theta}{\kappa}, \qquad \bar{z}(x) = \frac{\sqrt{2\kappa}}{\sigma}(\theta-x),\end{equation*}

and denote by $D_{\varsigma}(\!\cdot\!)$ , $E^{(0)}_{\varsigma}(\!\cdot\!)$ , and $E^{(1)}_{\varsigma}(\!\cdot\!)$ the parabolic cylinder functions [10] which are linked to the Hermite function via

\begin{align*} D_{\varsigma}(x) & = 2^{-{\varsigma}/{2}}{\textrm{e}}^{-{x^{2}}/{4}}H_{\varsigma}\bigg(\frac{x}{\sqrt{2}\,}\bigg), \\[5pt] & E^{(0)}_{\varsigma}(x) = \frac{\Gamma(({1-\varsigma})/{2})}{2^{({1}/{2})+\varsigma}\sqrt{\pi}\,}{\textrm{e}}^{-{x^{2}}/{4}} \bigg[H_{\varsigma}\bigg(\frac{x}{\sqrt{2}\,}\bigg) + H_{\varsigma}\bigg({-}\frac{x}{\sqrt{2}}\bigg)\bigg], \\[5pt] & E^{(1)}_{\varsigma}(x) = -\frac{\Gamma(\!-{\varsigma}/{2})}{2^{({1}/{2})+\varsigma}\sqrt{\pi}\,}{\textrm{e}}^{-{x^{2}}/{4}} \bigg[H_{\varsigma}\bigg(\frac{x}{\sqrt{2}\,}\bigg) - H_{\varsigma}\bigg({-}\frac{x}{\sqrt{2}\,}\bigg)\bigg].\end{align*}

Specializing the result (5.47) to the case when $f\equiv 1$ , we obtain the bond price given below.

Theorem 5.1. The zero-coupon bond price P(x,T) under the skew OU with sticky reflection model takes the series form (5.47) where the constants $\{\lambda_{n}\}_{n=1}^{\infty}$ , $\{l_{n}\}_{n=1}^{\infty}$ and the functions $\{\varphi_{n}(\!\cdot\!)\}_{n=1}^{\infty}$ are provided as follows. Write

(5.48) \begin{align} \bar{c}_{1} & = 2^{-({\varsigma_{\lambda}}/{2})-2}{\textrm{e}}^{-{\bar{z}^{2}(0)}/{4}}\bigg[2E^{(0)}_{\varsigma_{\lambda}-1}(\bar{z}(0)-\alpha) + \bigg(\frac{\sqrt{2}\alpha}{2} - \frac{\lambda\sigma}{\beta\sqrt{\kappa}\,}\bigg) E^{(1)}_{\varsigma_{\lambda}}(\bar{z}(0)-\alpha)\bigg], \quad\qquad \end{align}

(5.49) \begin{align} \bar{c}_{2} & = 2^{-({\varsigma_{\lambda}}/{2})-2}{\textrm{e}}^{-{\bar{z}^{2}(0)}/{4}} \bigg[\varsigma_{\lambda}E^{(1)}_{\varsigma_{\lambda}-1}(\bar{z}(0)-\alpha) + \bigg(\frac{\lambda\sigma}{\beta\sqrt{\kappa}\,} - \frac{\sqrt{2}\alpha}{2}\bigg) E^{(0)}_{\varsigma_{\lambda}}(\bar{z}(0)-\alpha)\bigg]; \quad\qquad \end{align}

then $\{\lambda_{n}\}_{n=0}^{\infty}$ are the discrete positive zeros of $\overline{w}_{\lambda}$ defined by

(5.50) \begin{align} \overline{w}_{\lambda} & = \frac{\sqrt{2\kappa}}{p(1-p)\sigma}{\textrm{e}}^{{\bar{z}^{2}(0)}/{2}} \bigg\{\bar{c}_{1}\bigg[(1-p)D_{\varsigma_{\lambda}}(\alpha-\bar{z}(a)) \bigg(\frac{\alpha}{2}E^{(0)}_{\varsigma_{\lambda}}(\bar{z}(a)-\alpha) - \frac{\varsigma_{\lambda}}{\sqrt{2}}E^{(1)}_{\varsigma_{\lambda}-1}(\bar{z}(a)-\alpha)\bigg) \nonumber \\[5pt] & \qquad\qquad - pE^{(0)}_{\varsigma_{\lambda}}(\bar{z}(a)-\alpha) \bigg(\frac{\alpha}{2}D_{\varsigma_{\lambda}}(\alpha-\bar{z}(a)) - \varsigma_{\lambda}D_{\varsigma_{\lambda}-1}(\alpha-\bar{z}(a))\bigg)\bigg] \nonumber \\[5pt] & \qquad\qquad + \bar{c}_{2}\bigg[(1-p)D_{\varsigma_{\lambda}}(\alpha-\bar{z}(a)) \bigg(\frac{\alpha}{2}E^{(1)}_{\varsigma_{\lambda}}(\bar{z}(a)-\alpha) + \sqrt{2}E^{(0)}_{\varsigma_{\lambda}-1}(\bar{z}(a)-\alpha)\bigg) \nonumber \\[5pt] & \qquad\qquad - pE^{(1)}_{\varsigma_{\lambda}}(\bar{z}(a)-\alpha) \bigg(\frac{\alpha}{2}D_{\varsigma_{\lambda}}(\alpha-\bar{z}(a)) - \varsigma_{\lambda}D_{\varsigma_{\lambda}-1}(\alpha-\bar{z}(a))\bigg)\bigg]\bigg\}, \end{align}

i.e. $\overline{w}_{\lambda_{n}}=0$ . Moreover,

(5.51) \begin{equation} \varphi_{n}(x) = \left\{\begin{array}{l@{\quad}l} \sqrt{\dfrac{\rho_{\lambda_{n}}^{(2)}(a)}{\rho_{\lambda_{n}}^{(1)}(a)\overline{w}'_{\lambda_{n}}}} \rho_{\lambda_{n}}^{(1)}(x), & 0\leq x\leq a, \\[5pt] \textrm{sgn}\big(\rho_{\lambda_{n}}^{(1)}(a)\rho_{\lambda_{n}}^{(2)}(a)\big) \sqrt{\dfrac{\rho_{\lambda_{n}}^{(1)}(a)}{\rho_{\lambda_{n}}^{(2)}(a)\overline{w}'_{\lambda_{n}}}} \rho_{\lambda_{n}}^{(2)}(x), & x>a, \end{array} \right. \end{equation}

where $\overline{w}'_{\lambda}$ denotes the first-order derivative of $\overline{w}_{\lambda}$ with respect to $\lambda$ , and

\begin{align*} \rho_{\lambda}^{(1)}(x) & = {\textrm{e}}^{{\bar{z}^{2}(x)}/{4}}\big[\bar{c}_{1}E^{(0)}_{\varsigma_{\lambda}}(\bar{z}(x)-\alpha) + \bar{c}_{2}E^{(1)}_{\varsigma_{\lambda}}(\bar{z}(x)-\alpha)\big], \\[5pt] \rho_{\lambda}^{(2)}(x) & = {\textrm{e}}^{{\bar{z}^{2}(x)}/{4}}D_{\varsigma_{\lambda}}(\alpha-\bar{z}(x)). \end{align*}

The coefficient $l_{n}$ is given by

(5.52) \begin{align} l_{n} & = \sqrt{\frac{\rho_{\lambda_{n}}^{(2)}(a)}{\rho_{\lambda_{n}}^{(1)}(a) \overline{w}'_{\lambda_{n}}}} \bigg[\int_{0}^{a}\rho_{\lambda_{n}}^{(1)}(x)\mathfrak{m}_X(x)\,{\textrm{d}}x + \frac{\rho_{\lambda_{n}}^{(1)}(0)}{p\beta}\bigg] \nonumber \\[5pt] & \quad + \textrm{sgn}\big(\rho_{\lambda_{n}}^{(1)}(a)\rho_{\lambda_{n}}^{(2)}(a)\big) \sqrt{\frac{\rho_{\lambda_{n}}^{(1)}(a)}{\rho_{\lambda_{n}}^{(2)}(a) \overline{w}'_{\lambda_{n}}}} \bigg[\int_{a}^{\infty}\rho_{\lambda_{n}}^{(2)}(x)\mathfrak{m}_X(x)\,{\textrm{d}}x\bigg], \end{align}

with $\mathfrak{m}_X(\!\cdot\!)$ defined in (3.16).

Proof. Consider the Sturm–Liouville problem

(5.53) \begin{equation} -\frac{\sigma^{2}}{2}f''(x)-\kappa(\theta-x)f'(x)+xf(x)=\lambda f(x) \end{equation}

on $[0,\infty)$ with the boundary condition $f'(0)=-({\lambda}/{\beta})f(0)$ . With the transformation

\begin{equation*} f(x) = {\textrm{e}}^{{\bar{z}^{2}(x)}/{4}}h(\bar{z}(x)), \end{equation*}

the ODE in (5.53) is converted into the Whittaker differential equation

\begin{equation*} h''_{\!\!\!\bar{z}\bar{z}} + \bigg[\frac{1}{2}+\varsigma_{\lambda}-\frac{(\alpha-\bar{z})^{2}}{4}\bigg]h = 0. \end{equation*}

As claimed by [10], $E^{(0)}_{\varsigma_{\lambda}}(\bar{z}-\alpha)$ and $E^{(1)}_{\varsigma_{\lambda}}(\bar{z}-\alpha)$ form a pair of linearly independent solutions to the above equation, and so do $D_{\varsigma_{\lambda}}(\alpha-\bar{z})$ and $D_{-\varsigma_{\lambda}-1}(\textrm{i}(\alpha-\bar{z}))$ . Now let $\Upsilon_{\lambda}(\!\cdot\!)$ be the unique (up to a constant multiple) continuous solution with continuous scale derivative of the ODE (5.53) on $[0,\infty)$ which is $\mathbb m{m}_X$ -square integrable near 0 and meets

(5.54) \begin{equation} \Upsilon'_{\!\!\lambda}(0) = -\frac{\lambda}{\beta}{\Upsilon_{\lambda}}(0). \end{equation}

Symmetrically, let $\Psi_{\lambda}(\!\cdot\!)$ be the unique (up to a constant multiple) continuous solution with continuous scale derivative of the ODE (5.53) on $[0,\infty)$ which is $\mathbb m{m}_X$ -square integrable near $\infty$ and meets

\begin{equation*} \lim_{x\rightarrow\infty}\frac{{\textrm{d}}\Psi_{\lambda}}{{\textrm{d}}\mathbb m{s}_{X}}(x) = 0. \end{equation*}

Then we have

\begin{align*} \Upsilon_{\lambda}(x) & = \left\{\begin{array}{l@{\quad}l} \bar{c}_{1}{\textrm{e}}^{{\bar{z}^{2}(x)}/{4}}E^{(0)}_{\varsigma_{\lambda}}(\bar{z}(x)-\alpha) + \bar{c}_{2}{\textrm{e}}^{{\bar{z}^{2}(x)}/{4}}E^{(1)}_{\varsigma_{\lambda}}(\bar{z}(x)-\alpha), & \qquad\ 0 \leq x \leq a, \\[5pt] \bar{c}_{3}{\textrm{e}}^{{\bar{z}^{2}(x)}/{4}}E^{(0)}_{\varsigma_{\lambda}}(\bar{z}(x)-\alpha) + \bar{c}_{4}{\textrm{e}}^{{\bar{z}^{2}(x)}/{4}}E^{(1)}_{\varsigma_{\lambda}}(\bar{z}(x)-\alpha), & \qquad\ x > a; \end{array} \right. \end{align*}

where the coefficients $\bar{c}_{k}$ , $k=1,2,3,4$ , could be derived by using the boundary condition (5.54) and the normalized condition $\Upsilon_{\lambda}(0)=1$ (to reduce one degree of freedom), as well as the continuity at a of both $\Upsilon_{\lambda}(\!\cdot\!)$ and ${{\textrm{d}}\Upsilon_{\lambda}}(\!\cdot\!)/{{\textrm{d}}\mathbb m{s}_{X}}$ , and where the coefficients $\bar{d}_{k}$ , $k=1,2$ , could be derived by using the continuity at a of both $\Psi_{\lambda}(\!\cdot\!)$ and ${{\textrm{d}}\Psi_{\lambda}}(\!\cdot\!)/{{\textrm{d}}\mathbb m{s}_{X}}$ . In particular, by virtue of the relations

\begin{equation*} {E^{(0)}_{\varsigma}}(x) = -\frac{x}{2}{E^{(0)}_{\varsigma}}(x) - \frac{\varsigma}{\sqrt{2}}{E^{(1)}_{\varsigma-1}}(x), \qquad {E^{(1)}_{\varsigma}}'(x) = -\frac{x}{2}{E^{(1)}_{\varsigma}}(x) + {\sqrt{2}}{E^{(0)}_{\varsigma-1}}(x), \end{equation*}

we obtain the expressions (5.48) and (5.49) for $\bar{c}_{1}$ and $\bar{c}_{2}$ .

The Wronskian defined by

\begin{equation*} \overline{w}_{\lambda} = \Upsilon_{\lambda}(x)\frac{{\textrm{d}}\Psi_{\lambda}}{{\textrm{d}}\mathbb m{s}_{X}}(x) - \Psi_{\lambda}(x)\frac{{\textrm{d}}\Upsilon_{\lambda}}{{\textrm{d}}\mathbb m{s}_{X}}(x) \end{equation*}

is independent of x and could be further expressed as

\begin{equation*} \overline{w}_{\lambda} = \Upsilon_{\lambda}(a)\frac{\Psi'_{\lambda}(a+\!)}{\mathfrak{s}_{X}(a+\!)} - \Psi_{\lambda}(a)\frac{\Upsilon'_{\lambda}(a-\!)}{\mathfrak{s}_{X}(a-\!)} = \rho_{\lambda}^{(1)}(a)\frac{{\rho_{\lambda}^{(2)}}'(a+\!)}{\mathfrak{s}_{X}(a+\!)} - \rho_{\lambda}^{(2)}(a)\frac{{\rho_{\lambda}^{(1)}}'(a-\!)}{\mathfrak{s}_{X}(a-\!)}, \end{equation*}

which, coupled with the property $D_{\varsigma}'(x) = \varsigma D_{\varsigma-1}(x) - ({x}/{2})D_{\varsigma}(x)$ , yields the formula in (5.50). According to the conclusions reached by [23, p. 351–352], the positive discrete roots of $\overline{w}_{\lambda}=0$ are exactly the eigenvalues $\{\lambda_{n}\}_{n=0}^{\infty}$ . Also, for $\lambda=\lambda_{n}$ , $\Upsilon_{\lambda}(\!\cdot\!)$ and $\Psi_{\lambda}(\!\cdot\!)$ become linearly dependent, and thus

\begin{equation*} \Upsilon_{\lambda_{n}}(x) = \begin{cases} \rho_{\lambda_{n}}^{(1)}(x), & 0\leq x\leq a, \\[5pt] \frac{\rho_{\lambda_{n}}^{(2)}(x)}{M_{n}}, & x>a, \end{cases} \qquad \text{where } M_{n} = \frac{\Psi_{\lambda_{n}}(a)}{\Upsilon_{\lambda_{n}}(a)} = \frac{\rho_{\lambda_{n}}^{(2)}(a)}{\rho_{\lambda_{n}}^{(1)}(a)}. \end{equation*}

The corresponding normalized eigenfunction $\varphi_{n}(\!\cdot\!)$ is given by

\begin{align*} \varphi_{n}(x) = \sqrt{\frac{M_{n}}{\overline{w}'_{\lambda_{n}}}}\Upsilon_{\lambda_{n}}(x) = \left\{\begin{array}{l@{\quad}l} \sqrt{\dfrac{M_{n}}{\overline{w}'_{\lambda_{n}}}}\rho_{\lambda_{n}}^{(1)}(x), & 0\leq x\leq a, \\[5pt] \textrm{sgn}\big(\rho_{\lambda_{n}}^{(1)}(a)\rho_{\lambda_{n}}^{(2)}(a)\big) \sqrt{\dfrac{1}{M_{n}\overline{w}'_{\lambda_{n}}}}\rho_{\lambda_{n}}^{(2)}(x), & x>a. \end{array}\right. \end{align*}

Finally, the integral expression in (5.52) for the expansion coefficient $l_{n}$ results from (5.51) and the fact that $l_{n} = \int_{[0,\infty)}\varphi_{n}(x)\,\mathbb m{m}_X({\textrm{d}}x)$ .

Next, we give numerical examples concerning the yield on the zero-coupon bond which refers to the return an investor realizes on a bond and is formally defined as

\begin{equation*} R(x,T) = -\frac{\ln P(x,T)}{T}.\end{equation*}

When evaluating bond prices, the built-in function FindRoot in Mathematica software is used to search for the positive zeros of $\overline{w}_{\lambda}$ and NIntegrate to compute the integrals appearing in (5.52) numerically. Also, all the special functions involved can be computed by employing their relationships with the Hermite function available as HermiteH in Mathematica.

Figure 1 depicts the yield curves as functions of time to maturity, where we examine the effects on bond yields of the skewness parameter p and the stickiness parameter $\beta$ . It can be seen from the left panel that the smaller the value of p, the lower the bond yield and the slope of the curve, because it is more likely for the interest rate to be constrained within the interval [0, a] for smaller p. On a related matter, the parameter a actually acts as a resistance level in the terminology of [4] if $p=0.2$ , as opposed to as a support level if $p=0.8$ . The right panel (under the setting $p=0.2$ ) reveals that the bond yield and the slope of the curve increase with larger values of $\beta$ , a phenomenon that was also found by [28] in the framework of the sticky reflecting OU interest rate model (i.e. the case when $p=0.5$ in our model). Moreover, both the values of the yield and the slope are bounded above by those of the instantaneous reflection case, which corresponds to the infinite $\beta$ value.

Figure 1: Yields on zero-coupon bonds against time to maturity with different skewness and stickiness parameters. In the left panel, the solid, dashed, and dotted lines correspond to the cases when $p=0.2$ , $0.5$ , and $0.8$ , respectively. In the right panel, the solid, dashed, dot-dashed, and dotted lines correspond to the cases when $\beta=0.01$ , $0.1$ , 1, and $\infty$ , respectively. The parameter values in the base case are: $X_{0}=0.02$ , $\kappa=0.1$ , $\theta=0.04$ , $\sigma=0.15$ , $a=0.01$ , $p=0.2$ , and $\beta=0.1$ .

To close this section, we remark that, aside from interest rate modeling, skew diffusions with sticky reflection may also find potential applications in particle dispersion problems and queuing systems. For example, in view of the close connections between skew diffusions (with reflecting boundaries) and random particle motions in layered media (with reflecting boundaries) [2, 33, 34], it is natural to consider which skew diffusion with sticky reflection can describe the particle motion in a layered medium with sticky reflecting boundary. On the other hand, noting that we can transform a skew diffusion with sticky reflection into a threshold diffusion (i.e. a diffusion with discontinuous coefficients) with sticky reflection, it is possible to construct a queuing system with the exceptional service policy whose heavy traffic limit is the latter process [9, 26, 32]. The above problems will be on the list of our future research topics.

6. Conclusion

Motivated by the applications to interest rate modeling, we considered the skew OU process with sticky reflecting boundary. Since this process has not been studied before, we started by establishing the existence and uniqueness of weak solutions to the SDE system satisfied by the process. Then, various important probabilistic properties, including the scale function, the speed measure, the stationary distribution, the Green function, and the distributions of first hitting times to constant and random jump boundaries were obtained in terms of analytical solutions. Finally, based on the skew OU sticky interest rate model, we evaluated zero-coupon bonds by using the spectral expansion approach, and numerical results showed the positive effects of skewness and stickiness parameters on bond yields.

The techniques employed in the current paper can be easily extended to deal with the case of multiple skew points, or with CIR and constant elasticity of variance (CEV) skew processes with sticky reflection. However, for processes with general drift and diffusion coefficients, the analytical tractability would be lost. To handle that case, we can first transform the skew diffusion with sticky reflection into the one without skewness, and then resort to the computational method developed by [25], where the continuous-time Markov chain (CTMC) is used to approximate a sticky diffusion. Under the CTMC model, the value functions (including the bond price as a special case) and the first hitting time probabilities can be obtained by computing matrix exponentials. The simulation problem for sticky diffusions can also be solved by means of the CTMC approximation. Similarly, we could further investigate the properties of multidimensional skew diffusions with sticky reflecting boundaries via the multidimensional CTMC method proposed by [26].