On extension of the Markov chain approximation method for computing Feynman--Kac type expectations

TL;DR

This paper extends a Markov chain approximation method with Brownian bridge correction to efficiently compute Feynman--Kac expectations, especially for option pricing, demonstrating high convergence rates in numerical experiments.

Contribution

The paper introduces an extension of the Markov chain approximation method to compute Feynman--Kac expectations for diffusion processes, applicable to option pricing.

Findings

Convergence rate of $O(n^{-2})$ for smooth integrands.

Effective extension to path-dependent functionals.

Numerical experiments confirm high accuracy.

Abstract

An efficient discrete time and space Markov chain approximation employing a Brownian bridge correction for computing curvilinear boundary crossing probabilities for general diffusion processes was recently proposed in Liang and Borovkov (2021). One of the advantages of that method over alternative approaches is that it can be readily extended to computing expectations of path-dependent functionals over the event of the process trajectory staying between two curvilinear boundaries. In the present paper, we extend the scheme to compute expectations of the Feynman--Kac type that frequently appear in option pricing. To illustrate our approximation scheme, we apply it in three special cases. For sufficiently smooth integrands, numerical experiments suggest that the proposed approximation converges at the rate $O(n^{-2})$, where $n$ is the number of steps on the uniform time grid used

Equations

$$\varphi (X(T))1\{g^{-}(s)<X(s)<g^{+}(s),s\in [0,T]\}.$$

$$\mathbf{E}\bigl{[}\phi(X(T))/N(T);g^{-}(s)<X(s)<g^{+}(s),\,s\in[0,T]\bigr{]}.$$

$$\mathbf{E}\Bigl{[}e^{-\int_{0}^{T}V(X(u))\,du}\phi(X(T));g^{-}(s)<X(s)<g^{+}(s),\,s\in[0,T]\Bigr{]}.$$

$$X(t)=x_0+{\int }_0^t\mu (s,X(s))ds+W(t),t\ge 0,$$

$$|\mu (t,x)|+|{\partial }_t\mu (t,x)|+|{\partial }_x\mu (t,x)|+|{\partial }_{xx}\mu (t,x)|\le K_r,t\in [0,1],|x|\le r.$$

$$\mathcal{G}:=\Bigl{\{}(f^{-},f^{+}):f^{\pm}\in C,\,f^{-}(0)<x_{0}<f^{+}(0),\,\min_{0\leq t\leq 1}(f^{+}(t)-f^{-}(t))>0\Bigr{\}}$$

$$S(f^{-},f^{+}):=\{v\in C:f^{-}(t)<v(t)<f^{+}(t),t\in [0,1]\},(f^{-},f^{+})\in \mathcal{G}$$

$$Q:=\mathbf{E}\Bigl{(}e^{-\int_{0}^{1}V(X(s))\,ds}\phi(X(1));X\in G\Bigr{)},\quad G:=S(g^{-},g^{+}),$$

$$T_{s,t}f(x)=\mathbf{E}^{s,x}\Bigl{(}e^{-\int_{s}^{t}V(X(u))\,du}f(X(t));g^{-}(u)<X(u)<g^{+}(u),u\in[s,t]\Bigr{)}$$

$$Q=(T_{0,t_{n,1}}{T}_{t_{n,2},t_{n,3}}\cdots T_{t_{n,n-1},t_{n,n}}\varphi )(x_0),$$

$$t_{n,k}:=k/n,k=0,1,\dots ,n,$$

