Alternative models for FX, arbitrage opportunities and efficient pricing of double barrier options in Lévy models

Alternative models for FX, arbitrage opportunities and efficient pricing of double barrier options in L\'evy models

{\sc Key words:} Heston model, L\'evy processes KoBoL, double barrier options, Wiener-Hopf factorization, Fourier transform, Laplace transform, Gaver-Wynn Rho algorithm, sinh-acceleration, double spiral method

{\sc MSC2020 codes:} 60-08,42A38,42B10,44A10,65R10,65G51,91G20,91G60

Introduction

There exists a large body of literature devoted to calculation of the expectation of a function of a L\'evy process and its running extremum, and related optimal stopping problems, standard examples being barrier and American options, and lookback options with barrier and/or American features. However, sufficiently accurate and fast algorithms are difficult to develop unless a process is of a rather simple nature. The errors can be especially large for prices of double barrier options in L\'evy and regime-switching L\'evy models. The importance of a reliable and fast pricing procedure for barrier options in regime-switching L\'evy models is apparent in the following situation. Typically, the best fit to empirically observed vanilla prices (or implied volatilities) is achieved with jump-diffusion models, and, among models with small number of parameters, with purely jump models. In fact, in many cases, finite variation L\'evy processes are calibrated better than the processes of infinite variation (see, e.g., examples in CGMY,ChenFengLin,levendorskii-xie-Asian,AsianGamma). However, as an extensive empirical study in dahlbokum07 demonstrates, the prices of barrier options provided by Markit agree better with prices in diffusion models. If the latter prices are the real prices, the market prices vanillas and barrier options using non-equivalent probability measures, hence, there may exist arbitrage opportunities. Using the well-known criterion of the equivalency of L\'evy measures (see., e.g., sato) one can easily prove that in several empirical examples documented in CGMY, the historic and risk-neutral measures inferred from the prices of a stock and vanillas are not equivalent. The differences among prices of vanillas in many models being small, the apparent violations of the no-arbitrage condition can be explained by the transaction cost.

In the case of barrier options, the differences in prices can be quite substantial due to boundary effects. The asymptotic analysis in NG-MBS,early-exercise,BIL,asymp-sens for single barrier options shows that, in many cases, close to the barrier, the prices in purely jump models are much larger than the prices in models with a sizable diffusion component. The numerical examples in BIL demonstrate that the difference is sizable at the distance of several percent from the barrier. Using the fairly complicated but exact analytical formula for prices of the double barrier options in terms of the inverse Laplace-Fourier transform (this formula is implicitly used in the method of the paper), one can derive similar but less explicit asymptotic formulas for double barrier options and conclude that, in the case of the double barrier options, the differences between prices in purely jump models and models with sizable diffusion components must be rather large, especially if the distance between the barriers is fairly small, e.g., 10%. Even if a high transaction cost is taken into account, the difference between the prices in the models used by the counterparties to price double barrier options may create arbitrage opportunities. The simplest way to realize the arbitrage opportunity in this setting is to buy a double-no-touch options (DNT) from a bank and keep it until maturity. If the discounted probability of the underlying remaining in the corridor between the two barriers until maturity, calculated in a jump model, is larger than the ask price, possibly, calculated at the bank using the Heston model, then one has an apparent arbitrage opportunity. The first aim of the paper is to analyze this opportunity using a small set of the real data; the second aim of the paper is to analyze analytical difficulties for accurate pricing double barrier options in jump models, and develop an efficient method for evaluation of double barrier options in one factor L\'evy models. The algorithm developed in the present paper is much more accurate and dozens of times faster than the one in BLdouble; the main block of the algorithm is the realization in the dual space of the main block in BLdouble, where the calculations are in the state space. We apply the Gaver-Wynn-Rho algorithm (GWR algorithm) and, for each value of the (positive) spectral parameter, apply the same algorithm as in BLdouble for perpetual double barrier options but do the calculations in the dual space. The conformal deformation technique (sinh-acceleration) developed in SINHregular,Contrarian,EfficientLevyExtremum,EfficientDoubleBarrier allows one to evaluate the prices of the perpetual options with the accuracy E-09 and better using arrays of the length of a hundred and shorter. The error tolerance of the order of E-14 requires arrays several hundreds long. We call the method in the paper the GWR-SINH method. Since the analytic properties of the characteristic exponent of L\'evy processes considered in the present paper are significantly worse that the properties of L\'evy processes considered in EfficientDoubleBarrier, the Laplace transform of the price of barrier options does not admit analytic continuation to the complex plane with the cut of the form $(-\infty, a]$ as in EfficientDoubleBarrier. In the result, the error of the GWR method in this paper is 1-2 orders of magnitude larger than that in EfficientDoubleBarrier, which agrees with the general characterization of the accuracy of the Gaver method in AbateValko04,AbValko04b.

