An Analytically Modified Finite Difference Scheme for Pricing Discretely Monitored Options

1. Introduction

Discretely monitored exotic options, including barrier options and autocallable structured products, are commonly traded throughout the world. Unlike continuously monitored barrier options, discrete barrier options cannot be priced using analytical formulas. Instead, these options are priced using numerical techniques such as Monte-Carlo simulations, finite difference methods, quadrature methods [1,2], and several other sophisticated mathematical techniques (e.g. [3,4,5,6,7]). Discrete barrier options can also be priced approximately using continuity correction techniques [8,9], but the accuracy of these methods is often quite limited.

To price options with low to moderate monitoring frequencies (e.g. monthly), quadrature methods can be quite effective [2]. However, these methods are not suitable for pricing options with high monitoring frequencies (e.g. daily), since they involve computations of convolution integrals at all monitoring points (or barriers), which can be very expensive when monitoring frequencies become large. Monte-Carlo simulations suffer from the same limitation, where significantly more random numbers need to be generated along a rapidly increasing number of paths as monitoring frequencies increase. Finite difference methods, in contrast, are not subject to such limitations, as long as monitoring points are placed on the computational grid, but they tend to converge slowly due to discontinuities contained in terminal conditions. Performance of finite difference methods can be improved using advanced numerical techniques such as adaptive mesh models [10], Rannacher time-stepping [11,12], and the TR-BDF2 method [13], but these methods either are formulated in a different setting (e.g. for tree-based pricing models [10]), or suffer from relatively low accuracies due to the incorporation of low-order approximations near terminal/boundary points (e.g. [11,12]).

We present an analytical modification to finite difference methods in the Black-Scholes framework to remove any discontinuities contained in the terminal conditions, which greatly improves the accuracy and efficiency of these methods. We will apply the modification to the popular Crank-Nicolson method, and introduce an effective adaptive mesh refinement technique to further improve the computational speed of the modified method. The modified finite difference method is valuable in practice thanks to its high accuracy and ease of implementation. It also complements our previous quadrature method [2] since its performance is unaffected by higher monitoring frequencies.

1.1. A Brief Review of Products Considered

The modified finite difference method is applicable to various types of discretely monitored options, including discrete barrier options, autocallable structured products, and snowball options.

Barrier options are among the most popular types of exotic options. A barrier option may be activated (knock-in option) or deactivated (knock-out option) when the underlying asset price crosses certain barrier levels. A single barrier option has one barrier at each observation date, while a double barrier option has two barriers at each observation date. The final payoff of a barrier option (if it is active at maturity) may be of the same type as that of a vanilla option or that of a digital option.

For an autocallable product, there is a pre-specified barrier level at each observation date. If the underlying asset price is greater (less) than or equal to the barrier level, the option is exercised and a pre-specified fixed-rate return is paid. If the asset price is below (above) the barriers at all observation dates, the option is never exercised and the investor receives a negative return at maturity.

A snowball option has an up-and-out barrier and a down-and-in barrier, and the down-and-in barrier is usually monitored daily. If the asset price reaches the up-and-out barrier, the option is exercised and a pre-specified fixed-rate return is paid. If the option expires without any barrier being reached, the investor receives a fixed coupon payment. If the up-and-out barrier is never reached but the down-and-in barrier is reached sometime before maturity, the investor automatically writes a vanilla put and may receive a negative return.

1.2. Organization of the Paper

The rest of the paper is organized as follows. Section 2.1 defines the class of (discrete) option pricing problems that our method is intended to solve, and the rest of Section 2 details the various ingredients of our method including, specifically, the analytical modification (Section 2.3) we introduce for standard finite difference methods (Section 2.2) in the context of discrete option pricing, and an adaptive mesh refinement technique for the analytically modified finite difference methods (Section 2.4). The efficacy and efficiency of our method is then demonstrated through extensive numerical studies in Section 3, and the paper is concluded in Section 4.

2. Materials and Methods

2.1. Basic Assumptions

We consider a general class of discretely monitored options with barriers. Note that an option with knock-in barriers is equivalent to the difference of an option without those knock-in barriers and another option with knock-out barriers at the same dates and same levels. Therefore, it suffices to consider options with only knock-out barriers. To this end, assume that

(A)

The option to be priced has two strike prices $K_m^{-},K_m^{+}\in [0,\infty ]$, with $K_m^{-}\le K_m^{+}$, at each observation date $T_m,m=1,2,\cdots ,M$. The expiration date is $T_M$.

(B)

The option is exercised if $S\le K_m^{-}$ or $S\ge K_m^{+}$ at some $T_m$, and the payoffs are given by $a_m^{-}S+b_m^{-}$ (if $S\le K_m^{-}$) and $a_m^{+}S+b_m^{+}$ (if $S\ge K_m^{+}$), respectively, for some $a_m^{\pm },b_m^{\pm }\in \mathbb{R}$.

(C)

The final payoff at maturity is

$$a_{M}S+b_M,\mathrm{for}{K}_M^{-}<S<K_M^{+},$$

for some $a_M,b_M\in \mathbb{R}$.

These assumptions are general enough to cover a wide class of discretely monitored options, such as the ones mentioned in the introduction. For instance, down-and-out put barrier options would have

$$\begin{array}{cc}\hfill 1\le m\le M-1:& K_m^{+}=\infty ,0<K_m^{-}<\infty ,a_m^{-}=0,b_m^{-}=0;\hfill \\ \hfill m=M:& K_M^{-}=0,0<K_M^{+}<\infty ,a_M=-1,b_M=K_M^{+},\hfill \\ \hfill & a_M^{+}=b_M^{+}=0.\hfill \end{array}$$

Double barrier knock-out call options would have

$$\begin{array}{cc}\hfill 1\le m\le M-1:& 0<K_m^{\pm }<\infty ,a_m^{\pm }=0,b_m^{\pm }=0;\hfill \\ \hfill m=M:& K_M^{+}=\infty ,0<K_M^{-}<\infty ,a_M=1,b_M=-K_M^{-},\hfill \\ \hfill & a_M^{-}=b_M^{-}=0.\hfill \end{array}$$

Common up-and-out autocallable products would have

$$\begin{array}{cc}\hfill 1\le m\le M:& 0<K_m^{+}<\infty ,K_m^{-}=0,a_m^{+}=0,b_m^{+}>0;\hfill \\ \hfill m=M:& a_M=0,b_M<0.\hfill \end{array}$$

We further assume that the price of the underlying asset $S\left(t\right)$ satisfies the following stochastic differential equation under the risk-neutral measure:

$$dS\left(t\right)=\left[r\left(t\right)-q\left(t\right)\right]S\left(t\right)dt+\sigma \left(t\right)S\left(t\right)dW\left(t\right),$$

(2.1)

where $r\left(t\right)$ is the risk-free forward interest rate, $q\left(t\right)$ the yield rate, $\sigma \left(t\right)$ the volatility, and $W\left(t\right)$ the Wiener process. Interest rates, yield rates, and volatilities are assumed to be deterministic functions of time.

To summarize, our basic assumptions are

• There are finitely many observation points, and two exercise levels (possibly ∞) at each observation point. If S is above the upper exercise level or below the lower exercise level at any observation point, the option is exercised and the payoff is a linear function in S.
• At maturity, if S is between the two exercise levels, a payoff is incurred which is also a linear function in S.
• The underlying asset price S follows a geometric Brownian motion with possibly time-dependent interest rates, yield rates, and volatilities.

Given the above assumptions, we shall begin the description of our method with a brief review of the finite difference methods commonly used in the pricing of discretely monitored options.

2.2. Finite Difference Methods

Suppose we wish to find an option’s value at some $T_0<T_M$ and the spot price of the underlying asset is $S_0$. Let $V(t,S)$ denote the value of the option (as a function of the asset price S) at any time $t\ge T_0$. Since S follows a geometric Brownian motion, the option’s value $V(t,S)$ satisfies the Black-Scholes PDE

$${\displaystyle \frac{\partial V}{\partial t}}(t,S)+\left[r\left(t\right)-q\left(t\right)\right]S{\displaystyle \frac{\partial V}{\partial S}}(t,S)+{\displaystyle \frac{{\sigma }^2\left(t\right)}{2}}{S}^2{\displaystyle \frac{{\partial }^{2}V}{\partial S^2}}(t,S)=r\left(t\right)V(t,S),$$

(2.2a)

for $T_0\le t\le T_M$, with the terminal condition

