Operators and Boundary Problems in Finance, Economics and Insurance: Peculiarities, Efficient Methods and Outstanding Problems

Abstract

The price V of a contingent claim in finance, insurance and economics is defined as an expectation of a stochastic expression. If the underlying uncertainty is modeled as a strong Markov process X, the Feynman–Kac theorem suggests that V is the unique solution of a boundary problem for a parabolic equation. In the case of PDO with constant symbols, simple probabilistic tools explained in this paper can be used to explicitly calculate expectations under very weak conditions on the process and study the regularity of the solution. Assuming that the Feynman–Kac theorem holds, and a more general boundary problem can be localized, the local results can be used to study the existence and regularity of solutions, and derive efficient numerical methods. In the paper, difficulties for the realization of this program are analyzed, several outstanding problems are listed, and several closely efficient methods are outlined.

1. Introduction

1.1. Expectations and Boundary Problems

The prices (values) of contingent claims in finance, insurance and economics are defined as expectations of certain stochastic expressions. In many cases, the underlying uncertainty is modeled as a strong Markov process X (with killing) on ${\mathbb{R}}^n$ or its subset $\mathcal{D}$, with the infinitesimal generator L. Let $V(t,x)$ denote the price of an option or other contingent claim in the market at time t and $X_t=x$. The Feynman–Kac theorem suggests that V is the unique solution of a boundary problem for the parabolic equation:

$$({\partial }_t+L)V(t,x)+g(t,x)=0.$$

(1)

In the case of diffusion models, L is a differential operator of order 2; in jump-diffusion models, L is an integro-differential operator, and hence, a pseudo-differential operator (PDO). We use the representation of L as a PDO because this facilitates the study of the regularity of solutions, and naturally leads to the construction of efficient numerical methods for option pricing. The function g represents the stream of payoffs ${\left\{g(t,X_t)\right\}}_{t\ge 0}$ that the owner of the contingent claim is entitled to. The value function V satisfies (1) in the open region U in the time-state space, where the derivative security remains alive. As the process $(t,X_t)$ leaves U, the owner of the contingent claim is entitled to an instantaneous payoff $G({\tau }_{U^c},X_{{\tau }_{U^c}})$, where ${\tau }_{U^c}$ is the hitting time of $U^c=:\mathcal{D}\backslash U$ by $\left\{(t,X_t)\right\}$.

The solution of the boundary problem can be used to calculate the expectation, and the properties of the expectation can be used to formulate appropriate boundary and regularity conditions. Unfortunately, this scheme is very difficult to realize in almost all cases, including many popular diffusion models, and the majority of the results in the literature are obtained assuming (without proof) that the boundary problem for the expectation has the unique solution. Furthermore, it is assumed (also without proof) that the expectation defines a sufficiently regular function so that Ito’s formula can be applied; then, the solution of the boundary problem is the expectation. The necessity of proving these crucial assumptions is brushed under the carpet. Even in many diffusion models, the proof of the Feynman–Kac theorem is lacking. See [1,2] for details; one of the main reason is a complicated degeneration of L at the boundary of the state space, which makes it difficult to apply standard tools for the analysis of degenerate elliptic operators. If L degenerates at the boundary of a half-space, then general tools [3,4,5] can be applied, with some modifications, but in the case of degeneration at two and more transversal hyperplanes, the study of the regularity of solutions is more difficult. Furthermore, even if a weak regularity of the solution of the boundary problem is established, the proof of the Feynman–Kac theorem remains quite non-trivial. In the case of jump-diffusion models, additional subtleties emerge, even in the simplest case of Lévy models (the infinitesimal generator is a PDO with the constant symbol). In the majority of empirical studies, when Lévy models are calibrated to the real data, L is of the form ${\partial }_t+\langle \mu ,{\partial }_x\rangle -a\left(D_x\right)$, where $a\left(D\right)$ is an elliptic PDO of order $\nu \in (0,1)$, with the symbol analytic in a tube domain. Hence, the operator of the boundary problem (1) becomes a parabolic operator of the standard form, with an elliptic stationary part, only after an appropriate change of variables in the $(t,x)$-space. This peculiarity leads to several non-standard properties of solutions, which we analyze in this paper. If one has in view applications of the general theory of fractional differential equations to finance, the class of operators of this kind deserves to be regarded as a model class rather than standard fractional differential operators ${\partial }_t-{|\Delta |}^{\nu }$, where ${|\Delta |}^{\nu }$ is the fractional Laplacian. (Note that ${|\Delta |}^{\nu }$, $\nu \ne 2$, cannot be used in the standard popular models in finance, economics and insurance.) An additional important feature of the infinitesimal generators of Lévy models used in finance is the existence of analytic continuation not only to tube domains but to conesas well. This observation is the basis of efficient numerical methods that evaluate integrals in pricing formulas, which we outline in this paper. These methods can be used in other situations, for instance, to evaluate special functions and stable distributions [6,7,8,9].

1.2. Black–Scholes Model and Diffusion Models

In the Black–Scholes model, the only source of uncertainty is the stock price $S_t=e^{X_t}$, where $X_t$ is the Brownian motion (BM) on $\mathbb{R}$, with drift. The infinitesimal generator L is the second-order elliptic differential operator. If the market with several stocks and/or several sources of uncertainty is considered, e.g., volatility and/or stochastic interest rate, then the underlying process X is of a more complicated nature.

Typically, the region U where (1) holds is of the form $\left\{\right(t,x\left)\right|t\in (0,T),x\in U\left(t\right)\}$, where $U\left(t\right)$ are open subsets of $\mathcal{D}$. If X is a diffusion, $X_{{\tau }_{U^c}}$ is at the boundary of U, and then V satisfies the terminal condition:

$$V(T,x)=G(T,x),x\in U\left(T\right),$$

(2)

and the boundary condition:

$$V(t,x)=G(t,x),x\in \partial U\left(t\right).$$

(3)

If U is unbounded, appropriate restrictions on the rate of growth of V at infinity are imposed. Boundary problems of this sort arise for the majority of contingent claims in finance and insurance, and for numerous value functions in economics, e.g., real options and stochastic games. In some important cases such as lookback or Asian options, the underlying process is not Markovian; hence, pricing is more involved. In this paper, we consider boundary problems of the form (1), (2), (3), including free boundary problems. In applications, the approximate pricing of contingent claims of a complicated structure is reduced to a sequence of embedded options, equivalently, to a sequence of boundary problems of the form (1), (2), (3). Hence, this study of the regularity of the solutions of the latter is important for studying other types of contingent claims as well. Efficient numerical procedures for the solutions of standard boundary problems can be (and are) used as building blocks to solve more complicated problems.

In the Black–Scholes model, a simple affine change of variables reduces (1) to the backward parabolic equation $V_t+V_{xx}-rV=0$; hence, the host of standard well-known methods can be used to price numerous contingent claims. Black and Scholes informally derived (1), constructing the perfect continuously changing hedging portfolio; later, Merton gave essentially equivalent proof constructing the perfect continuously changing replicating portfolio. Both proofs contain fundamental mathematical errors but can be made accurate for wide classes of diffusion models. Unfortunately, both proofs fail for jump-diffusion models—the statements and proofs of several popular imperfect substitutes for perfect hedging and replication are also incorrect (see [10]), and the reason for the failure is fundamental rather than technical.

1.3. The Case of Jump-Diffusions

In the Merton–Black–Scholes theory, the pricing equation is derived from the absence of frictions and no-arbitrage assumption. The no-arbitrage assumption means that it is impossible to construct a portfolio of securities traded in the market such that, at the expiration date, the value of the portfolio is non-negative with probability 1, and positive with positive probability. Leaving aside an important theoretical background from economics and finance, and the technical conditions necessary to make the statements mathematically accurate, the absence of frictions and the no-arbitrage assumption imply that the discounted prices of all securities traded in the market must be local martingales under a probability measure $\mathbb{Q}$ on the filtered probability space where the process X lives. The crucial point is that $\mathbb{Q}\ne \mathbb{P}$, where $\mathbb{P}$ is the historic or physical probability measure estimated using the time series for the prices of securities already traded on the market. The pricing measure $\mathbb{Q}$ may not assign non-zero probabilities to events of zero measure under $\mathbb{P}$ and vice versa; hence, $\mathbb{Q}$ must be equivalent to $\mathbb{P}$. This explains the name for $\mathbb{Q}$: an equivalent martingale measure (EMM; another name is risk-neutral measure). The discounting can be taken into account as the killing of the Markov process, and $L=L^{\mathbb{Q}}$ in (1) is the infinitesimal generator of the process with killing, under $\mathbb{Q}$. Then, $\left\{V(t,X_t)\right\}$, the discounted price process, is a local martingale if for any $t<{\tau }_{U^c}$ and stopping time $\tau \ge t$ s.t. $\tau \le {\tau }_{U^c}$, ${\mathbb{E}}^{\mathbb{Q}}\left[V(\tau ,X_{\tau })\right|X_t=x]=V(t,x)$.

The boundary problem (1), (2), (3) is used to calculate the prices of new securities. In (sufficiently regular) diffusion models, there exists a unique EMM; hence, one can calculate the price of any contingent claim. In financial markets and insurance, a great variety of contingent claims are constructed and sold, and the implicit assumption is that one can theoretically calculate the price of new securities.

The three conditions are as follows: (1) a perfect hedge is possible; (2) perfect replication is possible; (3) there exists a unique EMM are equivalent (naturally, under subtle technical conditions)—and markets satisfying these conditions are called complete markets. The Feynman–Kac theorem gives the representation of the price as the expectation of the payoff stream under $\mathbb{Q}$:

$$V(t,x)={\mathbb{E}}^{\mathbb{Q}}\left[{\int }_t^{{\tau }_{U^c}}g(s,X_s)ds+G({\tau }_{U^c},V_{{\tau }_{U^c}})\right].$$

(4)

If X is a jump process or jump-diffusion process, then the hedging/replicating argument is not applicable. Fortunately, according to the general economic theory, if the market does not admit arbitrage, there exists an EMM such that the price of each contingent claim is given by (4) (an accurate mathematical statement imposes additional subtle conditions on X); an EMM is not unique, and one may choose an EMM one believes in. In real financial markets and numerous empirical studies, the measures calibrated to prices of the underlying financial instruments and options are not only different; quite often, the measures are not equivalent. Several examples can be found in the well-known empirical study in [11]: prices of certain stocks are calibrated to processes of infinite variation, whereas the prices of options—to processes of finite variation, which contradicts the well-known conditions on equivalent probability measures in [12] (the authors of [11] did not notice this discrepancy of their calibration results). One of the important reasons for discrepancies of this kind is the inaccuracy of popular numerical methods used for pricing and calibration. Inaccurate methods are constructed and used because crucial qualitative properties of jump-diffusion models and prices $V(t,x)$ in these models are not taken into account. Examples of errors of popular methods, including the analysis of implications for risk management, can be found in [2,13,14,15,16,17,18,19,20]. The aims of this paper are to list and explain the irregularity of $V(t,x)$ (and the free boundary, in the case of options of American type) in several standard situations, and explain how to design efficient numerical methods that work well in wide classes of jump-diffusion models.

In incomplete markets and the jump-diffusion models which are models of incomplete markets, one cannot use the hedging or replication argument to derive the equation for contingent claims. However, one can start by choosing an EMM and defining the price by (4). In many cases of interest, the representation (4) and the Fourier/Laplace transform technique suffice to express the price in the form of an integral, and an efficient numerical procedure for the evaluation of the designed integral. In some popular classes of models, the derivation of pricing formulas is based on (1), although the rigorous justification is lacking. Informally, one can use Dynkin’s formula:

$$V(t,x)={\mathbb{E}}^{\mathbb{Q}}\left[{\int }_t^{{\tau }_{U^c}}(-{\partial }_s-L)V(s,X_s)ds+G({\tau }_{U^c},V_{{\tau }_{U^c}})\right]$$

(5)

to conclude that (1) must hold. Thus, (1) is the backward Kolmogorov equation. In [21,22], the formal derivation of (1) is given for Lévy processes satisfying the (ACP)-condition (absolute continuity of potential measures) in the infinite time horizon case (stationary problem) and the (ACT)-condition (absolute continuity of transition measures) in the non-stationary case—in ([12], Section 41), one can find equivalent conditions for the (ACP)- and (ACT)-conditions; the equation is understood in the sense of generalized functions. A similar proof can be given for wider classes of strong Markov processes satisfying the same conditions. However, the author is unaware of a published proof in the general setting. See [1,2] for a review of partial results, and [23,24,25] for proofs in several special cases.

The next important complication in the case of jump-diffusion processes is the form of the boundary condition which becomes non-local:

$$V(t,x)=G(t,x),x\in U{\left(t\right)}^c.$$

(6)

Note that the standard boundary problems for fractional differential equations are local although the form of the conditions is non-standard. The non-standard boundary conditions are formally invented in order for the Cauchy problem to be well defined. In the literature, one can find discussions about a proper choice of the boundary conditions. On the contrary, the non-local boundary condition (6) is a part of the definition of the value function, and cannot be replaced for convenience with a local condition.

One can use the boundary problem (1), (2), (6) as follows:

(1)

To prove the existence and uniqueness of the solution in the class of sufficiently regular functions; sufficiently regular means that Dynkin’s formula (5) is valid;

(2)

Applying (5) to conclude that V equals the RHS of (4).

Unfortunately, this scheme is very difficult to realize in many models, including many popular diffusion models, and the majority of the results in the literature are obtained assuming without proof that the boundary problem has the unique solution, and the solution is sufficiently regular so that Ito’s formula (or, more generally, Dynkin’s formula) can be applied; then, the solution satisfies (4) and (5), hence, the problem of calculation of the expectation (4) is solved. The necessity of proving this crucial assumption is brushed under the carpet.

The main difficulty for the proof stems from the irregularity of the value function V, which makes general theorems of Feynman–Kac type available in the literature not applicable; the irregularity of V is the artifact of the properties of the infinitesimal generators. However, for several classes of contingent claims and types of models, V given by (4) can be explicitly calculated using the Fourier/Laplace technique. In the case of Lévy processes (in terms of PDO, the case of operators with constant symbols), and problems of several basic types (boundary problems with flat boundaries), either no or very weak conditions on the process are needed. The result is an integral depending on $(t,x)$ as a parameter, in one or more dimensions. An explicit analytical expression can be used to study the regularity of solutions; in a number of important case, the integral can be efficiently and quickly calculated. Since the study of the regularity of solutions of boundary problems can be localized (see [26]), the author hopes that the results for operators with constant symbols and boundary problems with flat boundaries can be used to prove the correspondence between the stochastic expressions and boundary problems with curved boundaries, and operators with state- and time-dependent symbols.

1.4. Structure of This Paper

Lévy models are considered in Section 2, and the pricing of European options using the Fourier transform technique (solution of the Cauchy problem) and efficient numerical realizations are in Section 3. To price barrier options (boundary problems for parabolic operators), the Laplace transform, Wiener–Hopf factorization and maturity randomization (method of lines) are needed (Section 4). Explicit formulas involve multi-dimensional integrals, and efficient numerical calculations based on contour deformations become more involved. In Section 5, we analyze the peculiarities of the early exercise boundary and the prices of American options (solutions of free boundary problems), and present a general stable scheme based on the method of lines. The proof of convergence is based on the probabilistic interpretation of the operator form of the Wiener–Hopf factorization. In Section 6, we explain how the methods of Section 4 and Section 5 are modified to price options in regime-switching models (solve systems of boundary problems), and how to approximate more complicated stochastic volatility models and models with stochastic interest rates with regime-switching models. In Section 7, we outline the structure of several exactly solvable models with non-constant symbols and list several outstanding problems for these classes of models, which seem to be non-trivial and interesting from the point of view of the theory of boundary problems for PDE and PDO. Finally, in Section 8, we summarize the results presented in this paper, review several other groups of methods and, wherever possible, outline the relative advantages of different methods.

2. Lévy Models or PDO with Constant Symbols

2.1. Lévy Processes

Let X be the Lévy process on the filtered probability space $(\Omega ;\mathcal{F};{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0};\mathbb{Q})$. $\mathbb{E}={\mathbb{E}}^{\mathbb{Q}}$ denotes the expectation operator under $\mathbb{Q}$. We use the definition in [21,27] of the characteristic exponent $\psi \left(\xi \right)={\psi }^{\mathbb{Q}}\left(\xi \right)$ of a Lévy process X on ${\mathbb{R}}^n$, under $\mathbb{Q}$, which is marginally different from the definition in [12]. Namely, $\psi$ is definable from:

$$\mathbb{E}[e^{i\langle \xi ,X_t\rangle }]=e^{-t\psi \left(\xi \right)},\xi \in {\mathbb{R}}^n,t\ge 0,$$

