2.1. Existence and uniqueness of weak solutions

We make the following assumption.

The well-posedness of reflected SDEs has been well studied, and conditions can be found in, e.g., [Reference Lions and Sznitman37], [Reference Skorokhod47], [Reference SŁomínski48], [Reference Tanaka50], and [Reference Watanabe52]. Let $X^1$ denote the diffusion instantaneously reflected at l. We apply the time change

\begin{equation*}T_t=\bigg(t+\frac{1}{2\rho}L^l_t\big(X^1\big)\bigg)^{-1}\end{equation*}

to $X^1$ , which slows it down whenever it hits l. This introduces stickiness, which can be thought of as slow reflection. This construction enables us to show that there exists a unique weak solution to the sticky SDE system as long as this is so for the reflecting case.

The proof of this theorem can be found in the online supplement; it is similar to the proof of [Reference Nie42, Theorem 4.1]. We follow the standard approach to study the diffusion process with reflecting boundary condition and apply a time change to slow down the movement into the interior, which creates stickiness. However, the fact that the time change is continuous and strictly increasing is shown differently, as the arguments in [Reference Nie42] are only valid for the OU process and cannot be applied to general diffusions.

The following theorem shows that it is possible to add the local time given in (2.2) as an additional term in the drift.

Theorem 2. If $\rho\in(0,\infty)$ , then (2.1) and (2.2) are equivalent to the following SDE:

(2.3) \begin{equation}dX_t=\mu\left(X_t\right)I\left(X_t>l\right)dt+\rho I\left(X_t=l\right)dt+\sigma\left(X_t\right)I\left(X_t>l\right)dB_t .\end{equation}

Furthermore, X is a strong Markov process.

The proof of this theorem is omitted, as it is similar to that of Theorem 4.1 and Corollary 4.3 in [Reference Nie42]. Theorem 1 and Theorem 2 together show that X is a diffusion process. Its stickiness at l is measured by $\rho$ ; the smaller the value of $\rho$ , the more X sticks to l.

2.2. Feynman–Kac semigroup and infinitesimal generator

The Feynman–Kac operator $\mathcal{P}_t$ associated with the diffusion X is defined by

(2.4) \begin{equation}\mathcal{P}_tf(x):=\mathbb{E}_x\Big(\exp\left(-\int_0^tk\left(X_u\right)du\right)f\left(X_t\right)\Big), \qquad x\in\mathbb{S},\end{equation}

where the function k is assumed to be nonnegative. We can interpret k as the discount factor and f as the payoff function. Let $u(t,x)=\mathcal{P}_tf(x)$ ; this is called the value function. Financially, u(t,x) is the price for receiving payoff f at t given the current state x. The semigroup $${({{\mathcal P}_t})_{t \ge 0}}$$ is a strongly continuous contraction semigroup on $B_b(\mathbb{S})$ , the space of bounded Borel-measurable functions, endowed with the maximum norm.

Assume that $\mu$ , $\sigma$ , and k are continuous over $\mathbb{S}$ . Using the arguments in [Reference Nie42], one can easily show that the infinitesimal generator $\mathcal{G}$ of $${({{\mathcal P}_t})_{t \ge 0}}$$ acts on

\begin{equation*}\mathcal{D}\,{:}\,{\raise-1.5pt{=}}\,\Big\{f\in C^2(\mathbb{S})\cap C_b(\mathbb{S}): \mathcal{G}f\in C_b\left(\mathbb{S}\right),\mu\left(l\right)f^\prime\left(l\right)+\frac{1}{2}\sigma^2\left(l\right)f^{\prime\prime}\left(l\right)=\rho f^\prime\left(l\right)\Big\}\end{equation*}

as follows:

(2.5) \begin{equation}\mathcal{G}f(x)=\left(\mu(x)I\left(x>l\right)+\rho I\left(x=l\right)\right)f^\prime(x)+\frac{1}{2}\sigma^2(x)I\left(x>l\right)f^{\prime\prime}(x)-k(x)f(x).\end{equation}

In particular, for $f\in\mathcal{D}$ ,

(2.6) \begin{equation}\mathcal{G}f\left(l\right)=\rho f^\prime\left(l\right)-k\left(l\right)f\left(l\right).\end{equation}

Remark 1. The Feynman–Kac semigroup can also be seen as the transition semigroup of a process $\tilde{X}$ which is killed at rate k; i.e., $\tilde{X}$ has the same drift and diffusion coefficients as X but is killed at its lifetime $\tilde{\zeta}=\inf\{t\geq 0:\int_0^tk(X_u)du\geq e\}$ , where $e\sim\textrm{Exp}(1)$ is exponentially distributed with rate 1. At the lifetime $\tilde{\zeta}$ the process is sent to the cemetery state. It can be proved that

(2.7) \begin{equation}\mathcal{P}_tf(x)=\mathbb{E}_x\left(f(\tilde{X}_t)I(\tilde{\zeta}>t)\right)\!.\end{equation}

See e.g. [Reference Linetsky36, Section 1.1] for a detailed explanation.

We can write

\begin{equation*}\mathcal{P}_tf(x)=\int_{\mathbb{S}}P\left(t,x,dy\right)f\left(y\right),\end{equation*}

where $P(t,x,\cdot)$ is the transition probability measure of the process on $\mathbb{S}$ . It has a point mass at l and a density on $\mathbb{S}\setminus\{l\}$ , denoted by $p(t,x,\cdot)$ . We can further write

(2.8) \begin{equation}\mathcal{P}_tf(x)=P\left(t,x,l\right)f\left(l\right)+\int_\mathbb{S}p\left(t,x,y\right)f\left(y\right)dy.\end{equation}

Note that we simply write $P\left(t,x,l\right)$ for $P\left(t,x,\{l\}\right)$ .

2.3. Eigenfunction expansion representations

Let S and M be the scale and speed measure of X. Using the results in [7, Section II], we can obtain that the scale measure has a density, i.e., $S(dx)=s(x)dx$ with

\begin{equation*}s(x)=\exp\Big(-\int_l^x\frac{2\mu\left(y\right)}{\sigma^2\left(y\right)}dy\Big),\end{equation*}

and the speed measure is of a mixed type, with

