1. Introduction
Since the seminal works of [Reference Itô and McKean18] and [Reference Walsh41], 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, [Reference Harrison and Shepp17] 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, [Reference Appuhamillage2, Reference Atar and Budhiraja3, Reference Borodin and Salminen8, Reference Burdzy and Chen11, Reference Étoré and Martinez15]. There are also a number of papers devoted to the investigation of more general skew diffusions. For example, [Reference Blei5] introduced the skew Bessel process by determining the jumps of the local time in a certain way. By means of moving domains, [Reference Trutnau39] 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 [Reference Trutnau40]. The distributions of the skew Ornstein–Uhlenbeck (OU) process were probed in [Reference Wang, Song and Wang42], 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 [Reference Feller16]. 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 [Reference Amir1, Reference Engelbert and Peskir13, Reference Peskir30, Reference Rácz and Shkolnikov32, Reference Song and Wang38] 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 [Reference Longstaff24], use of the sticky reflecting boundary is consistent with empirical observations on treasury-bill yields when near zero. In view of that finding, [Reference Nie27] checked the empirical performance of a two-dimensional sticky interest rate model, and [Reference Nie and Linetsky28] 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, [Reference Meier, Li and Zhang25, Reference Meier, Li and Zhang26] obtained bond prices under single- and multi-factor sticky interest rate models, respectively. On the other hand, [Reference Decamps, Goovaerts and Schoutens12] 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 [Reference Bai, Wang, Zhang and Zhuo4] 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 [Reference Nie and Linetsky28, Reference Peskir30] 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 [Reference Bo, Wang and Yang6, Reference Perry, Stadje and Zacks29], 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 [Reference Linetsky23] to solve derivative-pricing problems. Our interest rate model contains the one used by [Reference Nie and Linetsky28] 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 [Reference Protter31, 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 [Reference Wang, Song and Wang42, Reference Wang, Xu and Song43] 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 [Reference Yamada44] 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 [Reference Karatzas and Shreve20] 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:
\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
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 [Reference Engelbert and Schmidt14, 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 [Reference Blei5, Lemma A.5], and the last equality holds almost surely because the occupation times formula [Reference Revuz and Yor35, 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)]
-
(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$
. -
(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. -
(iii)
$|\widetilde{{\sigma}}(x)|\geq \min\{p,1-p\}\sigma$
for any
$x\geq (1-p)a$
.
Hence, it follows from [Reference Semrau37, Theorem 4.2] and [Reference Rong36, 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
\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 [Reference Protter31, 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
\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
\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 [Reference Karatzas and Shreve20, 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,
\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 [Reference Karatzas and Shreve20, 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
\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$
,
\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 [Reference Engelbert and Schmidt14, p. 150] produces
\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 [Reference Blei5, Lemma A.5],
\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
\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 [Reference Protter31, 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 [Reference Karatzas and Shreve20, 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 [Reference Borodin and Salminen7] (see also [Reference Itô and McKean18]), 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
\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 [Reference Kohatsu-Higa, Taguchi and Zhong21, Reference Nie and Linetsky28], 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
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 [Reference Revuz and Yor35, 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 [Reference Revuz and Yor35, p. 301, Proposition 3.2], for
$0\leq b_{1}<y<b_{2}$
,
On the other hand, it follows from the relation
$\tau^{Y}_{b}=\tau^{Z}_{G(b)}$
that
As a consequence, the combination of (3.17) and (3.18) leads to
\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 [Reference Revuz and Yor35, 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
\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 [Reference Borodin and Salminen7, 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 [Reference Kallenberg19, 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 [Reference Borodin and Salminen7, 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
\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$
[Reference Lebedev22, p. 285], and
$\Gamma(\!\cdot\!)$
the gamma function [Reference Lebedev22, 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
\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}
\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
\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}
\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}
\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}
\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}
\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$
,
\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}
Proof. The general theory of linear diffusions [Reference Borodin and Salminen7, 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
\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
\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:
\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
\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
\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
\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$
,
\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:
\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 [Reference Perry, Stadje and Zacks29]. As pointed out in [Reference Bo, Wang and Yang6], 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
\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
\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}
\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 [Reference Engelbert and Schmidt14, 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
\begin{equation} (\mathscr{A}f)(x)=(q+r_{2}) f(x) \end{equation}
on
$[0,\infty)$
with the initial condition
\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
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.
\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}
\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 [Reference Karatzas and Shreve20, 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 [Reference Perry, Stadje and Zacks29, 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 [Reference Nie and Linetsky28]. 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 [Reference Linetsky23] 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 [Reference Nie and Linetsky28], 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)$
,
\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 [Reference Buchholz10] 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
\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}
\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
\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,
\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
\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
\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 [Reference Buchholz10],
$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
\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 [Reference Linetsky23, 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 [Reference Bai, Wang, Zhang and Zhuo4] 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 [Reference Nie and Linetsky28] 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) [Reference Appuhamillage2, Reference Ramirez33, Reference Ramirez34], 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 [Reference Browne, Whitt and Dshalalow9, Reference Meier, Li and Zhang26, Reference Rácz and Shkolnikov32]. 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 [Reference Meier, Li and Zhang25], 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 [Reference Meier, Li and Zhang26].
Acknowledgements
We wish to thank the anonymous reviewer of this paper for constructive comments.
Funding information
The first author was supported by the National Natural Science Foundation of China (Grant No. 11801407). The second author was supported by the National Natural Science Foundation of China (Grant No. 11701085) and by ‘the Fundamental Research Funds for the Central Universities’ in UIBE (22QN07).
Competing interests
There were no competing interests to declare which arose during the preparation or publication process of this article.
