1 Introduction

In this paper, we seek to price discrete options with two barriers. In particular, a double-barrier option is a derivative contract that is activated or extinguished when the underlying asset goes far beyond a pre-specified price interval at least once during the option life. Options with barrier features are cheaper than plain options and hence are common in structured products. As mentioned in Bernard et al. (2008), 25.7% of equity-linked notes have their payments driven by a triggered event based on the trajectory of the underlying stock.

The pricing of double-barrier options has been studied for quite a long time, and dates back to Kunitomo and Ikeda (1992). Since then, various studies have been conducted, such as Geman and Yor (1996); Pelsser (2000) under the Black–Scholes (BS) model and Cai et al. (2009) under a jump diffusion model. A recent development for Lévy-driven models is found in Boyarchenkoy and Levendorskiĭ (2012). These works all posit that the barriers are monitored continuously. However, in most real cases, barriers are imposed only at discrete points of time. Such discretely monitored path-dependent options are usually referred to as discrete options.

A pioneer and famous work in mathematical finance for discrete options is Broadie et al. (1997), which documents that under the BS setting, there are significant price errors between barrier options with and without continuous monitoring. More importantly, they further show such errors can be mitigated by a technique called continuity correction. Specifically, the price of a discrete barrier option can be approximately determined by the corresponding continuous-time pricing formula with a simple shift in the barrier level. The amount of shift is related to the expectation of “overshoot,” the excess price of the underlying asset over the barrier when such a breach is observed for the first time. As far as we know, this is the only technique that gives analytical approximations; otherwise, we resort to numerical methods; see (Milev & Tagliani, 2010; Golbabai et al., 2014; Ahmadian & Ballestra, 2015; Ahmadian et al., 2021) for references.Footnote 1

Other applications of continuity correction to option pricing include (Broadie et al., 1999) for lookback options, Lai et al. (2007) for American options, Fuh et al. (2013) for a jump-diffusion model, and Dia and Lamberton (2011) for exponential Lévy models. Nevertheless, these studies concern only single barriers. For two barriers at a time, it is generally expected that applying the same result in Broadie et al. (1997) to each barrier separately would be sufficient. This is also conjectured in Hörfelt (2003), but we still lack solid estimates of the approximation error. Hence, in this article we attempt to fill this gap by showing that the correction for two simultaneous barriers is different from that for a single barrier. The method involves not only the expected overshoot but also the probabilities of hitting the barrier.

Indeed, continuity correction for so-called chained options was established in Jun (2013), which deals with a typical barrier option but activated by the hitting of another primary barrier. For example, the \(DIC _u\) is a down-and-in call option activated when the price of the underlying asset hits an upper barrier. This case belongs to a type of option with two barriers, but is not exactly the same as what we here call double-barrier options. Other types of options with two barriers can be found in Li (1999). Notice that correction for chained options seems to be more straightforward because it allows the reflection principle to be sequentially applied at each barrier. As for the double barrier, we instead follow the PDE (partial differential equation) approach initiated by Keener (2013), who (implicitly) introduces the principle of smooth fitFootnote 2 from free-boundary problems (see Shepp & Shiryaev 1993 for example) to the study of boundary effects such as ours. The central idea of Keener’s approach is that only “non-smoothness” on boundaries necessitates correction for the overshoot effect in order to improve the convergence brought about by discretization in the time domain. We here extend his technique from one boundary to two simultaneous boundary conditions.

The idea of continuity correction can also be utilized in the reverse fashion. For instance, if in a continuous-time structural form model we seek to estimate a default probability via a simulation approach, which is clearly of a discrete nature, then we can shrink the distance to default to improve the efficiency of estimation. This practice is usually termed boundary correction. One relevant example can be found in Gobet and Menozzi (2010), where the authors begin with a multidimensional framework with a finite time domain and consider a general diffusion process. They also conclude that one should shift the boundary inward along a normal by the same amount as the aforementioned expected overshoot when adopting a Euler discretization scheme in a Monte-Carlo simulation to estimate an expectation of a continuous entity. In this study, by contrast, we borrow a continuous-time formula through a joint Laplace transform to approximate its discrete counterpart. As a result, we face a problem with an infinite time domain which then yields different correction results. More briefly, their boundary-correction technique in a numerical scheme is not totally equivalent to the continuity correction we seek. Here, our concern is more about how the correction should be incorporated into the existing formula, which in turn requires further properties with respect to asymptotic independence.

This work contributes to the literature in three main aspects. First, this is the first study to establish continuity correction for the pricing of double-barrier options. The correction result also shifts the barrier but may do so either farther or closer, relatively to the initial starting point. Additionally, the shifted amount for each barrier can be uneven. We do not observe both in a single-barrier case. Second, the boundary effect from the overshoot is no longer sufficiently captured by a linear expansion; now, to simultaneously account for non-smoothness at two points, it is instead a cubic polynomial; see (2.19). Nevertheless, the price approximation achieved between a discrete double-barrier option and its continuous correspondent is still a first-order approximation in light of renewal theory. Namely, the improved convergence rate is merely \(o\left( 1/\sqrt{m}\right)\) per the literature. Here m stands for the boundary monitoring frequency during the contract life. Last, we extend asymptotic properties of a typical renewal measure to the double-boundary case; see a battery of lemmas in the appendix.

The rest of the article is organized as follows. Section 2 introduces the working model and formulates our correction problem relevant to the pricing of double-barrier options. Section 3 then presents the main theoretical results, especially of continuity correction for the joint moment generating function of the discrete first passage time and its associated overshoot. Simple illustrative numerical examples are given in Sect. 4, and the final section concludes the paper. Some proofs are deferred to the appendices.

2 Problem formulation

Here, we will formulate the correction problem associated with the pricing of discrete double-barrier options, during which we will also illustrate the methodology and outline the key steps to achieve our final correction goal.

2.1 Model

We start with the classic Black–Scholes (BS) model; that is to say, we consider an underlying asset \(S_t\) whose price dynamics under the risk-neutral measure are

$$\begin{aligned} dS_t =r S_t\,dt+\sigma S_t\, dW_t, \end{aligned}$$

where r is the constant risk-free rate, \(\sigma\) is the constant volatility, and \(W_t\) is a one-dimensional standard Brownian motion. Let \(X_t=\log (S_t/S_0)\) denote the log-return process; then it is well-known that

$$\begin{aligned} X_t = \left( r-\sigma ^2/2\right) t+\sigma \,W_t. \end{aligned}$$
(2.1)

Also, by the celebrated Lévy–Khintchine formula, we have for \(\theta \in {\mathbb {R}}\), \(E[e^{\theta \,X_t}]=\exp \left\{ G(\theta )\,t\right\}\) with

$$\begin{aligned} G(\theta ) = \frac{1}{2}\,\sigma ^2\,\theta ^2+\left( r-\sigma ^2/2\right) \theta . \end{aligned}$$
(2.2)

Next, it is easy to see that for any \(\alpha >0\), the equation \(G(\theta )=\alpha\) has exactly two roots:

$$\begin{aligned} \begin{aligned} \beta _{1,\alpha }&=\frac{1}{\sigma ^2}\left[ -\left( r-\sigma ^2/2\right) -\sqrt{\left( r-\sigma ^2/2\right) ^2+2\,\sigma ^2\,\alpha }\,\right] ; \\ \beta _{2,\alpha }&=\frac{1}{\sigma ^2}\left[ -\left( r-\sigma ^2/2\right) +\sqrt{\left( r-\sigma ^2/2\right) ^2+2\,\sigma ^2\,\alpha }\,\right] . \end{aligned} \end{aligned}$$
(2.3)

Clearly, \(\beta _{1,\alpha }<0<\beta _{2,\alpha }\). In addition, we have \(\beta _{2,\alpha }>1\) if \(\alpha >r\) and