\begin{equation*} M\left(dx\right)=m(x)dx+\frac{1}{\rho}\delta_l\left(dx\right)=\frac{2}{\sigma^2(x)s(x)}dx+\frac{1}{\rho}\delta_l\left(dx\right),\end{equation*}

where $\delta_l$ is the Dirac delta measure at l. In future discussions, we will use M(l) instead of $M(\{l\})$ to simplify the notation.

It is possible to extend $${({{\mathcal P}_t})_{t \ge 0}}$$ to a self-adjoint semigroup on the Hilbert space $L^2(\mathbb{S},M)$ of square-integrable functions on $\mathbb{S}$ with respect to M. The inner product in this space is

\begin{equation*}\left(f,g\right)=\int_\mathbb{S}f(x)g(x)M\left(dx\right)=\frac{1}{\rho}f\left(l\right)g\left(l\right)+\int_{\mathbb{S}^\circ}\frac{2f(x)g(x)}{\sigma^2(x)s(x)}dx.\end{equation*}

The application of spectral methods to study diffusions dates back to [Reference McKean38]. When the spectrum of the generator of the diffusion is simple and purely discrete, a general spectral expansion reduces to an eigenfunction expansion. We make the following assumption for the existence of an eigenfunction expansion.

Under Assumption 2, X lives on $\mathbb{S}=[l,r)\cup\{\partial\}$ , and we extend the definition of k and the payoff function f to $\partial$ by setting $k(\partial)=f(\partial)=0$ in (2.4). Theorem 3.2 in [Reference Linetsky36] shows that the spectrum of $\mathcal{G}$ is purely discrete and simple, and we have the following eigenfunction expansions:

(2.9) \begin{align}P\left(t,x,l\right)&=M\left(l\right)\sum\limits_{k=1}^\infty\exp\left(-\lambda_kt\right)\varphi_k(x)\varphi_k\left(l\right)\qquad&& \textrm{for}\ x\in\mathbb{S}, \ t>0, \end{align}

(2.10) \begin{align}p\left(t,x,y\right)&=m\left(y\right)\sum\limits_{k=1}^\infty\exp\left(-\lambda_kt\right)\varphi_k(x)\varphi_k\left(y\right)\qquad&& \textrm{for}\ x,y\in\mathbb{S},\ y\neq l,\ t>0,\\u(t,x)&=\sum\limits_{k=1}^\infty\left(f,\varphi_k\right)\exp\left(-\lambda_kt\right)\varphi_k(x)\qquad&&\text{for}\ f\in L^2(\mathbb{S},M),\ x\in\mathbb{S},\ t>0,\nonumber \end{align}

where $$({\lambda _k},{\varphi _k})$$ is the kth eigenpair which solves the Sturm–Liouville problem

(2.11) \begin{align}\begin{split}&\mathcal{G}\varphi(x)=-\lambda\varphi(x) \qquad\qquad\qquad\ \textrm{for all}\ x\in\mathbb{S}^\circ,\\&\rho\varphi^\prime\left(l\right)-\left(k\left(l\right)-\lambda\right)\varphi\left(l\right)=0 \qquad\varphi\left(r\right)=0.\end{split}\end{align}

The eigenfunctions satisfy $\int_l^r\varphi_i(x)\varphi_j(x)M(dx)=\delta_{ij}$ (here $\delta_{ij}$ is the Kronecker delta), and hence the solutions of (2.11) form a complete orthonormal basis of $L^2(\mathbb{S},M)$ .

The eigenfunction expansion method has been extensively applied in finance for derivatives pricing. See for example [Reference Davydov and Linetsky14], [Reference Li and Linetsky31], [Reference Li and Linetsky32], [Reference Li and Linetsky33], and [Reference Linetsky36], among many others. The Sturm–Liouville problem associated with the sticky diffusion is different from those studied in the cited papers in that the eigenparameter also appears in the boundary condition, which makes the problem more difficult. The following proposition provides the asymptotic behavior of the eigenvalues and eigenfunctions of this type of Sturm–Liouville problem.

Proposition 1. Under Assumption 2, there exist constants $C_1,C_2,C_3>0$ such that for $k=1,2,\ldots$ ,

(2.12) \begin{align}{C_1}{k^2} \le {\lambda _k} \le {C_2}{k^2}, \end{align}

(2.13) $${\left\| {\varphi _k^{(i)}} \right\|_\infty } \le {C_3}{k^i},for\,{\mkern 1mu} i = 0,1,...,4.$$

Here $\varphi_k^{(i)}(x)$ is the ith-order derivative of $\varphi_k(x)$ , and $\varphi_k^{(0)}(x)=\varphi_k(x)$ . This result is proved by transforming the problem into the Liouville normal form and using the existing results in [Reference Altinisik, Kadakal and Mukhtarov1]. These results help to establish the bounds on the eigenvalues, which are then used to bound the norms of the eigenfunctions of the transformed problem. Finally, the bounds are preserved after transforming the problem back to its original formulation.

3.1. Transition rates of the CTMC

To obtain the transition rates of the CTMC that resembles the sticky diffusion $\tilde{X}$ , we apply finite differences to discretize its generator given by (2.5) and (2.6). Introduce the following difference operators acting on a general function g:

\begin{alignat*}{2}\nabla^+g(x)&=\frac{g\left(x^+\right)-g(x)}{x^+-x}, \qquad\qquad\qquad\quad\ \nabla^-g(x)&&=\frac{g(x)-g\left(x^-\right)}{x-x^-}, \\ \nabla g(x)&=\frac{\delta^-x}{2\delta x}\nabla^+g(x)+\frac{\delta^+x}{2\delta x}\nabla^-g(x), \qquad\Delta g(x)&&=\frac{1}{\delta x}\left(\nabla^+g(x)-\nabla^-g(x)\right).\end{alignat*}

For $x\in\mathbb{S}_n^\circ$ , we approximate $\mathcal{G}g(x)$ for $g\in\mathcal{D}$ by approximating $g^{\prime}$ by $\nabla g$ and $g^{\prime\prime}$ by $\Delta g$ , which gives us

(3.1) \begin{equation}G_ng(x)=\mu(x)\nabla g(x)+\frac{1}{2}\sigma^2(x)\Delta g(x)-k(x)g(x),\qquad\textrm{for}\ x\in\mathbb{S}_n^\circ. \end{equation}