(7)

where $\langle a,b\rangle ={\sum }_{j=1}^n{a}_j{b}_j$. This definition of the characteristic exponent implies that the infinitesimal generator $L_X$ of X is the pseudo-differential operator (pdo) $-\psi \left(D\right)$; equivalently, $L_X$ acts of oscillating exponents are as follows: $L_X{e}^{i\langle x,\xi \rangle }=-\psi \left(\xi \right)e^{ix\xi }$.

Denote $D=\left\{x\right|\left|x\right|\le 1\}$. The following theorem (Lévy–Khintchine theorem) can be found in many monographs, e.g., ([12], Thm. 8.1).

Theorem 1. 

(i) Let X be a Lévy process on ${\mathbb{R}}^n$. Then, its characteristic exponent admits the representation:

$$\psi \left(\xi \right)=\frac{1}{2}\langle A\xi ,\xi \rangle -i\langle \gamma ,\xi \rangle +{\int }_{{\mathbb{R}}^n}(1-e^{i\langle x,\xi \rangle }+i\langle x,\xi \rangle {\mathbf{1}}_D\left(x\right))F\left(dx\right),$$

(8)

where A is a symmetric non-negative-definite $n\times n$ matrix, $\gamma \in {\mathbb{R}}^n$, and $F\left(dx\right)$ is a measure on ${\mathbb{R}}^n$ satisfying:

$$F\left(\left\{0\right\}\right)=0,{\int }_{{\mathbb{R}}^n}{\left(\right|x|}^2\wedge 1)F\left(dx\right)<\infty .$$

(9)

(ii) The representation (8) is unique;

(iii) Conversely, if A is a symmetric non-negative-definite $n\times n$ matrix, $\gamma \in {\mathbb{R}}^n$, and F is a measure on ${\mathbb{R}}^n$ satisfying (9), then there exists a Lévy process X with the characteristic exponent (7).

The triple $(A,F,\gamma )$ is called the generating triplet of X. The A and F are called the Gaussian covariance matrix and Lévy measure of X. When $F=0$, X is Gaussian, and if $A=0$, X is called purely non-Gaussian.

Essentially, the term $-i\langle x,\xi \rangle {\mathbf{1}}_D\left(x\right)$ in (8) is needed to ensure the convergence of the integral, and hence different functions can be (and are) used instead of $c\left(x\right):={\mathbf{1}}_D\left(x\right)$, for instance, $c\left(x\right)=1/(1+|x{|}^2)$; the A and F are independent of the choice of c. If F satisfies the condition:

$$F\left(\left\{0\right\}\right)=0,{\int }_{{\mathbb{R}}^n}\left(\right|x|\wedge 1)F\left(dx\right)<\infty ,$$

(10)

which is stronger than (9), then (8) can be simplified:

$$\psi \left(\xi \right)=\frac{1}{2}\langle A\xi ,\xi \rangle -i\langle {\gamma }_0,\xi \rangle +{\int }_{{\mathbb{R}}^n}(1-e^{i\langle x,\xi \rangle })F\left(dx\right),$$

(11)

where ${\gamma }_0=\gamma -{\int }_{{\mathbb{R}}^n}x{\mathbf{1}}_D\left(x\right)F\left(dx\right).$

If the sample paths of a Lévy process have bounded variation on every compact time interval a. s., we say that the Lévy process has bounded variation. A Lévy process has bounded variation if and only if $A=0$ and (10) holds (see, e.g., [28], p. 15).

2.2. Examples of Lévy Processes on $\mathbb{R}$

Example 1.

The Lévy density of a (pure jump) one-dimensional stable Lévy process X of index $\alpha \in (0,2)$ is of the form

$$\begin{array}{ccc}\hfill F\left(dy\right)& =& {\left|y\right|}^{-\alpha -1}(c_{+}{\mathbf{1}}_{(0,+\infty )}\left(y\right)+c_{-}{\mathbf{1}}_{(-\infty ,0)}\left(y\right))dy,\hfill \end{array}$$

(12)

where $c_{\pm }\ge 0$ and $c_{+}+c_{-}>0$. If $\alpha \ne 1$, then, substituting (12) into the Lévy–Khintchine formula with ${\sigma }^2=0$, one easily derives ${\psi }_{st}\left(\xi \right)=-i\mu \xi +{\psi }_{st}^0(\alpha ,C_{+},\xi )$, where μ can be expressed in terms of α, $c_{\pm }$ and b:

$${\psi }_{st}^0(\alpha ,C_{+},\xi )=C_{+}{\xi }^{\alpha }{\mathbf{1}}_{(0,+\infty )}\left(\xi \right)+C_{-}{(-\xi )}^{\alpha }{\mathbf{1}}_{(-\infty ,0)}\left(\xi \right),$$

(13)

$C_{-}=\overline{C_{+}}$, and $C_{+}=C_{+}(\alpha ,c_{+},c_{-})=-c_{+}\Gamma (-\alpha )e^{-i\pi \alpha /2}-c_{-}\Gamma (-\alpha )e^{i\pi \alpha /2}.$ See [9]. For somewhat different (naturally, equivalent) formulas, see [29] and ([12], Thm.14.15). In the case $\alpha =1$, the formula for ${\psi }^0$ is different (see [9]):

$${\psi }^0\left(\xi \right)=\sigma \left|\xi \right|(1+i(2\beta /\pi )\mathrm{sign}\xi ln\left|\xi \right|),$$

(14)

where $\sigma =(c_{+}+c_{-})\pi /2$, $\beta =(c_{+}-c_{-})/(c_{+}+c_{-})$. This is a version of Zolotarev’s parametrizations [30] for stable processes of index 1. See [9,12] for further references. Note that the symbol of the infinitesimal generator L is non-smooth at zero.

The tails of the probability distributions of stable Lévy processes exhibit slow polynomial decay. Hence, the second moment of the distribution of $X_t$ is infinite. In real financial markets, the distributions have finite second moments, which makes stable processes unsuitable. The simplest way to ensure that the second moment is finite is to consider processes with exponentially decaying tails; the symbols of L are analytic in a strip. These equivalent properties were used as the basis for the definition of a general class of Regular Lévy processes of exponential type (RLPE) [21,22]. Below, we list several popular sub-classes of RLPE. The reader can easily observe that, in all these examples, the characteristic exponent $\psi \left(\xi \right)=-i\mu \xi +{\psi }^0\left(\xi \right)$ admits analytic continuation to the union of a strip and cone $\mathcal{C}$ around $\mathbb{R}$, and $\mathrm{Re}{\psi }^0\left(\xi \right)\to +\infty$ as $\xi \to \infty$ in a sub-cone ${\mathcal{C}}_{+}\subset \mathcal{C}$. In [31], the definition of the general class of processes enjoying these properties was given. The suggested name SINH-regular processes reflects the fact that the integrals in pricing formulas can be efficiently calculated using changes of variables of the form:

$$\xi ={\chi }_{{\omega }_1,b,\omega }\left(y\right)=i{\omega }_1+bsinh(i\omega +y)$$

(15)

(sinh-acceleration). Lévy processes used in quantitative finance are SINH-regular processes.

Example 2.

KoBoL processes. Modifying the Lévy density (12):

$$\begin{array}{ccc}\hfill F\left(dy\right)& =& (c_{+}{y}^{-{\nu }_{+}-1}{e}^{{\lambda }_{-}y}{\mathbf{1}}_{(0,+\infty )}\left(y\right)+c_{-}{(-y)}^{-{\nu }_{-}-1}{e}^{{\lambda }_{+}y}{\mathbf{1}}_{(-\infty ,0)}\left(y\right))dy,\hfill \end{array}$$

(16)

where ${\nu }_{\pm }\in (0,2)$, $c_{\pm }\ge 0$, $c_{+}+c_{-}>0$ and ${\lambda }_{-}<0<{\lambda }_{+}$, we obtain a class of Lévy processes with exponentially decaying tails. In the case ${\nu }_{-}\ne 1,{\nu }_{+}\ne 1$, the characteristic exponent of the KoBoL process is of the form $\psi \left(\xi \right)=-i\mu \xi +{\psi }^0\left(\xi \right)$, where:

$${\psi }^0\left(\xi \right)=c_{-}\Gamma (-{\nu }_{-})[{\lambda }_{+}^{{\nu }_{-}}-{({\lambda }_{+}+i\xi )}^{{\nu }_{-}}]+c_{+}\Gamma (-{\nu }_{+})[{(-{\lambda }_{-})}^{{\nu }_{+}}-{(-{\lambda }_{-}-i\xi )}^{{\nu }_{+}}];$$

(17)

if either ${\nu }_{+}=1$ or ${\nu }_{-}=1$, then the formula for ${\psi }^0$ is different. See [21,27]. In particular, if ${\nu }_{+}={\nu }_{-}=1$, then:

$$\begin{array}{ccc}\hfill {\psi }^0\left(\xi \right)& =& c_{+}((-{\lambda }_{-})ln(-{\lambda }_{-})-(-{\lambda }_{-}-i\xi )ln(-{\lambda }_{-}-i\xi ))\hfill \\ \hfill & & +c_{-}({\lambda }_{+}ln{\lambda }_{+}-({\lambda }_{+}+i\xi )ln({\lambda }_{+}+i\xi )).\hfill \end{array}$$

(18)

In the symmetric case ${\nu }_{-}={\nu }_{+}\in (0,2)\backslash \left\{1\right\},-{\lambda }_{-}={\lambda }_{+},c_{+}=c_{-}$, the class of processes with the Lévy density (16) was introduced by Koponen [32], who derived a rather inconvenient formula for the characteristic exponent (different from (17)) and suggested a somewhat misleading name truncated Lévy processes. The generalization (16) was introduced in [27], where the characteristic exponent was calculated and several option pricing problems were solved; the name Koponen’s family of truncated Lévy processes was used. Starting with [21], the name KoBoL processes is used. In [11], a subclass of KoBoL with ${\nu }_{-}={\nu }_{+}\ne 1$ and $c_{+}=c_{-}$ was called CGMY model and labels for the parameters of the KoBoL model were changed.

Rosinski [33] suggested the name exponentially tilted stable Lévy processes and gave a general definition of a class which is a subclass of the class RLPE.

Example 3.

Normal inverse Gaussian (NIG) processes, and the generalization: normal tempered stable (NTS) processes are constructed in [34,35], respectively. The characteristic exponent is given by

$${\psi }^0\left(\xi \right)=\delta [{({\alpha }^2+{(\xi +i\beta )}^2)}^{\nu /2}-{({\alpha }^2-{\beta }^2)}^{\nu /2}],$$

(19)

where $\nu \in (0,2)$, $\delta >0$, $\left|\beta \right|<\alpha$; NIG obtains with $\nu =1$.

Example 4.

Variance Gamma processes (VGPs) were introduced to finance in [36]. The characteristic exponent can be written in the form:

$${\psi }^0\left(\xi \right)=c[ln({\alpha }^2-{(\beta +i\xi )}^2)-ln({\alpha }^2-{\beta }^2)],$$

(20)

where $\alpha >\left|\beta \right|\ge 0$, $c>0$.

Example 5.

In the Merton model [37], the characteristic exponent is given by

$${\psi }^0\left(\xi \right)=\frac{{\sigma }^2}{2}{\xi }^2+\lambda ·\left(1-e^{im\xi -\frac{s^2}{2}{\xi }^2}\right),$$

(21)

where $\sigma ,s,\lambda >0$ and $\mu ,m\in \mathbb{R}$.

Example 6.

The hyper-exponential jump-diffusion processes (HEJD model) were introduced in [38,39] and studied in detail in [38,40]). The characteristic exponent is of the form:

$${\psi }^0\left(\xi \right)=\frac{{\sigma }^2}{2}{\xi }^2+{\lambda }^{+}·\sum _{j=1}^{n^{+}}\frac{i{p}_j^{+}\xi }{i\xi -{\alpha }_j^{+}}+{\lambda }^{-}·\sum _{k=1}^{n^{-}}\frac{i{p}_k^{-}\xi }{i\xi +{\alpha }_k^{-}},$$

(22)

where $n^{\pm }$ are positive integers and ${\alpha }_j^{\pm },{\lambda }^{\pm },p_j^{\pm }>0$ satisfy ${\sum }_{j=1}^{n^{\pm }}{p}_j^{\pm }=1$. A double-exponential jump diffusion model introduced to finance in [41] (and well known for decades) can be obtained as a special case of hyper-exponential jump-diffusion models by taking $n^{+}=n^{-}=1$.

Example 7.

Other classes of Lévy processes with rational characteristic exponents and non-trivial BM components [21,38,40,42,43].

Example 8.

The characteristic exponents of the processes of the β-class constructed in [44] are of the form:

$$\begin{array}{ccc}\hfill {\psi }^0\left(\xi \right)& =& \frac{{\sigma }^2}{2}{\xi }^2+\frac{c_1}{{\beta }_1}\left\{B({\alpha }_1,1-{\gamma }_1)-B({\alpha }_1-\frac{i\xi }{{\beta }_1},1-{\gamma }_1)\right\}\hfill \\ \hfill & & +\frac{c_2}{{\beta }_2}\left\{B({\alpha }_2,1-{\gamma }_2)-B({\alpha }_2+\frac{i\xi }{{\beta }_2},1-{\gamma }_2)\right\},\hfill \end{array}$$

(23)

where $c_j\ge 0,{\alpha }_j,{\beta }_j>0$ and ${\gamma }_j\in (0,3)\backslash \{1,2\}$, and $B(x,y)$ is the Beta-function. Evidently, all poles of ${\psi }^0$ are on $i{\mathbb{R}}_{+}\backslash 0$.

Example 9.

The Lévy density and characteristic exponent of the meromorphic Lévy processes introduced in [45] are defined by almost the same formulas as in the HEJD model. The difference is that the sums are infinite. A natural condition ${\alpha }_j^{\pm }\to +\infty$ as $j\to +\infty$ is imposed, and the requirement that an infinite sum defines a Lévy density is equivalent to ${\sum }_{j\ge 0}{p}_j^{\pm }{\left({\alpha }_j^{\pm }\right)}^{-2}<+\infty .$ The formula for ${\psi }^0$ is (22) with the infinite numbers of terms; the poles of $\psi \left(\xi \right)$ are on $i\mathbb{R}\backslash 0$.

Remark 1.

If one-factor Lévy models are calibrated to prices of stocks and indices in financial markets, then, in the majority of cases, KoBoL processes of order $\nu \in (0,1)$, without the diffusion component (or very small diffusion component) give the best fit. Typically, processes with a jump part of finite activity such as the Merton model, HEJD model and other models with rational characteristic exponents do not calibrate well to the real data but one may try to use simple models to approximate involved ones. There are publications where HEJD are used to approximate KoBoL dynamics. However, it is evident that the approximations of operators of order $\nu \in (0,1)$ by operators of order 2 cannot work well near the boundary. The reader can find examples in [46] which illustrate this point. A very flexible multi-parameter β-family and meromorphic processes can be used to approximate KoBoL processes, but these approximations are much less efficient than direct calculations in the KoBoL model.

For applications to qualitative problems in economics, processes with rational characteristic exponents are more suitable due to the triviality of the Wiener–Hopf factorization. See, e.g., [47,48,49,50,51,52].

2.3. Examples of Lévy Processes on ${\mathbb{R}}^n$

Example 10.

In the notation used in this paper, the characterization of non-trivial stable Lévy processes on ${\mathbb{R}}^n$ of index $\alpha \in (0,2)$ ([12], Thm. 14.10) is as follows. There exist $\mu \in {\mathbb{R}}^n$ and a finite non-zero measure $G_{st}\left(d\varphi \right)$ on the unit sphere $S_{n-1}:=\{\xi \in {\mathbb{R}}^n\left|\right|\xi |=1\}$ such that if $\alpha \ne 1$, then:

$$\begin{array}{ccc}\hfill \psi \left(\xi \right)& =& -i\langle \mu ,\xi \rangle +{\int }_{S_{n-1}}{\left|\langle \xi ,\varphi \rangle \right|}^{\alpha }\left(1+itan\frac{\alpha \pi }{2}\mathrm{sign}\langle \xi ,\varphi \rangle \right)G_{st}\left(d\varphi \right)\hfill \end{array}$$

(24)

is the characteristic exponent of a stable Lévy process of index α, and if $\alpha =1$, then:

$$\begin{array}{ccc}\hfill \psi \left(\xi \right)& =& -i\langle \mu ,\xi \rangle +{\int }_{S_{n-1}}\left|\langle \xi ,\varphi \rangle \right|\left(1-i\frac{2}{\pi }tan\frac{\alpha \pi }{2}\langle \xi ,\varphi \rangle ln\left|\langle \xi ,\varphi \rangle \right|\right)G_{st}\left(d\varphi \right)\hfill \end{array}$$