$$\begin{aligned} \begin{aligned} \beta _1^*:=\beta _{1,0}=\lim \limits _{\alpha \downarrow 0}\beta _{1,\alpha }&=\left\{ \begin{array}{ll} -2\gamma , &{} \text{ if }\, r-\sigma ^2/2\ge 0; \\ 0, &{} \text{ if }\, r-\sigma ^2/2<0; \end{array}\right. \\ \beta _2^*:=\beta _{2,0}=\lim \limits _{\alpha \downarrow 0}\beta _{2,\alpha }&=\left\{ \begin{array}{ll} 0, &{} \text{ if }\, r-\sigma ^2/2\ge 0; \\ -2\gamma , &{} \text{ if }\, r-\sigma ^2/2<0, \end{array}\right. \end{aligned} \end{aligned}$$
(2.4)

in which

$$\begin{aligned} \gamma =\left( r-\sigma ^2/2\right) /\sigma ^2. \end{aligned}$$
(2.5)

Note that \(\gamma\) has the same sign as the drift of \(X_t\). Finally, we will use \({\mathcal {F}}_t\) to denote the natural filtration generated by the process \(X_t\).

2.2 Option pricing

Throughout this paper, we will take the pricing of a knock-in double-barrier call (KIC) as our illustrative example to demonstrate the central idea and the correction results. The payoff of a such option at the maturity date T (\(>0\)) generally looks like

$$\begin{aligned} \chi ^{ kic }=\left( S_T-K\right) ^{+}\,1_{\{\tau \le T\}}, \end{aligned}$$

where K is the strike price, \(1_{\{\cdot \}}\) stands for the indicator function, and \(\tau\) is the associated first-hitting time. If the barriers are monitored continuously, \(\tau\) is

$$\begin{aligned} \tau _{\infty }=\inf \left\{ t\ge 0:\;S_t\notin \left( L,H\right) \right\} \end{aligned}$$

for two constants (barriers) LH such that \(L<S_0<H\). Otherwise, it is

$$\begin{aligned} \tau _m=\inf \left\{ n\in {\mathbb {N}}:\;S_{n\varDelta t}\notin \left( L,H\right) \right\} , \end{aligned}$$

where \(\varDelta t=T/m\) for some given monitoring frequency m. Namely, the barriers are imposed only at m discrete points of time: \(\varDelta t\), \(2\varDelta t\), \(\ldots\), \(m\varDelta t=T\).

Note that both \(\tau _{\infty }\) and \(\tau _m\) are finite almost surely. However, in order to assure finite expectations as well, from now on we further assume

$$\begin{aligned} r-\sigma ^2/2\ne 0\quad (\text{ i.e., }\;\;\gamma \ne 0) \end{aligned}$$
(2.6)

without loss of generality.Footnote 3

Let \(X>0\) be a scaling factor and \(k=\log (X/K)\); then the price of a KIC can be expressed as

$$\begin{aligned} KIC (T,k)=Xe^{-rT}E\left[ \left( S_T/X-e^{-k}\right) ^{+}\,1_{\{\tau \le T\}}\right] . \end{aligned}$$

Similar to Cai et al. (2009), we will compute this quantity in terms of a two-dimensional Laplace transform. Specifically, consider

$$\begin{aligned} \widehat{ KIC }(\alpha ,\theta ) =\int _0^\infty \int _{-\infty }^{\infty }\,e^{-\alpha \,T-\theta \,k}\, KIC (T,k)\,dk\,dT \end{aligned}$$

and the price then can be obtained via the inversion algorithm as proposed by Cai and Kou (2012). Actually, after some algebraic manipulation, the last expression can be further simplified asFootnote 4

$$\begin{aligned} \widehat{ KIC }(\alpha ,\theta )=\frac{X^{-\theta }}{\theta (\theta +1)}\, \frac{S_0^{\theta +1}}{r+\alpha -G(\theta +1)}\,E\left[ e^{-\left( r+\alpha \right) \tau +\left( \theta +1\right) X_\tau }\right] , \end{aligned}$$
(2.7)

provided that \(0<\theta <\left( \beta _{2,r+\alpha }-1\right)\) and \(\alpha >\max \left( G(\theta +1)-r,0\right)\). Here \(G(\cdot )\) is defined in (2.2) and X being introduced is merely to satisfy the parameter constraints behind (2.7) when applying the inversion algorithm.

Therefore, to study the price of a KIC, one needs only the joint moment generating function (MGF) of \(\tau\) and \(X_\tau\). This is true for both continuous (\(\tau _{\infty }\)) and discrete (\(\tau _m\)) options. Hence, from now on, let

$$\begin{aligned} f(t,x)=\exp \{-\alpha \,t+\theta \,x\}, \end{aligned}$$
(2.8)

which will serve as our working function in the subsequent analyses. We will also confine our discussion to the case with

$$\begin{aligned} \alpha >0\quad \text{ and }\quad \beta _{1,\alpha }<\theta <\beta _{2,\alpha } \end{aligned}$$
(2.9)

to assure the \(L^1\) finiteness to come. Finally, note that the events \(\left\{ S_t\notin \left( L,H\right) \right\}\) and \(\left\{ X_t\notin \left( a,b\right) \right\}\) are equivalent where \(a=\log (L/S_0)<0\) and \(b=\log (H/S_0)>0\).

2.3 Continuous-time joint distribution

For \(0\le t\le \tau _{\infty }\), consider

$$\begin{aligned} M_t=E\left[ f(\tau _\infty ,X_{\tau _\infty })\,|{\mathcal {F}}_t\right] . \end{aligned}$$

Then, \(M_t\) is a Doob’s martingale, and by the Markovian property there exists a function v(tx) such that \(M_t=v(t,X_t)\). In fact, the function v admits the form

$$\begin{aligned} v(t,x)=\left\{ \begin{array}{ll} u(t,x), &{} x\in \left( a,b\right) ; \\ f(t,x), &{} \text{ otherwise }, \end{array} \right. \end{aligned}$$

in which the function u satisfies the following PDE:

$$\begin{aligned} u_t+\left( r-\sigma ^2/2\right) u_x+\frac{1}{2}\,\sigma ^2\,u_{xx}=0 \end{aligned}$$
(2.10)

subject to the boundary conditions

$$\begin{aligned} u(t,a)=f(t,a)\quad \text{ and }\quad u(t,b)=f(t,b). \end{aligned}$$
(2.11)

For convenience, we will simply refer to such f as a terminal function to emphasize its role in the boundary conditions.

Recalling (2.8), u(0, 0) is clearly the desired joint MGF of \(\tau _{\infty }\) and \(X_{\tau _{\infty }}\). As suggested in Cai et al. (2009), one can easily check that the solution to (2.10) is given by

$$\begin{aligned} u(t,x)=e^{-\alpha \,t}\left[ A(\alpha ,\theta )\,e^{\theta \,a+\beta _{1,\alpha }\left( x-a\right) } +B(\alpha ,\theta )\,e^{\theta \,b+\beta _{2,\alpha }\left( x-b\right) }\right] , \end{aligned}$$
(2.12)

where \(\beta _{1,\alpha }\) and \(\beta _{2,\alpha }\) are given by (2.3), and

$$\begin{aligned} A(\alpha ,\theta )=\frac{1-e^{(b-a)\left( \theta -\beta _{2,\alpha }\right) }}{1-e^{(b-a) \left( \beta _{1,\alpha }-\beta _{2,\alpha }\right) }}\quad \text{ and }\quad B(\alpha ,\theta ) =\frac{1-e^{(b-a)\left( \beta _{1,\alpha }-\theta \right) }}{1-e^{(b-a)\left( \beta _{1,\alpha } -\beta _{2,\alpha }\right) }}. \end{aligned}$$
(2.13)

Thus, we have

$$\begin{aligned} E\left[ e^{-\alpha \tau _\infty +\theta X_{\tau _\infty }}\right] =u(0,0)=A(\alpha ,\theta )\, e^{a\left( \theta -\beta _{1,\alpha }\right) }+B(\alpha ,\theta )\,e^{b\left( \theta -\beta _{2,\alpha }\right) }. \end{aligned}$$
(2.14)

Now, as a final remark, let us list two implications of this joint MGF of the continuous version, which will be used greatly in the sequel. First, the MGF of the continuous stopping time \(\tau _{\infty }\) can then be obtained by directly setting \(\theta =0\) in (2.14), which leads to

$$\begin{aligned} \phi (\alpha ):=E\left[ e^{-\alpha \tau _{\infty }}\right] =A(\alpha ,0)\,e^{-a\beta _{1,\alpha }} +B(\alpha ,0)\,e^{-b\beta _{2,\alpha }} \end{aligned}$$
(2.15)

for \(\alpha >0\). On the other hand, when letting \(\alpha\) go to zero, together with the help of (2.4), we have

$$\begin{aligned} \xi (\theta ):=E\left[ e^{\theta X_{\tau _{\infty }}}\right] =A(0,\theta )\,e^{a\left( \theta -\beta _1^*\right) }+B(0,\theta )\,e^{b\left( \theta -\beta _2^*\right) }, \end{aligned}$$

from which we further achieve

$$\begin{aligned} E\left[ X_{\tau _{\infty }}\right] =\xi '(0)=a\cdot q_a+b\cdot q_b \end{aligned}$$

with

$$\begin{aligned} \begin{aligned} q_a&=\frac{e^{-a\beta _1^*}\left[ 1-e^{b\left( \beta _1^*-\beta _2^*\right) }\right] }{1-e^{(b-a)\left( \beta _1^*-\beta _2^*\right) }}=\frac{1-e^{2\gamma b}}{1-e^{2\gamma \left( b-a\right) }};\\ q_b&=\frac{e^{-b\beta _2^*}\left[ 1-e^{-a\left( \beta _1^*-\beta _2^*\right) }\right] }{1-e^{(b-a)\left( \beta _1^*-\beta _2^*\right) }}=\frac{1-e^{2\gamma a}}{1-e^{2\gamma \left( a-b\right) }}. \end{aligned} \end{aligned}$$
(2.16)

In fact, \(q_a=P(X_{\tau _{\infty }}=a)\) and \(q_b=P(X_{\tau _{\infty }}=b).\) Thus \(q_a\ge 0\), \(q_b\ge 0\) and \(q_a+q_b=1\), which can also be checked directly. Moreover, \(\lim _{b\rightarrow \infty }q_a=\lambda _a\) and \(\lim _{a\rightarrow -\infty }q_b=\lambda _b\), where

$$\begin{aligned} \lambda _y:=\lambda (y)=\min \left\{ 1,\,e^{2\gamma y}\right\} . \end{aligned}$$
(2.17)

Note that \(\lambda _y\) is merely the probability that a standard Brownian motion with drift \(\sigma \,\gamma\) hits the boundary of level y. Although these results are rather standard, we include them here for future use.

2.4 Continuity correction

With discrete monitoring, however, it is generally not possible to obtain a closed-form expression for \(E\left[ f(\tau _m\varDelta t,X_{\tau _m\varDelta t})\right]\). Thereby, we seek an analytical approximation. Given the property of weak convergence, a straightforward approach is to adopt

$$\begin{aligned} E\left[ f(\tau _m\varDelta t, X_{\tau _m\varDelta t})\right] =E\left[ f(\tau _\infty , X_{\tau _\infty })\right] +o(1)=u(0,0)+o(1) \end{aligned}$$
(2.18)

as \(m\rightarrow \infty\). However, in practice, the convergence rate of (2.18) is too slow for accurate approximation, even given daily monitoring.

Our goal here is to improve the convergence rate in (2.18) from o(1) to \(o(1/\sqrt{m})\) via adjustments to the levels of boundaries a and b before applying the formula (2.14). In sequential analysis, this is typically termed first-order approximation. The intuition is to shift away the boundary levels to account for a lower boundary-crossing probability with discrete monitoring. The amount of correction is highly related to the expectation of the so-called overshoot, the amount of the stopped level of the underlying process in excess of the boundary. Note that under the BS model, overshoot occurs only when monitoring discretely.

Keener (2013) observes that the boundary effect (due to overshoot) only matters (in the first order) when criteria like smooth pasting do not hold at the boundaries in the PDE problem (2.10). Namely, \(\frac{\partial }{\partial x}\,u\) and \(\frac{\partial }{\partial x}\,f\) do not always match on the boundaries. Thereupon, he further conjectures and verifies that the (path-wise) correction term, in our notation, takes the form

$$\begin{aligned} \left[ f_x(t,b)-u_x(t,b-)\right] (x-b) \end{aligned}$$

for the upper boundary, and

$$\begin{aligned} \left[ f_x(t,a)-u_x(t,a+)\right] (x-a) \end{aligned}$$

for the lower boundary. Note that these are linear in the state variable x, they vanish at the respective boundary, and they have a first-order derivative (with respect to x) equaling the difference between \(\frac{\partial }{\partial x}\,u\) and \(\frac{\partial }{\partial x}\,f\).

Nevertheless, this correction is only true for the single-barrier problem. To simultaneously account for the boundary effects at a and b, this form should not be linear. Indeed, a polynomial p(x) that satisfies \(p(a)=p(b)=0\) and also meets the two derivative differences at a and b (i.e., a total of four conditions) is at least cubic. Therefore, we consider

$$\begin{aligned} f^{1} (t,x) = & \frac{1}{{(b - a)^{2} }}{\mkern 1mu} {\mathcal{D}}_{b} (t)(x - a)^{2} (x - b) \\ & + \frac{1}{{(b - a)^{2} }}{\mkern 1mu} {\mathcal{D}}_{a} (t)(x - a)(x - b)^{2} , \\ \end{aligned}$$
(2.19)

where

$$\begin{aligned} {\mathcal {D}}_b(t)=f_x(t,b)-u_x(t,b-)\quad \textrm{and}\quad {\mathcal {D}}_a(t)=f_x(t,a)-u_x(t,a+) \end{aligned}$$
(2.20)

measure the first-order non-smoothness upon the boundaries b and a, respectively. Indeed, from (2.8) and (2.12) we have the general expression

$${\mathcal{D}}_{x} (t) = e^{{ - \alpha {\kern 1pt} t}} {\mkern 1mu} [\theta {\mkern 1mu} e^{{\theta {\kern 1pt} x}} - A(\alpha ,\theta ){\mkern 1mu} \beta _{{1,\alpha }} {\mkern 1mu} e^{{\theta {\kern 1pt} a + \beta _{{1,\alpha }} \left( {x - a} \right)}} - B(\alpha ,\theta ){\mkern 1mu} \beta _{{2,\alpha }} {\mkern 1mu} e^{{\theta {\kern 1pt} b + \beta _{{2,\alpha }} \left( {x - b} \right)}} ],$$
(2.21)

which is bounded for all \(x\in [a,b]\).

In this way,

$$\begin{aligned} f^1(t,a)=0=f^1(t,b),\quad f^1_x(t,a)={\mathcal {D}}_a(t),\quad \text{ and }\quad f^1_x(t,b)={\mathcal {D}}_b(t). \end{aligned}$$

As a result, if we let

$$\begin{aligned} f^0(t,x)=f(t,x)-f^1(t,x), \end{aligned}$$
(2.22)

then \(u(t,a)=f(t,a)=f^0(t,a)\), \(u(t,b)=f(t,b)=f^0(t,b)\), and

$$\begin{aligned} u_x(t,a)=f^0_x(t,a),\qquad u_x(t,b)=f^0_x(t,b). \end{aligned}$$
(2.23)

In other words, if we alter the terminal function from f to such a \(f^0\), then the continuity holds upon the boundaries in the PDE problem (2.10), as does the first-order smoothness.

In the next section, we will show that

$$\begin{aligned} E\left[ f^0(\tau _m\varDelta t,X_{\tau _m\varDelta t})\right] =u(0,0)+o\left( 1/\sqrt{m}\,\right) =E\left[ f^0(\tau _\infty ,X_{\tau _\infty })\right] +o\left( 1/\sqrt{m}\,\right) . \end{aligned}$$

Namely, there is no need (for order \(1/\sqrt{m}\)) for correction when smooth pasting holds upon the boundaries. This is exactly the core of Keener’s approach to extract the smoothing part that needs no correction from a general terminal function f via a decomposition like (2.22). Since clearly \(Ef=Ef^0+Ef^1\), the remaining task is to estimate \(E\left[ f^1(\tau _m\varDelta t,X_{\tau _m\varDelta t})\right]\) so that

$$\begin{aligned}&E\left[ f(\tau _m\varDelta t,X_{\tau _m\varDelta t})\right] \nonumber \\&\quad =u(0,0)+\lim \limits _{m\rightarrow \infty }E\left[ f^{1}(\tau _m\varDelta t, X_{\tau _m\varDelta t})\right] +o\left( 1/\sqrt{m}\,\right) . \end{aligned}$$
(2.24)

Thus, the improved convergence rate established for \(f^0\) still succeeds. More importantly, the limit of \(Ef^1\) is precisely the key we seek for correction to speed convergence in (2.18).

As a final remark, recall that \(\varDelta t=T/m\); hence \(m\rightarrow \infty\) is equivalent to \(\varDelta t\rightarrow 0\) here. To make our derivation more straightforward, we will simply take \(\varDelta t\) to describe the rate of convergence as m goes to infinity. Hence, for example, we write \(o(\sqrt{\varDelta t})\) rather than \(o(1/\sqrt{m})\) hereafter. The same logic applies for other situations.

2.5 Relevant random walks and the overshoot

Now, we define a discrete version of the continuous-time model (2.1), given that the monitoring frequency is m. We do so to facilitate the upcoming discussion. Let \(X_{m,0}=0\), and for \(n\in {\mathbb {N}}\), define

$$\begin{aligned} X_{m,n}=\sum _{k=1}^nV_{m,k}\quad \text{ with }\quad V_{m,k}=\left( r-\sigma ^2/2\right) \varDelta t +\sigma \sqrt{\varDelta t}\,Z_k^{(m)}. \end{aligned}$$
(2.25)

Here, for each m, the \(Z^{(m)}_k\)’s are independently and identically distributed (i.i.d.) standard normal random variables. Clearly, the random walk \(X_{m,n}\) follows the same law as \(X_{n\varDelta t}\); thus

$$\begin{aligned} \tau _m=\inf \left\{ n\in {\mathbb {N}}:\;X_{m,n}\notin \left( a,b\right) \right\} . \end{aligned}$$
(2.26)

Also, we will symbolize the filtration associated with \(X_{m,n}\) by \({\mathcal {F}}^{(m)}_n\).

Next, the associated overshoot can be defined by

$$\begin{aligned} \max \left\{ a-X_{m,\tau _m},\;X_{m,\tau _m}-b\right\} . \end{aligned}$$

To study its asymptotic property as \(m\rightarrow \infty\), we often rescale the random walk as \(X_{m,n}=\sigma \sqrt{\varDelta t}\,W_{m,n}\), where

$$\begin{aligned} W_{m,n}=\sum _{k=1}^nZ_{m,k}\quad \text{ with } \quad Z_{m,k}=\mu _m+Z_k^{(m)}\; \text {and}\; \mu _m=\frac{r-\sigma ^2/2}{\sigma }\,\sqrt{\varDelta t}. \end{aligned}$$
(2.27)

Thus, the overshoot is expressed as \(\max \left\{ \sigma \sqrt{\varDelta t}\,R_m(a),\;\sigma \sqrt{\varDelta t}\,R_m(b)\right\}\), in which

$$\begin{aligned} R_m(a):=a_m-W_{m,\tau _m}\quad \text{ and }\quad R_m(b):=W_{m,\tau _m}-b_m \end{aligned}$$
(2.28)

with \(a_m=\frac{a}{\sigma \sqrt{\varDelta t}}\) and \(b_m=\frac{b}{\sigma \sqrt{\varDelta t}}\).

The relevance of (2.28) to our problem can be observed from (2.19) via, for example, an expansion such as

$$\begin{aligned}&\left( \frac{x-a}{b-a}\right) ^2(x-b) \\&\quad =\left( \frac{x-b+b-a}{b-a}\right) ^2(x-b)\cdot 1_{\{x\ge b\}}+\left( \frac{x-a}{b-a}\right) ^2(x-a+a-b)\cdot 1_{\{x\le a\}}\\&\quad =\left\{ (x-b)+\frac{2}{b-a}\left( x-b\right) ^2+\frac{1}{(b-a)^2}\left( x-b\right) ^3\right\} \cdot 1_{\{x\ge b\}} \\&\qquad -\left\{ \frac{1}{b-a}\left( a-x\right) ^2+\frac{1}{(b-a)^2}\left( a-x\right) ^3\right\} \cdot 1_{\{x\le a\}} \end{aligned}$$

provided that \(x\notin (a,b)\). Thus, applying the same technique to both summands of \(f^1\), we finally obtain

$$\begin{aligned}&f^1(\tau _m\varDelta t, X_{m,\tau _m}) \nonumber \\&\quad =\Big \{{\mathcal {D}}_b(\tau _m\varDelta t)\left( X_{m,\tau _m}-b\right) +\frac{1}{b-a} \left[ 2{\mathcal {D}}_b(\tau _m\varDelta t)+{\mathcal {D}}_a(\tau _m\varDelta t)\right] \left( X_{m,\tau _m}-b\right) ^2 \nonumber \\&\qquad +\frac{1}{(b-a)^2}\left[ {\mathcal {D}}_b(\tau _m\varDelta t)+{\mathcal {D}}_a(\tau _m\varDelta t)\right] \left( X_{m,\tau _m}-b\right) ^3\Big \}\cdot 1_{\{X_{m,\tau _m}\ge b\}} \nonumber \\&\qquad -\Big \{{\mathcal {D}}_a(\tau _m\varDelta t)\left( a-X_{m,\tau _m}\right) +\frac{1}{b-a} \left[ 2{\mathcal {D}}_a(\tau _m\varDelta t)+{\mathcal {D}}_b(\tau _m\varDelta t)\right] \left( a-X_{m,\tau _m}\right) ^2 \nonumber \\&\qquad +\frac{1}{(b-a)^2}\left[ {\mathcal {D}}_a(\tau _m\varDelta t)+{\mathcal {D}}_b(\tau _m\varDelta t)\right] \left( a-X_{m,\tau _m}\right) ^3\Big \}\cdot 1_{\{X_{m,\tau _m}\le a\}} \nonumber \\&\quad =\Big \{{\mathcal {D}}_b(\tau _m\varDelta t)\left( \sigma \sqrt{\varDelta t}\,R_m(b)\right) \nonumber \\&\qquad +\frac{1}{b-a}\left[ 2{\mathcal {D}}_b(\tau _m\varDelta t)+{\mathcal {D}}_a(\tau _m\varDelta t)\right] \left( \sigma \sqrt{\varDelta t}\,R_m(b)\right) ^2 \nonumber \\&\qquad +\frac{1}{(b-a)^2}\left[ {\mathcal {D}}_b(\tau _m\varDelta t)+{\mathcal {D}}_a(\tau _m\varDelta t)\right] \left( \sigma \sqrt{\varDelta t}\,R_m(b)\right) ^3\Big \}\cdot 1_{\{W_{m,\tau _m}\ge b_m\}} \nonumber \\&\qquad -\Big \{{\mathcal {D}}_a(\tau _m\varDelta t)\left( \sigma \sqrt{\varDelta t}\,R_m(a)\right) \nonumber \\&\qquad +\frac{1}{b-a}\left[ 2{\mathcal {D}}_a(\tau _m\varDelta t)+{\mathcal {D}}_b(\tau _m\varDelta t)\right] \left( \sigma \sqrt{\varDelta t}\,R_m(a)\right) ^2 \nonumber \\&\qquad +\frac{1}{(b-a)^2}\left[ {\mathcal {D}}_a(\tau _m\varDelta t)+{\mathcal {D}}_b(\tau _m\varDelta t)\right] \left( \sigma \sqrt{\varDelta t}\,R_m(a)\right) ^3\Big \}\cdot 1_{\{W_{m,\tau _m}\le a_m\}}. \end{aligned}$$
(2.29)

Thus, to estimate \(Ef^1\), we must study the conditional moments of the overshoot; in particular,

$$\begin{aligned} E[R_m(b);\,W_{m,\tau _m}\ge b_m]\quad \text{ and }\quad E[R_m(a);\,W_{m,\tau _m}\le a_m]. \end{aligned}$$

Coincidentally, the situation here fits the framework considered in Siegmund (1979), where the incremental mean goes to zero and the boundary level goes to infinity synchronously. More precisely, as \(m\rightarrow \infty\), \(\mu _m\rightarrow 0\), \(|a_m|\rightarrow \infty\), and \(b_m\rightarrow \infty\), whereas both \(\mu _m|a_m|\) and \(\mu _mb_m\) remain constant. Therefore, we will choose to extend Siegmund’s results to cover our conditional version.

Finally, we emphasize that the expansion (2.29) from (2.19) is crucial for achieving our final correction formula. This is in fact not straightforward and requires some effort to figure it out. To the best of our knowledge, such a design is new in the related fields.

3 Main results

We here establish our correction result, under the two presumptions (2.6) and (2.9). According to the discussion in the previous section, there are actually three main steps. First, we will show there is no need to call for correction with regard to the \(f^0\) part in the decomposition (2.22). Second, we determine the estimation of the expected value of \(f^1\), and then combine these two results to achieve our final continuity correction formula for the pricing of discrete options with double barriers.

3.1 Estimation of \(Ef^0\)

First, recall the notation related to (2.12) and (2.22). Let

$$\begin{aligned} \bar{v}(t,x)=\left\{ \begin{array}{ll} u(t,x), &{} x\in (a,b); \\ f^0(t,x), &{} \text{ otherwise }, \end{array} \right. \end{aligned}$$

and

$$\begin{aligned} e_m(t,x)=\left\{ \begin{array}{lll} E\left[ \bar{v}\left( t+\varDelta t, x+V_{m,1}\right) \right] -\bar{v}(t,x), &{} x\in (a,b); \\ 0, &{} \text{ otherwise }. \end{array} \right. \end{aligned}$$
(3.1)

Then,

$$\begin{aligned}&E\left[ f^0(\tau _m\varDelta t, X_{m,\tau _m})\right] = E\left[ \bar{v}(\tau _m\varDelta t, X_{m,\tau _m})\right] \\&\quad =\bar{v}(0,0)+E\left[ \sum _{k=0}^{\tau _m-1}\bar{v}\left( (k+1)\varDelta t, X_{m,k+1}\right) -\bar{v}\left( k\varDelta t, X_{m,k}\right) \right] \\&\quad =u(0,0)+E\left[ \sum _{k=0}^{\tau _m-1}E\left[ \bar{v}\left( (k+1)\varDelta t, X_{m,k+1}\right) -\bar{v}\left( k\varDelta t, X_{m,k}\right) \big |\,{\mathcal {F}}^{(m)}_n\right] \right] . \end{aligned}$$

By the i.i.d. property of \(V_{m,k}\), the last equation becomes

$$\begin{aligned} E\left[ f^0(\tau _m\varDelta t, X_{m,\tau _m})\right] =u(0,0)+E\left[ \sum _{k=0}^{\tau _m-1}e_m \left( k\varDelta t, X_{m,k}\right) \right] . \end{aligned}$$
(3.2)

As a result, the key to estimating \(Ef^0\) lies in the sum of these \(e_m\) terms.

Theorem 3.1

Under the setting specified in Sect. 2, we have

$$\begin{aligned} E\left[ f^0(\tau _m\varDelta t, X_{m,\tau _m})\right] =E\left[ f(\tau _\infty , X_{\tau _\infty })\right] +o(\sqrt{\varDelta t}) \end{aligned}$$

as \(m\rightarrow \infty\).

Proof

First, by Lemma A.1 in the appendix,

$$\begin{aligned} e_m(k\varDelta t,X_{m,k})=O\left\{ (\varDelta t)^2\right\} +O\left\{ \frac{(\varDelta t)^2}{\sigma ^2\varDelta t +\min \left\{ (b-X_{m,k})^2, (a-X_{m,k})^2\right\} }\right\} , \end{aligned}$$

and note from (2.27) and (2.28) that the following three events are equivalent:

$$\begin{aligned} \left\{ X_{m,n}\notin (a,b)\right\} =\left\{ W_{m,n}\notin (a_m,b_m)\right\} =\left\{ W_{n}\notin (a_m(n),b_m(n))\right\} , \end{aligned}$$
(3.3)

whereFootnote 5

$$\begin{aligned} W_n=\sum _{k=1}^nZ_k^{(m)},\quad a_m(n)=a_m-n\mu _m,\quad \text{ and }\quad b_m(n)=b_m-n\mu _m. \end{aligned}$$
(3.4)

Therefore,

$$\begin{aligned}&\min \left\{ (b-X_{m,k})^2,\;(a-X_{m,k})^2\right\} \\&\quad =\sigma ^2\varDelta t\,\times \,\min \left\{ (b_m(k)-W_{k})^2,\;(a_m(k)-W_{k})^2\right\} , \end{aligned}$$

and hence

$$\begin{aligned} e_m(k\varDelta t,X_{m,k})=O\left\{ (\varDelta t)^2\right\} +O\left\{ \frac{\varDelta t}{1+\min \left\{ (b_m(k)-W_{k})^2,\,(a_m(k)-W_{k})^2\right\} }\right\} . \end{aligned}$$

Then, together with Lemmas A.2 and (A.10) in Appendix A, (3.2) becomes

$$\begin{aligned}&E\left[ f^0(\tau _m\varDelta t, X_{m,\tau _m})\right] \\&\quad =u(0,0)\\&\quad +O\left( \varDelta t\right) \cdot \left\{ E(\tau _m\varDelta t)+ E\left[ \sum _{k=0}^{\tau _m-1}\frac{1}{1+\min \left\{ (b_m(k)-W_k)^2,(a_m(k)-W_k)^2\right\} }\right] \right\} \\&\quad =E\left[ f(\tau _\infty , X_{\tau _\infty })\right] +O\left( \varDelta t\log ({\varDelta t}^{-1})\right) . \end{aligned}$$

Finally, note that

$$\begin{aligned} \lim _{m\rightarrow \infty }\frac{\varDelta t\log ({\varDelta t}^{-1})}{\sqrt{\varDelta t}}=\lim _{y\rightarrow \infty }\frac{2\log (y)}{y}=0 \qquad (\text{ letting }\, y=(\varDelta t)^{-\frac{1}{2}}) \end{aligned}$$

by L’Hôpital’s rule. The desired result thus follows. \(\square\)

3.2 Estimation of \(Ef^1\)

Recalling the notation in § 2.5, consider the following two hitting times of a single boundary,

$$\begin{aligned} \hat{\tau }_m=\inf \left\{ n\in {\mathbb {N}}:\;X_{m,n}\ge b\right\} \quad \text{ and }\quad \check{\tau }_m=\inf \left\{ n\in {\mathbb {N}}:\;X_{m,n}\le a\right\} , \end{aligned}$$
(3.5)

and define their associated overshoots as

$$\begin{aligned} \hat{R}_m(b)=W_{m,\hat{\tau }_m}-b_m\quad \text{ and }\quad \check{R}_m(a)=a_m-W_{m,\check{\tau }_m}, \end{aligned}$$
(3.6)

respectively. It is known that for any \(p\in {\mathbb {N}}\) both

$$\begin{aligned} \lim _{m\rightarrow \infty }E\left[ \hat{R}^p_m(b);\,\hat{\tau }_m<\infty \right]<\infty \;\;\text{ and }\;\; \lim _{m\rightarrow \infty }E\left[ \check{R}^p_m(a);\,\check{\tau }_m<\infty \right] <\infty ; \end{aligned}$$
(3.7)

see Siegmund (1985, pp. 186, 215) for example. With this preliminary knowledge, we have

Theorem 3.2

Under the setting specified in Sect. 2, we have as \(m\rightarrow \infty\)

$$\begin{aligned}&E\left[ f^1(\tau _m\varDelta t, X_{m,\tau _m})\right] \\&\quad =\sigma \sqrt{\varDelta t}\,\rho \cdot \Big (q_b\,E\left[ {\mathcal {D}}_b(\tau _{\infty })\right] -q_a\,E\left[ {\mathcal {D}}_a(\tau _{\infty })\right] \Big )+o(\sqrt{\varDelta t}), \end{aligned}$$

where \(q_a\) and \(q_b\) are given in (2.16), and \(\rho =-\zeta (1/2)/\sqrt{2\pi }\approx 0.5826\) which is defined in (A.13).

Proof

Because

$$\begin{aligned} \left\{ R^p_m(b)>x,\,W_{m,\tau _m}\ge b_m\right\}&\subseteq \left\{ \hat{R}^p_m(b)>x,\,\hat{\tau }_m<\infty \right\} \quad \text{ and }\\ \left\{ R^p_m(a)>x,\,W_{m,\tau _m}\le a_m\right]&\subseteq \left\{ \check{R}^p_m(a)>x,\,\check{\tau }_m<\infty \right\} , \end{aligned}$$

the fact (3.7) also implies

$$\begin{aligned} \lim _{m\rightarrow \infty }E\left[ R^p_m(b);\,W_{m,\tau _m}\ge b_m\right]<\infty \quad \text{ and } \quad \lim _{m\rightarrow \infty }E\left[ R^p_m(a);\,W_{m,\tau _m}\le a_m\right] <\infty . \end{aligned}$$

That is, both

$$\begin{aligned} E\left[ (X_{m,\tau _m}-b)^p\cdot 1_{\left\{ X_{m,\tau _m}\ge b\right\} }\right]&=O\left\{ (\sqrt{\varDelta t})^p\right\} \quad \text{ and }\\ E\left[ (a-X_{m,\tau _m})^p\cdot 1_{\left\{ X_{m,\tau _m}\le a\right\} }\right]&=O\left\{ (\sqrt{\varDelta t})^p\right\} . \end{aligned}$$

Also, it is easy to see from (2.21) that \(D_b(t)\) and \(D_a(t)\) are bounded uniformly in t; thus, we can conclude from (2.29) that

$$\begin{aligned}&{E\left[ f^1(\tau _m\varDelta t, X_{m,\tau _m})\right] } \\&\quad =\sigma \sqrt{\varDelta t}\cdot E\left[ {\mathcal {D}}_b\left( \tau _m\varDelta t\right) R_m(b);\,W_{m,\tau _m}\ge b_m\right] \\&\quad -\sigma \sqrt{\varDelta t}\cdot E\left[ {\mathcal {D}}_a\left( \tau _m\varDelta t\right) R_m(a);\,W_{m,\tau _m}\le a_m\right] +o(\sqrt{\varDelta t}). \end{aligned}$$

Then, the result follows from Lemma A.4. \(\square\)

Remark 3.1

It is possible to extend Theorems 3.1 and 3.2 to cases with time-dependent barriers; namely, consider \(a=a(t)=\log (L(t)/S_0)\) and \(b=b(t)=\log (H(t)/S_0)\) for some time-varying barriers L(t) and H(t) with \(L(0)<S_0<H(0)\). To this end, we must impose two further conditions on the barriers according to Keener (2013). First, a(t) and b(t) have bounded first derivatives. Second, for some \(\epsilon >0\),

$$\begin{aligned} b(t)-\left( r-\sigma ^2/2\right) t +\epsilon t \rightarrow -\infty \quad \text{ if }\, r-\sigma ^2/2>0 \end{aligned}$$

or

$$\begin{aligned} a(t)-\left( r-\sigma ^2/2\right) t -\epsilon t \rightarrow +\infty \quad \text{ if }\, r-\sigma ^2/2<0 \end{aligned}$$

such that the stopping times

$$\begin{aligned} \tau _m(t)=\inf \left\{ t=n\varDelta t:\;n\in {\mathbb {N}}\;\text{ and }\;X_{t}\notin \left( a(t),b(t)\right) \right\} ,\ \ m\ge 1, \end{aligned}$$

are uniformly integrable. However, the applicability of such a setting to real finance problems is still unknown or questionable.

3.3 Discrete option pricing

We are now ready to study the pricing of a discrete KIC. First, recall \(A(\alpha ,\theta )\) and \(B(\alpha ,\theta )\) in (2.13). Then, by Theorems 3.1 and 3.2, we have

Corollary 3.1

Under the setting specified in Sect. 2,

$$\begin{aligned} E\left[ e^{-\alpha \left( \tau _m\varDelta t\right) +\theta \,X_{m,\tau _m}}\right]&=A\left( \alpha ,\theta \right) e^{\left( a-\sigma \sqrt{\varDelta t}\,\rho \,\varphi _1(\alpha )\right) \left( \theta -\beta _{1,\alpha }\right) } \nonumber \\&\quad+B\left( \alpha ,\theta \right) e^{\left( b+\sigma \sqrt{\varDelta t}\,\rho \,\varphi _2(\alpha )\right) \left( \theta -\beta _{2,\alpha }\right) }+o(\sqrt{\varDelta t}) \end{aligned}$$
(3.8)

as \(m\rightarrow \infty\), in which

$$\begin{aligned} \varphi _1(\alpha )=\left( q_a\,e^{a\,\beta _{1,\alpha }}-q_b\,e^{b\,\beta _{1,\alpha }}\right) \phi (\alpha ) \end{aligned}$$

and

$$\begin{aligned} \varphi _2(\alpha )=\left( q_b\,e^{b\,\beta _{2,\alpha }}-q_a\,e^{a\,\beta _{2,\alpha }}\right) \phi (\alpha ) \end{aligned}$$

with \(\phi (\alpha )\) being defined in (2.15).

Proof

What remains is to work out \(E\left[ {\mathcal {D}}_a(\tau _{\infty })\right]\) and \(E\left[ {\mathcal {D}}_b(\tau _{\infty })\right]\) explicitly. The details are deferred to the appendix. \(\square\)

Remark 3.2

Note that from (2.3) we have \(\beta _{1,\alpha }+\beta _{2,\alpha }=-2\gamma\) for any \(\alpha \ge 0\). The hitting probabilities in (2.16) hence can be rewritten as

$$\begin{aligned} q_a=\frac{1-e^{-\left( \beta _{1,\alpha }+\beta _{2,\alpha }\right) b}}{1-e^{-\left( \beta _{1,\alpha } +\beta _{2,\alpha }\right) (b-a)}}\quad \text{ and }\quad q_b=\frac{1-e^{-\left( \beta _{1,\alpha } +\beta _{2,\alpha }\right) a}}{1-e^{\left( \beta _{1,\alpha }+\beta _{2,\alpha }\right) (b-a)}}. \end{aligned}$$

With the aid of this form, further computation (detailed but straightforward) leads to

$$\begin{aligned}&\lim _{a\rightarrow -\infty }\varphi _1(\alpha )=\infty ,\quad \lim _{a\rightarrow -\infty }\varphi _2(\alpha )=\lambda _b, \\&\quad \lim _{b\rightarrow \infty }\varphi _1(\alpha )=\lambda _a,\quad \text{ and }\quad \lim _{b\rightarrow \infty }\varphi _2(\alpha )=\infty , \end{aligned}$$

in which \(\lambda _a\) and \(\lambda _b\) are defined in (2.17). We then obtain the next two special cases of (3.8):

$$\begin{aligned}&\lim _{a\rightarrow -\infty }E\left[ e^{-\alpha \left( \tau _m\varDelta t\right) +\theta \,X_{m,\tau _m}}\right] =E\left[ e^{-\alpha \left( \hat{\tau }_m\varDelta t\right) +\theta \,X_{m,\hat{\tau }_m}};\,\hat{\tau }_m<\infty \right] \\&\quad =e^{\left( b+\sigma \sqrt{\varDelta t}\,\rho \,\lambda _b\right) \left( \theta -\beta _{2,\alpha }\right) }+o(\sqrt{\varDelta t}); \\&\lim _{b\rightarrow \infty }E\left[ e^{-\alpha \left( \tau _m\varDelta t\right) +\theta \,X_{m,\tau _m}}\right] =E\left[ e^{-\alpha \left( \check{\tau }_m\varDelta t\right) +\theta \,X_{m,\check{\tau }_m}};\,\check{\tau }_m<\infty \right] \\&\quad =e^{\left( a-\sigma \sqrt{\varDelta t}\,\rho \,\lambda _a\right) \left( \theta -\beta _{1,\alpha }\right) }+o(\sqrt{\varDelta t}). \end{aligned}$$

Note that when the drift is positive (negative), \(\lambda _b=1\) \((\lambda _a=1)\) indeed, and the above results are then in line with the literature on single barriers; see (Kou, 2008) for instance. Namely, our achievement here coincides with the single-barrier case.

Next, on the basis of (2.7), we re-symbolize the KIC price of a discrete version by

$$\begin{aligned} KIC _m(T,k)=Xe^{-rT}E\left[ \left( S_T/X-e^{-k}\right) ^{+}\,1_{\{\tau _{m}\varDelta t\le T\}}\right] \end{aligned}$$

and the associated Laplace transform by

$$\begin{aligned} \widehat{ KIC }_m(\alpha ,\theta )=\int _0^\infty \int _{-\infty }^{\infty }\, e^{-\alpha \,T-\theta \,k}\, KIC _m(T,k)\,dk\,dT. \end{aligned}$$

Recall that here \(k=\log (X/K)\) is the transformed strike price. Then, we have

Corollary 3.2

For \(0<\theta <\left( \beta _{2,r+\alpha }-1\right)\) and \(\alpha >\max \left( G(\theta +1)-r,\,0\right)\),

$$\begin{aligned} \widehat{ KIC }_m(\alpha ,\theta )&=\frac{X^{-\theta }}{\theta (\theta +1)} \times \frac{S_0^{\theta +1}}{r+\alpha -G(\theta +1)} \nonumber \\&\times \bigg \{A\left( r+\alpha ,\theta +1\right) e^{\left( a-\sigma \sqrt{\varDelta t}\, \rho \,\varphi _1(r+\alpha )\right) \left( \theta +1-\beta _{1,r+\alpha }\right) } \nonumber \\&+B\left( r+\alpha ,\theta +1\right) e^{\left( b+\sigma \sqrt{\varDelta t}\,\rho \, \varphi _2(r+\alpha )\right) \left( \theta +1-\beta _{2,r+\alpha }\right) }\bigg \} \nonumber \\&+o(\sqrt{\varDelta t}) \end{aligned}$$
(3.9)

as \(m\rightarrow \infty\).

Proof

See the appendix. \(\square\)

As mentioned in § 2.2, we will adopt the algorithm proposed by Cai and Kou (2012) to invert the Laplace transforms in Corollary 3.2. Generally speaking, their method is a combination of a two-sided Euler inversion and an Euler transformation for alternating series. The inversion parameters we choose are \(A_1=A_2=28.3\), \(n_1=38\), and \(n_2=11\) as suggested in Cai et al. (2009). On the other hand, the criteria for parameters \(\alpha\) and \(\theta\) in the corollary actually can be reduced to \(0<\theta <\beta _{2,r+\alpha }-1\), provided that \(r>0\). Hence, when implementing the inversion algorithm, we simply set the scaling factor to

$$\begin{aligned} X=K\,e^{A_2/(2k_c)}\quad \text{ with } \quad k_c=\frac{-\left( r+\sigma ^2/2\right) +\sqrt{\left( r+\sigma ^2/2\right) ^2+2\,\sigma ^2\,A_1/(2T)}}{c\,\sigma ^2} \end{aligned}$$
(3.10)

for some \(c>2\), where K is the original strike price and T is the time to maturity.

4 Numerical examples

In this section, we will present examples to show the performance of our approximation formulas (CC). The benchmark values are obtained using Monte-Carlo (MC) simulation. One million trials were conducted in each simulation. We also report values directly from continuous formulas without correction (NC). Comparisons were made in terms of relative error (R.E.) defined by (“X”−MC)/MC, where X is NC or CC.

Table 1 Approximation performance: \(E\left[ e^{-\alpha \left( \tau _m\varDelta t\right) +\theta \,X_{m,\tau _m}}\right]\) with positive drift

We first discuss the performance of (3.8). The outcomes are reported in Tables 1 and 2 for different considerations of incremental mean. Generally we see that our correction method does improve the rate of convergence, and that the improvement is quite stable across different boundary levels, regardless of whether the drift is positive or negative. Another important finding is that the signs of \(\varphi _1\) and \(\varphi _2\) in our correction formula are not necessarily positive. This means we do not always shift away each individual boundary to achieve continuity correction. In Table 2, for instance, when the lower boundary is \(a=-0.0175\) and the upper boundary is \(b=0.0225\), we must lower both boundary levels. By contrast, when \(a=-0.0225\), the upper barrier level is raised.

Additionally, the total amount adjusted at this point could be less than that in the single-boundary case. As hinted in Remark 3.2, the total adjustment is

$$\begin{aligned} {\mathcal {V}}:=\sigma \sqrt{\varDelta t}\,\rho \end{aligned}$$
(4.1)

when there is only one barrier. Here, interestingly, it is possible that \(|\varphi _1|+|\varphi _2|<1\) in the case of two simultaneous barriers. For example, in Table 1 with \(a=-0.0225\) and \(b=0.0175\), the total adjustment is just \(0.2781{\mathcal {V}}\). In other words, the shifted amount is only 28% of the single-barrier correction. However, the window length between barriers does widen after correction:

$$\begin{aligned}&[b+\sigma \sqrt{\varDelta t}\,\rho \,\varphi _2(\alpha )]-[a-\sigma \sqrt{\varDelta t}\,\rho \,\varphi _1(\alpha )] \\&\quad =(b-a)+\sigma \sqrt{\varDelta t}\,\rho \left[ \varphi _1(\alpha )+\varphi _2(\alpha )\right] \;>(b-a) \end{aligned}$$

thanks to the positiveness of \(\varphi _1+\varphi _2\). To see this, note that

$$\begin{aligned} \varphi _1(\alpha )+\varphi _2(\alpha )=\phi (\alpha )\left[ q_a\,(e^{a\,\beta _{1,\alpha }}-e^{a\,\beta _{2,\alpha }}) +q_b\,(e^{b\,\beta _{2,\alpha }}-e^{b\,\beta _{1,\alpha }})\right] >0. \end{aligned}$$

The outcome of the enlarged barrier window is consistent with the usual intuition that with discrete monitoring there is a smaller boundary-crossing probability than that in the continuous case. However, note that this argument is valid only given a fixed duration (time window).

Table 2 Approximation performance: \(E\left[ e^{-\alpha \left( \tau _m\varDelta t\right) +\theta \,X_{m,\tau _m}}\right]\) with negative drift
Table 3 Approximation performance: \(KIC _m\) with positive drift

On the other hand, for most cases in Tables 1 and 2, we see that the correction applied to each boundary is generally uneven. The only exception is when the lower and upper boundaries are equally far away (i.e., \(|a|=b\)). In this case, we see that \(\varphi _1=\varphi _2\). In fact, when taking \(a=-b\), we have

$$\begin{aligned} \varphi _1(\alpha )=\varphi _2(\alpha )=\frac{e^{b\beta _{2,\alpha }}-e^{b\beta _{1,\alpha }}}{e^{b\beta _{2,\alpha }} +e^{b\beta _{1,\alpha }}}:=\varphi (\alpha ). \end{aligned}$$
(4.2)

Then, it is easy to check \(0<\varphi (\alpha )<1\), and for each given \(\alpha\), \(\varphi (\alpha )\) increases in b. All observations from the tables with \(|a|=b\) meet these expectations. In fact, further graphical analysis shows that this seems to be the only situation for which \(\varphi _1(\alpha )=\varphi _2(\alpha )\). This special case also constitutes evidence against a common conjecture, that is, whether the correction tends to be greater for the boundary that is more likely to be reached. To illustrate, consider Table 1 where the drift is positive. In this case, we anticipate that the upper boundary \(b=0.0075\) is more likely to be hit earlier than the lower boundary \(a=-0.0075\). However, we see that the correction to each boundary requires the same level (\(\varphi _1=\varphi _2=0.0533\)). So far, these findings are specific to the case with two simultaneous boundaries.

Table 4 Approximation performance: \(KIC _m\) with negative drift

Now we consider the performance for the proposed KIC pricing. Note that in Tables 1 and 2 we consider only one parameter setting for \(\alpha\) and \(\theta\). Since the option price is computed via Laplace inversion, the oncoming analysis can be viewed as an integrated evaluation of the approximation performance over a wide range of parameter levels. The results are presented in Tables 3 and 4 for positive and negative drift, respectively. For simplicity, we consider only those cases in which the two barriers are of equal distance, roughly corresponding to the situation behind (4.2).

As evidenced by the tables, we see that pricing errors do significantly improve with our approximation method, in comparison with that without correction. In some situations, we see relative errors within 1% even with monthly monitoring. (Assuming 250 trading days in a year, a \(\varDelta t\) of 1/10 is roughly equivalent to monthly monitoring.) For instance, in Table 4 with \(L=82\) and \(H=118\), the error is a mere 0.51% with \(\varDelta t=1/10\) for at-the-money case (Panel B). Finally, note that we also report the outcomes when one simply extends the correction idea from single barrier to double barrier. Namely, both barriers are shifted away by the same amount of \({\mathcal {V}}\) (i.e., \(\varphi _1=\varphi _2=1\)). This naive approach echoes the boundary correction result of Gobet and Menozzi (2010). From the tables, we clearly see that this is not an accurate continuity correction. Sometimes, it even does worse than that without correction; see Panel B of Table 3 with \(L=92\) and \(H=108\). That is, continuity correction should involve not only the overshoot expectation (represented by \({\mathcal {V}}\)) but also the boundary hitting probability (represented by \(\varphi _1\) and \(\varphi _2\)). This latter analysis again confirms the merit and necessity of our efforts here.

Finally, recall from Sect. 1 that the improved rate of convergence \(o(\sqrt{\varDelta t})\) is of the first order.Footnote 6 Hence the numerical outcomes obtained up to this point seem to reveal that the idea of first-order correction, in terms of an enlarged boundary window, takes effect mainly for the concern of boundary-hitting probabilities. This partly explains why the overall improvement brought by our correction method in Tables 1 and 2 is not as significant as that in Tables 3 and 4. The moment generating functions studied in the former two tables actually do not presume a duration; sooner or later, the boundary will be reached, regardless of the monitoring frequency. Since correction is not always outward, it is possible to see the boundary-hitting event occur earlier after the correction, especially when the boundary window is corrected in a direction opposite to the mean direction. For example, in Table 2, where a negative drift is considered, when \(a=-0.0075\) and \(b=0.0025\), the correction raises both boundaries, making it more likely to hit the lower boundary. As \(\varDelta t=1/10\), the corrected lower boundary even becomes positive, which further overestimates the stopped process level. Such an observation mainly accounts for why we see the correction result even underperforming that without correction. Nevertheless, this odd phenomenon ultimately disappears when increasing the monitoring frequency (i.e., as \(\varDelta t\rightarrow 0\)), consistent with the theorem.

Note that these findings also suggest the need for higher-order correction, in order to improve the weak convergence result as in (2.18) for practical applications. To account for the delay of the first hitting time with discrete monitoring, perhaps a rate of \(o(\varDelta t)\) is needed, requiring the moment calculation of the overshoot to the second. Still, at what order would the improved rate of convergence be satisfactory?—this remains an open question. By contrast, neither the hitting time nor the stopped price level matters with respect to option payoff. Only the barrier hitting probability plays a role in option valuation. This may explain why simple first-order correction suffices for the pricing of discrete options, per the existing literature.

5 Concluding remarks

Continuity correction in option pricing has been a classic challenge since the pioneer work of Broadie et al. (1997). The general idea is to approximate the price of a path-dependent option whose underlying process is monitored discretely by its associated pricing formula, assuming continuous monitoring, with adjustments to relevant barriers. In particular, for a single barrier option, the original barrier level is shifted away by an amount \({\mathcal {V}}\) in (4.1) before directly applying the continuous-time pricing formula. The quantity \({\mathcal {V}}\) is associated with the expectation of a pertinent overshoot. However, the correction for a double barrier option is a longstanding issue which is unsolved even for the classical Black–Scholes model. Here we draw to a close in part by extending the PDE approach of Keener (2013).

It is generally conjectured that for the double-barrier case, the correction method is the same as for the single-barrier case: one simply shifts away the amount \({\mathcal {V}}\) simultaneously to both barriers. Now, more than twenty years later, we show that this is not true. Continuity correction involves not only the overshoot expectation but also the boundary hitting probability. As a result, the correction amounts and even directionsFootnote 7 are highly dependent on the distances to each barrier from the current level. This constitutes the main contribution of our work to the literature. What remains unchanged is that the region within the barriers does widen after correction, to account for a lower boundary-crossing probability with discrete monitoring during a given period of time, as in the single-barrier case.

One clear limitation of the current work is the sophistication of the model. It is widely accepted that the Black–Scholes model is too simple to address real price processes. Nevertheless, due to its tractability, the BS model remains popular both in academia and in the industry. More specifically, many other numerical techniques can be used to price a discrete option under the BS model. Thus one relevant question concerning our study is: What is the advantage of our correction method? A preliminary answer is that the correction result produced indicates how to adjust the analytical formula not only between discrete and continuous monitoring, but also between different discrete monitoring frequencies. The latter is of particular use when the frequency adopted for payoff calculation differs from that for barrier monitoring. This is one popular design in a structured product, where the payoff is typically calculated on a low-frequency basis whereas the barrier monitoring is at a high frequency. Note that the payoff and the barrier can be imposed on different assets. Then, it is plausible to develop an efficient simulation procedure based on a low frequency but further weighted by the boundary crossing probability calculated at a high frequency. Our correction achievement may be helpful in such a scenario.

Other promising future studies include the following. First, how general is Keener’s PDE approach with respect to the model setting? In particular, does the same idea apply to the mixed-exponential jump diffusion model proposed by Cai and Kou (2011), given that the model is flexible enough to compete with all exponential Lévy models of finite activity while maintaining tractability? Second, is the approach adequate to other type of path-dependent options? We are particularly interested in the pricing of American options. As mentioned, the overshoot effect only matters when the boundary is not passed smoothly. However, this is the usual condition we impose to determine the price for an American option of finite maturity. One exception is the real option, where smooth pasting conditions are not always assumed. Thereby, it is likely that continuity correction will be of assistance in the valuation of real options.