For $\mathcal{G}g(x_0)$ , we apply $\nabla^+g\left(x_0\right)$ to approximate $g^{\prime}(x_0)$ and obtain

\begin{equation*}G_ng\left(x_0\right)=\rho\nabla^+g\left(x_0\right)-k\left(x_0\right)g\left(x_0\right).\end{equation*}

We specify $x_{n+1}$ to be a killing boundary for the CTMC; i.e., the chain is sent to the cemetery state $\partial$ once it reaches $x_{n+1}$ . So the state space of the chain is $\{x_0,\cdots,x_n,\partial\}$ . Let $\mathbb{G}_n$ be an $(n+1)\times(n+1)$ matrix for the transition rates among the states $\{x_0,\cdots,x_n\}$ . Then, based on the expression for $G_n$ ,

(3.2) \begin{equation}\mathbb{G}_n=\begin{pmatrix} -\frac{\rho}{\delta^+ x_0}-k\left(x_0\right) &\quad \frac{\rho}{\delta^+ x_0} &\quad 0 &\quad \ldots \\ \\ \mathbb{G}_{n,1,0} &\quad -\mathbb{G}_{n,1,0}-\mathbb{G}_{n,1,2}-k\left(x_1\right) &\quad \mathbb{G}_{n,1,2} &\quad \ldots \\ \\ 0 &\quad \mathbb{G}_{n,2,1} &\quad -\mathbb{G}_{n,2,1}-\mathbb{G}_{n,2,3}-k\left(x_2\right) &\quad \ldots \\ \\ \vdots &\quad \vdots &\quad \vdots &\quad \ddots \end{pmatrix},\end{equation}

which is a tridiagonal matrix, and

\begin{alignat*}{2}\mathbb{G}_{n,i,i-1}&=\frac{-\mu\left(x_i\right)\delta^+x_i+\sigma^2\left(x_i\right)}{2\delta^-x_i\delta x_i},\qquad&& i=1,\ldots,n,\\ \mathbb{G}_{n,i,i+1}&=\frac{\mu\left(x_i\right)\delta^-x_i+\sigma^2\left(x_i\right)}{2\delta^+x_i\delta x_i},\qquad&& i=1,\ldots,n-1,\\ \mathbb{G}_{n,0,1}&=\frac{\rho}{\delta^+x_0}.\end{alignat*}

Note that in our calculation of (2.4) we do not need the transition rates involving $\partial$ . We refer to this construction of $\mathbb{G}_n$ in (3.2) as Scheme 1.

It turns out that Scheme 1 only converges at first order. Below we provide a better scheme that improves approximation for the sticky behavior. In Scheme 2, we try to match the expectation and variance of the sticky diffusion in a short time, given that it starts from $x_0$ , with those of the CTMC. Let $\hat{O}_t^{x_0}$ approximate the occupation time of $\tilde{X}$ at the left boundary $l=x_0$ up to time t, where t is understood to be very small. Then, after ignoring terms of higher order in t, we have the system of equations

(3.3) \begin{align}\mathbb{G}_{n,0,1}\delta^+x_0t&=\mu\left(x_0\right)\left(t-\hat{O}_t^{x_0}\right)+\rho\hat{O}_t^{x_0}, \end{align}

(3.4) \begin{align}\mathbb{G}_{n,0,1}\left(\delta^+x_0\right)^2t&=\sigma^2\left(x_0\right)\left(t-\hat{O}_t^{x_0}\right),\end{align}

together with the condition that $\mathbb{G}_{n,0,0}+\mathbb{G}_{n,0,1}+k(x_0)=0$ . Equation (3.3) matches the expected change from $x_0$ in t. The left-hand side is clearly the quantity for the CTMC, ignoring higher-order terms in t. The right-hand side can be explained as follows. Up to time t, the diffusion $\tilde{X}$ spent a length of time $\hat{O}^{x_0}_t$ at the boundary $x_0$ with drift $\rho$ and a length of time $t-\hat{O}^{x_0}_t$ near $x_0$ with drift approximately equal to $\mu(x_0)$ . Adding these up gives the right-hand side of (3.3). Equation (3.4) can be interpreted in an analogous way by noting that the change of $\tilde{X}$ does not have any variance while it is on the boundary.

Solving the above equations gives the following result:

\begin{align*}\mathbb{G}_{n,0,1}&= \frac{\rho}{\delta^+x_0+\frac{\rho-\mu\left(x_0\right)}{\sigma^2\left(x_0\right)}\left(\delta^+ x_0\right)^2},\qquad\mathbb{G}_{n,0,0}=-\mathbb{G}_{n,0,1}-k\left(x_0\right), \\ \hat{O}_t^{x_0}&=\frac{\sigma^2(x_0) - \mu(x_0) \delta^+ x_0}{\sigma^2(x_0) + (\rho - \mu(x_0)) \delta^+ x_0} t.\end{align*}

The above expression shows that $t-\hat{O}_t^{x_0}=O(\delta^+x_0)$ . Ignoring this difference in (3.3) yields the formula of $\mathbb{G}_{n,0,1}$ under Scheme 1. Now it is intuitively clear that Scheme 2 is better than Scheme 1, as it matches the first two moments, while Scheme 1 only cares about the first moment. In Section 4, we will prove that Scheme 2 achieves second-order convergence.

The operator $G_n$ acting on $g(x_0)$ can be written in a unified way for both schemes as follows:

\begin{equation*}G_ng\left(x_0\right)=\rho\beta\nabla^+g\left(x_0\right)-k\left(x_0\right)g\left(x_0\right),\end{equation*}

where

(3.5) \begin{equation}\beta=\begin{cases}1,&\quad\text{Scheme 1},\\ \frac{1}{1+\frac{\rho-\mu\left(x_0\right)}{\sigma^2\left(x_0\right)}\delta^+x_0},&\quad\text{Scheme 2}.\end{cases}\end{equation}

3.2. CTMC approximation of the Feynman–Kac semigroup and first passage probabilities

To approximate (2.4), we use $\mathbb{E}_x\left(f\left(Y_t\right)I\left(\zeta_Y>t\right)\right)$ , where $\zeta_Y$ is the lifetime of Y. This expectation is obtained in closed form in [Reference MijatoviĆ and Pistorius40]; it is the approximate value function given the payoff f and is calculated by