(25)

Conversely, for any $\mu \in {\mathbb{R}}^n$ and a finite non-zero measure $G_{st}\left(d\varphi \right)$ on $S_{n-1}$, (24) and (25) define the characteristic exponents of a stable Lévy process of index $\alpha \ne 1$ and $\alpha =1$, respectively.

A straightforward multi-dimensional analog of RLPE class was introduced in ([21], Sect. 9.1.4). The characteristic exponent of an RLPE admits analytic continuation to a tube domain U containing ${\mathbb{R}}^n$, and stabilizes to a positively homogeneous function as $\xi \to \infty$ remaining in U. As in the 1D-case, we impose conditions on ${\psi }^0$ in the representation:

$$\psi \left(\xi \right)=-i\langle \mu ,\xi \rangle +{\psi }^0\left(\xi \right),$$

(26)

where $\mu \in {\mathbb{R}}^n$.

Example 11.

Consider the multi-factor KoBoL family of pure jump Lévy processes constructed in ([21], Section 9.1.1). Let $\alpha \in (0,2)$, and let $G\left(dx\right)$ be a finite non-zero measure and λ a positive continuous function on the unit sphere $S_{n-1}$. Then:

$$F\left(dx\right)={\rho }^{-\alpha -1}exp(-\lambda \left(\varphi \right)\rho )d\rho G\left(d\varphi \right)$$

(27)

is a Lévy measure. If $\alpha \in (0,2),\alpha \ne 1,$ the characteristic exponent is:

$${\psi }^0\left(\xi \right)=\Gamma (-\alpha ){\int }_{S_{n-1}}[\lambda {\left(\varphi \right)}^{\alpha }-{(\lambda \left(\varphi \right)-i\langle \xi ,\varphi \rangle )}^{\alpha }]G\left(d\varphi \right)$$

(28)

and if $\alpha =1$:

$${\psi }^0\left(\xi \right)={\int }_{S_{n-1}}[\lambda \left(\varphi \right)ln\lambda \left(\varphi \right)-(\lambda \left(\varphi \right)-i\langle \xi ,\varphi \rangle )ln(\lambda \left(\varphi \right)-i\langle \xi ,\varphi \rangle )]G\left(d\varphi \right).$$

(29)

Passing to the limit $\lambda \downarrow 0$ in (28) and (29), one can easily derive alternative representations of the characteristic exponents of stable Lévy processes.

Example 12.

The class of multi-factor normal tempered stable (NTS) Lévy processes constructed in ([21], Section 9.1.2) can be defined by