The rest of the paper is organized as follows. In Sect. (ref), we collect the definitions of classes of L\'evy processes used in the paper and necessary basic facts of the Wiener-Hopf factorization. In Sect. (ref), we recall the scheme of evaluation of the Laplace transform BLdouble using the calculations in the state space. GWR-SINH method and the explicit algorithm are in Sect. (ref). In Sect. (ref), we discuss several popular approaches to pricing barrier options and analyze common sources of errors of various groups of methods. In Sect. (ref), we summarize the results of the paper and outline the extension of the GWR-SINH method to regime-switching models and models with the stochastic volatility and interest rates. Calibration results and prices of DNT options in the Heston and KoBoL models are in Section (ref). The GWR-SINH method for regime-switching L\'evy models is in Section (ref).

Auxilliary results

Let $X$ be a one-dimensional L\'evy process on the filtered probability space $(\Omega, {\mathcal F}, \{{\mathcal F}_t\}_{t\ge 0}, {\mathbb Q})$ satisfying the usual conditions; the riskless rate is constant and ${\mathbb Q}$ is an equivalent martingale measure. We denote the expectation operator under ${\mathbb Q}$ by ${\mathbb E}$. The underlying is modeled as $S_t=e^{X_t}$, the barriers are $H_-<H_+$, and the maturity date is denoted $T$. We set $h_\pm=\ln H_\pm$ and $x=\ln S$. The infimum process and supremum process ${\bar X}_t=\sup_{0\le s\le t}X_s$ and ${\underline X}_t=\inf_{0\le s\le t}X_s$ are defined pathwise, a.s. For $h\in {\mathbb R}$, $\tau^+_h$ and $\tau^-_h$ denote the first entrance time of $X$ into $[h,+\infty)$ and $(-\infty, h]$, respectively. For $q>0$, $T_q\sim\operatorname{Exp}q$ denotes an exponentially distributed random variable with mean $q^{-1}$, independent of $X$. $Q_{GWR}$ denotes the set of nodes in the GWR algorithm used in Sections (ref) and (ref).

Wiener-Hopf factorization

This subsection contains basic formulas and results systematically used in a number of publications, e.g., NG-MBS,barrier-RLPE,paired,EfficientDiscExtremum,EfficientDoubleBarrier. In probability, the Wiener-Hopf factors are defined as

Functions $\phi^\pm_q(\xi)$ appear in the Wiener-Hopf factorization formula

Define the expected present value operators (EPV-operators) under $X$, ${\bar X}$ and ${\underline X}$ (all three start at 0) by ${\mathcal E}_q u(x)={\mathbb E}\left[u(x+X_{T_q})\right]$, ${\mathcal E}^+_q u(x)={\mathbb E}\left[u(x+{\bar X}_{T_q})\right]$ and ${\mathcal E}^-_q u(x)={\mathbb E}\left[u(x+{\underline X}_{T_q})\right]$. The EPV operators are bounded operators in $L_\infty({\mathbb R})$. In the case of L\'evy processes with exponentially decaying L\'evy densities, the EPV operators are bounded operators in spaces with exponential weights. The operator version of (ref)