(3.6) \begin{equation}u_n(t,x)\,{:}\,{\raise-1.5pt{=}}\,\mathbb{E}_x\left(f\left(Y_t\right)I\left(\zeta_Y>t\right)\right)=\exp(\mathbb{G}_nt)f_n(x),\end{equation}

where $f_n=(f(x_0),\ldots,f(x_n))^T$ is an $(n+1)$ -dimensional column vector, and $\exp(\mathbb{G}_nt)f_n(x)$ in (3.6) is the entry corresponding to x in the vector $\exp(\mathbb{G}_nt)f_n$ .

One may also be interested in approximations of P(t,x,l) and p(t,x,y), where P is the distribution of $X_t$ given that the process starts at x, and p is the corresponding density on the interior of the state space. We now denote the transition probability of the CTMC from x to y in time t by $P_n(t,x,y)$ . It is given by

\begin{equation*}P_n(t,x,y)=\exp(\mathbb{G}_nt)(x,y),\qquad x,y\in\mathbb{S}_n^-,\end{equation*}

which is the entry of the matrix $\exp(\mathbb{G}_nt)$ with its row corresponding to x and its column to y. In particular, $P_n(t,x,l)$ approximates P(t,x,l). For $y\in\mathbb{S}_n^\circ$ , define

\begin{equation*}p_n(t,x,y):=P_n(t,x,y)/\delta y, \qquad x\in\mathbb{S}_n^-.\end{equation*}

Then we can approximate p(t,x,y) by $p_n(t,x,y)$ .

It is also of interest in practice to calculate the first passage probability $\mathbb{P}(\tau^X_z>t|X_0=x)$ for $l\leq x<z$ , where $\tau^X_z$ is the hitting time of z of the sticky diffusion X. The paper [Reference Nie and Linetsky43] obtains a semi-analytical formula for this probability under the sticky OU model. We use CTMC approximation to calculate the probability for general sticky diffusions. Let Y be the CTMC that approximates X (set $k\equiv 0$ in $\mathbb{G}_n$ ). Then, for Y, using [Reference MijatoviĆ and Pistorius40],

\begin{equation*}\mathbb{P}(\tau^Y_z>t|Y_0=x)=\exp(\mathbb{H}_nt)\textbf{1}(x),\end{equation*}

where $\tau^Y_z=\inf\{t: Y_t>z\}$ , $\mathbb{H}_n$ is the submatrix of $\mathbb{G}_n$ with entries corresponding to states smaller than z, and $\textbf{1}$ is a vector of ones with dimension compatible with $\mathbb{H}_n$ . In particular, setting $x=l$ , we obtain an approximation of the first passage probability when the diffusion starts from the sticky boundary.

To implement the proposed method, we need to compute the matrix exponential, for which a variety of algorithms exist. In this paper we evaluate the performance of the three approaches listed below:

• The scaling and squaring algorithm of [Reference Higham24] is used widely in the literature on Markov chain approximation. It is known that this algorithm can handle all kinds of matrices, and it is numerically stable. It requires $O(n^3)$ operations to compute $\exp(\mathbb{G}_nt)f_n$ .

• For the case when $\mathbb{G}_n$ is a tridiagonal matrix, [Reference Li and Zhang34] developed an algorithm based on the fast matrix eigendecomposition algorithm of [Reference Dhillon15], known as the MRRR algorithm. This algorithm takes only $O(n^2)$ operations to compute $\exp(\mathbb{G}_nt)f_n$ .

• The third approach is to numerically solve the ODE system that the matrix exponential $\mathbb{F}$ satisfies, which is

  $$d{\mathbb {F}}(t) = {\mathbb {G}}{_n}{\mathbb {F}}(t)dt,\qquad {\mathbb {F}}(0) = {I_{n + 1}}.$$
  We propose to apply the extrapolation approach of [Reference Feng and Linetsky19], which has not previously been considered in the literature on CTMC approximation.

Our experiment in Section 5.1 shows that the third approach performs the best, and we recommend using it for computing $\exp(\mathbb{G}_nt)$ , especially when t is large.

3.3. Eigenfunction expansions for the CTMC

We provide an eigenfunction expansion representation for $p_n(t,x,y)$ , which will be utilized later to analyze the convergence rate of CTMC approximation. Let

(3.7) \begin{align}M_n\left(x_0\right)&=M\left(x_0\right)\exp\Big(\frac{\alpha}{\sigma^2\left(x_0\right)}\delta^+x_0\Big), \end{align}

\begin{align}\frac{1}{s_n\left(x_0\right)}&=\rho M_n(x_0)=\exp\Big(\frac{\alpha}{\sigma^2\left(x_0\right)}\delta^+x_0\Big),\nonumber\end{align}

where

(3.8) \begin{equation}\alpha=\begin{cases}\mu(x_0),&\quad\text{Scheme 1},\\ \rho,&\quad\text{Scheme 2},\end{cases}\end{equation}

and for $x=x_1,\ldots,x_n$ ,

(3.9) \begin{align} m_n(x)&=M_n\left(x_0\right)\beta\frac{2\rho}{-\mu\left(x_1\right)\delta^+x_1+\sigma^2\left(x_1\right)}\prod_{y=x_1}^{x^-}\frac{\mu\left(y\right)\delta^-y+\sigma^2\left(y\right)}{-\mu\left(y^+\right)\delta^+y^++\sigma^2\left(y^+\right)},\end{align}

\begin{align}\frac{1}{s_n(x)}&=\frac{1}{2}m_n(x)(\mu(x)\delta^-x+\sigma^2(x)), \nonumber\end{align}

where $\alpha$ is as in (3.8), $\beta$ is as in (3.5), and the product term equals 1 if $x=x_1$ . We also define $M_n(x)=m_n(x)\delta x$ . As will be shown later, $m_n(x)\approx m(x)$ for $x\in\mathbb{S}_n^\circ$ , and (3.7) shows that $M_n(x_0)\approx M(x_0)$ . Furthermore, the last two terms in (3.9) are an approximation of $\int_{x_1}^x2\mu(y)/\sigma^2(y)dy$ , which appears in the definition of the scale density and speed density of the diffusion X. Hence, one can show that $M_n$ , $m_n$ , and $s_n$ are candidates to approximate the speed measure, speed density, and scale density of the diffusion.