$${\psi }^0\left(\xi \right)=\delta \left\{[{({\alpha }^2+\langle A(\xi -i\beta ),(\xi -i\beta )\rangle ]}^{\nu /2}-{[{\alpha }^2-\langle A\beta ,\beta \rangle ]}^{\nu /2}\right\},$$

(30)

where $\delta ,\alpha >0$, A is a positive-definite matrix, $\beta \in {\mathbb{R}}^n$, ${\alpha }^2-\langle A\beta ,\beta \rangle >0$.

The following is a straightforward definition of the class of SINH-regular processes in [31], in the simplest form.

Definition 1.

We say that X is a SINH-regular process of order $(\nu ,{\nu }^{\prime })$ if the characteristic exponent is of the form ${\psi }^0$ admits analytic continuation to a domain of the form $\mathcal{U}=iD+\mathcal{C}$, where D is an open set containing $\left\{0\right\}$, and $\mathcal{C}$ is an open cone $\mathcal{C}$ around ${\mathbb{R}}^n$ s.t.:

$$|{\psi }^0{\left(\xi \right)|\le C(1+\left|\xi \right|)}^{\nu },\xi \in \mathcal{U},$$

(31)

where C is independent of ξ, and there exists ${\mathcal{U}}_{+}\subset \mathcal{U}$, of the same form: ${\mathcal{U}}_{+}=i{D}_{+}+{\mathcal{C}}_{+}$, s.t.:

$$\mathrm{Re}{\psi }^0\left(\xi \right)\ge c{\left|\xi \right|}^{{\nu }^{\prime }}-C,\xi \in {\mathcal{U}}_{+},$$

(32)

where $C,c>0$ are independent of ξ.

2.4. General Remarks

• Stable Lévy processes can be characterized as the limiting case of SINH-regular processes when the tube domain $iD+{\mathbb{R}}^n$ of analyticity shrinks to $\left\{0\right\}$; the conditions are valid in $\mathcal{C}$ only. The calculation of expectations in stable Lévy models can be efficiently performed by either modifying the sinh-acceleration technique or approximating stable Lévy processes with SINH-regular ones [8,9].
• The definition in [31] allows for ${\mathcal{U}}_{+}$ to be adjacent to $\mathbb{R}$, and Definition 1 can be generalized in a similar fashion.
• In [31], the bounds (31) and (32) are formulated for $\psi$ rather than ${\psi }^0$. One can easily derive the bounds for $\psi$ using the bounds for ${\psi }^0$.
• It is easy to see that NTS processes are SINH-regular but a multi-factor KoBoL X is SINH-regular only if X is a mixture of independent KoBoL in 1D.
• VGP and their multi-factor generalizations are SINH-regular if we replace the weight ${\left|\xi \right|}^{\nu }$ with $ln(2+|\xi \left|\right)$.
• As in [31], more general weight functions can be used. It can be shown that, in the case of pure jump processes, ${\psi }^0\left(\xi \right)={o\left(\right|\xi |}^2)$ as $(\mathcal{C}\ni )\xi \to \infty$; hence, one can use the upper bound with $\nu =2$ in all cases.
• If $\nu ={\nu }^{\prime }$, then ${\psi }^0$ is an elliptic symbol. If ${\nu }^{\prime }<\nu$, then, typically, $\nu -{\nu }^{\prime }<1$, hence, ${\psi }^0$ is hypo-elliptic.
• If in the representation (26), $\mu \ne 0$, and $\nu <1$, then the principal symbol of $\psi \left(\xi \right)$ is $-i\langle \mu ,\xi \rangle$, which leads to the irregularity of the Wiener–Hopf factors and the solutions of the boundary problems.
• In the case of Cauchy problems in the whole space, the drift can be eliminated by the change of variables $x=x^{\prime }-t\mu$, and the same change of variables is implicit in efficient numerical methods for the Fourier inversion methods based on the conformal deformation of lines and hyperplanes of integration. Formally, the same change of variables can be applied in the case of more general boundary problems but the change makes a flat boundary (typical in pricing problems for standard barrier options) non-flat. When the boundary is flat, explicit pricing formulas can be derived and the study of the (ir)regularity of the solutions simplified.

3. Pricing European Options in Lévy Models and Cauchy Problems in ${\mathbb{R}}^n$ for Operators with Constant Symbols

3.1. Exact Formulas

Let $r\ge 0$ be the constant riskless rate, X a Lévy process in ${\mathbb{R}}^n$ under an EMM $\mathbb{Q}$ chosen for pricing, and $G\left(X_T\right)$ the payoff of the contingent claim at maturity date T. The price $V(t,x)$ of the claim at time t and $X_t=x$ is given by

$$V(t,x)={\mathbb{E}}^{\mathbb{Q}}\left[e^{-r(T-t)}G\left(X_T\right)|X_t=x\right].$$

(33)

Assume that the characteristic exponent $\psi$ of X is analytic in a tube domain $iD+{\mathbb{R}}^n$, where D is an open set containing $\left\{0\right\}$, and the Fourier transform:

$$\widehat{G}\left(\xi \right)={\left(2\pi \right)}^{-n}{\int }_{{\mathbb{R}}^n}{e}^{-i\langle x,\xi \rangle }G\left(x\right)dx$$

(34)

is analytic in a tube domain $i{D}^{\prime }+{\mathbb{R}}^n$, where $D^{\prime }\cap D$ is an non-empty open set. Let the product $e^{-(T-t)\psi \left(\xi \right)}\widehat{G}\left(\xi \right)$ be in $L_1(i\omega +{\mathbb{R}}^n)$ for $\omega \in D^{\prime }\cap D$. Decomposing G in the Fourier integral, substituting into (33), and using Fubini’s theorem to justify the change of the order of integration and taking expectation, we obtain, for any $\omega \in D\cap D^{\prime }$:

$$V(x,t)={\left(2\pi \right)}^{-n}{\int }_{\mathrm{Im}\xi =\omega }{e}^{i\langle x,\xi \rangle -(T-t)(r+\psi (\xi \left)\right)}\widehat{G}\left(\xi \right)d\xi .$$

(35)

The evident representation (35) was derived in [21,53,54]. Using the PDO notation, $V(t,x)=e^{-(T-t)(r+\psi \left(D_x\right))}G\left(x\right)$, where the operator on the RHS acts in appropriate spaces with exponential weights. It is easy to verify that V is a solution of the Cauchy problem

$$({\partial }_t-\psi \left(D_x\right)-r)V(t,x)=0,t<T,x\in {\mathbb{R}}^n,$$

(36)

subject to $V(T,x)=G\left(x\right),x\in {\mathbb{R}}^n,$ and appropriate bounds on V as $x\to \pm \infty$. Equation (36) is understood in the sense of generalized functions, and under certain regularity conditions, the solution is unique. See [53,54] and ([21], Chapt. 15) for details. Thus, in this case, the equivalence of the pricing problem and the Cauchy problem is established.

3.2. Efficient Numerical Realizations in 1D Case

In [21,53,54], the integral on the RHS is numerically realized using the trapezoid rule or Simpson rule, and the standard real-analytical error bounds are used to give prescriptions for the choice of the numerical scheme to satisfy a given error tolerance $ϵ>0$. Since the integrand is analytic in a strip around the line of integration, it is significantly more efficient to use the simplified trapezoid rule, the reason being that the error of the infinite trapezoid rule decays as $exp[-2\pi d/\zeta ]$, where d is the half-width of the strip of analyticity around the line of integration, and $\zeta$ is the step. See, e.g., Thm. 3.2.1 in [55]. Thus, if the strip of analyticity is not too narrow, it is relatively easy to satisfy a very small error tolerance for the discretization error. Note that popular variations of this straightforward approach such as the Carr–Madan method [56] and COS method [57,58,59] introduce additional errors which are difficult to control, and lead to systematic errors in practically important situations. See [2,13,14,16,17,18,19,20,31] for the analysis of errors of the Carr–Madan method and COS.

In many cases of interest, the integrand slowly decays at infinity, and a very large number of terms of the truncated sum (simplified trapezoid rule) are needed to satisfy even a moderate error tolerance. However, in the case of standard European options, and in the case of piece-wise polynomial approximations of complicated payoffs [16,46,60], $\widehat{G}$ is meromorphic with a finite number of simple poles; in [20], approximations with an infinite number of poles appear. If X is SINH-regular of order $(\nu ,{\nu }^{\prime })$ with ${\nu }^{\prime }>0$, one can use an appropriate conformal deformation and the corresponding change of variables to reduce calculations to the case of an integrand which is analytic in a strip around the line of integration and decays at infinity faster than exponentially. The complexity of the numerical scheme based on the sinh-acceleration (15) is of the order of $ElnE$, where $E=ln(1/ϵ)$; in the case of VGP with $\mu =0$, of the order of $O\left(E^2\right)$. See [31] for details. Note that the parameter $\omega$ in (15) is chosen so that the oscillating factor $e^{i(x+(T-t\left)\mu \right)\xi }$ becomes a quickly decaying one. The idea is similar to the idea of the saddle point method. However, simpler universal families of conformal deformations are easier to use, especially when the deformations of several lines of integration are needed, and the deformations must be in a certain agreement [10,61].

3.3. Stable Lévy Processes and Fractional Differential Operators

In the case of stable Lévy processes, the expectation can be finite only if G is bounded or increases at a sufficiently small polynomial rate at infinity. However, if $\widehat{G}$ is analytic in an open cone $\mathcal{C}$ around $\mathbb{R}\backslash \left\{0\right\}$, then the evaluation of probability distribution functions and expectations can be reduced to the evaluation of integrals over $(0,+\infty )$, and an exponential change of variables $\xi =e^{i\omega +y}$, where $\omega \in [-\pi /2,\pi /2]$ is used to efficiently evaluate integrals. In some cases, other conformal deformations are more efficient (see [9]). The same families of conformal deformations can be used to evaluate the solutions of the Cauchy problem for the fractional-parabolic equations of the form:

$$({\partial }_t-|D_x{|}^{\alpha })V(t,x)=0,t<T,x\in \mathbb{R},$$

where $\alpha >0$, and more general ones.

3.4. Calculation of Probability Distributions and Expectations (Prices) in Multi-Factor Lévy Models and Solution of the Cauchy Problems in ${\mathbb{R}}^n$

When n is moderate, X is SINH-regular, and G can be approximated by a function whose Fourier transform $\widehat{G}$ is analytic in the union of a tube domain and cone ${\mathcal{U}}_G\subset iD+\mathcal{C}$, and can be efficiently calculated; then one can construct appropriate conformal deformations of the form (15) for each variable. The complexity of the scheme is $O\left({(ElnE)}^n\right)$, which is not large if, e.g., $n=2,3,4$. For exchange or basket options with payoffs ${(e^{x_1}-e^{x_2}-K)}_{+}$, ${(e^{x_1}+e^{x_2}-K)}_{+}$, $\widehat{G}$ can be explicitly calculated in terms of the Beta function, and the values of the latter can be efficiently calculated using the same sinh-acceleration technique (see [6,7]).

4. Barrier Options in Lévy Models, and Boundary Problems for PDO with Constant Symbols

We consider in detail the 1D case, and in the end, outline extensions to the multi-factor case. We start with the basic notation and facts of the Wiener–Hopf factorization, and then formulate the results for barrier options.

4.1. Main Notation

• X: a Lévy process on $\mathbb{R}$;
• $(\Omega ,\mathcal{F},{\left\{\mathcal{F}\right\}}_{t\ge 0})$: the filtered measure space generated by X;
• $\mathcal{M}$: the set of all stopping times with respect to the filtration ${\left\{\mathcal{F}\right\}}_{t\ge 0}$;
• $\mathbb{Q}$: an EMM chosen for pricing;
• $\psi \left(\xi \right)=-i\mu \xi +{\psi }^0\left(\xi \right)$: the characteristic exponent of X;
• h, $h_{\pm }\in \mathbb{R}$, $h_{-}<h_{+}$: barriers;
• ${\tau }_h^{-}$ and ${\tau }_h^{+}$: first entrance time by X into $(-\infty ,h]$ and $[h,+\infty )$, respectively;
• $q>0$: the discount rate;
• The (expected) present value of the perpetual stream $g\left(X_t\right)$ which is lost at time ${\tau }_h^{-}$:

  $$V_{-}(g;q;h;x)={\mathbb{E}}^x\left[{\int }_0^{{\tau }_h^{-}}{e}^{-qt}g\left(X_t\right)dt\right];$$
• The (expected) present value of the perpetual stream $g\left(X_t\right)$ which is lost at time ${\tau }_h^{+}$:

  $$V_{+}(g;q;h;x)={\mathbb{E}}^x\left[{\int }_0^{{\tau }_h^{+}}{e}^{-qt}g\left(X_t\right)dt\right];$$
• the (expected present) value of the perpetual stream $g\left(X_t\right)$ which is lost at ${\tau }_{h_{+}}^{+}\wedge {\tau }_{h_{-}}^{-}$:

  $$V(g;q;h_{-},h_{+};x)={\mathbb{E}}^x\left[{\int }_0^{{\tau }_h^{+}\wedge {\tau }_{h_{-}}^{-}}{e}^{-qt}g\left(X_t\right)dt\right];$$
• $T_q$: exponentially distributed random variable of mean $1/q$, independent of X;
• ${\overline{X}}_t={sup}_{0\le s\le t}{X}_s$ and ${\underline{X}}_t={inf}_{0\le s\le t}{X}_s$—the supremum and infimum processes (defined pathwise, a.s.);
• Normalized EPV operators (normalized resolvents) under X, $\overline{X}$, and $\underline{X}$ calculate the (normalized) expected present value of the streams under $X,\overline{X}$ and $\underline{X}$:

  $$\begin{array}{ccc}\hfill \left({\mathcal{E}}_{q}g\right)\left(x\right)& :=& {\mathbb{E}}^{\mathbb{Q}}\left[g\left(X_{T_q}\right)\right]:=q{\mathbb{E}}^{\mathbb{Q}}\left[{\int }_0^{+\infty }{e}^{-qt}g\left(X_t\right)dt|X_0=x\right]\hfill \\ \hfill \left({\mathcal{E}}_q^{+}g\right)\left(x\right)& :=& {\mathbb{E}}^{\mathbb{Q}}\left[g\left({\overline{X}}_{T_q}\right)\right]:=q{\mathbb{E}}^{\mathbb{Q}}\left[{\int }_0^{+\infty }{e}^{-qt}g\left({\overline{X}}_t\right)dt|X_0={\overline{X}}_0=x\right]\hfill \\ \hfill \left({\mathcal{E}}_q^{-}g\right)\left(x\right)& :=& {\mathbb{E}}^{\mathbb{Q}}\left[g\left({\underline{X}}_{T_q}\right)\right]:=q{\mathbb{E}}^{\mathbb{Q}}\left[{\int }_0^{+\infty }{e}^{-qt}g\left({\underline{X}}_t\right)dt|X_0={\underline{X}}_0=x\right];\hfill \end{array}$$
• Wiener–Hopf factors: ${\varphi }_q^{+}\left(\xi \right)={\mathbb{E}}^{\mathbb{Q}}\left[e^{i{\overline{X}}_{T_q}}\right]$, ${\varphi }_q^{-}\left(\xi \right)={\mathbb{E}}^{\mathbb{Q}}\left[e^{i{\underline{X}}_{T_q}}\right]$

Basic facts:

(1)

${\varphi }_q^{+}\left(\xi \right)$ (resp., ${\varphi }_q^{-}\left(\xi \right)$) admits (uniformly bounded) analytic continuation to the upper (resp., lower) half-plane.

(2)

The EPV operators act in $L_{\infty }\left(\mathbb{R}\right)$. If $\mathrm{supp}g\subset (-\infty ,h]$, then $\mathrm{supp}{\mathcal{E}}_q^{+}g\subset (-\infty ,h]$. If $\mathrm{supp}g\subset [h,+\infty )$, then $\mathrm{supp}{\mathcal{E}}_q^{-}g\subset [h,+\infty )$.

(3)

${\mathcal{E}}_q=q{(q+\psi \left(D\right))}^{-1}$, ${\mathcal{E}}_q^{+}={\varphi }_q^{+}\left(D\right)$, ${\mathcal{E}}_q^{-}={\varphi }_q^{-}\left(D\right)$. If there exist ${\mu }_{-}\le 0\le {\mu }_{+}$, ${\mu }_{-}<{\mu }_{+}$, such that:

$${\mathbb{E}}^{\mathbb{Q}}\left[e^{\beta X_{T_q}}\right]<\infty ,\beta \in [-{\mu }_{-},-{\mu }_{+}],$$

(37)

then:

(a)

$\psi$ admits (uniformly bounded) analytic continuation to the strip $S_{[{\mu }_{-},{\mu }_{+}]}:=\left\{\xi \right|\mathrm{Im}\xi \in [{\mu }_{-},{\mu }_{+}]\}$ (meaning: analytic in the interior and continuous up to the boundary);

(b)

${\varphi }_q^{+}\left(\xi \right)$ admits analytic continuation to the half-plane $\{\mathrm{Im}\xi \ge {\mu }_{-}\}$;

(c)

${\varphi }_q^{-}\left(\xi \right)$ admits analytic continuation to the half-plane $\{\mathrm{Im}\xi \le {\mu }_{+}\}$;

(d)

the action of the EPV operators extends to $L_{\infty }$ and Sobolev spaces with exponential weights.

The properties (a)–(d) are used to study the asymptotics of prices of barrier options (solutions of the boundary problems) at the boundary [21,62,63], and develop efficient numerical methods for pricing barrier options, credit default swaps (CDSs) and lookbacks [10,61,64,65].

4.2. Wiener–Hopf Factorization

We use three equivalent versions of the Wiener–Hopf factorization

$$\begin{array}{ccc}\hfill \mathbb{E}\left[e^{i\xi X_{T_q}}\right]& =& \mathbb{E}\left[e^{i\xi {\overline{X}}_{T_q}}\right]\mathbb{E}\left[e^{i\xi {\underline{X}}_{T_q}}\right],\xi \in \mathbb{R};\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill \frac{q}{q+\psi \left(\xi \right)}& =& {\varphi }_q^{+}\left(\xi \right){\varphi }_q^{-}\left(\xi \right),\xi \in \mathbb{R};\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill {\mathcal{E}}_q& =& {\mathcal{E}}_q^{-}{\mathcal{E}}_q^{+}={\mathcal{E}}_q^{+}{\mathcal{E}}_q^{-}.\hfill \end{array}$$

(40)

Equation (39) is a special case of the initial form of the Wiener–Hopf factorization used in complex analysis since [66]; (40), with the interpretation of the EPV operators as PDO, is used in analysis (see, e.g., [26]). In both cases, the derivation is possible under certain regularity conditions on $\psi$. The probabilistic versions (38) and (40), hence, (39), hold for any Lévy process. In probability, the straightforward and short proof of (38) is based on

Lemma 1

([67], Lemma 2.1, and [68], p. 81). Let X and $T_q$ be as above. Then:

(a) 

the random variables ${\overline{X}}_{T_q}$ and $X_{T_q}-{\overline{X}}_{T_q}$ are independent; and

(b) 

the random variables ${\underline{X}}_{T_q}$ and $X_{T_q}-{\overline{X}}_{T_q}$ are identical in law.

Equation (40) is a special case of the next lemma derived in [21,22,49] under unnecessary restrictions on the process, and in [46], for any Lévy process. The proof below (in a shortened form) is borrowed from [46].

Theorem 2.

Let X be a Lévy process. Then, for any $g\in L_{\infty }\left(\mathbb{R}\right)$ and $h\in \mathbb{R}$,

$$\begin{array}{ccc}\hfill V_{-}(g;q;h;·)& =& q^{-1}{\mathcal{E}}_q^{-}{\mathbf{1}}_{(h,+\infty )}{\mathcal{E}}_q^{+}g,\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill V_{+}(g;q;h;·)& =& q^{-1}{\mathcal{E}}_q^{+}{\mathbf{1}}_{(-\infty ,h)}{\mathcal{E}}_q^{-}g.\hfill \end{array}$$

(42)

Proof. 

Let X, $\overline{X}$ and $\underline{X}$ start at 0. Then: ${\mathcal{E}}_{q}g\left(x\right)=\mathbb{E}\left[g(x+X_{T_q})\right],$

$${\mathcal{E}}^{+}g\left(x\right)=\mathbb{E}\left[g(x+{\overline{X}}_{T_q})\right],{\mathcal{E}}^{-}g\left(x\right)=\mathbb{E}\left[g(x+{\underline{X}}_{T_q})\right],$$

and, by definition:

$$\begin{array}{ccc}\hfill V_{-}(g;q;h;x)& :=& \mathbb{E}\left[{\int }_0^{{\tau }_h^{-}}{e}^{-qt}g(x+X_t)dt\right]\hfill \\ & =& \mathbb{E}\left[{\int }_0^{\infty }{\mathbf{1}}_{x+{\underline{X}}_t>h}{e}^{-qt}g(x+X_t)dt\right]\hfill \\ & =& q^{-1}\mathbb{E}\left[g(x+X_T){\mathbf{1}}_{x+{\underline{X}}_T>h}\right].\hfill \end{array}$$

Applying Lemma 1, we continue:

$$\begin{array}{ccc}\hfill V_{-}(q;h;x)& =& q^{-1}\mathbb{E}\left[g(x+{\underline{X}}_T+X_T-{\underline{X}}_T){\mathbf{1}}_{x+{\underline{X}}_T>h}\right]\hfill \\ & =& q^{-1}\mathbb{E}\left[{\mathbf{1}}_{x+{\underline{X}}_T>h}{\mathcal{E}}_q^{+}g(x+{\underline{X}}_T)\right]\hfill \\ & =& q^{-1}{\mathcal{E}}_q^{-}{\mathbf{1}}_{(h,+\infty )}{\mathcal{E}}_q^{+}g\left(x\right).\hfill \end{array}$$

This proves (41). Proof of (42) is by symmetry. □

Remark 2.

If (37) is satisfied, one can allow for exponentially increasing measurable $g\left(x\right)$ if:

$$\left|g\left(x\right)\right|\le C(e^{-{\mu }_{-}^{\prime }x}+e^{-{\mu }_{+}^{\prime }x}),$$

(43)

where ${\mu }_{-}\le {\mu }_{-}^{\prime }\le {\mu }_{+}^{\prime }\le {\mu }_{+}$, and

$$q+\psi (-i{\mu }_{\pm }^{\prime })>0.$$

(44)

For the proof of Theorem 2 in this case, it suffices to consider a non-negative measurable function g. We approximate g with the sequence $g_n\left(x\right)=min\{g\left(x\right),n\}$, apply Theorem 2 to $g_n$ and justify passage to the limit using the dominated convergence theorem.

The boundary problem for $V_{-}(q;x;h)$, in the PDO notation, is

$$(q+\psi \left(D_x\right))V_{-}(g;q;h;x)=g\left(x\right),x>h,$$

(45)

subject to $V_{-}(g;q;h;x)=0,x\le h$. Under certain regularity conditions on $\psi$ and g, the standard analytical technique [26] can be applied, and the existence and uniqueness of solutions in a class of bounded sufficiently regular functions proved. In [21], we derived

$$V_{-}(g;q;h;·)=q^{-1}{\varphi }_q^{-}\left(D\right){\mathbf{1}}_{(h,+\infty )}{\varphi }_q^{+}\left(D\right)g,$$

(46)

which is identical to (41); the (ACP)-condition was used. Thus, under a weak regularity condition, the expectation of the stochastic integral is the unique solution of the corresponding boundary problem, and vice versa.

The following theorem is a trivial corollary of Theorem 2 and (40).

Theorem 3.

Let X be a Lévy process satisfying (37), and let G admit the representation $G=q^{-1}{\mathcal{E}}_{q}g$, where g is a measurable function satisfying (43) and (44).

Then:

(a) 

${\left({\mathcal{E}}_q^{-}\right)}^{-1}G:={\mathcal{E}}_q^{+}g$ and ${\left({\mathcal{E}}_q^{+}\right)}^{-1}G:={\mathcal{E}}_q^{-}g$ are measurable functions satisfying (43);

(b) 

For any $h\in \mathbb{R}$:

$$\begin{array}{ccc}\hfill V_{-,\mathrm{inst}}(G;q;x;h)& :=& {\mathbb{E}}^{\mathbb{Q}}\left[e^{-q{\tau }_h^{-}}G\left(X_{{\tau }_h^{-}}\right)|X_0=x\right]=\left({\mathcal{E}}_q^{-}{\mathbf{1}}_{(-\infty ,h]}{\left({\mathcal{E}}_q^{-}\right)}^{-1}G\right)\left(x\right),\hfill \\ \hfill V_{+,\mathrm{inst}}(G;q;x;h)& :=& {\mathbb{E}}^{\mathbb{Q}}\left[e^{-q{\tau }_h^{+}}G\left(X_{{\tau }_h^{+}}\right)|X_0=x\right]=\left({\mathcal{E}}_q^{+}{\mathbf{1}}_{(-\infty ,h]}{\left({\mathcal{E}}_q^{+}\right)}^{-1}G\right)\left(x\right).\hfill \end{array}$$

4.3. Single Barrier Options

Let X be a Lévy process on $\mathbb{R}$ under an EMM $\mathbb{Q}$ chosen for pricing, $r\ge 0$ the riskless rate, and h the barrier. We consider the case of the down-and-out option, which expires if X enters $(-\infty ,h]$ before the maturity date T. If the barrier is not breached, the payoff is $G\left(X_T\right)$. Applying Fubini’s theorem and (41), we calculate the Laplace transform of the price $V_{n.t.d.}(G;h;T,x)$ with respect to T:

$$\begin{array}{ccc}\hfill {\tilde{V}}_{\mathrm{n}.\mathrm{t}.}(G;h;q,x)& =& {\int }_0^{+\infty }{e}^{-qT}{V}_{n.t.d.}(G;h;T,x)dT\hfill \\ & =& e^{-rT}{q}^{-1}{\mathbb{E}}^{\mathbb{Q}}\left[{\mathbf{1}}_{{\underline{X}}_{T_q}>h}G\left(X_{T_q}\right)|X_0=x\right]\hfill \\ & =& e^{-rT}{q}^{-1}\left({\mathcal{E}}_q^{-}{\mathbf{1}}_{(h,+\infty )}{\mathcal{E}}_q^{+}G\right)\left(x\right).\hfill \end{array}$$

Applying the inverse Laplace transform and using (41) and (46), we obtain, for any $\sigma >0$:

$$\begin{array}{ccc}\hfill V_{\mathrm{n}.\mathrm{t}.}(G;h;T;x)& =& \frac{e^{-rT}}{2\pi i}{\int }_{\mathrm{Re}q=\sigma }{e}^{qT}{q}^{-1}\left({\mathcal{E}}_q^{-}{\mathbf{1}}_{(h,+\infty )}{\mathcal{E}}_q^{+}G\right)\left(x\right)dq\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill & =& \frac{e^{-rT}}{2\pi i}{\int }_{\mathrm{Re}q=\sigma }{e}^{qT}{q}^{-1}\left({\varphi }_q^{-}\left(D\right){\mathbf{1}}_{(h,+\infty )}{\varphi }_q^{+}\left(D\right)G\right)\left(x\right)dq.\hfill \end{array}$$

(48)

The representation (47) is derived in [46] for arbitrary Lévy process. In [21,22], (48), is derived solving the corresponding boundary problem

$$\begin{array}{ccc}\hfill ({\partial }_t+L-q)V(t,x)& =& 0,x>h,t<T,\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill V(T,x)& =& G\left(x\right),\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill V(t,x)& =& G_b\left(x\right),x\le h,t<T,\hfill \end{array}$$

(51)

under certain regularity condition on V; X satisfies the (ACT)-condition—the ambiguity of the specification of the payoff at $(T,x)$, $x\le h$, is irrelevant because Lévy processes are stochastically continuous, and, therefore, $\mathbb{Q}[{\tau }_h^{-}=T]=0$. Thus, under weak regularity conditions (the condition on the process is stronger than in the time-independent case), the expectation of the stochastic expression defining the price is the unique solution of the corresponding boundary problem, and vice versa.

In [21,22,46,61,64], more general classes of single-barrier options, with payoff streams during the lifetime of the option, and non-zero payoffs $G_b\left(X_{{\tau }_h^{-}}\right)$ at time ${\tau }_h^{-}<T$ are considered and pricing formulas derived. In [10], similar general formulas are derived when the payoff $G_b$ depends on t and x.

4.4. Numerical Realizations

If $G_b$ depends on t and x, the explicit formula is a quadruple integral, the explicit form of (48) is a triple integral, and in addition, one needs to evaluate the Wiener–Hopf factors for all dual variables arising in the pricing formula. Thus, efficient calculations are difficult. When one uses the Gaver–Stehfest algorithm (see, e.g., [10,61,64]), only a positive q appears, but the coefficients in the approximate Laplace inversion formula are very large. In the result, in many cases of interest, high precision arithmetic is necessary to ensure the stability of results not to mention their accuracy (see [61,64] for examples). One can avoid large coefficients using the method of lines. In the probabilistic interpretation, time to maturity is approximated by a sum of exponentially distributed random variables, which results in a sequence of time-independent problems (maturity randomization or Carr’s randomization) (see [46,69] for efficient numerical realizations of this idea in applications to single and double barrier options, respectively). Appropriate conformal deformations of the lines of integration lead to faster and more accurate numerical algorithms [10,61]. We illustrate the choice of deformations with the simplest example of the no-touch digital option $V_{\mathrm{n}.\mathrm{t}.}(1;h;T;x)$ (the payoff at maturity is 1). Then, ${\mathcal{E}}_q^{+}G={\varphi }_q^{+}\left(D\right)G=1$, and (47)–(48) simplify. Assuming that (37) holds, we take a sufficiently small ${\omega }_0<0$, and write (48) as

$$V_{\mathrm{n}.\mathrm{t}.}(1;h;T;x)=\frac{e^{-rT}}{{\left(2\pi \right)}^{2}i}{\int }_{\mathrm{Re}q=\sigma }{e}^{qT}{q}^{-1}{\int }_{\mathrm{Im}\xi ={\omega }_0}{e}^{i(x-h)\xi }\frac{{\varphi }_q^{-}\left(\xi \right)}{-i\xi }d\xi dq.$$

(52)

The set of admissible ${\omega }_0<0$ is determined by the properties of the Wiener–Hopf factors; general formulas for the latter are well known. For efficient numerical realizations, the characteristic exponent should admit an analytic continuation to a strip $\{\mathrm{Im}\xi \in ({\lambda }_{-},{\lambda }_{+})\}$, ${\lambda }_{-}<0<{\lambda }_{+}$, around the real axis [21]. Then, (see [21,64]), for any $q>0$:

(I)

There exist ${\sigma }_{-}\left(q\right)<0<{\sigma }_{+}\left(q\right)$ such that:

$$q+\psi \left(\eta \right)\notin (-\infty ,0],\mathrm{Im}\eta \in ({\sigma }_{-}\left(q\right),{\sigma }_{+}\left(q\right));$$

(53)

(II)

The Wiener–Hopf factor ${\varphi }_q^{+}\left(\xi \right)$ admits analytic continuation to the half-plane $\{\mathrm{Im}\xi >{\sigma }_{-}\left(q\right)\}$, and can be calculated as follows—for any ${\omega }_{-}\in ({\sigma }_{-}\left(q\right),min\{\mathrm{Im}\xi ,{\sigma }_{+}\left(q\right)\})$:

$${\varphi }_q^{+}\left(\xi \right)=exp\left[\frac{1}{2\pi i}{\int }_{\mathrm{Im}\eta ={\omega }_{-}}\frac{\xi ln(q+\psi (\eta \left)\right)}{\eta (\xi -\eta )}d\eta \right];$$

(54)

(III)

The Wiener–Hopf factor ${\varphi }_q^{-}\left(\xi \right)$ admits analytic continuation to the half-plane $\{\mathrm{Im}\xi <{\sigma }_{+}\left(q\right)\}$, and can be calculated as follows—for any ${\omega }_{+}\in (max\{\mathrm{Im}\xi ,{\sigma }_{-}\left(q\right)\},{\sigma }_{+}\left(q\right))$:

$${\varphi }_q^{-}\left(\xi \right)=exp\left[-\frac{1}{2\pi i}{\int }_{\mathrm{Im}\eta ={\omega }_{+}}\frac{\xi ln(q+\psi (\eta \left)\right)}{\eta (\xi -\eta )}d\eta \right].$$

(55)

Analytic continuation with respect to $q,\xi$ is possible, and conformal deformation of the contours of integration are possible as well.

Assume that X is SINH-regular. Then, a numerical realization of (52) is designed choosing deformations ${\mathcal{L}}^{\left(1\right)}$, ${\mathcal{L}}^{\left(2\right)}$ and ${\mathcal{L}}^{\left(3\right)}$ of the lines of integration $\{\mathrm{Im}\xi ={\omega }_0\}$, $\{\mathrm{Re}q=\sigma \}$ in (52) and the line of integration $\{\mathrm{Im}\eta ={\omega }_{+}\}$ in (55). We use deformations of the form ${\mathcal{L}}^{\left(1\right)}={\chi }_{{\omega }_1,b,\omega }\left(\mathbb{R}\right)$, ${\mathcal{L}}^{\left(2\right)}={\chi }_{{\omega }_1^{\prime },b^{\prime },{\omega }^{\prime }}\left(\mathbb{R}\right)$, where the function ${\chi }_{{\omega }_1,b,\omega }$ is defined by (15), and ${\mathcal{L}}^{\left(3\right)}={\chi }_{L;\sigma ,b_l,{\omega }_l}\left(R\right)$, where ${\chi }_{L;\sigma ,b_l,{\omega }_l}$ is defined by

$${\chi }_{L;\sigma ,b_l,{\omega }_l}\left(y\right)=\sigma +i{b}_{l}sinh(i{\omega }_l+y).$$

(56)

Since $x-h>0$, the oscillating factor $e^{i(x-h)\xi }$ quickly decays if we choose $\omega >0$; hence, the wings of the contour point upward. The contour ${\mathcal{L}}^{\left(2\right)}$ must be above ${\mathcal{L}}^{\left(1\right)}$ so that $\xi -\eta$ is separated from 0 for all $\xi \in {\mathcal{L}}^{\left(1\right)}$ and $\eta \in {\mathcal{L}}^{\left(2\right)}$. In particular, ${\omega }^{\prime }\ge \omega$. The upper bound on ${\omega }^{\prime }$ is implied by the requirement that ${\mathcal{L}}^{\left(2\right)}$ be in the domain of analyticity of $\psi$. Finally, ${\mathcal{L}}^{\left(3\right)}$ must be chosen so that $q+\psi \left(\eta \right)\notin (-\infty ,0]$ for all $\eta \in {\mathcal{L}}^{\left(2\right)}$ and $q\in {\mathcal{L}}^{\left(3\right)}$.

4.5. Double Barrier Options

The following procedure is borrowed from [69]. Consider knock-out double-barrier options with barriers $h_{-}<h_{+}$; the option expires worthless if X leaves $(h_{-},h_{+})$ before maturity date T; if $X_t\in (h_{-},h_{+})$ until T, the option payoff at maturity is $G\left(X_T\right)$. Assuming the constant riskless rate, the option value at time 0 is:

$$V_{k.o.}(G;h_{-},h_{+};T,x)=e^{-rT}{\mathbb{E}}^{\mathbb{Q}}\left[{\mathbf{1}}_{{\tau }_{h_{-}}^{-}\wedge {\tau }_{h_{+}}^{+}>T}G\left(X_T\right)\right],$$

(57)

and the corresponding boundary problem is:

$$\begin{array}{ccc}\hfill ({\partial }_{\tau }+\psi \left(D_x\right)+r)V(\tau ,x)& =& 0,\tau >0,x\in (h_{-},h_{+}),\hfill \end{array}$$

(58)

$$\begin{array}{ccc}\hfill V(0,x)& =& {\mathbf{1}}_{(h_{-},h_{+})}\left(x\right)G\left(x\right),\hfill \end{array}$$

(59)

$$\begin{array}{ccc}\hfill V(\tau ,x)& =& 0,\tau \ge 0,x\notin (h_{-},h_{+}).\hfill \end{array}$$

(60)

As in the case of single barrier options, the Laplace transform (maturity randomization) reduces the calculation of the expectation on the RHS of (57) to the evaluation of the perpetual stream

$$W_{k.o.}(G;h_{-},h_{+};q,x)={\mathbb{E}}^{\mathbb{Q}}\left[{\int }_0^{{\tau }_{h_{-}}^{-}\wedge {\tau }_{h_{+}}^{+}}{e}^{-qt}G\left(X_t\right)dt\right],$$

(61)

which is abandoned at time ${\tau }_{h_{-}}^{-}\wedge {\tau }_{h_{+}}^{+}$. In [69], $W_{k.o.}(G;h_{-},h_{+};q,x)$ was evaluated as follows:

$$W_{k.o.}(G;h_{-},h_{+};q,x)=G^0\left(x\right)-G_{+}^1\left(x\right)-G_{-}^1\left(x\right)+G_{+}^2\left(x\right)+G_{-}^2\left(x\right)-G_{+}^3\left(x\right)-G_{-}^3\left(x\right)+\cdots$$

(62)

where:

$$G_{+}^0\left(x\right)=G^0\left(x\right){|}_{[h_{+},+\infty )},G_{-}^0\left(x\right)=G^0\left(x\right){|}_{(-\infty ,h_{-}]},$$

$$G_{+}^n\left(x\right)={\mathcal{E}}_q^{-}\left({\mathbf{1}}_{(-\infty ,h_{-}]}\left(x\right)·\left({\left({\mathcal{E}}_q^{-}\right)}^{-1}{G}_{-}^{n-1}\right)\left(x\right)\right)\forall n\ge 1,$$

$$G_{-}^n\left(x\right)={\mathcal{E}}_q^{+}\left({\mathbf{1}}_{[h_{+},+\infty )}\left(x\right)·\left({\left({\mathcal{E}}_q^{+}\right)}^{-1}{G}_{+}^{n-1}\right)\left(x\right)\right)\forall n\ge 1.$$

In the case of HEJD, the series can be explicitly calculated [70,71]. For general Lévy processes, [69] uses a general numerical method based on the piece-wise linear interpolation of functions $G_{\pm }^n$, fast convolution and the refined version of the fast Fourier transform (FFT) technique developed earlier in [46]. It is demonstrated that the standard version of FFT and fractional FFT are either inaccurate or inefficient for many classes of Lévy processes.

In [64], the calculations are in the dual space and fractional-parabolic deformations of the contours of integration are used; the more efficient sinh-acceleration technique can be applied in the same vein.

Under additional regularity conditions on $\psi$, it is possible to prove that the boundary problem (58)–(60) has a unique solution in the class of bounded functions continuous on $[0,+\infty )\times (h_{-},h_{+})$. Then, if X satisfies the (ACT)-condition, the expectation given by (57) is this unique solution.

4.6. Regularity of Solutions of Boundary Problems

Assume that ${\psi }^0$ is an elliptic symbol of order $\nu \in (0,2)$, hence, X is a pure jump process. Then, it is easy to prove (see [21]) that:

(a)

If $\nu \in (1,2)$ or $\nu \in (0,1]$ and $\mu =0$, then ${\varphi }_q^{\pm }$ are elliptic symbols of order ${\kappa }^{\pm }=-\nu /2$;

(b)

If $\nu =1$, then ${\varphi }_q^{\pm }$ are elliptic symbols of order ${\kappa }_q^{\pm }$, where ${\kappa }^{\pm }\in (-1,0)$ depend on $\mu$, and ${\kappa }_q^{+}+{\kappa }_q^{-}=-1$;

(c)

If $\nu \in (0,1)$ and $\mu >0$, then ${\varphi }_q^{+}\left(\xi \right)={(q-i\mu \xi )}^{-1}(1+{\alpha }_q^{+}\left(\xi \right))$ and ${\varphi }_q^{-}\left(\xi \right)=1+{\alpha }_q^{-}\left(\xi \right)$, where ${\alpha }_q^{\pm }$ are of order $-1+\nu +ϵ$, for any $ϵ>0$. Hence, ${\kappa }_q^{+}:=\mathrm{ord}{\varphi }_q^{+}=-1$, and ${\kappa }_q^{-}:=\mathrm{ord}{\varphi }_q^{-}=0$;

(d)

If $\nu \in (0,1)$ and $\mu <0$, then ${\varphi }_q^{-}\left(\xi \right)={(q-i\mu \xi )}^{-1}(1+{\alpha }_q^{-}\left(\xi \right))$ and ${\varphi }_q^{+}\left(\xi \right)=1+{\alpha }_q^{+}\left(\xi \right)$, where ${\alpha }_q^{\pm }$ are of order $-1+\nu +ϵ$, for any $ϵ>0$. Hence, ${\kappa }_q^{+}:=\mathrm{ord}{\varphi }_q^{+}=0$, and ${\kappa }_q^{-}:=\mathrm{ord}{\varphi }_q^{-}=-1$.

Straightforward calculations based on (a)–(d) show (see [21,62,63]) that the prices of down-and-out barrier options—hence, solutions of the corresponding boundary problems—are not smooth (and in the case of (c), discontinuous) at the barrier; the leading term of asymptotics as $x\downarrow h$ is calculated. In [63], the asymptotics of the derivatives of the price $V(t,x)$ (sensitivities) is calculated as well.

In cases (a), (b), (c), $V(\tau ,x)\sim c\left(\tau \right){(x-h)}^{{\kappa }_q^{-}}$, and $V_x(\tau ,x)\to +\infty$ as $x\downarrow h$. In case (d), $V_x$ is continuous up to the boundary but $V_{xx}$ is unbounded.

The case of the BM with embedded jumps was not explicitly studied but the study can be performed in a similar vein. If the jump component is of order $\nu \in [1,2)$, then the solution is continuous up to the boundary but $V_x$ is discontinuous; if $\nu \in (0,1)$, then $V_x$ is continuous but $V_{xx}$ is unbounded.

The case of the double-barrier options was not studied but since the study of the regularity at the boundary is easily localized (see [26]), one can easily prove that, for instance, in case (c), the solution is discontinuous at $h_{-}$ but smooth at $h_{+}$; $V_x$ is unbounded as $x\to h_{-}$, and $V_{xx}$ is unbounded as $x\to h_{+}$.

4.7. The Case of Time-Dependent Boundaries

The regularity of solutions is an open problem. In this section, we formulate several natural hypotheses. Assuming that the boundary (or two boundaries) are piece-wise smooth, one can try to study the regularity of solutions localizing the problem. In the case of processes of order $\nu \in [1,2)$, the operator $A(D_t,D_x)$ is quasi-elliptic (the symbol $\widehat{A}(\eta ,\xi )=i\eta -i\mu \xi +{\psi }^0\left(\xi \right)$ satisfies $|\widehat{A}{(\eta ,\xi )|\ge c(\left|\eta \right|+\left|\xi \right|}^{\nu })-C$, for $(\eta ,\xi )\in {\mathbb{R}}^2$), and elliptic if $\nu =1$; hence, localization can be performed as in the elliptic case [26]. If $\nu \in (1,2)$, one expects that, if t is fixed and $(t,x)$ tends to a point $(t,x_0)$ on a smooth part of the boundaries, the solution behaves as $c\left(t\right)|x-x_0{|}^{\nu /2}$. If $\nu =1$, then the study becomes complicated. One expects that the asymptotics at the boundary can be naturally described in the coordinates $x=x^{\prime }+f\left(t\right)$, which makes the boundary locally flat: $x^{\prime }=h^{\prime }$. For each point t in a small neighborhood of $t_0$, the asymptotics of the price $V(t,x^{\prime })$ in the new coordinate system is of the form $c\left(t\right)|x^{\prime }{|}^{{\kappa }^{\pm }\left(t\right)}$, where ${\kappa }^{\pm }\left(t\right)\in (0,1)$ (“-” for the lower boundary and “+” for the upper boundary) continuously depend on t and are defined by the “drift” ${\mu }^{\prime }(t,h^{\prime })$ in the new coordinate system and ${\psi }^0$ as in case (b) above.

If $\nu \in (0,1)$, then the localization itself becomes more difficult because the principal symbol $-i\eta -i\mu \xi$ is hyperbolic. Assuming that the localization is possible, one expects that the solution is smooth up to the boundary (or discontinuous at the boundary) when the vector field $-{\partial }_{\tau }+\mu {\partial }_x$ is transversal to the boundary and points to (or from, respectively) the boundary.

4.8. Multi-Factor Case

Consider first a boundary problem in $(0,+\infty )\times {\mathbb{R}}_{x^{\prime }}^{n-1}\times (0,+\infty )$. In the analytical setting, one makes the Fourier transform with respect to $x^{\prime }$ and solves the family of problems on $\{\tau >0,x_n>0\}$ using the results above. However, then the dependence of the Wiener–Hopf factors on ${\xi }^{\prime }$ does not allow for the interpretation of the Wiener–Hopf factors as the symbols of the EPV operators, and one is forced to impose additional potentially unnecessary conditions in order to justify the results.

A natural alternative is to represent the payoff functions as sums of functions supported on direct products of the half-axis. Then, after an appropriate change of variables, each new problem is a problem on ${(0,+\infty )}^{n+1}$. Then, we make the Laplace transform with respect to $\tau$ and $x_j,j=1,2,\dots ,n-1$, and solve the resulting family of problem on $(0,+\infty )$ depending on $q\in {(0,+\infty )}^n$ exactly as in the case of Lévy processes on $\mathbb{R}$. The characteristic exponent of the one-dimensional process is $\psi (-iq,{\xi }_n)$, and the resulting analytical expression admits analytic continuation to $\{q\in {\mathbb{C}}^n|\mathrm{Re}{q}_j>0\}$ (and wider regions, if the process is sufficiently regular).

In this way, one can derive explicit formulas for single barrier options, and the representation of the price of double-barrier options with flat parallel barriers. In the application to credit risk models (counterparty risk), important variations are problems with boundaries $x_j=h_j$, $j=n,n-1$, when $n=2$, and even $j=n,n-1,n-2$, when $n=3$. The case $n=2$ can be solved modifying in the straightforward fashion the iterative procedure used in the case of two parallel flat boundaries. The proof and design of efficient numerical procedures become more involved but, at the theoretical level, the tools remain essentially the same. The case $n=3$ becomes messier still but an iteration procedure can be designed, and convergence to the price proved as in the case of two boundaries.

The case of curved boundaries is similar to the case $n=1$, and naturally, the results are at the level of hypotheses so far.

5. American Options and Free Boundary Problems

5.1. Basic Example

The European options can be exercised only at expiry: if $S_T=e^{X_T}\ge K$, then it is optimal to exercise the European call option and receive $e^{X_T}-K$; otherwise the option expires and is worthless, and the payoff is 0. An American option can be exercised at any moment until the expiry date, T. If the process X is Markovian, it is natural to expect that there is a subset B of the half-space $\left\{\right(t,x\left)\right|t\le T\}$ such that the option is exercised the first time $(t,X_t)\in B$. B is called the exercise region and the part of the boundary $\partial B$ that is in the open half-plane $\left\{\right(t,x\left)\right|t<T\}$ is called the early exercise boundary. In the simplest cases of the American option put and call options, the early exercise regions are connected, and of the form $\left\{(t,x)\right|x\le h_{*}\left(t\right),t<T\}$ and $\left\{(t,x)\right|x\ge h^{*}\left(t\right),t<T\}$, respectively.

In the Brownian motion case, the early exercise boundary and option value can be found by solving the corresponding free boundary problem. In addition to the terminal condition at $t=T$, one adds the value matching and smooth pasting condition at the early exercise boundary. Let $r\ge 0$ be the riskless rate, $\delta \ge 0$ the dividend rate, and $X_t$ be a Lévy process under the risk-neutral measure chosen for pricing, with the infinitesimal generator L. For the American put option in the one-factor model, the free boundary problem is of the form

$$\begin{array}{ccc}\hfill ({\partial }_t+L-r+\delta )V(t,x)& =& 0,x>h\left(t\right),t<T,\hfill \end{array}$$

(63)

$$\begin{array}{ccc}\hfill V(T,x)& =& {(K-e^x)}_{+},\hfill \end{array}$$

(64)

$$\begin{array}{ccc}\hfill V(t,x)& =& K-e^x,x=h_{*}\left(t\right),t<T,\hfill \end{array}$$

(65)

$$\begin{array}{ccc}\hfill V_x(t,x)& =& -e^x,x=h_{*}\left(t\right),t<T,\hfill \end{array}$$

(66)

and V is sought in the class of continuous bounded functions, which are of the class $C^{1,2}$ above the early exercise boundary, with the first derivative $V_x$ continuous up to the boundary (these regularity conditions for the American option price are well known and proven). On the RHS of (65), we write $K-e^x$ instead of ${(K-e^x)}_{+}$ because it is non-optimal to exercise the option unless the payoff is non-negative. The free boundary problem for the American call option is similar. The equivalence of the pricing problems for the American put and call options and the corresponding free boundary problems is proven in the Brownian motion case and for Brownian motion with embedded compound Poisson process. For general Lévy processes, the equivalence is an open problem. Note that for processes with jumps, the value matching condition (65) becomes non-local:

$$V(t,x)=K-e^x,x\le h_{*}\left(t\right),t<T,$$

(67)

and the smooth pasting condition may fail. We will not use the smooth pasting condition. Instead, we look for the boundary which maximizes the solution of the problem (63), (64), (67).

5.2. Behavior of the Early Exercise Boundary near Maturity

The next subtle point is the limit $h^{*}(T-):={lim}_{t\to T}{h}_{*}\left(t\right)$. Designing numerical methods, one is tempted to assume that, as $t\to T$, the region in the state space $U_t$ where the American option with the payoff $g\left(X_t\right)$ remains alive, tends to $\left\{x\right|g\left(x\right)\ge 0\}$. In the case of the put option above, this assumption translates into the condition $h^{*}(T-0)=lnK$: the limit of the early exercise boundary at maturity (in the S-coordinate) is equal to the strike. This statement is correct if $r>0$ and $\delta =0$. If the stock pays dividends, then it is possible that $h_{*}(T-)<lnK$; in the case $r=0$, simple general no-arbitrage considerations spelled out by Merton prove that it is non-optimal to exercise the American put before maturity. Essentially the same no-arbitrage argument proves that it is non-optimal to exercise the American call on non-dividend paying stock before maturity. Starting with [72,73], the behavior of the critical stock price near maturity for American options in diffusion models has been studied in a number of publications (see the bibliography in [74]). It is proved, in particular, that if the stock pays dividends, then, depending on the parameters ${\sigma }^2,r,\delta$, the limit of the early exercise boundary for the American call in the Black–Scholes model is above the strike.

In [21,40], it is proven that in the presence of positive jumps, the early exercise boundary for the American put without dividends is separated from the strike by a margin. In [21], the result is obtained for KoBoL, NTS and NIG models using the Wiener–Hopf factorization technique. In [40], the proof is based on the calculation of $\theta (t,x):=V_t(t,x)$ of out-of-the-money options, at expiry (in the case of the put, in the region $x>lnK$, in the case of the call, in the region $x<lnK$). It was proved that, in the presence of jumps, ${\theta }_{\mathrm{call}}(T-,x)<0$ for $x<lnK$ (in the diffusion case, ${\theta }_{\mathrm{call}}(T-,x)=0$ for $x<lnK$). Using this result and the put-call parity, it is proven that it is not optimal to exercise the American put without dividends up to expiry in the region where $rK<-{\theta }_{\mathrm{call}}(T-,x)$. It is demonstrated that for parameters’ values documented in empirical studies of financial markets, the margin is several percent of the strike or even more than a dozen percent. Hence, a numerical method is based on the condition $h_{*}(T-)=lnK$ is expected to produce rather inaccurate results. In [24], a similar formula for $\theta$ of out-of-the-money options in one-factor Lévy driven quadratic term structure models (QTSMs) is presented but without the proof.

In [74], the results in [24,40] are generalized for wide classes of multi-factor Markov models with jumps, and the proofs are simplified. First, consider a European option with the (effective) payoff $g_{+}\left(X_T\right)=max\{0,g\left(X_T\right)\}$; for a European call on a stock with the price process $e^{X_t}$, and strike K, $g\left(x\right)=e^x-K$. Let $F(x,dy)$ be the density of jumps of the underlying Markov process $X_t$ and denote by $U_{-}\left(g\right):=\left\{x\right|g\left(x\right)<0\}$ the out-of-the-money region. Denote by $V(g_{+};x,\tau )$ the option price at time $\tau >0$ to expiry and $X_{T-\tau }=x$, and by $\mathcal{C}(g_{+};x)$ the limit of the $-\theta$ of the option, at expiry:

$$\mathcal{C}(g_{+};x)=\underset{\tau \to +0}{lim}\frac{V(g_{+};x,\tau )}{\tau }.$$

(68)

For wide classes of payoffs and jump-diffusion models, the limit exists for $x\in U_{-}\left(g\right)$, and it is given by

$$\mathcal{C}(g_{+};x)=\int g_{+}(x+y)F(x,dy),$$

(69)

for almost all $x\in U_{-}\left(g\right)$ (see [74]). If (69) is applied to digital options, one modifies the definition of the payoff: $g\left(x\right)=1$ in the in-the-money region, and $g\left(x\right)<0$ in the out-of-the-money region.

Let $U_{+}\left(g\right)=U_{-}(-g)$ be the in-the-money region of the option, and ${\Omega }_{\tau }$ be the optimal non-exercise region for the American option, in the $x-$space, at time $\tau$ to expiry. Clearly, ${\Omega }_{\tau }\supset U_{-}\left(g\right)$, hence, ${\cup }_{\tau >0}{\Omega }_{\tau }\supset U_{-}\left(g\right)$. However, it may be the case that it is non-optimal to exercise the option up to expiry even if x is in the in-the-money region $U_{+}\left(g\right)$. Let L be the infinitesimal generator of X. Define: ${\Omega }^J=\{x\in U_{+}\left(g\right)|Lg\left(x\right)+\mathcal{C}({(-g)}_{+},x)>0\}.$ The following two theorems from [74] give the “lower bound” for the limit of the no-exercise region at maturity and approximate formulas for the option price close to maturity, which can be used at the first step of backward induction procedures for pricing American options, in wide classes of Markov models.

Theorem 4

([74], Thm. 2.2). ${\Omega }_{+0}\supset {\Omega }^J\cup U_{-}\left(g\right).$

Theorem 5

([74], Thm. 2.3). Let $Lg\left(x\right)+\mathcal{C}({(-g)}_{+},x)\ne 0$. Then, as $\tau \to 0$:

1. 

For x in the out-of-the-money region, $U_{-}\left(g\right)$:

$$V_{\mathrm{am}}(g_{+};x,\tau )\sim \tau \mathcal{C}(g_{+};x);$$

(70)

2. 

For x in the in-the-money region $U_{+}\left(g\right)$:

$$V_{\mathrm{am}}(g_{+};x,\tau )\sim g\left(x\right)+\tau {(Lg\left(x\right)+\mathcal{C}({(-g)}_{+};x))}_{+}.$$

(71)

5.3. Perpetual American Options in One-Factor Models or Stationary Free Boundary Problems on $\mathbb{R}$

The solution of the optimal stopping problem in Lévy models with infinite time horizon (equivalently, stationary free boundary problems) is the main block for the solution of the problems with finite time horizon using maturity randomization (equivalently, method of lines), in regime-switching models, and the approximation of more general Markov models with regime-switching Lévy models.

The first simple but crucial trick which makes it possible to give simple proofs and design efficient procedures is the reduction in the option to acquire an instantaneous payoff $G\left(X_t\right)$ to the option to abandon the stream$-g\left(X_t\right)$. If q is the discount rate, then we assume that $G=q^{-1}{\mathcal{E}}_{q}g$, where g satisfies (43) and (44). This implies that G is sufficiently regular, does not increase too quickly at infinity, and the discount rate is not too small. In the case of options with finite time horizon, using sufficiently small time intervals in the method of lines, the latter condition can be satisfied in all cases. In the case of the American put with the infinite time horizon, on the non-dividend paying stock, the representation $G=q^{-1}{\mathcal{E}}_{q}g$ is impossible but a representation $G=q^{-1}{\left({\mathcal{E}}_q^{-}\right)}^{-1}g$ is possible and suffices for the proof.

Assuming that $G=q^{-1}{\mathcal{E}}_{q}g$, for any stopping time $\tau$:

$${\mathbb{E}}^x\left[e^{-q\tau }G\left(X_{\tau }\right)\right]=q^{-1}{\mathcal{E}}_{q}g\left(x\right)+V_{\mathrm{ex}}(-g;\tau ;x),$$

where:

$$V_{\mathrm{ex}}(f;\tau ;x)={\mathbb{E}}^x\left[{\int }_{\tau }^{+\infty }{e}^{-qt}f\left(X_t\right)dt\right].$$

Hence, an optimal time to exercise the American option with the payoff function G is an optimal time to abandon the stream ${\{-g\left(X_t\right)\}}_{t\ge 0}$, and vice versa.

Theorem 6.

Let the following conditions hold:

(i) 

X is a Lévy process with the non-trivial infimum process;

(ii) 

g is a non-decreasing stream that changes sign;

(iii) 

Bounds (37), (43) and (44) hold.

Then:

(a) 

There exists h such that:

$${\mathcal{E}}_q^{+}g\left(x\right)\le 0,x\le h,\mathrm{and}{\mathcal{E}}_q^{+}g\left(x\right)\ge 0,x\ge h.$$

(72)

(b) 

${\tau }_h^{-}$, the entry time into $(-\infty ,h]$, is an optimal exit time in class $\mathcal{M}$.

(c) 

The option value can be represented as the EPV of the stream $g\left(X_t\right)-U\left(X_t\right)$, where U is a non-decreasing function vanishing above h.

(d) 

If g is not monotone but (i) holds, then ${\tau }_h^{-}$ is an optimal exit time in the class of stopping times of the threshold type.

Proof. 

(a,d) are immediate from the representation

$$V_{\mathrm{ex}}(g;{\tau }_h^{-};x)=q^{-1}{\mathcal{E}}_q^{-}{\mathbf{1}}_{(h,+\infty )}{\mathcal{E}}_q^{+}g\left(x\right).$$

The expectation operator ${\mathcal{E}}_q^{-}$ being positive, $V_{\mathrm{ex}}(g;{\tau }_h^{-};x)$ is maximized when the function ${\mathbf{1}}_{(h,+\infty )}{\mathcal{E}}_q^{+}g$ is maximized. Hence, all negative values of ${\mathcal{E}}_q^{+}g$ must be set to 0.

(b,c) are directly derived (see [21,49,75]) verifying the conditions of the following general lemma; the proof is much shorter and simpler than purely probabilistic proofs. However, we need to assume that X satisfies the (ACP)-property. □

Lemma 2

([21,49,75]). Let ${\tau }_B$ be the first entrance time into a Borel set B. Let B be a Borel set such that:

$$\begin{array}{ccc}\hfill W_{\mathrm{ex}}(g;{\tau }_B;·)& :=& (q-L)V_{\mathrm{ex}}(g;{\tau }_B;·)\mathrm{is}\mathrm{universally}\mathrm{measurable}\hfill \end{array}$$

(73)

$$\begin{array}{ccc}\hfill W_{\mathrm{ex}}(g;{\tau }_B;f;x)& =& g\left(x\right),x\in \mathbb{R}\backslash B,\mathrm{a}.\mathrm{e}.\hfill \end{array}$$

(74)

$$\begin{array}{ccc}\hfill W_{\mathrm{ex}}(g;{\tau }_B;x)& \ge & g\left(x\right),x\in B,\mathrm{a}.\mathrm{e}.\hfill \end{array}$$

(75)

$$\begin{array}{ccc}\hfill W_{\mathrm{ex}}(g;{\tau }_B;x)& \ge & 0,\forall x.\hfill \end{array}$$

(76)

Then, ${\tau }_B$ maximizes $V_{\mathrm{ex}}(g;{\tau }_B;·)$ in the class $\mathcal{M}$.

Proof. 

Let $\tau$ be a stopping time. Then, using Dynkin’s formula and (74)–(76), we obtain:

$$\begin{array}{ccc}\hfill V_{\mathrm{ex}}(g;{\tau }_B;x)& =& {\mathbb{E}}^x\left[{\int }_0^{\tau }{e}^{-qt}(q-L)V_{\mathrm{ex}}(g;{\tau }_B;X_t)dt\right]\hfill \\ & & +{\mathbb{E}}^x\left[e^{-q\tau }{V}_{\mathrm{ex}}(g;{\tau }_B;X_{\tau })\right]\hfill \\ & \ge & {\mathbb{E}}^x\left[{\int }_0^{\tau }{e}^{-qt}g\left(X_t\right)dt\right].\hfill \end{array}$$

With $\tau ={\tau }_B$, we obtain the equality, which means that ${\tau }_B$ is optimal. □

Remark 3.

In [51], we proved optimality in the class of all stopping times for a wide class of non-monotone payoffs, under the assumption that the jump density is completely monotone. In [52], the results are applied to solve a game-theoretical problem. In [76], the American options with lookback features are studied in cases when the exercise region is discontinuous.

5.4. Good and Bad News Principles and the Failure of the Smooth Pasting Condition

In [77], for special cases, and in [47,49,78] in the general case, the following interpretation of the stopping rule (b) and its mirror reflection for the option to abandon a non-increasing stream were given:

Good news principle. Abandon an increasing stream (exit) when the EPV of the stream under the supremum process becomes negative.

Bad news principle. Acquire an increasing stream (entry) when the EPV of the stream under the infimum process becomes positive.

Failure of the smooth pasting condition [21,75]. If ${\varphi }_q^{-}\left(\xi \right)\sim c>0$ as $\xi \to \pm \infty$, ${\mathcal{E}}_q^{+}g$ is smooth at h, ${\mathcal{E}}_q^{+}g\left(h\right)=0$, and ${\left({\mathcal{E}}_q^{+}g\right)}^{\prime }\left(h\right)>0$, then h is the optimal exit boundary but the value function has a kink at h. For the option to abandon a non-increasing stream, replace ${\varphi }_q^{-},{\mathcal{E}}_q^{-}$ with ${\varphi }_q^{+},{\mathcal{E}}_q^{+}$.

In terms of atoms of the pdf of ${\overline{X}}_{T_q}$ and ${\underline{X}}_{T_q}$, the failure of the smooth pasting condition was reformulated and proven in [79].

5.5. American Options with Finite Time Horizon or Non-Stationary Boundary Problems on $\mathbb{R}$

We have an optimal stopping problem: for $t<T$, find:

$$V(t,x)=\underset{\tau \in \mathcal{M},\tau \le T}{sup}{\mathbb{E}}_t^{\mathbb{Q}}\left[e^{-r(\tau -t)}G\left(X_{\tau }\right)|X_t=x\right].$$

(77)

Carr’s randomization [80] or method of lines approximates $V(t,x)$ with a sequence of perpetual American options. The maturity period is divided into N subintervals, using points $0=t_0<t_1<\cdots <t_N=T$. Each sub-period $[t_s,t_{s+1}]$ is replaced with an exponentially distributed random maturity period $T_s$ with mean ${\Delta }_s=t_{s+1}-t_s$, $s=0,1,\dots ,N-1$. The random variables $T_s,s=0,\dots ,N-1$ and the process X are assumed to be independent. Typically, one takes ${\Delta }_s=T/N$ for all s. If no optimization decision is involved, then we can say that the deterministic maturity date T is replaced with $T^{\prime }=T_0+T_1+\cdots +T_{N-1}$. If all $T_s\sim \mathrm{Exp}N/T$, then $T^{\prime }$ is Erlang-distributed. The same idea can be applied to price barrier and lookback options, with flat boundaries. When optimizing decisions are involved, an accurate formulation of the approximate optimal stopping problem becomes more involved: a time-consistent exercise rule must take into account realizations of $T^s:=T_0+T_1+\cdots +T_s$, $s=0,1,\dots$, hence, the rule must be updated after the arrival of each $T^s$. The convergence of Carr’s randomization procedure for American options in wide classes of Lévy models is proven in [81]. The paper [80] solves the sequence in the BM model using explicit formulas for the boundary problem for second-order differential operators; this technique cannot be applied to other Lévy models. Furthermore, ref. [80] formulates an equivalent form of maturity randomization as the approximation of the American option with a finite-time horizon by the American option with the Erlang-distributed maturity date, and stated that Richardson’s extrapolation can be applied. Both statements are false for American options but hold for barrier options. The convergence of Carr’s randomization approximation is proven in [71] for wide classes of Markov processes (additional conditions on the process are necessary), and Richardson’s extrapolation of arbitrary order is justified in [63] for the value function and its derivatives.

The proof of optimality of the solution of the sequence and pricing procedure simplify if, at each step, the option is reduced to the option to abandon a stream. Then, the conditions of the basic theorems for perpetual options can be easily verified for each option in the sequence. This is the idea of the solution in [21,38,40,49]. For simplicity, consider equidistant dates, and set $\Delta =T/N$, $q=1/\Delta +r$.

The backward induction procedure is as follows.

• Set $V_0=G_{+}$ (the payoff at maturity).
• In the cycle $s=1,2,\dots ,N$, find $V_s$ as the solution to the optimal stopping problem:

  $$V_s\left(x\right)=\underset{\tau \in \mathcal{M}}{sup}{\mathbb{E}}^{\mathbb{Q}}\left[e^{-q\tau }G\left(X_{\tau }\right)+\frac{1}{\Delta }{\int }_0^{\tau }dt{e}^{-qt}{V}_{s-1}\left(X_t\right)|X_0=x\right].$$

  (78)
• $V_N$ is Carr’s randomization approximation to time-0 option price.

The reduction to the sequence of exit problems is as follows.

Set $W_s=V_s-G$, and notice that the maximization of $V_s$ is equivalent to the maximization of $W_s$. To reformulate the optimal stopping problem in terms of $W_s$, we find function $g_q$ such that:

$$G\left(x\right)=q^{-1}{\mathbb{E}}^{\mathbb{Q}}\left[{\int }_0^{+\infty }{e}^{-qt}{g}_q\left(X_t\right)dt|X_0=x\right].$$

For the put, $g_q\left(x\right)=qK-(q+\psi (-i))e^x$, where $\psi$ is the characteristic exponent of the Lévy process X. We have:

$$-G\left(x\right)+{\mathbb{E}}^{\mathbb{Q}}\left[e^{-q\tau }G\left(X_{\tau }\right)\right]=-{\mathbb{E}}^{\mathbb{Q},x}\left[{\int }_0^{\tau }{e}^{-qt}{g}_q\left(X_t\right)dt|X_0=x\right],$$

therefore, the optimal stopping problem in terms of $W_s$ is:

$$W_s\left(x\right)=\underset{\tau \in \mathcal{M}}{sup}{\mathbb{E}}^{\mathbb{Q}}\left[{\int }_0^{\tau }{e}^{-qt}(f\left(X_t\right)+\frac{1}{\Delta }{W}_{s-1}\left(X_t\right))dt|X_0=x\right],$$

(79)

where $f\left(x\right)=-g_q\left(x\right)+(1/\Delta )G\left(x\right)$. Assuming that $r+\psi (-i)\ge 0$, the function:

$$\begin{array}{ccc}\hfill f\left(x\right)& =& -(qK-(q+\psi (-i))e^x)+\frac{1}{\Delta }(K-e^x)=-rK+(r+\psi (-i))e^x\hfill \end{array}$$

is non-decreasing (the case of the put option). The function $W_0={(-G)}_{+}$ is also non-decreasing. Using Theorem 6, we find Carr’s randomization approximation $h_1$ to the early exercise boundary, and prove that $W_1$ is non-decreasing and vanishes below $h_1$.

In the cycle $s=1,2,\dots ,N-1$, using Theorem 6 and the induction assumptions: $W_s$ is non-decreasing and vanishes below $h_s$, we find the approximation to the early exercise boundary $h_{s+1}$ and $W_{s+1}$ and prove that $W_{s+1}$ is non-decreasing and vanishes below $h_{s+1}$. Finally, $V_s=W_s+G$, $s=1,2,\dots ,N$. See ([49], Chapt. 13) for details.

5.6. Shape of the Early Exercise Boundary and Smooth Pasting Condition

Assume that X is the process of finite variation, with non-zero drift $\mu$. Then, one can conjecture that far from maturity, the early exercise boundary is almost flat; hence, the smooth pasting condition fails if the vector field ${\partial }_t+\langle \mu ,{\partial }_x\rangle$ is transversal to the boundary and points from the boundary. One can also conjecture that the smooth pasting condition fails everywhere but this is far from evident in view of the fact that the behavior of the early exercise boundary at maturity can be very irregular, especially in the multi-factor case. One cannot exclude cases when, at some parts of the free boundary, the vector field ${\partial }_t+\langle \mu ,{\partial }_x\rangle$ is transversal to the boundary and points toward the boundary. At these parts, the smooth pasting condition would hold.

6. Barrier Options and American Options in Regime-Switching Lévy Models and Systems of Pseudo-Differential Equations, Approximation of Stochastic Volatility Models and Models with Stochastic Interest Rate

Let $\mathcal{M}$ be a finite state Markov chain with transition rates ${\lambda }_{jk}$, $j,k=1,2,\dots ,M$. Set ${\Lambda }_j={\sum }_{k\ne j}{\lambda }_{jk}$. For each j, let $X^{\left(j\right)}$ be a Lévy process on $\mathbb{R}$ with the characteristic exponent ${\psi }_j$ and infinitesimal generator $L_j$ under a measure ${\mathbb{Q}}_j$. The riskless rate $q_j$, instantaneous payoffs $G_j$, streams of payoffs $g_j$ and early exercise boundary $h_j$ depend on the state. In the case of barrier options, $h_j$ are given, and in the case of American options, $h_j$ are chosen, solving the optimal stopping problem. The infinitesimal generator of the model is a matrix PDO $\mathcal{L}=-\psi \left(D\right)$, where:

$$\psi \left(\xi \right)=\mathrm{diag}{({\psi }_j\left(\xi \right)+{\Lambda }_j)}_{j=1}^M-{\left[{\lambda }_{jk}\right]}_{j\ne k}).$$