$$w_{n,k}:=\{\begin{array}{ll}\frac{(g_{n,k}^{+}-g_{n,k}^{-})/{\Delta }_n^{1/2+\delta }}{\lfloor \gamma (g_{n,k}^{+}-g_{n,k}^{-})/{\Delta }_n^{1/2+\delta }\rfloor },& 1\le k<n,\\ \frac{(g^{+}(1)-g^{-}(1))/{\Delta }_n}{\lfloor \gamma (g^{+}(1)-g^{-}(1))/{\Delta }_n\rfloor },& k=n,\end{array}$$

$$h_{n,k}:=\{\begin{array}{ll}{w}_{n,k}{\Delta }_n^{1/2+\delta },& 1\le k<n,\\ w_{n,n}{\Delta }_n,& k=n.\end{array}$$

$$E_{n,k}:=\{g_{n,k}^{+}-j{h}_{n,k}:j\in \mathbb{Z}\},k=1,\dots ,n.$$

$$E_{n,k}^G:=\{x\in E_{n,k}:g_{n,k}^{-}<x<g_{n,k}^{+}\},k=1,\dots ,n.$$

$${\mu }_{n,k}(x)$$

$${\sigma }_{n,k}^2(x)$$

$$p_{n,k}(x,y):=\phi (y|x+{\mu }_{n,k}(x),{\sigma }_{n,k}^2(x))h_{n,k},(x,y)\in E_{n,k-1}\times E_{n,k},$$

$$\pi_{n,k}(x,y)=1-\mbox{$\exp\{\frac{-2}{\Delta_{n}}(g_{n,k-1}^{+}-x)(g_{n,k}^{+}-y)\}$}-\mbox{$\exp\{\frac{-2}{\Delta_{n}}(g_{n,k-1}^{-}-x)(g_{n,k}^{-}-y)\}$}.$$

$$e_{n,k}(x,y):=\exp\{\mbox{$\frac{-\Delta_{n}}{2}$}(V(x)+V(y))\}.$$

$$S_{n,k}:=\big{[}p_{n,k}(x,y)\pi_{n,k}(x,y)e_{n,k}(x,y)\big{]}_{(x,y)\in E_{n,k-1}^{G}\times E_{n,k}^{G}}.$$

$$Q_n:=S_{n,1}{S}_{n,2}\cdots S_{n,n}{\varphi }_n,n\ge 1,$$

$${\tilde{X}}_n(t):=B_{n,k}^{{\xi }_{n,k-1},{\xi }_{n,k}}(t),t\in [t_{n,k-1},t_{n,k}],k=1,\dots ,n,$$

$$B_{n,k}^{x,y}(t):=B_{n,k}^{\circ }(t)+x+n(t-t_{n,k-1})(y-x),x,y\in \mathbb{R},$$

$$Z(t):={\int }_0^{t}V(X(s))ds,t\in [0,1].$$

$$Q=\mathbf{E}\big{(}\Psi(\mathbf{X}(1));X\in G\big{)},$$

$$Q_{n}=\mathbf{E}\bigl{(}\Psi(\mathbf{X}_{n}(1));\widetilde{X}_{n}\in G\bigr{)},$$

$$Z_{n,k}:=\frac{\Delta_{n}}{2}\sum_{i=1}^{k}\bigl{[}V(\widetilde{X}_{n}(t_{n,i-1}))+V(\widetilde{X}_{n}(t_{n,i}))\bigr{]}.$$

$$\tau_{n}^{r}:=\inf\bigl{\{}t\geq 0:\|\widetilde{\mathbf{X}}_{n}(t)\|\vee\|\widetilde{\mathbf{X}}_{n}(t-)\|>r\bigr{\}},$$

$${\hat{Z}}_{n,k+1}={\hat{Z}}_{n,k}+V({\tilde{X}}_n(t_{n,k})){\Delta }_n,{\hat{Z}}_{n,0}=0.$$

$$\lim_{n\to\infty}\max_{1\leq k\leq\lfloor n\tau_{n}^{r}\rfloor}\Bigl{\lvert}\Delta_{n}^{-1}\mathbf{E}\bigl{(}\Delta\widehat{Z}_{n,k+1}\,|\,\widehat{Z}_{n,k}=z,\widetilde{X}_{n}(t_{n,k})=x\bigr{)}-V(x)\Bigr{\rvert}=0$$

Full text

1 Introduction

Barrier options are widely used path-dependent derivative securities, and there is extensive literature on their valuation (see e.g. Armstrong2001 , Carr1995 and the references therein). Most of the literature is focused on the case where the barriers are “flat”. However, curved boundaries also appear in the context of pricing options on underlying assets that have a volatility term structure. Using a deterministic time-change, the time dependence in the volatility can be transferred into the boundary, reducing the problem of pricing call and put options in the Black–Scholes type setting with deterministic interest rates to computing the curvilinear boundary crossing probability of a diffusion process (see e.g. Section 3 in Borovkov2005 ). Methods for solving the latter problem are numerous, and we refer the interested reader to Liang2021 for a literature review. To extend the aforementioned problem to the case of stochastic interest rates, one needs to introduce a numeraire process.

It is well-known that, under the no-arbitrage assumption, the fair price of a replicable derivative on an underlying asset $X$ with maturity $T$ and general payoff $\phi(X(T))$ is equal to $\mathbf{E}\phi(X(T))/N(T)$, where $\mathbf{E}$ denotes the expectation with respect to the martingale measure associated with the numeraire process $N$ (for details, see e.g. Chapter 6 in Bingham2004 ). If we introduce a barrier feature for the above option using two time-dependent barriers $g^{-}<g^{+}$, the payoff of the new option will be given by

[TABLE]

Hence the fair price of the barrier option is equal to

[TABLE]

The purpose of this paper is to extend the Markov chain approximation method proposed in Liang2021 to compute the above expectation in the case when the numeraire process takes the form $N(T)=\exp\{\int_{0}^{T}V(X(t))\,dt\}$ for some continuous function $V$, i.e. to evaluate

[TABLE]

The method suggested in Liang2021 was a development of the approach from FuWu2010 for computing boundary non-crossing probabilities of the Brownian motion process. That approach suggested to replace the continuous dynamics of the underlying Brownian motion process with those of a denumerable discrete time Markov chain, by discretising both time and space using uniform grids. The transition probabilities for these Markov chains were specified to be proportional to the values of the Brownian motion transition densities for the respective time and space increments. The desired approximation for the boundary non-crossing probability was then computed in FuWu2010 by multiplying finite-dimensional substochastic matrices obtained by retaining the entries corresponding to the space nodes located between the given boundaries at the respective times.

Our modification of this method developed in Liang2021 dramatically improved the convergence rate of the procedure by shifting and adjusting the uniform space grids so that they would have nodes lying exactly on the boundaries, and by implementing one-step Brownian bridge corrections taking into account the possibility that the process could cross a boundary at a time between the points on the discrete time grid. At the same time, the method was extended in Liang2021 to a much broader class of diffusions.

Roughly speaking, one can think about the workings of this method as summing up the probabilities of all the Markov chain trajectories lying between the boundaries (and applying appropriate corrections). As, for each trajectory, its probability is given by the product of one-step transition probabilities, one basically “tracks” each trajectory of the chain in this computational process. Unlike alternative approaches to computing boundary crossing probabilities, such as the method of integral equations (see e.g. GuRiRoTo1997 and further references therein), this provides one with an opportunity to compute not only probabilities but also expectations of path-dependent functionals, finding the values of expressions of the form (2) being of particular convenience.

Turning back to applications related to (1), we note that when the underlying $X$ is the spot interest rate, this includes the case when the numeraire is the bank account process, i.e. $V(x)=x$. This expression can also be used to price what we call a hybrid step-barrier option, an option that has a step option feature (see Linetsky1999 ) and a barrier above the step level (see Example 3 below).

It was reported in Liang2021 that numerical experiments had showed that the proposed method produces approximations for boundary crossing probabilities with convergence rate $O(n^{-2})$, where $n$ is the number of nodes on the regular time discretisation grid. It turns out that this convergence rate is preserved in the new extension of the method presented in this paper when $V$ is sufficiently smooth. Note that proving the above-mentioned very fast convergence rate for our scheme is a very difficult task already in the basic case of computing boundary crossing probabilities. Establishing this rate in the one-sided boundary case with $g^{+}\in C^{2}$ for the Brownian motion process is work in progress.

The paper is organised as follows. In Section 2 we describe our approximation procedure and demonstrate its convergence. In Section 3 we present the results of applying it in three special cases.

2 The method

We will begin by specifying the diffusion process model for the underlying asset $X$. We assume a unit diffusion coefficient since we can transform a large class of diffusion processes to this case using the well-known unit-diffusion transform (see e.g. Section 2 in Liang2021 ). Then, we will describe the sequence of Markov chain approximations and prove that the sequence converges weakly to the desired diffusion process. Using the weak convergence result, we will show that the corresponding approximations to the option price that are using the Brownian bridge correction techniques converge as well. Recall that employing this correction dramatically improves convergence rates in the basic problem of computing boundary crossing probabilities Liang2021 . As we will see this applies in the present paper setup as well.

Without loss of generality, we set the terminal time $T:=1$ since a deterministic time change can be applied to get the results for a general fixed $T<\infty$.

2.1 The setup

Suppose that the underlying asset’s price is modelled by the one-dimensional diffusion process

[TABLE]

where $\{W(t)\}_{t\geq 0}$ is a standard Brownian motion process, and $x_{0}$ is non-random. Assume the following condition is satisfied:

1. (C)

For any fixed $x\in\mathbb{R}$, one has $\mu(\,\cdot\,,x)\in C^{1}([0,1])$, and for any fixed $t\in[0,1]$, one has $\mu(t,\,\cdot\,)\in C^{2}(\mathbb{R})$. Moreover, for any $r>0$, there exists a $K_{r}<\infty$ such that one has

[TABLE]

We will also need some notations and conditions related to the boundaries $g^{\pm}$. Denote by $C=C([0,1])$ the space of continuous functions $f:[0,1]\to\mathbb{R}$ equipped with the uniform norm $\|f\|_{\infty}:=\sup_{t\in[0,1]}\lvert f(t)\rvert$. For a fixed $x_{0}\in\mathbb{R}$, consider the class

[TABLE]

of pairs of functions from $C$ and introduce the notation

[TABLE]

for the “bunch” of continuous functions whose graphs are entirely contained in the strip between the boundaries $f^{\pm}.$

Let $V:\mathbb{R}\to\mathbb{C}$ be a (possibly complex-valued) continuous function and $\phi:\mathbb{R}\to\mathbb{R}$ be some continuous function. The problem we deal with in this paper is how to compute

[TABLE]

for some $(g^{-},g^{+})\in\mathcal{G}$. To explain our approach, set $\mathbf{E}^{s,x}(\cdot):=\mathbf{E}(\cdot\,|\,X(s)=x)$ and introduce Feynman–Kac type time-dependent semigroup $\{T_{s,t}\}_{0\leq s<t\leq 1}$ defined by

[TABLE]

and represent (3) as

[TABLE]

for

[TABLE]

the uniform partition of $[0,1]$ of rank $\Delta_{n}:=1/n$, $n\geq 1$. The idea of our method is to approximate operators $T_{s,t}$ with their discrete versions. We will implement it in the next section, leading to approximation (8) to our $Q$.

2.2 Markov chain approximation

To specify our time-dependent Markov chain approximation, we first define the space grids $E_{n,k}$ and the transition probabilities $p_{n,k}$, and then introduce the corrective terms $\pi_{n,k}$ and $e_{n,k}$ which account for boundary correction and the presence of the term with $V$ respectively.

The grids $E_{n,k}$ are constructed as follows. Set $g_{n,k}^{\pm}:=g^{\pm}(t_{n,k})$, $k=1,\ldots,n,$ and, for fixed \delta\in(0,\mbox{\frac{1}{2}}] and $\gamma>0$, put

[TABLE]

assuming that $n$ is large enough such that the integer parts in all the denominators are positive. We set the time-dependent space lattice step sizes to be

[TABLE]

Next, we define

[TABLE]

We also put $E_{n,0}:=\{x_{0}\}$ and define the corresponding boundary restricted lattices

[TABLE]

Further, for $k=1,\ldots,n,$ we introduce the discrete time drift and diffusion coefficients

[TABLE]

and define one-step “pseudo-transition probabilities” (we use this expression because they do not sum to one, although the normalising constants converge to 1 so quickly as $n\to\infty$ that it makes no sense to normalise them, see Remark 7 in Liang2021 ) as

[TABLE]

where $\varphi(x\,|\,\mu,\sigma^{2}):=(2\pi\sigma^{2})^{-1/2}e^{-(x-\mu)^{2}/(2\sigma^{2})}$, $x,\mu\in\mathbb{R}$, $\sigma^{2}>0$. We will implement the Brownian bridge correction for two-sided boundary crossing probabilities by using the factor

[TABLE]

Note that we ignore here the highly unlikely event that the trajectory of the Brownian bridge process hits both the upper and lower boundaries in the small time interval $[t_{n,k-1},t_{n,k}]$, see Liang2021 . The following expression will be used as a trapezoidal approximation of the exponential term $\exp\{-\int_{(k-1)/n}^{k/n}V(X(s))\,ds\}$:

[TABLE]

Now we introduce the pseudo-transition matrices $S_{n,k}\in\mathbb{R}^{\lvert E_{n,k-1}^{G}\rvert\times\lvert E_{n,k}^{G}\rvert}$ involving both the boundary crossing correction terms $\pi_{n,k}$ and our adjustments $e_{n,k}$ for the presence of the exponential terms with $V$:

[TABLE]

Denoting by $\phi_{n}:=[\phi(x)]_{x\in E_{n,n}^{G}}$ the $n$-th payoff column vector, our approximation to (3) is given by the sequence of matrix products

[TABLE]

which are discrete analogues of (5).

2.3 Weak convergence

Theorem 2.1

Let condition (C) be met, $(g^{-},g^{+})\in\mathcal{G}$, and $V:\mathbb{R}\to\mathbb{C}$, $\phi:\mathbb{R}\to\mathbb{R}$ be continuous functions. Then $Q_{n}\to Q$ as $n\to\infty$.

Proof

We will use the method of weak convergence. First, we normalise the pseudo-transition probabilities $p_{n,k}$ and prove that the corresponding sequence of Markov chains converges weakly to the target diffusion process $X$ using convergence results from Liang2021 . Then we extend this convergence result to a sequence of auxiliary bivariate processes, where from the desired convergence $Q_{n}\to Q$ will follow as these quantities can be expressed as expectations of a suitable functional.

Let $\{\xi_{n,k}\}_{k=1}^{n}$ denote a Markov chain with transition probabilities $p_{n,k}(x,y)/C_{n,k}(x)$, $(x,y)\in E_{n,k-1}\times E_{n,k}$, where $C_{n,k}(x):=\sum_{y\in E_{n,k}}p_{n,k}(x,y)$. Furthermore, let $\widetilde{X}_{n}$ denote the Brownian bridge-interpolated version of $\{\xi_{n,k}\}_{k=1}^{n}$, i.e.

[TABLE]

where

[TABLE]

and $B_{n,k}^{\circ}(t)$, $t\in[t_{n,k-1},t_{n,k}]$, are independent Brownian motions “pinned” at the time-space points $(t_{n,k-1},0)$ and $(t_{n,k},0)$, these bridges being independent of our chain. From Corollary 1 in Liang2021 , we know that $\widetilde{X}_{n}\Rightarrow X$ as $n\to\infty$, where “$\Rightarrow$” denotes weak convergence of random elements in the respective functional space (in this case, in space $(C,\|\cdot\|_{\infty})$).

Define the two-dimensional process $\mathbf{X}(t):=(X(t),Z(t))$ by letting

[TABLE]

Then we can rewrite $Q$ from (3) as

[TABLE]

where $\Psi(\mathbf{X}(1)):=\phi(X(1))e^{Z(1)}$. Our approximation (8) can also be rewritten as

[TABLE]

where $\mathbf{X}_{n}(t):=(\widetilde{X}_{n}(t),Z_{n,\lfloor nt\rfloor})$, $Z_{n,0}=0$ and, for $k=0,1,\ldots,n,$

[TABLE]

We will now show that $(\widetilde{X}_{n},Z_{n})\Rightarrow(X,Z)$ in the respective Skorokhod space by verifying the consistency of infinitesimal moments and tightness. We first define the auxiliary process $\widehat{Z}_{n}(t):=\widehat{Z}_{n,\lfloor nt\rfloor}$, $t\in[0,1]$, with $\widehat{Z}_{n,k}=\sum_{i=0}^{k-1}V(\widetilde{X}_{n}(t_{n,i}))\Delta_{n}$ and show that $(\widetilde{X}_{n},\widehat{Z}_{n})\Rightarrow(X,Z)$.

Define the joint process $\widetilde{\mathbf{X}}_{n}(t):=(t,\widetilde{X}_{n}(t),\widehat{Z}_{n,\lfloor nt\rfloor})$ and the stopping time

[TABLE]

$\|\mathbf{u}\|:=\lvert u_{1}\rvert\vee\lvert u_{2}\rvert\vee\lvert u_{3}\rvert$ being the Chebyshev norm of $\mathbf{u}=(u_{1},u_{2},u_{3})\in\mathbb{R}^{3}$. Note that

[TABLE]

Let $\Delta\widehat{Z}_{n,k+1}:=\widehat{Z}_{n,k+1}-\widehat{Z}_{n,k}$. It follows immediately that, for all $r>0$,

[TABLE]

and

[TABLE]

Tightness of the sequence of the distributions of $\widetilde{\mathbf{X}}_{n}$ in the respective Skorokhod space is immediate since $V$ is bounded on compact sets and hence it follows from Corollary 4.2 in EthierKurtz that $(\widetilde{X}_{n},\widehat{Z}_{n})\Rightarrow(X,Z)$. For all $r>0,$

[TABLE]

and hence it follows that $(\widetilde{X}_{n},Z_{n})\Rightarrow(X,Z)$ as well. By the continuous mapping theorem applied to the function $(x,y)\mapsto\phi(x)e^{y}$, the Portmanteau theorem (see e.g. p. 24 in Billingsley1968 ), and Corollary 2 from Liang2021 , the claimed result follows from the established weak convergence of the processes and boundedness of the integrands in $Q$. $\Box$

3 Numerical examples

In this section, we will illustrate the efficiency of our approximation scheme (8). The run times for these computations using the programming language Julia run on a MacBook Pro 2020 laptop computer with an i5 processor (2 GHz, 16 GB RAM) are basically the same as the ones reported in Liang2021 for computing boundary-crossing probabilities employing the same hardware and software. To evaluate the partial derivatives in (6), we used the package HyperDualNumbers.jl

We apply our algorithm in three different examples, where $V(x)$ is equal to $-ix^{2}$, $x$, and $\kappa\mathbf{1}\{x>r\}$ respectively, for some $x,\kappa,r\in\mathbb{R}$. Due to the discontinuity of the function $V$ at $x=r$ in the last case, we modify the scheme slightly for that choice of $V$.

For each example, we compute $Q_{n}$ and plot the observed convergence rate of $Q_{n+1}-Q_{n}$ to zero as $n$ grows. Closed-form expressions for $Q$ are not available in any of these examples. Due to the iterative nature of our scheme (8), as a byproduct, we obtain approximations to $v(t,x):=(T_{t,1}\phi)(x)$, where $T_{s,t}$ appeared in (4). We will plot the surfaces $v(t,x)$ for each example as well.

3.1 The case $V(x)=-ix^{2}$

To demonstrate general methods for obtaining distributional properties of Wiener functionals of the form $\int_{0}^{1}V(W(s))\,ds$, the characteristic function of the random variable $\int_{0}^{1}W(s)^{2}\,ds$ was computed in Kac1949 and CameronMartin1945 yielding

[TABLE]

Using our scheme, we compute approximations to the expressions of this form restricted to the event of non-crossing given boundaries:

[TABLE]

with time-dependent boundaries

[TABLE]

In Fig 1, we plot our approximation of the form (8) to $v(t,x)$. We see from Fig 2, that $\lvert\operatorname{Re}Q_{n+1}-\operatorname{Re}Q_{n}\rvert$ converges to zero at the rate $O(n^{-3})$, suggesting that $\lvert\operatorname{Re}Q_{n}-\operatorname{Re}Q\rvert=O(n^{-2})$ as $n\to\infty$.

3.2 The case $V(x)=x$

In this example, we price a zero-coupon bond with a two-sided time-dependent knock-out barrier, using a one-factor Hull–White model for spot interest rates. This problem was also explored in Kuan2003 . For $\alpha,\sigma>0,$ the Hull–White model assumes that the spot interest rate $r$ has the following dynamics:

[TABLE]

where $\theta(t)$ depends on $\alpha$, $\sigma$ and the currently observed instantaneous forward rate curve $f(t)$ (see e.g. Proposition 10.1.6 in Piterbarg2010):

[TABLE]

assuming that $r(0)=f(0)$. To apply our scheme to this example, we scale all the space variables by $\sigma$ to transform $r(t)$ into a unit-diffusion process. Using the martingale pricing theorem, the fair price of a one-year maturity zero-coupon bond with a barrier option feature is equal to

[TABLE]

We used the time-dependent boundaries given by

[TABLE]

To demonstrate the performance of our scheme, we set $\alpha=0.01$, $\sigma=0.01$ and a flat instantaneous forward curve $f(t)\equiv 0.03$. In the left pane in Fig 3, we plot our approximation of

[TABLE]

for different initial values using our approximation algorithm (8). We see from the right pane in Fig 3 that $Q_{n+1}-Q_{n}$ converges to zero at the rate $O(n^{-3})$, suggesting that $\lvert Q_{n}-Q\rvert=O(n^{-2})$ as $n\to\infty$.

3.3 The case $V(x)=\kappa\mathbf{1}\{x>r\}$

The standard barrier option is one of the most popular first-generation exotic derivatives, due to its cheaper price compared to their vanilla counterparts. If an investor believes that a certain price is unlikely to fall below a certain level, they can choose to include a knock-out feature in the option, which reduces the price. These cost reductions can be substantial when volatility is high. However, standard barrier options come with a few disadvantages. Option buyers can lose their entire investment due to short time price spikes once the price is near the barrier. The delta sensitivity of such a barrier option is discontinuous near the boundary, making it difficult to hedge for options dealers.

The so-called “step option” was introduced in Linetsky1999 to address these issues. The payoff of an “up-and-out” step option with finite knock-out rate $\kappa\geq 0$ written on an option with payoff $\phi$ is given by the formula

[TABLE]

where $\kappa>0$ is a constant and $\tau_{r}(t)$ is the sojourn time of the underlying price process $X$ above the level $r>0$ until time $t$:

[TABLE]

We will assume $X$ to be the standard Brownian motion process $W$ in this example.

This option can be modified by adding barrier features that will further reduce its price. We introduce a time-dependent knock-out barrier $g^{+}$ above the step option barrier $r$ and a lower knock-out barrier $g^{-}$. We call this option a “hybrid step-barrier option”. Its payoff is given by

[TABLE]

The expectation of this functional is a special case of our $Q$ from (3) with $V(x)=\kappa\mathbf{1}\{x>r\}$. In this case, the trapezoidal approximation (7) converges very slowly due to the discontinuity of $V$ at the occupation level $r$. We now describe a modification of $e_{n,k}$ to address this issue. Let $\mathbf{E}^{t,x,y}(\cdot):=\mathbf{E}(\,\cdot\,|\,W(0)=x,W(t)=y)$. Conditioning $W$ on endpoints $(0,x)$ and $(t,y)$, one can expect that for small times $t$, the sojourn time $\tau_{r}(t)$ is “almost independent” of the boundary crossing event:

[TABLE]

There is a semi-explicit formula for $\mathbf{E}^{t,x,y}e^{-\kappa\tau_{r}(t)}$ given in (1.4.7) in Borodin2002 . However, efficient numerical computations of the double convolution present in that formula was challenging in the case when $x<r$ and $y>r$. Hence we decided to use the following simple approximation.

Lemma 1

For $x,y,\kappa\in\mathbb{R}$, as $t\downarrow 0$,

[TABLE]

where $\Phi(\cdot)$ is the standard normal distribution function and $\overline{\Phi}(x):=1-\Phi(x)$.

Proof

Since $\tau_{r}(t)\leq t$, a Taylor series expansion of $e^{-x}$ yields

[TABLE]

Changing the order of integration using Fubini’s theorem we get

[TABLE]

Since

[TABLE]

one has

[TABLE]

$\Box$

Remark 1

One could actually compute the second-order term in the expansion for $\mathbf{E}^{t,x,y}e^{-\kappa\tau_{r}(t)}$ using the following formula for the second moment of $\tau_{r}(t)$:

[TABLE]

Let $W^{t,x,y}:=\{W^{t,x,y}(s):s\in[0,t]\}$ denote the Brownian motion pinned at $(0,x)$ and $(t,y)$. For brevity, we also set $W^{t}:=W^{t,0,0}$. It is well-known that, for $s<t$,

[TABLE]

and hence

[TABLE]

where

[TABLE]

Now since $\mathbf{E}W^{t}(s)W^{t}(u)=(t-u)s/t$, $s<u$, we get

[TABLE]

Therefore

[TABLE]

where $\Phi_{\rho}$ denotes the bivariate standard normal cumulative distribution function with correlation $\rho$ given by (10). Hence we have

[TABLE]

However, numerically computing this integral turned out to be too computationally expensive and did not improve the efficiency of the algorithm.

To summarise, our modification of the term $e_{n,k}$ in this special case when $V$ is discontinuous amounts to replacing (7) with

[TABLE]

In our numerical implementation, we computed the time integral appearing in (11) using Gaussian quadratures. To demonstrate the performance of our scheme, we set $r=1/19$, $\kappa~{}=~{}2$ and $g^{\pm}(t)=\pm 4\mp t^{2}$. Without replacing $e_{n,k}(x,y)$ with (11), the convergence of $\lvert Q_{n+1}-Q_{n}\rvert$ to zero is highly non-smooth and slow. With the proposed modification, the scheme appears to converge at rate $O(n^{-1})$, as shown in the right pane of Fig 4. The left pane in Fig 4 shows our approximation for the values of

[TABLE]