Recall the expression (3.1) for $G_n$ , the generator of the CTMC. For $x\in\mathbb{S}_n^\circ$ , we can rewrite it in the form

\begin{equation*}G_ng(x)=\frac{1}{m_n(x)}\frac{\delta^-x}{\delta x}\nabla^-\Big(\frac{1}{s_n(x)}\nabla^+g(x)\Big)-k(x)g(x),\end{equation*}

which is crucial for the convergence rate analysis.

Consider the eigenvalue problem of $G_n$ :

(3.10) \begin{align}\begin{split}&G_n\varphi^n(x)=-\lambda\varphi^n(x),\qquad x=x_0,\ldots,x_n,\\&\varphi^n(x_{n+1})=0,\end{split}\end{align}

where $\varphi^n(x)$ is defined on $\mathbb{S}_n$ , i.e., it is a vector. Let $(\lambda^n_k,\varphi^n_k)$ be the kth eigenpair. The eigenvalue problem of the operator $G_n$ in (3.10) is essentially the eigenvalue problem of the matrix $\mathbb{G}_n$ , which is tridiagonal with negative diagonal elements and positive off-diagonal ones, and also diagonally dominant. Thus, it has $n+1$ simple and real eigenvalues with $0\leq\lambda_1^n<\lambda_2^n<\dots$ (see [Reference Li and Zhang34, Proposition 3.6]), and one expects $\lambda^n_k\approx\lambda_k$ . For any two functions $g_1,g_2$ defined on $\mathbb{S}_n$ , define their inner product $(g_1,g_2)_n\,{:}\,{\raise-1.5pt{=}}\,\sum_{x\in\mathbb{S}_n^-}g_1(x)g_2(x)M_n(x)$ . Since eigenfunctions are only unique up to a constant, we normalize them to satisfy $${(\varphi _j^n,\varphi _k^n)_n} = {\delta _{j,k}}$$ , so that $\varphi_k^n(x_i)\approx\varphi_k(x_i)$ .

We have the following bilinear eigenfunction expansion representations:

(3.11) \begin{align}{1}P_n(t,x,y)&=M_n(y)\sum_{k=1}^n\exp(-\lambda_k^nt)\varphi_k^n(x)\varphi_k^n(y),\qquad\quad&& x\in\mathbb{S}_n^-, y\in\mathbb{S}_n^-,\end{align}

(3.12) \begin{align}{1}p_n\left(t,x,y\right)&=m_n\left(y\right)\sum_{k=1}^n\exp\left(-\lambda_k^nt\right)\varphi_k^n(x)\varphi_k^n\left(y\right),\qquad&& x\in\mathbb{S}_n^-, y\in\mathbb{S}_n^\circ. \end{align}

3.4. CTMC approximation of 1D diffusions with a sticky interior point

We show how to construct a CTMC approximation for general diffusions with a sticky interior point $S^{*}\in\mathbb{S}^\circ$ (the boundary points are not sticky). To this end, one can generalize the steps in [Reference Jiang, Song and Wang26] and obtain the generator of this new diffusion as

\begin{equation*}\mathcal{G}f(x)=\begin{cases}\mu(x)f^\prime(x)+\frac{1}{2}\sigma^2(x)f^{\prime\prime}(x)-k(x)f(x), &\qquad x\neq S^{*}, \\ \rho\left(f^\prime\left(x+\right)-f^\prime\left(x-\right)\right)-k(x)f(x), & \qquad x=S^{*}, \end{cases}\end{equation*}

with domain

\begin{align*}\mathcal{D}&=\bigg\{f\in C^2\left(\mathbb{S}\setminus\{S^{*}\}\right)\cap C_b\left(\mathbb{S}\right):\mathcal{G}f\in C_b\left(\mathbb{S}\right), \\&\qquad\qquad\rho\left(f^\prime\left(S^{*}+\right)-f^\prime\left(S^{*}-\right)\right)-k\left(S^{*}\right)f\left(S^{*}\right)=\lim_{x\to S^{*}}\mathcal{G}f(x)\bigg\}.\end{align*}

By setting $S^{*}=l$ , we recover the sticky lower boundary case as the left limit $f^\prime(S^{*}-)=0$ . The CTMC that approximates this diffusion can be constructed as follows, and we call it Scheme 1. First, we put $S^{*}$ on the grid, and the index for this state is denoted by j; i.e., $x_j=S^{*}$ . Second, the transition rate matrix $\mathbb{G}_n$ is adjusted at $x_j$ as follows:

\begin{equation*}\mathbb{G}_{n,j,j-1}=\frac{\rho}{\delta^-S^{*}}, \qquad\mathbb{G}_{n,j,j}=-\rho\frac{2\delta S^{*}}{\delta^-S^{*}\delta^+S^{*}}-k\left(S^{*}\right), \qquad\mathbb{G}_{n,j,j+1}=\frac{\rho}{\delta^+S^{*}}.\end{equation*}

At the left boundary, the transition rate $\mathbb{G}_{n,0,1}=0$ is set according to the behavior of the boundary point. See Remark 3 for the absorbing and reflecting cases.

We can also develop a scheme that achieves second-order convergence, which we call Scheme 2. Applying Taylor expansion yields

(3.13) \begin{align}\nabla^+ f\left(S^{*}\right)&= f^\prime\left(S^{*}+\right) + \frac{1}{2} f^{\prime\prime}\left(S^{*}+\right) \delta^+S^{*} + O\left(h_n^2\right), \end{align}

(3.14) \begin{align}\nabla^- f\left(S^{*}\right)&= f^\prime\left(S^{*}-\right) - \frac{1}{2} f^{\prime\prime}\left(S^{*}-\right) \delta^-S^{*} + O\left(h_n^2\right).\end{align}

Since $f\in \mathcal{D}$ , we also have

\begin{align*}f^{\prime\prime}\left(S^{*}+\right)&=\frac{\rho\left(f^\prime\left(S^{*}+\right)-f^\prime\left(S^{*}-\right)\right)-\mu\left(S^{*}+\right)f^\prime\left(S^{*}+\right)}{\sigma^2\left(S^{*}+\right)/2},\\ f^{\prime\prime}\left(S^{*}-\right)&=\frac{\rho\left(f^\prime\left(S^{*}+\right)-f^\prime\left(S^{*}-\right)\right)-\mu\left(S^{*}-\right)f^\prime\left(S^{*}-\right)}{\sigma^2\left(S^{*}-\right)/2}.\end{align*}