If there are no barriers, then, evidently, one can solve the Cauchy problem as in the scalar case; the only difference is that the matrix exponential $exp[-(T-t)(\psi \left(\xi \right)+\mathrm{diag}{\left[q_j\right]}_{j=1}^M])$ appears, and the choice of appropriate contour deformations is more involved. Similarly, in the case of $h_1=\cdots =h_M$, one can use the matrix form of the Wiener–Hopf factorization, and repeat the calculations (in the PDO form) which we used in the no-regime switching case. However, if M is large, then even these theoretically straightforward and simple methods are very difficult for efficient numerical realization. Furthermore, the straightforward methods outlined above are not applicable if the boundaries are different in different states. This is the case when stochastic volatility models and models with stochastic interest rates are approximated by Markov modulated Lévy models and/or American options are priced.

The idea of the approximation is as follows [82,83,84,85,86]. The action of the infinitesimal generator of the Markov process is discretized—replaced by the infinitesimal generator ${\mathcal{L}}_{disc}$ of the Markov chain with an infinite number of states, and then truncated. At the boundary of the truncated discretized state space, it is necessary to impose appropriate conditions so that the truncated operator is the infinitesimal generator of a Markov chain, without killing. Hence, Dirichlet boundary conditions may not be used. Instead, in order to avoid killing, a kind of reflection condition must be imposed: transition rates from each point at the boundary to points in the truncated discretized state space must change. If ${\mathcal{L}}_M$ is a diffusion, the discretization and boundary conditions are straightforward; in the case of processes with jumps, constructions are more involved.

