The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
75,821 characters
Alternative models for FX, arbitrage opportunities and efficient pricing of double barrier options in L\'evy models
\title[Alternative models for FX and double barrier options in L\'evy models]
{Alternative models for FX, arbitrage opportunities and efficient pricing of double barrier options in L\'evy models}
\author[
Svetlana Boyarchenko and
Sergei Levendorski\u{i}]
{
Svetlana Boyarchenko and
Sergei Levendorski\u{i}}
\begin{abstract}
We analyze the qualitative differences between prices of double barrier no-touch options
in the Heston model and pure jump KoBoL model calibrated to
the same set of the empirical data, and discuss the potential for arbitrage opportunities if the correct model
is a pure jump model. We explain and demonstrate with numerical examples that accurate and fast calculations of prices of double barrier options in jump models are extremely difficult using the numerical methods available in the literature.
We develop a new efficient method (GWR-SINH method) based of the Gaver-Wynn-Rho acceleration applied to the Bromwich integral;
the SINH-acceleration and simplified trapezoid rule are used to evaluate perpetual double barrier options for each value of the spectral parameter in GWR-algorithm.
The program in Matlab running on a Mac with moderate characteristics achieves the precision of the order of E-5 and better in several several dozen of milliseconds; the precision E-07 is achievable in about 0.1 sec.
We outline the extension of GWR-SINH method to regime-switching models and models with stochastic parameters and stochastic interest rates.
\end{abstract}
\thanks{
\emph{S.B.:} Department of Economics, The
University of Texas at Austin, 2225 Speedway Stop C3100, Austin,
TX 78712--0301, {\tt [email removed]} \\
\emph{S.L.:}
Calico Science Consulting. Austin, TX.
Email address: {\tt
[email removed]}}
\maketitle
\noindent
{\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
\noindent
{\sc MSC2020 codes:} 60-08,42A38,42B10,44A10,65R10,65G51,91G20,91G60
\tableofcontents
\section{Introduction}\label{s:intro}
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 \cite{CGMY,ChenFengLin,levendorskii-xie-Asian,AsianGamma}). However,
as an extensive empirical study in \cite{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., \cite[Thm. 33.1]{sato})
one can easily prove that in several empirical examples documented in \cite{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 \cite{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 \cite{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 \cite{BLdouble}; the main block of the algorithm is the realization in the dual space of the main block in \cite{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
\cite{BLdouble} for perpetual double barrier options but do the calculations in the dual space. The conformal deformation technique
(sinh-acceleration) developed in \cite{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 \cite{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 \cite{EfficientDoubleBarrier}.
In the result, the error of the GWR method in this paper is 1-2 orders of magnitude larger than that in \cite{EfficientDoubleBarrier},
which agrees with the general characterization of the accuracy of the Gaver method in \cite{AbateValko04,AbValko04b}.
The rest of the paper is organized as follows.
In Sect. \ref{s:auxil}, 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{s:double_state}, we recall the scheme of evaluation of the Laplace transform
\cite{BLdouble} using the calculations in the state space. GWR-SINH method and the explicit algorithm are in Sect. \ref{s:GWR-SINH}.
In Sect. \ref{s:pricing_barr_Levy}, we discuss several popular approaches to pricing barrier options and analyze common sources of errors
of various groups of methods. In Sect. \ref{s:concl}, 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{s: calib_and_DNT}. The GWR-SINH method for regime-switching L\'evy models is in Section \ref{s:reg-switch}.
\section{Auxilliary results}\label{s:auxil}
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{s:double_state} and
\ref{s:GWR-SINH}.
\subsection{Wiener-Hopf factorization}\label{ss:WHF}
This subsection contains basic formulas and results systematically used in a number of publications,
e.g., \cite{NG-MBS,barrier-RLPE,paired,EfficientDiscExtremum,EfficientDoubleBarrier}.
In probability, the
Wiener-Hopf factors are defined as
\begin{equation}\label{defphipm}
\phi^+_q(\xi)={\mathbb E}[e^{i\xi {\bar X}_{T_q}}],\quad \phi^-_q(\xi)={\mathbb E}[e^{i\xi {\underline X}_{T_q}}].
\end{equation}
Functions $\phi^\pm_q(\xi)$ appear in the Wiener-Hopf factorization formula
\begin{equation}\label{whfprob}
\frac{q}{q+\psi(\xi)}=\phi^+_q(\xi)\phi^-_q(\xi),\quad \xi\in{\mathbb R}.
\end{equation}
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 \eqref{whfprob}
\begin{equation}\label{operWHF}
{\mathcal E}_q={\mathcal E}^+_q{\mathcal E}^-_q={\mathcal E}^-_q{\mathcal E}^+_q
\end{equation}
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., \cite{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 \cite{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., \cite{NG-MBS,barrier-RLPE,paired})
the following statements hold.
\medbreak\noindent
I. There exist
$\sigma_-(q)<0<\sigma_+(q)$ such that
\begin{equation}\label{crucial}
q+\psi(\eta)\not\in (-\infty,0],\quad \operatorname{\rm Im}\eta\in (\sigma_-(q),\sigma_+(q)).
\end{equation}
\medbreak\noindent
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)$,
\begin{eqnarray}\label{phip1}
\phi^+_q(\xi)&=&\exp\left[\frac{1}{2\pi i}\int_{\operatorname{\rm Im}\eta=\omega_-}\frac{\xi \ln (1+\psi(\eta)/q)}{\eta(\xi-\eta)}d\eta\right].
\end{eqnarray}
\medbreak\noindent
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))$,
\begin{eqnarray}\label{phim1}
\phi^-_q(\xi)&=&\exp\left[-\frac{1}{2\pi i}\int_{\operatorname{\rm Im}\eta=\omega_+}\frac{\xi \ln (1+\psi(\eta)/q)}{\eta(\xi-\eta)}d\eta\right].
\end{eqnarray}
We can (and will) use \eqref{whfprob} and (\ref{phip1}) to calculate ${\phi^-_q}(\xi)$ on $\{\operatorname{\rm Im}\xi\in (0,\sigma^+_q)\}$,
and \eqref{whfprob} and (\ref{phim1}) to calculate ${\phi^+_q}(\xi)$ on $\{\operatorname{\rm Im}\xi\in (\sigma^-_q,0)\}$.
The explicit algorithm formulated and used in \cite{EfficientDoubleBarrier} for processes of infinite variation
and finite variation processes without drift uses the representations (\ref{phip1})-(\ref{phim1}).
In the case of finite variation processes with positive drift, which we consider in the paper, the following set of formulas derived in \cite{NG-MBS,paired} is more efficient:
\begin{eqnarray}\label{phipq1f}
\phi^{+,0}_q(\xi)&=&\exp\left[-\frac{1}{2\pi i}\int_{\operatorname{\rm Im}\eta=\omega_-}\frac{\xi\ln(1+\psi^0(\eta)/(q-i\mu\eta))}{\eta(\eta-\xi)}d\eta\right],\\\label{phipq1ff}
{\phi^+_q}(\xi)&=&\frac{q}{q-i\mu\xi}\phi^{+,0}_q(\xi),
\\\label{phimq1f}
{\phi^-_q}(\xi)&=&\exp\left[\frac{1}{2\pi i}\int_{\operatorname{\rm Im}\eta=\omega_+}\frac{\xi\ln(1+\psi^0(\eta)/(q-i\mu\eta))}{\eta(\eta-\xi)}d\eta\right].
\end{eqnarray}
The formulas for the case of finite variation processes with negative drift are by symmetry.
\subsection{General classes of L\'evy processes amenable to efficient calculations}\label{ss:gen_eff_Levy}
Essentially all popular classes of L\'evy processes processes enjoy additional properties formalized in
\cite{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_+)\}$.
\begin{defin}\label{def:SINH_reg_proc_1D0}(\cite[Defin 2.1]{EfficientAmenable})
We say that $X$ is a SINH-regular L\'evy process (on ${\mathbb R}$) of order
$\nu$ and type $((\mu_-,\mu_+);{\mathcal C}; {\mathcal C}_+)$
iff
the following conditions are satisfied:
\begin{enumerate}[(i)]
\item
$\nu\in\{0+,1+\}\cup (0,2]$ and either $\mu_-<0\le \mu_+$ or $\mu_-\le 0<\mu_+$;
\item
${\mathcal C}={\mathcal C}_{\gamma_-,\gamma_+}, {\mathcal C}_+={\mathcal C}_{\gamma'_-,\gamma'_+}$, where $\gamma_-<0<\gamma_+$, $\gamma_-\le \gamma'_-\le 0\le \gamma'_+\le \gamma_+$,
and $|\gamma'_-|+\gamma'_+>0$;
\item
the characteristic exponent $\psi$ of $X$ can be represented in the form
\begin{equation}\label{eq:reprpsi}
\psi(\xi)=-i\mu\xi+\psi^0(\xi),
\end{equation}
where $\mu\in{\mathbb R}$, and
$\psi^0$ admits analytic continuation to $i(\mu_-,\mu_+)+ ({\mathcal C}\cup\{0\})$;
\item
for any $\varphi\in (\gamma_-,\gamma_+)$, there exists $c_\infty(\varphi)\in {\mathbb C}\setminus (-\infty,0]$ s.t.
\begin{equation}\label{asympsisRLPE}
\psi^0(\rho e^{i\varphi})\sim c_\infty(\varphi)\rho^\nu, \quad \rho\to+\infty;
\end{equation}
\item
the function $(\gamma_-,\gamma_+)\ni \varphi\mapsto c_\infty(\varphi)\in {\mathbb C}$ is continuous;
\item
for any $\varphi\in (\gamma'_-, \gamma'_+)$, $\operatorname{\rm Re} c_\infty(\varphi)>0$.
\end{enumerate}
\end{defin}
To simplify the constructions in the paper, we assume that $\mu_-<0<\mu_+$ and $\gamma'_-<0<\gamma'_+$.
\begin{example}\label{ex:KoBoL}{\rm In \cite{genBS,KoBoL}, we constructed a family of pure jump processes
generalizing the class of \cite{koponen}, with the L\'evy measure
\begin{equation}\label{KBLmeqdifnu}
F(dx)=c_+e^{\lambda_- x}x^{-\nu_+-1}{\bf 1}_{(0,+\infty)}(x)dx+
c_-e^{\lambda_+ x}|x|^{-\nu_--1}{\bf 1}_{(-\infty,0)}(x)dx,
\end{equation}
where $c_\pm>0, \nu_\pm\in [0,2), \lambda_-<0<\lambda_+$. Starting with \cite{NG-MBS}, we use the name KoBoL processes. If $\nu_\pm\in (0,2), \nu_\pm\neq 1$,
\begin{equation}\label{KBLnupnumneq01}
\psi^0(\xi)=c_+\Gamma(-\nu_+)((-\lambda_-)^{\nu_+}-(-\lambda_--i\xi)^{\nu_+})+c_-\Gamma(-\nu_-)(\lambda_+^{\nu_-}-(\lambda_++i\xi)^{\nu_-}).
\end{equation}
A specialization
$\nu_\pm=\nu\neq 1$, $c=c_\pm>0$, of KoBoL used in a series of numerical examples in \cite{genBS} was named CGMY model in \cite{CGMY} (and the labels were changed:
letters $C,G,M,Y$ replace the parameters $c,\nu,\lambda_-,\lambda_+$ of KoBoL):
\begin{equation}\label{KBLnuneq01}
\psi^0(\xi)= c\Gamma(-\nu)[(-\lambda_-)^{\nu}-(-\lambda_-- i\xi)^\nu+\lambda_+^\nu-(\lambda_++ i\xi)^\nu].
\end{equation}
Evidently, $\psi^0$ given by (\ref{KBLnuneq01}) is analytic in ${\mathbb C}\setminus i((-\infty,\lambda_-]\cup [\lambda_+,+\infty))$, and $\forall\ \varphi\in (-\pi/2,\pi/2)$, (\ref{asympsisRLPE}) holds with
\begin{equation}\label{ascofnupeqnumcc}
c_\infty(\varphi)=-2c\Gamma(-\nu)\cos(\nu\pi/2)e^{i\nu\varphi}.
\end{equation}
}
\end{example}
In \cite{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
\begin{equation}\label{eq:sSLrepr}
\psi(\xi)=(a^+_2\xi^2-ia^+_1\xi)ST({\mathcal G}^0_+)(-i\xi)+(a^-_2\xi^2+ia^-_1\xi)ST({\mathcal G}^0_-)(i\xi)+(\sigma^2/2)\xi^2-i\mu\xi,
\end{equation}
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 \cite{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 \cite{EfficientAmenable}, the reader can find a list of SINH-processes and SL-processes, with calculations of the order and type.
\section{Double barrier options: calculations in the state space}\label{s:double_state}
\subsection{Laplace transform and its inversion}
Let $r\in {\mathbb R}$, $T>0$, $h_-<x<h_+$, and $G\in L_\infty((h_-,h_+))$.
To evaluate
\begin{equation}\label{Vdoublent1}
V(G;r; h_-,h_+;T,x)={\mathbb E}^x[e^{-rT}{\bf 1}_{\tau^-_{h_+}\wedge \tau^+_{h_+}>T}G(X_T)],
\end{equation}
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
\begin{equation}\label{inv_tV0}
V(G; r+r_0; h_-,h_+;T,x)=\frac{1}{2\pi i}\int_{\operatorname{\rm Re} q'=\sigma}e^{q'T} {\tilde V}(G; r+r_0; h_-,h_+;q',x),
\end{equation}
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 \cite{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 \cite{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.
\subsection{Evaluation of the perpetual double barrier options}
The notation and scheme in this Section are borrowed from \cite{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
\begin{equation}\label{eq:repr1}
{\tilde V}(G; h_-,h_+;q,x)=q^{-1}{\mathbb E}^x[{\bf 1}_{\tau^-_{h_+}\wedge \tau^+_{h_+}>T_q}G(X_{T_q})].
\end{equation}
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
\begin{eqnarray}\label{def:tVp1}
{\tilde V}^+_1(lG;h_-,h_+;q,x)={\mathbb E}^x[e^{-q\tau^+_{h_+}}{\tilde V}^0(lG;q, X_{\tau^+_{h_+}})],\
\\
\label{def:tVm1}
{\tilde V}^-_1(lG;h_-,h_+;q,x)={\mathbb E}^x[e^{-q\tau^-_{h_-}}{\tilde V}^0(lG;q, X_{\tau^-_{h_-}})],
\end{eqnarray}
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
\begin{eqnarray}\label{def:tVpj}
{\tilde V}^+_j(lG;h_-,h_+;q,x)={\mathbb E}^x[e^{-q\tau^+_{h_+}}{\tilde V}^-_{j-1}(lG;h_-,h_+;q, X_{\tau^+_{h_+}})],\
\\
\label{def:tVmj}
{\tilde V}^-_j(lG;h_-,h_+;q;x)={\mathbb E}^x[e^{-q\tau^-_{h_-}}{\tilde V}^+_{j-1}(lG;h_-,h_+;q, X_{\tau^-_{h_-}})].
\end{eqnarray}
In \cite{BLdouble}, the following key theorem is derived from the stochastic continuity of $X$.
\begin{thm}\label{thm:convergence}
For any $\sigma>0$ and $h_-<h_+$, there exist $\delta_\pm=\delta_\pm(\sigma, h_+-h_-)\in (0,1)$ such that for all
$q\ge \sigma$, $lG\in L_\infty({\mathbb R})$ and $j=2,3,\ldots$,
\begin{eqnarray}\label{eq:boundtVp}
\sup_{x\ge h_+}|{\tilde V}^+_j(lG;h_-,h_+;q;x)|&\le&\delta_+\sup_{x\le h_-}|{\tilde V}^-_{j-1}(lG;h_-,h_+;q,x)|,\\
\label{eq:boundtVm}
\sup_{x\le h_-}|{\tilde V}^-_j(lG;h_-,h_+;q;x)|&\le&\delta_- \sup_{x\ge h_+}|{\tilde V}^+_{j-1}(lG;h_-,h_+;q,x)|,
\end{eqnarray}
and
\begin{eqnarray}\label{qtV}
{\tilde V}^1(lG;h_-,h_+;q,x)&=&\sum_{j=1}^{+\infty} (-1)^j ({\tilde V}^+_j(lG;h_-,h_+;q,x)+{\tilde V}^-_j(lG;h_-,h_+;q,x)).
\end{eqnarray}
The series on the RHS of (\ref{qtV}) exponentially converges in $L_\infty$-norm.
\end{thm}
Under additional weak conditions on $X$, Theorems 11.4.4 and 11.4.5 in \cite{IDUU} state that
\begin{eqnarray}\label{tVhm1}
{\tilde V}^-_1(lG; h_-,h_+; q,x)&=&q^{-1}({\mathcal E^-_q}{\bf 1}_{(-\infty,h_-]}{\mathcal E^+_q})lG)(x),\\
\label{tVhp1}
{\tilde V}^+_1(lG; h_-,h_+; q,x)&=&q^{-1}({\mathcal E^+_q}{\bf 1}_{[h_+,+\infty)}{\mathcal E^-_q})lG)(x);
\end{eqnarray}
in \cite{single}, (\ref{tVhm1})-(\ref{tVhp1}) are proved for any L\'evy process.
Similar representations for ${\tilde V}^\mp_j$, $j=2,3,\ldots$:
\begin{eqnarray}\label{tVhpj2}
{\tilde V}^+_j(lG; h_-,h_+; q,x)&=&({\mathcal E}^+_q{\bf 1}_{[h_+,+\infty)}({\mathcal E}^+_q)^{-1}){\tilde V}^-_{j-1}(lG;h_-,h_+; q,x),
\\\label{tVhmj2}
{\tilde V}^-_j(lG;h_-,h_+; q,x)&=&({\mathcal E}^-_q{\bf 1}_{(-\infty, h_-]}({\mathcal E}^-_q)^{-1}{\tilde V}^+_{j-1})(lG;h_-,h_+; q,x),
\end{eqnarray}
follow from Theorems 11.4.6 and 11.4.7 in \cite{IDUU}. See \cite{EfficientDoubleBarrier} for details.
\begin{rem}\label{rem:convergence}{\rm
\begin{enumerate}[(a)]
\item
For a numerical realization, the series on the RHS of (\ref{qtV}) are truncated,
and $\sum_{j=1}^\infty$ replaced with $\sum_{j=1}^{M_0}$; given the error tolerance,
$M_0$ can be chosen using the bounds in Theorem \ref{thm:convergence}.
\item If $G$ is bounded on $[h_-,h_+]$ but given by an analytical expression that defines an unbounded function on ${\mathbb R}$, we can use the scheme above replacing $G$ with a function $lG$ of class $L_\infty({\mathbb R})$, which coincides with $G$ on $[h_-,h_+]$. This simple consideration suffices for a numerical realization in the state space.
If the numerical realization is in the dual space, complex-analytical properties of the Fourier transform ${\hat G}$ are crucial.
If ${\hat G}$ has good properties as in the case of the call option, then a good replacement of $G$ with a bounded function
requires additional work. Instead, we will not replace $G$ but assume that the strip of analyticity of $\psi$ is sufficiently wide, and
$q>0$ is sufficiently large so that the series (\ref{qtV}) converges in a space with an appropriate exponential weight.
\item Alternatively, in the case of a call option, one can make an appropriate Esscher transform and reduce to the case of
a bounded $G$.
\item
The inverse Laplace transform of ${\tilde V}^0(q,x)=q^{-1}{\mathcal E_q} G(x)$ is the price $V_{\mathrm{euro}}(G; T,x)$ of the European option
with the payoff $G(x+X_T)$ at maturity date $T$. An explicit procedure for an efficient numerical evaluation of $V_{\mathrm{euro}}(G; T,x)$
can be found in \cite{SINHregular}. In the paper, we design an efficient numerical procedure for the evaluation of
$V^1(G; h_-,h_+;T,x)= V(G; h_-,h_+;T,x)-V_{\mathrm{euro}}(G; T,x)$. Note that $V^1(G; h_-,h_+;T,x)$ is the inverse Laplace transform of the series
on the RHS of (\ref{qtV}).
\item
To shorten the notation, we suppress the dependence of ${\tilde V}^0, {\tilde V}^1$ and ${\tilde V}^\pm_j$ on $lG$.
\end{enumerate}
}
\end{rem}
\section{GWR-SINH method}\label{s:GWR-SINH}
\subsection{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
\begin{enumerate}[(1)]
\item
the curve ${\mathcal L}^+_{\omega_{1,+}, b_+,\omega_+}$ intersects the imaginary axis at a point $ia_+$,
where $a_+\in (0, \lambda_+)$, and $q_0+\psi(ia_+)>0$;
\item
the curve ${\mathcal L}^-_{\omega_{1,-}, b_-,\omega_-}$ intersects the imaginary axis at a point $ia_-$,
where $a_-\in (\lambda_-,0)$, and $q_0+\psi(ia_-)>0$ and $q_0+a_-\mu>0;$
\item
for all $\eta\in {\mathcal L}^+_{\omega_{1,+},b_+,\omega_+}\cup{\mathcal L}^-_{\omega_{1,-},b_-,\omega_-}$, $q_0+\psi(\eta)\not\in (-\infty,0]$.
\end{enumerate}
We deform the contours of integration in (\ref{phipq1f}) and (\ref{phimq1f}) and calculate \begin{eqnarray*}\label{phipq1fsinh}
\phi^{+,0}_q(\xi)&=&\exp\left[-\frac{1}{2\pi i}\int_{{\mathcal L}^-_{\omega_{1,-}, b_-,\omega_-}}\frac{\xi\ln(1+\psi^0(\eta)/(q-i\mu\eta))},
\ \xi\in {\eta(\eta-\xi)}d\eta\right],\ {\mathcal L}^+_{\omega_{1,+}, b_+,\omega_+},
\\\label{phipq1fsinh2}
{\phi^+_q}(\xi)&=&(1-i\mu\xi/q)^{-1}\phi^{+,0}_q(\xi), \ {\mathcal L}^+_{\omega_{1,+}, b_+,\omega_+},
\\\label{phimq1fsinh}
{\phi^-_q}(\xi)&=&\exp\left[\frac{1}{2\pi i}\int_{{\mathcal L}^+_{\omega_{1,+}, b_+,\omega_+}}\frac{\xi\ln(1+\psi^0(\eta)/(q-i\mu\eta))}{\eta(\eta-\xi)}d\eta\right], \ \xi\in{\mathcal L}^-_{\omega_{1,-}, b_-,\omega_-}.
\end{eqnarray*}
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{whfprob}).
\begin{rem}\label{rem:choice}{\rm
In typical cases, $\mu$ and $c\Gamma(-\nu)$ are of the order of $1$ or even 0.1, and $\lambda_+$ and $-\lambda_-$ are greater than 3.
Therefore, it is possible to satisfy (1) with $a_+>2$ and (2) using $a_-<-1$. Typically, the choice of large $a_\pm$ (in absolute value)
is not optimal, and,
in this paper, we choose $a_\pm=\pm 1$, $\omega_{1,\pm}=0$, $\omega_\pm=\pm \pi/4$, $d_\pm=\pi/4$,
and set $b_\pm=1/\sin(\omega_++d_+)=1$ (or smaller multiplying by, e.g. $0.9$). In the program, we first check that, with this choice, the conditions (1) and (2) are satisfied, and then check condition (3).
}
\end{rem}
\subsection{General scheme for a fixed $q>0$}\label{ss:gen_q>0}
The numerical scheme in this section has small but crucial differences from
the scheme in \cite{EfficientDoubleBarrier}. Writing the program, the reader can check that the scheme formulated in
\cite{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 \cite{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,
\begin{equation}\label{htV1p}
{\hat{\tilde V}}^\pm_1(h_-,h_+;q,\xi)=\pm e^{-ih_\pm\xi}\frac{\phi^\pm_q(\xi)}{i\xi q}, \ \xi\in {\mathcal L}^\mp.
\end{equation}
For $j=1,2,\ldots, $ define
\begin{equation}\label{def:hWj}
{\hat W}_j^\pm(h_-,h_+; q,\xi)=qe^{ih_\pm\xi}\phi^\pm_q(\xi)^{-1}{\hat{\tilde V}}_j^\pm(h_-,h_+;q,\xi), \ \xi\in{\mathcal L}^\mp.
\end{equation}
Evidently,
\begin{equation}\label{tVtW}
{\hat{\tilde V}}_j^\pm(h_-,h_+;q,\xi)=q^{-1}e^{-ih_\pm\xi}\phi^\pm_q(\xi){\hat W}^\pm_j(h_-,h_+; q,\xi), \ \xi\in{\mathcal L}^\mp,
\end{equation}
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{tVhpj2})-(\ref{tVhmj2}) in terms of $ {\hat W}^\pm_j$ as follows.
For $\xi\in {\mathcal L}^-$, we have
\begin{eqnarray*}
{\hat W}^+_j(h_-,h_+; q,\xi)&=&qe^{ih_+\xi}{\mathcal F}_{x\to \xi}{\bf 1}_{[h_+,+\infty)}({\mathcal E}^+_q)^{-1}{\tilde V}^-_{j-1}(h_-,h_+;q,x)
\\
&=& e^{ih_+\xi}\int_{h_+}^{+\infty}dy\, e^{-iy\xi}
\frac{1}{2\pi}\int_{{\mathcal L}^+}e^{i(y-h_-)\eta}\frac{\phi^-_q(\eta)}{\phi^+_q(\eta)}{\hat W}^-_{j-1}(h_-,h_+;q,\eta)d\eta\\
&=&-\frac{e^{ih_+\xi}}{2\pi i}\int_{{\mathcal L}^+}\frac{e^{-ih_+(\xi-\eta)-ih_-\eta}}{\eta-\xi}
\frac{\phi^-_q(\eta)}{\phi^+_q(\eta)}{\hat W}^-_{j-1}(h_-,h_+;q,\eta)d\eta
\end{eqnarray*}
(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,
\begin{equation}\label{hWp}
{\hat W}^+_j(h_-,h_+; q,\xi)=-\frac{1}{2\pi i}\int_{{\mathcal L}^+}\frac{e^{i(h_+-h_-)\eta}}{\eta-\xi}
\frac{\phi^-_q(\eta)}{\phi^+_q(\eta)}{\hat W}^-_{j-1}(h_-,h_+;q,\eta)d\eta,\ \xi\in{\mathcal L}^-.
\end{equation}
Similarly, for $\xi\in {\mathcal L}^+$, we calculate
\begin{eqnarray*}
{\hat W}^-_{j}(h_-,h_+; q,\xi)&=& e^{ih_-\xi}\int_{-\infty}^{h_-}dy\, e^{-iy\xi}
\frac{1}{2\pi}\int_{{\mathcal L}^-}e^{i(y-h_+)\eta}\frac{\phi^+_q(\eta)}{\phi^-_q(\eta)}{\hat W}^+_{j-1}(h_-,h_+;q,\eta)d\eta\\
&=&\frac{e^{ih_-\xi}}{2\pi i}\int_{{\mathcal L}^-}\frac{e^{-ih_-(\xi-\eta)-ih_+\eta}}{\eta-\xi}
\frac{\phi^+_q(\eta)}{\phi^-_q(\eta)}{\hat W}^+_{j-1}(h_-,h_+;q,\eta)d\eta.
\end{eqnarray*}
Simplifying,
\begin{equation}\label{hWm}
{\hat W}^-_{j}(h_-,h_+; q,\xi)=\frac{1}{2\pi i}\int_{{\mathcal L}^-}\frac{e^{-i(h_+-h_-)\eta}}{\eta-\xi}
\frac{\phi^+_q(\eta)}{\phi^-_q(\eta)}{\hat W}^+_{j-1}(h_-,h_+;q,\eta)d\eta,\ \xi\in {\mathcal L}^+.
\end{equation}
In the cycle $j=2,3\ldots,$ we calculate ${\hat W}^\pm_{j}(h_-,h_+; q,\xi)$ and partials sums of the series
\begin{equation}\label{Winf}
{\hat W}^\pm(h_-,h_+; q,\xi)=\sum_{j=1}^\infty (-1)^j {\hat W}^\pm_j(h_-,h_+; q,\xi), \ \xi\in {\mathcal L}^\mp.
\end{equation}
Finally, for $x\in (h_-, h_+)$, we calculate
\begin{equation}\label{tVqpm}
{\tilde V}^1_\pm(lG;h_-,h_+;T,x)=\frac{1}{2\pi}\int_{{\mathcal L}^\mp}e^{i(x-h_\pm)\xi}\phi^\pm_q(\xi){\hat W}^\pm(h_-,h_+; q,\xi)d\xi,
\end{equation}
then
\begin{equation}\label{tVqntfin}
V^1(lG;h_-,h_+;T,x)=\frac{1}{2\pi i}\int_{\operatorname{\rm Re} q=\sigma}dq\,\frac{e^{qT}}{q}({\tilde V}^1_+(lG;h_-,h_+;T,x)+
{\tilde V}^1_-(lG;h_-,h_+;T,x)),
\end{equation}
and
\begin{equation}\label{eq:double_final}
V(G;h_-,h_+;T,x)=V_{\mathrm{euro}}(lG;T,x)+V^1(lG;h_-,h_+;T,x).
\end{equation}
\begin{rem}\label{rem:convergence2}{\rm Theorem \ref{thm:convergence} implies that if
the Gaver-Stehfest method or GWR algorithm is used to numerically evaluate
the Bromwich integral, then the inverse Fourier transform of the series on the RHS of (\ref{tVqntfin}) converges for any $h_-<x<h_+$
and any $q>0$ used in the algorithm. It can be proved that the infinite sums of the Fourier transforms converge as well.
}
\end{rem}
\subsection{
Efficient numerical evaluation of the series (\ref{Winf})
}\label{sss:effhWseriesI}
Define operators\\ ${\mathcal K}_{-+}(={\mathcal K}_{-+}(q; {\mathcal L}^+; h_-,h_+))$ and ${\mathcal K}_{+-}(={\mathcal K}_{+-}(q; {\mathcal L}^-; h_-,h_+))$ by
\begin{eqnarray}\label{defKmp}
{\mathcal K}_{-+}{\hat u}(\xi)&=&\frac{1}{2\pi }\int_{{\mathcal L}^+}\frac{e^{i(h_+-h_-)\eta}(1-i\mu\eta/q)}{\eta-\xi}
\frac{\phi^-_q(\eta)}{\phi^{+,0}_q(\eta)}{\hat u}(\eta)d\eta,\ \xi\in {\mathcal L}^-,\\
\label{defKpm}
{\mathcal K}_{+-}{\hat u}(\xi)&=&\frac{1}{2\pi }\int_{{\mathcal L}^-}\frac{e^{-i(h_+-h_-)\eta}}{\eta-\xi}
\frac{\phi^+_q(\eta)}{\phi^-_q(\eta)}{\hat u}(\eta)d\eta,\ \xi\in {\mathcal L}^+.
\end{eqnarray}
\begin{rem}\label{rem:Kmp_nu<1_mu>0} {\rm In \cite{EfficientDoubleBarrier}, the formula for ${\mathcal K}_{-+}$ is of the same form as
the one for ${\mathcal K}_{+-}$: just change the signs. However, in the case $\nu<1$ and $\mu>0$, the ratio $\phi^-_q(\eta)/\phi^+_q(\eta)$
increases as $\eta\to \infty$ very fast, and the algorithm (as in \cite{EfficientDoubleBarrier}) which explicitly uses this ratio
becomes unstable. The factor $e^{i(h_+-h_-)\eta}(1-i\mu\eta/q)$ is uniformly bounded and decays fast at infinity,
and the factor $\phi^-_q(\eta)/\phi^{+,0}_q(\eta)$ is bounded.}
\end{rem}
We write (\ref{hWp}) and (\ref{hWm}) as
\begin{equation}\label{hWpmj}
{\hat W}^+_{j+1}=i{\mathcal K}_{-+}{\hat W}^-_j,\ {\hat W}^-_{j+1}=-i{\mathcal K}_{+-}{\hat W}^+_j, j=1,2,\ldots
\end{equation}
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{rem:convergence}. 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 \cite{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:
\begin{eqnarray*}\label{Kmpjik}
{\mathcal K}_{-+}&=&\frac{\zeta^+b^+}{2\pi}\left[\frac{e^{i(h_+-h_-)\xi^+_k}(1-i\mu \xi^+_k/q)}{\xi^+_k-\xi^-_j}
\cdot\frac{{\phi^-_q}(\xi^+_k)}{\phi^{+,0}_q(\xi^+_k)}\cosh(i\omega^++y^+_k)\right]_{|j|\le N^-, |k|\le N^+},
\\\label{Kpmjik}
{\mathcal K}_{+-}&=&\frac{\zeta^-b^-}{2\pi}\left[\frac{e^{-i(h_+-h_-)\xi^-_k}}{\xi^-_k-\xi^+_j}
\cdot\frac{{\phi^+_q}(\xi^-_k)}{{\phi^-_q}(\xi^-_k)}\cosh(i\omega^-+y^-_k)\right]_{|j|\le N^+, |k|\le N^-}.
\end{eqnarray*}
\begin{rem}\label{rem:matrix_inv} {\rm If the sizes of matrices ${\mathcal K}_{-+}$ and ${\mathcal K}_{+-}$, hence,
the sizes of matrices ${\mathcal K}^+:={\mathcal K}_{-+}{\mathcal K}_{+-}, {\mathcal K}^-:={\mathcal K}_{+-}{\mathcal K}_{-+}$ are moderate so that the inverse matrices $(I-{\mathcal K}^\pm)^{-1}$
can be efficiently calculated, one can calculate infinite sums using standard matrix tools. See \cite{EfficientDoubleBarrier}.
}
\end{rem}
\subsection{Algorithm of GWR-SINH method for DNT options}\label{ss:algo_GWR-SINH_DNT}
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.
\begin{enumerate}[Step I.]
\item
Choose $r_0$ and the order of the GWR algorithm, and calculate the set of nodes $Q_{GWR}$ and weights of the GWR algorithm.
\item
{\em Grids for the approximations of ${\hat W}^\pm_m$ and operators ${\mathcal K}_{-+}, {\mathcal K}_{+-}, {\mathcal K}^\pm$}. Choose the sinh-deformations and grids for the simplified trapezoid rule on ${\mathcal L}^\pm$: $\vec{y^\pm}:=\zeta^\pm*(-N^\pm:1:N^\pm)$,
$\vec{\xi^\pm}:=i*\omega_1^\pm+ b^\pm*\sinh(i*\omega^\pm+i\vec{y^\pm})$. Calculate $\psi^\pm:=\psi(\vec{\xi^\pm})$ and
$\vec{der^\pm}:=b^\pm*\cosh(i*\omega^\pm+\vec{y^\pm}).
$
\item
{\em Grids for evaluation of the Wiener-Hopf factors $\phi^\pm_q(\xi)$.} Choose longer and finer grids for the simplified trapezoid rule on ${\mathcal L}^\pm_1$: $\vec{y^\pm_1}=\zeta_1^\pm*(-N^\pm_1:1:N^\pm_1)$,
$\vec{\xi^\pm_1}:=i*\omega^\pm_1+ b^\pm_1*\sinh(i*\omega_1^\pm+i\vec{y^\pm_1})$. Calculate $\psi^{0,\pm}_1:=\psi^0(\vec{\xi^\pm_1})$ and
$\vec{der^\pm_1}:=b^\pm_1*\cosh(i*\omega_1^\pm+\vec{y^\pm_1}).
$
\item
Calculate 2D arrays
\begin{eqnarray*}
D^{+-}_1&:=&1./(\mathrm{conj}(\vec{\xi^-_1})'*\mathrm{ones}(1,2*N^++1)-\mathrm{ones}(2*N^-_1+1,1)*\vec{\xi^+})),\\
D^{-+}_1&:=&1./(\mathrm{conj}(\vec{\xi^+_1})'*\mathrm{ones}(1,2*N^-+1)-\mathrm{ones}(2*N^+_1+1,1)*\vec{\xi^-})),\\
D^{+-}&:=&1./(\mathrm{conj}(\vec{\xi^-})'*\mathrm{ones}(1,2*N^++1)-\mathrm{ones}(2*N^-+1,1)*\vec{\xi^+})),\\
D^{-+}&:=&1./(\mathrm{conj}(\vec{\xi^+})'*\mathrm{ones}(1,2*N^-+1)-\mathrm{ones}(2*N^++1,1)*\vec{\xi^-})).
\end{eqnarray*}
\item
{\em Calculate ${\hat W}^+_1=-i./\vec{\xi^-}, {\hat W}^-_1=i./\vec{\xi^+}$.}
\item
{\em Calculate} $\vec{\phi^{+,0}_q}=\phi^{0,+}_q(\vec{\xi^+})$ and $\vec{{\phi^-_q}}={\phi^-_q}(\vec{\xi^-})$:
\begin{eqnarray*}
\Phi^\pm &=&1+\psi^{0,\pm}_1./(q-i*\mu*\vec{\xi^\pm_1}), \\
\vec{\phi^{0,+}_q}&:=&\exp((\zeta^-_1*i/(2*\pi))*\vec{\xi^+_1}./.*((\log(\Phi^-).*\vec{\xi^-_1}.*\vec{der^-_1})*D^{-+}_1)),\\
\vec{{\phi^-_q}}&:=&\exp(-(\zeta^+_1*i/(2*\pi))*\vec{\xi^-_1}.*((\log(\Phi^+)*\vec{\xi^-_1}.*\vec{der^+_1}.*D^{+-}_1)),
\end{eqnarray*}
and then
${\phi^-_q}(\vec{\xi^+}):=1./\Phi^+./\vec{{\phi^+_q}},\ \phi^{+,0}_q(\vec{\xi^-}):=1./\Phi^-./\vec{{\phi^-_q}},$\\ $ \phi^{+}_q(\vec{\xi^-})=\phi^{+,0}_q(\vec{\xi^-})./(1+\psi^{0,-}./(1-i*(\mu/q)*\vec{\xi^-})).
$
\item
{\em Calculate matrices ${\mathcal K}_{+-}, {\mathcal K}_{+-}$:}
\begin{eqnarray*}
{\mathcal K}_{-+}&:=&(\zeta^+/(2*\pi))*\mathrm{diag}(\vec{der^+}.*({\phi^-_q}(\vec{\xi^+})./\phi^{+,0}_q(\vec{\xi^+}))...\\
&& .*(\exp(i*(h_+-h_-)*\vec{\xi^+}).*(1-i*(\mu/q)*\vec{\xi^+}))*D^{-+};\\
{\mathcal K}_{+-}&:=&(\zeta^-/(2*\pi))*\mathrm{diag}(\vec{der^-}.*({\phi^+_q}(\vec{\xi^-})./{\phi^-_q}(\vec{\xi^-})).*\exp(-i*(h_+-h_-)*\vec{\xi^-}))*D^{+-}.\\
\end{eqnarray*}
\item
Use one of the following three blocks. Calculations using Block (1) and either Block (2) or Block (3) can be used to check the accuracy of the result.
\begin{enumerate}[(1)]
\item\begin{itemize}
\item
Assign ${\hat U}^\pm_1=-{\hat W}^\pm_1,$ ${\hat W}^\pm:={\hat U}^\pm_1$.
\item
In the cycle $j=1,2,\ldots,M_0$, calculate
$
{\hat U}^+_{2}:=-i*{\hat U}^-_1*{\mathcal K}_{-+},\\ {\hat U}^-_{2}:=i*{\hat U}^+_1*{\mathcal K}_{+-},\
{\hat U}^\pm:={\hat W}^\pm+{\hat U}^\pm_2, \ {\hat U}^\pm_1:={\hat U}^\pm_2$.
\item
Assign ${\hat W}^\pm={\hat U}^\pm$.
\end{itemize}
\item
Calculate
\begin{itemize}
\item
${\mathcal K}^-:={\mathcal K}_{-+}*{\mathcal K}_{+-},
{\mathcal K}^+:={\mathcal K}_{+-}*{\mathcal K}_{-+}$;
\item
${\hat W}^+_2:=i*{\hat W}^-_1*{\mathcal K}_{-+}, {\hat W}^-_2:=-i*{\hat W}^+_1*{\mathcal K}_{+-}$;
\item
$
{\hat W}^{\pm,0}:={\hat W}^\pm_2-{\hat W}^\pm_1$;
\item
inverse matrices
$(I-{\mathcal K}^\pm)^{-1}$;
\item
${\hat W}^\pm={\hat W}^{\pm,0}*(I-{\mathcal K}^\mp)^{-1}$.
\end{itemize}
\item
Replace the last two steps of Block (2) with
\[
{\hat W}^\pm=\mathrm{conj}(\mathrm{linsolv}(\mathrm{diag}(\mathrm{ones}(2*N^\mp+1))-\mathrm{conj}({\mathcal K}^\mp)', \mathrm{conj}({\hat W}^{\pm,0})'))'.
\]
Typically, the program with Block (1) achieves precision of the order of E-15 with $M_0=9$ or even $M_0=8$.
The CPU time is several times smaller than with Block (2). Program with Block (2) is faster if
the digitals or vanillas for many strikes need to be calculated. Block (3) is twice faster than Block (2)
if applied only once.
\end{enumerate}
\item
{\em For $x\in (h_-,h_+)$, calculate}
\begin{eqnarray*}
V^+&=&(\zeta^-/(2*\pi))*\sum({\hat W}^+.*\exp(i*(x-h_+)*\vec{\xi^-}).*{\phi^+_q}(\vec{\xi^-}).*\vec{der^-}), \\
V^-&=&(\zeta^+/(2*\pi))*\sum({\hat W}^-.*\exp(i*(x-h_-)*\vec{\xi^+}).*{\phi^-_q}(\vec{\xi^+}).*\vec{der^+});
\end{eqnarray*}
this step can be easily parallelized for a given array $\{x_j\}$.
\item
Use the GWR algorithm to evaluate $V^1$. \item
{\em Final step.} Set $V=e^{r_0T}(1+V^1)$.
\end{enumerate}
\subsection{A numerical example: Table~\ref{Table 3}}\label{ss:numer_example}
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.
\begin{table}
\caption{\small Prices of the double no-touch options, and errors (rounded) of Carr's randomization algorithm \cite{BLdouble}
and GWR-SINH method. Dependence on the spot. Barriers: $H_-=0.95, H_+=1.95$, time to maturity $T=0.25$.
KoBoL model with the parameters in Table \ref{Table 2}, MB: $\nu=0.445, c=1.125, \lambda_+= 27.93, \lambda_-= -51.66, \mu= 0.0940$.}
\begin{tabular}{c|ccccc|c}
\hline
$S$ & 0.96 & 0.98 & 1 & 1.02 & 1.04 & Time \\
BB & 0.4325056 & 0.6497429 & 0.6801758 & 0.5289720 & 0.224546 & 812
\\\hline
$ \epsilon_0$ &
-1.07E-08 & 1.82E-10 & -9.46E-09 & 3.68E-09 & -3.10E-09 & 767.2\\
$ \epsilon_1$ & 4.84E-06 & -7.86E-06 & -5.83E-06 & -4.09E-05 & 9.98E-05 & 615 \\
$\epsilon_2$ & 5.34E-06 & -3.82E-06 & -1.82E-05 & -6.20E-06 & -6.43E-05 & 70.1\\
$\epsilon_3$ & 5.62E-06 & -5.68E-06 & -2.01E-05 & -9.23E-06 & -3.44E-05 & 69.6\\
$\epsilon_4$ & 1.66E-04 & -4.30E-06& -4.66E-05 &
-7.44E-05 & -3.77E-05 & 29.5
\\\hline
$\epsilon_{C1}$ & -0.0015 & -0.0017 & -0.0027 & -0.0045 & -0.0074 & 1,564\\
$\epsilon_{C2}$ &-0.0024 & -0.0029& -0.00492& -0.0084& -0.013 & 546.3.3\\
$\epsilon_{C3}$ &-0.010& -0.0080& -0.01& -0.027& -0.0103 & 865.3\\\hline
$\epsilon_{Rich}$ & -0.00058 & -0.00050 & -0.00043 & -0.00063& -0.0015 & 2,156.2
\\\hline
\end{tabular}
\begin{flushleft}{\tiny
Time: CPU time in msec, average of 100 runs\\
BB: GWR-SINH(0;276,502): $r_0=0$, 276 points on ${\mathcal L}^\pm$ in the main cycle and 502 points to evaluate the Wiener-Hopf
factors
\\
$\epsilon_0$: differences with GWR-SINH(0;306,557) with different omegas\\
$\epsilon_1,\epsilon_2, \epsilon_3,\epsilon_4$: differences with GWR-SINH(5;276,502), GWR-SINH(0;108,188), GWR-SINH(0.5;88,155), GWR-SINH(-0.5;55,100)\\
$\epsilon_{C1}, \epsilon_{C2}, \epsilon_{C3}$: errors of Carr randomization method \cite{BLdouble} with $(N_T,\Delta_x)=(600,1/8000)$,
$=(300,1/4000)$, $=(150,1/2000)$.
\\ $\epsilon_{Rich}$: errors of Richardson's extrapolation of C1 and C2 prices.
}
\end{flushleft}
\label{Table 3}
\end{table}
\section{Pricing barrier options in L\'evy models: general discussion and examples}\label{s:pricing_barr_Levy}
\subsection{A short review of the literature} The general formulas for single barrier options with continuous monitoring were derived in \cite{KoBoL,barrier-RLPE,NG-MBS}
using the operator form of the Wiener-Hopf factorziation \cite{eskin}, under certain regularity conditions on the characteristic exponent.
In \cite{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 \cite{kou,KW1} and its generalization: Hyper-exponential jump-diffusion model (HEJD model) constructed independently in \cite{lipton-columbia,lipton-risk} (see also \cite{LiptonSelection}) and
\cite{amer-put-levy-maphysto,amer-put-levy}.
In \cite{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 \cite{lipton-sepp}. In
\cite{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 \cite{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 \cite{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 \cite{paraLaplace}). The performance of the GS-algorithm can be improved using Wynn-Rho acceleration (GWR algorithm) - see \cite{AbateValko04} and,
for examples in the context of pricing options of long maturity, \cite{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{ss:irr} 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 \cite{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., \cite{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.
\subsection{Analysis of errors of different methods and an example}\label{ss:irr}
The asymptotic analysis in \cite{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 \cite{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 \cite{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{fig:shortmaturity} 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 \cite{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{Table 3} 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{fig:shortmaturity} 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{fig:shortmaturity} (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{Table 3} 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 \cite{NG-MBS}). Convenient approximations of the kernels using cosine functions (COS method) or B-spline approximations (BPROJ method) lead to serious errors. See \cite{MarcoDiscBarr} for examples of very large errors of COS method used to price
single barrier options, and discussion in \cite{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 \cite{SINHregular,Contrarian,EfficientLevyExtremum} allows one to control errors efficiently,
and satisfy a small error tolerance using arrays of a moderate size. See Table \ref{Table 3}, where
\begin{itemize}
\item
the number of nodes in the GWR algorithm is 16;
\item GWR-SINH($r_0$, $a$,$b$) means that the grids on the curves ${\mathcal L}^\pm$ in the dual space are of length $a$ in the main cycle, and length $b$ for the evaluation of the Wiener-Hopf factors;
\item
$N$ is the number of time steps in Carr's randomization algorithm, and $\Delta$ is the step of the grid in the $x=\ln S$ space.
\end{itemize}
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{Table 3}. The justification of the Richardson extrapolation for single barrier options in \cite{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{Table 3} shows, the improvement in the accuracy is not large.
\begin{figure}
\centering
\begin{subfigure}[t]{0.45\textwidth}
{\includegraphics[width=\textwidth]{RuvimAAT0004.pdf}}\caption{$T=0.004$}
\end{subfigure}
\begin{subfigure}[t]{0.45\textwidth}
{\includegraphics[width=\textwidth]{GraphAAT001.pdf}}\caption{$T=0.01$}
\end{subfigure}
\caption{Prices of the double barrier option in KoBoL model. Parameters in Table \ref{Table 2} (AA). Barriers:
$H_-=0.95, H_+=1.05$. }\label{fig:shortmaturity}
\end{figure}
\section{Conclusion}\label{s:concl}
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 \cite{EfficientDoubleBarrier} for L\'evy processes of
either infinite variation or finite variation and zero drift. We explained why the method in \cite{EfficientDoubleBarrier}
required a modification. In the operator language, the infinitesimal generator $L$ of a process in \cite{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 \cite{SINHregular,Contrarian,EfficientLevyExtremum,EfficientDoubleBarrier}. The resulting GWR-SINH method is fast but less accurate
than the method in \cite{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 \cite{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 \cite{ExitRSw,stoch-int-rate-CF,amer-reg-sw-SIAM,SVolSSRN,BLHestonStIR08,BarrStIR}.
\begin{thebibliography}{10}
\bibitem{AbateValko04}
J.~Abate and P.P. Valko.
\newblock Multi-precision {L}aplace inversion.
\newblock {\em International Journal of Numerical Methods in Engineering},
60:979--993, 2004.
\bibitem{AsAndersenLipton}
L.B. Andersen and A.~Lipton.
\newblock Asymptotics for exponential {L}\'evy processes and their volatility
smile: Survey and new results.
\newblock {\em International Journal of Theoretical and Applied Finance}, 16,
February 2013.
\newblock Available at SSRN: http://ssrn.com/abstract=2095654.
\bibitem{AAP}
S.~Asmussen, F.~Avram, and M.R. Pistorius.
\newblock Russian and {A}merican put options under exponential phase-type
{L}\'evy models.
\newblock {\em Stochastic Processes and their Applications}, 109(1):79--111,
2004.
\bibitem{AKP}
F.~Avram, A.~Kyprianou, and M.R. Pistorius.
\newblock Exit problems for spectrally negative {L}\'evy processes and
applications to ({C}anadized) {R}ussian options.
\newblock {\em Annals of Applied Probability}, 14(2):215--238, 2004.
\bibitem{MSdouble}
M.~Boyarchenko and S.~Boyarchenko.
\newblock Double barrier options in regime-switching hyper-exponential
jump-diffusion models.
\newblock {\em International Journal of Theoretical and Applied Finance},
14(7):1005--1044, 2011.
\bibitem{BIL}
M.~Boyarchenko, M.~de~Innocentis, and S.~Levendorski\u{i}.
\newblock Prices of barrier and first-touch digital options in {L}\'evy-driven
models, near barrier.
\newblock {\em International Journal of Theoretical and Applied Finance},
14(7):1045--1090, 2011.
\newblock Available at SSRN: http://papers.ssrn.com/abstract=1514025.
\bibitem{single}
M.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Prices and sensitivities of barrier and first-touch digital options
in {L}\'evy-driven models.
\newblock {\em International Journal of Theoretical and Applied Finance},
12(8):1125--1170, December 2009.
\bibitem{BLdouble}
M.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Valuation of continuously monitored double barrier options and
related securities.
\newblock {\em Mathematical Finance}, 22(3):419--444, July 2012.
\bibitem{genBS}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Generalizations of the {B}lack-{S}choles equation for truncated
{L}\'evy processes.
\newblock Working Paper, University of Pennsylvania, April 1999.
\bibitem{KoBoL}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Option pricing for truncated {L}\'evy processes.
\newblock {\em International Journal of Theoretical and Applied Finance},
3(3):549--552, July 2000.
\bibitem{barrier-RLPE}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Barrier options and touch-and-out options under regular {L}\'evy
processes of exponential type.
\newblock {\em Annals of Applied Probability}, 12(4):1261--1298, 2002.
\bibitem{NG-MBS}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock {\em Non-{G}aussian {M}erton-{B}lack-{S}choles {T}heory}, volume~9 of
{\em Adv. Ser. Stat. Sci. Appl. Probab.}
\newblock World Scientific Publishing Co., River Edge, NJ, 2002.
\bibitem{BLSIAM02}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Perpetual {A}merican options under {L}\'evy processes.
\newblock {\em SIAM Journal on Control and Optimization}, 40(6):1663--1696,
2002.
\bibitem{SVolSSRN}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock American options in {L}\'evy models with stochastic volatility, 2007.
\newblock Available at SSRN: http://ssrn.com/abstract=1031280.
\bibitem{IDUU}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock {\em Irreversible {D}ecisions {U}nder {U}ncertainty ({O}ptimal
{S}topping {M}ade {E}asy)}.
\newblock Springer, Berlin, 2007.
\bibitem{ExitRSw}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Exit problems in regime-switching models.
\newblock {\em Journ. of Mathematical Economics}, 44(2):180--206, 2008.
\bibitem{stoch-int-rate-CF}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock American options in {L}\'evy models with stochastic interest rates.
\newblock {\em Journal of Computational Finance}, 12(4):1--30, Summer 2009.
\bibitem{amer-reg-sw-SIAM}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock American options in regime-switching models.
\newblock {\em SIAM Journal on Control and Optimization}, 48(3):1353--1376,
2009.
\bibitem{paraLaplace}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Efficient {L}aplace inversion, {W}iener-{H}opf factorization and
pricing lookbacks.
\newblock {\em International Journal of Theoretical and Applied Finance},
16(3):1350011 (40 pages), 2013.
\newblock Available at SSRN: http://ssrn.com/abstract=1979227.
\bibitem{BarrStIR}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Efficient pricing barrier options and {C}{D}{S} in {L}\'evy models
with stochastic interest rate.
\newblock {\em Mathematical Finance}, 27(4):1089--1123, 2017.
\newblock DOI: 10.1111/mafi.12121.
\bibitem{SINHregular}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Sinh-acceleration: Efficient evaluation of probability distributions,
option pricing, and {M}onte-{C}arlo simulations.
\newblock {\em International Journal of Theoretical and Applied Finance},
22(3):1950--011, 2019.
\newblock DOI: 10.1142/S0219024919500110. Available at SSRN:
https://ssrn.com/abstract=3129881 or http://dx.doi.org/10.2139/ssrn.3129881.
\bibitem{Contrarian}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Static and semi-static hedging as contrarian or conformist bets.
\newblock {\em Mathematical Finance}, 3(30):921--960, 2020.
\newblock Available at SSRN: https://ssrn.com/abstract=3329694 or
http://arxiv.org/abs/1902.02854.
\bibitem{EfficientDoubleBarrier}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Efficient evaluation of double barrier options and joint cpdf of a
{L}\'evy process and its two extrema.
\newblock Working paper, October 2022.
\newblock Available at SSRN: http://ssrn.com/abstract=4262396 or
http://arxiv.org/abs/2211.07765.
\bibitem{EfficientLevyExtremum}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Efficient evaluation of expectations of functions of a {L}\'evy
process and its extremum.
\newblock Working paper, June 2022.
\newblock Available at SSRN: https://ssrn.com/abstract=4140462 or
http://arXiv.org/abs/4362928.
\bibitem{EfficientDiscExtremum}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock Efficient inverse $z$-transform and pricing barrier and lookback
options with discrete monitoring.
\newblock Working paper, July 2022.
\newblock Available at SSRN: https://ssrn.com/abstract=4155587 or
https://doi.org/10.48550/arXiv.2207.02858.
\bibitem{EfficientAmenable}
S.~Boyarchenko and S.~Levendorski\u{i}.
\newblock {L}\'evy models amenable to efficient calculations.
\newblock Working paper, June 2022.
\newblock Available at SSRN: https://ssrn.com/abstract=4116959 or
http://arXiv.org/abs/4339862.
\bibitem{BSINH}
S.~Boyarchenko, S.~Levendorski\u{i}, J.L. Kirkby, and Z.~Cui.
\newblock {S}{I}{N}{H}-acceleration for {B}-spline projection with option
pricing applications.
\newblock {\em International Journal of Theoretical and Applied Finance},
08(24):2150042, 2021.
\newblock Available at SSRN: https://ssrn.com/abstract=3921840 or
arXiv:2109.08738.
\bibitem{BLHestonStIR08}
S.~Boyarchenko and S.Levendorski\u{i}.
\newblock American options in the {H}eston model with stochastic interest rate
and its generalizations.
\newblock {\em Appl. Mathem. Finance}, 20(1):26--49, 2013.
\bibitem{CGMY}
P.~Carr, H.~Geman, D.B. Madan, and M.~Yor.
\newblock The fine structure of asset returns: an empirical investigation.
\newblock {\em Journal of Business}, 75:305--332, 2002.
\bibitem{ChenFengLin}
Z.~Chen, L.~Feng, and X.~Lin.
\newblock Simulation of {L}\'evy processes from their characteristic functions
and financial applications.
\newblock {\em ACM Transactions on Modeling and Computer Simulation}, 22(3),
2011.
\newblock Available at: http://ssrn.com/abstract=1983134.
\bibitem{CV}
R.~Cont and E.~Voltchkova.
\newblock A finite difference scheme for option pricing in jump diffusion and
exponential {L}\'evy models.
\newblock {\em SIAM Journal on Numerical Analysis}, 43(4):1596--1626, October
2005.
\bibitem{dahlbokum07}
A.~Dahlbokum.
\newblock Empirischer {V}ergleich von {O}ptionspreismodellen auf {B}asis
zeitdeformierter {L}\'evy-prozesse: {K}alibrierung, {H}edging,
{M}odellrisiko., 2007.
\bibitem{SIAM2010}
M.~de~Innocentis and S.~Levendorski\u{i}.
\newblock A novel modelling of stochastic skewness with applications to fitting
and hedging.
\newblock SIAM Conference on Financial Mathematics and Engineering, November
2010.
\bibitem{MarcoDiscBarr}
M.~de~Innocentis and S.~Levendorski\u{i}.
\newblock Pricing discrete barrier options and credit default swaps under
{L}\'evy processes.
\newblock {\em Quantitative Finance}, 14(8):1337--1365, 2014.
\newblock Available at: DOI:10.1080/14697688.2013.826814.
\bibitem{eskin}
G.I. Eskin.
\newblock {\em Boundary {V}alue {P}roblems for {E}lliptic {P}seudodifferential
{E}quations}, volume~9 of {\em Transl. Math. Monogr.}
\newblock American Mathematical Society, Providence, RI, 1981.
\bibitem{feng-linetsky09}
L.~Feng and V.~Linetsky.
\newblock Computing exponential moments of the discrete maximum of a {L}\'evy
process and lookback options.
\newblock {\em Finance and Stochastics}, 13(4):501--529, 2009.
\bibitem{FusaiGermanoMarazzina}
G.~Fusai, G.~Germano, and D.~Marazzina.
\newblock Spitzer identity, {W}iener-{H}opf factorization and pricing of
discretely monitored exotic options.
\newblock {\em European Journal of Operational Research}, 251(1):124--134,
2016.
\newblock DOI:10.1016/j.ejor.2015.11.027.
\bibitem{GSh}
X.~Guo and L.A. Shepp.
\newblock Some optimal stopping problems with nontrivial boundaries for pricing
exotic options.
\newblock {\em J.Appl. Probability}, 38(3):647--658, 2001.
\bibitem{HaislipKaishev14}
G.G. Haislip and V.K. Kaishev.
\newblock Lookback option pricing using the {F}ourier transform {B}-spline
method.
\newblock {\em Quantitative Finance}, 14(5):789--803, 2014.
\bibitem{kirkbyJCompFinance18}
J.L. Kirkby.
\newblock American and {E}xotic {O}ption {P}ricing with {J}ump {D}iffusions and
other {L}\'evy processes.
\newblock {\em Journ. Comp. Fin.}, 22(3):13--47, 2018.
\bibitem{koponen}
I.~Koponen.
\newblock Analytic approach to the problem of convergence of truncated {L}\'evy
flights towards the {G}aussian stochastic process.
\newblock {\em Physics Review E}, 52:1197--1199, 1995.
\bibitem{kou}
S.G. Kou.
\newblock A jump-diffusion model for option pricing.
\newblock {\em Management Science}, 48(8):1086--1101, August 2002.
\bibitem{KW1}
S.G. Kou and H.~Wang.
\newblock First passage times of a jump diffusion process.
\newblock {\em Adv. Appl. Prob.}, 35(2):504--531, 2003.
\bibitem{KudrLev09}
O.~Kudryavtsev and S.Z. Levendorski\u{i}.
\newblock Fast and accurate pricing of barrier options under {L}\'evy
processes.
\newblock {\em Finance and Stochastics}, 13(4):531--562, 2009.
\bibitem{KudrLev11}
O.~Kudryavtsev and S.Z. Levendorski\u{i}.
\newblock Efficient pricing options with barrier and lookback features under
{L}\'evy processes.
\newblock Working paper, June 2011.
\newblock Available at SSRN: http://ssrn.com/abstract=1857943.
\bibitem{beta}
A.~Kuznetsov.
\newblock Wiener-{H}opf factorization and distribution of extrema for a family
of {L}\'evy processes.
\newblock {\em Ann.Appl.Prob.}, 20(5):1801--1830, 2010.
\bibitem{amer-put-levy-maphysto}
S.~Levendorski\u{i}.
\newblock Pricing of the {A}merican put under {L}\'evy processes.
\newblock Research Report MaPhySto, Aarhus, 2002.
\newblock Available at http://www.maphysto.dk/publications/MPS-RR/2002/44.pdf,
http://www.maphysto.dk/cgi-bin/gp.cgi?publ=441.
\bibitem{amer-put-levy}
S.~Levendorski\u{i}.
\newblock Pricing of the {A}merican put under {L}\'evy processes.
\newblock {\em International Journal of Theoretical and Applied Finance},
7(3):303--335, May 2004.
\bibitem{asymp-sens}
S.~Levendorski\u{i}.
\newblock Convergence of {C}arr's {R}andomization {A}pproximation {N}ear
{B}arrier.
\newblock {\em SIAM FM}, 2(1):79--111, 2011.
\bibitem{paired}
S.~Levendorski\u{i}.
\newblock Method of paired contours and pricing barrier options and {C}{D}{S}
of long maturities.
\newblock {\em International Journal of Theoretical and Applied Finance},
17(5):1--58, 2014.
\newblock 1450033 (58 pages).
\bibitem{AsianGamma}
S.~Levendorski\u{i}.
\newblock {D}ouble {S}piral method, {G}amma {T}ransform and pricing
{A}rithmetic {A}sian options.
\newblock Working paper, August 2016.
\newblock Available at SSRN: http://ssrn.com/abstract=2827138.
\bibitem{UltraFast}
S.~Levendorski\u{i}.
\newblock Ultra-{F}ast {P}ricing {B}arrier {O}ptions and {C}{D}{S}s.
\newblock {\em International Journal of Theoretical and Applied Finance},
20(5), 2017.
\newblock 1750033 (27 pages).
\bibitem{AsianGammaSIAMFM}
S.~Levendorski\u{i}.
\newblock Pricing arithmetic {A}sian options under {L}\'evy models by backward
induction in the dual space.
\newblock {\em SIAM FM}, 9(1):1--27, 2018.
\newblock Preliminary version available at SSRN:
http://papers.ssrn.com/abstract=2827138.
\bibitem{early-exercise}
S.Z. Levendorski\u{i}.
\newblock Early exercise boundary and option pricing in {L}\'evy driven models.
\newblock {\em Quantitative Finance}, 4(5):525--547, October 2004.
\bibitem{levendorskii-xie-Asian}
S.Z. Levendorski\u{i} and J.~Xie.
\newblock Pricing of {D}iscretely {S}ampled {A}sian {O}ptions {U}nder {L}\'evy
{P}rocesses.
\newblock Working paper, June 2012.
\newblock Available at SSRN: http://papers.ssrn.com/abstract=2088214.
\bibitem{LiLinetsky2015}
L.~Li and V.~Linetsky.
\newblock Discretely monitored first passage problems and barrier options: an
eigenfunction expansion approach.
\newblock {\em Finance and Stochastics}, 19(3):941--977, 2015.
\bibitem{Linetsky08}
V.~Linetsky.
\newblock Spectral methods in derivatives pricing.
\newblock In J.R. Birge and V.~Linetsky, editors, {\em Handbooks in OR \& MS,
Vol. 15}, pages 223--300. Elsevier, New York, 2008.
\bibitem{lipton-risk}
A.~Lipton.
\newblock Assets with jumps.
\newblock {\em Risk}, pages 149--153, September 2002.
\bibitem{lipton-columbia}
A.~Lipton.
\newblock Path-dependent options on assets with jumps.
\newblock 5$^{\textrm{th}}$ Columbia-Jaffe Conference, April 2002.
\newblock Available at http://www.math.columbia.edu/~lrb/columbia2002.pdf.
\bibitem{LiptonSelection}
A.~Lipton.
\newblock {\em Financial Engineering. Selected Works of Alexander Lipton}.
\newblock World Scientific, Singapore, 2018.
\bibitem{lipton-sepp}
A.~Lipton and A.~Sepp.
\newblock Credit value adjustment for credit default swaps via the structural
default model.
\newblock {\em Journal of Credit Risk}, 5(2):123--146, Summer 2009.
\bibitem{sato}
K.~Sato.
\newblock {\em L\'evy processes and infinitely divisible distributions},
volume~68 of {\em Cambridge Stud. Adv. Math.}
\newblock Cambridge University Press, Cambridge, 1999.
\bibitem{AbValko04b}
P.P. Valko and J.~Abate.
\newblock Comparison of sequence accelerators for the {G}aver {M}ethod of
{N}umerical {L}aplace {T}ransform inversion.
\newblock {\em Computers and Mathematics with Applications}, 48:629--636, 2004.
\bibitem{whystup_book_2e}
Uwe Wystup.
\newblock {\em {F}{X} {O}ptions and {S}tructured {P}roducts ({T}he {W}iley
{F}inance {S}eries)}.
\newblock John Wiley \& Sons, Chichester, West Sussex, UK, 2017.
\newblock 2d edition.
\end{thebibliography}