Substituting these two equations in (3.13) and (3.14), we obtain

\begin{align*}\nabla^+ f(S^{*}) &= \left( 1 + \frac{\rho - \mu(S^{*}+)}{\sigma^2(S^{*}+)} \delta^+S^{*} \right)f'(S^{*}+) - \frac{\rho\delta^+S^{*}}{\sigma^2(S^{*}+)} f'(S^{*}-) + O(h_n^2), \\ \nabla^- f(S^{*}) &= \left( 1 + \frac{\rho + \mu(S^{*}-)}{\sigma^2(S^{*}-)} \delta^-S^{*} \right)f'(S^{*}-) - \frac{\rho\delta^-S^{*}}{\sigma^2(S^{*}-)} f'(S^{*}+) + O(h_n^2).\end{align*}

Let

\begin{align*}A &= 1 + \frac{\rho - \mu(S^{*}+)}{\sigma^2(S^{*}+)} \delta^+S^{*},\\ B &= \frac{\rho\delta^+S^{*}}{\sigma^2(S^{*}+)},\\ C &= 1 + \frac{\rho + \mu(S^{*}-)}{\sigma^2(S^{*}-)} \delta^-S^{*},\\ D &= \frac{\rho\delta^-S^{*}}{\sigma^2(S^{*}-)};\end{align*}

then

\begin{align*}f'(S^{*}+)&= \frac{C\nabla^+f(S^{*}) + B\nabla^-f(S^{*})}{AC-BD} + O(h_n^2),\\ f'(S^{*}-)&= \frac{D\nabla^+f(S^{*}) + A\nabla^-f(S^{*})}{AC-BD} + O(h_n^2).\end{align*}

Therefore,

\begin{align*}\rho \left( f'(S^{*}+) - f'(S^{*}-) \right) - k(S^{*}) f(S^{*}) &= \frac{\rho (C - D)}{AC - BD} \nabla^+ f(S^{*}) + \frac{\rho(B - A)}{AC-BD} \nabla^- f(S^{*}) \\ &\qquad - k(S^{*})f(S^{*}) + O(h_n^2).\end{align*}

The corresponding elements of the generator matrix are then given by

\begin{align*}\mathbb{G}_{n,j,j-1} &= \frac{\rho (A - B)}{(AC - BD)\delta^-S^{*}}, \qquad\mathbb{G}_{n,j,j+1} = \frac{\rho (C - D)}{(AC - BD)\delta^+S^{*}}, \\ \mathbb{G}_{n,j,j} &= -\mathbb{G}_{n,j,j-1}-\mathbb{G}_{n,j,j+1} - k(S^{*}).\end{align*}

Remark 4. Now suppose we have p sticky points, and the set of such points is $I_{S^{*}}=\{S_1^{*},\ldots,S_p^{*}\}$ , with $\rho_i$ as the stickiness parameter at $S^*_i$ . The generator of the process is given by

\begin{equation*}\mathcal{G}f(x)=\begin{cases}\mu(x)f^\prime(x)+\frac{1}{2}\sigma^2(x)f^{\prime\prime}(x)-k(x)f(x), &\qquad x\notin I_{S^{*}}, \\ \rho_i\left(f^\prime\left(x+\right)-f^\prime\left(x-\right)\right)-k(x)f(x), & \qquad x=S_i^{*},\quad i=1,\ldots,p, \end{cases}\end{equation*}

with domain

\begin{align*}\mathcal{D}&=\bigg\{f\in C^2\left(\mathbb{S}\setminus I_{S^{*}}\right)\cap C_b\left(\mathbb{S}\right):\mathcal{G}f\in C_b\left(\mathbb{S}\right), \\ &\qquad\qquad\rho_i\left(f^\prime\left(S_i^{*}+\right)-f^\prime\left(S_i^{*}-\right)\right)-k\left(S_i^{*}\right)f\left(S_i^{*}\right)=\lim_{x\to S^*_i}\mathcal{G}f\left(S_i^{*}\right),\ i=1,\ldots,p\bigg\}.\end{align*}

A sticky left or right boundary can be included in this framework by setting $f^\prime\left(S_i^{*}-\right)=0$ or $f^\prime\left(S_i^{*}+\right)=0$ . One can then obtain the CTMC transition rates at each sticky point using finite differences in the same way as the single sticky point case.

5.1. Numerical results for the sticky OU short-rate model

The sticky OU process is given by the following SDE:

\begin{equation*}dX_t=\left(\kappa\left(\theta-X_t\right)I\left(X_t>0\right)+\rho I\left(X_t=0\right)\right)dt+\sigma I\left(X_t>0\right)dB_t,\end{equation*}

where $l=0$ is the left boundary. This process is used in [Reference Nie42] (and also in [Reference Nie and Linetsky43]) to model short rates, which are instantaneous interest rates. The advantages of using this model over standard short-rate models are explained in [Reference Nie42]. In particular, this model is able to produce various shapes of yield curves observed in the market, including the S-shaped yield curve in a low-interest environment seen after the 2008 financial crisis. To illustrate this point, we consider three different sets of values for the parameters of a sticky OU process, which are listed in Table 1. These sets are able to produce three different shapes of yield curves: upward-sloping (Model 1), inverted (Model 2), and S-shaped (Model 3). In general, the magnitude of the stickiness parameter plays an important role in controlling the shape of the curve. Calibration shows that the implied stickiness parameter from real market data is on similar scales with the values in the table (see [Reference Nie42]).

Table 1: Model parameters for the sticky OU process.

Throughout this section, we consider pricing zero-coupon bonds with unit payoff, i.e. $u(t,x)=\mathbb{E}_x(\exp(-\int_0^tX_udu))$ with $k(x)=x$ . The paper [Reference Nie and Linetsky43] obtained an eigenfunction expansion formula for the bond price under the sticky OU model. We use the method of [Reference Nie and Linetsky43] to obtain benchmarks for our method. Table 2 shows the absolute difference between the results of the CTMC method with 20,000 grid points and the eigenfunction expansion method using 50 to 100 terms in the expansion.