In a number of publications [50,82,83,84,85,86,87], we used the following simple iteration procedure, which we outline below in the case of down-and-out options: the option is exercised in state j when $X^{\left(j\right)}$ breaches the barrier $h_j$ from above. Numerical examples demonstrated the stability of the procedure even for $M>2000$ [86].

A problem with finite time horizon is reduced to a sequence of boundary problems with infinite time horizon using the method of lines (Carr’s randomization). The regularity conditions are as follows. For each j:

(i)

$X_t^{\left(j\right)}$ satisfies the (ACP)-property (needed in the case of American options only);

(ii)

$g_j$ is measurable, non-negative, and does not grow too fast at infinity;

(iiii)

$(q_j+{\Lambda }_j-L_j)G_j$ is continuous, monotone, and does not grow too fast at infinity.

(iv)

In each state, the rate of growth is controlled by conditions (37), (43), with $q_j+{\Lambda }_j$ in place of $q_j$.

Note that the general results for American options are obtained for non-negative $g_j$. However, we have optimal stopping results in the non-regime switching case when the payoff function g may assume negative values. The same technique can be used in regime-switching models, hence, the procedure can be generalized to the case of non-monotone payoff functions, under certain conditions on the payoffs and processes. The case of barrier options with payoff functions that change sign can be reduced to the case of non-negative payoff functions using the linearity of the expectation operator.