is a special case of the operator form of the Wiener-Hopf factorization in the theory of boundary problems for differential and pseudo-differential operators (pdo). Indeed, ${\mathcal E}^\pm_q e^{ix\xi}=\phi^\pm_q(\xi) e^{ix\xi}$. This means that ${\mathcal E}^\pm_q$ are pdo with symbols $\phi^\pm_q$, and ${\mathcal E}^\pm_qu(x)={\mathcal F}^{-1}_{\xi\to x}\phi^\pm_q(\xi){\mathcal F}_{x\to\xi}u(x)$ for sufficiently regular functions $u$. See e.g., eskin,NG-MBS. The Wiener-Hopf factor $\phi^+_q(\xi)$ (resp., $\phi^-_q(\xi)$) admits analytic continuation to the half-plane $\{\operatorname{\rm Im}\xi>0\}$ (resp., $\{\operatorname{\rm Im}\xi<0\}$).

The characteristic exponents of all popular classes of L\'evy processes bar stable L\'evy processes admit analytic continuation to a strip around the real axis. See NG-MBS,barrier-RLPE,BLSIAM02, where the general class of Regular L\'evy processes of exponential type (RLPE) is introduced. Let $X$ be a L\'evy process with the characteristic exponent admitting analytic continuation to a strip $\{\operatorname{\rm Im}\xi\in (\mu_-,\mu_+)\}$ around the real axis, and let $q>0$. Then (see., e.g., NG-MBS,barrier-RLPE,paired) the following statements hold. \medbreak I. There exist $\sigma_-(q)<0<\sigma_+(q)$ such that

\medbreak II. The Wiener-Hopf factor $\phi^+_q(\xi)$ admits analytic continuation to the half-plane $\{\operatorname{\rm Im}\xi>\sigma_-(q)\}$, and can be calculated as follows: for any $\omega_-\in (\sigma_-(q), \operatorname{\rm Im}\xi)$,

\medbreak III. The Wiener-Hopf factor $\phi^-_q(\xi)$ admits analytic continuation to the half-plane $\{\operatorname{\rm Im}\xi<\sigma_+(q)\}$, and can be calculated as follows: for any $\omega_+\in (\operatorname{\rm Im}\xi, \sigma_+(q))$,

We can (and will) use (ref) and ((ref)) to calculate ${\phi^-_q}(\xi)$ on $\{\operatorname{\rm Im}\xi\in (0,\sigma^+_q)\}$, and (ref) and ((ref)) to calculate ${\phi^+_q}(\xi)$ on $\{\operatorname{\rm Im}\xi\in (\sigma^-_q,0)\}$.

The explicit algorithm formulated and used in EfficientDoubleBarrier for processes of infinite variation and finite variation processes without drift uses the representations ((ref))-((ref)). In the case of finite variation processes with positive drift, which we consider in the paper, the following set of formulas derived in NG-MBS,paired is more efficient:

The formulas for the case of finite variation processes with negative drift are by symmetry.

General classes of L\'evy processes amenable to efficient calculations

Essentially all popular classes of L\'evy processes processes enjoy additional properties formalized in SINHregular,EfficientAmenable as follows. For $\nu=0+$ (resp., $\nu=1+$), set $|\xi|^\nu=\ln|\xi|$ (resp., $|\xi|^\nu=|\xi|\ln|\xi|$), and introduce the following complete ordering in the set $\{0+,1+\}\cup (0,2]$: the usual ordering in $(0,2]$; $\forall\ \nu>0, 0+<\nu$; $\forall\ \nu>1, 1<1+<\nu$. For $\gamma\in (0,\pi]$, $\gamma_+\in (0,\pi/2]$, $\gamma_-\in [-\pi/2,0)$ and $\mu_-<\mu_+$, define ${\mathcal C}_{\gamma_-,\gamma_+}=\{e^{i\varphi}\rho\ |\ \rho> 0, \varphi\in (\gamma_-,\gamma_+)\cup (\pi-\gamma_+,\pi-\gamma_-)\}$, ${\mathcal C}_{\gamma}=\{e^{i\varphi}\rho\ |\ \rho> 0, \varphi\in (-\gamma,\gamma)\}$, $S_{(\mu_-,\mu_+)}=\{\xi\ |\ \operatorname{\rm Im}\xi\in (\mu_-,\mu_+)\}$.

To simplify the constructions in the paper, we assume that $\mu_-<0<\mu_+$ and $\gamma'_-<0<\gamma'_+$.

In EfficientAmenable, we defined a class of Stieltjes-L\'evy processes (SL-processes). Essentially, $X$ is called a (signed) SL-process if $\psi$ is of the form