Table 2: Differences in prices for the CTMC method and the eigenfunction expansion method of [Reference Nie42, Section 6].

It can be seen that these methods yield consistent results for all maturities shown in the table. The eigenfunction expansion method can be slow for small maturities (e.g., 1 or 3 months). This is because the nth term in the expansion involves $e^{-\lambda_n t}$ (see [Reference Nie and Linetsky43, Proposition 4.2]), where $\lambda_n$ is the nth eigenvalue and t is the time to maturity. This decays slowly for small t, so more terms in the expansion are required to achieve good accuracy.

5.2. Convergence rates

We localize the state space $[0,\infty)$ to [0,1]. Simple calculations show that for the OU process with parameter values in Table 1, the probability of breaching this upper boundary is extremely small for all maturities under consideration, and hence the localization error is negligible. We then discretize [0,1] with $n=100, 200, 400, 800, 1600$ points using a uniform grid and apply the squaring and scaling algorithm to calculate the matrix exponential. The convergence rate is detected numerically. Figure 1 plots the results for three models with two maturities. It is clear that for Models 1 and 3, Scheme 2 converges at second order while Scheme 1 only attains first order; this can be seen by checking the slope of the convergence line against that of the small triangle in the lower left corner of each plot. This validates our theoretical estimates in Theorem 4. For Model 2, the convergence orders of Scheme 1 and Scheme 2 are both around two. This may at first seem unexpected, but the value of $\rho$ is very small in this model (much smaller than in the other two models), which makes $\mathbb{G}_{n,0,1}$ roughly zero under both schemes. So the two schemes perform similarly.

Figure 1: Absolute error vs. n on log-log scale for Models 1–3 (from top to bottom) with six-month and 30-year maturities.

We can further observe that in Models 2 and 3, Scheme 1 is more accurate than Scheme 2 for small n. The reason is that $\mathbb{G}_{n,0,1}$ is always positive in Scheme 1, whereas it can become negative in Scheme 2. In the latter case, the matrix $\mathbb{G}_n$ is not a proper transition rate matrix, which leads to larger errors near the boundary. Finally, for all models it is observed that the absolute error is greater when the maturity is longer.

A practical question is which scheme should be used. The choice depends on the parameter values. When $\rho$ is extremely small (like the value of $\rho$ in Model 2), we recommend Scheme 1 for the reasons mentioned above. For moderate and big values of $\rho$ , one can always implement Scheme 2 if highly accurate results are needed.

5.3. Some comparisons

The calculation of the value function under the CTMC model is based on the calculation of the matrix exponential. In this section, we compare three ways of computing the matrix exponential in our problem (the three ways listed in Section 3.2), and we also compare the CTMC approximation algorithm with a standard finite difference scheme.

We now briefly describe the extrapolation approach of [Reference Feng and Linetsky19] for numerically solving the ODE system satisfied by the matrix exponential. This approach uses a basic step size H, where $H=0.5$ if the maturity T is greater than $0.5$ and $H=T$ otherwise, to divide the interval [0,T] into smaller time periods. For each basic interval, $M_i=1,\ldots,s$ time steps are used to evolve the differential equation according to the implicit scheme to calculate the matrix exponential, where $s\geq 1$ denotes the extrapolation stage number. Let the approximation of the solution after one basic step be denoted by $A_{i,1}=u_n(H,x; M_i)$ where $M_i$ time points are used for the interval [0,H]. An extrapolation tableau is constructed by the equation

\begin{equation*}A_{i,j}=A_{i-1,j-1}+\frac{A_{i,j-1}-A_{i-1,j-1}}{\frac{M_i}{M_{i-j+1}}-1},\end{equation*}

for $i=2,\ldots,s$ and $j=2,\ldots,i$ . We then use the value $A_{s,s}$ after s extrapolation steps as the starting point of the approximation over the time interval [H,2H]. This calculation of approximations after basic steps is repeated until one obtains an approximation of $u_n(T,x)$ after M basic steps, where $T=MH$ . Finally, the value of $A_{s,s}$ after M steps is the approximation of the zero-coupon bond price with maturity T at time 0.

For the finite difference scheme, which discretizes both time and space in the PDE satisfied by the value function of the sticky diffusion, we use Crank–Nicolson time-stepping, which is a standard choice in the literature for numerical solutions of PDEs. It is also considered in [Reference Nie42]. We use equal time steps with 5 steps in a month. This allows for an adequate number of time steps even for longer maturities.

Figure 2: Comparison of four methods. ‘CTMC expm’ stands for using the scaling and squaring algorithm for computing the matrix exponential in the CTMC method. In the extrapolation approach, extrapolation is only applied once (i.e., $s=2$ ).

The MRRR algorithm, the extrapolation approach, and the finite difference method are implemented in C++, whereas the scaling and squaring method is already implemented in Matlab so we directly call it in Matlab. Although a Matlab implementation may generally be less efficient than an equivalent C++ implementation, it does not affect the conclusion that the scaling and squaring method is typically the slowest for obtaining similar levels of accuracy, as it involves the highest amount of computational complexity.

The results are obtained using a workstation running CentOS 7 with 3.2 GHz CPUs (32 cores) and memory of 256 GB. We did not make use of parallelization or other methods to speed up computations. The workstation was used because a larger amount of memory is needed to store the matrix exponential when calculating the benchmark prices with 20,000 grid points. All numerical calculations were run 10 times to obtain the average running time. Figure 2 displays a comparison of the three models for bonds with one-year and 30-year maturities.

The comparison clearly favors the extrapolation method, which defeats the other three methods in all cases, and its leading edge becomes greater as the bond’s maturity increases. It should be noted that the eigendecomposition method based on the MRRR algorithm does not always work. For Model 2, it cannot be applied when n is too large, and for Model 3, it fails for $n=100$ because of overflow/underflow errors in calculating a similar symmetric tridiagonal matrix required by the algorithm.

5.4. Simulation