At each step of the backward induction, we use the following block:

I.

Reduce the pricing problem to the problem of evaluation of a perpetual stream (in different states, the payoff streams are different).

II.

Assuming that the value functions in each state but state j are known, calculate the state-j option value.

III.

In the case of American options, calculate the approximation to the early exercise boundary in state j.

IV.

Using this conditional result as a guide, construct an iteration scheme for all states, and prove that the value functions converge to some limits.

The algorithm for perpetual American options (used as a block at each time step in the backward induction procedure) is as follows.

• Choose the grid $\overrightarrow{x}$ (it might be necessary to use different grids in different states).
• In the cycle $j=1,\dots ,M$, calculate the initial approximation $V_{j,0}$ to the option value:

  $$V_{j,0}\left(\overrightarrow{x}\right)=E^{{\mathbb{Q}}_j}\left[{\int }_0^{{\tau }_j^{-}}{e}^{-(q_j+{\Lambda }_j)t}{g}_j\left(X_t\right)dt|X_0=\overrightarrow{x}\right]$$
• In the cycle $\ell =0,1,\dots ,$ for each $j=1,2\dots ,M$, calculate:

  $$V_{j,\ell +1}\left(\overrightarrow{x}\right)=E^{{\mathbb{Q}}_j}\left[{\int }_0^{{\tau }_j^{-}}{e}^{-(q_j+{\Lambda }_j)t}\left(g_j\left(X_t\right)+\sum _{k\ne j}{\lambda }_{jk}{V}_{k,\ell }\left(X_t\right)\right)dt|X_0=\overrightarrow{x}\right]$$
  Stop when $\parallel V_{·,\ell +1}-V_{·,\ell }\parallel \le ϵ$, where $ϵ$ is the error tolerance.

The limit as $\ell \to \infty$ exists because the sequence of option values in each state is increasing (the non-negativity of $g_j$ is needed for the proof).

7. Affine and Quadratic Term Structure Models

7.1. Affine Processes

Affine processes are used in models with stochastic volatility and stochastic interest rates. The adjective affine is explained by two almost equivalent properties:

(1)

The characteristic function of the transition density is of the form of an exponential function of an affine function of factors of the models with the coefficients depending on time to maturity $\tau =T-t$ and spectral parameter $\xi$:

$$\mathbb{E}\left[e^{i\langle \xi ,X_T\rangle }\right|X_t=x]=exp[\langle A(\tau ,\xi ),x\rangle +B(\tau ,\xi )],$$

(80)

where $\langle a,b\rangle ={\sum }_j{a}_j{b}_j$ [88]. If the stochastic interest rate $r_t$ is modeled as an affine function of the factors of the model, the state space must be enlarged as in the probabilistic version of the Feynman–Kac formula, and an additional factor $Y_t={\int }_0^t{r}_{s}ds$ IS added. The extended model remains affine, and the representation (80) becomes possible.

(2)

The coefficients of the stochastic differential equation (SDE) defining the process are affine functions of the state variable.

Property (1) allows one to calculate expectations $V(t,x)=\mathbb{E}\left[G\left(X_T\right)\right|X_t=x]$ as in the case of Lévy models. However, given SDE, it is necessary to (1) prove that the representation of the form (80) exists; and (2) calculate the matrix function $A(\tau ,x)$ and vector-function $B(\tau ,x)$.

The formal proof is straightforward [89]. Assuming that the Feynman–Kac theorem holds, write down the Cauchy problem

$$\begin{array}{ccc}\hfill ({\partial }_t+\mathcal{L})V(t,x)& =& 0,t<T,x\in {\mathcal{D}}_0,\hfill \end{array}$$

(81)

$$\begin{array}{ccc}\hfill V(T,x)& =& e^{i\langle \xi ,x\rangle },\hfill \end{array}$$

(82)

and where ${\mathcal{D}}_0$ and $\mathcal{L}$ are the interior of the state space $\mathcal{D}$ of the process and infinitesimal generator, respectively, substitute anzatz (80) into (81)–(82) and reduce the calculation of $(A,B)$ to the solution of the generalized system of Riccati equations associated with the model. The proof in the general case is unknown (see [2] for the discussion and bibliography).

An incomplete list of outstanding problems for affine models, which can be regarded as mathematical problems for a general well-defined class of PDO, is as follows (each of the problems is solved for particular models or certain subclasses of affine models but the complete answers, in full generality, are unknown—for a short overview of the extant results, see [2]):

• Prove the equivalence of (1) and (2) in the general case or for as wide a class of SDE with affine coefficients as possible. The difficulty stems, in particular, from the fact that the state space is of the form ${\left({\mathbb{R}}_{+}\right)}^m\times {\mathbb{R}}^{n-m}$, the infinitesimal generator degenerates at the boundary and the term of order 0 (“electric potential”) is an bounded affine function;
• Prove the Feynman–Kac theorem for the (backward) Cauchy problem and more general boundary problems.
• Derive general conditions for the explosion of the solution of the boundary problem (81) and (82) as $t\to T$.
• Study the domain of analyticity of the characteristic function (80) and its behavior at infinity, hence, the applicability of the conformal deformation technique. For partial results, see [2,15].
• Derive conditions on the parameters of affine interest rate models which ensure that the solution of (81) and (82) with $\xi =0$ (price of the discount bond) is bounded by 1. For partial results, see [23,25].
• Study the asymptotics of the solution of the (backward) Cauchy problem with the terminal condition $V(T,x)=G\left(x\right)$, as $T\to +\infty$. For partial results based on the eigenfunction expansion technique, see [90,91,92,93]. The main block is the generalized eigenfunction expansion of the essentially non-self-adjoint quadratic Hamiltonian. This is a special case of the same procedure in quadratic term structure models [91].

7.2. Wishart Models

Factors are naturally organized as a matrix rather than vector. One may regard Wishart models as square root models where the positive scalar stochastic factor is replaced by a positive-definite matrix. See [94,95] for a list of references. The outstanding problems are the same as for affine models.

7.3. Quadratic Term Structure Models (QTSM)

In the pure diffusion case, the infinitesimal generator of QTSM is of the form:

$$L=\langle \kappa (\theta -x),{\partial }_x\rangle +\frac{1}{2}\langle A{\partial }_x,{\partial }_x\rangle +\langle Bx,x\rangle +\langle d,x\rangle +d_0,$$

(83)

where $\kappa$ is an anti-stable matrix with real entries, matrix A is positive definite, B is semi-definite (both with real entries), $\theta ,d\in {\mathbb{R}}^n$, $d_0\in \mathbb{R}$. The characteristic function is an exponential of a quadratic function of x, with the coefficients depending on time to maturity $\tau$ and the spectral parameter $\xi$ and can be calculated solving the associated system of generalized Riccati equations. The generalized eigenfunction expansion can be calculated as well [90,91]. It is interesting that, for the parameters of QTSM documented in real financial markets, L is essentially non-self-adjoint, hence, not diagonalizable.

If the diffusion part of the infinitesimal generator is replaced with a PDO, then the exact calculation of the characteristic function and generalized eigenfunction expansion is not known. For asymptotic approximations, see [24,92].

7.4. Systems of Affine and Quadratic Term Structure Models

To the best of my knowledge, there are no general results in this direction, although in view of a huge body of publications in quantitative finance, one expects that some special cases have been considered.