where $ST({\mathcal G})$ is the Stieltjes transform of the (signed) Stieltjes measure ${\mathcal G}$, $a^\pm_j\ge 0$, and $\sigma^2\ge0$, $\mu\in{\mathbb R}$. A (signed) SL-process is called regular if it is SINH-regular. We proved in EfficientAmenable that the characteristic exponent $\psi$ of a (signed) SL-process admits analytic continuation to the complex plane with two cuts along the imaginary axis. If $X$ is an SL-process, then, for any $q>0$, equation $q+\psi(\xi)=0$ has no solution on ${\mathbb C}\setminus i{\mathbb R}$. We proved that all popular classes of L\'evy processes bar the Merton model and Meixner processes are regular SL-processes, with $\gamma_\pm=\pm \pi/2$; the Merton model and Meixner processes are regular signed SL-processes, and $\gamma_\pm=\pm \pi/4$. In EfficientAmenable, the reader can find a list of SINH-processes and SL-processes, with calculations of the order and type.

Double barrier options: calculations in the state space

Laplace transform and its inversion

Let $r\in {\mathbb R}$, $T>0$, $h_-<x<h_+$, and $G\in L_\infty((h_-,h_+))$. To evaluate

we take $r_0\in {\mathbb R}$ and represent $V(G; r; h_-,h_+;T,x)$ in the form $V(G;h_-,h_+;T,x)=e^{rT} V(G; r+r_0; h_-,h_+;T,x)$. First, we apply the Laplace transform w.r.t. $T$ to $V(G; r+r_0; h_-,h_+;T,x)$. Then, assuming that $r+r_0>0$ is sufficiently large so that $ {\tilde V}(G; r+r_0; h_-,h_+;q',x)$ can be efficiently calculated for $q'$ in the right half-palne, we evaluate the Bromwich integral

where $\sigma>0$ is arbitrary, using an appropriate numerical method. Finally, we multiply the result by $e^{r_0T}$. In the case $\nu\in (0,1)$ and $\mu\neq 0$, the most efficient method of the Laplace inversion of complicated functions, namely, the sinh-acceleration used in EfficientDoubleBarrier, is not applicable. Therefore, we use the GWR algorithm. Let $Q_{GS}$ be the set of spectral parameters used in the GWR algorithm. We need to calculate the prices of perpetual options ${\tilde V}(G; h_-,h_+;q,x)$ for $q\in r+r_0+Q_{GS}$. The first advantage of using $r_0$ is that we can ensure that $q+\psi(\xi)\not\in (-\infty,0]$ for all $q,\xi$ of interest, which simplifies the verification of several important technical conditions in the paper, starting with conditions for an efficient evaluation of the Wiener-Hopf factors. The second advantage is that the condition $q+\psi(\xi)\not\in (-\infty,0]$ allows for deformations of the contours of integration in the formulas for the Wiener-Hopf factors and ${\tilde V}(G; h_-,h_+;q'+r+r_0,x)$ that make it possible to evaluate the prices of the perpetual options with the precision E-10 and better. The third advantage stemming from the second one is that if the differences between prices of the option with the finite time horizon calculated for different $r_0$'s are of the order of E-05, then E-05 is the order of the error of the GWR method itself. The disadvantage is that if $T$ is large, and the calculations are with double precision, then the factors $e^{-r_0T}$ and $e^{r_0T}$ may introduce large errors. This difficulty can be resolved using an additional trick in paired. However, in the case of double barrier options with small distances between barriers, the price of options of long maturity is negligible, hence, we may assume that $T$ is not large.

Evaluation of the perpetual double barrier options

The notation and scheme in this Section are borrowed from BLdouble,EfficientDoubleBarrier; we need the scheme as the starting point for the scheme where the calculations are in the dual space. Let $q\in r+r_0+Q_{GS}$ and $x\in (h_-,h_+)$. We have

Let $lG$ be a bounded measurable extension of $G$ to ${\mathbb R}$. Set ${\tilde V}^0(lG; q,\cdot)=q^{-1}{\mathcal E_q} lG$, and calculate ${\tilde V}^1(lG;h_-,h_+;q,x):={\tilde V}(G;h_-,h_+;q,x)-{\tilde V}^0(lG; q,x)$ as the sum of a series exponentially converging in $L_\infty$-norm. The terms of the series depend on the choice of the extension but ${\tilde V}(G;h_-,h_+;q,x)$ is independent of the choice. Define

and note that ${\tilde V}^+_1$ (resp., ${\tilde V}^-_1$) is the EPV of the stream $G(X_t)$ which starts to accrue the first time $X_t$ crosses $h_+$ from below (resp., $h_-$ from above). Inductively, for $j=2,3,\ldots,$ define

In BLdouble, the following key theorem is derived from the stochastic continuity of $X$.

Under additional weak conditions on $X$, Theorems 11.4.4 and 11.4.5 in IDUU state that

in single, ((ref))-((ref)) are proved for any L\'evy process. Similar representations for ${\tilde V}^\mp_j$, $j=2,3,\ldots$:

follow from Theorems 11.4.6 and 11.4.7 in IDUU. See EfficientDoubleBarrier for details.

GWR-SINH method

Calculation of the Wiener-Hopf factors for $q>0$. The case $\mu>0$ and $\nu\in (0,1)$

Set $q_0=\min Q_{GWR}$ . First, we construct the curves ${\mathcal L}^+_{\omega_{1,+}, b_+,\omega_+}$ and ${\mathcal L}^-_{\omega_{1,-}, b_-,\omega_-}$ in the upper and lower half-planes, respectively, such that

We deform the contours of integration in ((ref)) and ((ref)) and calculate

To calculate ${\phi^+_q}(\xi)$, $\xi\in{\mathcal L}^-_{\omega_{1,-}, b_-,\omega_-}$, and ${\phi^-_q}(\xi)$, $\xi\in{\mathcal L}^+_{\omega_{1,+}, b_+,\omega_+}$, we use ((ref)).

General scheme for a fixed $q>0$

The numerical scheme in this section has small but crucial differences from the scheme in EfficientDoubleBarrier. Writing the program, the reader can check that the scheme formulated in EfficientDoubleBarrier blows up for processes of finite variation with positive drift. We indicate the changes when they appear. For the same of brevity, we formulate the scheme for the case of DNT options. (The reader can easily adjust the schemes for call and put options and digitals using the first step of the scheme in EfficientDoubleBarrier.) With the exception of the last step, when the inverse Fourier transform is applied, the calculations are in the dual space. The Fourier transforms ${\hat{\tilde V}}^\pm_j(h_-,h_+;q,\xi)$ of functions ${\tilde V}^\pm_j(h_-,h_+;q,x)$ are evaluated on the curves ${\mathcal L}^\mp:={\mathcal L}^\mp_{\omega_\mp, b_\mp,\omega_{1,\mp}}$ used to evaluate the Wiener-Hopf factors; changing the curves, we can double-check the accuracy of calculations. In the case of DNT, ${\tilde V}^0(q,x)=1/q$, hence,

For $j=1,2,\ldots, $ define

Evidently,

and ${\hat W}^\pm_1(h_-,h_+,q,\xi)=\pm 1/(i\xi)=\mp i/\xi$, $\xi\in {\mathcal L}^\mp$. For $j=2,3,\ldots,$ ${\hat W}^\pm_j$ are calculated inductively. We rewrite ((ref))-((ref)) in terms of $ {\hat W}^\pm_j$ as follows. For $\xi\in {\mathcal L}^-$, we have

(we can apply Fubini's theorem because there exists $c>0$ such that $\operatorname{\rm Im}(\eta-\xi)>c|\eta|$ for $\eta\in {\mathcal L}^+$ and $\xi\in {\mathcal L}^-$). Simplifying,

Similarly, for $\xi\in {\mathcal L}^+$, we calculate

Simplifying,

In the cycle $j=2,3\ldots,$ we calculate ${\hat W}^\pm_{j}(h_-,h_+; q,\xi)$ and partials sums of the series

Finally, for $x\in (h_-, h_+)$, we calculate

then

and

Efficient numerical evaluation of the series ((ref))

Define operators\\ ${\mathcal K}_{-+}(={\mathcal K}_{-+}(q; {\mathcal L}^+; h_-,h_+))$ and ${\mathcal K}_{+-}(={\mathcal K}_{+-}(q; {\mathcal L}^-; h_-,h_+))$ by

We write ((ref)) and ((ref)) as

and calculate ${\hat W}^\pm_{j+1}$ and the partial sums in the cycle in $j=1,2,\ldots, M_0$. For the choice of the truncation parameter $M_0$ given the error tolerance, see Remark (ref). This choice is made assuming that the total error of calculation of the individual terms is sufficiently small and can be disregarded.

We calculate ${\hat W}^\pm_{j+1}$ at points of sinh-deformed uniform grids $\xi^\mp_k=i\omega^{1,\mp}+b^\mp \sinh(i\omega^\mp+y^\mp_k)$, $y^\mp_k=\zeta^\mp k$, $k\in {\mathbb Z}$, on ${\mathcal L}^\mp$, truncate the grids, and approximate the operators ${\mathcal K}_{-+}, {\mathcal K}_{+-}$ with the corresponding matrix operators. The parameters of the conformal deformations and corresponding changes of variables and steps $\zeta^\pm$ are chosen as in Contrarian,EfficientLevyExtremum,EfficientDoubleBarrier. The truncation parameters $N^\pm$ are chosen taking into account the exponential rate of decay of the kernels of the integral operators ${\mathcal K}_{-+}$ and ${\mathcal K}_{+-}$ w.r.t. the second argument, and exponential decay of $e^{i(x-h_-)\xi}$ as $\xi\to\infty$ along ${\mathcal L}^+$, and $e^{i(x-h_+)\xi}$ as $\xi\to\infty$ along ${\mathcal L}^-$. In the $y^\pm$-coordinates, the rate of decay is double-exponential, hence, the truncation parameters $\Lambda^\pm=N^\pm\zeta^\pm$ sufficient to satisfy a small error tolerance $\epsilon$ are moderately large. As a simple rule of thumb, we suggest to choose $\Lambda^\pm$ so that $ \exp[b^-(x-h_+)\kappa_-\sin|\omega^-|e^{\Lambda^-}]<\epsilon,\ \exp[b^+(h_--x)\kappa_+\sin(\omega^+)e^{\Lambda_+}]<\epsilon,$ where $\kappa_\pm\in (0,0.5)$, e.g., $\kappa_\pm=0.4$. Note that a choice of a smaller $\kappa_\pm$, e.g., $\kappa_\pm=0.3$, does not increase $N_\pm$ significantly but makes the prescription more reliable. Keeping the notation ${\mathcal K}_{-+}$ and ${\mathcal K}_{+-}$ for the matrices, we calculate the matrix elements as follows:

Algorithm of GWR-SINH method for DNT options

Steps I-V are preliminary ones; Steps VI-IX are performed for each $q\in r+r_0+Q_{GWR}$ used in the GWR method; the calculations can be easily parallelized.

A numerical example: Table (ref)

The calculations were performed in MATLAB 2017b-academic use, on a MacPro Chip Apple M1 Max Pro chip with 10-core CPU, 24-core GPU, 16-core Neural Engine 32GB unified memory, 1TB SSD storage. The CPU times shown can be significantly improved using parallelized calculations of the Wiener-Hopf factors and the main block, for each $q$ used in the Laplace inversion procedure. The parallelization w.r.t. $H_\pm, S$ is also possible.

Pricing barrier options in L\'evy models: general discussion and examples

A short review of the literature

The general formulas for single barrier options with continuous monitoring were derived in KoBoL,barrier-RLPE,NG-MBS using the operator form of the Wiener-Hopf factorziation eskin, under certain regularity conditions on the characteristic exponent. In single, the same formulas were proved for any L\'evy process. The pricing formulas are in terms of Laplace-Fourier inversion in dimensions 2 (first touch digitals and no-touch options), 3 (barrier puts and calls, and joint probability distributions of a L\'evy process and its extremum), and 4 (more general options with lookback and barrier features). Even marginally accurate realizations of these formulas are far from trivial unless the characteristic exponent $\psi$ of the process is a rational function, hence, the Wiener-Hopf factors are rational functions as well. The factors are especially simple in the Double exponential jump-diffusion model (DEJD model) used in kou,KW1 and its generalization: Hyper-exponential jump-diffusion model (HEJD model) constructed independently in lipton-columbia,lipton-risk (see also LiptonSelection) and amer-put-levy-maphysto,amer-put-levy. In lipton-columbia,lipton-risk, an explicit pricing formula for the joint distribution of the L\'evy process and its extremum was derived using the Gaver-Stehfest algorithm (GS algorithm); the formula can be used to price options with barrier-lookback features. Later, a variation of the same technique was used in structural default models lipton-sepp. In amer-put-levy-maphysto,amer-put-levy, American options with finite time horizon are priced using the maturity randomization technique (Carr's randomization). The Wiener-Hopf factors in HEJD model are represented as weighted sums of exponential functions, therefore, the calculations in the state space are reducible to application of convolution operators with exponentially decaying kernels. A straightforward numerical realization of the convolution operators with exponential kernels is very fast: see amer-put-levy-maphysto,amer-put-levy,ExitRSw,stoch-int-rate-CF for an explicit algorithm. An evident simplification of Carr's randomization can be applied to barrier options (see MSdouble, where double-barrier options in regime-switching models are priced): the early exercise boundary is fixed and it is unnecessary to fund an approximation to the boundary at each step of backward induction. In both cases (GS-algorithm and Carr's randomization), the main block is the evaluation of the perpetual options. If the GS-algorithm is used, it may be necessary to use high precision arithmetic because the weights are very large (see, e.g., examples in paraLaplace). The performance of the GS-algorithm can be improved using Wynn-Rho acceleration (GWR algorithm) - see AbateValko04 and, for examples in the context of pricing options of long maturity, paired. If the GWR-algorithm can be used with double precision arithmetic, then, typically, the CPU time is smaller than if Carr's randomization is applied. In the case of more general L\'evy processes, efficient calculations are much more difficult because the option price is very irregular at the barrier and maturity. See Section (ref) for details and examples. The irregular behavior makes it difficult to evaluate the prices of perpetual options sufficiently accurately so that the GWR-algorithm or Carr's randomization used in paraLaplace,paired,UltraFast can produce reasonably accurate results in regime-switchning models. The latter remark concerns other methods available in the literature - see, e.g., GSh,AAP,AKP,CV,beta,KudrLev09,KudrLev11,BIL,HaislipKaishev14,FusaiGermanoMarazzina,kirkbyJCompFinance18,Linetsky08,feng-linetsky09,LiLinetsky2015 and the bibliographies therein. The fundamental reasons for serious difficulties are as follows.

Analysis of errors of different methods and an example

The asymptotic analysis in NG-MBS,early-exercise,BIL,asymp-sens,AsAndersenLipton demonstrates that if the characteristic exponent $\psi$ of a L\'evy process $X$ of any popular class is not a rational function, the price of a single barrier option is irregular at the boundary.\footnote{The asymptotic formulas are derived using a couple of general properties of $\psi$; in EfficientAmenable we verify these properties for all popular classes of processes.} If $X$ is of finite variation and the drift points to (from) the boundary, then the option price is of class $C^1$ up to the boundary (is discontinuous at the boundary); option's gamma is unbounded in all cases. If $X$ is of infinite variation or the drift is 0, option's delta tend to infinity as the spot tends to the barrier. The representation of the Laplace transform of the double barrier options as a sum of the Laplace transform of an European option and two perpetual single barrier options derived in MSdouble,BLdouble and, in the dual form, used in the present paper, implies that in all cases, the price of the double barrier option is irregular at one of the barriers, at least. See Fig. (ref) for an example. In particular, in the vicinity of at least one barrier, the option price in a model with an irrational $\psi$ is larger than in a model with a rational characteristic exponent and in any model with non-negligible diffusion component.

If a numerical method uses the time discretization (and Carr's randomization can be interpreted as the time discretization), and the “true" process is replaced by a more convenient one so that the option price is of class $C^1$ up to the boundary, then the errors in the vicinity of a boundary accumulate and propagate. If, at the same time, the number of time steps is large, then the resulting error is quite large. See examples in single which demonstrate sizable errors of the approximation of KoBoL with HEJD: certainly, such an approximation will produces extremely large errors if applied in regime-switching models. An example in Table (ref) demonstrate, if Carr's randomization is used, it is difficult and time consuming to satisfy the moderate error tolerance of the order of 1% unless the time step is very small. Fig. (ref) clearly demonstrates that if the time step is small, then the grid in the state space must be very fine, and the interpolation of order greater than 2 cannot be applied. Note the following important fact evident from Fig. (ref) (the fact can be rigorously proved using the asymptotic analysis tools): the linear interpolation decreases the option price at each step of backward induction, which explains why the Carr's randomization prices in Table (ref) are smaller than the true prices. This effect dominates the well-known theoretical effect of the time discretization: if, at each time step, the price is calculated perfectly or very accurately, the resulting price must be larger than the true price. Depending on the ratio of the time step and step in the state space, either effect can dominate, and one can choose a sequence of refinements of the grids in on the time line and state space which produce a sequence of results converging to any “desirable" price - whether the latter is correct or not.

Errors of approximations of the continuous time model with the discrete time model are larger than the errors of Carr's randomization, and accurate numerical realization of the transition operators for small time steps is very difficult because their kernels have a very high peak or even discontinuous at 0 (see NG-MBS). Convenient approximations of the kernels using cosine functions (COS method) or B-spline approximations (BPROJ method) lead to serious errors. See MarcoDiscBarr for examples of very large errors of COS method used to price single barrier options, and discussion in BSINH on the range of applicability of BPROJ method: since the error of the approximation is in $H^2$-norm, the approximation cannot work for options of very short maturity, for processes close to Variance Gamma especially.

Calculations in the dual space SINHregular,Contrarian,EfficientLevyExtremum allows one to control errors efficiently, and satisfy a small error tolerance using arrays of a moderate size. See Table (ref), where

Note that, in the presence of jumps, the lengths of the grids in the state space needed for accurate calculations are larger than $\ln(H_+/H_-)/\Delta$. If the steepness parameters $\lambda_+,\lambda_-$ are not as large as in the examples (AA), (AB), (MA), (MB) in the paper, then one has to use the grids several times longer and then the CPU time can be much larger than in the example shown in Table (ref). The justification of the Richardson extrapolation for single barrier options in asymp-sens admits a straightforward modification to the case of double barrier options\footnote{The Richardson extrapolation is not applicable to American options although widely used.} but as the last line in Table (ref) shows, the improvement in the accuracy is not large.

Conclusion

In the paper, we constructed a new efficient method to price double barrier options in L\'evy models of finite variation and non-zero drift. The method is a variation of the method in EfficientDoubleBarrier for L\'evy processes of either infinite variation or finite variation and zero drift. We explained why the method in EfficientDoubleBarrier required a modification. In the operator language, the infinitesimal generator $L$ of a process in EfficientDoubleBarrier is a sectorial operator, and in the setting of the present paper $L$ is not sectorial, which makes it impossible to use efficient conformal deformations of all contours of integration in the pricing formula. We use the Gaver-Wynn-Rho algorithm to evaluate the Bromwich integral, and, for each value of the spectral parameter, calculate the price of the corresponding perpetual double barrier options making the calculations in the dual space. We calculate the resulting integrals using the sinh-acceleration technique SINHregular,Contrarian,EfficientLevyExtremum,EfficientDoubleBarrier. The resulting GWR-SINH method is fast but less accurate than the method in EfficientDoubleBarrier because the errors of the GWR algorithm for non-sectorial operators are larger. We discuss the sources of errors of popular groups of methods, which make it very difficult to price barrier options sufficiently accurately and fast so that arbitrage opportunity can be exploited, and present empirical examples, where arbitrage opportunities may present themselves and efficient methods for pricing double barrier options are necessary. A modification of the algorithm formulated in the paper can be used instead of the corresponding block in MSdouble to price double barrier options in the regime-switching models and approximations of stochastic volatility models and models with stochastic interest rates by regime-switching models that we developed earlier ExitRSw,stoch-int-rate-CF,amer-reg-sw-SIAM,SVolSSRN,BLHestonStIR08,BarrStIR.