We also consider how to perform simulations for sticky diffusions. Our CTMC approximation method offers a natural alternative to the standard Euler scheme. The simulation of sample paths from a CTMC is straightforward, and unlike the Euler scheme, it does not require time-discretization. We use 500 grid points for the CTMC Y constructed by Scheme 2, start with $Y_0=x$ , and then draw an exponentially distributed random variable with intensity $\vert\mathbb{G}_{n,x,x}\vert$ to determine the amount of time spent in the initial state. After determining the time point when the Markov chain is transitioning, we use the transition rates $\mathbb{G}_{n,x,x^-}$ , $\mathbb{G}_{n,x,x^+}$ , and k(x) to sample the new state. These steps are repeated until maturity is reached. Using the definition of the order of weak convergence given in [Reference Glasserman21] and [Reference Kloeden and Platen29] together with Theorem 4, we find that the weak convergence order of the CTMC simulation scheme using the transition rate matrix under Scheme 2 is two (it is one if Scheme 1 is used).

The Euler scheme is implemented in the following way. Time is discretized using 50 time steps per month. The process starts with $X_0=x$ , and subsequent values of the process are computed using the discretization of the SDE (2.3). In particular, for $t=0,\Delta t,\ldots,T-\Delta t$ ,

\begin{align*}Z_{t+\Delta t}&=\begin{cases} X_t+\mu\left(X_t\right)\Delta t+\sigma\left(X_t\right)\sqrt{\Delta t}\xi_{t+\Delta t} & \qquad\textrm{if}\ X_t>0, \\ \rho\Delta t & \qquad\textrm{if}\ X_t=0, \end{cases} \\e_{t+\Delta t}&\sim\textrm{Exp}\left(1\right),\end{align*}

where $\xi_{t+\Delta t}$ is standard normally distributed and $e_{t+\Delta t}$ is exponentially distributed with rate 1. The new value $Z_{t+\Delta t}$ is accepted as given by

\begin{equation*}X_{t+\Delta t}=\begin{cases} Z_{t+\Delta t} &\qquad\textrm{if}\ Z_{t+\Delta t}>0,\ e_{t+\Delta t}>k\left(X_t\right)\Delta t, \\ 0 &\qquad\textrm{if}\ Z_{t+\Delta t}\leq 0,\ e_{t+\Delta t}>k\left(X_t\right)\Delta t, \\ \partial &\qquad\textrm{if}\ e_{t+\Delta t}\leq k\left(X_t\right)\Delta t, \end{cases}\end{equation*}

where $\partial$ is the cemetery state and $f(\partial)=0$ .

Figure 3 shows the Monte Carlo simulation results when sample paths are simulated by a CTMC and the Euler scheme. In both cases, 1000 samples paths are simulated and the price estimator together with the 99% confidence intervals are displayed.

Figure 3: Monte Carlo simulation results for zero-coupon bond pricing. The blue line indicates the benchmark prices. The left panel shows the results for Model 1 with $\rho=4.0\times 10^{-3}$ . In the right panel, the results are obtained by setting $\rho=0.1$ with all other parameters remaining fixed.

The Euler scheme clearly fails when $\rho$ is small, i.e., when the process is very sticky. However, it becomes acceptable when $\rho$ is big enough. In contrast, the CTMC simulation scheme works well regardless of the degree of stickiness.

We provide some intuition about why the Euler scheme flops in very sticky cases. Note that this method only simulates a discrete-time process. If at some time point on the grid, say $i\Delta t$ , the process is at zero, it moves to $\rho \Delta t$ at $$(i + 1)\Delta t$$ . If $\Delta t$ is small enough, the interpolated path of the Euler scheme should resemble the path of the continuous-time process to which the Euler scheme converges in the limit. Thus, it is intuitively clear that the limiting continuous-time process is instantaneously reflected at zero and hence not sticky there. When the original diffusion is only mildly sticky, it is not very different from the reflected case, so the Euler scheme produces acceptable results. However, if the original diffusion is very sticky, the difference from the reflected case is large and the results of the Euler scheme become useless. Figure 4 shows one path simulated from the CTMC scheme and the other from the Euler scheme. The CTMC simulation scheme does not discretize time, so it can generate the phenomenon that the process sticks to zero, whereas the Euler scheme cannot.

Figure 4: Sample paths over 10 years generated by the CTMC and Euler scheme. The sample path for the Euler scheme was simulated using $\Delta t=1/600$ , i.e., with 50 time steps per month.

Figure 5: Absolute error vs. n on log-log scale for call option prices with one-year maturity obtained with the CTMC method: $S^{*}=110$ (left) and $S^{*}=90$ (right). The geometric Brownian motion is localized to [0,1000], and the scaling and squaring algorithm is used to compute matrix exponentials.

5.5. Geometric Brownian motion with a sticky interior point

We consider the model proposed in [Reference Jiang, Song and Wang26], where the asset price follows $dX_t=\mu(X_t)dt+\sigma(X_t)dB_t$ with $\mu(x)=(\mu_1I(x<S^*)+\mu_2I(x\geq S^*))x$ and $\sigma(x)=(\sigma_1I(x<S^*)+\sigma_2I(x\geq S^*))x$ . We set $\mu_1=0.1$ , $\mu_2=0.2$ , $\sigma_1=0.3$ , and $\sigma_2=0.4$ , and the process starts at $X_0=100$ . We consider two cases, $S^{*}=110$ and $S^{*}=90$ , and calculate the price of a European call option with strike price $\xi=100$ . The risk-free rate is 5%, and time to maturity equals one year. The stickiness at $S^{*}$ is given by $\rho=1/(2\times 0.02)=25$ . To obtain benchmark prices, we apply the method in [Reference Jiang, Song and Wang26] to numerically invert the Laplace transform of the option price which was derived by the authors. However, for more general diffusion models, explicit formulas for the Laplace transform may not be available.

To implement the CTMC approximation method, we use the grid design in [Reference Zhang and Li55] which places the strike price K midway between two adjacent grid points. This would make convergence smooth for call options. In our design, we also place the sticky point $S^*$ on the grid. Figure 5 shows convergence results for the two cases, from which the convergence orders for Scheme 1 and Scheme 2 are confirmed to be one and two. It should be noted that the initial price $X_0$ at which the option is priced is not on the grid, as $X_0=K$ . In this case, we used cubic splines to interpolate the prices on the neighboring grid points to obtain the price at $X_0$ . This does not reduce the convergence order of CTMC approximation, as a cubic spline interpolation has higher convergence order.