$$V(T_M,S)=\left\{\begin{array}{cc}{a}_{M}S+b_M,\hfill & K_M^{-}<S<K_M^{+}\hfill \\ a_M^{-}S+b_M^{-},\hfill & S\le K_M^{-}\hfill \\ a_M^{+}S+b_M^{+},\hfill & S\ge K_M^{+}\hfill \end{array}\right.,$$

(2.2b)

and boundary conditions

$$V(T_m,S)=\left\{\begin{array}{cc}{a}_m^{-}S+b_m^{-},\hfill & S\le K_m^{-}\hfill \\ a_m^{+}S+b_m^{+},\hfill & S\ge K_m^{+}\hfill \end{array}\right.,1\le m\le M-1,$$

(2.2c)

determined by Assumptions (B)–(C) from Section 2.1.

For convenience, we introduce the change of variable $Z=lnS$ to transform equation (2.2a) into a simpler form:

$${\displaystyle \frac{\partial V}{\partial t}}(t,e^Z)=\mathcal{L}\left[V(t,e^Z)\right],$$

(2.3a)

where the differential operator $\mathcal{L}$ is defined as

$$\mathcal{L}:=-{\displaystyle \frac{{\sigma }^2\left(t\right)}{2}}{\displaystyle \frac{{\partial }^2}{\partial Z^2}}-\rho \left(t\right){\displaystyle \frac{\partial }{\partial Z}}+r\left(t\right),\mathrm{with}\rho \left(t\right):=r\left(t\right)-q\left(t\right)-{\displaystyle \frac{{\sigma }^2\left(t\right)}{2}}.$$

(2.3b)

It is not convenient to solve the PDE (2.3a) directly due to the discrete boundary conditions (2.2c). Instead, for each $1\le m\le M$, we consider the PDE

$${\displaystyle \frac{\partial Y_m}{\partial t}}(t,Z)=\mathcal{L}{Y}_m(t,Z),T_{m-1}\le t\le T_m,$$

(2.4a)

supplemented with the terminal condition

$$Y_m(T_m,Z)=V(T_m,e^Z).$$

(2.4b)

Since the terminal value problem (2.4) has a unique solution for each m, we have

$$V(t,e^Z)=Y_m(t,Z),T_{m-1}<t\le T_m,$$

and by (2.2c),

$$V(T_{m-1},e^Z)=\left\{\begin{array}{cc}{Y}_m(T_{m-1},Z),\hfill & ln\left(K_{m-1}^{-}\right)<Z<ln\left(K_{m-1}^{+}\right)\hfill \\ a_{m-1}^{-}{e}^Z+b_{m-1}^{-},\hfill & Z\le ln\left(K_{m-1}^{-}\right)\hfill \\ a_{m-1}^{+}{e}^Z+b_{m-1}^{+},\hfill & Z\ge ln\left(K_{m-1}^{+}\right)\hfill \end{array}\right.,1\le m\le M.$$

(2.5)

The terminal value problems (2.4) are solved successively for $m=M,M-1,\cdots ,1$. Note that for $m=M$, the terminal condition $Y_M(T_M,Z)=V(T_M,e^Z)$ is known by (2.2b). Now for each $1\le m\le M$, assume the terminal condition $Y_m(T_m,Z)$ is known. The PDE (2.4a) is solved backward in t to yield $Y_m(T_{m-1},Z)$. Then we obtain $V(T_{m-1},e^Z)$ by (2.5), which by (2.4b) is equal to $Y_{m-1}(T_{m-1},Z)$, the terminal condition of the PDE (2.4a) for $m-1$. This process can then be repeated until we finally obtain $V(T_0,e^Z)$, from which the option’s value $V(T_0,S_0)$ can be easily deduced.

Now we give a brief review of the standard procedure of solving (2.4) numerically on a uniform grid using finite difference methods. Let $t_0,t_1,\cdots ,t_N$ denote the grid points in time, and $Z_0,Z_1,\cdots ,Z_{2P}$ the grid points in space where $Z_P=ln\left(S_0\right)$. Denote

$$\begin{array}{c}{Y}_{l,j}^m=Y_m(t_l,Z_j),r_l=r\left(t_l\right),{\rho }_l=\rho \left(t_l\right),{\sigma }_l=\sigma \left(t_l\right),\\ k={\displaystyle \frac{t_N-t_0}{N}},h={\displaystyle \frac{Z_{2P}-Z_0}{2P}},\lambda ={\displaystyle \frac{k}{h}},\mu ={\displaystyle \frac{k}{h^2}}.\end{array}$$

In practice, barriers are observed at the end of a business day, so it is convenient to find a uniform t-grid such that every monitoring point lies on the grid, i.e. for every $1\le m\le M$, there exists an $l_m$ such that $t_{l_m}=T_m$. The PDE (2.4a) can be discretized using the backward-time centered-space (BTCS), forward-time centered-space (FTCS), or Crank-Nicolson (C-N) finite difference method.

The BTCS (or explicit Euler) method consists of the recurrence equations

$$\begin{array}{cc}\hfill & Y_{l,j}^m-Y_{l-1,j}^m+\frac{1}{2}\lambda {\rho }_l(Y_{l,j+1}^m-Y_{l,j-1}^m)\hfill \\ \hfill & +\frac{1}{2}\mu {\sigma }_l^2(Y_{l,j+1}^m-2Y_{l,j}^m+Y_{l,j-1}^m)=k{r}_l{Y}_{l,j}^m,\hfill \end{array}$$

(2.6)

and is explicit (recall that we are solving (2.4a) backward in time). The FTCS (or implicit Euler) method, on the other hand, consists of the recurrence equations

$$\begin{array}{cc}\hfill & Y_{l,j}^m-Y_{l-1,j}^m+\frac{1}{2}\lambda {\rho }_{l-1}(Y_{l-1,j+1}^m-Y_{l-1,j-1}^m)\hfill \\ \hfill & +\frac{1}{2}\mu {\sigma }_{l-1}^2(Y_{l-1,j+1}^m-2Y_{l-1,j}^m+Y_{l-1,j-1}^m)=k{r}_{l-1}{Y}_{l-1,j}^m,\hfill \end{array}$$

(2.7)

and is implicit. The C-N method is a combination of the BTCS and FTCS methods and can be written as

$$\begin{array}{cc}\hfill & (\frac{1}{4}\mu {\sigma }_{l-1}^2-\frac{1}{4}\lambda {\rho }_{l-1})Y_{l-1,j-1}^m+(-1-\frac{1}{2}\mu {\sigma }_{l-1}^2-\frac{1}{2}k{r}_{l-1})Y_{l-1,j}^m\hfill \\ \hfill & +(\frac{1}{4}\mu {\sigma }_{l-1}^2+\frac{1}{4}\lambda {\rho }_{l-1})Y_{l-1,j+1}^m=-(\frac{1}{4}\mu {\sigma }_l^2-\frac{1}{4}\lambda {\rho }_l)Y_{l,j-1}^m\hfill \\ \hfill & +(-1+\frac{1}{2}\mu {\sigma }_l^2+\frac{1}{2}k{r}_l)Y_{l,j}^m-(\frac{1}{4}\mu {\sigma }_l^2+\frac{1}{4}\lambda {\rho }_l)Y_{l,j+1}^m.\hfill \end{array}$$

(2.8)

For $j\in \{0,2P\}$, we impose the extrapolation boundary condition1

$${\displaystyle \frac{{\partial }^{2}V}{\partial S^2}}(t,S)=0,$$

(2.9a)

which is equivalent to

$${\displaystyle \frac{{\partial }^2{Y}_m}{\partial Z^2}}(t,Z)-{\displaystyle \frac{\partial Y_m}{\partial Z}}(t,Z)=0,$$

(2.9b)

and can be discretized as

$$(1-\frac{1}{2}h)Y_{i,j+1}^m-2Y_{i,j}^m+(1+\frac{1}{2}h)Y_{i,j-1}^m=0,i=l-1,l,j=0,2P.$$

(2.9c)

The C-N method is generally more accurate than the BTCS and FTCS methods [14], although the FTCS method has better stability properties [13].

For convenience, we denote

$$\begin{array}{cc}\hfill {\alpha }_l=\frac{1}{4}\mu {\sigma }_l^2-\frac{1}{4}\lambda {\rho }_l,& {\beta }_l=-1-\frac{1}{2}\mu {\sigma }_l^2-\frac{1}{2}k{r}_l,\hfill \\ \hfill {\gamma }_l=\frac{1}{4}\mu {\sigma }_l^2+\frac{1}{4}\lambda {\rho }_l,& {\delta }_l=-1+\frac{1}{2}\mu {\sigma }_l^2+\frac{1}{2}k{r}_l.\hfill \end{array}$$

Equations (2.8) and (2.9c) can be written together as a matrix equation

$${\mathcal{A}}_{l-1}{\mathbf{Y}}_{l-1}^m={\mathcal{B}}_l{\mathbf{Y}}_l^m,$$

(2.10a)

where

$${\mathbf{Y}}_l^m={\langle Y_{l,0}^m,Y_{l,1}^m,\cdots ,Y_{l,2P}^m\rangle }^T,$$

is a $(2P+1)$-column vector and ${\mathcal{A}}_l,{\mathcal{B}}_l$ are $(2P+1)\times (2P+1)$ tridiagonal matrices:

$$\begin{array}{cc}\hfill {\mathcal{A}}_l& =\left(\begin{array}{cc}{\beta }_l+c_{1,+}{\alpha }_l& {\gamma }_l-c_{2,+}{\alpha }_l\\ {\alpha }_l& {\beta }_l& {\gamma }_l\\ & \ddots & \ddots & \ddots & \\ & & {\alpha }_l& {\beta }_l& {\gamma }_l\\ & & & {\alpha }_l-c_{2,-}{\gamma }_l& {\beta }_l+c_{1,-}{\gamma }_l\end{array}\right),\hfill \end{array}$$

(2.10b)

$$\begin{array}{cc}\hfill {\mathcal{B}}_l& =\left(\begin{array}{cc}{\delta }_l-c_{1,+}{\alpha }_l& -{\gamma }_l+c_{2,+}{\alpha }_l\\ -{\alpha }_l& {\delta }_l& -{\gamma }_l\\ & \ddots & \ddots & \ddots & \\ & & -{\alpha }_l& {\delta }_l& -{\gamma }_l\\ & & & -{\alpha }_l+c_{2,-}{\gamma }_l& {\delta }_l-c_{1,-}{\gamma }_l\end{array}\right),\hfill \end{array}$$

(2.10c)

where

$$c_{1,\pm }={\displaystyle \frac{2}{1\pm \frac{1}{2}h}},c_{2,\pm }={\displaystyle \frac{1\mp \frac{1}{2}h}{1\pm \frac{1}{2}h}}.$$

Note that the unknowns $Y_{i,-1}^m$ in the $j=0$ equation and $Y_{i,2P+1}^m$ in the $j=2P$ equation for $i=l-1,l$ can be eliminated using the boundary condition (2.9c). Equation (2.10a) can be solved for ${\mathbf{Y}}_{l-1}^m$ using the Thomas Algorithm, which has an $O\left(P\right)$ complexity.

Note that the terminal condition (2.4b) is typically discontinuous as can be seen from (2.5). This means that the C-N method, when applied as it is, cannot achieve a second-order convergence as for smooth terminal conditions [11].

2.3. Analytical Modification

We will remove the aforementioned discontinuities from terminal conditions by introducing an analytical modification to existing finite difference methods. First, we recall some classical results from the theory of binary options.

Lemma 2.1. 

Let $K>0$, and let ${\chi }_A$ denote the characteristic function of a set A. Consider a binary option with expiration time $T_m$ priced at $T_{m-1}$.

1. 

If the option has payoff ${\widehat{V}}_m\left(y\right)={\chi }_{[K,\infty )}y$, then ${\widehat{V}}_{m-1}\left(S\right)={\varphi }_m^a(S,K,1)$.

2. 

If ${\widehat{V}}_m\left(y\right)={\chi }_{(0,K]}y$, then ${\widehat{V}}_{m-1}\left(S\right)={\varphi }_m^a(S,K,-1)$.

3. 

If ${\widehat{V}}_m\left(y\right)={\chi }_{[K,\infty )}$, then ${\widehat{V}}_{m-1}\left(S\right)={\varphi }_m^b(S,K,1)$.

4. 

If ${\widehat{V}}_m\left(y\right)={\chi }_{(0,K]}$, then ${\widehat{V}}_{m-1}\left(S\right)={\varphi }_m^b(S,K,-1)$.

The functions ${\varphi }_m^a$ and ${\varphi }_m^b$ are defined as

$${\varphi }_m^a(S,K,ϵ)=e^{-q_m{\tau }_m}\mathcal{N}\left(ϵd_1\right)S,{\varphi }_m^b(S,K,ϵ)=e^{-r_m{\tau }_m}\mathcal{N}\left(ϵd_2\right),$$

where ${\tau }_m=T_m-T_{m-1}$,

$$r_m={\int }_{T_{m-1}}^{T_m}{\displaystyle \frac{r\left(t\right)}{{\tau }_m}}dt,q_m={\int }_{T_{m-1}}^{T_m}{\displaystyle \frac{q\left(t\right)}{{\tau }_m}}dt,{\sigma }_m^2={\int }_{T_{m-1}}^{T_m}{\displaystyle \frac{{\sigma }^2\left(t\right)}{{\tau }_m}}dt,$$

$\mathcal{N}$ is the cumulative normal distribution function, and

$$d_1={\displaystyle \frac{1}{{\sigma }_m\sqrt{{\tau }_m}}}\left\{log{\displaystyle \frac{S}{K}}+\left(r_m-q_m+{\displaystyle \frac{{\sigma }_m^2}{2}}\right){\tau }_m\right\},d_2=d_1-{\sigma }_m\sqrt{{\tau }_m}.$$

Proof. 

By definition, ${\varphi }_m^a$ is the value of an asset-or-nothing option, and ${\varphi }_m^b$ is the value of a cash-or-nothing option. The valuation formulas are just standard results for binary options [14]. □

Since $Y_m(T_m,Z)$ in the terminal condition (2.4b) generally has discontinuities at $ln\left(K_m^{\pm }\right)$, we do not solve the terminal value problem (2.4) directly, but instead seek to remove these discontinuities first. We describe our basic idea as follows.

2.3.1. Basic idea of the analytical modification

For each $1\le m\le M$, note that $V(T_m,S)$ is continuous for $K_m^{-}<S<K_m^{+}$. We define the one-sided limits

$$\begin{array}{cc}\hfill {\tilde{V}}_m^{+}=\underset{ϵ\to 0^{+}}{lim}V(T_m,K_m^{+}-ϵ),& {\tilde{V}}_m^{-}=\underset{ϵ\to 0^{+}}{lim}V(T_m,K_m^{-}+ϵ),\hfill \end{array}$$

(2.11a)

$$\begin{array}{cc}\hfill {\tilde{V}}_{S,m}^{+}=\underset{ϵ\to 0^{+}}{lim}{\displaystyle \frac{\partial V}{\partial S}}(T_m,K_m^{+}-ϵ),& {\tilde{V}}_{S,m}^{-}=\underset{ϵ\to 0^{+}}{lim}{\displaystyle \frac{\partial V}{\partial S}}(T_m,K_m^{-}+ϵ),\hfill \end{array}$$

(2.11b)

and consider $U_m(t,Z)$ defined as the solution to the PDE

$${\displaystyle \frac{\partial U_m}{\partial t}}(t,Z)=\mathcal{L}{U}_m(t,Z),T_{m-1}\le t\le T_m,$$

(2.12a)

supplemented with the (continuously differentiable) terminal condition

$$U_m(T_m,Z)=\left\{\begin{array}{cc}V(T_m,e^Z),\hfill & ln\left(K_m^{-}\right)<Z<ln\left(K_m^{+}\right)\hfill \\ {\tilde{V}}_{S,m}^{-}(e^Z-K_m^{-})+{\tilde{V}}_m^{-},\hfill & Z\le ln\left(K_m^{-}\right)\hfill \\ {\tilde{V}}_{S,m}^{+}(e^Z-K_m^{+})+{\tilde{V}}_m^{+},\hfill & Z\ge ln\left(K_m^{+}\right)\hfill \end{array}\right..$$

(2.12b)

Note that by (2.2b)-(2.2c), we have

$$V(T_m,e^Z)-U_m(T_m,Z)=\left\{\begin{array}{cc}0,\hfill & ln\left(K_m^{-}\right)<Z<ln\left(K_m^{+}\right)\hfill \\ {\tilde{a}}_m^{-}{e}^Z+{\tilde{b}}_m^{-},\hfill & Z\le ln\left(K_m^{-}\right)\hfill \\ {\tilde{a}}_m^{+}{e}^Z+{\tilde{b}}_m^{+},\hfill & Z\ge ln\left(K_m^{+}\right)\hfill \end{array}\right.,$$

(2.13)

where

$${\tilde{a}}_m^{\pm }=a_m^{\pm }-{\tilde{V}}_{S,m}^{\pm },{\tilde{b}}_m^{\pm }=b_m^{\pm }+{\tilde{V}}_{S,m}^{\pm }{K}_m^{\pm }-{\tilde{V}}_m^{\pm }.$$

This shows that $V(t,e^Z)-U_m(t,Z)$ is the sum of values of four binary options for $T_{m-1}\le t\le T_m$, and hence by Lemma 2.1,

$$\begin{array}{cc}\hfill & V(T_{m-1},e^Z)-U_m(T_{m-1},Z)={\tilde{a}}_m^{-}{\varphi }_m^a(e^Z,K_m^{-},-1)+{\tilde{b}}_m^{-}{\varphi }_m^b(e^Z,K_m^{-},-1)\hfill \\ \hfill & +{\tilde{a}}_m^{+}{\varphi }_m^a(e^Z,K_m^{+},1)+{\tilde{b}}_m^{+}{\varphi }_m^b(e^Z,K_m^{+},1),ln\left(K_{m-1}^{-}\right)<Z<ln\left(K_{m-1}^{+}\right),\hfill \end{array}$$

(2.14a)

or

$$\begin{array}{cc}\hfill & V(T_{m-1},e^Z)=U_m(T_{m-1},Z)+{\tilde{a}}_m^{-}{\varphi }_m^a(e^Z,K_m^{-},-1)+{\tilde{b}}_m^{-}{\varphi }_m^b(e^Z,K_m^{-},-1)\hfill \\ \hfill & +{\tilde{a}}_m^{+}{\varphi }_m^a(e^Z,K_m^{+},1)+{\tilde{b}}_m^{+}{\varphi }_m^b(e^Z,K_m^{+},1),ln\left(K_{m-1}^{-}\right)<Z<ln\left(K_{m-1}^{+}\right).\hfill \end{array}$$

(2.14b)

For each $1\le m\le M$, we solve (2.12) numerically for $U_m(T_{m-1},Z)$ and obtain $V(T_{m-1},e^Z)$ using (2.14b). Repeating this process for $m=M,M-1,\cdots ,1$, we finally obtain $V(T_0,e^Z)$, from which the option’s value $V(T_0,S_0)$ can be easily deduced.

Remark 2.2.

For the given terminal condition (2.2b), a closed-form expression can in fact be obtained for $U_M(t,Z)$; see, for example, [2]. In practice, however, M is often “reasonably large”, e.g., in the order of tens or hundreds, so that the incorporation of such exact formulas of $U_M$ makes little difference in the numerical approximations obtained for $V(T_0,S)$. Thus we do not specifically differentiate $U_M$ from other $U_m$’s ($m\le M-1$) in our study.

2.3.2. Implementation of the Analytical Modification

In practice, it may not be possible to place all barriers on the computational grid, which means we may not be able to compute the values of ${\tilde{V}}_m^{\pm },{\tilde{V}}_{S,m}^{\pm }$ for all $1\le m\le M-1$ as defined in (2.11). Note that the values of ${\tilde{V}}_M^{\pm },{\tilde{V}}_{S,M}^{\pm }$ for $m=M$ follow from the closed-form expression of $V(T_M,e^Z)$ as given in (2.2b). For $1\le m\le M-1$, however, we will have to find an approximation to these values. To this end, denote

$${\tilde{U}}_M(t,Z)=U_M(t,Z),$$

(2.15)

which is the smooth solution defined in (2.12) for $m=M$. For each $1\le m\le M-1$, assume we have already computed a three times continuously differentiable ($C^3$) function ${\tilde{U}}_{m+1}(T_m,Z)$ for $Z\in \{Z_0,Z_1,\cdots ,Z_{2P}\}$ and have obtained the smooth part of $V(T_m,e^Z)$, namely $Y_{m+1}(T_m,Z)$ (see (2.5)), at the grid points. Note that this is true for $m=M-1$, in which case we have already computed a $C^3$ function ${\tilde{U}}_M(T_{M-1},Z)$ by (2.15) and have obtained $Y_M(T_{M-1},Z)$ by (2.14b).

Denote

$$Y_{l,j}^{m+1}=Y_{m+1}(t_l,Z_j),Y_{S,l,j}^{m+1}={\displaystyle \frac{Y_{l,j+1}^{m+1}-Y_{l,j-1}^{m+1}}{2h{e}^{Z_j}}}\approx {\displaystyle \frac{\partial Y_{m+1}}{\partial S}}(t_l,Z_j),$$

and let

$$p_m^{-}=min\left\{j>0:Z_j>ln\left(K_m^{-}\right)\right\},p_m^{+}=max\left\{j<2P:Z_j<ln\left(K_m^{+}\right)\right\}.$$

Note that $ln\left(K_m^{\pm }\right)$ are located between $Z_{p_m^{\pm }}$ and $Z_{p_m^{\pm }\pm 1}$. We use linear interpolations of $Y_{m+1}(T_m,Z)$ on the grid to approximate ${\tilde{V}}_m^{\pm },{\tilde{V}}_{S,m}^{\pm }$ as follows:

$$\begin{array}{cc}\hfill c_m^{-}& =Y_{S,l_m,p_m^{-}}^{m+1}-{\displaystyle \frac{1}{h}}(Y_{S,l_m,p_m^{-}-1}^{m+1}-Y_{S,l_m,p_m^{-}}^{m+1})\left[ln\left(K_m^{-}\right)-Z_{p_m^{-}}\right],\hfill \end{array}$$

(2.16a)

$$\begin{array}{cc}\hfill c_m^{+}& =Y_{S,l_m,p_m^{+}}^{m+1}+{\displaystyle \frac{1}{h}}(Y_{S,l_m,p_m^{+}+1}^{m+1}-Y_{S,l_m,p_m^{+}}^{m+1})\left[ln\left(K_m^{+}\right)-Z_{p_m^{+}}\right],\hfill \end{array}$$

(2.16b)

$$\begin{array}{cc}\hfill d_m^{-}& =Y_{l_m,p_m^{-}}^{m+1}-{\displaystyle \frac{1}{h}}(Y_{l_m,p_m^{-}-1}^{m+1}-Y_{l_m,p_m^{-}}^{m+1})\left[ln\left(K_m^{-}\right)-Z_{p_m^{-}}\right],\hfill \end{array}$$

(2.16c)

$$\begin{array}{cc}\hfill d_m^{+}& =Y_{l_m,p_m^{+}}^{m+1}+{\displaystyle \frac{1}{h}}(Y_{l_m,p_m^{+}+1}^{m+1}-Y_{l_m,p_m^{+}}^{m+1})\left[ln\left(K_m^{+}\right)-Z_{p_m^{+}}\right].\hfill \end{array}$$

(2.16d)

Since ${\tilde{U}}_{m+1}(T_m,Z)$ and hence $Y_{m+1}(T_m,Z)$ is $C^3$ in Z, we have

$$c_m^{\pm }-{\tilde{V}}_{S,m}^{\pm }=O\left(h^2\right),d_m^{\pm }-{\tilde{V}}_m^{\pm }=O\left(h^2\right).$$

(2.17)

This means that the order of the interpolation error does not exceed that of the C-N method.

Consider ${\tilde{U}}_m(t,Z)$ defined as the solution to the PDE

$${\displaystyle \frac{\partial {\tilde{U}}_m}{\partial t}}(t,Z)=\mathcal{L}{\tilde{U}}_m(t,Z),T_{m-1}\le t\le T_m,$$

(2.18a)

supplemented with the terminal condition

$${\tilde{U}}_m(T_m,Z)=\left\{\begin{array}{cc}V(T_m,e^Z),\hfill & ln\left(K_m^{-}\right)<Z<ln\left(K_m^{+}\right)\hfill \\ c_m^{-}(e^Z-K_m^{-})+d_m^{-},\hfill & Z\le ln\left(K_m^{-}\right)\hfill \\ c_m^{+}(e^Z-K_m^{+})+d_m^{+},\hfill & Z\ge ln\left(K_m^{+}\right)\hfill \end{array}\right..$$

(2.18b)

Note that by (2.2c), we have

$$V(T_m,e^Z)-{\tilde{U}}_m(T_m,Z)=\left\{\begin{array}{cc}0,\hfill & ln\left(K_m^{-}\right)<Z<ln\left(K_m^{+}\right)\hfill \\ {\tilde{c}}_m^{-}{e}^Z+{\tilde{d}}_m^{-},\hfill & Z\le ln\left(K_m^{-}\right)\hfill \\ {\tilde{c}}_m^{+}{e}^Z+{\tilde{d}}_m^{+},\hfill & Z\ge ln\left(K_m^{+}\right)\hfill \end{array}\right.,$$

(2.19)

where

$${\tilde{c}}_m^{\pm }=a_m^{\pm }-c_m^{\pm },{\tilde{d}}_m^{\pm }=b_m^{\pm }+c_m^{\pm }{K}_m^{\pm }-d_m^{\pm }.$$

This shows that $V(t,e^Z)-{\tilde{U}}_m(t,Z)$ is the sum of values of four binary options for $T_{m-1}\le t\le T_m$, and hence by Lemma 2.1,

$$\begin{array}{cc}\hfill & V(T_{m-1},e^Z)-{\tilde{U}}_m(T_{m-1},Z)={\tilde{c}}_m^{-}{\varphi }_m^a(e^Z,K_m^{-},-1)+{\tilde{d}}_m^{-}{\varphi }_m^b(e^Z,K_m^{-},-1)\hfill \\ \hfill & +{\tilde{c}}_m^{+}{\varphi }_m^a(e^Z,K_m^{+},1)+{\tilde{d}}_m^{+}{\varphi }_m^b(e^Z,K_m^{+},1),ln\left(K_{m-1}^{-}\right)<Z<ln\left(K_{m-1}^{+}\right),\hfill \end{array}$$

(2.20a)

or

$$\begin{array}{cc}\hfill & V(T_{m-1},e^Z)={\tilde{U}}_m(T_{m-1},Z)+{\tilde{c}}_m^{-}{\varphi }_m^a(e^Z,K_m^{-},-1)+{\tilde{d}}_m^{-}{\varphi }_m^b(e^Z,K_m^{-},-1)\hfill \\ \hfill & +{\tilde{c}}_m^{+}{\varphi }_m^a(e^Z,K_m^{+},1)+{\tilde{d}}_m^{+}{\varphi }_m^b(e^Z,K_m^{+},1),ln\left(K_{m-1}^{-}\right)<Z<ln\left(K_{m-1}^{+}\right).\hfill \end{array}$$

(2.20b)

Now that we have computed ${\tilde{U}}_m(T_{m-1},Z)$ and the smooth part of $V(T_{m-1},e^Z)$, namely $Y_m(T_{m-1},Z)$, we continue to compute $c_{m-1}^{\pm },d_{m-1}^{\pm }$ as in (2.16), which means we have the terminal condition ${\tilde{U}}_{m-1}(T_{m-1},Z)$ as in (2.18b). Then we repeat this process for $m-1$.

We summarize the modified finite difference method as follows:

1:

for$m=M-1$ downto 1 do

2:

   Compute $c_m^{\pm },d_m^{\pm }$ using (2.16)

3:

   Solve (2.18) numerically for ${\tilde{U}}_m(T_{m-1},Z)$

4:

   Compute $Y_m(T_{m-1},Z)$, the smooth part of $V(T_{m-1},e^Z)$, using (2.20b)

5:

end for

6:

Obtain option value $V(T_0,S_0)=V(T_0,e^{Z_P})$

Note that the terminal value problem (2.18) is the same as (2.4) except that the terminal condition (2.18b) now has a higher degree of smoothness, and the problem (2.18) can be solved numerically using the same finite difference methods as the ones described in Section 2.2, in particular the C-N method. Since the terminal condition ${\tilde{U}}_m(T_m,Z)$ in (2.18b) is almost continuously differentiable (the jumps are of order $h^2$ by (2.17)), the C-N method converges quite well in practice (Section 3.1). In principle, Rannacher time-stepping can be used to improve the performance of the C-N method near monitoring points, but this is important only when accurate discrete approximations to option Greeks are needed (Section 3.4).

2.4. Adaptive Mesh Refinement

Note that the function ${\tilde{U}}_m(t,Z)$ in (2.18) tends to have large derivatives near the barriers and behaves relatively smoothly away from the barriers. This means that it makes sense to use finer grids near the barriers and coarser ones elsewhere.

Adaptive mesh refinement techniques are well-known and are widely used in numerical solutions of PDEs. In the field of quantitative finance, adaptive mesh models have been developed for the BTCS method [10]. Alternatively, nonuniform grids based on changes of variables have also been introduced in the literature [13]. Since here we are using an analytically modified C-N method, we introduce an effective adaptive mesh refinement technique which further improves the computational efficiency of our method.

We begin by assuming that the upper barriers $K_m^{+}$ (and lower barriers $K_m^{-}$, respectively) are equal or close to each other for all $1\le m\le M$, which is almost always satisfied in practice. Our basic idea is to construct fine grids near the barriers, and solve the Black-Scholes PDEs on both the original and the fine grids at each time step.

First, we find integers $0<q_1^{-}<q_2^{-}<q_1^{+}<q_2^{+}<2P$ such that

$$exp\left(Z_{q_1^{-}}\right)<K_m^{-}<exp\left(Z_{q_2^{-}}\right),\mathrm{and}exp\left(Z_{q_1^{+}}\right)<K_m^{+}<exp\left(Z_{q_2^{+}}\right),$$

for all $1\le m\le M$. Then we construct a fine grid for $[q_1^{-},q_2^{-}]$ and a similar one for $[q_1^{+},q_2^{+}]$ by dividing each interval of size h into n equal parts. We illustrate the fine grid near $K_m^{+}$ for $n=4$ in Figure 1.

Let

$$n^{\pm }=n(q_2^{\pm }-q_1^{\pm }),x_j^{\pm }=Z_{q_1^{\pm }}+jh/n,0\le j\le n^{\pm },$$

and define

$${\widehat{U}}_M(t,Z)=U_M(t,Z),$$

(2.21)

which is the smooth solution defined in (2.12) for $m=M$. For each $1\le m\le M-1$, assume we have already computed a $C^3$ function ${\widehat{U}}_{m+1}(T_m,Z)$ for

$$Z\in {\mathcal{G}}_h:=\left\{Z_0,Z_1,\cdots ,Z_{2P}\right\}\bigcup \left\{x_0^{\pm },x_1^{\pm },\cdots ,x_{n^{\pm }}^{\pm }\right\},$$

and have obtained the smooth part of $V(T_m,e^Z)$, namely $Y_{m+1}(T_m,Z)$, at the grid points. Note that this is true for $m=M-1$ by (2.21) and (2.14b).

Recall from (2.16) that $c_m^{\pm },d_m^{\pm }$ are approximations to ${\tilde{V}}_{S,m}^{\pm },{\tilde{V}}_m^{\pm }$ computed using the values of $Y_{m+1}(T_m,Z)$ on the original (coarse) grid. We shall now find approximations to ${\tilde{V}}_{S,m}^{\pm },{\tilde{V}}_m^{\pm }$ using the values of $Y_{m+1}(T_m,Z)$ on the fine grids. To this end, denote

$$\begin{array}{cc}\hfill Y_{l,j}^{\pm (m+1)}& =Y_{m+1}(t_l,x_j^{\pm }),\hfill \\ \hfill Y_{S,l,j}^{\pm (m+1)}& ={\displaystyle \frac{Y_{l,j+1}^{\pm (m+1)}-Y_{l,j-1}^{\pm (m+1)}}{2h{e}^{x_j^{\pm }}/n}}\approx {\displaystyle \frac{\partial Y_{m+1}}{\partial S}}(t_l,x_j^{\pm }),\hfill \end{array}$$

and let

$$u_m^{-}=min\left\{j>0:x_j^{-}>ln\left(K_m^{-}\right)\right\},u_m^{+}=max\left\{j<n^{+}:x_j^{+}<ln\left(K_m^{+}\right)\right\}.$$

Define

$$\begin{array}{cc}\hfill f_m^{-}& =Y_{S,l_m,u_m^{-}}^{-(m+1)}-{\displaystyle \frac{n}{h}}\left[Y_{S,l_m,u_m^{-}-1}^{-(m+1)}-Y_{S,l_m,u_m^{-}}^{-(m+1)}\right]\left[ln\left(K_m^{-}\right)-x_{u_m^{-}}^{-}\right],\hfill \end{array}$$

(2.22a)

$$\begin{array}{cc}\hfill f_m^{+}& =Y_{S,l_m,u_m^{+}}^{+(m+1)}+{\displaystyle \frac{n}{h}}\left[Y_{S,l_m,u_m^{+}+1}^{+(m+1)}-Y_{S,l_m,u_m^{+}}^{+(m+1)}\right]\left[ln\left(K_m^{+}\right)-x_{u_m^{+}}^{+}\right],\hfill \end{array}$$

(2.22b)

$$\begin{array}{cc}\hfill g_m^{-}& =Y_{l_m,u_m^{-}}^{-(m+1)}-{\displaystyle \frac{n}{h}}\left[Y_{l_m,u_m^{-}-1}^{-(m+1)}-Y_{l_m,u_m^{-}}^{-(m+1)}\right]\left[ln\left(K_m^{-}\right)-x_{u_m^{-}}^{-}\right],\hfill \end{array}$$

(2.22c)

$$\begin{array}{cc}\hfill g_m^{+}& =Y_{l_m,u_m^{+}}^{+(m+1)}+{\displaystyle \frac{n}{h}}\left[Y_{l_m,u_m^{+}+1}^{+(m+1)}-Y_{l_m,u_m^{+}}^{+(m+1)}\right]\left[ln\left(K_m^{+}\right)-x_{u_m^{+}}^{+}\right].\hfill \end{array}$$

(2.22d)

Note that, in general, $f_m^{\pm },g_m^{\pm }$ provide a better approximation to ${\tilde{V}}_{S,m}^{\pm },{\tilde{V}}_m^{\pm }$ than $c_m^{\pm },d_m^{\pm }$.

Consider ${\widehat{U}}_m(t,Z)$ defined as the solution to the PDE

$${\displaystyle \frac{\partial {\widehat{U}}_m}{\partial t}}(t,Z)=\mathcal{L}{\widehat{U}}_m(t,Z),T_{m-1}\le t\le T_m,$$

(2.23a)

supplemented with the terminal condition

$${\widehat{U}}_m(T_m,Z)=\left\{\begin{array}{cc}V(T_m,e^Z),\hfill & ln\left(K_m^{-}\right)<Z<ln\left(K_m^{+}\right)\hfill \\ f_m^{-}(e^Z-K_m^{-})+g_m^{-},\hfill & Z\le ln\left(K_m^{-}\right)\hfill \\ f_m^{+}(e^Z-K_m^{+})+g_m^{+},\hfill & Z\ge ln\left(K_m^{+}\right)\hfill \end{array}\right..$$

(2.23b)

The PDE (2.23a) is the same as (2.18a), but we do not solve it all the way back to $T_{m-1}$. Instead, for each $l_{m-1}<l\le l_m$, assume we have computed ${\widehat{U}}_m(t_l,Z)$ for $Z\in {\mathcal{G}}_h$, which is true for $l=l_m$ by (2.23b). We first solve the PDE (2.23a) numerically using the terminal condition ${\widehat{U}}_m(t_l,Z_k)$ on the coarse grid for just one time step from $t_l$ to $t_{l-1}$. More specifically, denote by ${\widehat{U}}_{l,k}^m={\widehat{U}}_m(t_l,Z_k)$ the solution computed on the coarse grid; we solve the matrix equation (2.10a) which in our new notation becomes

$${\mathcal{A}}_{l-1}{\widehat{\mathbf{U}}}_{l-1}^m={\mathcal{B}}_l{\widehat{\mathbf{U}}}_l^m,$$

(2.24)

where

$${\widehat{\mathbf{U}}}_l^m={\langle {\widehat{U}}_{l,0}^m,{\widehat{U}}_{l,1}^m,\cdots ,{\widehat{U}}_{l,2P}^m\rangle }^T,$$

is a $(2P+1)$-column vector and ${\mathcal{A}}_l,{\mathcal{B}}_l$ are $(2P+1)\times (2P+1)$ tridiagonal matrices as defined in (2.10). After solving for ${\widehat{\mathbf{U}}}_{l-1}^m$ from (2.24), we solve the PDE (2.23a) again on the fine grids for n time steps with reduced step size $k/n$, using the boundary conditions linearly interpolated (in time) from ${\widehat{\mathbf{U}}}_{l-1}^m,{\widehat{\mathbf{U}}}_l^m$ at $q_1^{\pm }$ and $q_2^{\pm }$. For convenience, we denote, for $0\le i\le n$,

$$\begin{array}{c}{r}_{l,i}=r(t_l+ik/n),{\rho }_{l,i}=\rho (t_l+ik/n),{\sigma }_{l,i}=\sigma (t_l+ik/n),\\ {\alpha }_{n,l,i}=\frac{1}{4}n\mu {\sigma }_{l,i}^2-\frac{1}{4}\lambda {\rho }_{l,i},{\beta }_{n,l,i}=-1-\frac{1}{2}n\mu {\sigma }_{l,i}^2-\frac{1}{2n}k{r}_{l,i},\\ {\gamma }_{n,l,i}=\frac{1}{4}n\mu {\sigma }_{l,i}^2+\frac{1}{4}\lambda {\rho }_{l,i},{\delta }_{n,l,i}=-1+\frac{1}{2}n\mu {\sigma }_{l,i}^2+\frac{1}{2n}k{r}_{l,i},\end{array}$$

and $W_{l,i,j}^{\pm m}={\widehat{U}}_m(t_{l,i},x_j^{\pm })$. The matrix equation on the fine grids then becomes

$${\mathcal{A}}_{n,l-1,i-1}{\mathbf{W}}_{l-1,i-1}^{\pm m}={\mathcal{B}}_{n,l-1,i}{\mathbf{W}}_{l-1,i}^{\pm m}+{\mathbf{w}}_{l-1,i}^{\pm },n\ge i\ge 1,$$

(2.25a)

where

$$\begin{array}{c}{\mathbf{W}}_{l-1,i}^{\pm m}={\langle W_{l-1,i,1}^{\pm m},W_{l-1,i,2}^{\pm m},\cdots ,W_{l-1,i,n^{\pm }-1}^{\pm m}\rangle }^T,\\ {\mathbf{w}}_{l-1,i}^{\pm }=\left(\begin{array}{c}-{\alpha }_{n,l-1,i}{W}_{l-1,i,0}^{\pm m}-{\alpha }_{n,l-1,i-1}{W}_{l-1,i-1,0}^{\pm m}\\ 0\\ ⋮\\ 0\\ -{\gamma }_{n,l-1,i}{W}_{l-1,i,n^{\pm }}^{\pm m}-{\gamma }_{n,l-1,i-1}{W}_{l-1,i-1,n^{\pm }}^{\pm m}\end{array}\right),\end{array}$$

(2.25b)

are $(n^{\pm }-1)$-column vectors and ${\mathcal{A}}_{n,l-1,i},{\mathcal{B}}_{n,l-1,i}$ are $(n^{\pm }-1)\times (n^{\pm }-1)$ tridiagonal Toeplitz matrices:

$$\begin{array}{cc}\hfill {\mathcal{A}}_{n,l-1,i}& =\left(\begin{array}{cc}{\beta }_{n,l-1,i}& {\gamma }_{n,l-1,i}\\ {\alpha }_{n,l-1,i}& {\beta }_{n,l-1,i}& {\gamma }_{n,l-1,i}\\ & \ddots & \ddots & \ddots & \\ & & {\alpha }_{n,l-1,i}& {\beta }_{n,l-1,i}& {\gamma }_{n,l-1,i}\\ & & & {\alpha }_{n,l-1,i}& {\beta }_{n,l-1,i}\end{array}\right),\hfill \end{array}$$

(2.25c)

$$\begin{array}{cc}\hfill {\mathcal{B}}_{n,l-1,i}& =\left(\begin{array}{cc}{\delta }_{n,l-1,i}& -{\gamma }_{n,l-1,i}\\ -{\alpha }_{n,l-1,i}& {\delta }_{n,l-1,i}& -{\gamma }_{n,l-1,i}\\ & \ddots & \ddots & \ddots & \\ & & -{\alpha }_{n,l-1,i}& {\delta }_{n,l-1,i}& -{\gamma }_{n,l-1,i}\\ & & & -{\alpha }_{n,l-1,i}& {\delta }_{n,l-1,i}\end{array}\right).\hfill \end{array}$$

(2.25d)

Note that, in (2.25b),

$$\begin{array}{cc}\hfill W_{l-1,i,0}^{\pm m}& ={\widehat{U}}_{l-1,q_1^{\pm }}^m+{\displaystyle \frac{i}{n}}({\widehat{U}}_{l,q_1^{\pm }}^m-{\widehat{U}}_{l-1,q_1^{\pm }}^m),\hfill \end{array}$$

(2.26a)

$$\begin{array}{cc}\hfill W_{l-1,i,n^{\pm }}^{\pm m}& ={\widehat{U}}_{l-1,q_2^{\pm }}^m+{\displaystyle \frac{i}{n}}({\widehat{U}}_{l,q_2^{\pm }}^m-{\widehat{U}}_{l-1,q_2^{\pm }}^m),\hfill \end{array}$$

(2.26b)

are boundary conditions linearly interpolated (in time) from coarse grid solutions ${\widehat{\mathbf{U}}}_{l-1}^m,{\widehat{\mathbf{U}}}_l^m$ at $q_1^{\pm }$ and $q_2^{\pm }$. After solving for ${\mathbf{W}}_{l-1,0}^{\pm m}$ from (2.25a), we set

$$\begin{array}{cc}\hfill {\widehat{U}}_m(t_{l-1},Z_j)={\widehat{U}}_{l-1,j}^m,& 0\le j\le 2P,j\notin \left\{q_1^{\pm }+1,q_1^{\pm }+2,\cdots ,q_2^{\pm }-1\right\},\hfill \end{array}$$

(2.27a)

and

$$\begin{array}{cc}\hfill 1-1{\widehat{U}}_m(t_{l-1},x_j^{\pm })=W_{l-1,0,j}^{\pm m},& 1\le j\le n^{\pm }-1.\hfill \end{array}$$

(2.27b)

Since $Z_{q_1^{\pm }+j}=x_{nj}^{\pm }$ for $0\le j\le q_2^{\pm }-q_1^{\pm }$, equation (2.27) completes the definition of ${\widehat{U}}_m(t_{l-1},Z)$ for $Z\in {\mathcal{G}}_h$. We then repeat this process for $l=l_m,l_m-1,\cdots ,l_{m-1}+1$ to obtain ${\widehat{U}}_m(T_{m-1},Z)$ for $Z\in {\mathcal{G}}_h$.

By (2.23b), $V(T_m,e^Z)-{\widehat{U}}_m(T_m,Z)$ is piecewise linear in $e^Z$, meaning $V(t,e^Z)-{\widehat{U}}_m(t,Z)$ is the sum of values of four binary options. It follows from Lemma 2.1 that

$$\begin{array}{cc}\hfill & V(T_{m-1},e^Z)={\widehat{U}}_m(T_{m-1},Z)+{\tilde{f}}_m^{-}{\varphi }_m^a(e^Z,K_m^{-},-1)+{\tilde{g}}_m^{-}{\varphi }_m^b(e^Z,K_m^{-},-1)\hfill \\ \hfill & +{\tilde{f}}_m^{+}{\varphi }_m^a(e^Z,K_m^{+},1)+{\tilde{g}}_m^{+}{\varphi }_m^b(e^Z,K_m^{+},1),ln\left(K_{m-1}^{-}\right)<Z<ln\left(K_{m-1}^{+}\right),\hfill \end{array}$$

(2.28)

where

$${\tilde{f}}_m^{\pm }=a_m^{\pm }-f_m^{\pm },{\tilde{g}}_m^{\pm }=b_m^{\pm }+f_m^{\pm }{K}_m^{\pm }-g_m^{\pm }.$$

Then we repeat this process for $m-1$.

We summarize the adaptive mesh refinement procedure as follows:

1:

Construct a fine grid for $[q_1^{-},q_2^{-}]$ and a similar one for $[q_1^{+},q_2^{+}]$ by dividing each interval of size h into n equal parts

2:

for$m=M-1$ downto 1 do

3:

    Compute $f_m^{\pm },g_m^{\pm }$ using (2.22)

4:

    Solve (2.23) numerically for one time step (with step size k) on coarse grid

5:

    Solve (2.23) numerically for n time steps (with step size $k/n$) on fine grids, using (2.26) as boundary conditions

6:

    Combine numerical solutions on coarse and fine grids to obtain ${\widehat{U}}_m(T_{m-1},Z)$

7:

    Compute $Y_m(T_{m-1},Z)$, the smooth part of $V(T_{m-1},e^Z)$, using (2.28)

8:

end for

9:

Obtain option value $V(T_0,S_0)=V(T_0,e^{Z_P})$

3. Results

3.1. Convergence Study

To study the convergence properties of the proposed method, either with or without mesh adaptation, we first apply the original and the analytically modified C-N method to the following five options:

• Option I: A cash-or-nothing call option with final payoff:

  $$V(T_M,S)=\left\{\begin{array}{cc}0,\hfill & S<K_M^{+}\hfill \\ b_M^{+},\hfill & S\ge K_M^{+}\hfill \end{array}\right..$$
• Option II: An asset-or-nothing call option with final payoff:

  $$V(T_M,S)=\left\{\begin{array}{cc}0,\hfill & S<K_M^{+}\hfill \\ a_M^{+}S,\hfill & S\ge K_M^{+}\hfill \end{array}\right..$$
• Option III: A vanilla call option with final payoff:

  $$V(T_M,S)=\left\{\begin{array}{cc}0,\hfill & S<K_M^{+}\hfill \\ a_M^{+}(S-K_M^{+}),\hfill & S\ge K_M^{+}\hfill \end{array}\right..$$
• Option IV: An exotic call option with final payoff:

  $$V(T_M,S)=\left\{\begin{array}{cc}0,\hfill & S<K_M^{+}\hfill \\ a_M^{+}{(S-K_M^{+})}^2,\hfill & S\ge K_M^{+}\hfill \end{array}\right..$$
  We remark that this option does not have a counterpart in real-life financial applications. It is designed here solely to test the convergence properties of the proposed method.

These four options are standard options without discretely monitored structures, and for all cases we set

$$\begin{array}{c}r=0.02,q=0,\sigma =0.2,\\ T_0=0,T_M=1,S_0=1.1,K_M^{+}=1.2,a_M^{+}=b_M^{+}=100.\end{array}$$

Note also that the final payoffs of the above four options have increasingly higher degree of smoothness (discontinuous, discontinuous, $C^0$, and $C^1$, respectively). The fifth option, on the other hand, is a discretely monitored option:

• Option V: A double-barrier knock-out option with different monitoring frequencies for the two barriers. Assume there are 250 business days in a year, and the option will expire one year from the valuation date which we denote as day 0. The spot price of the underlying asset at day 0 is 1.1. The upper barrier level is 1.2, and its observation dates are {15, 36, 57, 78, 99, 120, 141, 162, 183, 204, 225, 246}, resembling a monthly monitoring structure. The lower barrier is monitored daily, and the barrier level is 0.9. We further assume that the notional value of the option is 100, that the payoff at maturity is 10% of the notional value, namely 10, and that

  $$r=0.02,q=0,\sigma =0.2.$$

In practice, this type of barrier structure is common for snowball options and their variants. Note that the lower barrier of a snowball option is a knock-in barrier, but the option can be priced as the difference between an up-and-out barrier option and a double-barrier knock-out option using PDE methods. Table 1, Table 2, Table 3, Table 4, Table 5, Table 6, Table 7 and Table 8 present results from convergence studies of the four standard options (I-IV) priced using the original C-N method. Since our purpose here is to examine numerically the relationship between convergence rates and smoothness of terminal conditions, the analytically modified C-N method is not used which otherwise would produce (almost) exact pricing results for the first three options (I-III). All four options are priced on the interval2

$$I_0:=[S_0{e}^{-6\sigma },S_0{e}^{6\sigma }]\approx [0.33131363,3.65212862],$$

with a time step $k=(T_M-T_0)/N=1/N$ and a mesh size $h=k$, where N ranges from 10 to 5,120. For the first three options, exact option values are available from analytical formulas, while for the fourth option, the option value is estimated on a finer grid with $N=10,240$. For each numerical approximation $v_{0,j}$ to $V_{0,j}:=V(T_0,S_0{e}^{(j-P)h})$ calculated at $T_0$ on the computational grid, two types of errors are defined:

$$\begin{array}{cc}\hfill L^2-\mathrm{error}& ={\left(\sum _{j=0}^{2P}{\left|v_{0,j}-V_{0,j}\right|}^{2}h\right)}^{1/2},\hfill \\ \hfill L^{\infty }-\mathrm{error}& =\underset{j=0,1,\cdots ,2P}{max}\left|v_{0,j}-V_{0,j}\right|.\hfill \end{array}$$

The numerical order of an error ($L^2$ or $L^{\infty }$) ${ϵ}_N$ calculated on a grid of size N is estimated using the formula:

$$p_N={\displaystyle \frac{ln({ϵ}_{N/2}/{ϵ}_N)}{ln2}}.$$

It can be observed from these tables that:

• For the first two options (I-II), the C-N method converges at only first order due to the discontinuity at $S=K_M^{+}$ contained in the terminal conditions. On the other hand, mesh adaptation around the discontinuity $S=K_M^{+}$, where option prices are expected to undergo the most rapid variation, helps reduce the error by a factor of n where $n=2$ here is the refinement factor.
• For the next two options (III-IV), the C-N method converges at a higher order due to the improved smoothness of the terminal condition at $S=K_M^{+}$. The precise relationship between the convergence rate and the smoothness of the terminal condition is the subject of a future research, but these examples seem to suggest that a $C^1$ terminal condition is already sufficient to restore the second-order accuracy of the C-N method (see Remark 3.2).
• When mesh adaptation is applied to the third option, the error first improves by roughly a factor of $n=2$ up until $N=160$, after which no significant improvement in error is observed (Table 6). This is due to the fact that the terminal condition in this case has a higher degree of smoothness at $S=K_M^{+}$ (i.e., $C^0$), as a result of which the largest error on a sufficiently fine grid does not concentrate around $S=K_M^{+}$, and hence is not captured by the adaptive mesh which surrounds the strike price $S=K_M^{+}$. Similar observations apply to the fourth option, whose terminal condition has an even higher degree of smoothness at $S=K_M^{+}$ (i.e., $C^1$) and for which no improvement in error is observed for the adaptive mesh calculation at all.
• Despite these observations, the mesh adaptation technique still turns out to be useful when applied to discretely monitored options with high monitoring frequencies, where errors near the barriers caused by the non-smoothness of the terminal conditions dominate the calculations (Section 3.1, Section 3.3).

Remark 3.1.

Since the terminal conditions $V(T_M,S)$ for the first two options contain a discontinuity at $S=K_M^{+}$ and since the C-N method, with or without analytical modifications, is non-dissipative, the following modification is needed in discrete approximations $v_{N,j}$ to $V(T_M,S)$ so that the convergence of $v_{N,j}$ to $V(T_M,S)$ in a suitable weak sense (e.g., in the sense of $H^s$ for some $s>0$), and hence the convergence of $v_{0,j}$ to $V(T_0,S)$, can be ensured: let $j_0$ be chosen such that

$$ln\left(S_0\right)+(j_0-P)h\le ln\left(K_M^{+}\right)<ln\left(S_0\right)+(j_0+1-P)h,$$

or equivalently, such that

$$j_0\le P+{\displaystyle \frac{1}{h}}ln\left({\displaystyle \frac{K_M^{+}}{S_0}}\right)=:{\delta }_0<j_0+1.$$

Then for the first option, we set

$$v_{N,j}=\left\{\begin{array}{cc}0,\hfill & j<j_0\hfill \\ b_M^{+},\hfill & j>j_0\hfill \\ b_M^{+}(j_0+1-{\delta }_0),\hfill & j=j_0\hfill \end{array}\right.,$$

and for the second option, we set

$$v_{N,j}=\left\{\begin{array}{cc}0,\hfill & j<j_0\hfill \\ a_M^{+}{S}_0{e}^{(j-P)h},\hfill & j>j_0\hfill \\ {\displaystyle \frac{1}{2}}{a}_M^{+}\left[S_0{e}^{(j_0+1-P)h}+K_M^{+}\right](j_0+1-{\delta }_0),\hfill & j=j_0\hfill \end{array}\right..$$

For the next two options, the terminal conditions $V(T_M,S)$ are continuous throughout $I_0$ and hence such modification is not needed. One may simply take

$$v_{N,j}=V(T_M,S_0{e}^{(j-P)h}).$$

For dissipative schemes such as FTCS and TR-BDF2, the above modification is not needed either.

Remark 3.2.

We have carried out a preliminary error analysis for the proposed method and found that the convergence rate corresponding to a discontinuous, $C^0$, and $C^1$ terminal condition is $1/3,1$, and $5/3$, respectively. However, these results are likely suboptimal because they didn’t make full use of the dissipation properties of the Black-Scholes PDE.

Table 9, Table 10, Table 11, Table 12, Table 13 and Table 14 present results from convergence studies of the discretely monitored option (V) priced using both the original and the analytically modified C-N methods. To examine more closely the relationship between the convergence rate and the smoothness of the terminal conditions, we further distinguish between the following two cases for the analytically modified C-N method:

• A $C^0$-modification at the barriers:

  $$U_m(T_m,Z)=\left\{\begin{array}{cc}V(T_m,e^Z),\hfill & ln\left(K_m^{-}\right)<Z<ln\left(K_m^{+}\right)\hfill \\ {\tilde{V}}_m^{-},\hfill & Z\le ln\left(K_m^{-}\right)\hfill \\ {\tilde{V}}_m^{+},\hfill & Z\ge ln\left(K_m^{+}\right)\hfill \end{array}\right.,$$

  which is (2.12b) with ${\tilde{V}}_{S,m}^{+}={\tilde{V}}_{S,m}^{-}=0$ and which is only continuous at the barriers.
• A $C^1$-modification at the barriers:

  $$U_m(T_m,Z)=\left\{\begin{array}{cc}V(T_m,e^Z),\hfill & ln\left(K_m^{-}\right)<Z<ln\left(K_m^{+}\right)\hfill \\ {\tilde{V}}_{S,m}^{-}(e^Z-K_m^{-})+{\tilde{V}}_m^{-},\hfill & Z\le ln\left(K_m^{-}\right)\hfill \\ {\tilde{V}}_{S,m}^{+}(e^Z-K_m^{+})+{\tilde{V}}_m^{+},\hfill & Z\ge ln\left(K_m^{+}\right)\hfill \end{array}\right.,$$

  which is just (2.12b) and is continuously differentiable at the barriers.

For all three methods, namely, original C-N and analytically modified C-N with $C^0$- and $C^1$-modifications, the option is priced on the same interval as in the previous cases:

$$I_0:=[S_0{e}^{-6\sigma },S_0{e}^{6\sigma }]\approx [0.33131363,3.65212862].$$

The calculations are advanced with a time step $k=1/N_0$, i.e., with $N_0$ time steps per business day, and with a mesh size $h=k$ where $N_0$ ranges from 1 to 64. Since the exact option values $V(T_0,S)$ are not available, reference values $V_{0,j}$ against which all numerical approximations $v_{0,j}$ are compared are calculated on a uniform fine grid using the analytically modified C-N method with a $C^1$-modification at the barriers and with $N_0=128$ time steps per business day.

It can be observed from these tables that:

• For the original C-N method, the numerical approximations $v_{0,j}$ converge to the reference values $V_{0,j}$ in $L^2$ at an average rate of roughly 1 but do not converge in $L^{\infty }$ at all (Table 9). This reduced convergence rate or even lack of convergence is clearly a consequence of the discontinuities at the barriers contained in the terminal conditions. On the other hand, mesh adaptation around the barriers, where option prices are expected to undergo the most rapid variation, helps reduce the error by a factor of n where $n=2$ here is the refinement factor (Table 10).
• For the analytically modified C-N method, with either $C^0$- or $C^1$-modifications, the numerical approximations $v_{0,j}$ converge to the reference values $V_{0,j}$ in both $L^2$ and $L^{\infty }$ and do so at a higher rate due to the improved smoothness of the terminal conditions at the barriers. More specifically, with $C^0$-modifications, $v_{0,j}$ converge to $V_{0,j}$ in $L^2$ at an average rate of roughly $3/2$ and in $L^{\infty }$ at an average rate of roughly 1 (Table 11), while with $C^1$-modifications, the convergence in both $L^2$ and $L^{\infty }$ has an average rate of roughly 2 (Table 13). In other words, these examples seem to suggest that a $C^1$-modification to the terminal conditions is already sufficient to restore the second-order accuracy of the C-N method.
• When mesh adaptation is applied with $C^0$-modifications, both the $L^2$- and $L^{\infty }$-errors improve by roughly a factor of $n^{3/2}=2^{3/2}\approx 2.8284$ for all $N_0$, demonstrating the effectiveness of the adaptive mesh (Table 12).
• When mesh adaptation is applied with $C^1$-modifications, the numerical approximations $v_{0,j}$ calculated on the finest grid with $N_0=64$ are almost identical with the reference values $V_{0,j}$ calculated on the uniform fine grid with $N_0=128$, as suggested by the unusually small errors (last row, Table 14). On the other hand, for numerical approximations $v_{0,j}$ calculated on coarser grids with $N_0\le 32$, both the $L^2$- and $L^{\infty }$-errors improve by roughly a factor of $n^2=2^2=4$ for all $N_0$ when mesh adaptation is enabled. This provides another strong evidence for the effectiveness of the adaptive mesh, especially on discretely monitored options with high monitoring frequencies.
• We remark that for mesh adaptation applied with $C^1$-modifications, a sufficiently large area around the barriers needs to be refined in order for the adaptive mesh to be effective, in view of the improved smoothness of the terminal conditions near the barriers. Indeed, when the refinement ratio in Table 14 was reduced from 15% to 10%, meaning only the $(2P+1)\times 10\%$ grid points around each barrier are refined, the resulting adaptive mesh was only able to produce marginal improvement in error for $N_0\le 4$ and was not able to produce any meaningful improvement in error at all for all larger $N_0$. This is not an issue when mesh adaptation is applied with $C^0$-modifications or with no analytical modifications, where the errors near the barriers caused by the non-smoothness of the terminal conditions dominate the calculations.

3.2. Comparison with other Numerical Methods

To assess the efficacy and efficiency of the proposed method on options with discretely monitored structures, we next apply other popular, (formally) second-order numerical methods such as C-N with Rannacher time-stepping (CN-RAN) and trapezoidal rule with second-order backward differentiation formula (TR-BDF2) to the discretely monitored option introduced in the previous section. Both methods (CN-RAN and TR-BDF2) have been applied to the four standard options introduced above, and their convergence properties have been found to be very similar to those of the original C-N method (Table 1, Table 3, Table 5 and Table 7).

Table 15 and Table 16 present results from convergence studies of the discretely monitored option priced using the CN-RAN and the TR-BDF2 methods. It can be observed from these tables that:

• Both methods exhibit very similar convergence behaviors and converge in both $L^2$ and $L^{\infty }$ at an average rate of roughly 1. This shows, in particular, that both methods lose their second-order accuracy when applied to problems containing discontinuities.
• For a given time step k and mesh size h, both methods produce results that are more accurate than the original C-N method (Table 9) but less accurate than the analytically modified C-N method (Table 11 and Table 13).

As a more direct comparison of the different numerical methods studied in the literature (including in this work) and commonly used in practice, we list in Table 17 the pricing results of the following five methods applied to the discretely monitored option introduced in the previous section:

• Analytically modified C-N with a $C^1$-modification (CN-C1).
• C-N with Rannacher time-stepping (CN-RAN).
• Trapezoidal rule with second-order backward differentiation formula (TR-BDF2).
• Forward-time centered-space (or implicit Euler) (FTCS).
• Monte-Carlo simulations (MC).

Table 17. Comparison of different methods on a discretely monitored option.

Time Steps	CN-C1		CN-RAN
$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	CPU time	$V(T_0,S_0)\mathrm{error}$	CPU time
1	$6.7233\times {10}^{-4}$	$0.5618$	$1.1520\times {10}^{-1}$	$0.1195$
2	$3.1867\times {10}^{-5}$	$0.6440$	$5.5935\times {10}^{-2}$	$0.2298$
4	$1.9737\times {10}^{-5}$	$0.8517$	$2.4723\times {10}^{-2}$	$0.3972$
8	$6.3714\times {10}^{-6}$	$1.4262$	$1.1473\times {10}^{-2}$	$0.8562$
16	$1.8043\times {10}^{-6}$	$2.8770$	$5.5642\times {10}^{-3}$	$2.1088$
32	$4.5282\times {10}^{-7}$	$7.3977$	$2.7281\times {10}^{-3}$	$8.5549$
64	$9.4676\times {10}^{-8}$	$23.4999$	$1.3534\times {10}^{-3}$	$46.2784$
Reference option value at $S_0$: 1.83652751.
$\mathrm{Time}\mathrm{steps}$	TR-BDF2		FTCS
$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	CPU time	$V(T_0,S_0)\mathrm{error}$	CPU time
1	$8.5771\times {10}^{-2}$	$0.1145$	$1.3770\times {10}^{-1}$	$0.0648$
2	$4.4429\times {10}^{-2}$	$0.2203$	$6.7278\times {10}^{-2}$	$0.1141$
4	$2.1929\times {10}^{-2}$	$0.4887$	$3.3502\times {10}^{-2}$	$0.2415$
8	$1.0782\times {10}^{-2}$	$1.2014$	$1.6546\times {10}^{-2}$	$0.5837$
16	$5.3920\times {10}^{-3}$	$3.8852$	$8.2700\times {10}^{-3}$	$1.6694$
32	$2.6851\times {10}^{-3}$	$19.8665$	$4.1230\times {10}^{-3}$	$5.1475$
64	$1.3427\times {10}^{-3}$	$96.9560$	$2.0614\times {10}^{-3}$	$21.6645$
Reference option value at $S_0$: 1.83652751.
	Number	MC
	$\mathrm{of}\mathrm{paths}$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$
	40000	$2.6953\times {10}^{-2}$	$13.1037$
	60000	$5.6742\times {10}^{-3}$	$17.2020$
	80000	$6.3685\times {10}^{-3}$	$22.7554$
	100000	$3.2074\times {10}^{-3}$	$28.0916$
	120000	$1.6717\times {10}^{-3}$	$33.6178$
	140000	$3.2980\times {10}^{-4}$	$38.8197$
	160000	$3.2803\times {10}^{-3}$	$44.4004$
Reference option value at $S_0$: 1.83652751.

Table 17. Comparison of different methods on a discretely monitored option.

Time Steps	CN-C1		CN-RAN
$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	CPU time	$V(T_0,S_0)\mathrm{error}$	CPU time
1	$6.7233\times {10}^{-4}$	$0.5618$	$1.1520\times {10}^{-1}$	$0.1195$
2	$3.1867\times {10}^{-5}$	$0.6440$	$5.5935\times {10}^{-2}$	$0.2298$
4	$1.9737\times {10}^{-5}$	$0.8517$	$2.4723\times {10}^{-2}$	$0.3972$
8	$6.3714\times {10}^{-6}$	$1.4262$	$1.1473\times {10}^{-2}$	$0.8562$
16	$1.8043\times {10}^{-6}$	$2.8770$	$5.5642\times {10}^{-3}$	$2.1088$
32	$4.5282\times {10}^{-7}$	$7.3977$	$2.7281\times {10}^{-3}$	$8.5549$
64	$9.4676\times {10}^{-8}$	$23.4999$	$1.3534\times {10}^{-3}$	$46.2784$
Reference option value at $S_0$: 1.83652751.
$\mathrm{Time}\mathrm{steps}$	TR-BDF2		FTCS
$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	CPU time	$V(T_0,S_0)\mathrm{error}$	CPU time
1	$8.5771\times {10}^{-2}$	$0.1145$	$1.3770\times {10}^{-1}$	$0.0648$
2	$4.4429\times {10}^{-2}$	$0.2203$	$6.7278\times {10}^{-2}$	$0.1141$
4	$2.1929\times {10}^{-2}$	$0.4887$	$3.3502\times {10}^{-2}$	$0.2415$
8	$1.0782\times {10}^{-2}$	$1.2014$	$1.6546\times {10}^{-2}$	$0.5837$
16	$5.3920\times {10}^{-3}$	$3.8852$	$8.2700\times {10}^{-3}$	$1.6694$
32	$2.6851\times {10}^{-3}$	$19.8665$	$4.1230\times {10}^{-3}$	$5.1475$
64	$1.3427\times {10}^{-3}$	$96.9560$	$2.0614\times {10}^{-3}$	$21.6645$
Reference option value at $S_0$: 1.83652751.
	Number	MC
	$\mathrm{of}\mathrm{paths}$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$
	40000	$2.6953\times {10}^{-2}$	$13.1037$
	60000	$5.6742\times {10}^{-3}$	$17.2020$
	80000	$6.3685\times {10}^{-3}$	$22.7554$
	100000	$3.2074\times {10}^{-3}$	$28.0916$
	120000	$1.6717\times {10}^{-3}$	$33.6178$
	140000	$3.2980\times {10}^{-4}$	$38.8197$
	160000	$3.2803\times {10}^{-3}$	$44.4004$
Reference option value at $S_0$: 1.83652751.

For each method, the (pointwise) error of the option value at $S_0$, together with the CPU time (in seconds)3 used to obtain the pricing result, are shown for different time steps $N_0$. Each CPU time is measured by running the corresponding code 5 times and dividing the total CPU time by 5. It can be observed from the table that for a given level of accuracy, the CN-C1 method requires the least amount of CPU time while all other four methods typically require significantly more CPU time. As a specific example, consider the accuracy level of $2\times {10}^{-3}$. A more refined study shows that to achieve this level of accuracy, the CN-C1 method requires less than 0.5618 seconds while the other four methods require

• 16.7132 seconds (CN-RAN),
• 36.1549 seconds (TR-BDF2),
• 22.8074 seconds (FTCS), and
• 26.4529 seconds (MC),

respectively (Table 18). This confirms the efficiency of the CN-C1 method.

Table 18. CPU time (in seconds) required to achieve an accuracy level of $2\times {10}^{-3}$ by different methods on a discretely monitored option.

Time Steps	CN-C1		Time Steps	CN-RAN
$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$	$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$
1	$6.7233\times {10}^{-4}$	$0.5618$	43	$2.0210\times {10}^{-3}$	$16.2035$
−	−	−	44	$1.9749\times {10}^{-3}$	$17.3225$
−	−	−	$43.{4555}^{*}$	$2.0000\times {10}^{-3}$	$16.{7132}^{*}$
Reference option value at $S_0$: 1.83652751.
*: Estimates obtained from linear interpolation.
$\mathrm{Time}\mathrm{steps}$	TR-BDF2		$\mathrm{Time}\mathrm{steps}$	FTCS
$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$	$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$
42	$2.0449\times {10}^{-3}$	$34.4220$	65	$2.0309\times {10}^{-3}$	$21.8704$
43	$1.9971\times {10}^{-3}$	$36.2668$	66	$1.9978\times {10}^{-3}$	$22.8741$
$42.{9393}^{*}$	$2.0000\times {10}^{-3}$	$36.{1549}^{*}$	$65.{9335}^{*}$	$2.0000\times {10}^{-3}$	$22.{8074}^{*}$
Reference option value at $S_0$: 1.83652751.
*: Estimates obtained from linear interpolation.
	Number	MC
	$\mathrm{of}\mathrm{paths}$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$
	108300	$2.1369\times {10}^{-3}$	$26.7155$
	108400	$1.6615\times {10}^{-3}$	$25.8035$
	$108328.{80}^{*}$	$2.0000\times {10}^{-3}$	$26.{4529}^{*}$
	Reference option value at $S_0$: 1.83652751.
*: Estimates obtained from linear interpolation.

Table 18. CPU time (in seconds) required to achieve an accuracy level of $2\times {10}^{-3}$ by different methods on a discretely monitored option.

Time Steps	CN-C1		Time Steps	CN-RAN
$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$	$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$
1	$6.7233\times {10}^{-4}$	$0.5618$	43	$2.0210\times {10}^{-3}$	$16.2035$
−	−	−	44	$1.9749\times {10}^{-3}$	$17.3225$
−	−	−	$43.{4555}^{*}$	$2.0000\times {10}^{-3}$	$16.{7132}^{*}$
Reference option value at $S_0$: 1.83652751.
*: Estimates obtained from linear interpolation.
$\mathrm{Time}\mathrm{steps}$	TR-BDF2		$\mathrm{Time}\mathrm{steps}$	FTCS
$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$	$\mathrm{per}\mathrm{day}{N}_0$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$
42	$2.0449\times {10}^{-3}$	$34.4220$	65	$2.0309\times {10}^{-3}$	$21.8704$
43	$1.9971\times {10}^{-3}$	$36.2668$	66	$1.9978\times {10}^{-3}$	$22.8741$
$42.{9393}^{*}$	$2.0000\times {10}^{-3}$	$36.{1549}^{*}$	$65.{9335}^{*}$	$2.0000\times {10}^{-3}$	$22.{8074}^{*}$
Reference option value at $S_0$: 1.83652751.
*: Estimates obtained from linear interpolation.
	Number	MC
	$\mathrm{of}\mathrm{paths}$	$V(T_0,S_0)\mathrm{error}$	$\mathrm{CPU}\mathrm{time}$
	108300	$2.1369\times {10}^{-3}$	$26.7155$
	108400	$1.6615\times {10}^{-3}$	$25.8035$
	$108328.{80}^{*}$	$2.0000\times {10}^{-3}$	$26.{4529}^{*}$
	Reference option value at $S_0$: 1.83652751.
*: Estimates obtained from linear interpolation.

Remark 3.3.

It may appear from these results that MC is even more efficient that the TR-BDF2 method. The reality is that, if higher accuracy is demanded (higher than $2\times {10}^{-3}$), then the TR-BDF2 method will definitely outperform MC, even though it is relatively slow compared with other finite difference methods considered here.

3.3. Effectiveness of Adaptive Mesh Refinement

We also demonstrate the effectiveness of the mesh adaptation technique by applying the CN-C1 method to the discretely monitored option introduced in the previous sections, both with and without mesh adaptation. The result shows that, with a mesh refinement factor of $n=8$, an adaptive mesh calculation with $N_0$ time steps per business day is able to achieve an accuracy level comparable to that of a uniform grid calculation with $8N_0$ time steps per business day, with a more than 60% save in CPU time (Table 19). This confirms the efficiency of the mesh adaptation technique.

3.4. Approximations to Derivatives

It is well known that finite difference methods applied to options with non-smooth terminal conditions tend to produce poor approximations to option Greeks such as Delta $\Delta$ and Gamma $\Gamma$. In fact, as demonstrated in [11], when the original C-N method is applied to a European call option, the discrete approximations to $\Delta$ and $\Gamma$ both fail to converge, and the discrete approximation to $\Gamma$ even grows unboundedly when the computational grid is refined. To examine the quality of the discrete approximations to option Greeks produced by the proposed method, we calculate $\Delta$ and $\Gamma$ for the discretely monitored option introduced in the previous section, using the formulas

$$\begin{array}{cc}\hfill {\Delta }_{0,j}& \approx {\displaystyle \frac{V_{0,j+1}-V_{0,j-1}}{2h{S}_j}},\hfill \\ \hfill {\Gamma }_{0,j}& \approx {\displaystyle \frac{V_{0,j+1}-2V_{0,j}+V_{0,j-1}}{h^2{S}_j^2}}-{\displaystyle \frac{V_{0,j+1}-V_{0,j-1}}{2h{S}_j^2}},\hfill \end{array}$$

where $V_{0,j}$ are the same reference option values calculated on a uniform fine grid using the CN-C1 method with $N_0=128$ time steps per business day. The approximated Greeks are shown in Figure 2 and Figure 3 together with the option values, where it can be seen that while the discrete approximation to $\Delta$ is free of oscillations, the discrete approximation to $\Gamma$ contains a small amount of oscillations near the lower barrier 0.9 which is monitored daily. These oscillations are likely consequences of the $O\left(h^2\right)$ errors introduced by the linear interpolation approximations to $V,\partial V/\partial S$ at the barriers (see (2.17)) which have been significantly amplified by the high-frequency monitoring at the lower barrier. To remove these high-wavenumber oscillations, we replace the first C-N step after each monitoring point by two half-step dissipative steps such as FTCS (implicit Euler) steps and apply the resulting method, coded CN-C1-RAN1, again to the discretely monitored option introduced in the previous section. The results (Figure 4) demonstrate the effectiveness of the modified method in the elimination of high-wavenumber oscillations and also confirm the second-order convergence of both the option value V and option Greeks $\Delta ,\Gamma$ in both $L^2$ and $L^{\infty }$ (Table 20).

Remark 3.4.

When applied with mesh adaptation, the CN-C1-RAN1 method produces option Greeks $\Delta ,\Gamma$ converging at reduced convergence rates, most likely due to the small mismatches among different parts of ${\widehat{U}}_m(t_{l-1},Z_j)$ introduced by the direct transfer of the fine grid solutions ${\mathbf{W}}_{l-1,0}^{\pm m}$ to coarse grid (see (2.27)). While this issue may be addressed using, say, a high-order smoothing operator, we decide not to pursue this direction further but instead recommend the readers to consider CN-C1-RAN1 (without mesh adaptation) if accurate approximation to option Greeks is important and to consider CN-C1 (with or without mesh adaptation) if option Greeks are not needed or not important.

4. Discussion

We have introduced an analytical modification for finite difference methods used in the pricing of discretely monitored options in the Black-Scholes framework. The modified finite difference methods can be used to price single and double barrier options, autocallable structured products, and snowball options. We have demonstrated the efficacy and efficiency of the proposed method on the traditional Crank-Nicolson scheme which, in its analytically modified form CN-C1 or CN-C1-RAN1, is capable of producing much more accurate pricing results with much lower CPU cost than other popular option pricing methods such as CN-RAN, TR-BDF2, and Monte-Carlo simulations. We have also introduced a mesh adaptation technique for the analytically modified Crank-Nicolson method, which is capable of achieving an accuracy level comparable to that of a uniform fine grid calculation with a more than 60% improvement in computational speed. We strongly believe that the proposed method and the related ideas will prove a valuable tool for practitioners working with the aforementioned types of options.

This work can be extended in several directions. First, it is desirable to carry out a rigorous error analysis and obtain a precise quantification of the relationship between convergence rates and smoothness of terminal conditions for the proposed method. We have done some preliminary analysis along this direction (Remark 3.2) and the full analysis is currently a work in progress. Second, the proposed method has assumed a volatility that is independent of the asset price S, which may be too restrictive in some practical situations. It is desirable to relax this assumption and apply similar analytical modifications to option pricing problems with variable (in S) volatility. This is the subject of a current research and the results will be reported in a forthcoming paper.