8. Conclusions

In this paper, non-standard features of several basic problems arising in finance, insurance and economics are explained, and a group of related efficient methods (analytical and numerical) are outlined. The objective is to calculate the prices of contingent claims. In the general economic framework, the price V of a contingent claim is the expectation of a certain stochastic expression. In the case of Markov models, V is a function of time and the spot value $x=X_t$ of the underlying source of uncertainty. The formal application of Dynkin’s formula leads to a boundary problem for an integro-differential (pseudo-differential) equation of the form ${\partial }_t+L_X$, where $L_X$ is the infinitesimal generator of X (we allow for processes with killing/birth, hence, in the case of the Lévy model and the constant interest rate r, $L_X=-\psi \left(D_x\right)-r$, where $\psi \left(\xi \right)$ is the characteristic exponent of X). However, in the majority of cases of interest, the rigorous proof of the theorem of Feynman–Kac type is lacking. In the paper, several basic situations where proofs are available are considered in more detail; more general results are only outlined.

The following general features of the boundary problems discussed in this paper may be of general interest to specialists in PDE and PDO.

• In models with jumps, boundary conditions are non-local, whereas the standard boundary and co-boundary problems for PDE, PDO and fractional differential equations are local. See, e.g., [3,4,96,97]. Therefore, (1) one of the standard approaches to boundary problems, namely reduction to the boundary, cannot be used to reduce the dimension of the problem; (2) numerical methods which do not take the non-locality of the boundary conditions into account properly produce sizable errors; if the time horizon is large, the relative errors are, typically, very large.
• In popular pure diffusion models such as the Heston model, the operator degenerates at the boundary. The degeneration is sufficiently regular so that the generalization of the Boutet de Monvel calculus for degenerate elliptic operators [3] is applicable in a number of situations (interestingly, the infinitesimal generator in the Heston model is one of the basic examples in [3]). Formally, one can apply this calculus to boundary problems in models with jumps provided that the characteristic exponent of the jump part is a rational function. However, such an application would require the approximation of non-local boundary conditions by local ones. For a reduction to the boundary in applications to the Heston model and other basic diffusion models, see [98] and the bibliography therein.
• The degeneration and non-locality of the infinitesimal generators are the sources of fundamental difficulties for a rigorous proof of the Feynman–Kac theorem. One must establish certain regularity conditions of the solution for the proof. For applications of Dynkin’s formula, conditions are weaker than for applications of Ito’s formula, and, in some cases, general regularity results [3,5] can be used. However, the author is unaware of any general proof.
• In the case of Lévy models (PDO with constant symbols) and problems with flat boundaries, the probabilistic version of the Wiener–Hopf factorization technique can be used to derive an explicit formula for the price in the form of oscillatory integrals, and the analytic form of the same technique used to derive the same formula for the unique solution of the corresponding boundary problem thereby proving the Feynman–Kac theorem. In the case of problems with a curved boundary, the general regularity results and the proof of the Feynman–Kac theorem are unknown.
• The proof and study of regularity are especially non-trivial if the infinitesimal operator $L=-\psi \left(D_x\right)$ of a Lévy model is an elliptic PDO of order $\nu \in (0,1)$. Models of this kind are documented in the majority of empirical studies (if Lévy models are calibrated to the real data). We outlined approaches to study general boundary problems with operators of the form ${\partial }_t-\psi \left(D_x\right)-r$.
• In popular Lévy models, the solutions of the boundary problems are irregular at the boundary, hence numerical methods that (implicitly or explicitly) assume that the solution is more smooth than it is inevitably produce large errors.
• If the infinitesimal operator of a Lévy model is an elliptic PDO of order $\nu \in (0,1)$, then the smooth pasting principle for free boundary may fail.
• In many cases, the free boundary is discontinuous at the terminal date, which implies that a numerical method that assumes the continuity is bound to be inaccurate.

In the paper, proofs of facts 4–8 and several related efficient numerical methods for the solution are outlined. The main blocks are as follows.

• The interpretation of operators in the operator form from the Wiener–Hopf factorization as expectation operators under supremum and infimum processes, and explicit formulas for solutions of basic boundary problems under very mild restrictions on operators and boundary conditions, in the case of flat boundaries.
• The interpretation allows one to prove the convergence of general algorithms for pricing options:

  (a)

  With finite time horizon (stationary boundary problems) using maturity randomization (method of lines);

  (b)

  In regime-switching models;

  (c)

  Approximations of models with a stochastic interest rate and stochastic volatility by regime-switching models (systems of boundary problems);

  (d)

  Options with non-monotone and discontinuous payoffs, with applications to game-theoretical problems.

• In the interest of brevity, the refined version of the fast Fourier transform (FFT) and inverse Fourier transform (iFFT) introduced in [46,69] to price single and barrier options in Lévy models is not described in the paper. The refined version can be used to accurately evaluate options of various kinds. The advantage of refined FFT and iFFT as compared to existing versions, including the fractional FFT, stems from the flexibility of choices of grids of different length in the dual and state spaces, in each block of the numerical procedure. The final result is calculated using an almost optimal number of the standard FFT and iFFT blocks of the same (smaller) size.
• Instead, in the paper, we describe very fast and accurate methods for the numerical evaluation of integrals, in dimensions 1–4, based on the conformal deformation technique. The main idea is close to the idea of the saddle point method but the families of the contour deformations that we use allow one to relatively easily construct deformations with respect to several variables. The joint deformation of several contours above can be regarded as a further step in the realization of a general program of study of the efficiency of combinations of one-dimensional inverse transforms for high-dimensional inversions outlined in [99,100] with additional twists: the calculation of the Wiener–Hopf factors, which is necessary to price lookback and barrier options. The authors of [99,100] consider three main different one-dimensional algorithms for the numerical realization of the Bromwich integral (i) Fourier series expansions with the Euler summation; (ii) combinations of Gaver functionals; and (iii) deformation of the contour in the Bromwich integral, and discuss various methods of multi-dimensional inversion based on combinations of these three basic blocks. Our results imply that, for the purposes of multi-dimensional inversion, the class of deformations must be enlarged. In particular, in some practically important situations, deformations close to the steepest descent such as Talbot’s deformation $q=r\theta (cot\theta +i),-\pi <\theta <\pi$ [101] are not applicable, and one must resort to seemingly less efficient deformations. Note that the general conformal deformation technique that we develop is especially efficient in the case of highly oscillatory integrals.

We are grateful to the anonymous referee for the suggestion to include a short review of other methods. Naturally, essentially all methods developed to solve boundary problems for PDE and PDO can be used to price contingent claims, and many of these methods (if not all—the literature is huge) are used. We list and explain sources and types of errors of several groups of methods.

• “The Hilbert transform method” is used in backward induction procedures to price options of several types. Calculations are in the dual space. At each time step, operators of the form $\mathcal{F}{\mathbf{1}}_{(-\infty ,h)}{\mathcal{F}}^{-1}$ and $\mathcal{F}{\mathbf{1}}_{(h,+\infty )}{\mathcal{F}}^{-1}$ are expressed in terms of the Hilbert transform, and the latter is realized using the fast Hilbert transform. See [102,103,104] and the bibliographies therein. In these papers and other papers where the fast Hilbert transform is used, the grids of the same length in the state and dual spaces are used, which is presented as an advantage of the fast Hilbert transform approach. However, in many cases, the choice of grids of equal size leads to either very large errors or unnecessarily long grids and a large CPU time. See [46] for details and the explanation on how to use grids of various length in order to efficiently control discretization and truncation errors. Efficient methods of the numerical Fourier inversion described in the paper can be adjusted to the Hilbert transform. See [10] for an efficient numerical realization of operators $\mathcal{F}{\mathbf{1}}_{(-\infty ,h)}{\mathcal{F}}^{-1}$ and $\mathcal{F}{\mathbf{1}}_{(h,+\infty )}{\mathcal{F}}^{-1}$, applicable when these operators are applied only once. For an efficient numerical realization of operators generalizing the Hilbert transform, which is applicable in backward induction procedures (double spiral method), see [105].
• Variations of the straightforward application of the Fourier and inverse Fourier transform to European options, equivalently, the solution of Cauchy problems for parabolic PDO on the real line, namely, COS method and Carr–Madan method mentioned in Section 3.2, introduce additional unnecessary errors. The so-called Lewis–Lipton formula is the standard Fourier inversion formula with the prefixed line of integration. The choice of the line is non-optimal in the majority of cases; the conformal deformation method is much faster and more accurate. See [2] for numerical examples.
• In the CONV method [106,107], at each step of backward induction, an extremely inefficient interpolation procedure for the approximation of the value function at each time step is applied. In probabilistic terms, continuous distributions are approximated by discrete distributions supported on a uniform grid, a dual grid is chosen, and the calculations at each time step are reduced to the composition of FFT, multiplication by the array of values of the characteristic function, and iFFT. The procedure is simple but the errors are large. See [46] for the detailed analysis.
• The COS method is applicable (and has been applied) to price options of various kind. The essence of the method is an approximation of the kernel of the transition operator by a linear combination of cosines (hence the name). I find it difficult to find a sound mathematical argument in favor of this approximation. On the contrary, it is possible to indicate additional sources of errors and produce examples which demonstrate the inefficiency of COS. In backward induction procedures, the errors of COS accumulate very quickly, and pricing barrier options with even a moderate time horizon (0.5Y) is, essentially, impossible. See numerical examples in [2,13,14,16,17,18,19,20,31].
• In the PROJ method, the transition density of a random variable (equivalently, the kernel of the transition operator) is projected on a B-spline basis. See [20] for the bibliography, the discussion about the relative efficiency of COS, PROJ, the method in [46] which does not use an approximation of the transition kernel, and for an efficient procedure for the calculation of the projection coefficients using the sinh-acceleration. Note that the error of the approximation of the transition kernel is in the $H^2$-norm; hence, if the transition density has large derivatives or is non-smooth at the origin, which is the case of the VG model and Lévy models of order $\nu$ close to 0, then the errors of PRPJ can be very large.
• Approximations of value functions at each time step and filtering. In [16,20,46,69] (see also the bibliographies therein), the value function at each time step is approximated by piece-wise polynomials. In view of the irregularity of value functions near the boundary discussed in the paper, such an approximation introduces an error which can be controlled. See [16,60]. The approximation can be interpreted as a spectral filter, which is used in a number of publications to increase the speed of convergence. In [108], ad hoc spectral filters are used to increase the convergence of the integrals: “When Fourier techniques are employed to specific option pricing cases from computational finance with non-smooth functions, the so-called Gibbs phenomenon may become apparent. This seriously impacts the efficiency and accuracy of the pricing. For example, the Variance Gamma asset price process gives rise to algebraically decaying Fourier coefficients, resulting in a slowly converging Fourier series. We apply spectral filters to achieve faster convergence. Filtering is carried out in Fourier space; the series coefficients are pre-multiplied by a decreasing filter, which does not add significant computational cost. Tests with different filters show how the algebraic index of convergence is improved.” The quoted statement is correct. However, spectral filters are designed to regularize the results. The regularization of value functions results in serious errors in regions of paramount importance for risk management: near barrier and strike, close to maturity and for long dated options. For instance, close to the barrier or default boundary, the value can be overvalued or undervalued manifold. This remark is applicable to the applications of spectral filtering in [109] as well. Note that the conformal deformation technique allows one to eliminate the Gibbs phenomenon without sacrificing accuracy, and at a small CPU cost.
• Approximation of small jumps component by a diffusion. Cont and Volchkova [110] approximated the small jump component by a diffusion. In the result, a PDO of order $\nu <2$ is replaced with the sum $L^{ϵ}=ϵ{\partial }_x^2+\mu {\partial }_x+L_J$, where $ϵ>0$ is small, $\mu \in \mathbb{R}$ and $L_J$ is an integral operator with the kernel of class $L_1$; for an accurate approximation, the peak of the kernel has to be very high. After that, a standard implicit–explicit finite difference scheme is used to price barrier options. It is evident that if the infinitesimal generator L is a PDO of order $\nu <2$, hence, the derivative of the solution can be unbounded near the boundary, an approximation of L by $L^{ϵ}$ must lead to sizable errors near the boundary because the solution of the boundary problem becomes smooth up to the boundary. Numerous numerical examples in [111] have demonstrated the inaccuracy of the Cont–Voltchkova method. Note that the methods in [111] resemble but are less efficient than the method in [46,69].
• The approximation of KoBoL and other processes of infinite activity by HEJD model. In [38,40], the author constructed an HEJD model (without a special label attached) whilst keeping in mind to try such an approximation. For pricing American and barrier options, the advantage of HEJD is a simple explicit formula for the Wiener–Hopf factors in the case of positive values of the spectral parameter, derived in [38,40]. Unfortunately, for wide regions in the parameter space and large values of the spectral parameter which arise if a small time step in a backward induction procedure is used, an accurate approximation requires the use of HEJD with very large parameter values and high precision arithmetic is necessary. The reason is the same as in the case of the Cont–Voltchkova method. Due to this inefficiency, the author did not mention approximation of KoBoL by HEJD. Later, such an approximation was used in a number of publications, e.g., [112,113]. For a typical set of parameter values of KoBoL and moderate maturities, such an approximation can be very inaccurate at the distance of up to several percent of barrier. See [46] for numerical examples that illustrate the inefficiency of HEJD approximation. The problem of a large spectral parameter can be partially resolved using Richardson’s extrapolation. In applications to pricing barrier options, the convergence of Richardson’s extrapolation of arbitrary order is proven in [63]. Note that the technique of conformal deformations [10] is more efficient than approximation by HEJD, even in cases when the approximation is reasonably accurate.
• Approximation of underlying jump-diffusions with continuous time Markov chains. In [82,83,84,86], in the models with stochastic interest rates and/or stochastic volatility, the dynamics of additional factors is approximated by continuous time Markov chains. In the PDE language, a part of the infinitesimal generator is discretized; the result is a regime-switching Lévy model. At the first (discretization) step, a Markov chain with the infinite number of states (infinite grid) appears; at the second step, the infinite grid is truncated, and transition rates in a vicinity of the “boundary” of the truncated grid are adjusted so that the Markov chain remains the Markov chain without killing. Thus, the Dirichlet condition must be avoided. In the diffusion case, the adjustment can be interpreted as the discretization of the high contact condition ${\partial }_y^{2}V=0$; in the jump-diffusion case, the “discretized boundary condition” is non-local and more involved. Later, in a number of publications starting with [114], the approximation-by-continuous time Markov chain was used for more general Markov processes in 1D. In some publications, even the dynamics of Lévy factors was approximated by a continuous time Markov chain. Such an approximation is rather inefficient, especially if the tails decay slowly, and/or in the presence of barriers. Furthermore, in related publications, the discretized Dirichlet condition is used, which leads to significant errors of backward induction procedures with many steps.
• Eigenfunction expansion approach. In the case of diffusion models on the real line, there is a significant body of results obtained by V. Linetsky and their students (see, e.g., [115] and the bibliography therein). In [90,91,116], the generalized eigenfunction expansion is derived for solutions of the Cauchy problems in multi-factor models.
• Asymptotic methods. Due to the irregularity of solutions near the boundary, the asymptotic formulas are reasonably accurate only in a rather small vicinity of the boundary [62]: hence, they are rather useless for numerical purposes (although useful for the qualitative analysis). The same is true for asymptotics near maturity (in terms of time to maturity, short time asymptotics). The conformal deformation method can be used to calculate solutions close to maturity with high accuracy. For long time asymptotics, efficient methods can be derived using the eigenfunction expansion technique [92]. Note that there is a large body of the literature devoted to the study of the asymptotics of implied volatility close to maturity and far from maturity.
• Fast Gauss transform [117,118] can be efficiently used in certain diffusion models, and models with jumps of a special structure.
• For applications of finite elements to option pricing, see [119].
• Approximations based on purely probabilistic methods. The literature is huge. A typical feature is that the convergence of a method is proven without error bounds. A typical example is the Cont–Voltchkova method [110]. The proof of convergence is given; however, as the numerical examples in [111] demonstrate, in many cases, it is necessary to use extremely fine and long grids to satisfy the error tolerance of the order of one percent, at a very large CPU cost.
• Monte-Carlo simulations. The version of the Monte Carlo simulations that is closest to the methods of the present paper is based on the evaluation of the cumulative probability distribution function (cpdf) on an appropriate grid and interpolation. In applications to finance, the idea was used for the first time in [120]. Note that in [120] FFT and in a number of papers since FFT is used. As numerical examples in [31] demonstrate, the evaluation of the cpdf of Lévy processes using FFT leads to inaccurate results. The conformal deformations technique allows one to design much more accurate Monte-Carlo simulation procedures [9,31].